GLMM doesn't do what it is supposed to (wrong results).

Jake Westfall <jake987722 at ...> writes:

Joanie,

I see that you are looking at the simple effect of herring in the presence

of a herring*age interaction. Are

the age and herring predictors centered about their means?

Jake

Date: Mon, 30 Sep 2013 09:23:28 -0400
From: Joanie.VanDeWalle at ...
To: r-sig-mixed-models at ...
Subject: [R-sig-ME] GLMM doesn't do what it is supposed to (wrong results).

Hi, 

 

My problem is that the outcome of a GLMM I built doesn't make sense.
Since I always visually analyse my dataset prior applying any
statistical analysis, I know that the outcome is wrong.  

 

Basically, I want to test the effect of the food  availability on young
seals mass gain. I captured pups each season since 10 years. For each
capture, I noted pup weight and age.  Individual seals (identified with
their flipper tag or "fliptag") were captured from 1 to 8 times
throughout their lactation period, so "individual" was put as a random
factor.  I also put the intercept and the slope as random since I expect
them to vary from one individual to the other. I also have data on the
average availability of prey for each season (1 numeric value once a
year). Attached is my dataframe.

 

[snip]
 
 Jake Westfall is exactly right.
  It would be worth reading Schielzeth 2010 _Methods Ecol. Evol._ ...

database <- read.table("vandewalle.txt", header=TRUE)
## pups captured twise or more, before age 30, pos growth
data1 <- subset(database, numberofcaptures>=2 &
               age<30 & 
               growthrate>0)
data1$fliptag <- droplevels(data1$fliptag)
length(levels(data$fliptag)) ## 317

library(nlme)

## scale data
sdata <- data.frame(lapply(data1,
                function(x) if (is.numeric(x)) c(scale(x)) else x))
plyr::numcolwise(mean)(sdata,na.rm=TRUE)  ## double-check

## model 1 (unscaled)
mod1 <- lme(weight ~ age*herring,
            random=~1+age|fliptag,
            data=data1, method="REML")

anova(mod1, type="marginal")

## model 2 (scaled)
mod2 <- update(mod1,data=sdata)
pr <- function(x) {
    printCoefmat(summary(x)$tTable,digits=3,cs.ind=1:2,has.Pvalue=TRUE)
}
pr(mod1)
pr(mod2)
anova(mod2, type="marginal")
library(plyr)
hsum <- unlist(daply(data1,"fliptag",summarise,length(unique(herring))))
table(hsum)

library(ggplot2)
theme_set(theme_bw())
ggplot(data,aes(age,weight,colour=herring))+
    stat_sum(aes(size=factor(..n..)),alpha=0.5)

ggplot(data,aes(age,herring,colour=weight))+
    stat_sum(aes(size=factor(..n..)),alpha=0.5)+
    geom_smooth()