使用R軟件預測海藻數量李強強2013.112/3/20231R軟件R是一套完整的數據處理、計算和制圖軟件系統。其功能包括：數據存儲和處理系統；數組運算工具（其向量、矩陣運算方面功能尤其強大）；完整連貫的統計分析工具；優秀的統計制圖功能；簡便而強大的編程語言：可操縱數據的輸入和輸出，可實現分支、循環，用戶可自定義功能。R在語義上是函數設計語言。它允許在“語言上計算”。這使得它可以把表達式作為函數的輸入參數，而這種做法對統計模擬和繪圖非常有用。R是一個免費的自由軟件。本案例使用的是R的3.0.1版。2/3/20232背景描述某些高濃度的有害藻類嚴重破壞著河流的生態環境，因此，能夠監測并及早對海藻的繁殖進行預測對提高河流的質量是很有必要的。在約一年時間內，在不同的時間收集了多條不同河流的水樣。每個水樣測定了它們不同的化學性質和7種有害藻類的存在頻率。還記錄了如收集的季節、河流大小和水流速度。案例研究動機：1.化學監測相對人工檢測價格便宜，且易于自動化。2.更好地了解藻類的頻率和水樣的某些化學性質以及其他特性(如季節、河流類型等)是如何相關的。2/3/202330海藻數據第一個數據集：包含200個水樣。每條記錄是一條河流在該年的同一季節的三個月內收集的水樣的平均值。其中，每條記錄由11個變量組成。其中3個名義變量：水樣收集的季節、河流大小和河水速度。8個變量是水樣的不同化學參數：最大pH值、最小含氧量、平均氯化物含量、平均硝酸鹽含量、平均氨含量、平均正磷酸鹽含量、平均磷酸鹽含量和平均葉綠素含量。與之相關的是7種不同的有害藻類的頻率數目。第二個數據集：140個不含7種藻類頻率數目的測試集。2/3/202341數據加載1.點擊文件菜單下的"改變工作目錄"來設定當前工作目錄。2.輸入以下命令把文件中的數據讀入：

algae<-read.table('Analysis.txt',

s=c('season',

'size',

'speed',

'mxPH','mnO2',

'Cl',

'NO3',

'NH4',

'oPO4','PO4','Chla',

'a1',

'a2',

'a3',

'a4',

'a5',

'a6',

'a7'),

na.strings=c('XXXXXXX'))3.點擊文件菜單下的“保存工作空間”，輸入文件名，退出，下次打開R后可通過拖拽的方式直接打開。2/3/202352數據摘要鑒于沒有該問題領域足夠的信息，首先了解一些數據的統計特性是一種較好的方式，它方便我們更好地理解問題。獲取數據統計特性的方法之一是獲取其描述性的統計摘要。命令如下：

summary(algae)對于名義變量，它給出了每個可能取值的頻數。對于數值變量，它提供了均值、中位數、四分位數和極值等一系列統計信息。NA's表示缺失值的個數。通過觀察這些值，我們可以了解到數據分布的偏度和分散情況。2/3/202363數據可視化(1)1.繪制變量mxPH的直方圖的兩種方式：

hist(algae$mxPH)

hist(algae$mxPH,prob=T) 區別在于前者給出的是頻數，后者是區間的概率。2.繪制mxPH的Q-Q圖： library(car) qqPlot(algae$mxPH,main='NormalQQplotofmaximumpH')

Q-Q圖繪制變量值和正態分布的理論分位數的散點圖。同時，它給出正態分布的95%置信區間的帶狀圖(虛線)。main為設置圖形的標題。2/3/202373數據可視化(2)3.繪制變量oPO4的箱圖： boxplot(algae$oPO4,ylab="oPO4") rug(algae$oPO4,side=4) abline(h=mean(algae$oPO4,na.rm=T),lty=2) ylab為設置y軸標題；

rug函數繪制變量的實際值，side=4表示繪制在圖的右側(1在下方，2在左側，3在上方)；

