diff --git a/R/dynamic.R b/R/dynamic.R index c854af8e..a2c29222 100644 --- a/R/dynamic.R +++ b/R/dynamic.R @@ -135,6 +135,9 @@ interpret_mvgam = function(formula, N){ facs <- colnames(attr(terms.formula(formula), 'factors')) + # Check if formula has an intercept + keep_intercept <- attr(terms(formula), 'intercept') == 1 + # Re-arrange so that random effects always come last if(any(grepl('bs = \"re\"', facs, fixed = TRUE))){ newfacs <- facs[!grepl('bs = \"re\"', facs, fixed = TRUE)] @@ -263,5 +266,10 @@ interpret_mvgam = function(formula, N){ updated_formula <- newformula } + if(!keep_intercept){ + updated_formula <- update(updated_formula, . ~ . - 1) + attr(updated_formula, '.Environment') <- attr(newformula, '.Environment') + } + return(updated_formula) } diff --git a/R/gp.R b/R/gp.R index 974aa9b3..ae2b1637 100644 --- a/R/gp.R +++ b/R/gp.R @@ -93,6 +93,7 @@ make_gp_additions = function(gp_details, data, coefs_replace[[x]] <- which_replace } + # Replace basis functions with gp() eigenfunctions newX <- model_data$X # Training data eigenfunctions @@ -155,6 +156,22 @@ make_gp_additions = function(gp_details, data, # Add the GP attribute table to the mgcv_model attr(mgcv_model, 'gp_att_table') <- gp_att_table + # Assign GP labels to smooths + gp_assign <- data.frame(label = unlist(purrr::map(gp_att_table, 'name')), + first.para = unlist(purrr::map(gp_att_table, 'first_coef')), + last.para = unlist(purrr::map(gp_att_table, 'last_coef')), + by = unlist(purrr::map(gp_att_table, 'by'))) + for(i in seq_along(mgcv_model$smooth)){ + if(mgcv_model$smooth[[i]]$label %in% + gsub('gp(', 's(', gp_assign$label, fixed = TRUE) & + mgcv_model$smooth[[i]]$first.para %in% gp_assign$first.para){ + mgcv_model$smooth[[i]]$gp_term <- TRUE + } else { + mgcv_model$smooth[[i]]$gp_term <- FALSE + } + } + + # Return return(list(model_data = model_data, mgcv_model = mgcv_model, gp_stan_lines = gp_stan_lines, @@ -303,7 +320,7 @@ sim_hilbert_gp = function(alpha_gp, #' Mean-center and scale the particular covariate of interest #' so that the maximum Euclidean distance between any two points is 1 #' @noRd -scale_cov <- function(data, covariate, by, level, scale = TRUE, +scale_cov <- function(data, covariate, by, level, mean, max_dist){ Xgp <- data[[covariate]] if(!is.na(by) & @@ -311,26 +328,29 @@ scale_cov <- function(data, covariate, by, level, scale = TRUE, Xgp <- data[[covariate]][data[[by]] == level] } - if(is.na(mean)){ - Xgp_mean <- mean(Xgp, na.rm = TRUE) - } else { - Xgp_mean <- mean - } - + # Compute max Euclidean distance if not supplied if(is.na(max_dist)){ Xgp_max_dist <- sqrt(max(brms:::diff_quad(Xgp))) } else { Xgp_max_dist <- max_dist } - if(scale){ - # Mean center and divide by max euclidean distance - (Xgp - Xgp_mean) / Xgp_max_dist + # Scale + Xgp <- Xgp / Xgp_max_dist + # Compute mean if not supplied (after scaling) + if(is.na(mean)){ + Xgp_mean <- mean(Xgp, na.rm = TRUE) } else { - # Just mean center - Xgp - Xgp_mean + Xgp_mean <- mean } + + # Center + Xgp <- Xgp - Xgp_mean + + return(list(Xgp = Xgp, + Xgp_mean = Xgp_mean, + Xgp_max_dist = Xgp_max_dist)) } #' prep GP eigenfunctions @@ -355,8 +375,7 @@ prep_eigenfunctions = function(data, by = by, level = level, mean = mean, - max_dist = max_dist, - scale = scale) + max_dist = max_dist)$Xgp # Construct matrix of eigenfunctions eigenfunctions <- matrix(NA, nrow = length(covariate_cent), @@ -431,25 +450,24 @@ prep_gp_covariate = function(data, if(def_alpha == ''){ def_alpha<- 'student_t(3, 0, 2.5);' } + # Prepare the covariate + if(scale){ + max_dist <- NA + } else { + max_dist <- 1 + } + covariate_cent <- scale_cov(data = data, covariate = covariate, - scale = scale, by = by, mean = NA, - max_dist = NA, + max_dist = max_dist, level = level) - Xgp <- data[[covariate]] - if(!is.na(by) & - !is.na(level)){ - Xgp <- data[[covariate]][data[[by]] == level] - } - - covariate_mean <- mean(Xgp, na.rm = TRUE) - covariate_max_dist <- ifelse(scale, - sqrt(max(brms:::diff_quad(Xgp))), - 1) + covariate_mean <- covariate_cent$Xgp_mean + covariate_max_dist <- covariate_cent$Xgp_max_dist + covariate_cent <- covariate_cent$Xgp # Construct vector of eigenvalues for GP covariance matrix; the # same eigenvalues are always used in prediction, so we only need to @@ -469,9 +487,10 @@ prep_gp_covariate = function(data, covariate = covariate, by = by, level = level, + L = L, k = k, boundary = boundary, - mean = NA, + mean = covariate_mean, max_dist = covariate_max_dist, scale = scale, initial_setup = TRUE) diff --git a/R/mvgam.R b/R/mvgam.R index 8b46e6a1..6dae53fb 100644 --- a/R/mvgam.R +++ b/R/mvgam.R @@ -652,7 +652,7 @@ mvgam = function(formula, # in the model.frame formula <- gp_to_s(formula) if(!keep_intercept){ - formula <- update(formula, trend_y ~ . -1) + formula <- update(formula, . ~ . - 1) } } diff --git a/R/plot.mvgam.R b/R/plot.mvgam.R index c57e9e5f..51828d81 100644 --- a/R/plot.mvgam.R +++ b/R/plot.mvgam.R @@ -134,12 +134,6 @@ plot.mvgam = function(x, type = 'residuals', mgcv_plottable = object2$mgcv_model$smooth[[x]]$plot.me) })) - # Filter out any GP terms - if(!is.null(attr(object2$mgcv_model, 'gp_att_table'))){ - gp_names <- unlist(purrr::map(attr(object2$mgcv_model, 'gp_att_table'), 'name')) - smooth_labs %>% - dplyr::filter(!label %in% gsub('gp(', 's(', gp_names, fixed = TRUE)) -> smooth_labs - } n_smooths <- NROW(smooth_labs) if(n_smooths == 0) stop("No smooth terms to plot. Use plot_predictions() to visualise other effects", call. = FALSE) diff --git a/R/plot_mvgam_smooth.R b/R/plot_mvgam_smooth.R index d983f631..db2b2aa5 100644 --- a/R/plot_mvgam_smooth.R +++ b/R/plot_mvgam_smooth.R @@ -115,17 +115,6 @@ plot_mvgam_smooth = function(object, smooth_int <- smooth } - # Check whether this is actually a gp() term - if(!is.null(attr(object2$mgcv_model, 'gp_att_table'))){ - gp_names <- unlist(purrr::map(attr(object2$mgcv_model, 'gp_att_table'), 'name')) - if(any(grepl(object2$mgcv_model$smooth[[smooth_int]]$label, - gsub('gp(', 's(', gp_names, fixed = TRUE), - fixed = TRUE))){ - stop(smooth, ' is a gp() term. Use plot_predictions() instead to visualise', - call. = FALSE) - } - } - # Check whether this type of smooth is even plottable if(!object2$mgcv_model$smooth[[smooth_int]]$plot.me){ stop(paste0('unable to plot ', object2$mgcv_model$smooth[[smooth_int]]$label, @@ -290,14 +279,61 @@ plot_mvgam_smooth = function(object, # If this term has a by variable, need to use mgcv's plotting utilities if(object2$mgcv_model$smooth[[smooth_int]]$by != "NA"){ - # Deal with by variables - by <- rep(1,length(pred_vals)); dat <- data.frame(x = pred_vals, by = by) - names(dat) <- c(object2$mgcv_model$smooth[[smooth_int]]$term, - object2$mgcv_model$smooth[[smooth_int]]$by) + # Check if this is a gp() term + gp_term <- FALSE + if(!is.null(attr(object2$mgcv_model, 'gp_att_table'))){ + gp_term <- object2$mgcv_model$smooth[[smooth_int]]$gp_term + } + + if(gp_term){ + object2$mgcv_model$smooth[[smooth_int]]$label <- + gsub('s(', 'gp(', + object2$mgcv_model$smooth[[smooth_int]]$label, + fixed = TRUE) + # Check if this is a factor by variable + is_fac <- is.factor(object2$obs_data[[object2$mgcv_model$smooth[[smooth_int]]$by]]) + + if(is_fac){ + fac_levels <- levels(object2$obs_data[[object2$mgcv_model$smooth[[smooth_int]]$by]]) + whichlevel <- vector() + for(i in seq_along(fac_levels)){ + whichlevel[i] <- grepl(fac_levels[i], object2$mgcv_model$smooth[[smooth_int]]$label, + fixed = TRUE) + } + + pred_dat[[object2$mgcv_model$smooth[[smooth_int]]$by]] <- + rep(fac_levels[whichlevel], length(pred_dat$series)) + } + + if(!is_fac){ + pred_dat[[object2$mgcv_model$smooth[[smooth_int]]$by]] <- + rep(1, length(pred_dat$series)) + } + + if(trend_effects){ + Xp_term <- trend_Xp_matrix(newdata = pred_dat, + trend_map = object2$trend_map, + mgcv_model = object2$trend_mgcv_model) + } else { + Xp_term <- obs_Xp_matrix(newdata = pred_dat, + mgcv_model = object2$mgcv_model) + } + Xp[,object2$mgcv_model$smooth[[smooth_int]]$first.para: + object2$mgcv_model$smooth[[smooth_int]]$last.para] <- + Xp_term[,object2$mgcv_model$smooth[[smooth_int]]$first.para: + object2$mgcv_model$smooth[[smooth_int]]$last.para] + + } else { + # Deal with by variables in non-gp() smooths + by <- rep(1,length(pred_vals)); dat <- data.frame(x = pred_vals, by = by) + names(dat) <- c(object2$mgcv_model$smooth[[smooth_int]]$term, + object2$mgcv_model$smooth[[smooth_int]]$by) + + Xp_term <- mgcv::PredictMat(object2$mgcv_model$smooth[[smooth_int]], dat) + Xp[,object2$mgcv_model$smooth[[smooth_int]]$first.para: + object2$mgcv_model$smooth[[smooth_int]]$last.para] <- Xp_term + } - Xp_term <- mgcv::PredictMat(object2$mgcv_model$smooth[[smooth_int]], dat) - Xp[,object2$mgcv_model$smooth[[smooth_int]]$first.para: - object2$mgcv_model$smooth[[smooth_int]]$last.para] <- Xp_term } # Extract GAM coefficients diff --git a/docs/articles/SS_model.svg b/docs/articles/SS_model.svg index 415097bf..f6441719 100644 --- a/docs/articles/SS_model.svg +++ b/docs/articles/SS_model.svg @@ -5,12 +5,36 @@ xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#" xmlns:svg="http://www.w3.org/2000/svg" xmlns="http://www.w3.org/2000/svg" + xmlns:sodipodi="http://sodipodi.sourceforge.net/DTD/sodipodi-0.dtd" + xmlns:inkscape="http://www.inkscape.org/namespaces/inkscape" style="overflow:hidden" id="svg81" version="1.1" overflow="hidden" height="324.50705" - width="945.35211"> + width="945.35211" + sodipodi:docname="SS_model.svg" + inkscape:version="0.92.4 (5da689c313, 2019-01-14)"> + @@ -19,7 +43,7 @@ image/svg+xml - + @@ -45,8 +69,9 @@ d="m 850.01408,120.495 0.175,53.401 -3,0.01 -0.175,-53.401 z m 3.17,51.891 -4.47,9.015 -4.529,-8.986 z" id="path13" /> + d="m 753.51408,71 34.1,9e-5 v 3 l -34.1,-9e-5 z m 32.6,-2.99992 9,4.500025 -9,4.499975 z" + id="path15" + inkscape:connector-curvature="0" /> diff --git a/docs/articles/time_varying_effects.html b/docs/articles/time_varying_effects.html index 96704bdf..36fdce43 100644 --- a/docs/articles/time_varying_effects.html +++ b/docs/articles/time_varying_effects.html @@ -76,7 +76,7 @@

Nicholas J Clark

-

2023-10-15

+

2023-10-16

Source: vignettes/time_varying_effects.Rmd
time_varying_effects.Rmd
@@ -161,7 +161,7 @@

Simulating time-varying effectsThe dynamic() function

Time-varying coefficients can be fairly easily set up using the -s() or gp() wrapper functions in +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 @@ -170,18 +170,18 @@

The dynamic() functions() with bs = 'gp' and a +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 +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 +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 = 20)

