UA MATH571A QE練習(xí) R語(yǔ)言 單因子試驗(yàn)的回歸分析

2015年5月的第六題是單因子試驗(yàn)，因?yàn)闅v年只有這一道，所以單獨(dú)做一下。

土壤中的硅主要以硅酸鹽礦物的形式存在，受成土母質(zhì)和成土過(guò)程的影響，不同土類和同一土壤不同層次間的硅含量差異很大。母質(zhì)中的硅通過(guò)風(fēng)化、淋溶、淀積而變遷。在山地生態(tài)系統(tǒng)的風(fēng)化作用研究中，人們普遍認(rèn)為溫度升高會(huì)加劇土壤中硅的流失。為了驗(yàn)證這個(gè)猜測(cè)是否成立，有人設(shè)計(jì)了一個(gè)試驗(yàn)，試驗(yàn)數(shù)據(jù)如下：

這是四個(gè)地點(diǎn)在不同溫度下土壤中的硅含量。為這個(gè)數(shù)據(jù)建立的統(tǒng)計(jì)模型是
$$y_{ij}={\mu }_i+{?}_{ij},i=1,2,3,j=1,2,3,4,{?}_{ij}～N(0,{\sigma }^2)$$

下面用回歸的方法估計(jì)模型分析數(shù)據(jù)，并檢驗(yàn)上述猜測(cè)是否成立。雖然數(shù)據(jù)非常簡(jiǎn)單，我們還是畫一下散點(diǎn)圖看看大致的趨勢(shì)。數(shù)據(jù)可以在https://statistics.arizona.edu/sites/default/files/uagc_page/may_15_data_sets.xlsx下載。

silicon.df = read.csv( file.choose() ) attach( silicon.df ) Y = conc t = temp plot( Y ~ t, pch=19 )

從這個(gè)散點(diǎn)圖可以看出，溫度上升硅含量下降的趨勢(shì)是在的，因?yàn)槭窃谒膫€(gè)地點(diǎn)做的試驗(yàn)，所以每個(gè)溫度對(duì)應(yīng)四個(gè)點(diǎn)，線性回歸有可能會(huì)欠擬合。做線性回歸分析時(shí)，溫度的數(shù)值是有意義的，回歸模型其實(shí)是
$$y={\beta }_0+{\beta }_{1}t+?$$

siliconSLR.lm = lm( Y ~ t ) > summary(siliconSLR.lm)Call: lm(formula = Y ~ t)Residuals:Min 1Q Median 3Q Max -32.899 -13.924 5.246 14.726 28.803 Coefficients:Estimate Std. Error t value Pr(>|t|) (Intercept) 168.094 13.867 12.122 2.66e-07 *** t -16.699 2.075 -8.049 1.12e-05 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 21.16 on 10 degrees of freedom(12 observations deleted due to missingness) Multiple R-squared: 0.8663, Adjusted R-squared: 0.8529 F-statistic: 64.78 on 1 and 10 DF, p-value: 1.117e-05

線性回歸的結(jié)果顯示，硅含量與土壤溫度的確有顯著負(fù)相關(guān)，并且這個(gè)一元線性回歸的解釋力是比較強(qiáng)的，R方達(dá)到了0.85以上。但如果我們畫一下殘差圖：

plot( resid(siliconSLR.lm)~t, pch=19 ); abline( h=0 )

就會(huì)發(fā)現(xiàn)殘差其實(shí)是非常大的，而且殘差散點(diǎn)顯然不是均勻分布在橫軸上下的，而是有一些pattern，我們可以試試這種pattern是不是溫度的二次項(xiàng)，為此我們可以建立一個(gè)二次回歸模型
$$y={\beta }_0+{\beta }_{1}t+{\beta }_{11}{t}^2+?$$

siliconQR.lm = lm( Y ~ t + I(t^2) ) plot( resid(siliconQR.lm)~t, pch=19 ); abline( h=0 )

在上面的代碼中，函數(shù)lm()的formula輸入的是Y ~ t + I(t^2)，這里的I()是一定要加的，因?yàn)閒ormula會(huì)把所有的符號(hào)默認(rèn)為解釋變量的變量名，如果直接帶入運(yùn)算符，比如Y ~ t + t^2，輸出的結(jié)果就會(huì)與Y ~ t 一樣。除了加I()外，還可以定義一個(gè)新的變量來(lái)表示平方項(xiàng)，也可以避免這個(gè)問(wèn)題，比如

t_square = t^2 siliconQR.lm = lm( Y ~ t + t_square )

再說(shuō)上面的殘差圖，現(xiàn)在殘差關(guān)于溫度就沒(méi)有明顯的pattern了，并且殘差的數(shù)值也明顯變小了，所以加入二次項(xiàng)是能提升擬合效果的。我們可以把二次回歸曲線畫在數(shù)據(jù)散點(diǎn)圖中看看擬合效果：

bQR = coef( siliconQR.lm ) plot( Y ~ t, pch=19, xlim=c(3,10), ylim=c(0,150) ); par( new=T ) curve( bQR[1] + x*(bQR[2] + x*bQR[3]), xlim=c(3,10), ylim=c(0,150),ylab='', xlab='' )

上面代碼中，curve()里面的x不是我們定義的變量，它用在bQR[1] + x*(bQR[2] + x*bQR[3])中表示這是一個(gè)函數(shù)，注意在curve()中如果要自己寫表達(dá)式就只能用x表示curve的自變量。這個(gè)擬合看起來(lái)就比較完美了，但這個(gè)模型的缺點(diǎn)是不能得出溫度越高、硅含量越低這種簡(jiǎn)單明了的結(jié)論。答案的語(yǔ)言非常優(yōu)美，摘抄下來(lái)學(xué)習(xí)一下：
The ever-present danger with a quadratic fit is evident here: the very good fit also comes
with the unlikely implication that the mean response turns back up before we reach the
highest observed temperature. (Quick reflection suggests that this is hard to explain: it is
reasonable for the soil to lose silicon as temperature rises, but then how could it regain the
silicon as the temperature starts to rise even higher?)