abline函數繪制水平線，mean表示均值，na.rm=T指計算時不考慮NA值，lty=2設置線型為虛線。2/3/202384數據清理 數據缺失的情形在實際問題中非常普遍。處理含有缺失值的數據時常用的幾種策略：將含有缺失值的案例剔除。用中心趨勢值來填補缺失值。根據變量之間的相關關系填補缺失值。根據案例之間的相似性填補缺失值。使用能夠處理缺失值數據的工具(見下一節)。2/3/202394.1剔除缺失值(1)1.適用范圍：含缺失值的記錄在整個數據集中比例很小時。2.檢查含缺失值的記錄：

algae[!complete.cases(algae),

]3.剔除所有含缺失值的記錄： algae<-na.omit(algae)4.找出每個記錄中缺失值的個數： apply(algae,1,function(x)sum(is.na(x)))

函數apply()是元函數，可在某些條件下對對象應用其他函數。第二個參數“1”表示第一個參數algae中的對象的第一個維度，即行數據。第三個參數是臨時函數，功能是計算對象x中NA的數量。(注：R中有TRUE=1,FALSE=0)2/3/2023104.1剔除缺失值(2)5.可根據4編寫一個找出algae中含有給定數目缺失值的行。以下函數的功能是找出缺失值個數大于列數20%的行： library(DMwR) manyNAs(algae,0.2) 第二個參數如不指定，默認為0.2，下面的命令與上一條等價： manyNAs(algae)6.我們利用上面的函數來剔除缺失值較多的記錄： algae<-algae[-manyNAs(algae),]

這里第二個參數的默認值為0.2。2/3/2023114.2用中心趨勢值填補缺失值 代表中心趨勢的值反映了變量分布的最常見值，這種方法也最自然、簡便和快捷。對于接近正態的分布來說，均值是最佳選擇；對偏態分布或有離群值的分布而言，中位數通常是更好的代表數據中心趨勢的指標；對于名義變量，通常采用眾數。可用以下函數完成填補所有缺失值： data(algae)

algae<-algae[-manyNAs(algae),

] algae<-centralImputation(algae) 上述方法特別適用于大數據集，但是這種方法可能導致較大的數據偏差，影響后期的數據分析。但使用復雜的無偏方法尋找最佳數據填補對大型數據集可能也不適用。2/3/2023124.3通過變量的相關性填補缺失值(1)

1.用以下命令獲取變量之間的相關矩陣： data(algae) symnum(cor(algae[,4:18],use="complete.obs")) 其中，函數cor()產生相關值矩陣(忽略前3個名義變量)，use參數指計算相關值時忽略含有NA的記錄。2.結果顯示，有兩個相關性較大的值：NH4和NO3之間，PO4和oPO4之間。 前者相關性不是特別明顯(0.6~0.8)，考慮到只樣本62和樣本199含有過多的缺失值，若剔除它們，樣本中NH4和NO3就不存在缺失值了。后者相關值很高(大于0.9)*，可用變量的相關性填補缺失值。*根據領域專家的解釋，總的磷酸鹽值包含正磷酸鹽值。2/3/2023134.3通過變量的相關性填補缺失值(2)3.尋找PO4和oPO4之間的線性關系的方法： algae<-algae[-manyNAs(algae),] lm(formula=PO4~oPO4,data=algae) 我們得到線性模型：PO4=42.897+1.293*oPO4.4.剔除樣本62和樣本199后，僅樣本28在PO4上有缺失值，我們用上面的線性關系來填補：

algae[28,"PO4"]<-42.897+1.293*algae[28,"oPO4"]

查看填補的記錄： algae[28,]2/3/2023144.4通過案例的相關性填補缺失值1.度量相似性時，最常用的指標是歐式距離。我們可通過使用這種度量來尋找與任何含有缺失值的案例最相似的10個水樣，并用它們填補缺失值。 方法一：簡單計算這10個最近的案例的中位數并用中位數填補缺失值；若缺失值是名義變量則采用眾數。 方法二：采用這些最相似數據的加權均值。這里用高斯核函數從距離獲得權重。命令如下：

