2

我在设置先验对比时遇到问题，想寻求帮助。以下代码应给出与因子水平“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

但我希望 (Intercept) 的估计值为 5.00，因为它应该等于 d 级，对吧？其他估计也看起来很奇怪......

我知道我可以使用 relevel(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

我真的很感激对此的一些解释。

谢谢， 杰里米亚斯

4

1 回答 1

6

我认为问题在于理解contrasts（您可以?contrasts详细了解）。让我详细解释一下：

如果您使用默认方式factor 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

因此模型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 ***

所以对于 groupc来说，平均值就是截距2，对于 groupd来说，平均值是2 + 3 = 5相同的 group 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

于 2012-11-20T15:35:06.397 回答