Replacing Standard Errors in a Reg Model in R

5
votes

I am in search of a way to directly replace the standard errors in a regression model with my own standard errors in order to use the robust model in another R package that does not come with its own robust option and can only be fed particular types of models and not coeftest formats.

Let's say I have a linear model:

model <- lm(data=dat, Y ~ X1 + X2 + X3)

I then want robust standard errors:

robust <- coeftest(model, vcov=sandwich)

Next, I need to use this model in a specific package that cannot be fed a coeftest and does not have its own robust standard errors option. I would like to replace the original model's standard errors (and, for that matter, p-values, t-statistics, etc.) prior to feeding it into the package so that they are accounted for.

To access the standard errors in the original model, I use: 

summary(model)$coefficients[,2]

To extract the standard errors from the coeftest, I use:

coeftest.se <- robust[, 2]

The following method, however, returns an error when trying to replace the model's standard errors because it is treating "summary" itself as a command:

summary(model)$coefficients[,2] <- coeftest.se


Error in summary(M3)$coefficients[, 2] <- seM3 : could not find function "summary<-"

The Specifics

I am trying to run mediation analysis with the Mediation R package. The package will do one way clustered standard errors with the "mediate" function, but I would like two way clustered standard errors.

To get two way clustered standard errors, I am using Mahmood Arai's mclx function (code can be found here (p. 4).

My idea is to feed the package's mediate function with models that already report the correct standard errors. According the documentation, the mediation package accepts the following classes of models: lm, polr, bayespolr, glm, bayesglm, gam, rq, survreg, and merMod.

2
rregressionstandard-error
So, summary(model) is a function of the model, not a part of the model. The standard errors aren't part of the model, the summary function calculates them on the fly (see summary.lm to see how). – Gregor Thomas
So the easiest solution is probably to modify your "specific package" to use robust standard errors, or maybe it could be fed an rlm object. Either way, can't help much without knowing what you're trying to use this in. – Gregor Thomas
It really depends on the function you'd be massing the model to. The lm object is a "living" object. It knows the data used to fit it and the original call. Modifying it's internal properties isn't guaranteed to result in the behavior you want. you really should provide a more complete reproducible example that shows what you are trying to do with the model down stream. – MrFlick
I have included some specific details. I was mostly thinking about this from a theoretical perspective, but I'd be happy to include a reproducible example soon. – firebird17139

2 Answers

3
votes

What worked for me was:

    model <- lm(data=dat, Y ~ X1 + X2 + X3)

    library(sandwich)

    SE_robust <- sqrt(diag(vcovHC(model, type="HC2")))

    model2 <- summary(model)

    model2$coefficients[,2] <- SE_robust
0
votes

@MySeppo's answer is great, but I think it is worth giving a function here that changes everything to be robust (SEs, t-stats, p-values) for more generality (I'm also using HC1 to be consistent with Stata):

model <- lm(mpg~cyl,data=mtcars)

robustsummary <- function(model) { 
  library(sandwich)
  library(lmtest)
  coeftest <- coeftest(model, vcov = vcovHC(model, type="HC1"))
  summ <- summary(model)
  summ$coefficients[,2] <- coeftest[,2]
  summ$coefficients[,3] <- coeftest[,3]
  summ$coefficients[,4] <- coeftest[,4]
  summ
}
summary(model)
robustsummary(model)

Or if you don't want to use coeftest

robustsummary <- function(model) { 
  library(sandwich)
  SE_robust <- sqrt(diag(vcovHC(model, type="HC1")))
  summ <- summary(model)
  summ$coefficients[,2] <- SE_robust
  summ$coefficients[,3] <- summ$coefficients[,1]/SE_robust # get t-stats
  summ$coefficients[,4] <- 2*(1-pt(abs(summ$coefficients[,3]),summ$df[2])) # apply t distribution with the right degrees of freedom to get p-values
  summ
}