-mod <- mvgam(out ~ dynamic(temp, rho = 8, stationary = TRUE, k = 20),
+mod <- mvgam(out ~ dynamic(temp, rho = 8, stationary = TRUE, k = 40),
              family = gaussian(),
              data = data_train)

Inspect the model summary, which shows how the dynamic() @@ -190,7 +190,7 @@

The dynamic() function
 summary(mod, include_betas = FALSE)
## GAM formula:
-## out ~ s(time, by = temp, bs = "gp", m = c(-2, 8, 2), k = 20)
+## out ~ s(time, by = temp, bs = "gp", m = c(-2, 8, 2), k = 40)
 ## 
 ## Family:
 ## gaussian
@@ -212,11 +212,11 @@ 
The dynamic() function## 
 ## Observation error parameter estimates:
 ##              2.5%  50% 97.5% Rhat n.eff
-## sigma_obs[1] 0.23 0.26  0.29    1  2146
+## sigma_obs[1] 0.23 0.25  0.28    1  2773
 ## 
 ## GAM coefficient (beta) estimates:
 ##             2.5% 50% 97.5% Rhat n.eff
-## (Intercept)    4   4   4.1    1  2852
+## (Intercept)    4   4   4.1    1  2558
 ## 
 ## Stan MCMC diagnostics:
 ## n_eff / iter looks reasonable for all parameters
@@ -249,8 +249,8 @@ 
The dynamic() function\(temperature\):

 
-plot_predictions(mod, 
-                 newdata = datagrid(time = unique,
+plot_predictions(mod, 
+                 newdata = datagrid(time = unique,
                                     temp = fivenum),
                  by = c('time', 'temp', 'temp'),
                  type = 'link')

@@ -260,14 +260,14 @@ 
The dynamic() functionplot(mod, type = 'forecast', newdata = data_test)

## Out of sample CRPS:
-## [1] 1.453761
+## [1] 1.298432

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 = 20.

+similar to gp(time, by = 'temp', c = 5/4), k = 20).

-mod <- mvgam(out ~ dynamic(temp, k = 20),
+mod <- mvgam(out ~ dynamic(temp, k = 40),
              family = gaussian(),
              data = data_train)

This model summary now contains estimates for the marginal deviation @@ -276,7 +276,7 @@

The dynamic() function
 summary(mod, include_betas = FALSE)
## GAM formula:
-## out ~ gp(time, by = temp, c = 5/4, k = 20, scale = TRUE)
+## out ~ gp(time, by = temp, c = 5/4, k = 40, scale = TRUE)
 ## 
 ## Family:
 ## gaussian
@@ -298,43 +298,52 @@ 
The dynamic() function## 
 ## Observation error parameter estimates:
 ##              2.5%  50% 97.5% Rhat n.eff
-## sigma_obs[1] 0.33 0.37  0.41    1  1730
+## sigma_obs[1] 0.24 0.26   0.3    1  2254
 ## 
 ## GAM coefficient (beta) estimates:
 ##             2.5% 50% 97.5% Rhat n.eff
-## (Intercept)    4 4.1   4.1    1  2941
+## (Intercept)    4   4   4.1    1  3513
 ## 
 ## GAM gp term marginal deviation (alpha) and length scale (rho) estimates:
-##                       2.5%   50% 97.5% Rhat n.eff
-## alpha_gp(time):temp 1.2000 2.100 4.400 1.01   374
-## rho_gp(time):temp   0.0095 0.027 0.067 1.00  1471
+##                      2.5%   50% 97.5% Rhat n.eff
+## alpha_gp(time):temp 0.630 0.880 1.400    1   710
+## rho_gp(time):temp   0.029 0.053 0.069    1   629
 ## 
 ## Approximate significance of GAM observation smooths:
-##               edf    F p-value
-## s(time):temp 1.98 3605    0.59
+##               edf   F p-value
+## s(time):temp 2.33 699     0.9
 ## 
 ## Stan MCMC diagnostics:
 ## n_eff / iter looks reasonable for all parameters
 ## Rhat looks reasonable for all parameters
-## 3 of 2000 iterations ended with a divergence (0.15%)
+## 1 of 2000 iterations ended with a divergence (0.05%)
 ## *Try running with larger adapt_delta to remove the divergences
 ## 0 of 2000 iterations saturated the maximum tree depth of 12 (0%)
 ## E-FMI indicated no pathological behavior
-

The plot_predictions call shows that the effect in this -case is similar to what we estimated above:

+

Effects for gp() terms can also be plotted as +smooths:

-plot_predictions(mod, 
-                 newdata = datagrid(time = unique,
+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)
+

+

Both the above plot and the below plot_predictions() +call show that the effect in this case is similar to what we estimated +in the approximate GP smooth model above:

+
+plot_predictions(mod, 
+                 newdata = datagrid(time = unique,
                                     temp = fivenum),
                  by = c('time', 'temp', 'temp'),
                  type = 'link')
-

+

Forecasts are also similar:

-
+
 plot(mod, type = 'forecast', newdata = data_test)

-

+

 
## Out of sample CRPS:
-## [1] 6.21817

+## [1] 1.677305
@@ -352,7 +361,7 @@

Salmon survival exampleMARSS package:

-
+
 load(url('https://github.com/atsa-es/MARSS/raw/master/data/SalmonSurvCUI.rda'))
 dplyr::glimpse(SalmonSurvCUI)

 
## Rows: 42
@@ -368,7 +377,7 @@ 
Salmon survival exampletime indicator and a series
 indicator for working in mvgam:

-
+
 SalmonSurvCUI %>%
   # create a time variable
   dplyr::mutate(time = dplyr::row_number()) %>%
@@ -382,7 +391,7 @@ 
Salmon survival example  # convert logit-transformed survival back to proportional
   dplyr::mutate(survival = plogis(logit.s)) -> model_data

 
Inspect the data

-
+
 dplyr::glimpse(model_data)

 
## Rows: 42
 ## Columns: 6
@@ -395,26 +404,26 @@ 
Salmon survival examplePlot features of the outcome variable, which shows that it is a
 proportional variable with particular restrictions that we want to
 model:

-
+
 plot_mvgam_series(data = model_data, y = 'survival')

-

+

 

 
A State-Space Beta regression
 

 
mvgam can easily handle data that are bounded at 0 and 1
 with a Beta observation model (using the mgcv function
-betar(), see ?mgcv::betar for details). First
+betar(), see ?mgcv::betar for details). First
 we will fit a simple State-Space model that uses a Random Walk dynamic
 process model with no predictors and a Beta observation model:

-
+
 mod0 <- mvgam(formula = survival ~ 1,
              trend_model = 'RW',
-             family = betar(),
+             family = betar(),
              data = model_data)

 
The summary of this model shows good behaviour of the Hamiltonian
 Monte Carlo sampler and provides useful summaries on the Beta
 observation model parameters:

-
+
 summary(mod0)

 
## GAM formula:
 ## survival ~ 1
@@ -439,33 +448,32 @@ 
A State-Space Beta regression## 
 ## Observation precision parameter estimates:
 ##        2.5% 50% 97.5% Rhat n.eff
-## phi[1]  150 310   550 1.01   386
+## phi[1]  160 310   570    1   617
 ## 
 ## GAM coefficient (beta) estimates:
 ##             2.5%  50% 97.5% Rhat n.eff
-## (Intercept) -4.4 -3.4  -2.4 1.05    81
+## (Intercept) -4.2 -3.4  -2.5 1.04    87
 ## 
 ## Latent trend variance estimates:
 ##          2.5%  50% 97.5% Rhat n.eff
-## sigma[1] 0.15 0.33  0.56 1.04   125
+## sigma[1] 0.18 0.33  0.54 1.01   199
 ## 
 ## Stan MCMC diagnostics:
 ## n_eff / iter looks reasonable for all parameters
-## Rhats above 1.05 found for 14 parameters
-## *Diagnose further to investigate why the chains have not mixed
+## 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

 
A plot of the underlying dynamic component shows how it has easily
 handled the temporal evolution of the time series:

-
+
 plot(mod0, type = 'trend')

-

+

 
Posterior hindcasts are also good and will automatically respect the
 bounding at 0 and 1:

-
+
 plot(mod0, type = 'forecast')

-

+

 

 

 
Including time-varying upwelling effects
@@ -477,26 +485,22 @@ 
Including time-varying upwelli
 observations, but this time will include a dynamic() effect
 of CUI.apr in the latent process model. We do not specify
 the \(\rho\) parameter, instead opting
-to estimate it using a Hilbert space approximate GP. Note also that we
-no longer need an intercept in the observation model as this may be hard
-to identify well alongside the flexible process model.
-mvgam allows us to fix the coefficient for the intercept at
-0 by specifying a -1 in the observation formula:

-
-mod1 <- mvgam(formula = survival ~ -1,
-              trend_formula = ~ dynamic(CUI.apr, k = 20),
+to estimate it using a Hilbert space approximate GP:

+
+mod1 <- mvgam(formula = survival ~ 1,
+              trend_formula = ~ dynamic(CUI.apr, k = 25),
               trend_model = 'RW',
-              family = betar(),
+              family = betar(),
               data = model_data)

 
The summary for this model now includes estimates for the
 time-varying GP parameters:

-
+
 summary(mod1, include_betas = FALSE)

 
## GAM observation formula:
 ## survival ~ 1
 ## 
 ## GAM process formula:
-## ~dynamic(CUI.apr, k = 20)
+## ~dynamic(CUI.apr, k = 25)
 ## 
 ## Family:
 ## beta
@@ -521,40 +525,41 @@ 
Including time-varying upwelli
 ## 
 ## Observation precision parameter estimates:
 ##        2.5% 50% 97.5% Rhat n.eff
-## phi[1]  190 360   680    1   504
+## phi[1]  180 360   650 1.01   832
 ## 
 ## GAM observation model coefficient (beta) estimates:
-##             2.5% 50% 97.5% Rhat n.eff
-## (Intercept)    0   0     0  NaN   NaN
+##             2.5%  50% 97.5% Rhat n.eff
+## (Intercept)   -4 -3.2  -2.4  1.1    43
 ## 
 ## Process error parameter estimates:
-##          2.5%  50% 97.5% Rhat n.eff
-## sigma[1] 0.18 0.31   0.5 1.01   334
+##          2.5% 50% 97.5% Rhat n.eff
+## sigma[1] 0.16 0.3   0.5 1.04   134
 ## 
 ## GAM process model coefficient (beta) estimates:
-##                   2.5%  50% 97.5% Rhat n.eff
-## (Intercept)_trend -4.3 -3.3  -2.4    1  1246
+##                            2.5%  50% 97.5% Rhat n.eff
+## gp(time):CUI.apr.1_trend -0.081 0.12  0.51    1   687
 ## 
 ## GAM process model gp term marginal deviation (alpha) and length scale (rho) estimates:
 ##                                2.5%  50% 97.5% Rhat n.eff
-## alpha_gp_time_byCUI_apr_trend 0.023 0.33  1.00 1.02   276
-## rho_gp_time_byCUI_apr_trend   0.031 0.13  0.72 1.01   317
+## alpha_gp_time_byCUI_apr_trend 0.015 0.31  1.00 1.01   421
+## rho_gp_time_byCUI_apr_trend   0.033 0.15  0.81 1.00   354
 ## 
 ## Stan MCMC diagnostics:
 ## n_eff / iter looks reasonable for all parameters
-## Rhat looks reasonable for all parameters
-## 81 of 2000 iterations ended with a divergence (4.05%)
+## Rhats above 1.05 found for 42 parameters
+## *Diagnose further to investigate why the chains have not mixed
+## 119 of 2000 iterations ended with a divergence (5.95%)
 ## *Try running with larger adapt_delta to remove the divergences
 ## 0 of 2000 iterations saturated the maximum tree depth of 12 (0%)
 ## E-FMI indicated no pathological behavior

 
The estimates for the underlying dynamic process haven’t changed
 much:

-
+
 plot(mod1, type = 'trend')

-

+

 
But the process error parameter \(\sigma\) is slightly smaller for this model
 than for the first model:

-
+
 # Extract estimates of the process error 'sigma' for each model
 mod0_sigma <- as.data.frame(mod0, variable = 'sigma', regex = TRUE) %>%
   dplyr::mutate(model = 'Mod0')
@@ -567,22 +572,27 @@ 
Including time-varying upwelli
 ggplot(sigmas, aes(y = `sigma[1]`, fill = model)) +
   geom_density(alpha = 0.3, colour = NA) +
   coord_flip()

-

+

 
Why does the process error not need to be as flexible in the second
 model? Because the estimates of this dynamic process are now informed
 partly by the time-varying effect of upwelling, which we can visualise
-using plot_predictions:

-
-plot_predictions(mod1, newdata = datagrid(CUI.apr = fivenum,
+on the link scale using plot() with
+trend_effects = TRUE:

+
+plot(mod1, type = 'smooth', trend_effects = TRUE)

+

+
Or on the outcome scale, at a range of possible CUI.apr
+values, using plot_predictions():

+
+plot_predictions(mod1, newdata = datagrid(CUI.apr = fivenum,
                                           time = unique),
                  by = c('time', 'CUI.apr', 'CUI.apr'),
-                 type = 'link',
                  process_error = FALSE)

 
## Warning: These arguments are not supported for models of class `mvgam`:
 ## process_error. Please file a request on Github if you believe that additional
 ## arguments should be supported:
 ## https://github.com/vincentarelbundock/marginaleffects/issues

-

+

 

 

 
Comparing model predictive performances
@@ -593,14 +603,14 @@ 
Comparing model predictive perf
 First, we can compare models based on in-sample approximate
 leave-one-out cross-validation as implemented in the popular
 loo package:

-
-loo_compare(mod0, mod1)

+
+loo_compare(mod0, mod1)

 
## Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.
 
 ## Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.

 
##      elpd_diff se_diff
 ## mod1  0.0       0.0   
-## mod0 -2.8       1.3

+## mod0 -1.7       1.7

 
The second model has the larger Expected Log Predictive Density
 (ELPD), meaning that it is slightly favoured over the simpler model that
 did not include the time-varying upwelling effect. However, the two
@@ -615,21 +625,21 @@ 
Comparing model predictive perf
 sampling to reweight posterior predictions, acting as a kind of particle
 filter so that we don’t need to refit the model too often (you can read
 more about how this process works in Bürkner et al. 2020).

-
+
 lfo_mod0 <- lfo_cv(mod0, min_t = 30)
 lfo_mod1 <- lfo_cv(mod1, min_t = 30)

 
The model with the time-varying upwelling effect tends to provides
 better 1-step ahead forecasts, with a higher total forecast ELPD

-
-sum(lfo_mod0$elpds)

-
## [1] 34.99668

 
+sum(lfo_mod0$elpds)

+
## [1] 34.80785

+
 sum(lfo_mod1$elpds)

-
## [1] 36.44244

+
## [1] 36.43776

 
We can also plot the ELPDs for each model as a contrast. Here, values
 less than zero suggest the time-varying predictor model (Mod1) gives
 better 1-step ahead forecasts:

-
+
 plot(x = 1:length(lfo_mod0$elpds) + 30,
      y = lfo_mod0$elpds - lfo_mod1$elpds,
      ylab = 'ELPDmod0 - ELPDmod1',
@@ -638,7 +648,7 @@ 
Comparing model predictive perf
      col = 'darkred',
      bty = 'l')
 abline(h = 0, lty = 'dashed')

-

+

 
A useful exercise to further expand this model would be to think
 about what kinds of predictors might impact measurement error, which
 could easily be implemented into the observation formula in
diff --git a/docs/articles/time_varying_effects_files/figure-html/unnamed-chunk-10-1.png b/docs/articles/time_varying_effects_files/figure-html/unnamed-chunk-10-1.png
index 482b8442..bed85c1c 100644
Binary files a/docs/articles/time_varying_effects_files/figure-html/unnamed-chunk-10-1.png and b/docs/articles/time_varying_effects_files/figure-html/unnamed-chunk-10-1.png differ
diff --git a/docs/articles/time_varying_effects_files/figure-html/unnamed-chunk-11-1.png b/docs/articles/time_varying_effects_files/figure-html/unnamed-chunk-11-1.png
index 95ca3667..7a8e2dbd 100644
Binary files a/docs/articles/time_varying_effects_files/figure-html/unnamed-chunk-11-1.png and b/docs/articles/time_varying_effects_files/figure-html/unnamed-chunk-11-1.png differ
diff --git a/docs/articles/time_varying_effects_files/figure-html/unnamed-chunk-12-1.png b/docs/articles/time_varying_effects_files/figure-html/unnamed-chunk-12-1.png
index 8547c620..e27b1a0e 100644
Binary files a/docs/articles/time_varying_effects_files/figure-html/unnamed-chunk-12-1.png and b/docs/articles/time_varying_effects_files/figure-html/unnamed-chunk-12-1.png differ
diff --git a/docs/articles/time_varying_effects_files/figure-html/unnamed-chunk-13-1.png b/docs/articles/time_varying_effects_files/figure-html/unnamed-chunk-13-1.png
index 95d32755..6c8cc0c2 100644
Binary files a/docs/articles/time_varying_effects_files/figure-html/unnamed-chunk-13-1.png and b/docs/articles/time_varying_effects_files/figure-html/unnamed-chunk-13-1.png differ
diff --git a/docs/articles/time_varying_effects_files/figure-html/unnamed-chunk-17-1.png b/docs/articles/time_varying_effects_files/figure-html/unnamed-chunk-17-1.png
index a3f45dfb..91a5ded6 100644
Binary files a/docs/articles/time_varying_effects_files/figure-html/unnamed-chunk-17-1.png and b/docs/articles/time_varying_effects_files/figure-html/unnamed-chunk-17-1.png differ
diff --git a/docs/articles/time_varying_effects_files/figure-html/unnamed-chunk-18-1.png b/docs/articles/time_varying_effects_files/figure-html/unnamed-chunk-18-1.png
index 45dcd1e3..9527d4e1 100644
Binary files a/docs/articles/time_varying_effects_files/figure-html/unnamed-chunk-18-1.png and b/docs/articles/time_varying_effects_files/figure-html/unnamed-chunk-18-1.png differ
diff --git a/docs/articles/time_varying_effects_files/figure-html/unnamed-chunk-19-1.png b/docs/articles/time_varying_effects_files/figure-html/unnamed-chunk-19-1.png
new file mode 100644
index 00000000..b8c3f6ee
Binary files /dev/null and b/docs/articles/time_varying_effects_files/figure-html/unnamed-chunk-19-1.png differ
diff --git a/docs/articles/time_varying_effects_files/figure-html/unnamed-chunk-23-1.png b/docs/articles/time_varying_effects_files/figure-html/unnamed-chunk-23-1.png
new file mode 100644
index 00000000..3d281c36
Binary files /dev/null and b/docs/articles/time_varying_effects_files/figure-html/unnamed-chunk-23-1.png differ
diff --git a/docs/articles/time_varying_effects_files/figure-html/unnamed-chunk-27-1.png b/docs/articles/time_varying_effects_files/figure-html/unnamed-chunk-27-1.png
index e330f19d..e3a28fc8 100644
Binary files a/docs/articles/time_varying_effects_files/figure-html/unnamed-chunk-27-1.png and b/docs/articles/time_varying_effects_files/figure-html/unnamed-chunk-27-1.png differ
diff --git a/docs/articles/time_varying_effects_files/figure-html/unnamed-chunk-28-1.png b/docs/articles/time_varying_effects_files/figure-html/unnamed-chunk-28-1.png
new file mode 100644
index 00000000..00c0ae6c
Binary files /dev/null and b/docs/articles/time_varying_effects_files/figure-html/unnamed-chunk-28-1.png differ
diff --git a/docs/articles/time_varying_effects_files/figure-html/unnamed-chunk-32-1.png b/docs/articles/time_varying_effects_files/figure-html/unnamed-chunk-32-1.png
index c131dfe7..46282632 100644
Binary files a/docs/articles/time_varying_effects_files/figure-html/unnamed-chunk-32-1.png and b/docs/articles/time_varying_effects_files/figure-html/unnamed-chunk-32-1.png differ
diff --git a/docs/articles/time_varying_effects_files/figure-html/unnamed-chunk-33-1.png b/docs/articles/time_varying_effects_files/figure-html/unnamed-chunk-33-1.png
index 75846af3..97ee6f8c 100644
Binary files a/docs/articles/time_varying_effects_files/figure-html/unnamed-chunk-33-1.png and b/docs/articles/time_varying_effects_files/figure-html/unnamed-chunk-33-1.png differ
diff --git a/docs/articles/time_varying_effects_files/figure-html/unnamed-chunk-34-1.png b/docs/articles/time_varying_effects_files/figure-html/unnamed-chunk-34-1.png
new file mode 100644
index 00000000..8871e518
Binary files /dev/null and b/docs/articles/time_varying_effects_files/figure-html/unnamed-chunk-34-1.png differ
diff --git a/docs/articles/time_varying_effects_files/figure-html/unnamed-chunk-35-1.png b/docs/articles/time_varying_effects_files/figure-html/unnamed-chunk-35-1.png
new file mode 100644
index 00000000..1a0ad537
Binary files /dev/null and b/docs/articles/time_varying_effects_files/figure-html/unnamed-chunk-35-1.png differ
diff --git a/docs/articles/time_varying_effects_files/figure-html/unnamed-chunk-40-1.png b/docs/articles/time_varying_effects_files/figure-html/unnamed-chunk-40-1.png
new file mode 100644
index 00000000..bbf06662
Binary files /dev/null and b/docs/articles/time_varying_effects_files/figure-html/unnamed-chunk-40-1.png differ
diff --git a/src/RcppExports.o b/src/RcppExports.o
index c1d41286..4acd929b 100644
Binary files a/src/RcppExports.o and b/src/RcppExports.o differ
diff --git a/src/mvgam.dll b/src/mvgam.dll
index a5b9a478..882ad816 100644
Binary files a/src/mvgam.dll and b/src/mvgam.dll differ
diff --git a/src/trend_funs.o b/src/trend_funs.o
index ad153f68..30ae531e 100644
Binary files a/src/trend_funs.o and b/src/trend_funs.o differ
diff --git a/tests/testthat/Rplots.pdf b/tests/testthat/Rplots.pdf
index cf61d416..a984bf68 100644
Binary files a/tests/testthat/Rplots.pdf and b/tests/testthat/Rplots.pdf differ
diff --git a/vignettes/time_varying_effects.Rmd b/vignettes/time_varying_effects.Rmd
index f2741ab6..6d428d2d 100644
--- a/vignettes/time_varying_effects.Rmd
+++ b/vignettes/time_varying_effects.Rmd
@@ -80,13 +80,13 @@ plot_mvgam_series(data = data_train, newdata = data_test, y = 'out')
 ### 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 = 20)`
 ```{r, include=FALSE}
-mod <- mvgam(out ~ dynamic(temp, rho = 8, stationary = TRUE, k = 20),
+mod <- mvgam(out ~ dynamic(temp, rho = 8, stationary = TRUE, k = 40),
              family = gaussian(),
              data = data_train)
 ```
 
 ```{r, eval=FALSE}
-mod <- mvgam(out ~ dynamic(temp, rho = 8, stationary = TRUE, k = 20),
+mod <- mvgam(out ~ dynamic(temp, rho = 8, stationary = TRUE, k = 40),
              family = gaussian(),
              data = data_train)
 ```
@@ -123,15 +123,15 @@ This results in sensible forecasts of the observations as well
 plot(mod, type = 'forecast', newdata = data_test)
 ```
 
-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 = 20`.
+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 = 20)`.
 ```{r include=FALSE}
-mod <- mvgam(out ~ dynamic(temp, k = 20),
+mod <- mvgam(out ~ dynamic(temp, k = 40),
              family = gaussian(),
              data = data_train)
 ```
 
 ```{r eval=FALSE}
-mod <- mvgam(out ~ dynamic(temp, k = 20),
+mod <- mvgam(out ~ dynamic(temp, k = 40),
              family = gaussian(),
              data = data_train)
 ```
@@ -141,7 +141,15 @@ This model summary now contains estimates for the marginal deviation and length
 summary(mod, include_betas = FALSE)
 ```
 
-The `plot_predictions` call shows that the effect in this case is similar to what we estimated above:
+Effects for `gp()` terms can also be plotted as smooths:
+```{r}
+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)
+```
+
+Both the above plot and the below `plot_predictions()` call show that the effect in this case is similar to what we estimated in the approximate GP smooth model above:
 ```{r}
 plot_predictions(mod, 
                  newdata = datagrid(time = unique,
@@ -221,18 +229,18 @@ plot(mod0, type = 'forecast')
 
 
 ### Including time-varying upwelling effects
-Now we can increase the complexity of our model by constructing and fitting a State-Space model with a time-varying effect of the coastal upwelling index in addition to the autoregressive dynamics. We again use a Beta observation model to capture the restrictions of our proportional observations, but this time will include a `dynamic()` effect of `CUI.apr` in the latent process model. We do not specify the $\rho$ parameter, instead opting to estimate it using a Hilbert space approximate GP. Note also that we no longer need an intercept in the observation model as this may be hard to identify well alongside the flexible process model. `mvgam` allows us to fix the coefficient for the intercept at 0 by specifying a `-1` in the observation formula:
+Now we can increase the complexity of our model by constructing and fitting a State-Space model with a time-varying effect of the coastal upwelling index in addition to the autoregressive dynamics. We again use a Beta observation model to capture the restrictions of our proportional observations, but this time will include a `dynamic()` effect of `CUI.apr` in the latent process model. We do not specify the $\rho$ parameter, instead opting to estimate it using a Hilbert space approximate GP:
 ```{r include=FALSE}
-mod1 <- mvgam(formula = survival ~ -1,
-              trend_formula = ~ dynamic(CUI.apr, k = 20),
+mod1 <- mvgam(formula = survival ~ 1,
+              trend_formula = ~ dynamic(CUI.apr, k = 25),
               trend_model = 'RW',
               family = betar(),
               data = model_data)
 ```
 
 ```{r eval=FALSE}
-mod1 <- mvgam(formula = survival ~ -1,
-              trend_formula = ~ dynamic(CUI.apr, k = 20),
+mod1 <- mvgam(formula = survival ~ 1,
+              trend_formula = ~ dynamic(CUI.apr, k = 25),
               trend_model = 'RW',
               family = betar(),
               data = model_data)
@@ -264,12 +272,16 @@ ggplot(sigmas, aes(y = `sigma[1]`, fill = model)) +
   coord_flip()
 ```
 
-Why does the process error not need to be as flexible in the second model? Because the estimates of this dynamic process are now informed partly by the time-varying effect of upwelling, which we can visualise using `plot_predictions`:
+Why does the process error not need to be as flexible in the second model? Because the estimates of this dynamic process are now informed partly by the time-varying effect of upwelling, which we can visualise on the link scale using `plot()` with `trend_effects = TRUE`:
+```{r}
+plot(mod1, type = 'smooth', trend_effects = TRUE)
+```
+
+Or on the outcome scale, at a range of possible `CUI.apr` values, using `plot_predictions()`:
 ```{r}
 plot_predictions(mod1, newdata = datagrid(CUI.apr = fivenum,
                                           time = unique),
                  by = c('time', 'CUI.apr', 'CUI.apr'),
-                 type = 'link',
                  process_error = FALSE)
 ```