为什么 R 中的用户定义对比度需要作为权重的逆矩阵提供?

问题描述 投票:0回答:1

除了默认对比代码(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() 函数是一个鬼鬼祟祟的小混蛋,如上所述。要应用对比权重,你需要给它你的权重矩阵的逆矩阵。

都没有解释为什么。此外,计算逆矩阵会弄乱对比度系数,因此它们变得难以解释,因为它们不再是标准单位,所以我想知道为什么有必要。

r categorical-data anova
1个回答
0
投票

让我试试看一个简单的例子会发生什么。

我们从数据开始:

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://stats.stackexchange.com/questions/535123/custom-contrasts-in-r-why-should-i-take-the-generalized-inverse-of-transpose-of/612043#612043

一个很好的参考是:

https://cran.r-project.org/web/packages/codingMatrices/vignettes/codingMatrices.pdf

© www.soinside.com 2019 - 2024. All rights reserved.