R中的lm()的先验对比

I'm having trouble with setting a priori contrasts and would like to ask for some help. The following code should give two orthogonal contrasts to the factor level "d".

我在设置先验对比方面有困难，想寻求一些帮助。下面的代码应该给因子级“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))

What I get is:

我得到的是:

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

But I expect the Estimate of (Intercept) to be 5.00, as it should be equal to the level d, right? Also the other estimates look strange...

但是我期望(截距)的估计值是5，因为它应该等于d，对吧?另一个估算看起来很奇怪……

I know I can get the correct values easier using relevel(A, ref="d") (where they are displayed correctly), but I am interested in learning the correct formulation to test own hypotheses. If I run a similar example with the folowing code (from a website), it works as expected:

我知道使用relevel(A, ref="d")可以更容易地得到正确的值(在这里它们被正确显示)，但我对学习正确的公式来测试自己的假设很感兴趣。如果我使用folowing代码运行一个类似的示例(来自一个网站)，它会像预期的那样工作:

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

I would really appreciate some explanation for this.

我很想知道为什么。

Thanks, Jeremias

耶利米亚,谢谢你

1 个解决方案

#1

I think the problem is in the understanding of contrasts (You may ?contrasts for detail). Let me explain in detail:

我认为问题在于对对比的理解(你可以用对比来描述细节)。让我详细解释一下:

If you use the default way for factor A,

如果你用A的默认方法，

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

thus the model lm gives you are

因此，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 ***

So for group c, the mean is just intercept 2, for group d , the mean is 2 + 3 = 5, same for group h.

c组的均值是2,d组的均值是2 + 3 = 5,h组也是一样。

What if you use your own contrast:

如果你用自己的对比:

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

The model you fit turns out to be

你适合的模型是

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 ***

So for group c, the mean is 4 - 2 = 2, for group d, the mean is 4 - 1 + 2 = 5, for group h, the mean is 4 - (-1) = 5.

c组的均值是4 - 2 = 2,d组的均值是4 -1 + 2 = 5,h组的均值是4 -(-1)= 5。

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

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

Update:

更新:

The easiest way to do your contrast is to set the base (reference) level to be d.

最简单的方法是将基准(引用)级别设置为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

If you want to use your contrast:

如果你想用你的对比:

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 ***

Reference: http://www.ats.ucla.edu/stat/r/library/contrast_coding.htm

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

#1