# Functional ANOVA using INLA – update

February 29, 2012
By

(This article was first published on dahtah » R, and kindly contributed to R-bloggers)

INLA author Håvard Rue wrote me to point out a problem in the Functional ANOVA code given in this post. I made a mistake in setting the precision of the fixed effects (I used “default” instead of “prec”). I’ve put Håvard’s corrected version of the code below.


##Usage:
##Run INLA
##res <- inla.fanova.temperature(cw)
##
##Plot "fitted" model
##plot.fitted(cw,res)
##
##etc.

require(fda)
require(INLA)
require(ggplot2)
require(stringr)
require(mgcv)
require(directlabels)

##Use the same colour scheme Ramsay&Silverman did
col.scale <- c('red','blue','darkgreen','black')

{
##We need to create multiple copies of the time index because we need multiple functions of time
cw <- mutate(cw,day.region=day,day.place=day,
region.ind=as.numeric(region),
place.ind=as.numeric(place))

##INLA apparently doesn't like the original factor levels, we modify them
levels(cw$place) <- str_replace(levels(cw$place),'. ','_')
cw
}

inla.fanova.temperature <- function(cw)
{
##The formula for the model
formula <- temp ~ f(day.region,model="rw2",replicate=region.ind,cyclic=T, diagonal=1e-5)+
f(day,model="rw2",cyclic=T, diagonal=1e-5)+
f(day.place,model="rw2",cyclic=T,replicate=place.ind, diagonal=1e-5) +
region
##Note that 'region' comes in as a factor, and inla treats factors
##in the same way as the 'lm' function, i.e., using contrasts This
##is not strictly necessary in a Bayesian analysis and here
##complicates things a bit The default "treatment" contrasts used
##here mean that the "place" factor has the Pacific region as the
##default level to which other regions are compared

##Call inla We use control.fixed to impose a proper gaussian prior
##on the fixed effects and control.predictor to make INLA compute
##marginal distributions for each value of the linear predictor
##The call takes 160 sec. on my machine
inla(formula,data=cw,family="gaussian",
control.predictor=list(compute=T),
control.fixed=list(prec=0.01, prec.intercept = 0.01),verbose=T)
}

plot.raw.data <- function(cw)
{
p <- ggplot(cw,aes(day,temp,group=place,colour=region))+geom_point(alpha=.1)+geom_smooth(method="gam",form = y ~ s(x))+labs(x='Day of the year',y='Temp. (C)')+scale_colour_manual(values=col.scale)+facet_wrap(~ region)
p + theme_bw() + opts(legend.position="none")
}

plot.fitted <- function(cw,res)
{
cw$fitted <- (res$summary.fitted.values$mean) p <- ggplot(cw,aes(day,temp,group=place,colour=region))+geom_point(alpha=.1)+facet_wrap(~ region)+scale_colour_manual(values=col.scale)+geom_path(aes(y=fitted)) p + theme_bw() + opts(legend.position="none") } extract.regional.effects <- function(cw,res) { reg.effect <- reff(res)$day.reg
names(reg.effect)[1] <- 'day'
reg.effect$region <- gl(4,365,lab=levels(cw$region))

##Total regional effect is equal to smooth component + intercept + regional effect
feff.region <- feff(res)[str_detect(attributes(feff(res))$row.names,"region"),] feff.inter <- feff(res)[str_detect(attributes(feff(res))$row.names,"Inter"),]
offsets <- c(feff.inter$mean,feff.inter$mean+feff.region$mean) reg.effect$total.effect <- reg.effect$mean+rep(offsets,each=365) reg.effect } plot.regional.effects <- function(cw,res,smooth=F) { reg.effect <- extract.regional.effects(cw,res) p <- ggplot(reg.effect,aes(day,total.effect,colour=region))+geom_line()+scale_colour_manual(values=c('red','blue','darkgreen','black'))+scale_y_continuous(lim=c(-18,15)) if (smooth) p <- p + geom_smooth(method="gam",form = y ~ s(x),lty=2) p <- p + geom_abline(slope=0,intercept=0,lty=3,col="lightblue") p <- p + geom_dl(aes(label=region),method="top.qp") p + theme_bw() +opts(panel.grid.minor = theme_blank(),panel.grid.major = theme_blank(),legend.position="none") + labs(x='\n Day of the year',y='Temp. (C)') } plot.regional.avg <- function(cw,res) { reg.effect <- extract.regional.effects(cw,res) regionavg <- data.frame(day=1:365, glob.effect=reff(res)$day$mean, reg.avg=reff(res)$day$mean+reg.effect$total.effect,
region=reg.effect$region ) p <- ggplot(regionavg,aes(day,reg.avg,colour=region))+geom_line()+facet_wrap(~ region)+scale_colour_manual(values=col.scale)+geom_line(aes(y=glob.effect),col="darkgrey",lty=2) p <- p + geom_abline(slope=0,intercept=0,lty=3,col="lightblue") p + theme_bw() +opts(panel.grid.minor = theme_blank(),panel.grid.major = theme_blank(),legend.position="none") + labs(x='\n Day of the year',y='Temp. (C)') } ##Extract "random effects" summary from an inla object reff <- function(res.inla) { res.inla$summary.random
}

##Extract "fixed effects"
feff <- function(res.inla)
{
as.data.frame(res.inla\$summary.fixed)
}