My basic question is, how do I specify random effects in gamm4 properly?
Reproducible example:
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)
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")
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?
Thank you very much in advance for any insight.
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.
