[R-sig-ME] lme4 heteroscedasticity???

[Paul Johnson]

I did this a few years ago for ~250 patients (nested within ~10 hospitals) being assessed by 7 adjudicators, so roughly a 250 x 6 table of crossed random effects. The aim was to look for differences in precision between clinicians rating patients on a symptoms scale. It did take a several minutes to fit on a server with 24 GB of RAM, and required a little fiddling with the control options, but was nevertheless worth the time and trouble. Looking back at my code (below), I find that your advice helped me to fit it, via a post of yours from 2004, for which thanks!

The model syntax is a little fiddly and easy to get wrong, so I’d advise anyone doing this for the first time to check that the results of the homoscedastic model match those from lme4, where the random effects syntax is straightforward, before fitting the heteroscedastic model.

       # fit models
  
        # extend number of iterations allowed
        
          control<-lmeControl(maxIter=500,msMaxIter=500,niterEM=250,msMaxEval=2000)
        
        # I want to fit three random effects, id crossed with adj and id nested within site
        # how to model this is  described by douglas bates here: http://tolstoy.newcastle.edu.au/R/help/04/02/0830.html
        # beware - lme takes several minutes to fit these models - lmer is much faster but doesn't allow heteroscedasticity to be modelled
    
          mod.sep.adj.rancrossed<-
            lme(mrs~1,random=list(Block=pdBlocked(list(pdIdent(~adj-1),pdIdent(~site-1),pdIdent(~id-1)))),
              data=cars.rep,control=control)
    
        # fit heteroscedastic model
    
          mod.sep.adj.rancrossed.w<-update(mod.sep.adj.rancrossed,weights=varIdent(form=~1|adj))

        # get p-value for heteroscedastic model vs homoscedastic model 
    
          (precision.p.val<-anova(mod.sep.adj.rancrossed.w,mod.sep.adj.rancrossed)[2,"p-value"])


On 13 Oct 2014, at 16:35, Douglas Bates <bates at stat.wisc.edu> wrote:

> On Mon, Oct 13, 2014 at 10:24 AM, Paul Johnson <paul.johnson at glasgow.ac.uk> wrote:
> You can model both heteroscedasticity and crossed random effects in nlme. Crossed random effects can be handled in nlme using pdBlocked and pdIdent. See p163-6 of Pinheiro, J.C. & Bates, D.M. (2000). Mixed-Effects Models in S and S-PLUS. Springer, New York.
> 
> But it will be a very slow and memory-intensive process unless one of the grouping factors has very few levels.  The nlme package does not use sparse matrix methods, nor does it use a profiled objective criterion.
> 
> In other words, yes it is possible but no, it is not easy to do so this way.
>  
> On 13 Oct 2014, at 16:12, Douglas Bates <bates at stat.wisc.edu> wrote:
> 
> > It is best to send questions like this to the
> > R-SIG-Mixed-Models at R-Project.org mailing list, which I am cc:ing on this
> > reply.  Several of those who read that list can respond more quickly than I
> > am able to.
> >
> > As far as I know there is not yet the capability in lme4 to model
> > heteroscedasticity in the distribution of the response given the random
> > effects.
> >
> > On Mon, Oct 13, 2014 at 6:01 AM, Košmelj, Katarina <
> > Katarina.Kosmelj at bf.uni-lj.si> wrote:
> >
> >> Hello,
> >>
> >> I am analyzing a mixed model with three crossed factors, two random
> >> (sample, laboratory) and one fixed (method); the response variable is the
> >> number of somatic cells in milk.  The main question is: is the precision of
> >> the means of the three method is comparable? Therefore, I would like to
> >> compare a model with different variances for the methods with the model
> >> considering the same variance for the methods.
> >>
> >>
> >>
> >> In nlme, this is feasible, however, two crossed random factors can not be
> >> tackled, this can be analyzed with lme4.
> >>
> >>
> >>
> >> In nlme, the problem of heteroscedasticity if solved. Is this problem
> >> solved in lme4 yet?
> >>
> >> Do you have any suggestion how to deal with this problem?
> >>
> >>
> >>
> >> Regards,
> >>
> >> Katarina
> >>
> >>
> >>
> >
> >       [[alternative HTML version deleted]]
> >
> > _______________________________________________
> > R-sig-mixed-models at r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> 
>