|
但是,，事實(shí)上，他們是同一個(gè)模型的特例而已,。
這個(gè)部分包含一些復(fù)雜模型以及使用lm()構(gòu)造模型的過程以及在這個(gè)過程中經(jīng)常出現(xiàn)的問題的處理,。
在多元回歸里有的時(shí)候不像它看起來那么簡(jiǎn)單,，有時(shí)可以在多元回歸分析中納入變量的二次和高次冪，盡管這個(gè)看似是非線性關(guān)系的模型依舊屬于線性模型的范疇,，重點(diǎn)在于參數(shù)和預(yù)期的觀測(cè)值是線性關(guān)系,。
下面以之前使用的cystic fibrosis 數(shù)據(jù)集為例。
首先我們展示一下數(shù)據(jù)：
> cystfibr
age sex height weight bmp fev1 rv frc tlc pemax
1 7 0 109 13.1 68 32 258 183 137 95
2 7 1 112 12.9 65 19 449 245 134 85
…
24 23 0 175 51.1 71 33 224 131 113 95
25 23 0 179 71.5 95 52 225 127 101 195
另外,，我們展示一下之前在線性回歸中,，所給出的相關(guān)圖：
展示一下變量pemax和變量height之間的關(guān)系?？梢钥吹蕉邔儆诜蔷€性關(guān)系,。通過在模型公式中增加一個(gè)變量height的平方可對(duì)此進(jìn)行改造。
> summary(lm(pemax~height+I(height^2)))
Call:
lm(formula = pemax ~ height + I(height^2))
Residuals:
Min 1Q Median 3Q Max
-51.411 -14.932 -2.288 12.787 44.933
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 615.36248 240.95580 2.554 0.0181 *
height -8.08324 3.32052 -2.434 0.0235 *
I(height^2) 0.03064 0.01126 2.721 0.0125 *
—
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 24.18 on 22 degrees of freedom
Multiple R-squared: 0.5205, Adjusted R-squared: 0.4769
F-statistic: 11.94 on 2 and 22 DF, p-value: 0.0003081
#Tips：需要注意的是,，模型公式中的height^2需要用I()函數(shù)進(jìn)行保護(hù),。這個(gè)技術(shù)常常用來防止模型公式中的操作符被特殊解釋。這種解釋不作用于函數(shù)命令內(nèi)部,，I()是反身函數(shù),，原封不動(dòng)地返回自身的輸入?yún)?shù)。
使用predict()函數(shù)可以繪制帶預(yù)測(cè)值和置信帶的擬合曲線,。我們采用幾個(gè)有序的的身高點(diǎn)來做圖：
> pred.frame<-data.frame(height=seq(110,180,2))
> lm.pemax.hq<-lm(pemax~height+I(height^2))
> predict(lm.pemax.hq,interval=”pred”,newdata=pred.frame)
fit lwr upr
1 96.90026 37.94461 155.8559
2 94.33611 36.82985 151.8424
…
35 147.21294 93.51117 200.9147
36 152.98174 98.36718 207.5963
> pp<-predict(lm.pemax.hq,newdata=pred.frame,interval=”pred”)
> pc<-predict(lm.pemax.hq,newdata=pred.frame,interval=”conf”)
> plot(height,pemax,ylim=c(0,200))
> matlines(pred.frame$height,pp,lty=c(1,2,2),col=”black”)
> matlines(pred.frame$height,pc,lty=c(1,3,3),col=”black”)
#Tips：我們可以從圖里看出來,，身高較小時(shí)，擬合曲線隨著身高的增加而減少,，而后隨著身高的增加而快速增加,。這是為擬合數(shù)據(jù)選取了二次多項(xiàng)式形式的模型導(dǎo)致的。而身高和pemax變量之間的關(guān)系也可以很好地使用二次項(xiàng)來解釋,。
有時(shí)候,，數(shù)據(jù)會(huì)根據(jù)某個(gè)連續(xù)尺度的分段進(jìn)行分組，或者試驗(yàn)設(shè)計(jì)過程中令變量取幾個(gè)特定的x值的集合,。這兩種情況都跟比較線性回歸的結(jié)果和方差分析的結(jié)果相關(guān),。
對(duì)于同樣的數(shù)據(jù)，我們有兩種可供選擇的數(shù)據(jù)模型,。兩者都屬于線性模型的范疇,，且都能通過lm()函數(shù)擬合。線性回歸模型是單因素方差分析模型的子模型,，因?yàn)榍罢呖梢酝ㄟ^向后者的參數(shù)添加約束來獲得,。
可以通過檢驗(yàn)較大模型對(duì)應(yīng)的殘差方差來決定模型降價(jià)是否通過，這就是F檢驗(yàn),。
我們以trypsin（胰蛋白酶）濃度的例子里,，數(shù)據(jù)按照年齡進(jìn)行分組，數(shù)據(jù)以6個(gè)組對(duì)應(yīng)的均值和標(biāo)準(zhǔn)差的形式給出,。R沒法直接處理這種形式的數(shù)據(jù),，因此,，有必要?jiǎng)?chuàng)造一個(gè)均值和標(biāo)準(zhǔn)差都相同的偽數(shù)據(jù)，對(duì)應(yīng)的代碼如下：
方差分析的真是結(jié)果僅依賴于均值和標(biāo)準(zhǔn)差,，對(duì)生成偽數(shù)據(jù)的過程感興趣的可以看一下ISwR文件夾中的fake.trypsin.R文件,。
數(shù)據(jù)框fake.trypsin共包含3個(gè)變量，可以運(yùn)行下面代碼查看：
> summary(fake.trypsin)
trypsin grp grpf
Min. :-39.96 Min. :1.000 1: 32
1st Qu.:119.52 1st Qu.:2.000 2:137
Median :167.59 Median :2.000 3: 38
Mean :168.68 Mean :2.583 4: 44
3rd Qu.:213.98 3rd Qu.:3.000 5: 16
Max. :390.13 Max. :6.000 6: 4
#Tips：數(shù)據(jù)框既包括數(shù)值型向量grp,，也包括一個(gè)具有6個(gè)水平的因子型變量grpf,。
對(duì)生成的數(shù)據(jù)執(zhí)行單因素方差分析得到ANOVA表如下：
> anova(lm(trypsin~grpf))
Analysis of Variance Table
Response: trypsin
Df Sum Sq Mean Sq F value Pr(>F)
grpf 5 224103 44821 13.508 9.592e-12 ***
Residuals 265 879272 3318
—
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
如果用變量grp代替模型公式中的grpfin變量，將得到一個(gè)關(guān)于分組號(hào)的線性回歸模型,。得到的ANOVA表格如下：
> anova(lm(trypsin~grp))
Analysis of Variance Table
Response: trypsin
Df Sum Sq Mean Sq F value Pr(>F)
grp 1 206698 206698 62.009 8.451e-14 ***
Residuals 269 896677 3333
—
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Tips：殘差均方的變化不大,，這意味著這兩個(gè)模型對(duì)數(shù)據(jù)的刻畫效果基本上一樣。
如果想做一個(gè)正規(guī)的檢驗(yàn)來比較簡(jiǎn)單線性模型和各組具有獨(dú)立均值的模型的話,，可以直接運(yùn)行下面代碼：
> anova(lm(trypsin~grp+grpf))
Analysis of Variance Table
Response: trypsin
Df Sum Sq Mean Sq F value Pr(>F)
grp 1 206698 206698 62.2959 7.833e-14 ***
grpf 4 17405 4351 1.3114 0.2661
Residuals 265 879272 3318
—
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Tips：從上面結(jié)果可以看出,，grpf對(duì)模型的擬合沒有grp好。而之前他們對(duì)整個(gè)模型的解釋都是有意義的,，而grpf在這個(gè)模型里的p值突然無意義了,，說明了兩個(gè)變量之間有明顯的共線問題。
我們可以繪制帶擬合線和經(jīng)驗(yàn)均值的偽數(shù)據(jù)的散點(diǎn)圖：
> xbar.trypsin<-tapply(trypsin,grpf,mean)
> stripchart(trypsin~grp,method=”jitter”,jitter=0.1,vertical=T,pch=20)
> lines(1:6,xbar.trypsin,type=”b”,pch=4,cex=2,lty=2)
> abline(lm(trypsin~grp))
#Tips：這個(gè)圖的技巧本質(zhì)上同之前的圖相同,，所以這里不再對(duì)繪圖細(xì)節(jié)進(jìn)行詳細(xì)描述,。圖中有一個(gè)點(diǎn)表示胰島素密度是負(fù)的，這個(gè)數(shù)據(jù)是偽造的,，原始數(shù)據(jù)是沒法看到,。
多元回歸模型的一個(gè)基本假設(shè)就是模型中的各變量對(duì)響應(yīng)變量的影響具有疊加效應(yīng)。然而,，這并不是說線性模型模型沒法刻畫非疊加效應(yīng),。我們可以通過添加特殊交互項(xiàng)來指定一個(gè)變量受另一個(gè)變量水平變動(dòng)的影響程度。在R中的模型公式里,，交互項(xiàng)可以使用“：”來生成,，比如a：b。通常,，我們還會(huì)在模型中包含a和b這兩項(xiàng)，同時(shí),，R的模型里允許a*b或者a+b+a:b這種公式,，這兩個(gè)公式是等效的。
當(dāng)然在模型建立的過程中還有很多需要注意很多事項(xiàng),，我們這里就不一一列舉了,。
參考資料：
1.《R語言統(tǒng)計(jì)入門（第二版）》人民郵電出版社 Peter Dalgaard著
2.《R語言初學(xué)者指南》人民郵電出版社 Brian Dennis著
3.Vicky的小筆記本《blooming for you》by Vicky
|
