Skip to content

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

3 messages · Van De Walle, Joanie, Jake Westfall, Ben Bolker

Original
Van De Walle, Joanie
# 
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.

 

 

database <- read.table("DataGLMMquestion.txt", header=T)

data1 <- subset(database, database$numberofcaptures>=2) # Only pups
captured twice or more

datalac <- subset(data1, data1$age<30) # I'm only interested in captured
made before age 30 days

data <- subset(datalac, datalac$growthrate>0) # Keep only the
individuals with positive growth rate

 

My model:

 

library(nlme)

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

anova(mod1, type="marginal")

 

So, the output of my model says that herring abundance has a huge effect
on individual mass.

 

However, if I look at the average of the average individual mass each
year and I plot it against the abundance of herring, there is no pattern
at all .

 

What am I doing wrong?

 

 

Thanks in advance,

 

 

Joanie Van de Walle

 

-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: DataGLMMquestion.txt
URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20130930/2aea8878/attachment-0001.txt>
Ben Bolker
# 
Jake Westfall <jake987722 at ...> writes:
of a herring*age interaction. Are
[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()