Example of Markov Random Field with temporal interaction

mrf-time.R

# Markov Random Fields with temporal interactions

# David L Miller 2015
# Released under MIT license, YMMV


# example from ?mgcv::smooth.construct.mrf.smooth.spec
library(mgcv)
## Load Columbus Ohio crime data (see ?columbus for details and credits)
data(columb)       ## data frame
data(columb.polys) ## district shapes list
xt <- list(polys=columb.polys) ## neighbourhood structure info for MRF


## First a full rank MRF...
b <- gam(crime ~ s(district,bs="mrf",xt=xt),data=columb,method="REML")

#plot(b)

# for 10 time periods
times <- 10
dat <- c()

# build a sine wave into the response
# duplicate the data 10 times and create a time covar (tt)
for(i in 1:times){
  newdat <- columb
  newdat$crime <- newdat$crime * sin(i/2)
  newdat$tt <- i

  dat <- rbind(dat, newdat)
}


# this dies
#b <- gam(crime ~ te(district, tt, bs=c("mrf","cr"), d=c(1,1), xt=xt),
#         data=dat, method="REML")

# additive
b.add <- gam(crime ~ s(district, bs="mrf", xt=xt) + s(tt),
             data=dat, method="REML")

# with ti interaction
b.add.ti <- gam(crime ~ s(district, bs="mrf", xt=xt) + s(tt)+
                ti(district, tt, bs=c("mrf","cr"), xt=xt),
                data=dat, method="REML")

# make predictions
preds <- predict(b.add.ti)

# plot per time period
# this code is goofy 💩
par(mfrow=c(2,5))
for(i in 1:times){
  pp <- preds[1:49]
  names(pp) <- 0:48
  polys.plot(columb.polys, pp)
  preds <- preds[-c(1:49)]
}