Dear all,
I was wondering if anyone could help solve a problem of a missing interaction
effect!!
I carried out a 2 x 2 factorial experiment to see if eggs from 2 different
locations (Origin = 1 or 2) had different hatching success under 2 different
incubation schedules (Treat = 1 or 2). Six eggs were taken from 10 females
(random = Female) at each location and split between the treatments, giving 30
eggs from each location in each treatment.
Overall proportions hatching were as follows:
Treat
1 2
Origin
1 29/30 5/30
2 29/30 16/30
I made a binomial response in which hatching was a success and not-hatching was
a failure, and analysed as a binomial GLMM. I'm particularly interested in the
interaction between the two factors. An expression reproducing the raw data is
attached at the end of the post in case it is helpful.
hatch.frame$success<-cbind(hatch.frame$Hatched,hatch.frame$Nothatched)
model<-lmer(success~Origin*Treat+(1|Female),family=binomial,method="ML",data=hatch.frame)
model2<-update(model,~.-Origin:Treat)
anova(model,model2)
Data:
Models:
model2: success ~ Origin + Treat + (1 | Female)
model: success ~ Origin * Treat + (1 | Female)
Df AIC BIC logLik Chisq Chi Df Pr(>Chisq)
model2 4 94.707 105.857 -43.353
model 5 95.350 109.287 -42.675 1.3572 1 0.244
model3<-update(model2,~.-Origin)
anova(model2,model3)
Data:
Models:
model3: success ~ Treat + (1 | Female)
model2: success ~ Origin + Treat + (1 | Female)
Df AIC BIC logLik Chisq Chi Df Pr(>Chisq)
model3 3 98.863 107.225 -46.431
model2 4 94.707 105.857 -43.353 6.1558 1 0.01310 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
model4<-update(model2,~.-Treat)
anova(model2,model4)
Data:
Models:
model4: success ~ Origin + (1 | Female)
model2: success ~ Origin + Treat + (1 | Female)
Df AIC BIC logLik Chisq Chi Df Pr(>Chisq)
model4 3 155.592 163.954 -74.796
model2 4 94.707 105.857 -43.353 62.885 1 2.191e-15 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
So the model implies that there is a very significant effect of treatment
(reduced hatching at treatment 2) with a small effect of origin (improved
hatching from origin 2). However the lack of interaction effect implies
hatching was better for Origin 2 at both treatments, which looking at the raw
values above does not seem to be the case. Identical numbers of eggs hatched
from both Origins in Treatment 1, but much more from Origin 2 hatched at
Treatment 2.
If you divide the analysis by treatments, Origin only has a significant effect
on hatching under Treatment 2 and not with Treatment 1
Hot<-data.frame(hatch.frame[hatch.frame$Treat==2,])
Cold<-data.frame(hatch.frame[hatch.frame$Treat==1,])
#2
model<-lmer(success~Origin+(1|Female),family=binomial,method="ML",data=Hot)
model2<-update(model,~.-Origin)
anova(model,model2)
Data: Hot
Models:
model2: incubate ~ (1 | Code)
model: incubate ~ Origin + (1 | Code)
Df AIC BIC logLik Chisq Chi Df Pr(>Chisq)
model2 2 78.633 82.821 -37.316
model 3 73.697 79.980 -33.848 6.9357 1 0.008449 **
#1
model<-lmer(success~Origin+(1|Female),family=binomial,method="ML",data=Cold)
model2<-update(model,~.-Origin)
anova(model,model2)
Data: Cold
Models:
model2: incubate ~ (1 | Code)
model: incubate ~ Origin + (1 | Code)
Df AIC BIC logLik Chisq Chi Df Pr(>Chisq)
model2 2 21.5086 25.6973 -8.7543
model 3 23.3472 29.6302 -8.6736 0.1615 1 0.6878
So I can't understand where the interaction effect has gone in the full model?!
I get the same result in a binomial GLM, without the random effect of Female
i.e. a small effect of origin but no interaction with treatment. I'm sure I
must be missing something here so I would be very grateful to anyone who can
point out my mistakes. I've read previous R Help posts that suggest binomial
GLM(M) can create problems when estimated probabilities are close to 0 or 1. In
Treatment 1 hatching probability was 0.97 for both Origins, so could this be
the source of the problem?
Thanks for your help
Sam Weber
----------------------------------------------------------------------------
University of Exeter
Centre for Ecology and Conservation
Tremough Campus
Penryn
Cornwall TR10 9EZ
UK
hatch.frame <-
structure(list(Female = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 2L,
2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L,
4L, 5L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L, 6L, 7L, 7L, 7L,
7L, 7L, 7L, 8L, 8L, 8L, 8L, 8L, 8L, 9L, 9L, 9L, 9L, 9L, 9L, 10L,
10L, 10L, 10L, 10L, 10L, 11L, 11L, 11L, 11L, 11L, 11L, 12L, 12L,
12L, 12L, 12L, 12L, 13L, 13L, 13L, 13L, 13L, 13L, 14L, 14L, 14L,
14L, 14L, 14L, 15L, 15L, 15L, 15L, 15L, 15L, 16L, 16L, 16L, 16L,
16L, 16L, 17L, 17L, 17L, 17L, 17L, 17L, 18L, 18L, 18L, 18L, 18L,
18L, 19L, 19L, 19L, 19L, 19L, 19L, 20L, 20L, 20L, 20L, 20L, 20L
), .Label = c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10",
"11", "12", "13", "14", "15", "16", "17", "18", "19", "20"), class = "factor"),
Origin = structure(c(2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L,
2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L,
2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L,
1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L), .Label = c("1", "2"), class = "factor"),
Treat = structure(c(1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L,
1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L,
2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L,
1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L,
2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L,
1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L,
2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L,
1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L,
2L, 1L, 2L, 1L, 2L), .Label = c("1", "2"), class = "factor"),
Hatched = c(1L, 1L, 0L, 1L, 1L, 0L, 1L, 0L, 1L, 0L, 1L, 0L,
1L, 0L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 0L, 1L, 0L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 0L, 1L, 0L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 0L, 1L,
0L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 0L,
1L, 0L, 1L, 1L, 1L, 0L, 1L, 0L, 1L, 0L, 1L, 0L, 1L, 0L, 1L,
0L, 1L, 0L, 1L, 0L, 1L, 0L, 1L, 0L, 1L, 0L, 1L, 0L, 1L, 0L,
1L, 0L, 1L, 0L, 1L, 0L, 1L, 0L, 1L, 0L, 1L, 1L, 1L, 0L, 1L,
0L, 1L, 1L), Nothatched = c(0, 0, 1, 0, 0, 1, 0, 1, 0, 1,
0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,
0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0,
0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 0, 0, 1, 0,
0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1,
0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0,
1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0)), .Names = c("Female",
"Origin", "Treat", "Hatched", "Nothatched"), row.names = c(NA,
-120L), class = "data.frame")
[[alternative HTML version deleted]]
______________________________________________
[email protected] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
