除了默认对比代码(contr.treatment / contr.sum / contr.helmert)之外,我还想在 R 中使用一些用户定义的对比。但是,我阅读的指南表明这些需要在 inverse 矩阵中提供。谁能解释为什么?
即这里的指导:https://stats.idre.ucla.edu/r/library/r-library-contrast-coding-systems-for-categorical-variables/ 状态:
最终的对比矩阵(或编码方案)原来是 mat 转置的逆矩阵。
同样,这个网站:https://rstudio-pubs-static.s3.amazonaws.com/65059_586f394d8eb84f84b1baaf56ffb6b47f.html 写道:
虽然 R 中有一些自动对比功能(到目前为止我一直在使用的功能),但您有时会发现自己想要运行其中未包含的比较。当发生这种情况时,您可以自己指定它们。但是,您需要小心,因为 contrasts() 函数是一个鬼鬼祟祟的小混蛋,如上所述。要应用对比权重,你需要给它你的权重矩阵的逆矩阵。
都没有解释为什么。此外,计算逆矩阵会弄乱对比度系数,因此它们变得难以解释,因为它们不再是标准单位,所以我想知道为什么有必要。
让我试试看一个简单的例子会发生什么。
我们从数据开始:
x = sample(letters[1:3], size = 300, replace = T)
y = rnorm(300)
head(x, 10)
[1] "c" "c" "c" "c" "a" "a" "a" "c" "b" "a"
假设您的自定义对比是:
contr_mat = matrix(c(1, -0.5, -0.5,
1, 0, -1),
nrow = 2, byrow = TRUE)
contr_mat
[,1] [,2] [,3]
[1,] 1 -0.5 -0.5
[2,] 1 0.0 -1.0
因此,例如,您正在比较 a 组中 y 的平均值与 c 组中 y 的平均值(第二个对比)。
让我们使用对比矩阵的伪逆构建对比编码矩阵。我们还添加了一个截距项:
C1 <- cbind(1, ginv(contr_mat))
round(C1, 3)
[,1] [,2] [,3]
[1,] 1 0.667 0
[2,] 1 -1.333 1
[3,] 1 0.667 -1
这个矩阵是方阵和满秩的,通过反转它我们得到了对比,有一个截距项:
C <- round(solve(C1), 2)
C
[,1] [,2] [,3]
[1,] 0.33 0.33 0.33
[2,] 1.00 -0.50 -0.50
[3,] 1.00 0.00 -1.00
因为 C 和 C1 是逆的,您可以将它们包含在支持您的比较的模型中:
y = Xu + e = XIu + e= X(C1C)u + e = (XC1)(Cu) + e
其中 X 是设计矩阵(对于此单元均值模型),u 是均值向量,I 是单位矩阵,e 是具有正态分布的随机分量的向量。 (Cu)代表你想要的反差
现在我们可以使用最小二乘方程来评估对比度 Cu,使用修改后的设计矩阵 (XC1)。所以对比矩阵的逆被用来评估实际的对比 Cu:
X <- model.matrix(~x + 0) # model matrix for the cell means model
X1 <- X %*% C1 # this is the modified model matrix with the inverse
现在我们用最小二乘法或者lm()求解。 lm() 将对比度矩阵的逆作为参数并添加截距项,然后通过最小二乘法评估对比度:
# least-squares equation:
solve ( t(X1) %*% X1 ) %*% t(X1) %*% y
[,1]
[1,] -0.06524626
[2,] 0.03627632
[3,] 0.04820965
# lm() with modified design matrix
summary(lm(y ~ X1 +0 ) )
Call:
lm(formula = y ~ X1 + 0)
Residuals:
Min 1Q Median 3Q Max
-3.1953 -0.7073 0.0250 0.7382 2.7251
Coefficients:
Estimate Std. Error t value Pr(>|t|)
X11 -0.06525 0.06049 -1.079 0.282
X12 0.03628 0.12596 0.288 0.774
X13 0.04821 0.14336 0.336 0.737
# lm() with your contrasts
summary( lm(y ~ x, contrasts = list(x = ginv(contr_mat))) )
Call:
lm(formula = y ~ x, contrasts = list(x = ginv(contr_mat)))
Residuals:
Min 1Q Median 3Q Max
-3.1953 -0.7073 0.0250 0.7382 2.7251
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.06525 0.06049 -1.079 0.282
x1 0.03628 0.12596 0.288 0.774
x2 0.04821 0.14336 0.336 0.737
最后两行包含您想要的比较:
means <- tapply(X = y, INDEX = factor(x), FUN = mean)
C %*% means
[,1]
[1,] -0.06524626
[2,] 0.03627632
[3,] 0.04820965
(means[1] + means[2] + means[3])/3
-0.06524626
means[1] - (means[2] + means[3])/2
0.03627632
means[1] - means[3]
0.04820965
我在这里对矩阵项进行了更好的格式化,做出了类似的回应:
一个很好的参考是:
https://cran.r-project.org/web/packages/codingMatrices/vignettes/codingMatrices.pdf