[R] dealing with heteroscedasticity in lmer: problem with the method weights

[Ben Bolker]

doris gomez <dodogomez <at> yahoo.fr> writes:

> 
> Dear lmer users,

  You should probably redirect this sort of question to
r-sig-mixed-models at r-project.org , rather than the generic
R help list.

> 
> The experiment includes 15 groups of (3 males and 1 female). 
> The female is characterized by its quality Q1
> and Q2. Each male of a group is characterized by the number of 
> MatingAttempts (with Poisson
> distribution). I want to examine if male mating attempts 
> depend on female quality.  I can see from graphic
> exploration  that the within-group heterogeneity of 
> male attempts increases with female quality Q1.

  Is this true above and beyond the expected increase in
variance due to the Poisson distribution of the data?
 
> When including the method weights in the function lmer, 
>  I get the message that variables' length varies and
> the model does not run.
> lmer(MatingAttempts~Q1+Q2+(1|Group),data=file,
> family="poisson",na.action=na.omit,
> REML=FALSE, weights=varExp(form=~Q1))

  The REML and weights arguments do not make sense in the context
of a GLMM fitted with [g]lmer.

> If I run the same model (fixed effects and random effects) with lme, it 
> works properly, which shows that there is no problem with data structure.
> lme(MatingAttempts~Q1+Q2,random=~1|Group,
> data=file,na.action=na.omit, method="ML", weights=varExp(form=~Q1))

  As you presumably know, this model does not incorporate the Poisson
distribution of the data.

> 
> I saw on the forum that lmer had problems in taking into account 
> variance heterogeneity. Yet, the messages
> were old and there are maybe new solutions.
> How can I correct the analyses for this problem of heteroscedasticity? 
> Should I normalise the within group variance before implementing the model? 
> And deal with the variance (as
> a new variable to explain) in another model?
> Is there another way to solve this problem?

  (1) lmer still does not allow for models of heteroscedasticity
as implemented in the weights() argument of lme, and may not for some
time ...
  (2) You _could_ use glmmPQL() from the MASS package, along the
lines of 

MASS::glmmPQL(MatingAttempts~Q1+Q2,random=~1|Group,
  family=poisson,
  data=file,na.action=na.omit, weights=varExp(form=~Q1))

 **HOWEVER**: it's not entirely clear what statistical model
this represents, or whether it makes sense.
 (3) I would recommend fitting the GLMM, with lmer, without the
heteroscedasticity, and then see whether the heteroscedasticity
remains in the Pearson residuals or not ... (don't forget to
check for overdispersion too).