我在设置先验对比方面遇到了麻烦,并希望寻求一些帮助.以下代码应给出与因子水平"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级,对吧？另外估计看起来很奇怪......
我知道使用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
我真的很感激对此的一些解释.
谢谢,Jeremias
我认为问题在于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 ***
因此对于组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