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