rregressionlinear-regressionlmmlm

Get confidence intervals for regression coefficients of "mlm" object returned by `lm()`


I'm running a multivariate regression with 2 outcome variables and 5 predictors. I would like to obtain the confidence intervals for all regression coefficients. Usually I use the function lm but it doesn't seem to work for a multivariate regression model (object mlm).

Here's a reproducible example.

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

Any alternative way to do it? (I could just use the value of the standard error and multiply by the right critical t value, but I was wondering if there was an easier way to do it).


Solution

  • confint won't return you anything, because there is no "mlm" method supported:

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

    As you said, we can just plus / minus some multiple of standard error to get upper / lower bound of confidence interval. You were probably going to do this via coef(summary(mod)), then use some *apply method to extract standard errors. But my answer to Obtain standard errors of regression coefficients for an “mlm” object returned by lm() gives you a supper efficient way to get standard errors without going through summary. Applying std_mlm to your example model gives:

    se <- std_mlm(mod)
    #                 income   prestige
    #(Intercept) 1162.299027 3.54212524
    #education    103.731410 0.31612316
    #women          8.921229 0.02718759
    

    Now, we define another small function to compute lower and upper bound:

    ## 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
    

    It is easy to interpret this. For example, for response income, the 95%-confidence interval for all variables are

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