nicholasjclark
commented
10 months ago Thanks for the reprex @matt-w-rees. These are challenging to implement, and unfortunately the mgcv documentation doesn't make it easy to learn how to do them. You are correct that the by variable needs to be a matrix, but it must be a numeric matrix. If you look at the help for linear functionals (?mgcv::linear.functional.terms and ?mgcv::gam.models), you will see that the by variable is used as a weighting matrix, and the resulting smooth is constructed as rowSums(L*F), where F is your distributed lag smooth and L is the weighting matrix. I spent a lot of time on this for a few recent projects, and finally realised that I needed to supply group-level distributed lag terms for each group in the formula (this is reasonably well-described here: https://github.com/nicholasjclark/portal_VAR. Below is some code that should work for you. I'd also recommend you reinstall mvgam as the most recent version does a better job of passing along error messages from mgcv:
library(mvgam)
lagard <- function(x, n.lag = 6) {
n <- length(x); X <- matrix(NA, n, n.lag)
for (i in 1:n.lag) X[i:n, i] <- x[i:n - i + 1]
X
}
# simulate time series data at 3 sites
dat <- sim_mvgam(T = 100, n_series = 3, n_lv = 2,
family = 'poisson',
mu = runif(6, -1, 2),
trend_rel = 0.6)
# our data is
data_train <- dat$data_train
# add fake precip covariate
data_train$precip <- rnorm(nrow(data_train))
# make training data a list with lag
data_train_list <- list(lag=matrix(0:5, nrow(data_train), 6, byrow = TRUE),
y = data_train$y,
time = data_train$time,
series = data_train$series)
# calculate lagged values
data_train_list$precip <- lagard(data_train$precip)
# remove first 5 rows due to NA's (from lag)
data_train_list_trim <- lapply(data_train_list, tail, -6)
# fit SS model - same lag response across series
mod1 <- mvgam(y ~ -1,
trend_formula = ~ te(precip, lag, k = c(6, 6)),
data = data_train_list_trim,
family = poisson(),
trend_model = 'AR1')
# now try to fit SS model - different lag responses across series
# with the most recent version of {mvgam}, you will get a more informative error
# here telling you that factor 'by' variables can't be used with matrix covariates
mod2 <- mvgam(y ~ -1,
trend_formula = ~ te(precip, lag, by = trend, k = c(6, 6)),
data = data_train_list_trim,
family = poisson(),
trend_model = 'AR1')
# the best way I've found to do this, which is unforunately a bit cumbersone:
# create weight matrices for each level of the grouping facotr for setting up hierarchical
# distributed lag terms
weights_s1 <- weights_s2 <-
weights_s3 <- matrix(1, ncol = ncol(data_train_list_trim$lag),
nrow = nrow(data_train_list_trim$lag))
# for each series, the rows in its weighting matrix need to have '1s' when the corresponding
# observation in data matches that series, and '0s' otherwise
weights_s1[!(data_train_list_trim$series == 'series_1'), ] <- 0
data_train_list_trim$weights_s1 <- weights_s1
weights_s2[!(data_train_list_trim$series == 'series_2'), ] <- 0
data_train_list_trim$weights_s2 <- weights_s2
weights_s3[!(data_train_list_trim$series == 'series_3'), ] <- 0
data_train_list_trim$weights_s3 <- weights_s3
# independent series-level distributed lag terms
mod3 <- mvgam(y ~ -1,
trend_formula = ~
te(precip, lag, by = weights_s1, k = c(5, 5)) +
te(precip, lag, by = weights_s2, k = c(5, 5)) +
te(precip, lag, by = weights_s3, k = c(5, 5)),
data = data_train_list_trim,
family = poisson(),
trend_model = 'AR1')
summary(mod3)
plot(mod3, type = 'smooths', trend_effects = TRUE)
# unfortunately plot_predictions() and friends will fail because {marginaleffects}
# only currently works with data.frames, not lists (I'm working on this)
conditional_effects(mod3)
# alternatively you can set these up as a hierarchical GAM, where each series
# shares a 'global' distributed lag term and then is allowed to deviate from
# that term
mod3 <- mvgam(y ~ -1,
trend_formula = ~
# the shared distributed lag term
te(precip, lag,k = c(6, 6)) +
# series-level deviation terms
te(precip, lag, by = weights_s1, k = c(5, 5)) +
te(precip, lag, by = weights_s2, k = c(5, 5)) +
te(precip, lag, by = weights_s3, k = c(5, 5)),
data = data_train_list_trim,
family = poisson(),
trend_model = 'AR1')
matt-w-rees
Hi Nick,
I am trying to fit a multivariate model with trend-specific distributed lag effects (e.g., ' te(precip, lag, by = trend, k = c(6, 6))'). However, this appears to be failing due to internal formatting errors in mvgam. I assume this is because 'trend' needs to be in the same matrix format as the two variables ('precip', 'lag' in the previous e.g.,), but when I try to do this, the model fits but errors right at the end with "Error in PredictMat(object$smooth[[k]], data) : Can't find by variable".
I've attached a reproducible example below, any help or ideas would be greatly appreciated!
Cheers, Matt.
Created on 2023-11-16 with reprex v2.0.2