我正在尝试调整线性模型上的对比度编码，我想知道因子的每个水平是否与总平均值显着不同。
假设因子具有水平"A"、"B"和"C"。最常见的对照-治疗对比显然将"A"设置为参考水平，并将"B"和"C"与此进行比较。这不是我想要的，因为"A"级不会出现在模型摘要中。
偏差编码似乎也没有给我我想要的东西，因为它将级别"C"的对比度矩阵设置为[-1,-1,-1]，现在这个级别没有显示在模型摘要中。
set.seed(1)
y <- rnorm(6, 0, 1)
x <- factor(rep(LETTERS[1:3], each = 2))
fit <- lm(y ~ x, contrasts = list(x = contr.sum))
summary(fit)
此外，报告的级别名称已从"A"、"B"更改为"1"和"2"。
Call:
lm(formula = y ~ x, contrasts = list(x = contr.sum))
Residuals:
1 2 3 4 5 6
-0.405 0.405 -1.215 1.215 0.575 -0.575
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.02902 0.46809 -0.062 0.954
x1 -0.19239 0.66198 -0.291 0.790
x2 0.40885 0.66198 0.618 0.581
Residual standard error: 1.147 on 3 degrees of freedom
Multiple R-squared: 0.1129, Adjusted R-squared: -0.4785
F-statistic: 0.1909 on 2 and 3 DF, p-value: 0.8355
我错过了什么吗？我是否应该添加一个等于总平均值的虚拟变量，以便我可以将其用作参考水平？
我去年看到一个类似的问题(但可能要求更高)，但没有解决方案(尚未)：分类变量上所有系数的"平均值差异"模型; 得到"对比编码"来做到这一点？我已经悬赏了它，希望引起注意。
这里接受的答案有效，但作者没有提供解释。我已经在统计数据SE上询问了它，并增加了一个赏金：https://stats.stackexchange.com/questions/600798/understanding-the-process-of-tweaking-contrasts-in-linear-model-fitting-to-show
这个答案告诉你如何获得下面的系数表：
# Estimate Std. Error t value Pr(>|t|)
#(Intercept) -0.02901982 0.4680937 -0.06199574 0.9544655
#A -0.19238543 0.6619845 -0.29061922 0.7902750
#B 0.40884591 0.6619845 0.61760645 0.5805485
#C -0.21646049 0.6619845 -0.32698723 0.7651640
太神奇了，不是吗？它模仿你从summary(fit)中看到的东西，即
# Estimate Std. Error t value Pr(>|t|)
#(Intercept) -0.02901982 0.4680937 -0.06199574 0.9544655
#x1 -0.19238543 0.6619845 -0.29061922 0.7902750
#x2 0.40884591 0.6619845 0.61760645 0.5805485
但现在我们显示了所有因子水平。
为什么lm汇总不显示所有因子水平？
2016 年，我回答了这个堆栈溢出问题："lm"摘要不显示所有因子水平，从那时起，它已成为标记类似主题重复问题的目标。
回顾一下，基本思想是，为了获得最小二乘拟合的全秩设计矩阵，我们必须对因子变量应用对比。假设因子有 N 个水平，那么无论我们选择哪种类型的对比(有关列表，请参见?contrasts)，它都会将原始 N 个虚拟变量减少为一组新的N- 1个变量。因此，只有 N- 1个系数与N水平因子相关联。
但是，我们可以使用对比矩阵将 N- 1系数转换回原始N系数。变换使我们能够获得所有因子水平的系数表。我现在将演示如何根据 OP 的可重现示例执行此操作：
set.seed(1)
y <- rnorm(6, 0, 1)
x <- factor(rep(LETTERS[1:3], each = 2))
fit <- lm(y ~ x, contrasts = list(x = contr.sum))
在此示例中，总和到零的对比度应用于因子x。要了解有关如何控制模型拟合的对比度的更多信息，请参阅如何使用 R 在回归分析中为变量设置对比度？中的回答。
R 代码演练
对于受和零对比影响的 N 个水平的因子变量，我们可以使用以下函数来获取 N x (N - 1) 转换矩阵，该矩阵将lm估计的(N - 1)系数映射回所有水平的N系数。
ContrSumMat <- function (fctr, sparse = FALSE) {
if (!is.factor(fctr)) stop("'fctr' is not a factor variable!")
N <- nlevels(fctr)
Cmat <- contr.sum(N, sparse = sparse)
dimnames(Cmat) <- list(levels(fctr), seq_len(N - 1))
Cmat
}
对于示例 3 水平因子x，此矩阵为：
Cmat <- ContrSumMat(x)
# 1 2
#A 1 0
#B 0 1
#C -1 -1
拟合模型fit报告此因子的系数为 3 - 1 = 2。我们可以将它们提取为：
## coefficients After Contrasts
coef_ac <- coef(fit)[2:3]
# x1 x2
#-0.1923854 0.4088459
因此，特定于水平的系数为：
## coefficients Before Contrasts
coef_bc <- (Cmat %*% coef_ac)[, 1]
# A B C
#-0.1923854 0.4088459 -0.2164605
## note that they sum to zero as expected
sum(coef_bc)
#[1] 0
类似地，我们可以在对比后得到协方差矩阵：
var_ac <- vcov(fit)[2:3, 2:3]
# x1 x2
#x1 0.4382235 -0.2191118
#x2 -0.2191118 0.4382235
并将其转换为对比前的那个：
var_bc <- Cmat %*% var_ac %*% t(Cmat)
# A B C
#A 0.4382235 -0.2191118 -0.2191118
#B -0.2191118 0.4382235 -0.2191118
#C -0.2191118 -0.2191118 0.4382235
解释：
模型截距
coef(fit)[1]是大平均值。计算
coef_bc是每个水平的均值与总均值的偏差。var_bc的对角线条目给出了这些偏差的估计方差。
然后，我们可以继续计算这些系数的 t 统计量和 p 值，如下所示。
## standard error of point estimate `coef_bc`
std.error_bc <- sqrt(diag(var_bc))
# A B C
#0.6619845 0.6619845 0.6619845
## t-statistics (Null Hypothesis: coef_bc = 0)
t.stats_bc <- coef_bc / std.error_bc
# A B C
#-0.2906192 0.6176065 -0.3269872
## p-values of the t-statistics
p.value_bc <- 2 * pt(abs(t.stats_bc), df = fit$df.residual, lower.tail = FALSE)
# A B C
#0.7902750 0.5805485 0.7651640
## construct a coefficient table that mimics `coef(summary(fit))`
stats.tab_bc <- cbind("Estimate" = coef_bc,
"Std. Error" = std.error_bc,
"t value" = t.stats_bc,
"Pr(>|t|)" = p.value_bc)
# Estimate Std. Error t value Pr(>|t|)
#A -0.1923854 0.6619845 -0.2906192 0.7902750
#B 0.4088459 0.6619845 0.6176065 0.5805485
#C -0.2164605 0.6619845 -0.3269872 0.7651640
我们还可以通过包含总平均值的结果(即模型截距)来增强它。
## extract statistics of the intercept
intercept.stats <- coef(summary(fit))[1, , drop = FALSE]
# Estimate Std. Error t value Pr(>|t|)
#(Intercept) -0.02901982 0.4680937 -0.06199574 0.9544655
## augment the coefficient table
stats.tab <- rbind(intercept.stats, stats.tab_bc)
# Estimate Std. Error t value Pr(>|t|)
#(Intercept) -0.02901982 0.4680937 -0.06199574 0.9544655
#A -0.19238543 0.6619845 -0.29061922 0.7902750
#B 0.40884591 0.6619845 0.61760645 0.5805485
#C -0.21646049 0.6619845 -0.32698723 0.7651640
我们也可以用意义星星打印这个表格。
printCoefmat(stats.tab)
# Estimate Std. Error t value Pr(>|t|)
#(Intercept) -0.02902 0.46809 -0.0620 0.9545
#A -0.19239 0.66199 -0.2906 0.7903
#B 0.40885 0.66199 0.6176 0.5805
#C -0.21646 0.66199 -0.3270 0.7652
嗯？为什么没有星星？在这个例子中，所有的p值都非常大。如果 p 值较小，则星星将显示。这是一个令人信服的演示：
fake.tab <- stats.tab
fake.tab[, 4] <- fake.tab[, 4] / 100
printCoefmat(fake.tab)
# Estimate Std. Error t value Pr(>|t|)
#(Intercept) -0.02902 0.46809 -0.0620 0.009545 **
#A -0.19239 0.66199 -0.2906 0.007903 **
#B 0.40885 0.66199 0.6176 0.005805 **
#C -0.21646 0.66199 -0.3270 0.007652 **
#---
#Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
哦，这太美了。有关这些星星的含义，请参阅我的回答：方差分析表的交互R显著性代码？
闭幕辞
应该可以编写一个函数(甚至R包)来执行这样的表转换。但是，要使此类功能足够灵活，可能需要付出很大的努力才能处理：
所有类型的对比(这很容易做到);
复杂的模型项，例如因子与其他数值/因子变量之间的交互作用(这似乎真的很涉及！！
所以，我暂时到此为止。
>杂项答复我从 lm 的摘要中获得的模型分数是否仍然准确，即使它没有显示因子的所有级别？
是的，他们是。lm进行精确的最小二乘拟合。
此外，系数表的变换不会影响 R 平方、自由度、残差、拟合值、F 统计量、方差分析表等。
这只是一个初步的功能化版本@Zheyuan惊人的答案。正如@Zheyuan也提到的那样，这个答案有局限性，包括具有对比的复杂模型等。我不相信如果因子是有序的，该函数会起作用。
奶嘴数据：
set.seed(1)
test.df <- data.frame(y = rnorm(10, 0, 1),
x = factor(rep(LETTERS[1:5], each = 2)))
test.fit <- lm(y ~ x, contrasts = list(x=contr.sum), data= test.df)
返回关卡统计信息的函数：
# Output deviation coded factor matrix, used within `DevContrStats()`
ContrSumMat <- function (fctr, sparse = FALSE) {
if (!is.factor(fctr)) stop("'fctr' is not a factor variable!")
N <- nlevels(fctr)
Cmat <- contr.sum(N, sparse = sparse)
dimnames(Cmat) <- list(levels(fctr), seq_len(N - 1))
Cmat
}
## `dat` is the original data
## `fit.model` is the linear model created from the data
## `fctr_col` is the colum you'd like to get all factor levels from
DevContrStats <- function(dat, fit.model, fctr_col) {
fctr.mtx <- ContrSumMat(factor(dat[, fctr_col]))
## coefficients After Contrasts
coef_a_contr <- coef(fit.model)[-1]
## coefficients Before Contrasts
## coef_bc tells how each factor level deviates from the grand mean.
coef_b_contr <- (fctr.mtx %*% coef_a_contr)[, 1]
## Covariance matrix after contrasts:
var_a_contr <- vcov(fit.model)[-1,-1]
## Transform covariance matrix (after contrasts) to the one before contrasts:
## The diagonal of var_bc gives the estimated variance of factor-level deviations.
var_b_contr <- fctr.mtx %*% var_a_contr %*% t(fctr.mtx)
## standard error of point estimate `coef_bc`
std.err_b_contr <- sqrt(diag(var_b_contr))
## t-statistics (Null Hypothesis: coef_bc = 0)
t.stats_b_contr <- coef_b_contr / std.err_b_contr
## p-values of the t-statistics
p.value_b_contr <- 2 * pt(abs(t.stats_b_contr),
df = fit.model$df.residual,
lower.tail = FALSE)
## construct a coefficient table that mimics `coef(summary(fit))`
stats.tab_b_contr <- cbind("Estimate" = coef_b_contr,
"Std. Error" = std.err_b_contr,
"t value" = t.stats_b_contr,
"Pr(>|t|)" = p.value_b_contr)
## extract statistics of the intercept, which = grand mean
intercept.stats <- coef(summary(fit.model))[1, , drop = FALSE]
## augment the coefficient table with the intercept stats
stats.tab <- rbind(intercept.stats, stats.tab_b_contr)
## print stats table with significance stars
printCoefmat(stats.tab)
}
17:20:20> DevContrStats(test.df, test.fit, 'x')
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.164183 0.245068 0.6699 0.53258
A 0.448694 0.490137 0.9154 0.40195
B -0.028986 0.490137 -0.0591 0.95513
C 0.786629 0.490137 1.6049 0.16942
D -1.582153 0.490137 -3.2280 0.02326 *
E 0.375816 0.490137 0.7668 0.47785
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
做一系列 t 检验，每个检验将子集的均值与总均值进行比较怎么样？
lapply(levels(x), (i) t.test(y,y[x==i]))
对比度是模型参数的线性组合。计算它们需要一点线性代数。您可以按照其他答案中的说明手动完成。或者，您可以使用几个软件包之一为您进行数学计算。
我将使用对比度包。在此处阅读其小插曲。
library("contrast")
set.seed(1)
xlevels <- LETTERS[1:3]
y <- rnorm(6, 0, 1)
x <- factor(rep(xlevels, each = 2))
铌。您无需更改模型的参数化(通过在lm中使用contrasts参数进行设置)。您可以计算任何感兴趣的对比度与任何参数化，那么为什么不保留默认参数化。
使用默认参数化拟合线性回归。
fit <- lm(y ~ x)
默认参数化将 A 作为参考水平。
broom::tidy(fit)
#> # A tibble: 3 × 5
#> term estimate std.error statistic p.value
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 (Intercept) -0.221 0.811 -0.273 0.803
#> 2 xB 0.601 1.15 0.524 0.636
#> 3 xC -0.0241 1.15 -0.0210 0.985
依次将每个级别与三个级别的平均值进行比较。
for (xi in xlevels) {
print(
paste("Contrast", xi, "with the average")
)
print(
contrast(
fit,
list(x = xi),
list(x = xlevels),
type = "average"
)
)
}
#> [1] "Contrast A with the average"
#> lm model parameter contrast
#>
#> Contrast S.E. Lower Upper t df Pr(>|t|)
#> 1 -0.1923854 0.6619845 -2.299116 1.914345 -0.29 3 0.7903
#> [1] "Contrast B with the average"
#> lm model parameter contrast
#>
#> Contrast S.E. Lower Upper t df Pr(>|t|)
#> 1 0.4088459 0.6619845 -1.697884 2.515576 0.62 3 0.5805
#> [1] "Contrast C with the average"
#> lm model parameter contrast
#>
#> Contrast S.E. Lower Upper t df Pr(>|t|)
#> 1 -0.2164605 0.6619845 -2.323191 1.89027 -0.33 3 0.7652
创建于 2023-01-07 带有 reprex v2.0.2