clean.algae<-knnImputation(algae,k=10)2.這種方法看起來更合理，但也可能存在不相關的變量扭曲相似性，甚至造成大型數據集的計算特別復雜等問題。因此，填補缺失值時，大多應根據分析領域的知識來確定。2/3/2023155獲取預測模型用于預測海藻的兩種模型：多元線性回歸模型和回歸樹模型。線性回歸不能使用有缺失值的數據集，而回歸樹模型可以很自然地處理帶缺失值的數據。多元線性回歸模型是最常用的統計數據分析方法，該方法給出了一個有關目標變量與一組解釋變量關系的線性函數。由于多元線性回歸模型中沒有處理缺失值的方法，因此，我們可以做如下的數據預處理：data(algae)algae<-algae[-manyNAs(algae),]clean.algae<-knnImputation(algae,k=10)這里還是先移除缺失值較多的記錄，然后根據訓練集數據個案的相似性來填補缺失值。2/3/2023165.1線性回歸模型(1) 建立用于預測海藻頻率的線性回歸模型：

lm.a1<-lm(a1~.,data=clean.algae[,1:12]) 函數lm()建立一個線性回歸模型，其中，第一個參數給出了模型的函數形式。這個函數的形式是用數據中的其他所有變量來預測變量a1，第一個參數中的點“.”代表數據框中的所有除a1外的變量。如需要用變量mxPH和NH4來預測變量a1，就要定義模型為"a1~mxPH+NH4"。參數data是用來設定建模所用的數據集。2/3/2023175.1線性回歸模型(2) 通過下面的代碼，我們可以獲取更多線性模型的信息：

summary(lm.a1) 首先，給出數據擬合的殘差(residuals)。殘差應該是均值為0且為正態分布的。 其次，對于每個多元線性回歸方程的系數(變量)，給出它的估計值和標準誤差，并使用t檢驗來驗證系數為0的假設檢驗。 再者，給出模型與數據的吻合度，即模型所能解釋的數據變差的比例。R-squared越接近于1說明模型擬合得越好，越小則代表模型擬合得越差。 最后，還可以檢驗任何解釋變量與目標變量的依賴關系。

2/3/2023185.2回歸樹模型(1)因模型解釋的方差比例太低，才約32%，故實際驗證結果表明：對海藻案例應用線性模型是不合適的。用線性思維去考慮非線性問題，得不到理想的結果。我們考慮使用回歸樹預測。建立回歸樹： library(rpart) data(algae) algae<-algae[-manyNAs(algae),] rt.a1<-rpart(a1~.,data=algae[,1:12])函數的形式是用數據中其他所有變量來預測a1，data是用來設定建模所用的數據集。2/3/2023195.2回歸樹模型(2) 回歸樹rt.a1的圖形表示的兩種方法： plot(rt.a1) text(rt.a1) 或

prettyTree(rt.a1) 建立回歸樹通常分兩步。最初，生成一棵較大的樹，然后通過統計估計刪除底部的一些結點來對樹進行修剪。這樣是為了防止過度擬合。2/3/2023206模型評價和選擇 使用已有的訓練數據獲得模型的性能指標是不可靠的，因為這些計算是有偏的。實際上，有的模型可以很容易獲得訓練數據的零誤差預測。然而，這一優秀性能很難推廣到目標變量值未知的新樣本上。這種現象我們通常稱為過度擬合訓練數據。我們需要一個模型，使它在未知數據上有可靠的預測性能。 k折交叉驗證是獲得模型性能可靠估計的一種常用方法，它適用于小數據集。2/3/202321k折交叉驗證方法

