預(yù)測海藻數(shù)量

1、R語言實(shí)驗(yàn)設(shè)計(jì)題目：海藻數(shù)量數(shù)據(jù)預(yù)處理與數(shù)據(jù)探索性分析姓名：樊 海 濤 班級：12-1 學(xué)號(hào)：1207060106專業(yè)：統(tǒng) 計(jì) 學(xué) 預(yù)測海藻數(shù)量1.1基本任務(wù)：數(shù)據(jù)預(yù)處理、探索性數(shù)據(jù)分析和預(yù)測模型的構(gòu)建1.2問題描述與目標(biāo)某些高濃度的有害海藻對河流生態(tài)環(huán)境的強(qiáng)大破壞是一個(gè)嚴(yán)重問題，她們不僅僅破壞河流的生物，也破壞水質(zhì)。能夠檢測并在早期對海藻的繁殖進(jìn)行預(yù)測對提高河流質(zhì)量是很有必要的。針對這一問題的預(yù)測目標(biāo)，在大約一年的時(shí)間里，在不同時(shí)間內(nèi)收集了歐洲多條不同河流的水樣。對于每個(gè)水樣，測定了它們的不同化學(xué)性質(zhì)以及7種有害藻類的存在頻率，在水樣的收集過程中，也記錄了一些其他特征，如收集的季節(jié)、河流的大

2、小和水流的速度。本案例研究的動(dòng)機(jī)之一是化學(xué)檢測價(jià)格便宜，并且易于自動(dòng)化。而通過分析生物樣品來識(shí)別水中的藻類要涉及顯微鏡檢驗(yàn)，需要訓(xùn)練有素的工作人員，因此即昂貴又緩慢，因此，構(gòu)建一個(gè)可以基于化學(xué)性質(zhì)來準(zhǔn)確預(yù)測藻類的模型將有助于建立監(jiān)測有害藻類的廉價(jià)的自動(dòng)化系統(tǒng)。另一個(gè)目的是更好地了解影響藻類頻率的因素。也就是說，我們要了解藻類的頻率和水樣的某些化學(xué)性質(zhì)以及其他性質(zhì)（季節(jié)、河流類型等）是如何相關(guān)的。1.3數(shù)據(jù)說明來自于ERUDIT研究網(wǎng)絡(luò)，并被用于1999年的COIL國際數(shù)據(jù)分析競賽。在R軟件中可以直接獲得。本案例有兩個(gè)數(shù)據(jù)集，第一個(gè)數(shù)據(jù)集有200個(gè)水樣。名為：algae，更精確的說，該數(shù)據(jù)集的每

3、一條記錄是同一河流在該年的同一個(gè)季節(jié)的三個(gè)月內(nèi)收集的水樣的平均值。每條記錄由11個(gè)變量構(gòu)成。其中3個(gè)變量是名義變量，它們分別描述水樣收集的季節(jié)、收集樣品的河流大小和河水速度。余下的8個(gè)變量是所觀測水樣的不同化學(xué)參數(shù)，即最大PH值最小含氧量O2平均氯化物含量CI平均硝酸鹽含量NO3-平均氨含量NH4+平均正磷酸鹽含量PO43-平均磷酸鹽含量PO4平均葉綠素含量與這些參數(shù)相關(guān)的是7種不同有害藻類在相應(yīng)水樣中的頻率數(shù)目。并未提供所觀察藻類的名稱的相關(guān)信息。第二個(gè)數(shù)據(jù)集由140個(gè)額外觀測值構(gòu)成。名為：algae.sols，它們的基本結(jié)構(gòu)和第一個(gè)數(shù)據(jù)集一樣，但是它不包含7種藻類的頻率數(shù)目?？梢暈闇y試集。

4、本案例的主要目標(biāo)是預(yù)測140個(gè)水樣中7種海藻的頻率。在這種問題中，任務(wù)是建立預(yù)測模型，并預(yù)測在給定預(yù)測變量的取值時(shí)相應(yīng)的目標(biāo)變量的值。說明哪一個(gè)預(yù)測變量對目標(biāo)變量有較大的影響。1.4加載數(shù)據(jù)利用R軟件載入DMwR添加包，里面有我們需要的名為algae的數(shù)據(jù)框。這個(gè)數(shù)據(jù)框里含有前面提到的200個(gè)觀測值：> library(DMwR)> head(algae) season size speed mxPH mnO2 Cl NO3 NH4 oPO41 winter small medium 8.00 9.8 60.800 6.238 578.000 105.0002 spring sma

5、ll medium 8.35 8.0 57.750 1.288 370.000 428.7503 autumn small medium 8.10 11.4 40.020 5.330 346.667 125.6674 spring small medium 8.07 4.8 77.364 2.302 98.182 61.1825 autumn small medium 8.06 9.0 55.350 10.416 233.700 58.2226 winter small high 8.25 13.1 65.750 9.248 430.000 18.250 PO4 Chla a1 a2 a3 a

