# NOT RUN {
## NOTE: most examples are flagged as 'do not run' simply to
## save time in package checking on CRAN.
###################################
## A simple additive mixed model...
###################################
library(gamm4)
set.seed(0)
dat <- gamSim(1,n=400,scale=2) ## simulate 4 term additive truth
## Now add 20 level random effect `fac'...
dat$fac <- fac <- as.factor(sample(1:20,400,replace=TRUE))
dat$y <- dat$y + model.matrix(~fac-1)%*%rnorm(20)*.5
br <- gamm4(y~s(x0)+x1+s(x2),data=dat,random=~(1|fac))
plot(br$gam,pages=1)
summary(br$gam) ## summary of gam
summary(br$mer) ## underlying mixed model
anova(br$gam)
## compare gam fit of the same
bg <- gam(y~s(x0)+x1+s(x2)+s(fac,bs="re"),
data=dat,method="REML")
plot(bg,pages=1)
gam.vcomp(bg)
##########################
## Poisson example GAMM...
##########################
## simulate data...
x <- runif(100)
fac <- sample(1:20,100,replace=TRUE)
eta <- x^2*3 + fac/20; fac <- as.factor(fac)
y <- rpois(100,exp(eta))
## fit model and examine it...
bp <- gamm4(y~s(x),family=poisson,random=~(1|fac))
plot(bp$gam)
bp$mer
# }
# NOT RUN {
#################################################################
## Add a factor to the linear predictor, to be modelled as random
## and make response Poisson. Again compare `gamm' and `gamm4'
#################################################################
set.seed(6)
dat <- gamSim(1,n=400,scale=2) ## simulate 4 term additive truth
## add random effect...
g <- as.factor(sample(1:20,400,replace=TRUE))
dat$f <- dat$f + model.matrix(~ g-1)%*%rnorm(20)*2
dat$y <- rpois(400,exp(dat$f/7+1))
b2<-gamm(y~s(x0)+s(x1)+s(x2)+s(x3),family=poisson,
data=dat,random=list(g=~1))
plot(b2$gam,pages=1)
b2r<-gamm4(y~s(x0)+s(x1)+s(x2)+s(x3),family=poisson,
data=dat,random = ~ (1|g))
plot(b2r$gam,pages=1)
rm(dat)
vis.gam(b2r$gam,theta=35)
##################################
# Multivariate varying coefficient
# With crossed and nested random
# effects.
##################################
## Start by simulating data...
f0 <- function(x, z, sx = 0.3, sz = 0.4) {
(pi^sx * sz) * (1.2 * exp(-(x - 0.2)^2/sx^2 - (z -
0.3)^2/sz^2) + 0.8 * exp(-(x - 0.7)^2/sx^2 -
(z - 0.8)^2/sz^2))
}
f1 <- function(x2) 2 * sin(pi * x2)
f2 <- function(x2) exp(2 * x2) - 3.75887
f3 <- function (x2) 0.2 * x2^11 * (10 * (1 - x2))^6 + 10 * (10 * x2)^3 *
(1 - x2)^10
n <- 1000
## first set up a continuous-within-group effect...
g <- factor(sample(1:50,n,replace=TRUE)) ## grouping factor
x <- runif(n) ## continuous covariate
X <- model.matrix(~g-1)
mu <- X%*%rnorm(50)*.5 + (x*X)%*%rnorm(50)
## now add nested factors...
a <- factor(rep(1:20,rep(50,20)))
b <- factor(rep(rep(1:25,rep(2,25)),rep(20,50)))
Xa <- model.matrix(~a-1)
Xb <- model.matrix(~a/b-a-1)
mu <- mu + Xa%*%rnorm(20) + Xb%*%rnorm(500)*.5
## finally simulate the smooth terms
v <- runif(n);w <- runif(n);z <- runif(n)
r <- runif(n)
mu <- mu + f0(v,w)*z*10 + f3(r)
y <- mu + rnorm(n)*2 ## response data
## First compare gamm and gamm4 on a reduced model
br <- gamm4(y ~ s(v,w,by=z) + s(r,k=20,bs="cr"),random = ~ (1|a/b))
ba <- gamm(y ~ s(v,w,by=z) + s(r,k=20,bs="cr"),random = list(a=~1,b=~1),method="REML")
par(mfrow=c(2,2))
plot(br$gam)
plot(ba$gam)
## now fit the full model
br <- gamm4(y ~ s(v,w,by=z) + s(r,k=20,bs="cr"),random = ~ (x+0|g) + (1|g) + (1|a/b))
br$mer
br$gam
plot(br$gam)
## try a Poisson example, based on the same linear predictor...
lp <- mu/5
y <- rpois(exp(lp),exp(lp)) ## simulated response
## again compare gamm and gamm4 on reduced model
br <- gamm4(y ~ s(v,w,by=z) + s(r,k=20,bs="cr"),family=poisson,random = ~ (1|a/b))
ba <- gamm(y ~ s(v,w,by=z) + s(r,k=20,bs="cr"),family=poisson,random = list(a=~1,b=~1))
par(mfrow=c(2,2))
plot(br$gam)
plot(ba$gam)
## and now fit full version (very slow)...
br <- gamm4(y ~ s(v,w,by=z) + s(r,k=20,bs="cr"),family=poisson,random = ~ (x|g) + (1|a/b))
br$mer
br$gam
plot(br$gam)
####################################
# Different smooths of x2 depending
# on factor `fac'...
####################################
dat <- gamSim(4)
br <- gamm4(y ~ fac+s(x2,by=fac)+s(x0),data=dat)
plot(br$gam,pages=1)
summary(br$gam)
####################################
# Timing comparison with `gam'... #
####################################
dat <- gamSim(1,n=600,dist="binary",scale=.33)
system.time(lr.fit0 <- gam(y~s(x0)+s(x1)+s(x2),
family=binomial,data=dat,method="ML"))
system.time(lr.fit <- gamm4(y~s(x0)+s(x1)+s(x2),
family=binomial,data=dat))
lr.fit0;lr.fit$gam
cor(fitted(lr.fit0),fitted(lr.fit$gam))
## plot model components with truth overlaid in red
op <- par(mfrow=c(2,2))
fn <- c("f0","f1","f2","f3");xn <- c("x0","x1","x2","x3")
for (k in 1:3) {
plot(lr.fit$gam,select=k)
ff <- dat[[fn[k]]];xx <- dat[[xn[k]]]
ind <- sort.int(xx,index.return=TRUE)$ix
lines(xx[ind],(ff-mean(ff))[ind]*.33,col=2)
}
par(op)
# }
# NOT RUN {
######################################
## A "signal" regression example, in
## which a univariate response depends
## on functional predictors.
######################################
## simulate data first....
rf <- function(x=seq(0,1,length=100)) {
## generates random functions...
m <- ceiling(runif(1)*5) ## number of components
f <- x*0;
mu <- runif(m,min(x),max(x));sig <- (runif(m)+.5)*(max(x)-min(x))/10
for (i in 1:m) f <- f+ dnorm(x,mu[i],sig[i])
f
}
x <- seq(0,1,length=100) ## evaluation points
## example functional predictors...
par(mfrow=c(3,3));for (i in 1:9) plot(x,rf(x),type="l",xlab="x")
## simulate 200 functions and store in rows of L...
L <- matrix(NA,200,100)
for (i in 1:200) L[i,] <- rf() ## simulate the functional predictors
f2 <- function(x) { ## the coefficient function
(0.2*x^11*(10*(1-x))^6+10*(10*x)^3*(1-x)^10)/10
}
f <- f2(x) ## the true coefficient function
y <- L%*%f + rnorm(200)*20 ## simulated response data
## Now fit the model E(y) = L%*%f(x) where f is a smooth function.
## The summation convention is used to evaluate smooth at each value
## in matrix X to get matrix F, say. Then rowSum(L*F) gives E(y).
## create matrix of eval points for each function. Note that
## `smoothCon' is smart and will recognize the duplication...
X <- matrix(x,200,100,byrow=TRUE)
## compare `gam' and `gamm4' this time
b <- gam(y~s(X,by=L,k=20),method="REML")
br <- gamm4(y~s(X,by=L,k=20))
par(mfrow=c(2,1))
plot(b,shade=TRUE);lines(x,f,col=2)
plot(br$gam,shade=TRUE);lines(x,f,col=2)
# }
Run the code above in your browser using DataLab