Skip to content

sandwich variance estimation using glmer?

17 messages · T.D. Rudolph, Doran, Harold, Dimitris Rizopoulos +6 more

Original
Doran, Harold
# 
Out of curiosity, why would you want a sandwich estimator from lmer? That estimator is typically used when the likelihood is misspecified, but you still want standard errors that account for correlations among units within a cluster.

Since this is what lmer standard errors already account for, is there a need for the sandwich?
Dimitris Rizopoulos
#
On 11/3/2010 1:57 PM, Doran, Harold wrote:
well, it is possible that the random-effects structure that you have 
specified is not the correct one (i.e., in order to fully account for 
the correlations), and in this case it makes sense to use the sandwich 
estimator (of course, the sandwich estimator has its own problems, but 
this is probably another discussion...)

Best,
Dimitris

  
    

      
      
    

  


  


      

    

    

      
      

        
      

    

    

      
      

        


  

        

      David Atkins
    

    

      
      #
    

  


  

      
Tyler--

This question about "robust" (usually, heteroscedastically consistent) 
SE has come up before with respect to lme4 and lmer.  For example, see 
the following thread:

http://tolstoy.newcastle.edu.au/R/e3/help/07/10/0754.html

I have a vague memory of seeing a post from Achim Zeileis (the author of 
the "sandwich" package) mentioning why it might be challenging to extend 
his methods to lmer, but I couldn't (quickly) find it.

[BTW, you might be cautious about comparisons to SAS... R is a wonderful 
resource, but not universally guaranteed to do everything SAS can, and 
sometimes for good reasons.  Moreover, it is the result of an 
incredible, collective generosity of time and expertise (ie, nobody is 
making money here...).]

Hope that helps.

cheers, Dave

        
Tyler wrote:

        
Indeed, in this case the correlation structure of the random effects is not
fully appreciated or known, in which case the standard errors are likely
underestimated.  The use of sandwich estimators should render variance
estimates, and therefore inference, somewhat more realistic.  While this is
currently possible with GEEs, that approach does not ask the same question
as a GLMM (i.e. marginal or "population" estimates vs. conditional or
"subject-specific" estimates).

I used to think R updates were ahead of SAS upgrades in terms of new
approaches but apparently that is often not the case.  Does anyone have the
know-how required to implement this in R, or is there something I'm still
missing?

Best,
Tyler


On Wed, Nov 3, 2010 at 9:12 AM, Dimitris Rizopoulos <

        
d.rizopoulos at erasmusmc.nl> wrote:

        

        
> On 11/3/2010 1:57 PM, Doran, Harold wrote:

        
>
 >> Out of curiosity, why would you want a sandwich estimator from lmer? 
That
 >> estimator is typically used when the likelihood is misspecified, but you
 >> still want standard errors that account for correlations among units 
within
 >> a cluster.
 >>
 >> Since this is what lmer standard errors already account for, is there a
 >> need for the sandwich?
 >>
 >
 > well, it is possible that the random-effects structure that you have
 > specified is not the correct one (i.e., in order to fully account for the
 > correlations), and in this case it makes sense to use the sandwich 
estimator
 > (of course, the sandwich estimator has its own problems, but this is
 > probably another discussion...)
 >
 > Best,
 > Dimitris
 >
 >
 >
 >  -----Original Message-----
 >>> From: r-sig-mixed-models-bounces at r-project.org [mailto:
 >>> r-sig-mixed-models-
 >>> bounces at r-project.org] On Behalf Of Tyler Dean Rudolph
 >>> Sent: Tuesday, November 02, 2010 4:41 PM
 >>> To: r-sig-mixed-models at r-project.org
 >>> Subject: [R-sig-ME] sandwich variance estimation using glmer?
 >>>
 >>> Are there any current functionalities in R that permit estimation of
 >>> robust
 >>> sandwich variances based on lmer (mixed model) objects??  I'm aware of
 >>> the
 >>> sandwich package and gee implementations but to my knowledge these are
 >>> not
 >>> yet compatible with mixed model objects.
 >>>
 >>> Apparently these are already implemented in SAS....
 >>>
 >>> Tyler
 >>>
 >>>        [[alternative HTML version deleted]]
 >>>
 >>> _______________________________________________
 >>> 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
 >>
 >>
 > --
 > Dimitris Rizopoulos
 > Assistant Professor
 > Department of Biostatistics
 > Erasmus University Medical Center
 >
 > Address: PO Box 2040, 3000 CA Rotterdam, the Netherlands
 > Tel: +31/(0)10/7043478
 > Fax: +31/(0)10/7043014
 > Web: http://www.erasmusmc.nl/biostatistiek/
 >

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Andrew Robinson
    

    

      
      #
    

  


  

      
Hi Tyler,