6、4 a5 a6 a71 170.000 50.0 0.0 0.0 0.0 0.0 34.2 8.3 0.02 558.750 1.3 1.4 7.6 4.8 1.9 6.7 0.0 2.13 187.057 15.6 3.3 53.6 1.9 0.0 0.0 0.0 9.74 138.700 1.4 3.1 41.0 18.9 0.0 1.4 0.0 1.45 97.580 10.5 9.2 2.9 7.5 0.0 7.5 4.1 1.06 56.667 28.4 15.1 14.6 1.4 0.0 22.5 12.6 2.9函數(shù)head()將顯示數(shù)據(jù)框的前6行。數(shù)據(jù)框的每一行代表一個(gè)觀測值。

7、1.5數(shù)據(jù)可視化和摘要鑒于開始我們對該領(lǐng)域一無所知，首先我們要了解一些數(shù)據(jù)的統(tǒng)計(jì)特性，為后面的數(shù)據(jù)處理與建模提供更多的信息。獲取數(shù)據(jù)統(tǒng)計(jì)特性的一個(gè)方法是獲取數(shù)據(jù)的描述性統(tǒng)計(jì)摘要。> summary(algae) season size speed mxPH autumn:40 large :45 high :84 Min. :5.600 spring:53 medium:84 low :33 1st Qu.:7.700 summer:45 small :71 medium:83 Median :8.060 winter:62 Mean :8.012 3rd Qu.:8.400 Max.

8、:9.700 NA's :1 mnO2 Cl NO3 Min. : 1.500 Min. : 0.222 Min. : 0.050 1st Qu.: 7.725 1st Qu.: 10.981 1st Qu.: 1.296 Median : 9.800 Median : 32.730 Median : 2.675 Mean : 9.118 Mean : 43.636 Mean : 3.282 3rd Qu.:10.800 3rd Qu.: 57.824 3rd Qu.: 4.446 Max. :13.400 Max. :391.500 Max. :45.650 NA's :2

9、NA's :10 NA's :2 NH4 oPO4 PO4 Min. : 5.00 Min. : 1.00 Min. : 1.00 1st Qu.: 38.33 1st Qu.: 15.70 1st Qu.: 41.38 Median : 103.17 Median : 40.15 Median :103.29 Mean : 501.30 Mean : 73.59 Mean :137.88 3rd Qu.: 226.95 3rd Qu.: 99.33 3rd Qu.:213.75 Max. :24064.00 Max. :564.60 Max. :771.60 NA's

10、 :2 NA's :2 NA's :2 Chla a1 a2 Min. : 0.200 Min. : 0.00 Min. : 0.000 1st Qu.: 2.000 1st Qu.: 1.50 1st Qu.: 0.000 Median : 5.475 Median : 6.95 Median : 3.000 Mean : 13.971 Mean :16.92 Mean : 7.458 3rd Qu.: 18.308 3rd Qu.:24.80 3rd Qu.:11.375 Max. :110.456 Max. :89.80 Max. :72.600 NA's :12

11、 a3 a4 a5 Min. : 0.000 Min. : 0.000 Min. : 0.000 1st Qu.: 0.000 1st Qu.: 0.000 1st Qu.: 0.000 Median : 1.550 Median : 0.000 Median : 1.900 Mean : 4.309 Mean : 1.992 Mean : 5.064 3rd Qu.: 4.925 3rd Qu.: 2.400 3rd Qu.: 7.500 Max. :42.800 Max. :44.600 Max. :44.400 a6 a7 Min. : 0.000 Min. : 0.000 1st Qu

12、.: 0.000 1st Qu.: 0.000 Median : 0.000 Median : 1.000 Mean : 5.964 Mean : 2.495 3rd Qu.: 6.925 3rd Qu.: 2.400 Max. :77.600 Max. :31.600 這個(gè)簡單函數(shù)立即給出了數(shù)據(jù)的的統(tǒng)計(jì)特征概括，對于名義變量，他給出了每個(gè)取值的變量的頻數(shù)。例如，從結(jié)果中可知冬季采集的水樣比其他季節(jié)更多，河流size為最大的有45個(gè)，河流的流速low比較少。對于數(shù)值型變量，R為我們提供了四分之一位數(shù)，中位數(shù)，均值，四分之三位數(shù)，極值等一系列信息。這些統(tǒng)計(jì)信息提供了變量值分布的初步信息，在變量有

