分类变量的所有系数均具有“与平均值的差异”的模型；得到“对比编码”来做到这一点？
daa*_*onr 10 r linear-regression lm
假设我们想做一个简单的“收入描述模型”。假设我们有三个组，北部、中部和南部（想想美国地区）。\n比较其他相似的组，假设北部的平均收入为 130，中部为 80，南部为 60。假设组规模相等，则平均值为90.
\n在（线性回归）模型中，应该有一种方法将系数报告为与总体平均值的差异（在多变量上下文中，\xe2\x80\x98所有其他都等于\xe2\x80\x99）并为每个系数获取一个：
\n$\\beta_{北} = 40$
\n$\\beta_{中心} = -10$
\n$\\beta_{南} = -30$
\n显然，要跳过拦截。
\n这对我来说似乎是最直观的。但我可以\xe2\x80\x99t在我的一生中弄清楚如何使用R\xe2\x80\x99s\xe2\x80\x98对比编码\xe2\x80\x99来得到这个。（而且，这似乎弄乱了变量名称）。
\n为我的模拟/mwe 设置参数
\nm_inc <- 90\nb_n <- 40\nb_c <- -10\nb_s <- -30\n\nsd_prop <- 0.5 #sd as share of mean\npop_per <- 1000 \n模拟数据
\n\nset.seed(100)\n\nn_income <- rnorm(pop_per, m_inc + b_n, (m_inc + b_n)*sd_prop)\nc_income <- rnorm(pop_per, m_inc + b_c, (m_inc + b_s)*sd_prop)\ns_income <- rnorm(pop_per, m_inc + b_s, (m_inc + b_s)*sd_prop)\n\nnoise_var <- rnorm(pop_per*3, 0, (m_inc + b_s)*sd_prop)\n\ni_df <- tibble(\n region = rep( c("n", "c", "s"), c(pop_per, pop_per, pop_per) ),\n income = c(n_income, c_income, s_income),\n noise_var\n) %>% \n mutate(region = as.factor(region))\n\n\ni_df %>% # Summary by group using purrr\n split(.$region) %>%\n purrr::map(summary)\n\n看起来足够接近。
\n现在我想“建立收入模型”来检查“控制其他因素”按地区划分的差异。为了说明这个问题，让我们将南方作为基地组。我设置了默认值contr.treatment，以防您重置它。
\ni_df <- i_df %>% mutate(region = relevel(region, ref="s"))\noptions(contrasts = rep ("contr.treatment", 2))\n\n\n(\n basic_lm <- i_df %>% lm(income ~ region + noise_var, .)\n)\n标准的东西：截距（大约）是“基地组”、南方和系数的平均值regionc，regionn代表这些的相对调整，大约为+20和+70。
这是标准的“虚拟编码”，即“处理编码”，默认 i R。
\n我们可以将这个默认值（对于无序变量）调整为针对无序和有序变量的“总和对比编码”
\noptions(contrasts = rep ("contr.sum", 2))\n\n(\n basic_lm_cc <- i_df %>% lm(income ~ region + noise_var, .)\n)\n现在这似乎让我们得到了我们正在寻找的调整系数，但是
\n- \n
- 地区名称丢失；我怎么知道哪个是哪个？ \n
- 它显然报告了 s（南）和 c（中）的调整系数。不是很直观。 \n
无论我们如何重新调整区域以设置特定的基础组（我尝试过），情况似乎都是如此......系数不会改变。
\n我找到了解决方案，但这不是“正确的方法”。我对结果（收入）变量进行去均值处理，并强制截距为 0：
\n\ni_df %>% \n mutate(m_inc = mean(income)) %>% \n lm(income - m_inc ~ 0 + region + noise_var, .)\n\n耶！这就是我想要的，并且奇迹般地保留了变量名称。但这似乎是一种奇怪的方法。另请注意，使用上面的代码，无论我强制使用哪个对比矩阵，无论是求和还是处理，都会出现这组系数。
\n如何通过对比编码或其他工具以“正确的方式”完成此操作？
\n我们无法使用对比编码来实现此目的。在单向方差分析中，对比用于将N级因子编码为N - 1 个变量加上一个截距，因此仍然是N 个变量对N个变量。但尝试在模型中包含组均值的总均值和偏差是从N个变量到N + 1 个变量的重新参数化。这（即使我们找到一种方法来做到这一点）使得设计矩阵秩亏，并且lm/ aov/glm等将把一个变量放入NA。
一般来说，我们要做后续的统计分析。在这个答案中，我将总结对比编码的作用，并展示如何以四种方式比较组均值和总均值：我们自己编码、使用multcomp、emmeans和car。
\nlibrary(ggplot2)\nlibrary(car)\nlibrary(multcomp)\nlibrary(emmeans)\n设置
\n我将使用一个与您的示例类似的示例。
\n
SimData <- function (group.size, group.mean, group.variance) {\n ## number of groups\n ng <- length(group.size)\n if (ng > 5) stop("There is no need to experiment that many groups!")\n ## number of observations per group\n n <- sum(group.size)\n ## generate a factor variable \'f\' for these groups\n f <- rep.int(factor(sprintf("G%d", 1:ng)), group.size)\n ## simulate samples from each group\n mu <- rep.int(group.mean, group.size)\n se <- rep.int(sqrt(group.variance), group.size)\n y <- rnorm(n, mu, se)\n ## numerical covariate \'x\' with slope = 1\n lim <- sd(y)\n br <- seq.int(-lim, lim, length.out = ng + 1)\n interval <- cbind(br[-(ng + 1)], br[-1])\n interval <- interval[sample.int(ng), ]\n a <- rep.int(interval[, 1], group.size)\n b <- rep.int(interval[, 2], group.size)\n x <- runif(n, a, b)\n ## create data.frame\n data.frame(y = y + x, f = f, x = x)\n}\n\nset.seed(4891738) ## my Stack Overflow ID\ngroup.size <- c(100, 125, 150)\ngroup.mean <- c(130, 80, 60)\ngroup.variance <- 0.25 * group.mean\ndat <- SimData(group.size, group.mean, group.variance)\n\nggplot(data = dat, mapping = aes(x = x, y = y, colour = f)) + geom_point()\n
对每组均值和方差的天真计算具有误导性，因为我们与事实相差甚远！
\n## true values are 130, 80, 60\nwith(dat, tapply(y, f, mean))\n## G1 G2 G3 \n## 148.36985 60.38273 59.45486\n\n## true values are 32.5, 20 and 15\nwith(dat, tapply(y, f, var))\n## G1 G2 G3 \n## 67.37108 55.25185 41.69867\n对比编码
\n
## treatment coding\ncontr.treatment(3)\n## 2 3\n## 1 0 0\n## 2 1 0\n## 3 0 1\n\n## sum-to-zero coding\ncontr.sum(3)\n## [,1] [,2]\n## 1 1 0\n## 2 0 1\n## 3 -1 -1\n
fit.treatment <- lm(y ~ f + x, dat, contrasts = list(f = "contr.treatment"))\ncoef(fit.treatment)\n#(Intercept) fG2 fG3 x \n# 128.94609 -48.71525 -69.03433 1.03121 \n\nsummary(fit.treatment)\nanova(fit.treatment)\n\nfit.sum <- lm(y ~ f + x, dat, contrasts = list(f = "contr.sum"))\ncoef(fit.sum)\n#(Intercept) f1 f2 x \n# 89.696234 39.249860 -9.465391 1.031210 \n\nsummary(fit.sum)\nanova(fit.sum)\n请注意，尽管使用不同的对比编码会给出不同的回归系数，但它们实际上通过产生相同的拟合值而等效。
\nall.equal(fit.treatment$fitted.values, fit.sum$fitted.values)\n## [1] TRUE\n比较群体均值与总均值
\n
R 中的线性假设检验
\n1. 没有精美软件包的普通方法
\n
## 3 x 2 linear combination matrix\nwt.treatment <- matrix(c(-1, 2, -1, -1, -1, 2), nrow = 3) / 3\nwt.sum <- matrix(c(1, 0, -1, 0, 1, -1), nrow = 3)\n\nvanilla <- function (wt, beta.ind, lmfit) {\n ## beta coefficients and their covariance matrix\n beta <- coef(lmfit)[beta.ind]\n V <- vcov(lmfit)[beta.ind, beta.ind]\n ## linear combination and their standard errors\n MEAN <- c(wt %*% beta)\n ## get standard errors for sum of `LinearComb`\n SE <- sqrt(diag(wt %*% tcrossprod(V, wt)))\n ## perform t-test\n tscore <- MEAN / SE\n pvalue <- 2 * pt(abs(tscore), lmfit$df.residual, lower.tail = FALSE)\n ## return a matrix\n ans <- matrix(c(MEAN, SE, tscore, pvalue), ncol = 4L)\n colnames(ans) <- c("Estimate", "Std. Error", "t value", "Pr(>|t|)")\n printCoefmat(ans)\n}\n\nvanilla(wt.treatment, 2:3, fit.treatment)\n## Estimate Std. Error t value Pr(>|t|) \n## [1,] 39.24986 0.90825 43.215 < 2.2e-16 ***\n## [2,] -9.46539 0.89404 -10.587 < 2.2e-16 ***\n## [3,] -29.78447 0.32466 -91.741 < 2.2e-16 ***\n## ---\n## Signif. codes: 0 \'***\' 0.001 \'**\' 0.01 \'*\' 0.05 \'.\' 0.1 \' \' 1\n\nvanilla(wt.sum, 2:3, fit.sum) ## identical to above\n2. 使用包\'multcomp\'
\n
## pad columns of zeros to zero out the effect of alpha and gamma\nwt.treatment0 <- cbind(0, wt.treatment, 0)\nwt.sum0 <- cbind(0, wt.sum, 0)\n\nsummary(glht(fit.treatment, linfct = wt.treatment0))\n## Linear Hypotheses:\n## Estimate Std. Error t value Pr(>|t|) \n## 1 == 0 39.2499 0.9083 43.22 <2e-16 ***\n## 2 == 0 -9.4654 0.8940 -10.59 <2e-16 ***\n## 3 == 0 -29.7845 0.3247 -91.74 <2e-16 ***\n## ---\n## Signif. codes: 0 \'***\' 0.001 \'**\' 0.01 \'*\' 0.05 \'.\' 0.1 \' \' 1\n## (Adjusted p values reported -- single-step method)\n\nsummary(glht(fit.sum, linfct = wt.sum0)) ## identical to above\n3.使用包\xe2\x80\x98emmeans\xe2\x80\x99
\n
emmeans(fit.treatment, specs = eff ~ f)\n## $emmeans\n## f emmean SE df lower.CL upper.CL\n## G1 127.3 1.002 371 125.4 129.3\n## G2 78.6 0.874 371 76.9 80.3\n## G3 58.3 0.379 371 57.5 59.0\n## \n## Confidence level used: 0.95 \n## \n## $contrasts\n## contrast estimate SE df t.ratio p.value\n## G1 effect 39.25 0.908 371 43.215 <.0001\n## G2 effect -9.47 0.894 371 -10.587 <.0001\n## G3 effect -29.78 0.325 371 -91.741 <.0001\n## \n## P value adjustment: fdr method for 3 tests\n\nemmeans(fit.sum, specs = eff ~ f) ## identical to above\n这里，使用specs = ~f或specs = "f"，仅$emmeans报告（边际均值）分量。左侧的“eff”意味着eff.emmc()要调用它来应用对比，并且该$contrasts组件呈现这样的结果。
4.使用包\xe2\x80\x98car\xe2\x80\x99
\ncarslinearHypothesis()的函数执行F 检验来测试所有线性组合是否同时为 0。因此它与上面演示的 t 检验不同。此外，它还给出错误：
linearHypothesis(fit.treatment, hypothesis.matrix = wt.treatment0)\n#Error in solve.default(vcov.hyp) : \n# system is computationally singular: reciprocal condition number = 1.12154e-17\n\nlinearHypothesis(fit.sum, hypothesis.matrix = wt.sum0) ## identical to above\n\n
这个答案从某种意义上来说是我之前几个答案的总结。
\n- \n
- \n \n
- \n \n
- \n \n
- \n \n
- \n \n
- \n \n
| 归档时间： |
|
| 查看次数： |
767 次 |
| 最近记录： |