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:
##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)
}

 


To leave a comment for the author, please follow the link and comment on his blog: dahtah » R.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.