From: Pedro J. Aphalo <pedro.aphalo_at_cc.jyu.fi>

Date: Thu 18 Aug 2005 - 19:20:02 EST

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) ~~~~~Received on Thu Aug 18 19:28:06 2005

> 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

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