3

我正在运行具有 2 个结果变量和 5 个预测变量的多元回归。我想获得所有回归系数的置信区间。通常我使用该函数lm，但它似乎不适用于多元回归模型（对象mlm）。

这是一个可重现的示例。

library(car)
mod <- lm(cbind(income, prestige) ~ education + women, data=Prestige)
confint(mod) # doesn't return anything.

有什么替代方法吗？（我可以只使用标准误差的值并乘以正确的临界 t 值，但我想知道是否有更简单的方法来做到这一点）。

4

3 回答 3

1

confint不会返回任何东西，因为不支持“mlm”方法：

methods(confint)
#[1] confint.default confint.glm*    confint.lm      confint.nls*

正如你所说，我们可以通过加/减标准误差的倍数来获得置信区间的上限/下限。您可能会通过 来执行此操作coef(summary(mod))，然后使用某种*apply方法提取标准错误。但是我对“获取返回的“mlm”对象的回归系数的标准误差的回答为您提供了一种超级有效的方法来获得标准误差而无需通过。应用于您的示例模型给出：lm()summarystd_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

这很容易解释。例如，对于 response income，所有变量的 95% 置信区间为

#(intercept)    (-3798.25140, 814.25546)
#  education    (739.05564, 1150.70689)
#      women    (-81.75738, -46.35407)
于 2016-10-07T23:13:50.937 回答
0

这来自 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")
于 2015-02-10T21:46:41.107 回答
0

最近（2018 年 7 月）在R-devel 列表上似乎已经讨论了这个问题，所以希望在 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

就个人而言，我喜欢它以一种干净的方式（使用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
于 2018-10-11T16:38:05.690 回答