Cross posted on Stackoverflow with a bounty of 200.
EDIT:
I think I have to clarify this question a little bit more. So what I am looking for, is a function in which I can provide both the vcov matrix (the vcov2sls
), and have robust and clustered standard errors. However it seems that they both pertain to the same vcov matrix option. So if I supply one, I already have to make sure the se's are clustered and robust.. So I guess I am essentially asking how I can make the vcov2sls
function have robust and clustered errors. Obviously any other solution leading to the same practical outcome would be great as well.
END OF EDIT
A while ago, I asked this question, which was about correcting the standard errors when using IV/2SLS and the first stage has a tobit distribution, on which I got an amazing answer from jay.sf (example data at the bottom). He provided me with the following function:
vcov2sls <- function(s1, s2, data, type=2) {
## turn factor variables into dummies
DATA <- as.data.frame(model.matrix(phantom ~ ., transform(data, phantom=0)))
## list variable names
vn <- lapply(list(s1=s1, s2=s2), function(s)
c(all.vars(s$call)[1], colnames(model.matrix(s))[-1]))
## auxilliary model matrix
X <- cbind(`(Intercept)`=1, DATA[, c(vn$s1[1], vn$s2[-(1:2)]), F])
## get y
y <- DATA[, vn$s2[1]]
## betas second stage
b <- s2$coefficients
## calculate corrected sums of squares
sse <- sum((y - b %*% t(X))^2)
rmse <- sqrt(mean(s2$residuals^2)) ## RMSE 2nd stage
V0 <- vcov(s2) ## biased vcov 2nd stage
dof <- s2$df.residual ## degrees of freedom 2nd stage
## calculate corrected RMSE
rmse.c <- sqrt(sse/dof)
## calculate corrected vcov
V <- (rmse.c/rmse)^2 * V0
return(V)
}
It works great. Even when I want to use robust/clustered standard errors, that is not a problem, because AER::tobit
, calculates the robust/clustered standard errors within the function:
tobit(y~x, left=12, right=33, data=DT, robust=robust, cluster=cluster)
However I want to use jay.sf's function, when the first stage is an lm
, the clustering takes part in the summary (source), for example:
first_stage_ols <- lm(y~x, data=DT)
summary(first_stage_ols, robust=T)
Is there either, a way to correct the standard errors from within the lm function, or (replaced them in the result), or adapt the vcov2sls
function to also account for robust/clustered standard errors?
EDIT: I know that also lmtest:coeftest
exists, but I want to able to use weights
. See this link. I am having trouble figuring out if this is possible in lmtest:coeftest
.
Example Data
DF <- structure(list(country = c("C", "C", "C", "C", "J", "J", "B",
"B", "F", "F", "E", "E", "D", "D", "F", "F", "I", "I", "J", "J",
"E", "E", "C", "C", "I", "I", "I", "I", "I", "I", "C", "C", "H",
"H", "J", "J", "G", "G", "J", "J", "I", "I", "C", "C", "D", "D",
"A", "A", "G", "G", "E", "E", "J", "J", "G", "G", "I", "I", "I",
"I", "J", "J", "G", "G", "E", "E", "G", "G", "E", "E", "F", "F",
"I", "I", "B", "B", "E", "E", "H", "H", "B", "B", "A", "A", "I",
"I", "I", "I", "F", "F", "E", "E", "I", "I", "J", "J", "D", "D",
"F", "F"), year = c(2005, 2010, 2010, 2005, 2005, 2010, 2010,
2005, 2010, 2005, 2005, 2010, 2010, 2005, 2005, 2010, 2005, 2010,
2005, 2010, 2010, 2005, 2010, 2005, 2005, 2010, 2005, 2010, 2010,
2005, 2010, 2005, 2005, 2010, 2010, 2005, 2005, 2010, 2005, 2010,
2005, 2010, 2005, 2010, 2010, 2005, 2005, 2010, 2010, 2005, 2010,
2005, 2010, 2005, 2010, 2005, 2010, 2005, 2010, 2005, 2010, 2005,
2010, 2005, 2010, 2005, 2010, 2005, 2005, 2010, 2005, 2010, 2005,
2010, 2005, 2010, 2005, 2010, 2005, 2010, 2010, 2005, 2005, 2010,
2005, 2010, 2010, 2005, 2010, 2005, 2010, 2005, 2005, 2010, 2005,
2010, 2010, 2005, 2010, 2005), sales = c(15.48, 12.39, 3.72,
23.61, 4, 31.87, 25.33, 7.64, -0.26, 2.9, 15.48, 12.39, 3.72,
23.61, 4, 31.87, 25.33, 7.64, -0.26, 2.9, 15.48, 12.39, 3.72,
23.61, 4, 31.87, 25.33, 7.64, -0.26, 2.9, 15.48, 12.39, 3.72,
23.61, 4, 31.87, 25.33, 7.64, -0.26, 2.9, 15.48, 12.39, 3.72,
23.61, 4, 31.87, 25.33, 7.64, -0.26, 2.9, 15.48, 12.39, 3.72,
23.61, 4, 31.87, 25.33, 7.64, -0.26, 2.9, 15.48, 12.39, 3.72,
23.61, 4, 31.87, 25.33, 7.64, -0.26, 2.9, 15.48, 12.39, 3.72,
23.61, 4, 31.87, 25.33, 7.64, -0.26, 2.9, 15.48, 12.39, 3.72,
23.61, 4, 31.87, 25.33, 7.64, -0.26, 2.9, 15.48, 12.39, 3.72,
23.61, 4, 31.87, 25.33, 7.64, -0.26, 2.9), industry = c("D",
"D", "E", "E", "F", "F", "F", "F", "D", "D", "E", "E", "D", "D",
"E", "E", "F", "F", "F", "F", "D", "D", "F", "F", "E", "E", "D",
"D", "D", "D", "E", "E", "F", "F", "D", "D", "E", "E", "E", "E",
"D", "D", "E", "E", "D", "D", "D", "D", "E", "E", "D", "D", "F",
"F", "D", "D", "D", "D", "E", "E", "D", "D", "E", "E", "D", "D",
"D", "D", "D", "D", "F", "F", "F", "F", "E", "E", "D", "D", "E",
"E", "F", "F", "E", "E", "F", "F", "E", "E", "F", "F", "D", "D",
"D", "D", "D", "D", "D", "D", "F", "F"), urbanisation = c("B",
"B", "A", "A", "B", "B", "A", "A", "C", "C", "C", "C", "A", "A",
"B", "B", "C", "C", "A", "A", "C", "C", "B", "B", "A", "A", "A",
"A", "A", "A", "A", "A", "A", "A", "C", "C", "B", "B", "B", "B",
"B", "B", "C", "C", "A", "A", "B", "B", "B", "B", "A", "A", "B",
"B", "A", "A", "A", "A", "B", "B", "C", "C", "A", "A", "C", "C",
"A", "A", "B", "B", "A", "A", "B", "B", "B", "B", "B", "B", "C",
"C", "A", "A", "A", "A", "A", "A", "A", "A", "C", "C", "A", "A",
"B", "B", "A", "A", "B", "B", "B", "B"), size = c(1, 1, 5, 5,
5, 5, 1, 1, 1, 1, 5, 5, 5, 5, 2, 2, 2, 2, 5, 5, 1, 1, 1, 1, 5,
5, 5, 5, 4, 4, 5, 5, 5, 5, 4, 4, 2, 2, 5, 5, 1, 1, 1, 1, 2, 2,
1, 1, 2, 2, 5, 5, 1, 1, 3, 3, 2, 2, 2, 2, 5, 5, 4, 4, 1, 1, 5,
5, 2, 2, 5, 5, 2, 2, 2, 2, 4, 4, 3, 3, 4, 4, 3, 3, 3, 3, 3, 3,
5, 5, 3, 3, 2, 2, 3, 3, 1, 1, 5, 5), base_rate = c(14L, 14L,
14L, 14L, 19L, 19L, 30L, 30L, 20L, 20L, 29L, 29L, 20L, 20L, 20L,
20L, 24L, 24L, 19L, 19L, 29L, 29L, 14L, 14L, 24L, 24L, 24L, 24L,
24L, 24L, 14L, 14L, 17L, 17L, 19L, 19L, 33L, 33L, 19L, 19L, 24L,
24L, 14L, 14L, 20L, 20L, 23L, 23L, 33L, 33L, 29L, 29L, 19L, 19L,
33L, 33L, 24L, 24L, 24L, 24L, 19L, 19L, 33L, 33L, 29L, 29L, 33L,
33L, 29L, 29L, 20L, 20L, 24L, 24L, 30L, 30L, 29L, 29L, 17L, 17L,
30L, 30L, 23L, 23L, 24L, 24L, 24L, 24L, 20L, 20L, 29L, 29L, 24L,
24L, 19L, 19L, 20L, 20L, 20L, 20L), taxrate = c(12L, 14L, 14L,
12L, 21L, 18L, 30L, 30L, 20L, 20L, 29L, 30L, 20L, 20L, 20L, 20L,
24L, 24L, 21L, 18L, 30L, 29L, 14L, 12L, 24L, 24L, 24L, 24L, 24L,
24L, 14L, 12L, 18L, 19L, 18L, 21L, 33L, 32L, 21L, 18L, 24L, 24L,
12L, 14L, 20L, 20L, 22L, 25L, 32L, 33L, 30L, 29L, 18L, 21L, 32L,
33L, 24L, 24L, 24L, 24L, 18L, 21L, 32L, 33L, 30L, 29L, 32L, 33L,
29L, 30L, 20L, 20L, 24L, 24L, 30L, 30L, 29L, 30L, 18L, 19L, 30L,
30L, 22L, 25L, 24L, 24L, 24L, 24L, 20L, 20L, 30L, 29L, 24L, 24L,
21L, 18L, 20L, 20L, 20L, 20L), vote = c(0, 0, 0, 0, 1, 1, 1,
0, 1, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 0, 0, 0, 1, 0, 1,
1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 1, 1,
1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0,
1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 0, 0, 1, 0, 1, 1, 0, 1, 1, 0,
1, 0, 1, 1, 1, 1, 0, 1, 1), votewon = c(0, 0, 0, 0, 1, 0, 1,
0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1,
1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1,
0, 0, 0, 0, 1, 1, 0, 1, 0, 1, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0,
1, 1, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 1, 1, 0,
0, 0, 1, 1, 0, 1, 0, 1, 1)), class = "data.frame", row.names = c(NA,
-100L))
## convert variables to factors beforehand
DF[c(1, 2, 4, 5, 6, 9, 10)] <- lapply(DF[c(1, 2, 4, 5, 6, 9, 10)], factor)
s1.tobit <- AER::tobit(taxrate ~ votewon + industry + size + urbanisation + vote,
left=12, right=33, data=DF)
yhat <- fitted(s1.tobit)
s2.tobit <- lm(sales ~ yhat + industry + size + urbanisation + vote, data=DF)
lmtest::coeftest(s2.tobit, vcov.=vcov2sls(s1.tobit, s2.tobit, DF))