3

I am missing something when trying to specifiy random effects on Gamm4.

Reproducible example:

Consider the following simulated data:

library(gamm4) 
library(ggplot2)
residualNoise <- 20 
blockNoise<- 10 
nBlocks <- 3 
n<-100 
df <- data.frame(x=rep(1:n,nBlocks),y=rep(0.01*(1:n)^2,nBlocks)) 
df$block <- rep(1:nBlocks,each=n) 
blockNoise <- rnorm(nBlocks,mean=0,sd=blockNoise) 
df$y <- df$y+blockNoise[df$block]
df$y <- df$y+rnorm(1:nrow(df),mean=0,sd=residualNoise) 
df$block <- factor(df$block)

We do a fit using gamm4 in the way that I think should be done but I get the same fit on each block:

fit1<-gamm4(y~s(x),random=~(1|block),data=df) 
df$predicted1 <- predict(fit1$gam)

ggplot(df,aes(x,y))+geom_point()+facet_wrap(~block)+
      geom_line(aes(x=x,y=predicted1),colour="red")

The following does give me different fits per block, as I imagine the random effect one should

fit2<-gamm4(y~s(x)+s(block,bs="re"),data=df)
df$predicted2 <- predict(fit2$gam)

ggplot(df,aes(x,y))+
  geom_point()+
  facet_wrap(~block)+
  geom_line(aes(x=x,y=predicted2),colour="red")

With fit1 I am getting the same fit on every block, whereas on fit2 I am getting a different fit on every block, although in both I am specifying random effects on block. I would expect different fits per block when specifying random effects on them. Why doesn't it happen on fit1?

ramiro
  • 101
  • 3

1 Answers1

1

Even though these are essentially the same model, you're seeing different behaviour because of how these models are actually fitted.

In fit1, the $gam part contains only the smooth of x, the random effects are in the mixed effects version of the model only (which is where the wiggly parts of the smooth of x are also put in that representation). When you predict from the $gam component therefore you get predictions where the random effects are set to 0.

In fit2, the random effect exists as a "smooth" in the main model matrix, not the random effects matrix. Here the model doesn't know the difference between a smooth smooth and a random effect smooth; from the point of view of the model these are all columns of "basis functions" with associated penalty matrices and coefficients. When you predict from this model, the predictions are performed using the observed values of the random variable.

To make fit2 like fit1 you can exclude the block random effect using:

predict(fit2$gam, exclude = "s(block)")

To make fit1 like fit2, you would need to predict from the $mer component and look at how to get predictions that include the random effects.

Gavin Simpson
  • 37,567
  • 5
  • 110
  • 153