2

Is there any possible way to get 95% CI for regression coefficients from the robust regression, as implemented in MASS::rlm?

# libraries needed
library(MASS)
library(stats)
library(datasets)

# running robust regression
(x <-
  MASS::rlm(formula = scale(Sepal.Length) ~ scale(Sepal.Width),
            data = iris))
#> Call:
#> rlm(formula = scale(Sepal.Length) ~ scale(Sepal.Width), data = iris)
#> Converged in 5 iterations
#> 
#> Coefficients:
#>        (Intercept) scale(Sepal.Width) 
#>        -0.03728607        -0.14343268 
#> 
#> Degrees of freedom: 150 total; 148 residual
#> Scale estimate: 1.06

# getting confidence interval for the regression coefficient
stats::confint(object = x,
               parm = "scale(Sepal.Width)",
               level = 0.95)
#>                    2.5 % 97.5 %
#> scale(Sepal.Width)    NA     NA
Spacedman
  • 92,590
  • 12
  • 140
  • 224
Indrajeet Patil
  • 4,673
  • 2
  • 20
  • 51
  • 1
    `confint.default(object = x, parm = "scale(Sepal.Width)", level = 0.95)`? – Dan Mar 07 '18 at 16:43
  • Woot! That worked! Had no idea there are two different functions called `stats::confint` and `stats::confint.default`. Can you post this as answer so that I can accept it? Thanks! – Indrajeet Patil Mar 07 '18 at 16:45

2 Answers2

4

Explicitly calling confint.default seems to provide good results, like this:

confint.default(object = x, parm = "scale(Sepal.Width)", level = 0.95)

#                        2.5 %     97.5 %
#scale(Sepal.Width) -0.3058138 0.01894847

Edit

confint uses method confint.lm when it is passed x because x is of class lm (as well as rlm). Calling confint.default explicitly avoids this. These two functions only differ in one line of code, as shown below:

confint.lm

fac <- qt(a, object$df.residual)

confint.default

fac <- qnorm(a)

The issue is that x$df.residual is NA and, consequently, qt(a, object$df.residual) produces an NA whereas qnorm(a) doesn't have this problem.

Dan
  • 11,370
  • 4
  • 43
  • 68
1

Late to the party, but note that CIs from the normal distribution will have lower than expected coverage for small sample sizes.

The residual degrees of freedom for the rlm object can be gotten from,

summary(x)$df[2]  # see code for MASS:::summary.rlm

To write your own confint method for rlm results, assign the df to the df.residual slot, and then call confint.lm:

confint.rlm <- function(object, ...){
  object$df.residual <- MASS:::summary.rlm(object)$df[2]
  confint.lm(object, ...)
}

Now confint behaves as expected, and is also based on Student's t:

confint(x)
                        2.5 %     97.5 %
(Intercept)        -0.2004593 0.12588715
scale(Sepal.Width) -0.3071526 0.02028719
abcdefg
  • 11
  • 1