Ada*_*m_G 5 regression r linear-regression lm categorical-data
我正在尝试调整线性模型上的对比编码，我想知道因子的每个级别是否与总平均值显着不同。
\n让\xe2\x80\x99s 表示该因子具有级别“A”、“B”和“C”。最常见的对照治疗对比显然将“A”设置为参考水平，并将“B”和“C”与其进行比较。这不是我想要的，因为“A”级没有出现在模型摘要中。
\n偏差编码似乎也没有给我想要的东西，因为它将级别“C”的对比矩阵设置为[-1,-1,-1]，现在这个级别没有显示在模型摘要中。
set.seed(1)\ny <- rnorm(6, 0, 1)\nx <- factor(rep(LETTERS[1:3], each = 2))\nfit <- lm(y ~ x, contrasts = list(x = contr.sum))\nsummary(fit)\n此外，报告的级别名称已从“A”、“B”更改为“1”和“2”。
\nCall:\nlm(formula = y ~ x, contrasts = list(x = contr.sum))\n\nResiduals:\n 1 2 3 4 5 6 \n-0.405 0.405 -1.215 1.215 0.575 -0.575 \n\nCoefficients:\n Estimate Std. Error t value Pr(>|t|)\n(Intercept) -0.02902 0.46809 -0.062 0.954\nx1 -0.19239 0.66198 -0.291 0.790\nx2 0.40885 0.66198 0.618 0.581\n\nResidual standard error: 1.147 on 3 degrees of freedom\nMultiple R-squared: 0.1129, Adjusted R-squared: -0.4785 \nF-statistic: 0.1909 on 2 and 3 DF, p-value: 0.8355\n我错过了什么吗？我是否应该添加一个等于总平均值的虚拟变量，以便我可以将其用作参考水平？
\n我去年看到了一个类似的问题（但可能要求更高），但还没有解决方案：分类变量的所有系数具有“与平均值的差异”的模型；得到“对比编码”来做到这一点？。
\n这里接受的答案有效，但作者没有提供解释。我在统计SE上询问过这个问题：https ://stats.stackexchange.com/questions/600798/understanding-the-process-of-tweaking-contrasts-in-linear-model-fitting-to-show
\n这个答案向您展示如何获得以下系数表：
\n# Estimate Std. Error t value Pr(>|t|)\n#(Intercept) -0.02901982 0.4680937 -0.06199574 0.9544655\n#A -0.19238543 0.6619845 -0.29061922 0.7902750\n#B 0.40884591 0.6619845 0.61760645 0.5805485\n#C -0.21646049 0.6619845 -0.32698723 0.7651640\n太棒了，不是吗？它模仿你所看到的summary(fit)，即
# Estimate Std. Error t value Pr(>|t|)\n#(Intercept) -0.02901982 0.4680937 -0.06199574 0.9544655\n#x1 -0.19238543 0.6619845 -0.29061922 0.7902750\n#x2 0.40884591 0.6619845 0.61760645 0.5805485\n但现在我们显示了所有因素水平。
\nlm汇总不显示所有因素水平？2016年，我回答了这个Stack Overflow问题：“lm”摘要不显示所有因素级别，从那时起，它就成为标记类似主题的重复问题的目标。
\n回顾一下，基本思想是，为了获得用于最小二乘拟合的满秩设计矩阵，我们必须对因子变量应用对比。假设该因子有N 个水平，那么无论我们选择哪种类型的对比（请参阅?contrasts参考资料 中的列表），它都会将原始N 个虚拟变量减少为一组新的N - 1 个变量。因此，只有N - 1 个系数与N级因子相关。
然而，我们可以使用对比矩阵将N-1系数变换回原始N系数。该变换使我们能够获得所有因子水平的系数表。我现在将根据 OP 的可重现示例演示如何做到这一点：
\nset.seed(1)\ny <- rnorm(6, 0, 1)\nx <- factor(rep(LETTERS[1:3], each = 2))\nfit <- lm(y ~ x, contrasts = list(x = contr.sum))\n在此示例中，将和与零对比应用于因子x。要了解有关如何控制模型拟合对比度的更多信息，请参阅我的答案：如何使用 R 进行回归分析中的变量设置对比度？。
对于受归零对比影响的N 个级别的因子变量，我们可以使用以下函数来获取N x (N - 1)变换矩阵，该矩阵将估计的(N - 1)个系数映射lm回N个系数适合所有级别。
ContrSumMat <- function (fctr, sparse = FALSE) {\n if (!is.factor(fctr)) stop("\'fctr\' is not a factor variable!")\n N <- nlevels(fctr)\n Cmat <- contr.sum(N, sparse = sparse)\n dimnames(Cmat) <- list(levels(fctr), seq_len(N - 1))\n Cmat\n}\n对于示例 3 水平因子x，该矩阵为：
Cmat <- ContrSumMat(x)\n# 1 2\n#A 1 0\n#B 0 1\n#C -1 -1\n拟合模型fit报告该因子的 3 - 1 = 2 个系数。我们可以将它们提取为：
## coefficients After Contrasts\ncoef_ac <- coef(fit)[2:3]\n# x1 x2 \n#-0.1923854 0.4088459 \n因此，特定级别的系数为：
\n## coefficients Before Contrasts\ncoef_bc <- (Cmat %*% coef_ac)[, 1]\n# A B C \n#-0.1923854 0.4088459 -0.2164605 \n\n## note that they sum to zero as expected\nsum(coef_bc)\n#[1] 0\n同理，对比后可以得到协方差矩阵：
\nvar_ac <- vcov(fit)[2:3, 2:3]\n# x1 x2\n#x1 0.4382235 -0.2191118\n#x2 -0.2191118 0.4382235\n并将其转换为对比之前的：
\nvar_bc <- Cmat %*% var_ac %*% t(Cmat)\n# A B C\n#A 0.4382235 -0.2191118 -0.2191118\n#B -0.2191118 0.4382235 -0.2191118\n#C -0.2191118 -0.2191118 0.4382235\n解释：
\n模型截距coef(fit)[1]是总平均值。
计算出的coef_bc是每个级别的平均值与总平均值的偏差。
的对角线条目var_bc给出了这些偏差的估计方差。
然后我们可以继续计算这些系数的 t 统计量和 p 值，如下所示。
\n## standard error of point estimate `coef_bc`\nstd.error_bc <- sqrt(diag(var_bc))\n# A B C \n#0.6619845 0.6619845 0.6619845 \n\n## t-statistics (Null Hypothesis: coef_bc = 0)\nt.stats_bc <- coef_bc / std.error_bc\n# A B C \n#-0.2906192 0.6176065 -0.3269872 \n\n## p-values of the t-statistics\np.value_bc <- 2 * pt(abs(t.stats_bc), df = fit$df.residual, lower.tail = FALSE)\n# A B C \n#0.7902750 0.5805485 0.7651640 \n\n## construct a coefficient table that mimics `coef(summary(fit))`\nstats.tab_bc <- cbind("Estimate" = coef_bc,\n "Std. Error" = std.error_bc,\n "t value" = t.stats_bc,\n "Pr(>|t|)" = p.value_bc)\n# Estimate Std. Error t value Pr(>|t|)\n#A -0.1923854 0.6619845 -0.2906192 0.7902750\n#B 0.4088459 0.6619845 0.6176065 0.5805485\n#C -0.2164605 0.6619845 -0.3269872 0.7651640\n我们还可以通过包含总平均值（即模型截距）的结果来增强它。
\n## extract statistics of the intercept\nintercept.stats <- coef(summary(fit))[1, , drop = FALSE]\n# Estimate Std. Error t value Pr(>|t|)\n#(Intercept) -0.02901982 0.4680937 -0.06199574 0.9544655\n\n## augment the coefficient table\nstats.tab <- rbind(intercept.stats, stats.tab_bc)\n# Estimate Std. Error t value Pr(>|t|)\n#(Intercept) -0.02901982 0.4680937 -0.06199574 0.9544655\n#A -0.19238543 0.6619845 -0.29061922 0.7902750\n#B 0.40884591 0.6619845 0.61760645 0.5805485\n#C -0.21646049 0.6619845 -0.32698723 0.7651640\n我们还可以打印带有重要星号的表格。
\nprintCoefmat(stats.tab)\n# Estimate Std. Error t value Pr(>|t|)\n#(Intercept) -0.02902 0.46809 -0.0620 0.9545\n#A -0.19239 0.66199 -0.2906 0.7903\n#B 0.40885 0.66199 0.6176 0.5805\n#C -0.21646 0.66199 -0.3270 0.7652\n嗯？为什么没有星星？嗯，在这个例子中，所有 p 值都非常大。如果 p 值很小，星星就会出现。这是一个令人信服的演示：
\nfake.tab <- stats.tab\nfake.tab[, 4] <- fake.tab[, 4] / 100\nprintCoefmat(fake.tab)\n# Estimate Std. Error t value Pr(>|t|) \n#(Intercept) -0.02902 0.46809 -0.0620 0.009545 **\n#A -0.19239 0.66199 -0.2906 0.007903 **\n#B 0.40885 0.66199 0.6176 0.005805 **\n#C -0.21646 0.66199 -0.3270 0.007652 **\n#---\n#Signif. codes: 0 \xe2\x80\x98***\xe2\x80\x99 0.001 \xe2\x80\x98**\xe2\x80\x99 0.01 \xe2\x80\x98*\xe2\x80\x99 0.05 \xe2\x80\x98.\xe2\x80\x99 0.1 \xe2\x80\x98 \xe2\x80\x99 1\n哦，这太漂亮了。有关这些星星的含义，请参阅我的回答：Interpeting Rsignificancecodes for ANOVA table?
\n应该可以编写一个函数（甚至R包）来执行此类表转换。然而，可能需要付出很大的努力才能使此类函数足够灵活，以处理：
\n所有类型的对比（这很容易做到）；
\n复杂的模型术语，例如因子与其他数字/因子变量之间的相互作用（这似乎确实涉及！！）。
\n所以，我暂时就停在这里。
\n\n\n我从 lm 的摘要中获得的模型分数是否仍然准确，即使它没有显示该因素的所有级别？
\n
是的，他们是。lm进行精确的最小二乘拟合。
另外，系数表的变换不影响R平方、自由度、残差、拟合值、F统计量、ANOVA表等。
\n| 归档时间： |
|
| 查看次数： |
1185 次 |
| 最近记录： |