10
$\begingroup$

I am trying to replicate a colleague's work and am moving the analysis from Stata to R. The models she employs invoke the "cluster" option within the nbreg function to cluster the standard errors.

See http://repec.org/usug2007/crse.pdf for a fairly complete description of the what and why of this option

My question is how to invoke this same option for negative binomial regression within R?

The primary model in our paper is specified in Stata as follows

 xi: nbreg cntpd09 logpop08 pcbnkthft07 pccrunion07 urbanpop pov00 pov002 edu4yr ///
 black04 hispanic04 respop i.pdpolicy i.maxloan rollover i.region if isser4 != 1,   
 cluster(state)

and I have replaced this with

pday<-glm.nb(cntpd09~logpop08+pcbnkthft07+pccrunion07+urbanpop+pov00+pov002+edu4yr+
black04+hispanic04+respop+as.factor(pdpolicy)+as.factor(maxloan)+rollover+
as.factor(region),data=data[which(data$isser4 != 1),])

which obviously lacks the clustered errors piece.

Is it possible to do an exact replication? If so how? If not, what are some reasonable alternatives?

Thanks

[Edit] As noted in the comments, I was hoping for a solution that didn't take me into the realm of multilevel models. While my training allows me to see that these things should be related, it is more of a leap than I am comfortable taking on my own. As such I kept digging and found this link: http://landroni.wordpress.com/2012/06/02/fama-macbeth-and-cluster-robust-by-firm-and-time-standard-errors-in-r/

that points to some fairly straightforward code to do what I want:

library(lmtest)
pday<-glm.nb(cntpd09~logpop08+pcbnkthft07+pccrunion07+urbanpop+pov00+pov002+edu4yr+
 black04+hispanic04+respop+as.factor(pdpolicy)+as.factor(maxloan)+rollover+
 as.factor(region),data=data[which(data$isser4 != 1),])
summary(pday)

coeftest(pday, vcov=function(x) vcovHC(x, cluster="state", type="HC1"))

This doesn't replicate the results from the analysis in Stata though, probably because it is designed to work on OLS not negative binomial. So the search goes on. Any pointers on where I am going wrong would be much appreciated

$\endgroup$
  • 3
    $\begingroup$ You might find Ben Bolker's notes useful here. $\endgroup$ – fmark Jun 14 '12 at 6:50
  • $\begingroup$ And see also this earlier question $\endgroup$ – fmark Jun 14 '12 at 7:40
  • $\begingroup$ FYI here is a definition of Stata's robust clustered standard errors. They don't seem that arduous to implement. IMO you may be better off with bootstrapped or jackknifed standard errors anyway (see the help on vce). I can't suggest any R packages though. Good luck on finding a replacement! $\endgroup$ – Andy W Jun 14 '12 at 17:53
  • $\begingroup$ Thanks @fmark - very useful comments, much better than my "answer" and I've updated it accordingly. $\endgroup$ – Peter Ellis Jun 14 '12 at 20:01
  • $\begingroup$ Thanks to all. I think the short answer to my question is that there is no straightforward replacement (eg, a pre-made function that exactly replaces the cluster option). Clearly someone with more experience can see the path through Ben Bolker's notes, but it takes me into new territory where I could not be sure I was getting the formula statements correct. I am not sure what the appropriate way to say "Thank you" without accepting an answer is, but you do have my thanks, and the shortcomings are my own. $\endgroup$ – csfowler Jun 18 '12 at 19:33
4
$\begingroup$

This document shows how to get cluster SEs for a glm regression:

http://dynaman.net/R/clrob.pdf

$\endgroup$
  • $\begingroup$ I will have to give it a test comparison with the stata results, but this looks exactly like what I was hoping for. $\endgroup$ – csfowler Aug 15 '12 at 3:11
1
$\begingroup$

This is not a fully satisfactory answer...

I haven't tried it myself but it looks like the glmmADMB package might do what you want.

I will shamelessly pinch from @fmark's comment on the question and agree with him that Ben Bolker's notes are useful, as is this earlier question, which is not quite an exact duplicate but covers very similar issues.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.