Using RW penalties in {mvgam}

[nicholasjclark]

Hi all, as I mentioned recently on X I think this set of tools is potentially very useful. In my package {mvgam} I have a wrapper function for easily setting up time-varying coefficient functions (https://nicholasjclark.github.io/mvgam/articles/time_varying_effects.html). This uses low-rank Gaussian Processes, either as a gp basis in {mgcv} where the user must specify the length scale, or with Hilbert space approximate GPs where the length scale is estimated, as in {brms}. These work very well, particularly the Hilbert space option, but I cannot include them easily in tensor products. This is problematic as many types of series I deal with might show subtle time-varying seasonality, and I'd like to have a tensor product of a cyclic smooth and a temporal smooth. But I don't want the temporal smooth to be a standard spline as these often don't extrapolate very well, as you all know. I think the Random Walk MRF basis might be a great option to get around this. Can I please ask whether you have plans to submit to CRAN anytime soon? And also, is my code below the correct way to implement these? It all seems to work very well, which is great to see!

# Simulate from a Gaussian Random Walk with T = 25 timesteps
y <- cumsum(rnorm(25))

# Create training data (first 20 observations)
df <- data.frame(y = y[1:20],
                 # we need a time factor variable in data for specifying the
                 # MRF basis; for this factor, the factor levels need to go 
                 # as high as we would potentially want to forecast
                 time_factor = factor(1:20, levels = as.character(1:40)))

# The penalty matrix needs to also go as high as we'd like to forecast
rw_penalty <- mrf_penalty.numeric(object = 1:40, type = 'linear')

# Fit a GAM using an MRF basis with the RW penalty; note that k can still be 
# smaller than T for a reduced rank MRF
mod <- gam(y ~ s(time_factor, bs = 'mrf', 
                 xt = list(penalty = rw_penalty),
                 k = 18),
           data = df,
           drop.unused.levels = FALSE)

summary(mod)

# Predict ahead 5 timesteps and plot against truth
preds <- predict(mod, 
                 newdata = data.frame(time_factor = as.factor(1:25)),
                 se.fit = TRUE)
plot(preds$fit, type = 'l', lwd = 2,
     ylim = c(min(preds$fit - 1.96*preds$se.fit),
              max(preds$fit + 1.96*preds$se.fit)),
     bty = 'l',
     ylab = 'Predictions',
     xlab = 'Time')
box(bty = 'l', lwd = 3)
lines(preds$fit - 1.96*preds$se.fit, col = 'darkgrey', lwd = 3)
lines(preds$fit + 1.96*preds$se.fit, col = 'darkgrey', lwd = 3)
abline(v = 20, lty = 'dashed', lwd = 2)
points(y, col = 'white', pch = 16, cex = 1)
points(y, col = 'darkred', pch = 16, cex = 0.8)

[eric-pedersen]

Hi Nick,

Glad to hear you've found this useful so far! Also, I love that example random walk model; that's exactly the kind of application we were hoping to see this used for. The forecasting behaviour is nice to see too. The code looks correct to me, although I don't think you need mrf_penalty.numeric; you should be able to use mrf_penalty(object = 1:40, type = 'linear') and it should automatically use the mrf_penalty.numeric method to create the penalty.

I still want to get this out on CRAN, but right now I don't think @gavinsimpson or I have the time it would take to finish the final steps on this to submit to CRAN in the next month or so; I hope to get back to this this coming winter, and a PhD student of mine (@firvine) has expressed interest in re-starting development on this package, so I'm cautiously optimistic that we can get a version out sooner rather than later, but I don't want to promise anything right now.

However, if you're wanting to use the code we've developed, I'm okay with you including a copy of the linear MRF penalty building code in mvgam, as long as we get credited in it. @gavinsimpson, does that work for you? We based that function off the specification of a first-order continuous-time random walk (RW1) model from Rue and Held 2006 (equation 3.26). You can view it as a Weiner process prior distribution, observed at fixed points in time.

[nicholasjclark]

Thanks very much for the reply @eric-pedersen, I'm pleased to hear you have plans for further development. Ok wonderful, I'll include a modified version of the penalty builder that only uses the linear version, and I'll make sure you all receive credit in the code and help documentation. I appreciate all the support

[gavinsimpson]

I actually don't think it would take much effort (and I'm willing to do it) to get this to the point so you can submit it to CRAN, @eric-pedersen. There's some trivial stuff in the output of R CMD check that needs fixing, but what's there works, it's documented, it just needs a little tidying.

I don't mind @nicholasjclark just taking that numeric method, but I can't imagine this would take me more than a couple of hours to get it to the point it will pass CRAN's checks. We can always polish later. Only thing really needed from you @eric-pedersen would be 10-15 mins to do the upload and submit (it requires some manual steps, like responding to an email).

How about I put my money where my mouth is and get this so it will pass CRAN's checks on all OSes and then we reevaluate?

[nicholasjclark]

Ok thanks @gavinsimpson for that update. I'll hold off then and wait for your decision about a CRAN submission. Cheers

[eric-pedersen]

Okay, let's submit it to CRAN! I had thought the process was more involved than that, but I think there's already enough here to be useful to people working with MRFs, so we might as well go ahead and publish it. There's a lingering issue from Noam regarding phylogenetic MRFs that I want to tackle before final publication, but I can tackle that this week.

[gavinsimpson]

@eric-pedersen I think Noam's Issue re adding eps was fixed in a PR that has already been merged: c116ba3

[gavinsimpson]

The package now passes CRAN checks on many OSes and R versions.

If you modify anything now @eric-pedersen be sure to check against CRAN's winbuilder system, which is very easy if you are using the devtools package: devtools::check_win_xxx() will upload to winbuilder, where xxx is the R version to check against.

[nicholasjclark]

Hi @eric-pedersen and @gavinsimpson. Happy New Year! Just checking in to see how you're going with the CRAN submission. I'd love to implement these before I submit mvgam to CRAN, but of course I'm happy to wait if there are too many other priorities right now. All the best,
Nick