[R-sig-ME] Beginner help for mixed effects model

[Ben Bolker]

  Comments inline marked with BMB>

On Fri, Oct 21, 2016 at 10:52 AM, Isabella Mandl <i.mandl at gmx.at> wrote:
> Dear all,
> I apologise for the rather simplistic questions I am about to ask but I am
> very much at the beginning of my analysis and have only just worked my way
> through to mixed effect models in R. I collected some ecological data over
> the past years and have been advised by my PhD supervisors to use a GLMM. It
> should have been rather straightforward:
> I have measured "looking time" in response to different playback stimuli
> (0,1,2,3) in a group of 13 individuals - measures were repeated four times
> over the year. I now wanted to look at the effect of Stimulus Type, Season
> and Sex on "looking time" as well as look at whether different sexes have
> different looking times in different seasons (sex-season interaction).
> Random effects are ID of the animal and Trial.
> This is the model I built (using lme4):
>
> glmer(LookingTime~Stimulus+Season*Sex+(1|ID)+(1|Trial),
> family=Gamma(link="log"), data=playbacks)
> which gives me the following output:
>
> Generalized linear mixed model fit by maximum likelihood (Laplace
> Approximation) ['glmerMod']
>  Family: Gamma  ( log )
> Formula: VigilanceTowardsAdjus ~ Stimulus + Season * Sex + (1 | ID) +
> (1 | Trial)
>    Data: playbacks
>
>      AIC      BIC   logLik deviance df.resid
>   2731.9   2775.1  -1354.0   2707.9      258
>
> Scaled residuals:
>     Min      1Q  Median      3Q     Max
> -1.4529 -0.5723 -0.0321  0.6412  2.1003
>
> Random effects:
>  Groups   Name        Variance  Std.Dev.
>  Trial    (Intercept) 3.985e-01 0.6313040
>  ID       (Intercept) 6.185e-08 0.0002487
>  Residual             4.464e-01 0.6681655
> Number of obs: 270, groups:  Trial, 135; ID, 13

BMB>  this reflects essentially no among-ID variation once other variation is
taken into account.  (Not a big problem/surprising but worth noting)

>
> Fixed effects:
>                                        Estimate      Std. Error   t value
> Pr(>|z|)
> (Intercept)                       3.425875   0.184784  18.540 < 2e-16 ***
> Stimulus                           0.181798   0.036725   4.950 7.41e-07 ***
> Season[T.ED]                   0.382235   0.259361   1.474 0.1405
> Season[T.EW]                  0.338762   0.259361   1.306 0.1915
> Season[T.W]                    0.138404   0.259609   0.533 0.5939
> Sex[T.M]                           0.207001   0.272354   0.760 0.4472
> Season[T.ED]:Sex[T.M] -0.649771   0.387518  -1.677   0.0936 .
> Season[T.EW]:Sex[T.M] -0.428328   0.378500  -1.132   0.2578
> Season[T.W]:Sex[T.M]  -0.008209   0.385262  -0.021   0.9830
> ---
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> convergence code: 0
> unable to evaluate scaled gradient
> Model failed to converge: degenerate  Hessian with 1 negative eigenvalues
>
> The model fails to converge and gives me the following residual plot (which
> to me looks bad):

BMB>  The mailing list strips out graphical attachments.  Any chance you
can post the graph somewhere, e.g. imgur.com ?

  For failure to converge, see ?convergence

  I would strongly recommend that you try a log-Normal model as a
complement/backup for the Gamma model (i.e.
lmer(log(VigilanceTowardsAdjus) ~ ....) ) - in general log-LMMs
and Gamma-GLMMs give fairly similar results, and the former are a
little more stable.

>
> I don't really know what to look for to make it fit better. There is an
> inbalance in trial numbers (double the amount of 0-Stimulus trials than any
> of the others) and a slight inbalance in subjects (not all tested all four
> times) - could that have something to do with it?

  That shouldn't be a big deal.

> Should I be looking at PQL
> instead of ML?

Probably not ...

> Or is it the error distribution that I'm getting wrong?
> I am grateful for any help that points me into the right direction because I
> feel like I'm missing obvious things here.
> Kindest regards,
> Isabella

   Just guessing without the diagnostic plots, but I often look at the raw data
for things that might be violations of the model  (nonlinearities on
the log scale of the response to Stimulus - at the moment you're
assuming log-linear response), outliers ...)

  It's conceivable that the convergence warning is just wrong in the
singular fit case (which yours is very close to) - if you like, you
can send me your data (in our never-ending struggle to get the
convergence warnings for lme4 to be just right)