I think that there's something that you're missing.

R is not motivated by comparisons with SAS or any package.  So, your
impression that R was ahead of SAS or behind SAS is mistaken, or at
least, it's your impression, so you are responsible for it.  R
responds exactly to the community's needs because the community
supports it.  If the functionality that you want isn't there, it's
because noone else has wanted it badly enough to

a) code it up, or

b) pay someone else to code it up.

If you want that function, and you know that SAS has it, then use SAS.
If you want to use that function in R, then see the above two points.

Good luck,

Andrew

        
On Wed, Nov 03, 2010 at 04:04:23PM -0400, Tyler Dean Rudolph wrote:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        
      

    

    

      
      

        


  

        

      Rolf Turner
    

    

      
      #
    

  


  

      
On 4/11/2010, at 10:49 AM, Andrew Robinson wrote:

        
<SNIP>

            

      
      
      
    

        
<SNIP>

There is another possibility that should be taken into consideration.
The person wanting that functionality may lack both the requisite skills
and insight to code up the functionality his or herself and the financial
resources to pay someone else to do this.

Yet another possibility is that the person responsible for the original
code simply never thought of including the requested functionality.  I am
sure I've seen instances where people have asked (on the R-help list) for
features in specific packages, and the maintainers of those packages have
responded ``Yes, that's a good idea.  I'll include that in the next release.''
Can't give specific examples off the top of my head.  However I know for
certain that I have accommodated such feature requests from users in respect
of packages that I maintain (specifically Iso and deldir).

	cheers,

		Rolf Turner

  
  


  


      

    

    

      
      

        


  

        

      Andrew Robinson
    

    

      
      #
    

  


  

      
Hi Tyler,

I guess that it is a sensitive topic.  This is a community of
volunteers.  Specious comparisons with commercial software products
aren't helpful, unless they're *positive* specious comparisons ;).  

The other problem with working in a community is that it's very
difficult for any one person to definitively state that functionality
does not exist and is not presently being worked on.  

So, I repeat, good luck ....

Andrew

        
On Wed, Nov 03, 2010 at 06:06:43PM -0400, Tyler Dean Rudolph wrote:

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Doran, Harold
    

    

      
      #
    

  


  

      
Let me push on this just a bit to spark further discussion. The OP was interested in robust standard errors given misspecification in the likelihood. So, one possible avenue was to explore Huber-White standard errors, or the sandwich estimator, to account for this misspecification and obtain "better" standard errors, but still use the point estimates of the fixed effects as given.

Some discussion on this has noted that the misspecification occurs in many ways, sometimes given that distributional assumptions were not met. Let's assume someone was willing and skilled to code up the HW as a possible solution within lmer to account for not meeting certain distributional assumptions.

My question is now why not directly code up models that permit for different distributional assumptions, such as t-distributions of residuals (random effects) or whatever the case might be? In other words, why not write code that addresses the problems directly (misspecification of the likelihood) rather than focusing on HW estimates.

Isn't it a better use of time and energy to focus on properly specifying the likelihood and estimating parameters from that model rather than HW?

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Ben Bolker
    

    

      
      #
    

  


  

      
On 10-11-04 10:18 AM, Doran, Harold wrote:

            

      
      
      
    

        
I would say there are two issues here, one philosophical and one
technical/pragmatic.

  Philosophical: how should one divide one's time and effort between
trying to come up with better parametric models vs. admitting that 'all
models are wrong' and coming up instead with robust alternatives that
perform reasonably well across a broad range of (unknown/unspecified)
parametric models?

  Technical/pragmatic: taking your case of t-distributions of residuals
-- this could be reasonably straightforward to implement in a framework
that is closer to a brute force computational solution (MCMC, AD Model
Builder), but I wouldn't know where to begin trying to implement it in
the elegant/computationally efficient lme4 framework ...

  cheers
    Ben Bolker

  
  


  


      

    

    

      
      

        


  

        

      David Atkins
    

    

      
      #
    

  


  

      
One slightly different perspective on robust SE in mixed models:

The place where I have seen these used regularly is in the HLM software 
(popular in education and psychology circles).  HLM *always* reports 
both "standard" and robust SE.

What I find interesting is that if you read Raudenbush and Bryk (and the 
HLM manual), they suggest using the robust SE as a model diagnostic (my 
term).  That is, when there is a discrepancy between SE, they rightly 
note that something is amiss, and you should do further detective work 
related to the random-effects specification.

