获取组均值差的 p 值，无需使用新的参考水平重新拟合线性模型

当我们有一个带有因子变量的线性模型时X（有等级A, B, and C)

y ~ factor(X) + Var2 + Var3 

结果显示估计值XB and XC这就是差异B - A and C - A。 （假设参考文献是A).

如果我们想知道之间差异的 p 值B and C: C - B, 我们应该指定 B 或 C 作为参考组并重新运行模型。

我们能得到效应的 p 值吗B - A, C - A, and C - B一次？

您正在通过检查回归系数的某些线性组合的 p 值来寻找线性假设检验。根据我的回答：如何使用聚类协方差矩阵对回归系数进行线性假设检验？，我们只考虑系数之和，我将扩展该函数LinearCombTest处理更一般的情况，假设alpha作为变量的一些组合系数vars:

LinearCombTest <- function (lmObject, vars, alpha, .vcov = NULL) {
  ## if `.vcov` missing, use the one returned by `lm`
  if (is.null(.vcov)) .vcov <- vcov(lmObject)
  ## estimated coefficients
  beta <- coef(lmObject)
  ## linear combination of `vars` with combination coefficients `alpha`
  LinearComb <- sum(beta[vars] * alpha)
  ## get standard errors for sum of `LinearComb`
  LinearComb_se <- sum(alpha * crossprod(.vcov[vars, vars], alpha)) ^ 0.5
  ## perform t-test on `sumvars`
  tscore <- LinearComb / LinearComb_se
  pvalue <- 2 * pt(abs(tscore), lmObject$df.residual, lower.tail = FALSE)
  ## return a matrix
  form <- paste0("(", paste(alpha, vars, sep = " * "), ")")
  form <- paste0(paste0(form, collapse = " + "), " = 0")
  matrix(c(LinearComb, LinearComb_se, tscore, pvalue), nrow = 1L,
         dimnames = list(form, c("Estimate", "Std. Error", "t value", "Pr(>|t|)")))
  }

考虑一个简单的例子，我们对三个组进行了平衡设计A, B and C，组均值分别为 0、1、2。

x <- gl(3,100,labels = LETTERS[1:3])
set.seed(0)
y <- c(rnorm(100, 0), rnorm(100, 1), rnorm(100, 2)) + 0.1

fit <- lm(y ~ x)
coef(summary(fit))

#             Estimate Std. Error   t value     Pr(>|t|)
#(Intercept) 0.1226684 0.09692277  1.265631 2.066372e-01
#xB          0.9317800 0.13706949  6.797866 5.823987e-11
#xC          2.0445528 0.13706949 14.916177 6.141008e-38

Since A是参考电平，xB正在给予B - A while xC正在给予C - A。假设我们现在对组之间的差异感兴趣B and C, i.e., C - B， 我们可以用

LinearCombTest(fit, c("xC", "xB"), c(1, -1))

#                         Estimate Std. Error  t value     Pr(>|t|)
#(1 * xC) + (-1 * xB) = 0 1.112773  0.1370695 8.118312 1.270686e-14

注意，这个函数也可以方便地计算出组平均值B and C， 那是(Intercept) + xB and (Intercept) + xC:

LinearCombTest(fit, c("(Intercept)", "xB"), c(1, 1))

#                                 Estimate Std. Error  t value     Pr(>|t|)
#(1 * (Intercept)) + (1 * xB) = 0 1.054448 0.09692277 10.87926 2.007956e-23

LinearCombTest(fit, c("(Intercept)", "xC"), c(1, 1))

#                                 Estimate Std. Error  t value     Pr(>|t|)
#(1 * (Intercept)) + (1 * xC) = 0 2.167221 0.09692277 22.36029 1.272811e-65

替代解决方案lsmeans

再次考虑上面的玩具示例：

library(lsmeans)
lsmeans(fit, spec = "x", contr = "revpairwise")

#$lsmeans
# x    lsmean         SE  df    lower.CL  upper.CL
# A 0.1226684 0.09692277 297 -0.06807396 0.3134109
# B 1.0544484 0.09692277 297  0.86370603 1.2451909
# C 2.1672213 0.09692277 297  1.97647888 2.3579637
#
#Confidence level used: 0.95 
#
#$contrasts
# contrast estimate        SE  df t.ratio p.value
# B - A    0.931780 0.1370695 297   6.798  <.0001
# C - A    2.044553 0.1370695 297  14.916  <.0001
# C - B    1.112773 0.1370695 297   8.118  <.0001
#
#P value adjustment: tukey method for comparing a family of 3 estimates

The $lsmeans域返回边缘组均值，而$contrasts返回成对组平均差，因为我们使用了“revpairwise”对比。阅读第 32 页lsmeans对于之间的差异"pairwise" and "revpairwise".

这当然很有趣，因为我们可以与以下结果进行比较LinearCombTest。我们看到LinearCombTest做得正确。