📈 fitting categorical covariates as 1D Markov random fields

is_this_a_good_idea.R

# example using an MRF for ordered categorical predictors

library(mgcv)


# simulate some data...
set.seed(2)
dat <- gamSim(1,n=400,dist="normal",scale=2)

# fit a simple model
b <- gam(y~s(x0)+s(x1)+s(x2)+s(x3),data=dat, method="REML")
summary(b)
plot(b,pages=1)

# okay now let's pretend we have ordinal data
# we'll just chop-up the x0 covariate and pretend that's the way things are
dat$x0c <- cut(dat$x0, seq(0,1,by=0.2), labels=1:5)
table(dat$x0c)

# okay how can we fit this model?
## this bit is the bit that's important for you the bit above is just faff

# first setup the Markov random field
# see the ?mrf manual page

# we need to setup a neighbourhood structure, so we have a list of the
# categories and then say in each element what it's neighbours are

nb <- list()
nb[[1]] <- c(2)
nb[[2]] <- c(1,3)
nb[[3]] <- c(2,4)
nb[[4]] <- c(3,5)
nb[[5]] <- c(4)

names(nb) <- as.character(1:5)

b1 <- gam(y~s(x0c, xt=list(nb=nb), bs="mrf", k=5)+s(x1)+s(x2)+s(x3),data=dat, method="REML")
summary(b1)

# compare to the original smooth version
# this might be just hieroglyphics, but it's just for our purposes to test here
plot(b, select=1)

pred_grid <- data.frame(x0=seq(0.1,1,by=0.2), x0c=1:5, x1=0, x2=0, x3=0)
pred_grid$p <- predict(b1, pred_grid, type="iterms")[,1]
points(pred_grid$x0, pred_grid$p)

# not a bad approximation!!