t(mrf, time) but trying to not have a ti(time) term

mrf_ti.R

# Markov Random Fields with temporal interactions
# based on https://gist.github.com/dill/ecfe7d2f0e542bb274ff

# David L Miller 2020
# 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)
}



# full ti interaction
b.full.ti <- gam(crime ~ ti(district, bs="mrf", xt=xt) +
                         ti(tt, bs="bs")+
                         ti(district, tt, bs=c("mrf","bs"), xt=xt),
                 data=dat, method="REML")
# drop temporal smooth
# need to specify mc here to constrain tt
b.drop <- gam(crime ~ ti(district, bs="mrf", xt=xt) +
                      ti(district, tt, bs=c("mrf","bs"), xt=xt, mc=c(0,1)),
              data=dat, method="REML")
# big penalty on time instead
b.pen <- gam(crime ~ ti(district, bs="mrf", xt=xt) +
                     ti(tt, sp=1e10, bs="bs")+
                     ti(district, tt, bs=c("mrf","bs"), xt=xt),
             data=dat, method="REML")


# plot per time period
# this code is goofy 💩
plot_all <- function(model, main){
  # make predictions
  preds <- predict(model)
  for(i in 1:times){
    pp <- preds[1:49]
    names(pp) <- 0:48
    polys.plot(columb.polys, pp, lab=paste(main, i))
    preds <- preds[-c(1:49)]
  }
}

# avoid slow quartz plotting on mac
pdf(file="boop.pdf", width=12, height=4)
par(mfrow=c(3,10))

plot_all(b.full.ti, "full")
plot_all(b.drop, "drop")
plot_all(b.pen, "big penalty")

dev.off()

preds.full.ti <- predict(b.full.ti)
preds.drop <- predict(b.drop)
preds.pen <- predict(b.pen)

dev.new()
par(mfrow=c(1,3))
# these are basically the same (slightly different interaction effects)
plot(preds.full.ti, preds.drop, xlab="full", ylab="drop")
plot(preds.full.ti, preds.pen, xlab="full", ylab="pen")
plot(preds.drop, preds.pen, ylab="pen", xlab="drop")


# how do temporal smooth look?

pred_tt <- expand.grid(district=unique(dat$district),
                       tt = 1:10)

# rough base plot
par(mfrow=c(1,3))
plot(pred_tt$tt, predict(b.full.ti))
plot(pred_tt$tt, predict(b.drop))
plot(pred_tt$tt, predict(b.pen))

library(ggplot2)

pp <- rbind(pred_tt, pred_tt, pred_tt)
pp$model <- c(rep("full", nrow(pred_tt)),
              rep("drop", nrow(pred_tt)),
              rep("big penalty", nrow(pred_tt)))
pp$pred <- c(predict(b.full.ti),
             predict(b.drop),
             predict(b.pen))

ggplot(pp) +
  geom_line(aes(x=tt, y=pred, group=district, colour=district)) +
  facet_wrap(~model) +
  theme_minimal()