UA MATH571A QE練習 R語言 單因子試驗的回歸分析

2015年5月的第六題是單因子試驗，因為歷年只有這一道，所以單獨做一下。

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

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

下面用回歸的方法估計模型分析數據，并檢驗上述猜測是否成立。雖然數據非常簡單，我們還是畫一下散點圖看看大致的趨勢。數據可以在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 )

從這個散點圖可以看出，溫度上升硅含量下降的趨勢是在的，因為是在四個地點做的試驗，所以每個溫度對應四個點，線性回歸有可能會欠擬合。做線性回歸分析時，溫度的數值是有意義的，回歸模型其實是
$$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

線性回歸的結果顯示，硅含量與土壤溫度的確有顯著負相關，并且這個一元線性回歸的解釋力是比較強的，R方達到了0.85以上。但如果我們畫一下殘差圖：

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

就會發現殘差其實是非常大的，而且殘差散點顯然不是均勻分布在橫軸上下的，而是有一些pattern，我們可以試試這種pattern是不是溫度的二次項，為此我們可以建立一個二次回歸模型
$$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 )

在上面的代碼中，函數lm()的formula輸入的是Y ~ t + I(t^2)，這里的I()是一定要加的，因為formula會把所有的符號默認為解釋變量的變量名，如果直接帶入運算符，比如Y ~ t + t^2，輸出的結果就會與Y ~ t 一樣。除了加I()外，還可以定義一個新的變量來表示平方項，也可以避免這個問題，比如

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

再說上面的殘差圖，現在殘差關于溫度就沒有明顯的pattern了，并且殘差的數值也明顯變小了，所以加入二次項是能提升擬合效果的。我們可以把二次回歸曲線畫在數據散點圖中看看擬合效果：

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])中表示這是一個函數，注意在curve()中如果要自己寫表達式就只能用x表示curve的自變量。這個擬合看起來就比較完美了，但這個模型的缺點是不能得出溫度越高、硅含量越低這種簡單明了的結論。答案的語言非常優美，摘抄下來學習一下：
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?)

于是我們排除掉二次回歸。但線性回歸發現的問題還沒有解決，模型是有非線性pattern的，最常規的處理非線性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 )

如果上面的代碼報錯，可以另外定義一個變量，把as.numeric(t)賦值給它，重新做一個線性回歸，然后把模型對象輸入給boxcox()。上圖的結果說明$\lambda$取0.2-0.7時Box-Cox在95%置信水平下的結果是相同的，最優的$\lambda$為0.4。下面我們估計這個模型：
$$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

現在的模型擬合曲線的形狀和二次回歸差不多，但看起來要差一點，它的優點在于擬合曲線是單調的，結合回歸的結果可以判斷溫度越高，硅含量越低，并且這個關系是顯著的。然而，畫出這個模型的殘差圖

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

會發現這種表示凸性的pattern依然存在，并且似乎還多一個異方差的問題。這個就簡直了，為了同時解決這兩個問題，我們對Box-Cox變換后的模型加上二次項做WLS，權重的選擇答案解釋的是可以選擇溫度對應的數據的組內方差，

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 )

然而可怕的事情發生了，殘差圖顯示異方差的問題根本沒有得到解決！

我們再來回想一下整個過程，首先我們嘗試了線性回歸，但線性回歸殘差比較大，并且殘差圖存在非線性的pattern，為了解決非線性pattern，我們先是把模型修改為二次回歸，但二次回歸擬合曲線不單調，很難做解釋；于是我們又試圖用Box-Cox變換解決非線性pattern，但Box-Cox變換不但沒有解決非線性pattern，還引入了異方差，為了解決這兩個問題，我們最后又用了WLS估計Box-Cox變換后加入二次項的模型，殘差圖顯示非線性pattern解決了，但異方差問題仍然存在。下圖可以看出擬合效果已經非常完美了，擬合曲線也是單調的，雖然有異方差的pattern，但整體殘差都是比較小的，在試驗溫度范圍內下結論問題不大了。

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分析試驗結果，可以參考UA MATH571B 試驗設計III 單因素試驗設計3中的Tukey檢驗、Fisher Least Significant Difference方法、Dunnett方法。

總結

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

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。