13、缺失值的情況下，字符串NA后面的數(shù)值即為缺失值的個(gè)數(shù)，通過中位數(shù)，均值，四分位數(shù)的信息，我們可以了解數(shù)據(jù)分布的偏度和分散情況。且這些信息大多數(shù)都可以通過圖形來表達(dá)出來。我們看一個(gè)例子。> hist(algae$mxPH,probability = T) 上面程序繪制變量mxPH的直方圖。如下圖，變量mxPH的分布非常接近正態(tài)分布，它的值大多數(shù)都集中在變量的均值附近，我們做Q-Q來檢驗(yàn)該變量是否是正態(tài)分布，加載包c(diǎn)ar，使用里面的qq.plot()繪制Q-Q圖> library(car)> par(mfrow=c(1,2)> hist(algae$mxPH,prob=T,

14、xlab = "",ylim = 0:1,main = "Histogram of maximum pH value")> lines(density(algae$mxPH,na.rm=T)> rug(jitter(algae$mxPH)> qq.plot(algae$mxPH,main="Normal QQ plot of maximum pH")上圖的左邊是添加核密度曲線的直方圖，而在X軸附近繪制變量的實(shí)際值，從而容易辨別離群點(diǎn)，例如我們可以發(fā)現(xiàn)有兩個(gè)字顯著低于其它所有值。這種數(shù)據(jù)檢查是很重要的，因?yàn)樗艽_定數(shù)據(jù)

15、樣本中可能出現(xiàn)的錯(cuò)誤。右邊的Q-Q圖，他繪制了變量值與正態(tài)分布的理論分位數(shù)的散點(diǎn)圖，同時(shí)他給出正態(tài)分布的95%的置信區(qū)間的帶狀圖，從右上圖知，變量有幾個(gè)小的值明顯在95%置信區(qū)間之外，它們不服從正態(tài)分布。下面給出另一個(gè)數(shù)據(jù)檢驗(yàn)的例子，用它來檢驗(yàn)變量OP04：> boxplot(algae$oPO4,ylab="Orthophosphate")> rug(jitter(algae$oPO4),side = 2)> abline(h=mean(algae$oPO4,na.rm = T),lty=2)箱圖上方小橫線上面的小圓圈表示與其他值比較特別大的值，通常認(rèn)為

16、是離群值，這意味著箱線圖給出大量的信息，它不僅給出變量的中心趨勢，也給出了變量的發(fā)散情況和離群值。上圖中與X軸平行的直線，是變量的均值所在的位置，將均值線與中位數(shù)線進(jìn)行比較，就可以知道變量的多個(gè)離群值使得作為變量中心的均值產(chǎn)生了扭曲。變量oPO4的分布集中在較小的觀測值范圍，因此分布為正偏。大多數(shù)水樣的oPO4值都比較低，也有幾個(gè)水樣的觀測值較高，甚至特別高。有時(shí)，當(dāng)有離群值時(shí)，需要確定哪些有離群值的觀測，可以使用圖形法，例如繪制NH4的值，將會(huì)注意到一個(gè)特別大的值，我們可以使用一下方法識(shí)別特大值相應(yīng)的水樣：> plot(algae$NH4,xlab = "")>

17、; abline(h=mean(algae$NH4,na.rm = T),lty=1)> abline(h=mean(algae$NH4,na.rm = T)+sd(algae$NH4,na.rm = T),lty=2)> abline(h=median(algae$NH4,na.rm = T),lty=3)> identify(algae$NH4)上面函數(shù)第一條繪制變量的所有值，調(diào)用函數(shù)abline繪制三條有用的直線，第一條為均值，第二條為均值加一個(gè)標(biāo)準(zhǔn)差，第三條為中位數(shù)，對于離群值的識(shí)別盡管這三條不是必須的，但它們能提供有用的信息。 從上圖中知道，水樣153的NH4值是一

18、個(gè)極大異常值，20，34,35,88,89,153等水樣中NH4值都異常高，可以把這些水樣選擇出來，看看他們與別的水樣中海藻有什么種類有什么差異。下面給出幾種其他類型的數(shù)據(jù)檢驗(yàn)的例子，這些例子應(yīng)用R的添加包lattice，做Trellis圖形。假設(shè)需要研究海藻變量a1的值得分布，然而，這里需要分布如何依賴于其他變量，就需要新的變量，新的工具。條件繪圖是依賴于某個(gè)特定因子的圖形表示，因子是一個(gè)為一個(gè)取值為有限集合的名義變量。例如，對于size的不同取值，可以繪制變量a1的一組箱圖。每個(gè)箱圖是對應(yīng)于變量size的某個(gè)特定值的水樣子集。通過這些箱線圖可以研究名義變量size如何影響變量a1值得分布。

19、bwplot(sizea1,data=algae,ylab = "River Size",xlab = "Algal A1")從圖知，在規(guī)模較小的河流中，海藻a1的頻率越高。這種箱圖的另外一個(gè)類型是分為箱圖，它可以給出繪制變量的更多信息，R添加包Hmisc可以繪制分位箱圖，下面繪制上面例子中a1變量的條件分位箱圖：> library(Hmisc)> bwplot(sizea1,data=algae,panel = panel.bpplot,probs=seq(.01,.49,by=.01),datadensity=T,ylab="R

20、iver Size",xlab = "Algal A1")上圖中的點(diǎn)代表代表不同大小的河流中海藻頻數(shù)的均值，而圖中的豎線分別代表變量的第一份位數(shù)，中位數(shù)和第三分位數(shù)。圖中的小豎線代表數(shù)據(jù)的真實(shí)取值，這些值分布信息由分位數(shù)圖來體現(xiàn)。分位數(shù)箱圖提供的信息要多于傳統(tǒng)的箱圖。例如我們可以確認(rèn)上面的觀測結(jié)論，小型河流有更高頻率的海藻，但我們也觀察到小型河流的海藻頻率的分布比其他類型的河流的海藻頻率的分布分散。這種類型的條件繪圖不局限于名義變量，也不局限于單個(gè)因子。只需把連續(xù)變量離散化，也同樣可以進(jìn)行條件繪圖，下面給出兩因子條件繪圖的例子?？紤]變量a3在給定變量season和

21、變量mnO2下的條件繪圖，變量mnO2是一個(gè)連續(xù)變量，代碼如下：> m1<-equal.count(na.omit(algae$mnO2),number=4,overlap=1/5)> stripplot(seasona3|m1,data = algae!is.na(algae$mnO2),)得到的圖形如下，上面的代碼第一行是調(diào)用函數(shù)equal.Count()對連續(xù)變量mnO2離散化，把該變量轉(zhuǎn)化為因子類型。參數(shù)number設(shè)置需要的區(qū)間的個(gè)數(shù)。每個(gè)區(qū)間的觀測值的個(gè)數(shù)相等。 第二個(gè)調(diào)用函數(shù)stripplot()，該函數(shù)是lattice包中的一個(gè)繪圖函數(shù)，他根據(jù)另一個(gè)變量sea

22、son把變量的實(shí)際值繪制到不同的圖形中。然后對變量mnO2的每個(gè)不同的區(qū)間繪制不同的圖形。這些區(qū)間按照從左到右，從上到下的順序來排列。左下方對應(yīng)的是較小的mnO2值。1.6數(shù)據(jù)缺失處理在許多水樣中，一些變量含有缺失值。這種情形在現(xiàn)實(shí)問題中非常普遍，這會(huì)導(dǎo)致一些不能處理缺失值的分析方法無法應(yīng)用。當(dāng)我們處理含有缺失值的數(shù)據(jù)時(shí)，可以運(yùn)用以下幾種最常見的策略：將含有缺失值得案例剔除根據(jù)變量之間的相關(guān)關(guān)系填補(bǔ)缺失值根據(jù)案例之間的相似性填補(bǔ)缺失值使用能夠處理缺失值數(shù)據(jù)的工具最后一種方式是最嚴(yán)格的，因?yàn)樗薅晌覀兛梢允褂玫墓ぞ?。然而，?dāng)對這些可以處理缺失值的數(shù)據(jù)挖掘工具有信心時(shí)，它們可能是一個(gè)好的選擇。在

23、下面的章節(jié)中我們將通過舉例來演示如何在R中實(shí)現(xiàn)上述處理缺失值的方法。如果你決定試一試本書中的代碼，你需要了解它們不是互相補(bǔ)充的，由于每節(jié)都用不同的方法來處理這些數(shù)據(jù)，這意味著當(dāng)你選擇其他方式來處理這些缺失時(shí)，需要重新讀取含有所有缺失值的原始數(shù)據(jù)。重新讀取數(shù)據(jù)的代碼如下：library(DMwR)data(algae)16.1 將缺失部分剔除剔除含有缺失數(shù)據(jù)的記錄非常容易實(shí)現(xiàn)，尤其是當(dāng)這些記錄所占的比例在可用數(shù)據(jù)集中非常小的時(shí)候，這個(gè)選擇就比較合理。在剔除某些變量中至少含有一個(gè)缺失數(shù)據(jù)的所有觀測值時(shí)，最好先檢查觀測值，或者至少得到這些觀測值的個(gè)數(shù)，例如：algae!complete.cases(

24、algae),nrow(algae!complete.cases(algae),)1 16 16個(gè)含缺失的樣本函數(shù)complete.cases() 產(chǎn)生一個(gè)布爾值向量，該向量的元素個(gè)數(shù)與algae數(shù)據(jù)框中的行數(shù)相同，如果數(shù)據(jù)框的相應(yīng)行中不含NA 值（即為一個(gè)完整的觀測值），函數(shù)返回值就是TRUE。前面提到過的“!”運(yùn)算符，它是取邏輯否，因此上述指令顯示了含有缺失值得水樣記錄。為了從數(shù)據(jù)框中剔除16個(gè)樣本，我們可以簡單的輸入：algae <- na.omit(algae)即使我們決定不使用剔除所有包含缺失值記錄的極端方法，我們也可以剔除某些觀測值。因?yàn)檫@些樣本的缺失值太多，所以它們幾乎是無

