[R-sig-ME] GLMM doesn't do what it is supposed to (wrong results).

[Ben Bolker]

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()