于是我們排除掉二次回歸。但線性回歸發(fā)現(xiàn)的問(wèn)題還沒(méi)有解決，模型是有非線性pattern的，最常規(guī)的處理非線性pattern的方法是Box-Cox變換：
$$y^{\lambda }={\beta }_0+{\beta }_{1}t+{\beta }_{11}{t}^2+?$$

require( MASS ) siliconSLR.bc = boxcox( siliconSLR.lm,lambda=seq(-1, 1, 0.1), interp=F )

如果上面的代碼報(bào)錯(cuò)，可以另外定義一個(gè)變量，把a(bǔ)s.numeric(t)賦值給它，重新做一個(gè)線性回歸，然后把模型對(duì)象輸入給boxcox()。上圖的結(jié)果說(shuō)明$\lambda$取0.2-0.7時(shí)Box-Cox在95%置信水平下的結(jié)果是相同的，最優(yōu)的$\lambda$為0.4。下面我們估計(jì)這個(gè)模型：
$$y^{0.4}={\beta }_0+{\beta }_{1}t+{\beta }_{11}{t}^2+?$$

Y_boxcox = Y^0.4 siliconBC.lm = lm( Y_boxcox~t ) bQR.bc = coef( siliconBC.lm ) plot( Y ~ t, pch=19, xlim=c(3,10), ylim=c(0,150) ); par( new=T ) curve( (bQR.bc[1] + x*bQR.bc[2])^2.5, xlim=c(3,10), ylim=c(0,150),ylab='', xlab='' )> summary(siliconBC.lm)Call: lm(formula = Y_boxcox ~ t)Residuals:Min 1Q Median 3Q Max -0.73947 -0.43646 0.06689 0.30905 1.12419 Coefficients:Estimate Std. Error t value Pr(>|t|) (Intercept) 8.94390 0.37660 23.75 3.98e-10 *** t -0.69117 0.05635 -12.27 2.38e-07 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 0.5747 on 10 degrees of freedom Multiple R-squared: 0.9377, Adjusted R-squared: 0.9314 F-statistic: 150.5 on 1 and 10 DF, p-value: 2.377e-07

現(xiàn)在的模型擬合曲線的形狀和二次回歸差不多，但看起來(lái)要差一點(diǎn)，它的優(yōu)點(diǎn)在于擬合曲線是單調(diào)的，結(jié)合回歸的結(jié)果可以判斷溫度越高，硅含量越低，并且這個(gè)關(guān)系是顯著的。然而，畫出這個(gè)模型的殘差圖

resids = Y - (bQR.bc[1] + t*bQR.bc[2])^2.5 plot( resids~t, pch=19 ); abline( h=0 )

會(huì)發(fā)現(xiàn)這種表示凸性的pattern依然存在，并且似乎還多一個(gè)異方差的問(wèn)題。這個(gè)就簡(jiǎn)直了，為了同時(shí)解決這兩個(gè)問(wèn)題，我們對(duì)Box-Cox變換后的模型加上二次項(xiàng)做WLS，權(quán)重的選擇答案解釋的是可以選擇溫度對(duì)應(yīng)的數(shù)據(jù)的組內(nèi)方差，

s2 = by(data=silicon.df$conc, INDICES=factor(silicon.df$temp), FUN=var) w = rep( 1/s2, each=4 ) siliconWBC.lm = lm( Y_boxcox~t+I(t^2), weight=w ) plot( resid(siliconWBC.lm)~t, pch=19 ); abline( h=0 )

然而可怕的事情發(fā)生了，殘差圖顯示異方差的問(wèn)題根本沒(méi)有得到解決！

我們?cè)賮?lái)回想一下整個(gè)過(guò)程，首先我們嘗試了線性回歸，但線性回歸殘差比較大，并且殘差圖存在非線性的pattern，為了解決非線性pattern，我們先是把模型修改為二次回歸，但二次回歸擬合曲線不單調(diào)，很難做解釋；于是我們又試圖用Box-Cox變換解決非線性pattern，但Box-Cox變換不但沒(méi)有解決非線性pattern，還引入了異方差，為了解決這兩個(gè)問(wèn)題，我們最后又用了WLS估計(jì)Box-Cox變換后加入二次項(xiàng)的模型，殘差圖顯示非線性pattern解決了，但異方差問(wèn)題仍然存在。下圖可以看出擬合效果已經(jīng)非常完美了，擬合曲線也是單調(diào)的，雖然有異方差的pattern，但整體殘差都是比較小的，在試驗(yàn)溫度范圍內(nèi)下結(jié)論問(wèn)題不大了。

bQR.bc = coef( siliconWBC.lm ) plot( Y ~ t, pch=19, xlim=c(3,10), ylim=c(0,150) ); par( new=T ) curve( (bQR.bc[1] + x*(bQR.bc[2]+x*bQR.bc[3]))^2.5, xlim=c(3,10), ylim=c(0,150),ylab='', xlab='' )

如果不用回歸方法，改用ANOVA分析試驗(yàn)結(jié)果，可以參考UA MATH571B 試驗(yàn)設(shè)計(jì)III 單因素試驗(yàn)設(shè)計(jì)3中的Tukey檢驗(yàn)、Fisher Least Significant Difference方法、Dunnett方法。

總結(jié)

以上是生活随笔為你收集整理的UA MATH571A QE练习 R语言 单因子试验的回归分析的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。