K折交叉驗證：初始采樣分割成K個子樣本，一個單獨的子樣本被保留作為驗證模型的數據，其他K-1個樣本用來訓練。交叉驗證重復K次，每個子樣本驗證一次，平均K次的結果或者使用其它結合方式，最終得到一個單一估測。常見的選擇是k=10。有時還會重復進行多次K折交叉驗證以獲得更加可靠的估計。 總之，在做一項預測任務時要做出以下決策： 1)為預測任務選擇模型； 2)選擇比較模型性能的評估指標； 3)選擇獲取評估指標的可靠估計的實驗方法。2/3/202322進行模型比較的函數 在R中，有函數experimentalComparison()可用來進行模型之間的選擇和比較。它有三個參數：1)用于比較的數據集；2)需要比較的可選模型；3)實驗過程中的系數。此函數適用于任何模型和任何數據，從這個意義上說，它是一個泛型函數。 使用者提供一組實現待比較的模型的函數，其中每一個函數應該對訓練集和測試集實現一個完整的“訓練+測試+評估”周期。在評估過程的每一次迭代中，調用這些函數。這些函數返回一個向量，其元素為交叉驗證中用戶需要的性能評估指標量。2/3/202323兩個目標模型的函數cv.rpart<-function(form,train,test,...){ m<-rpartXse(form,train,...) p<-predict(m,test) mse<-mean((p-resp(form,test))^2) c(nmse=mse/mean((mean(resp(form,train))-resp(form,test))^2))}cv.lm<-function(form,train,test,...){ m<-lm(form,train,...) p<-predict(m,test) p<-ifelse(p<0,0,p) mse<-mean((p-resp(form,test))^2) c(nmse=mse/mean((mean(resp(form,train))-resp(form,test))^2))}2/3/202324函數的說明11.我們用NMSE(標準化后的平均絕對誤差)作為回歸樹模型和線性回歸模型的性能評估指標。2.這些函數的前三個參數是公式、訓練數據和測試數據。3.特殊參數"..."可用在任意的R函數中，它允許一個特定函數具有可變的參數，它用來給實際模型傳遞所需要的額外參數。4.函數resp()用于根據公式獲得數據集的目標變量值。2/3/202325模型的交叉驗證比較

定義好用于模型學習和測試的函數后，我們可按下列代碼進行模型的交叉驗證比較：

res<-experimentalComparison(

c(dataset(a1~.,clean.algae[,1:12],'a1')),

c(variants('cv.lm'),

variants('cv.rpart',se=c(0,0.5,1))),

cvSettings(3,10,1234)

)

參數一的形式為dataset(<formula>,<dataframe>,<label>);

參數二包含可選的模型方法并通過variants()指定，該函數第一個參數為目標模型的函數名，第二個作為可選參數；

參數三設定交叉驗證試驗的參數(repetition,fold,seed)。2/3/202326查看比較結果1.查看比較結果的摘要：

summary(res) 從結果中可知，apart.v1有最優的平均NMSE值。2.驗證結果也可轉化成可視化圖形：

plot(res)3.可通過以下代碼查看模型所對應的參數： getVariant("cv.rpart.v1",res)2/3/202327多個預測任務同時進行DSs

<-

sapply(

names(clean.algae)[12:18], function(x,names.attrs){ f<-as.formula(paste(x,"~.")) dataset(f,clean.algae[,c(names.attrs,x)],x) }, names(clean.algae)[1:11])res.all<-experimentalComparison(

DSs, c(variants('cv.lm'),

variants('cv.rpart',se=c(0,0.5,1)) ), cvSettings(5,10,1234))2/3/202328函數的說明21.該代碼首先創建用于比較7個預測任務的數據集向量。對每一個預測問題需構建一個公式，該公式由一個字符串構成，它是數據集中相應的需要預測的目標變量和符號"~."連接而成的。然后，該字符串通過函數as.formula()轉換為一個R公式。2.這次采用重復5次10折交叉驗證以提高統計結果的顯著性。3.本條指令運行的時間稍長(處理器：Intel(R)Core(TM)i5-3210MCPU@2.50GHz，2501Mhz，2個內核，4個邏輯處理器，約需1分鐘)。2/3/202329模型對不同海藻的結果1.所有海藻交叉驗證結果的可視化：

