先驗對比R中的lm（）

Question

我在設置先驗對比時遇到了問題，並且希望尋求一些幫助。下面的代碼應該給出因子級別「d」的兩個正交對比。

Response <- c(1,3,2,2,2,2,2,2,4,6,5,5,5,5,5,5,4,6,5,5,5,5,5,5) 
A <- factor(c(rep("c",8),rep("d",8),rep("h",8))) 
contrasts(A) <- cbind("d vs h"=c(0,1,-1),"d vs c"=c(-1,1,0)) 
summary.lm(aov(Response~A)) 

我得到的是：

Call: 
aov(formula = Response ~ A) 

Residuals: 
    Min   1Q  Median   3Q  Max 
-1.000e+00 -3.136e-16 -8.281e-18 -8.281e-18 1.000e+00 

Coefficients: 
     Estimate Std. Error t value Pr(>|t|)  
(Intercept) 4.0000  0.1091 36.661 < 2e-16 *** 
Ad vs h  -1.0000  0.1543 -6.481 2.02e-06 *** 
Ad vs c  2.0000  0.1543 12.961 1.74e-11 *** 
--- 
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.5345 on 21 degrees of freedom 
Multiple R-squared: 0.8889,  Adjusted R-squared: 0.8783 
F-statistic: 84 on 2 and 21 DF, p-value: 9.56e-11 

但我希望（攔截）的估計是5.00，因爲它應該是相等的水平d，對不對？另外其他估計看起來很奇怪...

我知道我可以更容易地使用更高版本（A，ref =「d」）（正確顯示它們的位置）獲得正確的值，但我有興趣學習正確的表述測試自己的假設。如果我運行與如下因素代碼（從網站）一個類似的例子，它按預期工作：

irrigation<-factor(c(rep("Control",10),rep("Irrigated 10mm",10),rep("Irrigated20mm",10))) 
biomass<-1:30 
contrastmatrix<-cbind("10 vs 20"=c(0,1,-1),"c vs 10"=c(-1,1,0)) 
contrasts(irrigation)<-contrastmatrix 
summary.lm(aov(biomass~irrigation)) 


Call: 
aov(formula = biomass ~ irrigation) 

Residuals: 
     Min   1Q  Median   3Q  Max 
-4.500e+00 -2.500e+00 3.608e-16 2.500e+00 4.500e+00 

Coefficients: 
        Estimate Std. Error t value Pr(>|t|)  
(Intercept)   15.5000  0.5528 28.04 < 2e-16 *** 
irrigation10 vs 20 -10.0000  0.7817 -12.79 5.67e-13 *** 
irrigationc vs 10 10.0000  0.7817 12.79 5.67e-13 *** 
--- 
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 3.028 on 27 degrees of freedom 
Multiple R-squared: 0.8899,  Adjusted R-squared: 0.8817 
F-statistic: 109.1 on 2 and 27 DF, p-value: 1.162e-13 

我真的很感激這一些解釋。

感謝，赫雷米亞斯

Answer 1

我認爲這個問題是在contrasts理解（您可能?contrasts查看詳細）。讓我詳細解釋：

如果使用默認方式爲factorA，

A <- factor(c(rep("c",8),rep("d",8),rep("h",8))) 
> contrasts(A) 
    d h 
c 0 0 
d 1 0 
h 0 1 

因此模型lm給你

Mean(Response) = Intercept + beta_1 * I(d = 1) + beta_2 * I(h = 1) 

summary.lm(aov(Response~A)) 
Coefficients: 
    Estimate Std. Error t value Pr(>|t|)  
(Intercept) 2.000  0.189 10.6 7.1e-10 *** 
Ad    3.000  0.267 11.2 2.5e-10 *** 
Ah    3.000  0.267 11.2 2.5e-10 *** 

因此，對於組c，平均只是攔截2，對於組d，平均值爲2 + 3 = 5，對於組h也是如此。

如果你使用自己的對比：

contrasts(A) <- cbind("d vs h"=c(0,1,-1),"d vs c"=c(-1,1,0)) 

A 
[1] c c c c c c c c d d d d d d d d h h h h h h h h 
attr(,"contrasts") 
    d vs h d vs c 
c  0  -1 
d  1  1 
h  -1  0 

適合你原來的模式是

Mean(Response) = Intercept + beta_1 * (I(d = 1) - I(h = 1)) + beta_2 * (I(d = 1) - I(c = 1)) 
    = Intercept + (beta_1 + beta_2) * I(d = 1) - beta_2 * I(c = 1) - beta_1 * I(h = 1) 

Coefficients: 
Estimate Std. Error t value Pr(>|t|)  
(Intercept) 4.000  0.109 36.66 < 2e-16 *** 
Ad vs h  -1.000  0.154 -6.48 2.0e-06 *** 
Ad vs c  2.000  0.154 12.96 1.7e-11 *** 

因此，對於組c，均值爲4 - 2 = 2，爲組d，平均是4 - 1 + 2 = 5，對於組h，平均值是4 - (-1) = 5。

==========================

更新：

做你的對比度，最簡單的方法是設置基準（參考）等級爲d。

contrasts(A) <- contr.treatment(3, base = 2) 
Coefficients: 
      Estimate Std. Error t value Pr(>|t|)  
(Intercept) 5.00e+00 1.89e-01 26.5 < 2e-16 *** 
A1   -3.00e+00 2.67e-01 -11.2 2.5e-10 *** 
A3   -4.86e-17 2.67e-01  0.0  1  

如果你想用你的對比：

Response <- c(1,3,2,2,2,2,2,2,4,6,5,5,5,5,5,5,4,6,5,5,5,5,5,5) 
A <- factor(c(rep("c",8),rep("d",8),rep("h",8))) 
mat<- cbind(rep(1/3, 3), "d vs h"=c(0,1,-1),"d vs c"=c(-1,1,0)) 
mymat <- solve(t(mat)) 
my.contrast <- mymat[,2:3] 
contrasts(A) <- my.contrast 
summary.lm(aov(Response~A)) 

Coefficients: 
Estimate Std. Error t value Pr(>|t|)  
(Intercept) 4.00e+00 1.09e-01 36.7 < 2e-16 *** 
Ad vs h  7.69e-16 2.67e-01  0.0  1  
Ad vs c  3.00e+00 2.67e-01 11.2 2.5e-10 *** 

參考：http://www.ats.ucla.edu/stat/r/library/contrast_coding.htm