25、用的樣本，如果采用復(fù)雜的方法來填補(bǔ)缺失值，就會(huì)導(dǎo)致較大偏差。需要注意的是，如果執(zhí)行了前面的命令，你需要重新讀取數(shù)據(jù)，因?yàn)檫@個(gè)指令已經(jīng)剔除了所有缺失數(shù)據(jù)，所以接下來的命令就沒有意義! 觀察這個(gè)樣本中的數(shù)據(jù)，我們可以看到第62條和第199條記錄中的11個(gè)解釋變量有6個(gè)是缺失值。這種情況下，最好是剔除它們：algae <- algae-c(62,199),在有些問題中，由于大量記錄中含有缺失值，用上面的觀察方法來檢查數(shù)據(jù)的缺失值是不可行的，所以需要找到缺失值較多的樣本所在的行。下面的代碼可以找到海藻數(shù)據(jù)集中每行數(shù)據(jù)的缺失值個(gè)數(shù)：apply(algae,1,function(x) sum(is.

26、na(x)函數(shù)apply（）屬于R中功能非常強(qiáng)大的一類函數(shù)。這類函數(shù)又稱為元函數(shù)，它們可以在某些條件下對對象應(yīng)用其他函數(shù)。對函數(shù)apply（）而言，它把一個(gè)函數(shù)應(yīng)用到數(shù)據(jù)框的每一行。這個(gè)被應(yīng)用的函數(shù)在apply（）函數(shù)的第三個(gè)參數(shù)中給出，對數(shù)據(jù)框的每一行都分別調(diào)用該函數(shù)。在這個(gè)案例中我們使用一個(gè)臨時(shí)函數(shù)。它只有在調(diào)用apply（）函數(shù)時(shí)才存在。另外，函數(shù)apply（）的第三個(gè)參數(shù)也可以是一個(gè)正常的函數(shù)的函數(shù)名。臨時(shí)函數(shù)的功能是計(jì)算對象中NA的數(shù)量。在R中邏輯值TRUE 等于數(shù)值1，邏輯值FALSE 等于數(shù)值0，這意味著當(dāng)加一個(gè)布爾值向量時(shí)，得到向量中取值為TRUE的元素的個(gè)數(shù)。根據(jù)以上代碼，可

27、以編寫一個(gè)程序找到algae中含有給定數(shù)目缺失值的行。可以如下應(yīng)用該函數(shù)：data(algae)manyNAs(algae,0.2)1 62 199只有在前面操作中剔除了缺失值較多的幾行數(shù)據(jù)后才需要調(diào)用函數(shù)data（）。函數(shù)manyNA是（）的功能是找到缺失值個(gè)數(shù)大于列數(shù)20%的行。在第二個(gè)參數(shù)中可以設(shè)置一個(gè)精確的列數(shù)作為界限。因此，用下面的代碼就無須知道含有缺失值較多的行的具體數(shù)量。algae <- algae-manyNAs(algae),以上操作可以把樣本中62，199號(hào)水樣數(shù)據(jù)剔除。1.6.2 用高頻率值來填補(bǔ)缺失值填補(bǔ)含有缺失值記錄的另一種方法是嘗試找到這些缺失值最可能的值，同

28、樣，這里有多種策略可供選擇，不同策略對逼近程度和算法復(fù)雜度的權(quán)衡不同，填補(bǔ)缺失數(shù)據(jù)最簡單和快捷的方法是使用一些代表中心趨勢的值，代表中心趨勢的值反映了變量分布的最常見的值，因此中心趨勢值是最自然的選擇。有多個(gè)代表數(shù)據(jù)中心趨勢的指標(biāo)，例如平均值，中位數(shù)，眾數(shù)等。最合適的選擇由變量的分布決定。對于接近正態(tài)的分布來說，所有的觀測值都較好地聚集在平均值周圍，平均值數(shù)就是最佳選擇。然而，對于偏態(tài)分布，或者離群值的變量來說，選擇平均值就不好。偏態(tài)分布的大部分值都聚集在變量分布的一側(cè)，因此平均值不能作為最常見的代表。另一方面，離群值（極值）的存在會(huì)扭曲平均值，這就導(dǎo)致了平均值不具有代表性的問題。因此，在對變

29、量分布進(jìn)行檢查之前選擇平均值作為中心趨勢的代表是不明智的，例如某些R的繪制工具。對偏態(tài)分布或者有離群值的分布而言，中位數(shù)是更好的代表數(shù)據(jù)中心趨勢的指標(biāo)。比如，樣本algae48,中變量mxPH有缺失值。由于該變量分布近似正態(tài)分布我們可以選用平均值來填補(bǔ)這個(gè)“洞”，計(jì)算方法如下：algae48,'mxPH' <- mean(algae$mxPH,na.rm=T)這里，函數(shù)mean（）計(jì)算數(shù)值向量的平均值，參數(shù)na.rm=T 使計(jì)算時(shí)忽略缺失數(shù)據(jù)。大多數(shù)時(shí)候采用一次填補(bǔ)一列中的所有缺失值而不是像上面那樣一行一行地逐個(gè)填補(bǔ)。以變量Chla為例，這個(gè)變量在第12行上有缺失值。另外

30、，這也是平均值不能代表大多數(shù)變量值的一種情況。事實(shí)上，Chla的分布偏向于較低的數(shù)值，并且它有幾個(gè)極端值，這些都使得平均值不能代表大多數(shù)的變量值。因此，我們使用中位數(shù)來填補(bǔ)這一類的缺失值：algaeis.na(algae$Chla),'Chla' <- median(algae$Chla,na.rm=T)由于缺失值的存在會(huì)導(dǎo)致某些方法不能使用，所以使用上面的方法填補(bǔ)缺失值通常也認(rèn)為不是很好的方法。雖然上述的方法速度快，特別適用于大數(shù)據(jù)集，但是它可能導(dǎo)致較大的數(shù)據(jù)偏差，影響后期的數(shù)據(jù)分析工作。然而，使用無偏差方法來尋找最佳數(shù)據(jù)調(diào)補(bǔ)值復(fù)雜，對于大型數(shù)據(jù)挖掘問題可能并不適用。1

31、.6.3 通過變量的相關(guān)關(guān)系來調(diào)補(bǔ)缺失值另一種獲取缺失值較少偏差估計(jì)值的方法是探尋變量之間的相關(guān)關(guān)系。比如，通過變量值之間的相關(guān)關(guān)系，能夠發(fā)現(xiàn)某變量與mxPH高度相關(guān)。這可以使我們得到含有缺失值的第48條樣本的更可能的填補(bǔ)值。這比之前使用平均值的方法將更勝一籌。> options(digits = 1)> cor(algae,4:18,use="complete.obs") mxPH mnO2 Cl NO3 NH4 oPO4 PO4 Chla a1mxPH 1.00 -0.10 0.15 -0.17 -0.15 0.090 0.10 0.43 -0.16mnO2

32、 -0.10 1.00 -0.26 0.12 -0.08 -0.394 -0.46 -0.13 0.25Cl 0.15 -0.26 1.00 0.21 0.07 0.379 0.45 0.14 -0.36NO3 -0.17 0.12 0.21 1.00 0.72 0.133 0.16 0.15 -0.25NH4 -0.15 -0.08 0.07 0.72 1.00 0.219 0.20 0.09 -0.12oPO4 0.09 -0.39 0.38 0.13 0.22 1.000 0.91 0.11 -0.39PO4 0.10 -0.46 0.45 0.16 0.20 0.912 1.00 0.

33、25 -0.46Chla 0.43 -0.13 0.14 0.15 0.09 0.107 0.25 1.00 -0.27a1 -0.16 0.25 -0.36 -0.25 -0.12 -0.395 -0.46 -0.27 1.00a2 0.34 -0.07 0.08 0.02 -0.04 0.124 0.13 0.37 -0.26a3 -0.03 -0.24 0.08 -0.09 -0.11 0.006 0.03 -0.06 -0.11a4 -0.18 -0.38 0.14 -0.01 0.27 0.382 0.41 -0.09 -0.09a5 -0.11 0.21 0.15 0.21 0.0

34、2 0.122 0.16 -0.07 -0.27a6 -0.17 0.19 0.17 0.54 0.40 0.003 0.05 0.01 -0.26a7 -0.17 -0.10 -0.04 0.08 -0.03 0.026 0.08 0.02 -0.19 a2 a3 a4 a5 a6 a7mxPH 0.34 -0.027 -0.18 -0.11 -0.173 -0.17mnO2 -0.07 -0.235 -0.38 0.21 0.189 -0.10Cl 0.08 0.077 0.14 0.15 0.169 -0.04NO3 0.02 -0.092 -0.01 0.21 0.544 0.08NH

35、4 -0.04 -0.113 0.27 0.02 0.401 -0.03oPO4 0.12 0.006 0.38 0.12 0.003 0.03PO4 0.13 0.032 0.41 0.16 0.053 0.08Chla 0.37 -0.063 -0.09 -0.07 0.010 0.02a1 -0.26 -0.108 -0.09 -0.27 -0.262 -0.19a2 1.00 0.010 -0.18 -0.19 -0.134 0.04a3 0.01 1.000 0.03 -0.14 -0.197 0.04a4 -0.18 0.033 1.00 -0.10 -0.085 0.07a5 -

36、0.19 -0.142 -0.10 1.00 0.389 -0.05a6 -0.13 -0.197 -0.08 0.39 1.000 -0.03a7 0.04 0.039 0.07 -0.05 -0.030 1.00函數(shù)cor（）的功能是產(chǎn)生變量之間的相關(guān)值矩陣（因?yàn)榍?個(gè)變量時(shí)名義變量，所以計(jì)算相關(guān)值時(shí)不考慮它們）。設(shè)定參數(shù)use=“complete.obs”時(shí)，R在計(jì)算相關(guān)值時(shí)忽略含有NA的記錄。相關(guān)值在1（或-1）周圍表示相應(yīng)的兩個(gè)變量之間有強(qiáng)正（或負(fù)）線性相關(guān)關(guān)系。然后其他R函數(shù)可以得到變量間線性相關(guān)的近似函數(shù)形式，它可以讓我們通過一個(gè)變量的值計(jì)算出另一個(gè)變量的值。函數(shù)cor（）的輸出

37、結(jié)果并不是很清晰，可用symnum來改善結(jié)果的輸出形式.> symnum(cor(algae,4:18,use="complete.obs") mP mO Cl NO NH o P Ch a1 a2 a3 a4 a5 a6 a7mxPH 1 mnO2 1 Cl 1 NO3 1 NH4 , 1 oPO4 . . 1 PO4 . . * 1 Chla . 1 a1 . . . 1 a2 . . 1 a3 1 a4 . . . 1 a5 1 a6 . . . 1 a7 1 attr(,"legend")1 0 0.3 . 0.6 , 0.8 + 0.9

38、* 0.9 B 1這種用符號(hào)表示相關(guān)值得方法更為清晰，特別是對于大的相關(guān)矩陣。在本案例的數(shù)據(jù)中，大多數(shù)變量之間是不相關(guān)的，然而，有兩個(gè)例外：變量NH4和NO3之間，變量PO4和oPO4之間。后兩個(gè)變量之間的相關(guān)值很高（大于0.9）。變量NH4和NO3之間的相關(guān)性不是特別明顯（為0.72），因此，根據(jù)它們來確定缺失數(shù)據(jù)是危險(xiǎn)的。此外，因?yàn)闃颖?2和樣本199有太多的變量含有缺失值，所以如果剔除它們，樣本中的變量NH4和NO3就沒有缺失值了。至于變量PO4和 oPO4，它們之間相關(guān)性可以幫助填補(bǔ)這兩個(gè)變量的缺失值。為了達(dá)到這個(gè)目標(biāo)，我們需要找到這兩個(gè)變量之間的線性關(guān)系，方法如下：> data

39、(algae)> algae<-algae-manyNAs(algae),> lm(PO4oPO4,data = algae)Call:lm(formula = PO4 oPO4, data = algae)Coefficients:(Intercept) oPO4 42.90 1.29 述該函數(shù)。線性模型是 PO4=1.29oPO4+42.90 如果這兩個(gè)變量不是同時(shí)有缺失值，那么可以用這個(gè)公式計(jì)算這些變量的缺失值。在剔除樣本62和樣本199后，還剩下一個(gè)樣本28 在PO4上有缺失值，可以簡單地使用上面的線性關(guān)系計(jì)算缺失值的填補(bǔ)值：> algae28,'PO4

40、' <- 42.897 + 1.293 * algae28,'oPO4'> algae28,"PO4"1 48通過變量之間的相關(guān)關(guān)系求出水樣28，PO4的缺失值填補(bǔ)為48。對線性關(guān)系的研究使我們能夠填充一些新的缺失值。然而還有幾個(gè)觀測值含有缺失值?？梢栽囍剿靼咐龜?shù)據(jù)中含有缺失值的變量和名義變量之間的關(guān)系，這可以通過R添加包lattice中的函數(shù)來繪制條件直方圖來進(jìn)行。繪圖代碼如下：> histogram(mxPH|season,data=algae)上圖是繪制在不同的季節(jié)變量mxPH的直方圖。mxPH直方圖在不同季節(jié)沒有顯著的不同

41、，下面對河流的大小size進(jìn)行上面類似的分析，執(zhí)行命令：> histogram(mxPH|season,data=algae)那么從得到的直方圖知，較小的河流有較小的mxPH值。對這種相關(guān)性研究，可以擴(kuò)展到多個(gè)名義變量：> stripplot(sizemxPH|speed,data=algae,jitter=T)上面的指令結(jié)果如下圖，參數(shù)jitter=T說明Y的變量值進(jìn)行小范圍的隨機(jī)排列，這些可以避免相同值之間的互相重疊而導(dǎo)致導(dǎo)致失去具有某些特定值的觀察值集中度的信息。這種類型的分析可以應(yīng)用到其它含有缺失值的變量中，然而這種分析是一種繁瑣的過程，有大量的變量組合需要分析。適合名義變

42、量較少的數(shù)據(jù)集的分析。1.6.4通過探索案例之間的相似性來填補(bǔ)缺失值不同于探索數(shù)據(jù)集列之間的相關(guān)性，嘗試使用行之間的相似性來填充缺失值。我們可以使用這種方法來填補(bǔ)剔除了那兩個(gè)含有太多的NA值得樣本外的其它缺失數(shù)據(jù)。如上所描述的方法假如兩個(gè)水樣是相似的，其中某些變量有缺失值，那么該缺失值很可能與另一個(gè)水樣的值是相似的。為了使用這種直觀的方法，首先定義相似性的概念。相似性經(jīng)常由描述觀察值的多元度量空間的變量所定義。在文獻(xiàn)中有許多度量相似性的指標(biāo)，常用的是歐式距離。這個(gè)距離可以非正式的定義為任何兩個(gè)案例的觀察值之差的平方和，計(jì)算公式如下：dx,y= 1p(xi-yi)2下面要描述的方法是使用這種度量

43、來尋找與任何含有缺失值的案例最相似的10個(gè)水樣，并用它們來填充缺失值。我們考慮兩種應(yīng)用這些值的方法。第一種方法簡單地計(jì)算這10個(gè)最相近的案例的中位數(shù)并用這個(gè)中位數(shù)來填充這些缺失值。第二種方法是采用這些相似數(shù)據(jù)的加權(quán)均值。權(quán)重的大小隨著距待填補(bǔ)缺失值的個(gè)案的距離為d，則它的值在加權(quán)平均中的權(quán)重為：wd=e-d上面的方法可以通過添加包函數(shù)knnImputation()來實(shí)現(xiàn)，這個(gè)函數(shù)用一個(gè)歐式距離的變種來找到距離任何個(gè)案最近的k個(gè)鄰居。在計(jì)算距離時(shí)都要對數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化。下面介紹如何使用knnImputation()函數(shù)> algae<-algae-manyNAs(algae),>

44、 algae<-knnImputation(algae,k=10)> algae如果使用中位數(shù)來填補(bǔ)缺失值，可以使用下面的代碼> algae<-knnImputation(algae,k=10,meth = "median")總之使用這些簡單地操作，數(shù)據(jù)中不再含有缺失值NA，未使用R的其他函數(shù)進(jìn)行分析做好充分的準(zhǔn)備。當(dāng)決定使用哪種方法來填充缺失值時(shí)，大多數(shù)時(shí)候根據(jù)你所分析的領(lǐng)域知識(shí)來確定，根據(jù)個(gè)案之間的相似性來填補(bǔ)缺失看起來合理，但這種方法也存在其他問題，例如可能存在不相關(guān)的變量扭曲相似性，甚至造成大型數(shù)據(jù)集的計(jì)算特別復(fù)雜的問題。另外，對于這些大數(shù)據(jù)

45、集問題，可以通過隨機(jī)抽取樣本的方法來計(jì)算他們之間的相關(guān)性。1.6獲取預(yù)測模型本案例主要研究的目的是預(yù)測140個(gè)水樣中7種海藻的出現(xiàn)頻率。假設(shè)海藻頻率是數(shù)值型數(shù)據(jù)因此可以考慮進(jìn)行回歸分析，簡單地說，預(yù)測任務(wù)是建立一個(gè)模型來找到一個(gè)數(shù)值變量和一組解釋變量之間的關(guān)系。這個(gè)模型既可以根據(jù)未來解釋變量的值來預(yù)測目標(biāo)變量，也可以更好的理解問題中各個(gè)變量之間的聯(lián)系。下面將研究適用于海藻的兩種的不同的模型，多元線性回歸模型和回歸樹模型。這兩種方法是應(yīng)用于回歸的較好的模型，他們對逼近的回歸函數(shù)的形式有完全不同的假設(shè)，兩類模型都易于解釋。這兩個(gè)問題以不同的方式解決缺失值問題。R中的線性回歸不能使用帶有缺失值的數(shù)據(jù)

46、集，而回歸樹模型可以很自然的處理這些帶缺失值得數(shù)據(jù)。對線性回歸模型我們運(yùn)用1.5節(jié)中描述的方法進(jìn)行數(shù)據(jù)的預(yù)處理，然后應(yīng)用線性回歸模型。在回歸樹模型中，我們將直接使用這些帶缺失值的數(shù)據(jù)。在下面的分析中，假定140個(gè)測試水樣的目標(biāo)變量的真實(shí)值是未知的，等我們處理完之后再與這些測試水樣的正確數(shù)據(jù)對照，看我們處理的正確率。2.6.1多元線性回歸多元線性分析是最常用的統(tǒng)計(jì)分析方法，該模型給出了一個(gè)有關(guān)目標(biāo)變量與一組解釋變量關(guān)系的線性函數(shù)。在多元線性回歸中沒有處理缺失值的方法。因此，這里將應(yīng)用根據(jù)訓(xùn)練集數(shù)據(jù)個(gè)案的相似性方法來填補(bǔ)缺失值。要注意的是，在使用這些方法之前，首先應(yīng)該移除第62條和第199條水樣的

47、記錄，因?yàn)檫@兩天記錄中11個(gè)變量6個(gè)是缺失值的。一下代碼是獲取一個(gè)不含缺失值的數(shù)據(jù)框：> algae<-algae-manyNAs(algae),> clean.algae<-knnImputation(algae,k=10)> clean.algae在運(yùn)行以上的代碼后，得到的數(shù)據(jù)框clean.algae將不含缺失值。接下來，將建立一個(gè)用于預(yù)測海藻頻率的線性回歸模型：> lm.a1<-lm(a1.,data=clean.algae,1:12)> summary(lm.a1)Call:lm(formula = a1 ., data = clean.

48、algae, 1:12)Residuals: Min 1Q Median 3Q Max -37.68 -11.89 -2.57 7.41 62.19 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 42.94206 24.01088 1.79 0.0754 . seasonspring 3.72698 4.13774 0.90 0.3689 seasonsummer 0.74760 4.02071 0.19 0.8527 seasonwinter 3.69295 3.86539 0.96 0.3406 size

