Получить доверительные интервалы для коэффициентов регрессии объекта "mlm", возвращаемых функцией `lm()`

Question

Я управляю многомерной регрессией с 2 ​​исходными переменными и 5 предикторами. Я хотел бы получить доверительные интервалы для всех коэффициентов регрессии. Обычно я использую функцию lm но это не похоже на многомерную модель регрессии (объект mlm).

Вот воспроизводимый пример.

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

Есть ли альтернативный способ сделать это? (Я мог бы просто использовать значение стандартной ошибки и умножить на правильное критическое значение t, но мне было интересно, есть ли более простой способ сделать это).

Answer 1

confint не вернет вам ничего, потому что не поддерживается метод "mlm":

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

Как вы сказали, мы можем просто плюс / минус несколько кратных стандартной ошибки, чтобы получить верхнюю / нижнюю границу доверительного интервала. Вы, вероятно, собирались сделать это через coef(summary(mod)), а затем использовать некоторые *apply метод извлечения стандартных ошибок. Но мой ответ на Получение стандартных ошибок коэффициентов регрессии для объекта "mlm", возвращаемого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

Это легко интерпретировать. Например, для ответа income95-процентный доверительный интервал для всех переменных

#(intercept)    (-3798.25140, 814.25546)
#  education    (739.05564, 1150.70689)
#      women    (-81.75738, -46.35407)

Answer 2

Это происходит из примера Foregnet. Вы хотите 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")

Answer 3

Похоже, что это недавно обсуждалось (июль 2018 г.) в списке 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