[R-sig-ME] Huber-White standard errors for glmer Poisson fits

[Steve Walker]

If it helps, I have just pushed a new lme4:::weights.merMod to the 
development version of lme4, that does allow a "working" option:

https://github.com/lme4/lme4/commit/9953d93d1837

Cheers,
Steve

On 2/17/2014, 12:08 PM, Ben Bolker wrote:
>    This looks reasonable to me (at least following Hardin and Hilbe's
> definition of the score residual).  For the Poisson distribution, it
> turns out that the score residuals and the raw residuals are the same
> (because (1/d(eta)/d(mu))) = Var(mu) in this case).
>
>    Here's a function (not thoroughly tested) that computes score
> residuals more generally, at the price of using some information about
> the internal structure of the object -- it's a fairly little translation
> of Hardin and Hilbe's definition.
>
> scoreResids <- function(object,...) {
>      ff <- object at resp$family
>      eta.mu <- 1/ff$mu.eta(predict(object))
>      y <- object at resp$y
>      mu <- object at resp$mu
>      (y-mu)/ff$var(mu)/eta.mu
> }
>
>   Otherwise, though, I haven't had the time to work through to see why
> working residuals are used in estfun.glm.  I'm guessing that the issue
> is actually a difference in the lme4:::weights.merMod function, which
> doesn't actually have a "working" option (but the option is silently
> ignored !!) -- I think that if one implemented 'working' rather than
> 'prior' weights, that this would probably work in parallel with estfun.glm.
>
>    I agree that your functions will be very useful.
>
>   Ben Bolker
>
>
> On 14-02-17 05:02 AM, Theodore Lytras wrote:
>> Dear list,
>>
>> I am a PhD student in epidemiology, trying to analyze a multi-site cohort
>> study with R. I have a series of mixed-effects Poisson model objects fitted
>> with glmer(), and I am trying to find a way to calculate Huber-White robust
>> standard errors for these.
>>
>> The only program I found that calculates Huber-White SEs for mixed-model fits
>> is GLLAMM in Stata, so I would like to reproduce this functionality in R.
>>
>> I took as a starting point the package "sandwich" by Achim Zeileis, which
>> unfortunately provides a sandwich covariance matrix only for glm objects. So I
>> tried to modify his code (functions estfun.glm(), bread.glm() and sandwich())
>> and the function robust.se() found in [1] to make it work with glmer objects,
>> simplifying things as much as possible (focusing only on Poisson fits).
>>
>> Thus I came up with the code in [2]. Function robustSEglmm() takes a glmer fit
>> with family="Poisson" and the associated clustering variable and tries to
>> compute Huber-White standard errors.
>>
>> What I found was the following: if estfunGlmm() uses the "working" residuals
>> (provided by residuals.merMod()), as is the case in the original estfun.glm()
>> and explained in the "sandwich" package documentation, the resulting SEs are
>> far too large, 6-10 times as those produced by GLLAMM.
>>
>> However, after reading a few things about the various types of residuals in
>> [3], I tried replacing the working residuals in estfunGlmm() with "score"
>> residuals, and to my surprise I obtained SEs very close to those returned by
>> GLLAMM (±1%).
>>
>> So my question is: is this approach statistically correct and valid?
>>
>> Since I can't say I really understand the maths involved (I'm a medical doctor
>> by profession), I turned to Achim Zeileis for help. Cutting and pasting from
>> his answer:
>>
>> [quote]
>>
>> Well, one needs to work out the "right" residuals: The estfun() needs to
>> be the derivative of the log-likelihood with respect to the parameters.
>> And the bread() essentially the second derivative.
>>
>> Given that the derivative of the log-likelihood is also called the score
>> function, it may well be that the "score residuals" that you mention above
>> are what you are looking for. But I'm not familiar enough with the model
>> classes and their implementation in R to know for sure.
>>
>> [/quote]
>>
>> So, could someone more familiar with the internals of lme4 shed more light on
>> the issue? If this approach to calculate Huber-White standard error is indeed
>> correct, I guess it would be useful to many people beside myself (I've seen
>> this discussed before on this list) and would extend the capabilities of R.
>>
>> Thank you all in advance,
>>
>> Theodore Lytras
>>
>>
>> [1] http://diffuseprior.wordpress.com/2013/01/12/the-cluster-bootstrap/
>> [2] https://gist.github.com/anonymous/9046934
>> [3]
>> http://books.google.com/books?id=tOeqO6Hs-6gC&lpg=PA53&vq=residuals&hl=el&pg=PA55#v=onepage&q&f=false
>>
>> _______________________________________________
>> R-sig-mixed-models at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>