plot(res.all) 圖中顯示有幾個結果很差，即其NMSE值明顯大于1。這意味著測試結果比簡單采用目標變量的均值這一基準模型還要差！2.查看每個問題對應的最優模型的代碼：

bestScores(res.all)

其結果說明，只有海藻1的預測結果尚可。3.考慮用組合法進行模型構建。通過產生大量可選模型并把其進行組合，這樣得到的模型可以克服單個模型的局限性。2/3/202330隨機森林1.隨機森林作為組合模型的代表，它由大量的樹模型(回歸樹或分類樹)構成。每個樹是完全生長的(沒有事后剪枝)，在樹生長的每一步，最好的結點分割方法將從變量集合的一個隨機子集中選取。回歸任務的預測采用組合中預測結果的平均值。2.以下代碼是包含三個版本的隨機森林模型的交叉驗證，在組合中每個模型有不同數目的數。2/3/202331隨機森林的交叉驗證library(randomForest)cv.rf<-function(form,train,test,...){ m<-randomForest(form,train,...) p<-predict(m,test) mse<-mean((p-resp(form,test))^2) c(nmse=mse/mean((mean(resp(form,train))-resp(form,test))^2))}res.all<-experimentalComparison(

DSs, c(variants('cv.lm'),

variants('cv.rpart',se=c(0,0.5,1)),

variants('cv.rf',ntree=c(200,500,700))), cvSettings(5,10,1234))2/3/202332模型應用的結果1.采用下面的函數來證實組合方法的優勢：

bestScores(res.all) 對于某些問題，隨機森林給出很好的結果。但像海藻7，結果還不能令人滿意。2.還不知道這些最佳模型和其它模型之間的區別是否顯著。即，采用別的隨機數據我們能得到相似結果的可能性有多大？ 我們用Wilcoxon(威爾科克森)檢驗來判斷顯著性。上述結果表明，對海藻1、2、4和6，模型"cv.rf.v3"最好，檢驗代碼為：

compAnalysis(res.all,against='cv.rf.v3',

datasets=c('a1','a2','a4','a6'))2/3/202333顯著性分析1.結果中的“sig.X”列是我們需要的信息。沒有任何標識符則意味著相應的模型和"cv.rf.v3"模型之間有顯著差異的可能性低于95%。加號意味著相應模型的平均性能估計指標顯著高于模型"cv.rf.v3"。由于好的模型對應較低的NMSE值，所以該模型的性能比模型"cv.rf.v3"差。減號的含義相反。2.從結果可知，不同版本之間的隨機森林模型的差異在統計上通常不顯著。與其他模型相比，在大部分情況下，隨機森林具有顯著的優勢。3.參數against和datasets取不同的值，可對在其他海藻上有最優性能的模型進行類似分析。2/3/2023347預測海藻頻率 本案例的目的：預測140個水樣的7個海藻的頻率值。 我們已經采用交叉驗證的過程給出了最佳的預測模型，下面應該應用所有可得的訓練數據來構建模型，并將得到的模型應用到測試數據集。為了避免回歸樹采用它自身的缺失值的處理方法，我們采用了k近鄰法填補數據框clean.algae的NA值。隨機森林本身沒有處理缺失值的方法，我們把數據框clean.algae作為它的訓練集數據。1.為每種藻類選擇最優的預測模型：

bestModelsNames<-sapply(bestScores(res.all),

function(x)x['nmse','system'])

learners<-c(rf='randomForest',rpart='rpartXse')

funcs<-learners[sapply(strsplit(bestModelsNames,'\\.'),

function(x)x[2])]2/3/202335parSetts<-lapply(bestModelsNames, function(x)getVariant(x,res.all)@pars)