Time-varying effects
Dynamic fixed-effect coefficients (often referred to as dynamic
linear models) can be readily incorporated into GAMs / DGAMs. In
mvgam, the dynamic() formula wrapper offers a
convenient interface to set these up. The plan is to incorporate a range
of dynamic options (such as random walk, AR1 etc…) but for the moment
only low-rank Gaussian Process (GP) smooths are allowed (making use
either of the gp basis in mgcv of of Hilbert
space approximate GPs). These are advantageous over splines or random
walk effects for several reasons. First, GPs will force the time-varying
effect to be smooth. This often makes sense in reality, where we would
not expect a regression coefficient to change rapidly from one time
point to the next. Second, GPs provide information on the ‘global’
dynamics of a time-varying effect through their length-scale parameters.
This means we can use them to provide accurate forecasts of how an
effect is expected to change in the future, something that we couldn’t
do well if we used splines to estimate the effect. An example below
illustrates.
Simulating time-varying effects
Simulate a time-varying coefficient using a squared exponential
Gaussian Process function with length scale \(\rho\)=10. We will do this using an
internal function from mvgam (the sim_gp
function):
set.seed(1111)
N <- 200
beta_temp <- mvgam:::sim_gp(rnorm(1),
alpha_gp = 0.75,
rho_gp = 10,
h = N) + 0.5A plot of the time-varying coefficient shows that it changes smoothly through time:
plot(beta_temp, type = 'l', lwd = 3,
bty = 'l', xlab = 'Time', ylab = 'Coefficient',
col = 'darkred')
box(bty = 'l', lwd = 2)
Next we need to simulate the values of the covariate, which we will
call temp (to represent \(temperature\)). In this case we just use a
standard normal distribution to simulate this covariate:
Finally, simulate the outcome variable, which is a Gaussian observation process (with observation error) over the time-varying effect of \(temperature\)
out <- rnorm(N, mean = 4 + beta_temp * temp,
sd = 0.25)
time <- seq_along(temp)
plot(out, type = 'l', lwd = 3,
bty = 'l', xlab = 'Time', ylab = 'Outcome',
col = 'darkred')
box(bty = 'l', lwd = 2)
Gather the data into a data.frame for fitting models,
and split the data into training and testing folds.
The dynamic() function
Time-varying coefficients can be fairly easily set up using the
s() or gp() wrapper functions in
mvgam formulae by fitting a nonlinear effect of
time and using the covariate of interest as the numeric
by variable (see ?mgcv::s or
?brms::gp for more details). The dynamic()
formula wrapper offers a way to automate this process, and will
eventually allow for a broader variety of time-varying effects (such as
random walk or AR processes). Depending on the arguments that are
specified to dynamic, it will either set up a low-rank GP
smooth function using s() with bs = 'gp' and a
fixed value of the length scale parameter \(\rho\), or it will set up a Hilbert space
approximate GP using the gp() function with
c=5/4 so that \(\rho\) is
estimated (see ?dynamic for more details). In this first
example we will use the s() option, and will mis-specify
the \(\rho\) parameter here as, in
practice, it is never known. This call to dynamic() will
set up the following smooth:
s(time, by = temp, bs = "gp", m = c(-2, 8, 2), k = 40)
mod <- mvgam(out ~ dynamic(temp, rho = 8, stationary = TRUE, k = 40),
family = gaussian(),
data = data_train,
silent = 2)Inspect the model summary, which shows how the dynamic()
wrapper was used to construct a low-rank Gaussian Process smooth
function:
summary(mod, include_betas = FALSE)
#> GAM formula:
#> out ~ s(time, by = temp, bs = "gp", m = c(-2, 8, 2), k = 40)
#> <environment: 0x0000026c107b3fd8>
#>
#> Family:
#> gaussian
#>
#> Link function:
#> identity
#>
#> Trend model:
#> None
#>
#> N series:
#> 1
#>
#> N timepoints:
#> 190
#>
#> Status:
#> Fitted using Stan
#> 4 chains, each with iter = 1000; warmup = 500; thin = 1
#> Total post-warmup draws = 2000
#>
#>
#> Observation error parameter estimates:
#> 2.5% 50% 97.5% Rhat n_eff
#> sigma_obs[1] 0.23 0.25 0.28 1 1709
#>
#> GAM coefficient (beta) estimates:
#> 2.5% 50% 97.5% Rhat n_eff
#> (Intercept) 4 4 4.1 1 1960
#>
#> Approximate significance of GAM smooths:
#> edf Ref.df Chi.sq p-value
#> s(time):temp 17 40 164 <2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Stan MCMC diagnostics:
#> n_eff / iter looks reasonable for all parameters
#> Rhat looks reasonable for all parameters
#> 0 of 2000 iterations ended with a divergence (0%)
#> 0 of 2000 iterations saturated the maximum tree depth of 12 (0%)
#> E-FMI indicated no pathological behavior
#>
#> Samples were drawn using NUTS(diag_e) at Wed Sep 04 11:52:26 AM 2024.
#> For each parameter, n_eff is a crude measure of effective sample size,
#> and Rhat is the potential scale reduction factor on split MCMC chains
#> (at convergence, Rhat = 1)Because this model used a spline with a gp basis, it’s
smooths can be visualised just like any other gam. We can
plot the estimates for the in-sample and out-of-sample periods to see
how the Gaussian Process function produces sensible smooth forecasts.
Here we supply the full dataset to the newdata argument in
plot_mvgam_smooth to inspect posterior forecasts of the
time-varying smooth function. Overlay the true simulated function to see
that the model has adequately estimated it’s dynamics in both the
training and testing data partitions
plot_mvgam_smooth(mod, smooth = 1, newdata = data)
abline(v = 190, lty = 'dashed', lwd = 2)
lines(beta_temp, lwd = 2.5, col = 'white')
lines(beta_temp, lwd = 2)
We can also use plot_predictions from the
marginaleffects package to visualise the time-varying
coefficient for what the effect would be estimated to be at different
values of \(temperature\):
require(marginaleffects)
range_round = function(x){
round(range(x, na.rm = TRUE), 2)
}
plot_predictions(mod,
newdata = datagrid(time = unique,
temp = range_round),
by = c('time', 'temp', 'temp'),
type = 'link')
This results in sensible forecasts of the observations as well
The syntax is very similar if we wish to estimate the parameters of
the underlying Gaussian Process, this time using a Hilbert space
approximation. We simply omit the rho argument in
dynamic to make this happen. This will set up a call
similar to gp(time, by = 'temp', c = 5/4, k = 40).
This model summary now contains estimates for the marginal deviation and length scale parameters of the underlying Gaussian Process function:
summary(mod, include_betas = FALSE)
#> GAM formula:
#> out ~ gp(time, by = temp, c = 5/4, k = 40, scale = TRUE)
#> <environment: 0x0000026c107b3fd8>
#>
#> Family:
#> gaussian
#>
#> Link function:
#> identity
#>
#> Trend model:
#> None
#>
#> N series:
#> 1
#>
#> N timepoints:
#> 190
#>
#> Status:
#> Fitted using Stan
#> 4 chains, each with iter = 1000; warmup = 500; thin = 1
#> Total post-warmup draws = 2000
#>
#>
#> Observation error parameter estimates:
#> 2.5% 50% 97.5% Rhat n_eff
#> sigma_obs[1] 0.24 0.26 0.3 1 2285
#>
#> GAM coefficient (beta) estimates:
#> 2.5% 50% 97.5% Rhat n_eff
#> (Intercept) 4 4 4.1 1 3056
#>
#> GAM gp term marginal deviation (alpha) and length scale (rho) estimates:
#> 2.5% 50% 97.5% Rhat n_eff
#> alpha_gp(time):temp 0.630 0.880 1.400 1.00 619
#> rho_gp(time):temp 0.026 0.053 0.069 1.01 487
#>
#> Stan MCMC diagnostics:
#> n_eff / iter looks reasonable for all parameters
#> Rhat looks reasonable for all parameters
#> 0 of 2000 iterations ended with a divergence (0%)
#> 0 of 2000 iterations saturated the maximum tree depth of 12 (0%)
#> E-FMI indicated no pathological behavior
#>
#> Samples were drawn using NUTS(diag_e) at Wed Sep 04 11:53:17 AM 2024.
#> For each parameter, n_eff is a crude measure of effective sample size,
#> and Rhat is the potential scale reduction factor on split MCMC chains
#> (at convergence, Rhat = 1)Effects for gp() terms can also be plotted as
smooths:
plot_mvgam_smooth(mod, smooth = 1, newdata = data)
abline(v = 190, lty = 'dashed', lwd = 2)
lines(beta_temp, lwd = 2.5, col = 'white')
lines(beta_temp, lwd = 2)
