Does "coeftest" correctly use weights from "svydesign" in "svyglm" object

Question

I am using data from the European Social Survey (ESS) and I would like to calculate country-level cluster-robust standard errors for a regression model in R that includes country fixed effects and employs the design weights that come with the ESS.

To correctly use the weights, I use the 'survey' package and the functions 'svydesign' and 'svyglm'. This step looks like this:

design_1 <- svydesign(id=~1, weights=~dweight, data=ESS)

m1 <- svyglm(y ~ cntry + x, design = design_1)

My question is: when I now apply the functions 'cluster.vcov' and 'coeftest' from the packages 'lmtest' and 'multiwayvcov' to the model m1, do the resulting standard errors correctly account for the design weights? This step looks like this:

vcov_m1 <- cluster.vcov(m1, ESS$cntry)

coeftest(m1, vcov_m1)

Note that I do not use 'cntry' as an id variable in the svydesign function, because then I cannot include country dummies in the regression model.

Thanks in advance for your feedback!

Answer 1

You should include cntry as an id variable in the svydesign function; it won't stop you having country dummies in your model.

However, if you want to do it the hard way, fit the model as you have done, extract the influence functions, create the correct survey design object, and compute the variance of the total of the influence functions

efun <- model.matrix(m1) * resid(m1,"response")
infl <- efun %*% m1$naive.cov/mean(weights(design_1))
design_2 <- svydesign(id=~cntry, weights=~dweight, data=ESS)
svytotal(infl, design_2)

That's basically what's happening under the hood anyway. You should get numbers very close to zero for the point estimates, and the standard errors will be the standard errors for the model, under the new, correct design.

The idea is that $\hat\beta=\beta+\sum IF+o(n^{-1/2})$, so that inference about $\hat\beta$ can be replaced by inference about the influence functions, where we've got a population total as the estimand and we understand population totals. For a linear model, the remainder term is actually zero, not just small.