49、medium 3.26373 3.80205 0.86 0.3918 sizesmall 9.68214 4.17997 2.32 0.0217 * speedlow 3.92208 4.70631 0.83 0.4057 speedmedium 0.24676 3.24187 0.08 0.9394 mxPH -3.58912 2.70353 -1.33 0.1860 mnO2 1.05264 0.70502 1.49 0.1372 Cl -0.04017 0.03366 -1.19 0.2343 NO3 -1.51124 0.55134 -2.74 0.0067 *NH4 0.00163

50、0.00100 1.63 0.1052 oPO4 -0.00543 0.03988 -0.14 0.8918 PO4 -0.05224 0.03075 -1.70 0.0911 . Chla -0.08802 0.08000 -1.10 0.2727 -Signif. codes: 0 * 0.001 * 0.01 * 0.05 . 0.1 1Residual standard error: 18 on 182 degrees of freedomMultiple R-squared: 0.373,Adjusted R-squared: 0.321 F-statistic: 7.22 on 15 and 182 DF, p-value: 2.44e-12在解釋函數(shù)summary()應(yīng)用到線性模型對象所得到的信息之前，先介紹R處理3個(gè)名義變量。當(dāng)像上面一樣進(jìn)行模型構(gòu)建時(shí)R會(huì)生成一組鋪助變量，即對每個(gè)有k個(gè)水平的因子變量，R會(huì)生成k-1個(gè)鋪助變量。這些鋪助變量的值為0或者為1，,當(dāng)鋪助變量為1表明該因子值出現(xiàn)，同時(shí)表明其他鋪助變量的值為0.如果所有這k-1個(gè)鋪助變量的取值都為0，則表明因子變量的取值為第k個(gè)剩余的值。對得到的線性模型對象應(yīng)用函數(shù)summary()，將給出所建模型的一些診斷信息。首先是有關(guān)線性模型中數(shù)據(jù)擬合殘差的。殘差應(yīng)該是均值為0并且為這個(gè)態(tài)分布?？梢允褂胮lot(l