Functional ANOVA using INLA – update

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:
##Load the data
##cw = load.data()
##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')

load.data <- function()
{
cw <- with(CanadianWeather,ldply(1:length(place),function(ind) data.frame(temp=dailyAv[,ind,1],place=place[ind],region=region[ind],day=1:365)))
##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)
}