That seems like a very valid use of robust SE, though I fully 
acknowledge such info (ie, model isn't fitting well) could be got other 
ways.

[BTW, I'd love to see other robust approaches, such as t-distributed 
error and/or priors, but as Ben notes that's an awfully high bar to 
implement -- either in lmer or MCMCglmm.  The heavy package is an 
initial attempt, but seems to be "stalled out" at the moment.]

For what it's worth.

cheers, Dave

        
Harold wrote:

        
Let me push on this just a bit to spark further discussion. The OP was 
interested in robust standard errors given misspecification in the 
likelihood. So, one possible avenue was to explore Huber-White standard 
errors, or the sandwich estimator, to account for this misspecification 
and obtain "better" standard errors, but still use the point estimates 
of the fixed effects as given.

Some discussion on this has noted that the misspecification occurs in 
many ways, sometimes given that distributional assumptions were not met. 
Let's assume someone was willing and skilled to code up the HW as a 
possible solution within lmer to account for not meeting certain 
distributional assumptions.

My question is now why not directly code up models that permit for 
different distributional assumptions, such as t-distributions of 
residuals (random effects) or whatever the case might be? In other 
words, why not write code that addresses the problems directly 
(misspecification of the likelihood) rather than focusing on HW estimates.

Isn't it a better use of time and energy to focus on properly specifying 
the likelihood and estimating parameters from that model rather than HW?

  
    

      
      
    

  


  


      

    

    

      
      

        
      

    

    

      
      

        
      

    

    

      
      

        


  

        

      David Atkins
    

    

      
      #
    

  


  

      
On 11/4/10 3:06 PM, Andrew Robinson wrote:

            

      
      
      
    

        
Andrew--

You are, of course, absolutely correct.  I use R almost exclusively, and 
for the most part, not having robust SE doesn't really cause any 
problems for me (in general, or with respect to checking diagnostics... 
actually, I love it!).

 From my own experience, I have seen the HLM's software use of both be 
helpful in the following way: I have several times had colleagues / 
students come to me with output from HLM, noting the discrepancy between 
SE and wondering what this means?  In essence, b/c that software puts 
both of them in front of you, it's hard to miss when they are different 
(esp. since p-values will then be radically different).

These times tend to be "teachable moments" in terms of what it is saying 
about the model and data (and fit).  More often that not, the outcome 
has been some form of count variable (which lurk quite commonly in 
psychological/psychiatric waters), and thus we can talk about 
distributional properties of the outcome vs. model, etc. (And I think 
the mean-variance relationship of the Poisson tends to lead to large 
adjustments with the robust SE.)

I would in no way suggest that lme4 should follow suit with the HLM 
software, and I honestly doubt whether it would serve a similar purpose, 
as more novice statistical users tend to be intimidated by R (given lack 
of menus, GUI, etc.).

As I mentioned before, what I *would* be interested in are robust 
approaches a la incorporating t-distribution in prior or likelihood.  I 
am currently working with a colleague on an analysis of a small number 
of groups, where we have discused this -- he's exploring WinBUGS as an 
option at present (to follow-up on our analyses with glmer and MCMCglmm).

For what it's worth (and realizing it ain't all that much about R...).

cheers, Dave

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Ben Bolker
    

    

      
      #
    

  


  

      
-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

        
On 10-11-05 12:12 PM, David Atkins wrote:

            

      
      
      
    

        
Channeling Dave Fournier here for a few moments ...

  For that particular project, I would recommend that you give AD Model
Builder and (to plug my own project and make things a little more R-ish)
the R2ADMB package a try. AD Model Builder is fast, flexible, uses
Laplace approximation for mixed modeling, and allows you to do post-hoc
MCMC sampling around the estimated modes; changing from a Gaussian to
t-distributed residual/likelihood calculations is in principle just a
matter of adding a parameter and changing the line of code that does the
likelihood calculation.  And ADMB is now open source.
  The downside is that you have to learn a new system (objective
functions in ADMB are written in a superset of C++). Overall I would say
that the effort, and the payoff, of learning ADMB are on the same order
of magnitude as learning WinBUGS/JAGS ...

  Ben Bolker
-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.4.10 (GNU/Linux)
Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org/

iEYEARECAAYFAkzUL+wACgkQc5UpGjwzenPStQCglIPjBsjlEFsea5CjUFugaIFW
q7sAn1M78DGVTEctM/VETxBu/HdrrDj6
=fTTC
-----END PGP SIGNATURE-----

  
  


  


      

    

    

      
      

        


  

        

      Whit Armstrong
    

    

      
      #
    

  


  

      
Do you have a winbugs example of the model you would like to fit?

If it is as simple as using a t-dist prior for your coefficients, then
it's pretty simple to do that in winbugs or PyMC.

Of course I'd like to plug my own project which is a c++ equivalent of
winbugs, now directly callable from R via inline.

I'm happy to translate a winbugs model and post the cppbugs equivalent.

-Whit

        
On Fri, Nov 5, 2010 at 12:25 PM, Ben Bolker <bbolker at gmail.com> wrote: