[R] GLMM - Am I trying the impossible?

From: Pedro J. Aphalo <pedro.aphalo_at_cc.jyu.fi>
Date: Thu 18 Aug 2005 - 19:20:02 EST


Dear all,

I have tried to calculate a GLMM fit with lmer (lme4) and glmmPQL (MASS), I also used glm for comparison.

I am getting very different results from different functions, and I suspect that the problem is with our dataset rather than the functions, but I would appreciate help in deciding whether my suspicions are right. If indeed we are attempting the wrong type of analysis, some guidance about what would be the right thing to do would be greatly appreciated.

The details:
The data:
The data are from the end-point of a survival experiment with fish. The design of the experiment is a 2 x 2 factorial, with each factor (Bacteria and Parasite) at two levels (yes and no). There were 16 fish in each tank, and the treatment was applied to the whole tank. There were in all 10 tanks (160 fish), with 2 tanks for controls (no/no), 2 tanks for (Parasite:yes/Bacteria:no) and 3 tanks for each of the other 2 treatments. A dead fish was considered a success, and a binomial family with the default logit link was used in the fits. No fish died in the control treatment (Is this the problem?).

Overall "probabilities" as dead/total for the four treatments were: Paras Bact Prob
no no 0

yes   no    0.0625
no    yes   0.208
yes   yes   0.458

We are interested in testing main effects and interaction, but the interaction is of special interest to us.

The data for "dead" are coded as 0/1 with 1 indicating a dead fish, and the file has one row per fish.

Some results:

lme4 (ver 0.98-1, R 2.1.1, Windows XP)


> fish1.lmerPQL <- lmer(dead ~ Parasite * Bacteria + (1|Tank),
data=fish.data, family=binomial)
Error in lmer(dead ~ Parasite * Bacteria + (1 | Tank), data = fish.data, :

         Unable to invert singular factor of downdated X'X In addition: Warning message:
Leading minor of size 4 of downdated X'X is indefinite
>

without the interaction:
> fish3.lmerPQL <- lmer(dead ~ Parasite + Bacteria + (1|Tank),
data=fish.data, family=binomial)
> anova(fish3.lmerPQL)

Analysis of Variance Table

          Df Sum Sq Mean Sq Denom F value Pr(>F) Parasite 1 0.012 0.012 157.000 0.0103 0.9192 Bacteria 1 0.058 0.058 157.000 0.0524 0.8192
> summary(fish3.lmerPQL)

Generalized linear mixed model fit using PQL Formula: dead ~ Parasite + Bacteria + (1 | Tank)

    Data: fish.data
  Family: binomial(logit link)

       AIC BIC logLik deviance
  141.3818 156.7577 -65.69091 131.3818
Random effects:

      Groups        Name    Variance    Std.Dev.
        Tank (Intercept)       5e-10  2.2361e-05
# of obs: 160, groups: Tank, 10

Estimated scale (compare to 1) 0.9318747

Fixed effects:

             Estimate Std. Error z value  Pr(>|z|)
(Intercept) -4.24380    0.79924 -5.3098 1.098e-07 ***
Parasiteyes  1.26407    0.44694  2.8283 0.0046801 **
Bacteriayes  2.85062    0.75970  3.7523 0.0001752 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
             (Intr) Prstys
Parasiteyes -0.429
Bacteriayes -0.898  0.093

Very different P-values from anova and summary.

MASS: (ver 7.2-17, R 2.1.1, Windows XP)
~~~~~


> summary(fish1.glmmpql)
Linear mixed-effects model fit by maximum likelihood Data: fish.data AIC BIC logLik 1236.652 1255.103 -612.326 Random effects: Formula: ~1 | Tank (Intercept) Residual StdDev: 0.02001341 0.8944214 Variance function: Structure: fixed weights Formula: ~invwt Fixed effects: dead ~ Parasite * Bacteria Value Std.Error DF t-value p-value (Intercept) -18.56607 1044.451 150 -0.01777591 0.9858 Parasiteyes 15.85802 1044.451 6 0.01518311 0.9884 Bacteriayes 17.23107 1044.451 6 0.01649772 0.9874 Parasiteyes:Bacteriayes -14.69007 1044.452 6 -0.01406487 0.9892 Correlation: (Intr) Prstys Bctrys Parasiteyes -1 Bacteriayes -1 1 Parasiteyes:Bacteriayes 1 -1 -1 Standardized Within-Group Residuals: Min Q1 Med Q3 Max -1.028634e+00 -5.734674e-01 -2.886770e-01 -4.224474e-14 4.330155e+00 Number of Observations: 160 Number of Groups: 10
> anova(fish1.glmmpql)
numDF denDF F-value p-value (Intercept) 1 150 17.452150 <.0001 Parasite 1 6 4.136142 0.0882 Bacteria 1 6 12.740212 0.0118 Parasite:Bacteria 1 6 0.000198 0.9892
>

> anova(glmmPQL(dead~Bacteria*Parasite, random=~1|Tank,
family=binomial, data=fish.data)) iteration 1 numDF denDF F-value p-value (Intercept) 1 150 17.452150 <.0001 Bacteria 1 6 8.980833 0.0241 Parasite 1 6 7.895521 0.0308 Bacteria:Parasite 1 6 0.000198 0.9892
>
Now anova indicates significance, but summary gives huge P-values. I have looked in MASS, ISwR, Fox's and Crawley's book, for hints, but I have probably missed the right spots/books. Hints about what to read will be also appreciated. Many thanks in advance, and sorry about the long message. The report on this research is obviously not yet published, but if the dataframe would be of help, it can be found as saved from R at: http://people.cc.jyu.fi/~aphalo/R_fish/fish.Rda. (Use load to read it into R). Pedro. -- ================================================================== Pedro J. Aphalo Department of Biological and Environmental Science University of Jyväskylä P.O. Box 35, 40351 JYVÄSKYLÄ, Finland Phone +358 14 260 2339 Mobile +358 50 3721504 Fax +358 14 260 2321 mailto:pedro.aphalo@cc.jyu.fi http://www.jyu.fi/~aphalo/ ,,,^..^,,, ______________________________________________ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Received on Thu Aug 18 19:28:06 2005

This archive was generated by hypermail 2.1.8 : Fri 03 Mar 2006 - 03:39:51 EST