Functional ANOVA using INLA – update

[This article was first published on dahtah » R, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

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.

##Load the data
##cw =
##Run INLA
##res <- inla.fanova.temperature(cw)
##Plot "fitted" model


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

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

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(,model="rw2",cyclic=T,replicate=place.ind, diagonal=1e-5) +
##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
control.fixed=list(prec=0.01, prec.intercept = 0.01),verbose=T)
} <- 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)

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,
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)

##Extract "fixed effects"
feff <- function(res.inla)


To leave a comment for the author, please follow the link and comment on their blog: dahtah » R. offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)