我正在运行具有 2 个结果变量和 5 个预测变量的多元回归。我想获得所有回归系数的置信区间。通常我使用函数lm，但它似乎不适用于多元回归模型（对象mlm）。
这是一个可重现的例子。
library(car)
mod <- lm(cbind(income, prestige) ~ education + women, data=Prestige)
confint(mod) # doesn't return anything.
还有其他方法吗？ （我可以只使用标准误差的值并乘以正确的临界 t 值，但我想知道是否有更简单的方法。
【问题讨论】：
标签： r regression linear-regression lm mlm
这来自 predict.lm 示例。你想要interval = 'confidence' 选项。
x <- rnorm(15)
y <- x + rnorm(15)
predict(lm(y ~ x))
new <- data.frame(x = seq(-3, 3, 0.5))
predict(lm(y ~ x), new, se.fit = TRUE)
pred.w.clim <- predict(lm(y ~ x), new, interval = "confidence")
matplot(new$x, pred.w.clim,
lty = c(1,2,2,3,3), type = "l", ylab = "predicted y")
【讨论】：
confint 不会返回任何东西，因为不支持“mlm”方法：
methods(confint)
#[1] confint.default confint.glm* confint.lm confint.nls*
正如你所说，我们可以通过加/减标准误差的倍数来获得置信区间的上限/下限。您可能会通过coef(summary(mod)) 来执行此操作，然后使用一些*apply 方法来提取标准错误。但是my answer 到Obtain standard errors of regression coefficients for an “mlm” object returned by lm() 为您提供了一种超级有效的方法来获得标准错误，而无需通过summary。将std_mlm 应用于您的示例模型给出：
se <- std_mlm(mod)
# income prestige
#(Intercept) 1162.299027 3.54212524
#education 103.731410 0.31612316
#women 8.921229 0.02718759
现在，我们定义另一个小函数来计算下限和上限：
## add "mlm" method to generic function "confint"
confint.mlm <- function (model, level = 0.95) {
beta <- coef(model)
se <- std_mlm (model)
alpha <- qt((1 - level) / 2, df = model$df.residual)
list(lower = beta + alpha * se, upper = beta - alpha * se)
}
## call "confint"
confint(mod)
#$lower
# income prestige
#(Intercept) -3798.25140 -15.7825086
#education 739.05564 4.8005390
#women -81.75738 -0.1469923
#
#$upper
# income prestige
#(Intercept) 814.25546 -1.72581876
#education 1150.70689 6.05505285
#women -46.35407 -0.03910015
这很容易解释。例如，对于响应income，所有变量的 95% 置信区间为
#(intercept) (-3798.25140, 814.25546)
# education (739.05564, 1150.70689)
# women (-81.75738, -46.35407)
【讨论】：
confint.mlm 方法并不完全正确。相反，修补confint.lm 方法相对容易——从lm() 调用的方法是 也适用于多变量情况。我（作为 R 核心成员）现在已经这样做了，用于 R 的开发版本和“R 3.5.1 patched”。
这似乎已经在最近（2018 年 7 月）在 R-devel list 上讨论过，所以希望在 R 的下一个版本中它会得到修复。该列表中建议的解决方法是使用：
confint.mlm <- function (object, level = 0.95, ...) {
cf <- coef(object)
ncfs <- as.numeric(cf)
a <- (1 - level)/2
a <- c(a, 1 - a)
fac <- qt(a, object$df.residual)
pct <- stats:::format.perc(a, 3)
ses <- sqrt(diag(vcov(object)))
ci <- ncfs + ses %o% fac
setNames(data.frame(ci),pct)
}
测试：
fit_mlm <- lm(cbind(mpg, disp) ~ wt, mtcars)
confint(fit_mlm)
给予：
2.5 % 97.5 %
mpg:(Intercept) 33.450500 41.119753
mpg:wt -6.486308 -4.202635
disp:(Intercept) -204.091436 -58.205395
disp:wt 90.757897 134.198380
就我个人而言，我喜欢简洁的 tibble 方式（使用 broom::tidy 会更好，但目前有一个问题）
library(tidyverse)
confint(fit_mlm) %>%
rownames_to_column() %>%
separate(rowname, c("response", "term"), sep=":")
给予：
response term 2.5 % 97.5 %
1 mpg (Intercept) 33.450500 41.119753
2 mpg wt -6.486308 -4.202635
3 disp (Intercept) -204.091436 -58.205395
4 disp wt 90.757897 134.198380
【讨论】：