0

I need to estimate a fractional (response taking values between 0 and 1) model with R. I also want to cluster the standard errors. I have found several examples in SO and elsewhere and I built this function based on my findings:

require(sandwich)
require(lmtest)
clrobustse <- function(fit, cl){
    M <- length(unique(cl))
    N <- length(cl)
    K <- fit$rank
    dfc <- (M/(M - 1))*((N - 1)/(N - K))
    uj <- apply(estfun(fit), 2, function(x) tapply(x, cl, sum))
    vcovCL <- dfc*sandwich(fit, meat = crossprod(uj)/N)
    coeftest(fit, vcovCL)
}

I estimate my model like this:

fit <- glm(dep ~ exp1 + exp2 + exp3, data = df, fam = quasibinomial("probit"))
clrobustse(fit, df$cluster)

Everything works fine and I get the results. However, I suspect that something is not right as the non-clustered version:

coeftest(fit)

gives the exact same standard errors. I checked that Stata reports and that displays different clustered errors. I suspect that I have misspecified the function clrobustse but I just don't know how. Any ideas about what could be going wrong here?

Antti
  • 1,263
  • 2
  • 16
  • 28
  • 1
    You probably don't need to write your own function. Try the `multiwayvcov` package, or `plm::vcovHC`. This question on [stats.se] may help: https://stats.stackexchange.com/questions/10017/standard-error-clustering-in-r-either-manually-or-in-plm –  May 11 '18 at 08:52
  • Thank you! `multiwayvcov` was a good find. There's some difference to Stata std. errors but only around fourth or fifth digit. Somehow `vcovHC` did not replicate Stata results very closely at all. – Antti May 11 '18 at 11:26

0 Answers0