diff --git a/R/add_nmixture.R b/R/add_nmixture.R index 8a5ad4ef..969d8638 100644 --- a/R/add_nmixture.R +++ b/R/add_nmixture.R @@ -15,13 +15,36 @@ add_nmixture = function(model_file, stop('Max abundances must be supplied as a variable named "cap" for N-mixture models', call. = FALSE) } - cap <- data_train$cap + + if(inherits(data_train, 'data.frame')){ + cap = data_train %>% + dplyr::arrange(series, time) %>% + dplyr::pull(cap) + } else { + cap = data.frame(series = data_train$series, + cap = data_train$cap, + time = data_train$time)%>% + dplyr::arrange(series, time) %>% + dplyr::pull(cap) + } + if(!is.null(data_test)){ if(!(exists('cap', where = data_test))) { stop('Max abundances must be supplied in test data as a variable named "cap" for N-mixture models', call. = FALSE) } - cap <- c(cap, data_test$cap) + if(inherits(data_test, 'data.frame')){ + captest = data_test %>% + dplyr::arrange(series, time) %>% + dplyr::pull(cap) + } else { + captest = data.frame(series = data_test$series, + cap = data_test$cap, + time = data_test$time)%>% + dplyr::arrange(series, time) %>% + dplyr::pull(cap) + } + cap <- c(cap, captest) } validate_pos_integers(cap) @@ -33,7 +56,7 @@ add_nmixture = function(model_file, model_data$cap <- as.vector(cap) - if(any(model_data$cap < model_data$y)){ + if(any(model_data$cap[model_data$obs_ind] < model_data$flat_ys)){ stop(paste0('Some "cap" terms are < the observed counts. This is not allowed'), call. = FALSE) } @@ -243,7 +266,7 @@ add_nmixture = function(model_file, 'array[n, n_series] int latent_ypred;\n', 'array[total_obs] int latent_truncpred;\n', 'vector[n_nonmissing] flat_ps;\n', - 'int flat_caps[n_nonmissing];', + 'int flat_caps[n_nonmissing];\n', 'vector[total_obs] flat_trends;\n', 'vector[n_nonmissing] flat_trends_nonmis;\n', 'vector[total_obs] detprob;\n', diff --git a/R/get_mvgam_priors.R b/R/get_mvgam_priors.R index 60ee0e2e..228a5cd7 100644 --- a/R/get_mvgam_priors.R +++ b/R/get_mvgam_priors.R @@ -462,7 +462,7 @@ get_mvgam_priors = function(formula, } # Remove sigma prior if this is an N-mixture with no dynamics - if(add_nmix){ + if(add_nmix & trend_model == 'None'){ out <- out[-grep('vector[n_lv] sigma;', out$param_name, fixed = TRUE),] diff --git a/R/gp.R b/R/gp.R index d257d049..04ec030f 100644 --- a/R/gp.R +++ b/R/gp.R @@ -158,10 +158,15 @@ make_gp_additions = function(gp_details, data, # Add coefficient indices to attribute table and to Stan data for(covariate in seq_along(gp_att_table)){ - # coef_indices <- grep(gp_att_table[[covariate]]$name, - # names(coef(mgcv_model)), fixed = TRUE) - coef_indices <- which(grepl(gp_att_table[[covariate]]$name, - names(coef(mgcv_model)), fixed = TRUE) & + # coef_indices <- which(grepl(gp_att_table[[covariate]]$name, + # names(coef(mgcv_model)), fixed = TRUE) & + # !grepl(paste0(gp_att_table[[covariate]]$name,':'), + # names(coef(mgcv_model)), fixed = TRUE) == TRUE) + + coef_indices <- which(grepl(paste0(gsub("([()])","\\\\\\1", + gp_att_table[[covariate]]$name), + '\\.+[0-9]'), + names(coef(mgcv_model)), fixed = FALSE) & !grepl(paste0(gp_att_table[[covariate]]$name,':'), names(coef(mgcv_model)), fixed = TRUE) == TRUE) diff --git a/R/mvgam.R b/R/mvgam.R index baaafd2a..32ff9743 100644 --- a/R/mvgam.R +++ b/R/mvgam.R @@ -1718,8 +1718,7 @@ mvgam = function(formula, # Auto-format the model file if(autoformat){ if(requireNamespace('cmdstanr') & cmdstanr::cmdstan_version() >= "2.29.0") { - tmp_file <- cmdstanr::write_stan_file(vectorised$model_file) - vectorised$model_file <- .autoformat(tmp_file, + vectorised$model_file <- .autoformat(vectorised$model_file, overwrite_file = FALSE) } vectorised$model_file <- readLines(textConnection(vectorised$model_file), @@ -1908,7 +1907,7 @@ mvgam = function(formula, cpp_options = list(stan_threads = TRUE)) } else { cmd_mod <- cmdstanr::cmdstan_model(cmdstanr::write_stan_file(vectorised$model_file), - stanc_options = list('O1')) + stanc_options = list('O1'),) } } else { diff --git a/R/stan_utils.R b/R/stan_utils.R index f6d6efac..d295cbe8 100644 --- a/R/stan_utils.R +++ b/R/stan_utils.R @@ -36,12 +36,313 @@ remove_likelihood = function(model_file){ } #' @noRd -.autoformat <- function(stan_file, overwrite_file = TRUE) { +.autoformat <- function(stan_file, overwrite_file = TRUE){ + + # Old ways of specifying arrays have been converted to errors in + # the latest version of Cmdstan (2.34.0); this coincides with + # a decision to stop automatically replacing these depracations with + # the canonicalizer, so we have no choice but to replace the old + # syntax with this ugly bit of code: + if(requireNamespace('cmdstanr') & cmdstanr::cmdstan_version() >= "2.33.0"){ + + # Data modifications + stan_file[grep("int ytimes[n, n_series]; // time-ordered matrix (which col in X belongs to each [time, series] observation?)", + stan_file, fixed = TRUE)] <- + 'array[n, n_series] int ytimes; // time-ordered matrix (which col in X belongs to each [time, series] observation?)' + + stan_file[grep("int flat_ys[n_nonmissing]; // flattened nonmissing observations", + stan_file, fixed = TRUE)] <- + 'array[n_nonmissing] int flat_ys; // flattened nonmissing observations' + + stan_file[grep("int obs_ind[n_nonmissing]; // indices of nonmissing observations", + stan_file, fixed = TRUE)] <- + "array[n_nonmissing] int obs_ind; // indices of nonmissing observations" + + if(any(grepl('int ytimes_trend[n, n_lv]; // time-ordered matrix for latent states', + stan_file, fixed = TRUE))){ + stan_file[grep("int ytimes_trend[n, n_lv]; // time-ordered matrix for latent states", + stan_file, fixed = TRUE)] <- + "array[n, n_lv] int ytimes_trend;" + } + + if(any(grepl('int idx', stan_file) & + grepl('// discontiguous index values', + stan_file, fixed = TRUE))){ + lines_replace <- which(grepl('int idx', stan_file) & + grepl('// discontiguous index values', + stan_file, fixed = TRUE)) + for(i in lines_replace){ + split_line <- strsplit(stan_file[i], ' ')[[1]] + + idxnum <- gsub(';', '', + gsub("\\s*\\[[^\\]+\\]", "", + as.character(split_line[2]))) + idx_length <- gsub("\\]", "", gsub("\\[", "", + regmatches(split_line[2], + gregexpr("\\[.*?\\]", split_line[2]))[[1]])) + + stan_file[i] <- + paste0('array[', + idx_length, + '] int ', + idxnum, + '; // discontiguous index values') + } + } + + if(any(grepl('int cap[total_obs]; // upper limits of latent abundances', + stan_file, fixed = TRUE))){ + stan_file[grep('int cap[total_obs]; // upper limits of latent abundances', + stan_file, fixed = TRUE)] <- + 'array[total_obs] int cap; // upper limits of latent abundances' + + stan_file[grep('int flat_caps[n_nonmissing];', + stan_file, fixed = TRUE)] <- + 'array[n_nonmissing] int flat_caps;' + } + + # Model modifications + if(any(grepl('real flat_phis[n_nonmissing];', + stan_file, fixed = TRUE))){ + stan_file[grep("real flat_phis[n_nonmissing];", + stan_file, fixed = TRUE)] <- + "array[n_nonmissing] real flat_phis;" + } + + # n-mixture modifications + if(any(grepl('real p_ub = poisson_cdf(max_k, lambda);', + stan_file, fixed = TRUE))){ + stan_file[grep('real p_ub = poisson_cdf(max_k, lambda);', + stan_file, fixed = TRUE)] <- + 'real p_ub = poisson_cdf(max_k | lambda);' + } + + # trend_formula modifications + if(any(grepl('int trend_idx', stan_file) & + grepl('// discontiguous index values', + stan_file, fixed = TRUE))){ + lines_replace <- which(grepl('int trend_idx', stan_file) & + grepl('// discontiguous index values', + stan_file, fixed = TRUE)) + for(i in lines_replace){ + split_line <- strsplit(stan_file[i], ' ')[[1]] + + trend_idxnum <- gsub(';', '', + gsub("\\s*\\[[^\\]+\\]", "", + as.character(split_line[2]))) + idx_length <- gsub("\\]", "", gsub("\\[", "", + regmatches(split_line[2], + gregexpr("\\[.*?\\]", split_line[2]))[[1]])) + + stan_file[i] <- + paste0('array[', + idx_length, + '] int ', + trend_idxnum, + '; // discontiguous index values') + } + } + + if(any(grepl('vector[n_series] trend_raw[n];', + stan_file, fixed = TRUE))){ + stan_file[grep("vector[n_series] trend_raw[n];", + stan_file, fixed = TRUE)] <- + "array[n] vector[n_series] trend_raw;" + } + + if(any(grepl('vector[n_lv] error[n];', + stan_file, fixed = TRUE))){ + stan_file[grep("vector[n_lv] error[n];", + stan_file, fixed = TRUE)] <- + "array[n] vector[n_lv] error;" + } + + if(any(grepl('vector[n_series] error[n];', + stan_file, fixed = TRUE))){ + stan_file[grep("vector[n_series] error[n];", + stan_file, fixed = TRUE)] <- + "array[n] vector[n_series] error;" + } + + if(any(grepl('vector[n_lv] LV[n];', + stan_file, fixed = TRUE))){ + stan_file[grep("vector[n_lv] LV[n];", + stan_file, fixed = TRUE)] <- + "array[n] vector[n_lv] LV;" + } + + if(any(grepl('vector[n_series] mu[n - 1];', + stan_file, fixed = TRUE))){ + stan_file[grep("vector[n_series] mu[n - 1];", + stan_file, fixed = TRUE)] <- + "array[n - 1] vector[n_series] mu;" + } + + if(any(grepl('vector[n_lv] mu[n - 1];', + stan_file, fixed = TRUE))){ + stan_file[grep("vector[n_lv] mu[n - 1];", + stan_file, fixed = TRUE)] <- + "array[n - 1] vector[n_lv] mu;" + } + + if(any(grepl('vector[n_series] mu[n];', + stan_file, fixed = TRUE))){ + stan_file[grep("vector[n_series] mu[n];", + stan_file, fixed = TRUE)] <- + "array[n] vector[n_series] mu;" + } + + if(any(grepl('vector[n_lv] mu[n];', + stan_file, fixed = TRUE))){ + stan_file[grep("vector[n_lv] mu[n];", + stan_file, fixed = TRUE)] <- + "array[n] vector[n_lv] mu;" + } + # Generated quantity modifications + if(any(grepl('real ypred[n, n_series];', + stan_file, fixed = TRUE))){ + stan_file[grep("real ypred[n, n_series];", + stan_file, fixed = TRUE)] <- + "array[n, n_series] real ypred;" + } + + if(any(grepl('real ypred[n, n_series];', + stan_file, fixed = TRUE))){ + stan_file[grep("real ypred[n, n_series];", + stan_file, fixed = TRUE)] <- + "array[n, n_series] real ypred;" + } + + # ARMA model modifications + if(any(grepl('vector[n_series] epsilon[n];', + stan_file, fixed = TRUE))){ + stan_file[grep("vector[n_series] epsilon[n];", + stan_file, fixed = TRUE)] <- + "array[n] vector[n_series] epsilon;" + } + + if(any(grepl('vector[n_lv] epsilon[n];', + stan_file, fixed = TRUE))){ + stan_file[grep("vector[n_lv] epsilon[n];", + stan_file, fixed = TRUE)] <- + "array[n] vector[n_lv] epsilon;" + } + + # VARMA model modifications + if(any(grepl('matrix[n_series, n_series] P[1];', + stan_file, fixed = TRUE))){ + stan_file[grep("matrix[n_series, n_series] P[1];", + stan_file, fixed = TRUE)] <- + "array[1] matrix[n_series, n_series] P;" + + stan_file[grep("matrix[n_series, n_series] phiGamma[2, 1];", + stan_file, fixed = TRUE)] <- + "array[2, 1] matrix[n_series, n_series] phiGamma;" + } + + if(any(grepl('matrix initial_joint_var(matrix Sigma, matrix[] phi, matrix[] theta) {', + stan_file, fixed = TRUE))){ + stan_file[grep("matrix initial_joint_var(matrix Sigma, matrix[] phi, matrix[] theta) {", + stan_file, fixed = TRUE)] <- + "matrix initial_joint_var(matrix Sigma, array[] matrix phi, array[] matrix theta) {" + } + + if(any(grepl('matrix[n_lv, n_lv] P[1];', + stan_file, fixed = TRUE))){ + stan_file[grep("matrix[n_lv, n_lv] P[1];", + stan_file, fixed = TRUE)] <- + "array[1] matrix[n_lv, n_lv] P;" + + stan_file[grep("matrix[n_lv, n_lv] R[1];", + stan_file, fixed = TRUE)] <- + "array[1] matrix[n_lv, n_lv] R;" + + stan_file[grep("matrix[n_lv, n_lv] A_init[1];", + stan_file, fixed = TRUE)] <- + "array[1] matrix[n_lv, n_lv] A_init;" + + stan_file[grep("matrix[n_lv, n_lv] theta_init[1];", + stan_file, fixed = TRUE)] <- + "array[1] matrix[n_lv, n_lv] theta_init;" + } + + if(any(grepl('matrix[n_series, n_series] R[1];', + stan_file, fixed = TRUE))){ + + stan_file[grep("matrix[n_series, n_series] R[1];", + stan_file, fixed = TRUE)] <- + "array[1] matrix[n_series, n_series] R;" + + stan_file[grep("matrix[n_series, n_series] A_init[1];", + stan_file, fixed = TRUE)] <- + "array[1] matrix[n_series, n_series] A_init;" + + stan_file[grep("matrix[n_series, n_series] theta_init[1];", + stan_file, fixed = TRUE)] <- + "array[1] matrix[n_series, n_series] theta_init;" + } + + if(any(grepl('matrix[] rev_mapping(matrix[] P, matrix Sigma) {', + stan_file, fixed = TRUE))){ + stan_file[grep("matrix[] rev_mapping(matrix[] P, matrix Sigma) {", + stan_file, fixed = TRUE)] <- + "array[] matrix rev_mapping(array[] matrix P, matrix Sigma) {" + + stan_file[grep("matrix[m, m] phi_for[p, p]; matrix[m, m] phi_rev[p, p];", + stan_file, fixed = TRUE)] <- + 'array[p, p] matrix[m, m] phi_for; array[p, p] matrix[m, m] phi_rev;' + + stan_file[grep("matrix[m, m] Sigma_for[p+1]; matrix[m, m] Sigma_rev[p+1];", + stan_file, fixed = TRUE)] <- + 'array[p+1] matrix[m, m] Sigma_for; array[p+1] matrix[m, m] Sigma_rev;' + + stan_file[grep("matrix[m, m] S_for_list[p+1];", + stan_file, fixed = TRUE)] <- + 'array[p+1] matrix[m, m] S_for_list;' + } + + # VAR model modifications + if(any(grepl('matrix[n_lv, n_lv] phiGamma[2, 1];', + stan_file, fixed = TRUE))){ + stan_file[grep('matrix[n_lv, n_lv] phiGamma[2, 1];', + stan_file, fixed = TRUE)] <- + 'array[2, 1] matrix[n_lv, n_lv] phiGamma;' + } + + if(any(grepl('matrix[,] rev_mapping(matrix[] P, matrix Sigma) {', + stan_file, fixed = TRUE))){ + stan_file[grep("matrix[,] rev_mapping(matrix[] P, matrix Sigma) {", + stan_file, fixed = TRUE)] <- + "array[,] matrix rev_mapping(array[] matrix P, matrix Sigma) {" + + stan_file[grep("matrix[m, m] phi_for[p, p]; matrix[m, m] phi_rev[p, p];", + stan_file, fixed = TRUE)] <- + 'array[p, p] matrix[m, m] phi_for; array[p, p] matrix[m, m] phi_rev;' + + stan_file[grep("matrix[m, m] Sigma_for[p+1]; matrix[m, m] Sigma_rev[p+1];", + stan_file, fixed = TRUE)] <- + 'array[p+1] matrix[m, m] Sigma_for; array[p+1] matrix[m, m] Sigma_rev;' + + stan_file[grep("matrix[m, m] S_for_list[p+1];", + stan_file, fixed = TRUE)] <- + 'array[p+1] matrix[m, m] S_for_list;' + + stan_file[grep("matrix[m, m] Gamma_trans[p+1];", + stan_file, fixed = TRUE)] <- + 'array[p+1] matrix[m, m] Gamma_trans;' + + stan_file[grep("matrix[m, m] phiGamma[2, p];", + stan_file, fixed = TRUE)] <- + 'array[2, p] matrix[m, m] phiGamma;' + } + } + + stan_file <- cmdstanr::write_stan_file(stan_file) cmdstan_mod <- cmdstanr::cmdstan_model(stan_file, compile = FALSE) out <- utils::capture.output( cmdstan_mod$format( max_line_length = 80, - canonicalize = list("deprecations", "parentheses"), + canonicalize = TRUE, overwrite_file = overwrite_file, backup = FALSE)) paste0(out, collapse = "\n") } diff --git a/R/summary.mvgam.R b/R/summary.mvgam.R index 1d93e2d4..cdbb66ca 100644 --- a/R/summary.mvgam.R +++ b/R/summary.mvgam.R @@ -647,9 +647,9 @@ if(!is.null(object$trend_call)){ if(!is.null(attr(object$trend_mgcv_model, 'gp_att_table'))){ gp_names <- clean_gpnames(unlist(purrr::map(attr(object$trend_mgcv_model, 'gp_att_table'), 'name'))) - alpha_params <- gsub('gp_', 'gp_trend_', gsub(':', 'by', gsub(')', '_', + alpha_params <- gsub('series', 'trend', gsub('gp_', 'gp_trend_', gsub(':', 'by', gsub(')', '_', gsub('(', '_', paste0('alpha_', gp_names), - fixed = TRUE), fixed = TRUE))) + fixed = TRUE), fixed = TRUE)))) alpha_summary <- mcmc_summary(object$model_output, alpha_params, ISB = TRUE, digits = digits, variational = object$algorithm %in% c('fullrank', 'meanfield', 'laplace', 'pathfinder'))[,c(3:7)] @@ -657,9 +657,9 @@ if(!is.null(object$trend_call)){ rownames(alpha_summary) <- paste0(gsub('series', 'trend', paste0('alpha_', gp_names)), '_trend') - rho_params <- gsub('gp_', 'gp_trend_', gsub(':', 'by', gsub(')', '_', + rho_params <- gsub('series', 'trend', gsub('gp_', 'gp_trend_', gsub(':', 'by', gsub(')', '_', gsub('(', '_', paste0('rho_', gp_names), - fixed = TRUE), fixed = TRUE))) + fixed = TRUE), fixed = TRUE)))) rho_summary <- mcmc_summary(object$model_output, rho_params, ISB = TRUE, digits = digits, variational = object$algorithm %in% c('fullrank', 'meanfield', 'laplace', 'pathfinder'))[,c(3:7)] diff --git a/R/validations.R b/R/validations.R index 06cb04ba..8b2d868a 100644 --- a/R/validations.R +++ b/R/validations.R @@ -6,6 +6,8 @@ validate_series_time = function(data, name = 'data'){ # Series label must be present as a factor and # must contain a time variable if(inherits(data, 'data.frame')){ + data %>% + dplyr::ungroup() -> data if(!'series' %in% colnames(data)){ data$series <- factor('series1') } diff --git a/README.md b/README.md index b035b155..dd47756c 100644 --- a/README.md +++ b/README.md @@ -317,29 +317,29 @@ summary(lynx_mvgam) #> #> #> GAM coefficient (beta) estimates: -#> 2.5% 50% 97.5% Rhat n_eff -#> (Intercept) 6.200 6.6000 7.000 1 1020 -#> s(season).1 -0.590 0.0500 0.700 1 1094 -#> s(season).2 -0.240 0.8100 1.800 1 417 -#> s(season).3 -0.057 1.2000 2.500 1 365 -#> s(season).4 -0.510 0.4200 1.400 1 951 -#> s(season).5 -1.200 -0.1400 1.000 1 538 -#> s(season).6 -1.100 0.0074 1.100 1 632 -#> s(season).7 -0.790 0.3800 1.400 1 847 -#> s(season).8 -1.000 0.2900 1.800 1 413 -#> s(season).9 -1.100 -0.2700 0.670 1 574 -#> s(season).10 -1.400 -0.6900 -0.025 1 641 +#> 2.5% 50% 97.5% Rhat n_eff +#> (Intercept) 6.10 6.600 7.000 1.00 625 +#> s(season).1 -0.63 0.040 0.760 1.00 788 +#> s(season).2 -0.25 0.800 1.800 1.01 520 +#> s(season).3 -0.10 1.200 2.500 1.01 459 +#> s(season).4 -0.53 0.420 1.400 1.01 907 +#> s(season).5 -1.20 -0.130 0.960 1.01 622 +#> s(season).6 -1.10 0.022 1.000 1.01 616 +#> s(season).7 -0.77 0.320 1.500 1.00 900 +#> s(season).8 -1.00 0.200 1.900 1.01 453 +#> s(season).9 -1.20 -0.320 0.680 1.00 470 +#> s(season).10 -1.40 -0.670 -0.022 1.00 850 #> #> Approximate significance of GAM observation smooths: -#> edf Chi.sq p-value -#> s(season) 5 17851 0.28 +#> edf Chi.sq p-value +#> s(season) 3.66 18756 0.28 #> #> Latent trend AR parameter estimates: #> 2.5% 50% 97.5% Rhat n_eff -#> ar1[1] 0.74 1.10 1.400 1 695 -#> ar2[1] -0.82 -0.41 0.045 1 1630 -#> ar3[1] -0.47 -0.12 0.280 1 453 -#> sigma[1] 0.40 0.50 0.630 1 1101 +#> ar1[1] 0.74 1.10 1.400 1.00 829 +#> ar2[1] -0.84 -0.41 0.033 1.00 1598 +#> ar3[1] -0.47 -0.13 0.310 1.01 517 +#> sigma[1] 0.40 0.50 0.630 1.00 1278 #> #> Stan MCMC diagnostics: #> n_eff / iter looks reasonable for all parameters @@ -348,7 +348,7 @@ summary(lynx_mvgam) #> 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 Tue Jan 16 8:42:08 PM 2024. +#> Samples were drawn using NUTS(diag_e) at Fri Jan 19 12:32:58 PM 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) @@ -470,7 +470,7 @@ plot(lynx_mvgam, type = 'forecast', newdata = lynx_test) #> Out of sample CRPS: - #> [1] 2856.97 + #> [1] 2831.125 And the estimated latent trend component, again using the more flexible `plot_mvgam_...()` option to show first derivatives of the estimated @@ -592,31 +592,31 @@ summary(mod, include_betas = FALSE) #> #> Observation precision parameter estimates: #> 2.5% 50% 97.5% Rhat n_eff -#> phi[1] 5.4 8.3 12 1 1692 -#> phi[2] 5.6 8.7 13 1 1512 -#> phi[3] 5.6 8.4 12 1 1786 +#> phi[1] 5.5 8.4 12 1 1516 +#> phi[2] 5.8 8.7 13 1 1333 +#> phi[3] 5.6 8.5 12 1 1690 #> #> GAM coefficient (beta) estimates: -#> 2.5% 50% 97.5% Rhat n_eff -#> (Intercept) -0.15 0.2 0.46 1.01 891 +#> 2.5% 50% 97.5% Rhat n_eff +#> (Intercept) -0.19 0.19 0.44 1 822 #> #> Approximate significance of GAM observation smooths: #> edf Chi.sq p-value -#> s(season) 5.00 15.00 1.2e-05 *** -#> s(season):seriesseries_1 3.82 0.12 0.98 -#> s(season):seriesseries_2 3.94 0.09 0.98 -#> s(season):seriesseries_3 3.99 0.74 0.51 +#> s(season) 5.00 16.91 7.4e-06 *** +#> s(season):seriesseries_1 3.90 0.14 0.97 +#> s(season):seriesseries_2 3.90 0.10 0.98 +#> s(season):seriesseries_3 3.96 0.82 0.58 #> --- #> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 #> #> Latent trend marginal deviation (alpha) and length scale (rho) estimates: -#> 2.5% 50% 97.5% Rhat n_eff -#> alpha_gp[1] 0.095 0.42 0.94 1.00 793 -#> alpha_gp[2] 0.370 0.73 1.30 1.00 901 -#> alpha_gp[3] 0.150 0.46 0.94 1.00 899 -#> rho_gp[1] 1.200 3.80 13.00 1.01 451 -#> rho_gp[2] 1.800 7.50 33.00 1.03 254 -#> rho_gp[3] 1.200 4.60 19.00 1.00 1281 +#> 2.5% 50% 97.5% Rhat n_eff +#> alpha_gp[1] 0.10 0.44 0.98 1.00 770 +#> alpha_gp[2] 0.37 0.74 1.30 1.00 1129 +#> alpha_gp[3] 0.16 0.46 0.98 1.00 824 +#> rho_gp[1] 1.20 3.80 14.00 1.00 705 +#> rho_gp[2] 1.80 7.30 33.00 1.01 411 +#> rho_gp[3] 1.30 4.90 23.00 1.01 658 #> #> Stan MCMC diagnostics: #> n_eff / iter looks reasonable for all parameters @@ -626,7 +626,7 @@ summary(mod, include_betas = FALSE) #> 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 Tue Jan 16 8:44:20 PM 2024. +#> Samples were drawn using NUTS(diag_e) at Fri Jan 19 12:34:12 PM 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) @@ -642,11 +642,11 @@ for(i in 1:3){ ``` #> Out of sample CRPS: - #> [1] 2.079879 + #> [1] 2.0731 #> Out of sample CRPS: - #> [1] 1.843941 + #> [1] 1.833157 #> Out of sample CRPS: - #> [1] 1.728574 + #> [1] 1.736085 diff --git a/docs/articles/data_in_mvgam.html b/docs/articles/data_in_mvgam.html index d9fa9525..42f6f2fc 100644 --- a/docs/articles/data_in_mvgam.html +++ b/docs/articles/data_in_mvgam.html @@ -77,7 +77,7 @@

Nicholas J Clark

-

2024-01-12

+

2024-01-19

Source: vignettes/data_in_mvgam.Rmd
data_in_mvgam.Rmd
@@ -106,22 +106,22 @@

Required long data formatsimdat <- sim_mvgam(n_series = 4, T = 24, prop_missing = 0.2) head(simdat$data_train, 16) 
##     y season year   series time
-## 1   0      1    1 series_1    1
-## 2  NA      1    1 series_2    1
-## 3   1      1    1 series_3    1
-## 4  NA      1    1 series_4    1
+## 1   2      1    1 series_1    1
+## 2   3      1    1 series_2    1
+## 3   2      1    1 series_3    1
+## 4   2      1    1 series_4    1
 ## 5   0      2    1 series_1    2
-## 6   0      2    1 series_2    2
-## 7   2      2    1 series_3    2
+## 6   3      2    1 series_2    2
+## 7   3      2    1 series_3    2
 ## 8   0      2    1 series_4    2
-## 9  NA      3    1 series_1    3
+## 9   1      3    1 series_1    3
 ## 10  1      3    1 series_2    3
-## 11 NA      3    1 series_3    3
+## 11  4      3    1 series_3    3
 ## 12  0      3    1 series_4    3
-## 13  0      4    1 series_1    4
+## 13  3      4    1 series_1    4
 ## 14 NA      4    1 series_2    4
 ## 15 NA      4    1 series_3    4
-## 16  0      4    1 series_4    4
+## 16 3 4 1 series_4 4

series as a factor variable @@ -173,22 +173,22 @@

A single outcome variable## ## Coefficients: ## Estimate Std. Error z value Pr(>|z|) -## (Intercept) -1.32292 0.52914 -2.500 0.01241 * -## seriesseries_2 0.93776 0.52643 1.781 0.07486 . -## seriesseries_3 1.30099 0.50364 2.583 0.00979 ** -## seriesseries_4 0.76958 0.53337 1.443 0.14906 -## time 0.03096 0.02863 1.081 0.27959 +## (Intercept) 0.24304 0.27119 0.896 0.37015 +## seriesseries_2 0.38948 0.27346 1.424 0.15438 +## seriesseries_3 0.78087 0.25731 3.035 0.00241 ** +## seriesseries_4 -0.36133 0.32553 -1.110 0.26701 +## time 0.01246 0.01768 0.705 0.48104 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## (Dispersion parameter for poisson family taken to be 1) ## -## Null deviance: 103.369 on 56 degrees of freedom -## Residual deviance: 93.339 on 52 degrees of freedom -## (15 observations deleted due to missingness) -## AIC: 162.91 +## Null deviance: 161.62 on 60 degrees of freedom +## Residual deviance: 140.99 on 56 degrees of freedom +## (11 observations deleted due to missingness) +## AIC: 258.79 ## -## Number of Fisher Scoring iterations: 6 +## Number of Fisher Scoring iterations: 5 
 summary(gam(y ~ series + s(time, by = series),
             data = simdat$data_train,
@@ -202,22 +202,24 @@ 
A single outcome variable## 
 ## Parametric coefficients:
 ##                Estimate Std. Error z value Pr(>|z|)  
-## (Intercept)     -1.0260     0.4503  -2.279   0.0227 *
-## seriesseries_2  -1.1485     1.4713  -0.781   0.4350  
-## seriesseries_3   1.1179     0.5259   2.126   0.0335 *
-## seriesseries_4  -8.3112    24.6052  -0.338   0.7355  
+## (Intercept)      0.3629     0.2086   1.740   0.0819 .
+## seriesseries_2  -1.5109     1.4864  -1.016   0.3094  
+## seriesseries_3  -0.7010     0.6697  -1.047   0.2953  
+## seriesseries_4  -0.3985     0.3330  -1.197   0.2314  
 ## ---
 ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
 ## 
 ## Approximate significance of smooth terms:
-##                          edf Ref.df Chi.sq p-value
-## s(time):seriesseries_1 1.000  1.000  0.336   0.562
-## s(time):seriesseries_2 7.460  8.210  8.495   0.425
-## s(time):seriesseries_3 4.586  5.577  7.662   0.225
-## s(time):seriesseries_4 6.926  7.305  7.461   0.399
+##                          edf Ref.df Chi.sq  p-value    
+## s(time):seriesseries_1 1.000  1.000  0.036 0.850572    
+## s(time):seriesseries_2 8.577  8.877 15.670 0.092017 .  
+## s(time):seriesseries_3 8.543  8.933 30.274 0.000374 ***
+## s(time):seriesseries_4 1.298  1.537  1.230 0.539553    
+## ---
+## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
 ## 
-## R-sq.(adj) =  0.715   Deviance explained =   74%
-## UBRE = 0.31255  Scale est. = 1         n = 57
+## R-sq.(adj) = 0.715 Deviance explained = 65.5% +## UBRE = 0.68285 Scale est. = 1 n = 61

Depending on the observation families you plan to use when building models, there may be some restrictions that need to be satisfied within the outcome variable. For example, a Beta regression can only handle @@ -239,17 +241,17 @@

A single outcome variable levels = 'series1'), time = 1:10) gauss_dat

-
##       outcome  series time
-## 1  -0.5283394 series1    1
-## 2   1.0519794 series1    2
-## 3   0.5998988 series1    3
-## 4   0.6118495 series1    4
-## 5   1.1667473 series1    5
-## 6   0.4938658 series1    6
-## 7  -0.1106668 series1    7
-## 8  -0.7518892 series1    8
-## 9  -2.3121463 series1    9
-## 10 -0.2919469 series1   10
+
##          outcome  series time
+## 1   0.1660605580 series1    1
+## 2   1.5241908906 series1    2
+## 3   0.2700075110 series1    3
+## 4   1.8622153994 series1    4
+## 5  -1.2011344802 series1    5
+## 6  -1.1331587806 series1    6
+## 7   0.0005641071 series1    7
+## 8  -0.6097032999 series1    8
+## 9   1.2220741287 series1    9
+## 10  1.8483330116 series1   10

A call to gam using the mgcv package leads to a model that actually fits (though it does give an unhelpful warning message):

@@ -260,14 +262,14 @@

A single outcome variable## Warning in family$saturated.ll(y, prior.weights, theta): saturated likelihood ## may be inaccurate 
## 
-## Family: Beta regression(0.116) 
+## Family: Beta regression(0.087) 
 ## Link function: logit 
 ## 
 ## Formula:
 ## outcome ~ time
 ## Total model degrees of freedom 2 
 ## 
-## REML score: -176.7579
+## REML score: -182.0205

But the same call to mvgam gives us something more useful:

 
##   time   series    outcome
-## 1    1 series_1 -0.3103526
-## 2    3 series_1 -1.0459039
-## 3    5 series_1 -1.3013471
-## 4    7 series_1 -0.5744845
-## 5    9 series_1  0.1146766
-## 6   11 series_1 -0.4771703
-## 7   13 series_1  1.4805867
-## 8   15 series_1  0.1176757
+## 1 1 series_1 -0.6853424 +## 2 3 series_1 -1.4364053 +## 3 5 series_1 0.2405592 +## 4 7 series_1 -0.7739931 +## 5 9 series_1 0.3248164 +## 6 11 series_1 -0.2119142 +## 7 13 series_1 0.9151188 +## 8 15 series_1 1.0588870

Next we call get_mvgam_priors by simply specifying an intercept-only model, which is enough to trigger all the checks:

@@ -375,21 +377,21 @@ 
Checking data with get_mvgam_
 
 good_times

 
##    time   series    outcome
-## 1     1 series_1 -0.3103526
+## 1     1 series_1 -0.6853424
 ## 2     2 series_1         NA
-## 3     3 series_1 -1.0459039
+## 3     3 series_1 -1.4364053
 ## 4     4 series_1         NA
-## 5     5 series_1 -1.3013471
+## 5     5 series_1  0.2405592
 ## 6     6 series_1         NA
-## 7     7 series_1 -0.5744845
+## 7     7 series_1 -0.7739931
 ## 8     8 series_1         NA
-## 9     9 series_1  0.1146766
+## 9     9 series_1  0.3248164
 ## 10   10 series_1         NA
-## 11   11 series_1 -0.4771703
+## 11   11 series_1 -0.2119142
 ## 12   12 series_1         NA
-## 13   13 series_1  1.4805867
+## 13   13 series_1  0.9151188
 ## 14   14 series_1         NA
-## 15   15 series_1  0.1176757

+## 15   15 series_1  1.0588870

Now the call to get_mvgam_priors, using our filled in data, should work:

@@ -399,9 +401,9 @@ 
Checking data with get_mvgam_
 
##                             param_name param_length           param_info
 ## 1                          (Intercept)            1          (Intercept)
 ## 2 vector<lower=0>[n_series] sigma_obs;            1 observation error sd
-##                                    prior                  example_change
-## 1 (Intercept) ~ student_t(3, -0.4, 2.5);     (Intercept) ~ normal(0, 1);
-## 2      sigma_obs ~ student_t(3, 0, 2.5); sigma_obs ~ normal(0.71, 0.86);
+##                                 prior                  example_change
+## 1 (Intercept) ~ student_t(3, 0, 2.5);     (Intercept) ~ normal(0, 1);
+## 2   sigma_obs ~ student_t(3, 0, 2.5); sigma_obs ~ normal(0.93, 0.33);
 ##   new_lowerbound new_upperbound
 ## 1             NA             NA
 ## 2             NA             NA

@@ -449,9 +451,9 @@ 
Checking data with get_mvgam_
 
##                             param_name param_length           param_info
 ## 1                          (Intercept)            1          (Intercept)
 ## 2 vector<lower=0>[n_series] sigma_obs;            1 observation error sd
-##                                   prior                   example_change
-## 1 (Intercept) ~ student_t(3, 0.7, 2.5);      (Intercept) ~ normal(0, 1);
-## 2     sigma_obs ~ student_t(3, 0, 2.5); sigma_obs ~ normal(-0.41, 0.22);
+##                                    prior               example_change
+## 1 (Intercept) ~ student_t(3, -0.8, 2.5);  (Intercept) ~ normal(0, 1);
+## 2      sigma_obs ~ student_t(3, 0, 2.5); sigma_obs ~ normal(0, 0.19);
 ##   new_lowerbound new_upperbound
 ## 1             NA             NA
 ## 2             NA             NA

@@ -474,23 +476,23 @@ 
Covariates with no NAs= 'series1'),
                        time = 1:10)
 miss_dat
-
##       outcome          cov  series time
-## 1  -0.8073610           NA series1    1
-## 2   1.4759469 -1.479350061 series1    2
-## 3  -2.3326864  0.623521555 series1    3
-## 4  -0.1135510 -1.072257096 series1    4
-## 5  -1.7327234 -0.631452996 series1    5
-## 6  -1.2263422  0.508258203 series1    6
-## 7  -1.9541538  0.000188825 series1    7
-## 8   0.8113813  0.595749976 series1    8
-## 9   1.0077533  0.438257074 series1    9
-## 10 -0.5593047  1.096516533 series1   10
+
##        outcome         cov  series time
+## 1   0.52094669          NA series1    1
+## 2  -0.58644030  0.73867095 series1    2
+## 3   1.28310150 -0.58577764 series1    3
+## 4  -0.35030777 -1.07774912 series1    4
+## 5   0.82824135 -0.26907960 series1    5
+## 6   0.18857106  1.65804013 series1    6
+## 7  -0.36209231 -0.67063650 series1    7
+## 8  -0.02767463  2.52337077 series1    8
+## 9   0.63317969  0.02399786 series1    9
+## 10  1.97227522 -0.41125516 series1   10
 get_mvgam_priors(outcome ~ cov,
                  data = miss_dat,
                  family = gaussian())
## Error: Missing values found in data predictors:
-##  Error in na.fail.default(structure(list(outcome = c(-0.807361010516451, : missing values in object
+## Error in na.fail.default(structure(list(outcome = c(0.520946693584384, : missing values in object

Just like with the mgcv package, mvgam can also accept data as a list object. This is useful if you want to set up linear @@ -510,7 +512,7 @@

Covariates with no NAs= miss_dat, family = gaussian())

## Error: Missing values found in data predictors:
-##  Error in na.fail.default(structure(list(outcome = c(-0.548669693991903, : missing values in object
+## Error in na.fail.default(structure(list(outcome = c(-2.90806279836404, : missing values in object

@@ -711,12 +713,12 @@

Example with NEON tick data## $ S6 : num [1:8, 1:8] 1.037 -0.416 0.419 0.117 0.188 ... ## $ S7 : num [1:8, 1:8] 1.037 -0.416 0.419 0.117 0.188 ... ## $ S8 : num [1:8, 1:8] 1.037 -0.416 0.419 0.117 0.188 ... -## $ p_coefs : Named num 0.809 +## $ p_coefs : Named num 0.659 ## ..- attr(*, "names")= chr "(Intercept)" -## $ p_taus : num 301 +## $ p_taus : num 387 ## $ ytimes : int [1:416, 1:8] 1 9 17 25 33 41 49 57 65 73 ... ## $ n_series : int 8 -## $ sp : Named num [1:9] 12.38 4.42 400.88 1.51 11.13 ... +## $ sp : Named num [1:9] 1.19 4.84 4.21 16.59 2.61 ... ## ..- attr(*, "names")= chr [1:9] "s(epiWeek):seriesBLAN_005" "s(epiWeek):seriesBLAN_012" "s(epiWeek):seriesSCBI_002" "s(epiWeek):seriesSCBI_013" ... ## $ y_observed : num [1:416, 1:8] 0 0 0 0 0 0 0 0 0 0 ... ## $ total_obs : int 3328 diff --git a/docs/articles/data_in_mvgam_files/figure-html/unnamed-chunk-23-1.png b/docs/articles/data_in_mvgam_files/figure-html/unnamed-chunk-23-1.png index 2f138a66..3847029f 100644 Binary files a/docs/articles/data_in_mvgam_files/figure-html/unnamed-chunk-23-1.png and b/docs/articles/data_in_mvgam_files/figure-html/unnamed-chunk-23-1.png differ diff --git a/docs/articles/data_in_mvgam_files/figure-html/unnamed-chunk-24-1.png b/docs/articles/data_in_mvgam_files/figure-html/unnamed-chunk-24-1.png index bddc1d01..c5e649dd 100644 Binary files a/docs/articles/data_in_mvgam_files/figure-html/unnamed-chunk-24-1.png and b/docs/articles/data_in_mvgam_files/figure-html/unnamed-chunk-24-1.png differ diff --git a/docs/articles/data_in_mvgam_files/figure-html/unnamed-chunk-25-1.png b/docs/articles/data_in_mvgam_files/figure-html/unnamed-chunk-25-1.png index 00e1e949..aa99b673 100644 Binary files a/docs/articles/data_in_mvgam_files/figure-html/unnamed-chunk-25-1.png and b/docs/articles/data_in_mvgam_files/figure-html/unnamed-chunk-25-1.png differ diff --git a/docs/articles/forecast_evaluation.html b/docs/articles/forecast_evaluation.html index f356addf..2eaccb45 100644 --- a/docs/articles/forecast_evaluation.html +++ b/docs/articles/forecast_evaluation.html @@ -77,7 +77,7 @@

Nicholas J Clark

-

2024-01-12

+

2024-01-19

Source: vignettes/forecast_evaluation.Rmd
forecast_evaluation.Rmd
@@ -211,15 +211,15 @@

Modelling dynamics with splines## ## ## GAM coefficient (beta) estimates: -## 2.5% 50% 97.5% Rhat n_eff -## (Intercept) -0.42 -0.22 -0.03 1.01 847 +## 2.5% 50% 97.5% Rhat n_eff +## (Intercept) -0.41 -0.21 -0.039 1 813 ## ## Approximate significance of GAM observation smooths: -## edf Chi.sq p-value -## s(season) 3.71 9.25 0.0191 * -## s(time):seriesseries_1 7.55 13.18 0.0919 . -## s(time):seriesseries_2 11.05 346.28 0.0023 ** -## s(time):seriesseries_3 6.24 16.78 0.0384 * +## edf Chi.sq p-value +## s(season) 3.77 9.48 0.01603 * +## s(time):seriesseries_1 6.50 13.64 0.09218 . +## s(time):seriesseries_2 9.49 256.09 0.00021 *** +## s(time):seriesseries_3 5.93 16.79 0.04680 * ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## @@ -230,7 +230,7 @@

Modelling dynamics with splines## 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 Fri Jan 12 3:49:18 PM 2024. +## Samples were drawn using NUTS(diag_e) at Fri Jan 19 11:52:33 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) @@ -287,32 +287,32 @@

Modelling dynamics with GPs## ## GAM coefficient (beta) estimates: ## 2.5% 50% 97.5% Rhat n_eff -## (Intercept) -1.1 -0.53 0.19 1 700 +## (Intercept) -1.1 -0.52 0.31 1 768 ## ## GAM gp term marginal deviation (alpha) and length scale (rho) estimates: -## 2.5% 50% 97.5% Rhat n_eff -## alpha_gp(time):seriesseries_1 0.24 0.77 2.1 1.00 1090 -## alpha_gp(time):seriesseries_2 0.75 1.40 3.0 1.01 900 -## alpha_gp(time):seriesseries_3 0.48 1.10 2.5 1.00 900 -## rho_gp(time):seriesseries_1 1.20 4.90 20.0 1.00 750 -## rho_gp(time):seriesseries_2 2.20 11.00 17.0 1.00 512 -## rho_gp(time):seriesseries_3 1.40 8.20 23.0 1.00 825 +## 2.5% 50% 97.5% Rhat n_eff +## alpha_gp(time):seriesseries_1 0.21 0.8 2.1 1.01 763 +## alpha_gp(time):seriesseries_2 0.74 1.4 2.9 1.00 1028 +## alpha_gp(time):seriesseries_3 0.50 1.1 2.8 1.00 1026 +## rho_gp(time):seriesseries_1 1.20 5.1 23.0 1.00 681 +## rho_gp(time):seriesseries_2 2.20 10.0 17.0 1.00 644 +## rho_gp(time):seriesseries_3 1.50 8.8 23.0 1.00 819 ## ## Approximate significance of GAM observation smooths: ## edf Chi.sq p-value -## s(season) 6 24.8 0.00028 *** +## s(season) 6 25 0.00016 *** ## --- ## 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 -## 5 of 2000 iterations ended with a divergence (0.25%) +## 4 of 2000 iterations ended with a divergence (0.2%) ## *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 ## -## Samples were drawn using NUTS(diag_e) at Fri Jan 12 3:50:20 PM 2024. +## Samples were drawn using NUTS(diag_e) at Fri Jan 19 11:53:19 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) @@ -395,22 +395,22 @@

Forecasting with the forec ## ..$ series_3: int [1:25] 1 0 0 1 0 0 1 0 1 0 ... ## $ test_times : int [1:25] 76 77 78 79 80 81 82 83 84 85 ... ## $ hindcasts :List of 3 -## ..$ series_1: num [1:2000, 1:75] 1 2 1 1 2 0 0 1 0 0 ... +## ..$ series_1: num [1:2000, 1:75] 1 1 0 0 0 1 1 1 0 0 ... ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : NULL ## .. .. ..$ : chr [1:75] "ypred[1,1]" "ypred[2,1]" "ypred[3,1]" "ypred[4,1]" ... -## ..$ series_2: num [1:2000, 1:75] 0 0 0 0 0 0 0 0 0 0 ... +## ..$ series_2: num [1:2000, 1:75] 0 0 0 0 0 0 0 1 0 0 ... ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : NULL ## .. .. ..$ : chr [1:75] "ypred[1,2]" "ypred[2,2]" "ypred[3,2]" "ypred[4,2]" ... -## ..$ series_3: num [1:2000, 1:75] 0 1 3 1 0 0 0 0 3 5 ... +## ..$ series_3: num [1:2000, 1:75] 3 0 2 1 0 1 2 1 5 1 ... ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : NULL ## .. .. ..$ : chr [1:75] "ypred[1,3]" "ypred[2,3]" "ypred[3,3]" "ypred[4,3]" ... ## $ forecasts :List of 3 -## ..$ series_1: num [1:2000, 1:25] 1 0 0 0 0 1 1 0 1 0 ... -## ..$ series_2: num [1:2000, 1:25] 3 1 1 1 1 0 1 0 5 0 ... -## ..$ series_3: num [1:2000, 1:25] 0 3 2 1 0 0 2 2 1 0 ... +## ..$ series_1: num [1:2000, 1:25] 1 3 2 1 0 0 1 1 0 0 ... +## ..$ series_2: num [1:2000, 1:25] 6 0 0 0 0 2 0 0 0 0 ... +## ..$ series_3: num [1:2000, 1:25] 0 1 1 3 3 1 3 2 4 2 ... ## - attr(*, "class")= chr "mvgam_forecast"

We can plot the forecasts for each series from each model using the S3 plot method for objects of this class:

@@ -418,32 +418,32 @@

Forecasting with the forec plot(fc_mod1, series = 1)

## Out of sample CRPS:
-## [1] 13.96002
+## [1] 14.62964 
 plot(fc_mod2, series = 1)

## Out of sample DRPS:
-## [1] 10.79944
+## [1] 10.92516 
 plot(fc_mod1, series = 2)

## Out of sample CRPS:
-## [1] 4.209141e+12
+## [1] 84201962708 
 plot(fc_mod2, series = 2)

## Out of sample DRPS:
-## [1] 14.30626
+## [1] 14.31168 
 plot(fc_mod1, series = 3)

## Out of sample CRPS:
-## [1] 41.10631
+## [1] 32.44136 
 plot(fc_mod2, series = 3)

## Out of sample DRPS:
-## [1] 15.62831
+## [1] 15.44332

Clearly the two models do not produce equivalent forecasts. We will come back to scoring these forecasts in a moment.

@@ -477,7 +477,7 @@

Forecasting with newdata plot(fc_mod2, series = 1) 

## Out of sample DRPS:
-## [1] 10.75464
+## [1] 10.85762

Scoring forecast distributions @@ -494,55 +494,55 @@

Scoring forecast distributionsstr(crps_mod1)

## List of 4
 ##  $ series_1  :'data.frame':  25 obs. of  5 variables:
-##   ..$ score         : num [1:25] 0.1975 0.1398 1.3519 NA 0.0359 ...
+##   ..$ score         : num [1:25] 0.1938 0.1366 1.355 NA 0.0348 ...
 ##   ..$ in_interval   : num [1:25] 1 1 1 NA 1 1 1 1 1 1 ...
 ##   ..$ interval_width: num [1:25] 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 ...
 ##   ..$ eval_horizon  : int [1:25] 1 2 3 4 5 6 7 8 9 10 ...
 ##   ..$ score_type    : chr [1:25] "crps" "crps" "crps" "crps" ...
 ##  $ series_2  :'data.frame':  25 obs. of  5 variables:
-##   ..$ score         : num [1:25] 0.362 0.349 0.941 0.492 0.558 ...
+##   ..$ score         : num [1:25] 0.379 0.306 0.941 0.5 0.573 ...
 ##   ..$ in_interval   : num [1:25] 1 1 1 1 1 1 1 1 1 NA ...
 ##   ..$ interval_width: num [1:25] 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 ...
 ##   ..$ eval_horizon  : int [1:25] 1 2 3 4 5 6 7 8 9 10 ...
 ##   ..$ score_type    : chr [1:25] "crps" "crps" "crps" "crps" ...
 ##  $ series_3  :'data.frame':  25 obs. of  5 variables:
-##   ..$ score         : num [1:25] 0.329 0.597 0.386 0.349 0.235 ...
+##   ..$ score         : num [1:25] 0.32 0.556 0.379 0.362 0.219 ...
 ##   ..$ in_interval   : num [1:25] 1 1 1 1 1 1 1 1 1 1 ...
 ##   ..$ interval_width: num [1:25] 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 ...
 ##   ..$ eval_horizon  : int [1:25] 1 2 3 4 5 6 7 8 9 10 ...
 ##   ..$ score_type    : chr [1:25] "crps" "crps" "crps" "crps" ...
 ##  $ all_series:'data.frame':  25 obs. of  3 variables:
-##   ..$ score       : num [1:25] 0.888 1.086 2.68 NA 0.829 ...
+##   ..$ score       : num [1:25] 0.892 0.999 2.675 NA 0.827 ...
 ##   ..$ eval_horizon: int [1:25] 1 2 3 4 5 6 7 8 9 10 ...
 ##   ..$ score_type  : chr [1:25] "sum_crps" "sum_crps" "sum_crps" "sum_crps" ...
 crps_mod1$series_1
##         score in_interval interval_width eval_horizon score_type
-## 1  0.19749300           1            0.9            1       crps
-## 2  0.13975150           1            0.9            2       crps
-## 3  1.35193650           1            0.9            3       crps
+## 1  0.19375525           1            0.9            1       crps
+## 2  0.13663925           1            0.9            2       crps
+## 3  1.35502175           1            0.9            3       crps
 ## 4          NA          NA            0.9            4       crps
-## 5  0.03594650           1            0.9            5       crps
-## 6  1.55512000           1            0.9            6       crps
-## 7  1.52746400           1            0.9            7       crps
-## 8  0.62147725           1            0.9            8       crps
-## 9  0.60050300           1            0.9            9       crps
-## 10 0.58537625           1            0.9           10       crps
-## 11 1.27717750           1            0.9           11       crps
-## 12 2.06627175           1            0.9           12       crps
-## 13 0.58131925           1            0.9           13       crps
-## 14 0.13444200           1            0.9           14       crps
-## 15 0.63533750           1            0.9           15       crps
-## 16 0.06961475           1            0.9           16       crps
-## 17 0.06879775           1            0.9           17       crps
-## 18 0.07628275           1            0.9           18       crps
-## 19 0.10143425           1            0.9           19       crps
+## 5  0.03482775           1            0.9            5       crps
+## 6  1.55416700           1            0.9            6       crps
+## 7  1.51028900           1            0.9            7       crps
+## 8  0.62121225           1            0.9            8       crps
+## 9  0.62630125           1            0.9            9       crps
+## 10 0.59853100           1            0.9           10       crps
+## 11 1.30998625           1            0.9           11       crps
+## 12 2.04829775           1            0.9           12       crps
+## 13 0.61251800           1            0.9           13       crps
+## 14 0.14052300           1            0.9           14       crps
+## 15 0.65110800           1            0.9           15       crps
+## 16 0.07973125           1            0.9           16       crps
+## 17 0.07675600           1            0.9           17       crps
+## 18 0.09382375           1            0.9           18       crps
+## 19 0.12356725           1            0.9           19       crps
 ## 20         NA          NA            0.9           20       crps
-## 21 0.16723125           1            0.9           21       crps
-## 22 0.76396600           1            0.9           22       crps
+## 21 0.20173600           1            0.9           21       crps
+## 22 0.84066825           1            0.9           22       crps
 ## 23         NA          NA            0.9           23       crps
-## 24 0.90052575           1            0.9           24       crps
-## 25 0.50255325           1            0.9           25       crps
+## 24 1.06489225 1 0.9 24 crps +## 25 0.75528825 1 0.9 25 crps

The returned list contains a data.frame for each series in the data that shows the CRPS score for each evaluation in the testing data, along with some other useful information about the fit of the @@ -557,31 +557,31 @@

Scoring forecast distributionscrps_mod1 <- score(fc_mod1, score = 'crps', interval_width = 0.6) crps_mod1$series_1 
##         score in_interval interval_width eval_horizon score_type
-## 1  0.19749300           1            0.6            1       crps
-## 2  0.13975150           1            0.6            2       crps
-## 3  1.35193650           0            0.6            3       crps
+## 1  0.19375525           1            0.6            1       crps
+## 2  0.13663925           1            0.6            2       crps
+## 3  1.35502175           0            0.6            3       crps
 ## 4          NA          NA            0.6            4       crps
-## 5  0.03594650           1            0.6            5       crps
-## 6  1.55512000           0            0.6            6       crps
-## 7  1.52746400           0            0.6            7       crps
-## 8  0.62147725           1            0.6            8       crps
-## 9  0.60050300           1            0.6            9       crps
-## 10 0.58537625           1            0.6           10       crps
-## 11 1.27717750           0            0.6           11       crps
-## 12 2.06627175           0            0.6           12       crps
-## 13 0.58131925           1            0.6           13       crps
-## 14 0.13444200           1            0.6           14       crps
-## 15 0.63533750           1            0.6           15       crps
-## 16 0.06961475           1            0.6           16       crps
-## 17 0.06879775           1            0.6           17       crps
-## 18 0.07628275           1            0.6           18       crps
-## 19 0.10143425           1            0.6           19       crps
+## 5  0.03482775           1            0.6            5       crps
+## 6  1.55416700           0            0.6            6       crps
+## 7  1.51028900           0            0.6            7       crps
+## 8  0.62121225           1            0.6            8       crps
+## 9  0.62630125           1            0.6            9       crps
+## 10 0.59853100           1            0.6           10       crps
+## 11 1.30998625           0            0.6           11       crps
+## 12 2.04829775           0            0.6           12       crps
+## 13 0.61251800           1            0.6           13       crps
+## 14 0.14052300           1            0.6           14       crps
+## 15 0.65110800           1            0.6           15       crps
+## 16 0.07973125           1            0.6           16       crps
+## 17 0.07675600           1            0.6           17       crps
+## 18 0.09382375           1            0.6           18       crps
+## 19 0.12356725           1            0.6           19       crps
 ## 20         NA          NA            0.6           20       crps
-## 21 0.16723125           1            0.6           21       crps
-## 22 0.76396600           1            0.6           22       crps
+## 21 0.20173600           1            0.6           21       crps
+## 22 0.84066825           1            0.6           22       crps
 ## 23         NA          NA            0.6           23       crps
-## 24 0.90052575           1            0.6           24       crps
-## 25 0.50255325           1            0.6           25       crps
+## 24 1.06489225 1 0.6 24 crps +## 25 0.75528825 1 0.6 25 crps

We can also compare forecasts against out of sample observations using the Expected Log Predictive Density (ELPD; also known as the log score). The ELPD is a strictly proper scoring rule that can be @@ -593,31 +593,31 @@

Scoring forecast distributionslink_mod1 <- forecast(mod1, newdata = simdat$data_test, type = 'link') score(link_mod1, score = 'elpd')$series_1 
##         score eval_horizon score_type
-## 1  -0.5415389            1       elpd
-## 2  -0.4394250            2       elpd
-## 3  -2.9516572            3       elpd
+## 1  -0.5304414            1       elpd
+## 2  -0.4298955            2       elpd
+## 3  -2.9617583            3       elpd
 ## 4          NA            4       elpd
-## 5  -0.2039540            5       elpd
-## 6  -3.3591450            6       elpd
-## 7  -3.2434347            7       elpd
-## 8  -1.9910208            8       elpd
-## 9  -2.0099446            9       elpd
-## 10 -2.0168619           10       elpd
-## 11 -3.0178822           11       elpd
-## 12 -3.6389118           12       elpd
-## 13 -2.0934546           13       elpd
-## 14 -0.3015793           14       elpd
-## 15 -2.2982746           15       elpd
-## 16 -0.2153677           16       elpd
-## 17 -0.2021458           17       elpd
-## 18 -0.2102088           18       elpd
-## 19 -0.2230817           19       elpd
+## 5  -0.2007644            5       elpd
+## 6  -3.3781408            6       elpd
+## 7  -3.2729088            7       elpd
+## 8  -2.0363750            8       elpd
+## 9  -2.0670612            9       elpd
+## 10 -2.0844818           10       elpd
+## 11 -3.0576463           11       elpd
+## 12 -3.6291058           12       elpd
+## 13 -2.1692669           13       elpd
+## 14 -0.2960899           14       elpd
+## 15 -2.3738851           15       elpd
+## 16 -0.2160804           16       elpd
+## 17 -0.2036782           17       elpd
+## 18 -0.2115539           18       elpd
+## 19 -0.2235072           19       elpd
 ## 20         NA           20       elpd
-## 21 -0.2440297           21       elpd
-## 22 -2.5577404           22       elpd
+## 21 -0.2413680           21       elpd
+## 22 -2.6791984           22       elpd
 ## 23         NA           23       elpd
-## 24 -2.5796575           24       elpd
-## 25 -0.2890367           25       elpd
+## 24 -2.6851981 24 elpd +## 25 -0.2836901 25 elpd

Finally, when we have multiple time series it may also make sense to use a multivariate proper scoring rule. mvgam offers two such options: the Energy score and the Variogram score. The first @@ -642,7 +642,7 @@

Scoring forecast distributions## ..$ interval_width: num [1:25] 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 ... ## ..$ eval_horizon : int [1:25] 1 2 3 4 5 6 7 8 9 10 ... ## $ all_series:'data.frame': 25 obs. of 3 variables: -## ..$ score : num [1:25] 0.77 1.115 1.25 NA 0.453 ... +## ..$ score : num [1:25] 0.773 1.147 1.226 NA 0.458 ... ## ..$ eval_horizon: int [1:25] 1 2 3 4 5 6 7 8 9 10 ... ## ..$ score_type : chr [1:25] "energy" "energy" "energy" "energy" ...

The returned object still provides information on interval coverage @@ -651,31 +651,31 @@

Scoring forecast distributions
 energy_mod2$all_series
##        score eval_horizon score_type
-## 1  0.7702021            1     energy
-## 2  1.1146337            2     energy
-## 3  1.2502919            3     energy
+## 1  0.7728517            1     energy
+## 2  1.1469836            2     energy
+## 3  1.2258781            3     energy
 ## 4         NA            4     energy
-## 5  0.4534938            5     energy
-## 6  1.7972077            6     energy
-## 7  1.4989659            7     energy
-## 8  0.7437805            8     energy
-## 9  1.0982204            9     energy
+## 5  0.4577536            5     energy
+## 6  1.8094487            6     energy
+## 7  1.4887317            7     energy
+## 8  0.7651593            8     energy
+## 9  1.1180634            9     energy
 ## 10        NA           10     energy
-## 11 1.4463652           11     energy
-## 12 3.3035827           12     energy
-## 13 1.5265317           13     energy
-## 14 1.1512861           14     energy
-## 15 1.0331594           15     energy
-## 16 1.8705236           16     energy
+## 11 1.5008324           11     energy
+## 12 3.2142460           12     energy
+## 13 1.6129732           13     energy
+## 14 1.2704438           14     energy
+## 15 1.1335958           15     energy
+## 16 1.8717420           16     energy
 ## 17        NA           17     energy
-## 18 0.7118636           18     energy
-## 19 0.8919367           19     energy
+## 18 0.7953392           18     energy
+## 19 0.9919119           19     energy
 ## 20        NA           20     energy
-## 21 1.0178829           21     energy
-## 22 1.3557390           22     energy
+## 21 1.2461964           21     energy
+## 22 1.5170615           22     energy
 ## 23        NA           23     energy
-## 24 2.1700795           24     energy
-## 25 1.2258707           25     energy
+## 24 2.3824552 24 energy +## 25 1.5314557 25 energy

You can use your score(s) of choice to compare different models. For example, we can compute and plot the difference in CRPS scores for each series in data. Here, a negative value means the Gaussian Process model diff --git a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-13-1.png b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-13-1.png index 6a14c23c..edac75d5 100644 Binary files a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-13-1.png and b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-13-1.png differ diff --git a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-14-1.png b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-14-1.png index 367dcbe0..db1be0c9 100644 Binary files a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-14-1.png and b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-14-1.png differ diff --git a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-15-1.png b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-15-1.png index 01340e0d..7b709ac3 100644 Binary files a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-15-1.png and b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-15-1.png differ diff --git a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-16-1.png b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-16-1.png index 52168bc2..59cffcc2 100644 Binary files a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-16-1.png and b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-16-1.png differ diff --git a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-19-1.png b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-19-1.png index 5c7f7c1a..0bf0db70 100644 Binary files a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-19-1.png and b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-19-1.png differ diff --git a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-19-2.png b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-19-2.png index 1b252cd5..b27b2b90 100644 Binary files a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-19-2.png and b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-19-2.png differ diff --git a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-19-3.png b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-19-3.png index aaad97ad..a70e489b 100644 Binary files a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-19-3.png and b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-19-3.png differ diff --git a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-19-4.png b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-19-4.png index 7f236548..d471c24e 100644 Binary files a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-19-4.png and b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-19-4.png differ diff --git a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-19-5.png b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-19-5.png index 2d5cc59f..eb6f009f 100644 Binary files a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-19-5.png and b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-19-5.png differ diff --git a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-19-6.png b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-19-6.png index 3ee70081..4be6fa4f 100644 Binary files a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-19-6.png and b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-19-6.png differ diff --git a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-23-1.png b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-23-1.png index 0463aaf0..3769cdfb 100644 Binary files a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-23-1.png and b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-23-1.png differ diff --git a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-29-1.png b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-29-1.png index 31a7f559..3bb928b4 100644 Binary files a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-29-1.png and b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-29-1.png differ diff --git a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-29-2.png b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-29-2.png index 6b846d85..2c74644e 100644 Binary files a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-29-2.png and b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-29-2.png differ diff --git a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-29-3.png b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-29-3.png index 3ef5dcd8..6b9a196a 100644 Binary files a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-29-3.png and b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-29-3.png differ diff --git a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-9-1.png b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-9-1.png index f299083b..47909451 100644 Binary files a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-9-1.png and b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-9-1.png differ diff --git a/docs/articles/mvgam_overview.html b/docs/articles/mvgam_overview.html index 71b41a49..dee0648c 100644 --- a/docs/articles/mvgam_overview.html +++ b/docs/articles/mvgam_overview.html @@ -77,7 +77,7 @@

Nicholas J Clark

-

2024-01-15

+

2024-01-19

Source: vignettes/mvgam_overview.Rmd
mvgam_overview.Rmd
@@ -444,19 +444,19 @@

Manipulating data for modelling
 data <- sim_mvgam(n_series = 4, T = 24)
 head(data$data_train, 12)
-
##    y season year   series time
-## 1  1      1    1 series_1    1
-## 2  7      1    1 series_2    1
-## 3  8      1    1 series_3    1
-## 4  1      1    1 series_4    1
-## 5  2      2    1 series_1    2
-## 6  5      2    1 series_2    2
-## 7  1      2    1 series_3    2
-## 8  2      2    1 series_4    2
-## 9  1      3    1 series_1    3
-## 10 0      3    1 series_2    3
-## 11 1      3    1 series_3    3
-## 12 0      3    1 series_4    3
+
##     y season year   series time
+## 1   1      1    1 series_1    1
+## 2   6      1    1 series_2    1
+## 3  12      1    1 series_3    1
+## 4   0      1    1 series_4    1
+## 5   3      2    1 series_1    2
+## 6   5      2    1 series_2    2
+## 7   3      2    1 series_3    2
+## 8   0      2    1 series_4    2
+## 9   2      3    1 series_1    3
+## 10  4      3    1 series_2    3
+## 11  3      3    1 series_3    3
+## 12  0      3    1 series_4    3

Notice how we have four different time series in these simulated data, but we do not spread the outcome values into different columns. Rather, there is only a single column for the outcome variable, labelled @@ -649,8 +649,8 @@

GLMs with temporal random effects## 1 vector[1] mu_raw; 1 s(year_fac) pop mean ## 2 vector<lower=0>[1] sigma_raw; 1 s(year_fac) pop sd ## prior example_change -## 1 mu_raw ~ std_normal(); mu_raw ~ normal(-0.6, 0.26); -## 2 sigma_raw ~ student_t(3, 0, 2.5); sigma_raw ~ exponential(0.52); +## 1 mu_raw ~ std_normal(); mu_raw ~ normal(0.21, 0.23); +## 2 sigma_raw ~ student_t(3, 0, 2.5); sigma_raw ~ exponential(0.29); ## new_lowerbound new_upperbound ## 1 NA NA ## 2 NA NA @@ -686,33 +686,33 @@

GLMs with temporal random effects## ## ## GAM coefficient (beta) estimates: -## 2.5% 50% 97.5% Rhat n_eff -## s(year_fac).1 1.8 2.1 2.3 1.00 2439 -## s(year_fac).2 2.5 2.7 2.8 1.00 3112 -## s(year_fac).3 3.0 3.1 3.2 1.00 3182 -## s(year_fac).4 3.1 3.3 3.4 1.00 2476 -## s(year_fac).5 1.9 2.1 2.3 1.00 2944 -## s(year_fac).6 1.5 1.8 2.0 1.00 3386 -## s(year_fac).7 1.8 2.0 2.3 1.00 2740 -## s(year_fac).8 2.8 3.0 3.1 1.00 2897 -## s(year_fac).9 3.1 3.2 3.4 1.00 2916 -## s(year_fac).10 2.6 2.8 2.9 1.00 2911 -## s(year_fac).11 2.9 3.1 3.2 1.00 2772 -## s(year_fac).12 3.1 3.2 3.3 1.00 3060 -## s(year_fac).13 2.0 2.2 2.5 1.00 2228 -## s(year_fac).14 2.5 2.6 2.8 1.00 3100 -## s(year_fac).15 1.9 2.2 2.4 1.00 2161 -## s(year_fac).16 1.9 2.1 2.3 1.00 3331 -## s(year_fac).17 -0.3 1.1 1.9 1.01 532 +## 2.5% 50% 97.5% Rhat n_eff +## s(year_fac).1 1.80 2.1 2.3 1.00 2678 +## s(year_fac).2 2.50 2.7 2.8 1.00 2668 +## s(year_fac).3 3.00 3.1 3.2 1.00 3090 +## s(year_fac).4 3.10 3.3 3.4 1.00 2636 +## s(year_fac).5 1.90 2.1 2.3 1.00 2405 +## s(year_fac).6 1.50 1.8 2.0 1.00 2637 +## s(year_fac).7 1.80 2.0 2.3 1.00 2381 +## s(year_fac).8 2.80 3.0 3.1 1.00 2866 +## s(year_fac).9 3.10 3.2 3.4 1.00 2801 +## s(year_fac).10 2.60 2.8 2.9 1.00 2844 +## s(year_fac).11 3.00 3.1 3.2 1.00 3150 +## s(year_fac).12 3.10 3.2 3.3 1.00 2307 +## s(year_fac).13 2.00 2.2 2.5 1.00 2857 +## s(year_fac).14 2.50 2.6 2.8 1.00 2939 +## s(year_fac).15 1.90 2.2 2.4 1.00 2610 +## s(year_fac).16 1.90 2.1 2.3 1.00 3129 +## s(year_fac).17 -0.28 1.1 1.9 1.02 367 ## ## GAM group-level estimates: ## 2.5% 50% 97.5% Rhat n_eff -## mean(s(year_fac)) 2.10 2.40 2.8 1.03 198 -## sd(s(year_fac)) 0.45 0.66 1.1 1.01 238 +## mean(s(year_fac)) 2.00 2.40 2.7 1.02 190 +## sd(s(year_fac)) 0.46 0.68 1.1 1.03 203 ## ## Approximate significance of GAM observation smooths: ## edf Chi.sq p-value -## s(year_fac) 13.2 24457 <2e-16 *** +## s(year_fac) 12.8 23511 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## @@ -723,7 +723,7 @@

GLMs with temporal random effects## 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 Mon Jan 15 9:06:57 AM 2024. +## Samples were drawn using NUTS(diag_e) at Fri Jan 19 11:54:28 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) @@ -740,23 +740,23 @@

GLMs with temporal random effectsdplyr::glimpse(beta_post) 
## Rows: 2,000
 ## Columns: 17
-## $ `s(year_fac).1`  <dbl> 2.30080, 2.19798, 2.04816, 2.01919, 2.19020, 1.91288,…
-## $ `s(year_fac).2`  <dbl> 2.72378, 2.69498, 2.64144, 2.73021, 2.80479, 2.63429,…
-## $ `s(year_fac).3`  <dbl> 3.04055, 2.94148, 3.27943, 3.05446, 3.13977, 3.14182,…
-## $ `s(year_fac).4`  <dbl> 3.31494, 3.28852, 3.19607, 3.15786, 3.12430, 3.15365,…
-## $ `s(year_fac).5`  <dbl> 2.24893, 2.16990, 1.99415, 2.11878, 2.20520, 2.09616,…
-## $ `s(year_fac).6`  <dbl> 1.69947, 1.60053, 1.92989, 1.49806, 1.58634, 1.98859,…
-## $ `s(year_fac).7`  <dbl> 2.05134, 2.03602, 2.13280, 2.11045, 2.26266, 1.83901,…
-## $ `s(year_fac).8`  <dbl> 2.99787, 2.97008, 2.99874, 2.91971, 2.91366, 3.02957,…
-## $ `s(year_fac).9`  <dbl> 3.33622, 3.20377, 3.32396, 3.24947, 3.20354, 3.23998,…
-## $ `s(year_fac).10` <dbl> 2.82800, 2.76225, 2.76828, 2.72482, 2.75613, 2.79170,…
-## $ `s(year_fac).11` <dbl> 3.10896, 3.01869, 3.16489, 2.99758, 3.12221, 3.17723,…
-## $ `s(year_fac).12` <dbl> 3.21659, 3.04013, 3.37160, 3.30859, 3.26822, 3.30016,…
-## $ `s(year_fac).13` <dbl> 2.36270, 2.33583, 2.26894, 2.07640, 2.18187, 2.36848,…
-## $ `s(year_fac).14` <dbl> 2.69838, 2.64463, 2.66172, 2.59959, 2.59452, 2.64214,…
-## $ `s(year_fac).15` <dbl> 2.25933, 2.25514, 2.07515, 2.15506, 2.17237, 2.20526,…
-## $ `s(year_fac).16` <dbl> 2.09584, 2.10991, 1.83435, 2.20681, 2.26212, 1.92478,…
-## $ `s(year_fac).17` <dbl> 0.4150130, 0.4578770, 0.2064460, 0.9136730, 1.0479400…
+## $ `s(year_fac).1` <dbl> 1.86190, 2.11388, 2.11600, 2.10553, 2.23891, 1.86875,… +## $ `s(year_fac).2` <dbl> 2.64166, 2.68429, 2.72846, 2.64397, 2.77831, 2.60079,… +## $ `s(year_fac).3` <dbl> 3.16581, 3.08324, 3.08527, 2.97360, 3.14883, 3.12684,… +## $ `s(year_fac).4` <dbl> 3.27785, 3.18822, 3.23306, 3.22290, 3.19447, 3.38818,… +## $ `s(year_fac).5` <dbl> 2.09716, 2.03221, 2.16466, 2.03947, 2.17555, 2.10229,… +## $ `s(year_fac).6` <dbl> 1.71619, 1.58724, 1.86737, 1.90242, 1.56511, 1.75312,… +## $ `s(year_fac).7` <dbl> 1.96664, 2.09183, 2.15583, 2.10924, 2.00211, 2.11828,… +## $ `s(year_fac).8` <dbl> 2.89847, 2.80166, 2.97784, 2.90876, 3.02493, 2.90480,… +## $ `s(year_fac).9` <dbl> 3.24516, 3.23472, 3.17024, 3.28808, 3.10487, 3.24116,… +## $ `s(year_fac).10` <dbl> 2.70289, 2.75808, 2.63608, 2.81192, 2.77074, 2.76922,… +## $ `s(year_fac).11` <dbl> 3.03255, 3.04433, 3.03362, 3.01215, 3.01482, 3.19735,… +## $ `s(year_fac).12` <dbl> 3.21954, 3.16691, 3.18440, 3.17227, 3.14342, 3.22443,… +## $ `s(year_fac).13` <dbl> 2.18136, 2.38175, 2.18528, 2.29116, 2.27636, 2.24354,… +## $ `s(year_fac).14` <dbl> 2.68183, 2.68278, 2.65206, 2.68147, 2.67352, 2.62832,… +## $ `s(year_fac).15` <dbl> 1.94873, 2.16423, 2.07559, 2.12836, 2.03424, 2.32574,… +## $ `s(year_fac).16` <dbl> 2.16423, 2.37222, 1.92845, 2.11019, 2.03066, 2.17096,… +## $ `s(year_fac).17` <dbl> 0.9960880, 0.9922320, 0.9462070, 0.6951830, 0.8883960…

With any model fitted in mvgam, the underlying Stan code can be viewed using the code function:

@@ -877,7 +877,7 @@

## $ test_observations : NULL ## $ test_times : NULL ## $ hindcasts :List of 1 -## ..$ PP: num [1:2000, 1:199] 6 9 9 5 4 8 8 5 6 9 ... +## ..$ PP: num [1:2000, 1:199] 7 6 7 5 9 7 4 10 9 8 ... ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : NULL ## .. .. ..$ : chr [1:199] "ypred[1,1]" "ypred[2,1]" "ypred[3,1]" "ypred[4,1]" ... @@ -890,7 +890,7 @@ 

 hc <- hindcast(model1, type = 'link')
 range(hc$hindcasts$PP)
-
## [1] -1.65123  3.47762
+
## [1] -1.79937  3.44340

Objects of class mvgam_forecast have an associated plot function as well:

@@ -942,14 +942,14 @@ 
Automatic forecasting for new dataplot(model1b, type = 'forecast')

## Out of sample DRPS:
-## [1] 182.1551
+## [1] 184.8855

We can also view the test data in the forecast plot to see that the predictions do not capture the temporal variation in the test set

 plot(model1b, type = 'forecast', newdata = data_test)

## Out of sample DRPS:
-## [1] 182.1551
+## [1] 184.8855

As with the hindcast function, we can use the forecast function to automatically extract the posterior distributions for these predictions. This also returns an object of @@ -977,12 +977,12 @@

Automatic forecasting for new data## ..$ PP: int [1:39] NA 0 0 10 3 14 18 NA 28 46 ... ## $ test_times : num [1:39] 161 162 163 164 165 166 167 168 169 170 ... ## $ hindcasts :List of 1 -## ..$ PP: num [1:2000, 1:160] 11 10 4 6 7 3 12 5 4 8 ... +## ..$ PP: num [1:2000, 1:160] 11 9 8 6 9 7 9 7 10 9 ... ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : NULL ## .. .. ..$ : chr [1:160] "ypred[1,1]" "ypred[2,1]" "ypred[3,1]" "ypred[4,1]" ... ## $ forecasts :List of 1 -## ..$ PP: num [1:2000, 1:39] 9 9 16 3 9 8 11 20 9 10 ... +## ..$ PP: num [1:2000, 1:39] 12 7 11 10 5 12 9 10 6 6 ... ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : NULL ## .. .. ..$ : chr [1:39] "ypred[161,1]" "ypred[162,1]" "ypred[163,1]" "ypred[164,1]" ... @@ -1039,33 +1039,33 @@

Adding predictors as “fixed” eff ## ## GAM coefficient (beta) estimates: ## 2.5% 50% 97.5% Rhat n_eff -## ndvi 0.32 0.39 0.46 1 1786 -## s(year_fac).1 1.10 1.40 1.70 1 2685 -## s(year_fac).2 1.80 2.00 2.20 1 2607 -## s(year_fac).3 2.20 2.40 2.60 1 2057 -## s(year_fac).4 2.30 2.50 2.70 1 2129 -## s(year_fac).5 1.20 1.40 1.60 1 2567 +## ndvi 0.32 0.39 0.46 1 1720 +## s(year_fac).1 1.10 1.40 1.70 1 2505 +## s(year_fac).2 1.80 2.00 2.20 1 2281 +## s(year_fac).3 2.20 2.40 2.60 1 1897 +## s(year_fac).4 2.30 2.50 2.70 1 2232 +## s(year_fac).5 1.20 1.40 1.60 1 2288 ## s(year_fac).6 1.00 1.30 1.50 1 2267 -## s(year_fac).7 1.10 1.40 1.70 1 2636 -## s(year_fac).8 2.10 2.30 2.50 1 2082 -## s(year_fac).9 2.70 2.90 3.00 1 2173 -## s(year_fac).10 2.00 2.20 2.40 1 2669 -## s(year_fac).11 2.30 2.40 2.60 1 2050 -## s(year_fac).12 2.60 2.70 2.80 1 2114 -## s(year_fac).13 1.40 1.60 1.80 1 2484 -## s(year_fac).14 0.68 2.00 3.30 1 1519 -## s(year_fac).15 0.76 2.00 3.40 1 1648 -## s(year_fac).16 0.64 2.00 3.20 1 1658 -## s(year_fac).17 0.61 2.00 3.30 1 1370 +## s(year_fac).7 1.10 1.40 1.70 1 2556 +## s(year_fac).8 2.10 2.30 2.50 1 2217 +## s(year_fac).9 2.70 2.90 3.00 1 1968 +## s(year_fac).10 2.00 2.20 2.40 1 2603 +## s(year_fac).11 2.30 2.40 2.60 1 2064 +## s(year_fac).12 2.50 2.70 2.80 1 1902 +## s(year_fac).13 1.40 1.60 1.80 1 2325 +## s(year_fac).14 0.63 2.00 3.30 1 1454 +## s(year_fac).15 0.73 2.00 3.20 1 1673 +## s(year_fac).16 0.74 1.90 3.20 1 1613 +## s(year_fac).17 0.73 2.00 3.20 1 1416 ## ## GAM group-level estimates: ## 2.5% 50% 97.5% Rhat n_eff -## mean(s(year_fac)) 1.6 2.00 2.30 1.00 456 -## sd(s(year_fac)) 0.4 0.59 0.99 1.02 352 +## mean(s(year_fac)) 1.6 2.00 2.30 1.02 451 +## sd(s(year_fac)) 0.4 0.59 0.96 1.00 546 ## ## Approximate significance of GAM observation smooths: ## edf Chi.sq p-value -## s(year_fac) 10.3 2802 <2e-16 *** +## s(year_fac) 10.9 2832 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## @@ -1076,7 +1076,7 @@

Adding predictors as “fixed” eff ## 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 Mon Jan 15 9:07:58 AM 2024. +## Samples were drawn using NUTS(diag_e) at Fri Jan 19 11:55:20 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) @@ -1087,24 +1087,24 @@

Adding predictors as “fixed” eff 
 coef(model2)
##                     2.5%       50%     97.5% Rhat n_eff
-## ndvi           0.3214333 0.3900135 0.4563159    1  1786
-## s(year_fac).1  1.1140800 1.4079600 1.6670795    1  2685
-## s(year_fac).2  1.7910845 2.0025550 2.1936545    1  2607
-## s(year_fac).3  2.1944393 2.3821600 2.5643877    1  2057
-## s(year_fac).4  2.3301122 2.5071400 2.6851682    1  2129
-## s(year_fac).5  1.1836370 1.4268450 1.6446277    1  2567
-## s(year_fac).6  1.0164710 1.2687300 1.5169110    1  2267
-## s(year_fac).7  1.1292298 1.4117000 1.6742532    1  2636
-## s(year_fac).8  2.0858535 2.2694550 2.4555455    1  2082
-## s(year_fac).9  2.7214870 2.8538250 2.9875630    1  2173
-## s(year_fac).10 1.9630840 2.1858000 2.3866162    1  2669
-## s(year_fac).11 2.2738888 2.4383000 2.6029485    1  2050
-## s(year_fac).12 2.5519082 2.6920100 2.8434890    1  2114
-## s(year_fac).13 1.4001263 1.6197550 1.8443220    1  2484
-## s(year_fac).14 0.6821051 1.9962300 3.3376105    1  1519
-## s(year_fac).15 0.7618126 2.0004150 3.3627562    1  1648
-## s(year_fac).16 0.6390487 1.9976450 3.2277897    1  1658
-## s(year_fac).17 0.6054341 1.9780000 3.3044415    1  1370
+## ndvi 0.3238703 0.3902365 0.4581023 1 1720 +## s(year_fac).1 1.1365188 1.4058450 1.6577235 1 2505 +## s(year_fac).2 1.8010670 2.0009350 2.2021112 1 2281 +## s(year_fac).3 2.1889132 2.3762450 2.5703225 1 1897 +## s(year_fac).4 2.3193515 2.5045900 2.6835790 1 2232 +## s(year_fac).5 1.1944500 1.4269400 1.6385940 1 2288 +## s(year_fac).6 1.0194370 1.2684500 1.5165563 1 2267 +## s(year_fac).7 1.1363357 1.4122750 1.6805735 1 2556 +## s(year_fac).8 2.0830938 2.2729100 2.4619788 1 2217 +## s(year_fac).9 2.7169883 2.8535000 2.9868788 1 1968 +## s(year_fac).10 1.9698398 2.1826100 2.3704190 1 2603 +## s(year_fac).11 2.2599978 2.4376100 2.5998005 1 2064 +## s(year_fac).12 2.5364615 2.6917000 2.8401345 1 1902 +## s(year_fac).13 1.3762790 1.6143250 1.8473850 1 2325 +## s(year_fac).14 0.6335216 2.0001150 3.3246795 1 1454 +## s(year_fac).15 0.7257230 1.9997600 3.1933355 1 1673 +## s(year_fac).16 0.7395432 1.9439600 3.2223760 1 1613 +## s(year_fac).17 0.7292550 1.9868700 3.2186612 1 1416

Look at the estimated effect of ndvi using plot.mvgam with type = 'pterms'

@@ -1120,24 +1120,24 @@ 
Adding predictors as “fixed” eff
 dplyr::glimpse(beta_post)
## Rows: 2,000
 ## Columns: 18
-## $ ndvi             <dbl> 0.365016, 0.383852, 0.382832, 0.374714, 0.385262, 0.3…
-## $ `s(year_fac).1`  <dbl> 1.26080, 1.58275, 1.32614, 1.34134, 1.36449, 1.43590,…
-## $ `s(year_fac).2`  <dbl> 2.11456, 1.97718, 1.83386, 1.89789, 1.95921, 2.05343,…
-## $ `s(year_fac).3`  <dbl> 2.30024, 2.37424, 2.42140, 2.45382, 2.38137, 2.43344,…
-## $ `s(year_fac).4`  <dbl> 2.58320, 2.48406, 2.51264, 2.51145, 2.41215, 2.49424,…
-## $ `s(year_fac).5`  <dbl> 1.65731, 1.21655, 1.27151, 1.29860, 1.40799, 1.44363,…
-## $ `s(year_fac).6`  <dbl> 1.11182, 1.40885, 1.23147, 1.32212, 1.30500, 1.24335,…
-## $ `s(year_fac).7`  <dbl> 1.55856, 1.20594, 1.60187, 1.65128, 1.18561, 1.65570,…
-## $ `s(year_fac).8`  <dbl> 2.31500, 2.17933, 2.19752, 2.27384, 2.21278, 2.26591,…
-## $ `s(year_fac).9`  <dbl> 2.86616, 2.85139, 2.92218, 2.84696, 2.80555, 2.81140,…
-## $ `s(year_fac).10` <dbl> 2.25933, 2.12217, 2.14634, 2.24707, 2.19890, 2.17169,…
-## $ `s(year_fac).11` <dbl> 2.52557, 2.49356, 2.62267, 2.52679, 2.43452, 2.40733,…
-## $ `s(year_fac).12` <dbl> 2.72165, 2.70566, 2.67655, 2.70252, 2.68610, 2.71396,…
-## $ `s(year_fac).13` <dbl> 1.55743, 1.70503, 1.48361, 1.55631, 1.58906, 1.62060,…
-## $ `s(year_fac).14` <dbl> 2.566090, 1.447150, 1.984190, 2.192650, 1.231120, 1.9…
-## $ `s(year_fac).15` <dbl> 1.05225, 2.12175, 1.72275, 1.93314, 1.75620, 1.66417,…
-## $ `s(year_fac).16` <dbl> 1.447540, 2.158620, 1.625380, 1.515920, 2.205820, 0.7…
-## $ `s(year_fac).17` <dbl> 0.828068, 0.831492, 2.645610, 2.822590, 0.700252, 2.3…
+## $ ndvi <dbl> 0.425996, 0.413570, 0.419061, 0.422454, 0.398450, 0.3… +## $ `s(year_fac).1` <dbl> 1.57912, 1.35922, 1.35987, 1.26014, 1.35486, 1.34965,… +## $ `s(year_fac).2` <dbl> 1.99289, 1.93975, 1.90113, 2.01752, 1.89284, 2.00161,… +## $ `s(year_fac).3` <dbl> 2.39844, 2.43849, 2.17660, 2.43780, 2.47060, 2.40514,… +## $ `s(year_fac).4` <dbl> 2.42520, 2.43230, 2.38353, 2.50616, 2.49209, 2.51539,… +## $ `s(year_fac).5` <dbl> 1.41111, 1.50256, 1.20099, 1.48075, 1.47780, 1.43706,… +## $ `s(year_fac).6` <dbl> 1.18568, 1.01049, 1.22573, 1.09919, 1.05730, 1.33812,… +## $ `s(year_fac).7` <dbl> 1.38444, 1.50095, 1.34319, 1.39171, 1.28743, 1.50470,… +## $ `s(year_fac).8` <dbl> 2.24025, 2.19628, 2.22380, 2.18654, 2.29103, 2.33180,… +## $ `s(year_fac).9` <dbl> 2.85087, 2.82355, 2.77289, 2.91273, 2.70009, 2.73278,… +## $ `s(year_fac).10` <dbl> 2.26549, 2.02703, 2.24575, 2.09158, 2.15786, 2.19188,… +## $ `s(year_fac).11` <dbl> 2.27807, 2.21506, 2.36005, 2.40034, 2.35533, 2.38627,… +## $ `s(year_fac).12` <dbl> 2.63560, 2.58363, 2.65288, 2.62733, 2.62621, 2.62781,… +## $ `s(year_fac).13` <dbl> 1.67680, 1.55301, 1.66027, 1.45560, 1.65408, 1.56315,… +## $ `s(year_fac).14` <dbl> 2.898330, 4.072930, 3.201270, 3.381070, 1.957030, 2.3… +## $ `s(year_fac).15` <dbl> 3.015740, 2.935880, 3.108230, 3.318080, 1.149890, 2.9… +## $ `s(year_fac).16` <dbl> 3.347350, 2.120660, 1.605260, 1.202070, 1.962580, 2.2… +## $ `s(year_fac).17` <dbl> 2.893970, 2.884200, 2.368160, 1.820480, 1.562180, 2.3…

The posterior distribution for the effect of ndvi is stored in the ndvi column. A quick histogram confirms our inference that log(counts) respond positively to increases @@ -1274,26 +1274,26 @@

Adding predictors as smooths## ## GAM coefficient (beta) estimates: ## 2.5% 50% 97.5% Rhat n_eff -## (Intercept) 2.00 2.10 2.200 1.00 1043 -## ndvi 0.26 0.33 0.390 1.00 1059 -## s(time).1 -2.10 -1.10 -0.031 1.01 337 -## s(time).2 0.44 1.30 2.300 1.02 283 -## s(time).3 -0.55 0.44 1.500 1.02 243 -## s(time).4 1.60 2.50 3.600 1.02 262 -## s(time).5 -1.20 -0.21 0.910 1.02 237 -## s(time).6 -0.62 0.37 1.600 1.02 281 -## s(time).7 -1.50 -0.54 0.590 1.02 249 -## s(time).8 0.60 1.40 2.600 1.02 248 -## s(time).9 1.20 2.10 3.200 1.02 225 -## s(time).10 -0.39 0.54 1.700 1.02 259 -## s(time).11 0.83 1.70 2.900 1.02 240 -## s(time).12 0.67 1.50 2.500 1.02 254 -## s(time).13 -1.20 -0.32 0.670 1.01 322 -## s(time).14 -7.60 -4.10 -1.200 1.01 351 +## (Intercept) 2.00 2.10 2.200 1.01 1231 +## ndvi 0.26 0.33 0.390 1.01 1215 +## s(time).1 -2.10 -1.10 0.078 1.01 517 +## s(time).2 0.50 1.30 2.400 1.00 395 +## s(time).3 -0.46 0.46 1.500 1.01 359 +## s(time).4 1.60 2.50 3.600 1.01 340 +## s(time).5 -1.10 -0.20 0.840 1.01 359 +## s(time).6 -0.54 0.39 1.600 1.00 380 +## s(time).7 -1.50 -0.51 0.550 1.01 381 +## s(time).8 0.60 1.50 2.600 1.01 350 +## s(time).9 1.20 2.10 3.200 1.01 325 +## s(time).10 -0.32 0.55 1.600 1.01 355 +## s(time).11 0.88 1.80 2.900 1.01 324 +## s(time).12 0.71 1.50 2.500 1.01 368 +## s(time).13 -1.10 -0.29 0.650 1.01 428 +## s(time).14 -7.60 -4.30 -1.200 1.00 417 ## ## Approximate significance of GAM observation smooths: ## edf Chi.sq p-value -## s(time) 8.57 797 <2e-16 *** +## s(time) 9.44 793 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## @@ -1304,7 +1304,7 @@

Adding predictors as smooths## 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 Mon Jan 15 9:08:44 AM 2024. +## Samples were drawn using NUTS(diag_e) at Fri Jan 19 11:56:04 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) @@ -1431,7 +1431,7 @@

Latent dynamics in mvgamplot(model3, type = 'forecast', newdata = data_test) 

## Out of sample DRPS:
-## [1] 286.5788
+## [1] 288.607

Why is this happening? The forecasts are driven almost entirely by variation in the temporal spline, which is extrapolating linearly forever beyond the edge of the training data. Any slight @@ -1514,32 +1514,33 @@

Latent dynamics in mvgam## ## GAM coefficient (beta) estimates: ## 2.5% 50% 97.5% Rhat n_eff -## (Intercept) 0.890 1.9000 2.500 1.03 61 -## s(ndvi).1 -0.082 0.0110 0.220 1.00 533 -## s(ndvi).2 -0.160 0.0140 0.400 1.00 454 -## s(ndvi).3 -0.054 -0.0018 0.061 1.01 760 -## s(ndvi).4 -0.300 0.1200 1.700 1.01 268 -## s(ndvi).5 -0.100 0.1500 0.340 1.00 462 +## (Intercept) 0.640 1.9000 2.400 1.08 72 +## s(ndvi).1 -0.083 0.0110 0.170 1.00 599 +## s(ndvi).2 -0.160 0.0180 0.320 1.00 465 +## s(ndvi).3 -0.062 -0.0012 0.052 1.00 755 +## s(ndvi).4 -0.280 0.1300 1.300 1.01 266 +## s(ndvi).5 -0.073 0.1500 0.340 1.00 643 ## ## Approximate significance of GAM observation smooths: -## edf Chi.sq p-value -## s(ndvi) 0.978 95.5 0.088 . +## edf Chi.sq p-value +## s(ndvi) 1 84.1 0.078 . ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Latent trend parameter AR estimates: ## 2.5% 50% 97.5% Rhat n_eff -## ar1[1] 0.69 0.81 0.94 1.01 173 -## sigma[1] 0.67 0.80 0.96 1.01 518 +## ar1[1] 0.70 0.81 0.94 1.03 164 +## sigma[1] 0.67 0.79 0.95 1.01 498 ## ## Stan MCMC diagnostics: ## n_eff / iter looks reasonable for all parameters -## Rhat looks reasonable for all parameters +## Rhats above 1.05 found for 79 parameters +## *Diagnose further to investigate why the chains have not mixed ## 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 Mon Jan 15 9:10:06 AM 2024. +## Samples were drawn using NUTS(diag_e) at Fri Jan 19 11:57:03 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) @@ -1555,7 +1556,7 @@

Latent dynamics in mvgamplot(model4, type = 'forecast', newdata = data_test) 

## Out of sample DRPS:
-## [1] 150.9452
+## [1] 150.1069

The trend is evolving as an AR1 process, which we can also view:

 plot(model4, type = 'trend', newdata = data_test)
@@ -1570,7 +1571,7 @@

Latent dynamics in mvgam## Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details. 
##        elpd_diff se_diff
 ## model4    0.0       0.0 
-## model3 -560.2      66.5
+## model3 -557.6 66.2

The higher estimated log predictive density (ELPD) value for the dynamic model suggests it provides a better fit to the in-sample data.

@@ -1585,7 +1586,7 @@

Latent dynamics in mvgamscore_mod3 <- score(fc_mod3, score = 'drps') score_mod4 <- score(fc_mod4, score = 'drps') sum(score_mod4$PP$score, na.rm = TRUE) - sum(score_mod3$PP$score, na.rm = TRUE) -
## [1] -135.6337
+
## [1] -138.5

A strongly negative value here suggests the score for the dynamic model (model 4) is much smaller than the score for the model with a smooth function of time (model 3)

diff --git a/docs/articles/mvgam_overview_files/figure-html/Histogram of NDVI effects-1.png b/docs/articles/mvgam_overview_files/figure-html/Histogram of NDVI effects-1.png index 56b72cb6..f9e3d96c 100644 Binary files a/docs/articles/mvgam_overview_files/figure-html/Histogram of NDVI effects-1.png and b/docs/articles/mvgam_overview_files/figure-html/Histogram of NDVI effects-1.png differ diff --git a/docs/articles/mvgam_overview_files/figure-html/Plot NDVI effect-1.png b/docs/articles/mvgam_overview_files/figure-html/Plot NDVI effect-1.png index bcdffa60..ce70d3e6 100644 Binary files a/docs/articles/mvgam_overview_files/figure-html/Plot NDVI effect-1.png and b/docs/articles/mvgam_overview_files/figure-html/Plot NDVI effect-1.png differ diff --git a/docs/articles/mvgam_overview_files/figure-html/Plot extrapolated temporal functions using newdata-1.png b/docs/articles/mvgam_overview_files/figure-html/Plot extrapolated temporal functions using newdata-1.png index 1109c200..d39cc498 100644 Binary files a/docs/articles/mvgam_overview_files/figure-html/Plot extrapolated temporal functions using newdata-1.png and b/docs/articles/mvgam_overview_files/figure-html/Plot extrapolated temporal functions using newdata-1.png differ diff --git a/docs/articles/mvgam_overview_files/figure-html/Plot hindcasts on the linear predictor scale-1.png b/docs/articles/mvgam_overview_files/figure-html/Plot hindcasts on the linear predictor scale-1.png index a001b44e..e1ba65a7 100644 Binary files a/docs/articles/mvgam_overview_files/figure-html/Plot hindcasts on the linear predictor scale-1.png and b/docs/articles/mvgam_overview_files/figure-html/Plot hindcasts on the linear predictor scale-1.png differ diff --git a/docs/articles/mvgam_overview_files/figure-html/Plot posterior hindcasts-1.png b/docs/articles/mvgam_overview_files/figure-html/Plot posterior hindcasts-1.png index a0d871f8..20a70abb 100644 Binary files a/docs/articles/mvgam_overview_files/figure-html/Plot posterior hindcasts-1.png and b/docs/articles/mvgam_overview_files/figure-html/Plot posterior hindcasts-1.png differ diff --git a/docs/articles/mvgam_overview_files/figure-html/Plot posterior residuals-1.png b/docs/articles/mvgam_overview_files/figure-html/Plot posterior residuals-1.png index 08683df5..82e9b96a 100644 Binary files a/docs/articles/mvgam_overview_files/figure-html/Plot posterior residuals-1.png and b/docs/articles/mvgam_overview_files/figure-html/Plot posterior residuals-1.png differ diff --git a/docs/articles/mvgam_overview_files/figure-html/Plot random effect estimates-1.png b/docs/articles/mvgam_overview_files/figure-html/Plot random effect estimates-1.png index d794301d..50698111 100644 Binary files a/docs/articles/mvgam_overview_files/figure-html/Plot random effect estimates-1.png and b/docs/articles/mvgam_overview_files/figure-html/Plot random effect estimates-1.png differ diff --git a/docs/articles/mvgam_overview_files/figure-html/Plot smooth term derivatives-1.png b/docs/articles/mvgam_overview_files/figure-html/Plot smooth term derivatives-1.png index 7c90e2d2..0f46b68f 100644 Binary files a/docs/articles/mvgam_overview_files/figure-html/Plot smooth term derivatives-1.png and b/docs/articles/mvgam_overview_files/figure-html/Plot smooth term derivatives-1.png differ diff --git a/docs/articles/mvgam_overview_files/figure-html/Plotting predictions against test data-1.png b/docs/articles/mvgam_overview_files/figure-html/Plotting predictions against test data-1.png index 3a4bedc4..0bf4b031 100644 Binary files a/docs/articles/mvgam_overview_files/figure-html/Plotting predictions against test data-1.png and b/docs/articles/mvgam_overview_files/figure-html/Plotting predictions against test data-1.png differ diff --git a/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-13-1.png b/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-13-1.png index d23650fc..f26b8f48 100644 Binary files a/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-13-1.png and b/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-13-1.png differ diff --git a/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-17-1.png b/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-17-1.png index 90fc8860..7c773fe9 100644 Binary files a/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-17-1.png and b/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-17-1.png differ diff --git a/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-18-1.png b/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-18-1.png index 3a4bedc4..0bf4b031 100644 Binary files a/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-18-1.png and b/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-18-1.png differ diff --git a/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-22-1.png b/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-22-1.png index b8408084..8d543b8f 100644 Binary files a/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-22-1.png and b/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-22-1.png differ diff --git a/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-23-1.png b/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-23-1.png index 52974248..c297f088 100644 Binary files a/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-23-1.png and b/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-23-1.png differ diff --git a/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-23-2.png b/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-23-2.png index ae02aec2..935a8cf2 100644 Binary files a/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-23-2.png and b/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-23-2.png differ diff --git a/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-26-1.png b/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-26-1.png index 036dd03e..75cc413c 100644 Binary files a/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-26-1.png and b/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-26-1.png differ diff --git a/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-27-1.png b/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-27-1.png index a244e649..cee5fa21 100644 Binary files a/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-27-1.png and b/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-27-1.png differ diff --git a/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-28-1.png b/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-28-1.png index 66a4b17b..29efc6e3 100644 Binary files a/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-28-1.png and b/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-28-1.png differ diff --git a/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-28-2.png b/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-28-2.png index 84cd11b3..a4995c7f 100644 Binary files a/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-28-2.png and b/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-28-2.png differ diff --git a/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-29-1.png b/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-29-1.png index b289cefe..08a531b8 100644 Binary files a/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-29-1.png and b/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-29-1.png differ diff --git a/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-29-2.png b/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-29-2.png index f5219b5c..3dd8f4e4 100644 Binary files a/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-29-2.png and b/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-29-2.png differ diff --git a/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-31-1.png b/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-31-1.png index 275429b3..10725ddf 100644 Binary files a/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-31-1.png and b/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-31-1.png differ diff --git a/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-33-1.png b/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-33-1.png index 29cbf193..9c1ce98e 100644 Binary files a/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-33-1.png and b/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-33-1.png differ diff --git a/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-34-1.png b/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-34-1.png index 1cb85882..65c0fae1 100644 Binary files a/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-34-1.png and b/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-34-1.png differ diff --git a/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-35-1.png b/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-35-1.png index c03916e2..c02d880c 100644 Binary files a/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-35-1.png and b/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-35-1.png differ diff --git a/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-6-1.png b/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-6-1.png index a8a340cd..eed90f62 100644 Binary files a/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-6-1.png and b/docs/articles/mvgam_overview_files/figure-html/unnamed-chunk-6-1.png differ diff --git a/docs/articles/shared_states.html b/docs/articles/shared_states.html index c8b68c3d..d78c7cca 100644 --- a/docs/articles/shared_states.html +++ b/docs/articles/shared_states.html @@ -77,7 +77,7 @@

Nicholas J Clark

-

2024-01-12

+

2024-01-19

Source: vignettes/shared_states.Rmd
shared_states.Rmd
@@ -185,7 +185,7 @@

Checking trend_map ## matrix[total_obs, num_basis] X; // mgcv GAM design matrix ## matrix[n * n_lv, num_basis_trend] X_trend; // trend model design matrix ## array[n, n_series] int<lower=0> ytimes; // time-ordered matrix (which col in X belongs to each [time, series] observation?) -## array[n, n_lv] int<lower=0> ytimes_trend; // time-ordered matrix for latent states +## array[n, n_lv] int ytimes_trend; ## int<lower=0> n_nonmissing; // number of nonmissing observations ## matrix[4, 4] S_trend1; // mgcv smooth penalty matrix S_trend1 ## array[n_nonmissing] int<lower=0> flat_ys; // flattened nonmissing observations @@ -357,30 +357,30 @@

Fitting and inspecting the model## ## GAM observation model coefficient (beta) estimates: ## 2.5% 50% 97.5% Rhat n_eff -## seriesseries_1 -0.17 0.091 0.32 1.00 1357 -## seriesseries_2 0.92 1.100 1.20 1.00 1024 -## seriesseries_3 1.90 2.100 2.30 1.02 275 +## seriesseries_1 -0.15 0.088 0.31 1.00 1895 +## seriesseries_2 0.92 1.100 1.20 1.00 1267 +## seriesseries_3 1.90 2.100 2.30 1.02 256 ## ## Process model AR parameter estimates: -## 2.5% 50% 97.5% Rhat n_eff -## ar1[1] -0.73 -0.4300 -0.065 1 605 -## ar1[2] -0.30 -0.0047 0.270 1 1457 +## 2.5% 50% 97.5% Rhat n_eff +## ar1[1] -0.72 -0.420 -0.056 1 676 +## ar1[2] -0.28 -0.011 0.280 1 1433 ## ## Process error parameter estimates: ## 2.5% 50% 97.5% Rhat n_eff -## sigma[1] 0.32 0.49 0.67 1.01 296 -## sigma[2] 0.60 0.73 0.91 1.00 1144 +## sigma[1] 0.33 0.49 0.67 1 487 +## sigma[2] 0.59 0.73 0.91 1 948 ## ## GAM process model coefficient (beta) estimates: -## 2.5% 50% 97.5% Rhat n_eff -## s(season).1_trend -0.22 -0.0093 0.22 1 1636 -## s(season).2_trend -0.30 -0.0500 0.16 1 1427 -## s(season).3_trend -0.15 0.0730 0.29 1 1408 -## s(season).4_trend -0.14 0.0650 0.27 1 1654 +## 2.5% 50% 97.5% Rhat n_eff +## s(season).1_trend -0.22 -0.011 0.20 1 1612 +## s(season).2_trend -0.27 -0.045 0.18 1 1745 +## s(season).3_trend -0.15 0.074 0.29 1 1347 +## s(season).4_trend -0.15 0.067 0.28 1 1561 ## ## Approximate significance of GAM process smooths: ## edf F p-value -## s(season) 1.91 0.1 0.91 +## s(season) 1.52 0.1 0.91 ## ## Stan MCMC diagnostics: ## n_eff / iter looks reasonable for all parameters @@ -389,7 +389,7 @@

Fitting and inspecting the model## 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 Fri Jan 12 3:56:28 PM 2024. +## Samples were drawn using NUTS(diag_e) at Fri Jan 19 11:57:59 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) @@ -598,45 +598,46 @@

The shared signal model## ## Observation error parameter estimates: ## 2.5% 50% 97.5% Rhat n_eff -## sigma_obs[1] 0.57 1.2 1.9 1.02 122 -## sigma_obs[2] 1.80 2.1 2.4 1.00 2326 -## sigma_obs[3] 1.80 2.1 2.3 1.00 2364 +## sigma_obs[1] 1.6 1.9 2.2 1 1757 +## sigma_obs[2] 1.4 1.7 2.0 1 1090 +## sigma_obs[3] 1.3 1.5 1.8 1 1339 ## ## GAM observation model coefficient (beta) estimates: ## 2.5% 50% 97.5% Rhat n_eff -## (Intercept) 0.76 1.70 2.500 1.01 702 -## seriessensor_2 -2.10 -0.94 0.300 1.00 940 -## seriessensor_3 -3.50 -2.00 0.016 1.01 849 +## (Intercept) 0.72 1.70 2.50 1.01 360 +## seriessensor_2 -2.10 -0.96 0.32 1.00 1068 +## seriessensor_3 -3.40 -2.00 -0.39 1.00 1154 ## ## Approximate significance of GAM observation smooths: -## edf F p-value -## s(temperature) 7.25 8.25 <2e-16 *** -## s(series,temperature) 2.35 0.95 0.001 ** +## edf F p-value +## s(temperature) 1.22 12.66 < 2e-16 *** +## s(series,temperature) 1.92 0.95 6.9e-06 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Process model AR parameter estimates: -## 2.5% 50% 97.5% Rhat n_eff -## ar1[1] -0.19 0.12 0.44 1 878 +## 2.5% 50% 97.5% Rhat n_eff +## ar1[1] 0.33 0.59 0.83 1 492 ## ## Process error parameter estimates: ## 2.5% 50% 97.5% Rhat n_eff -## sigma[1] 1.1 1.8 2.3 1.04 133 +## sigma[1] 0.72 1 1.3 1.01 392 ## ## Approximate significance of GAM process smooths: -## edf F p-value -## s(productivity) 1.47 0.31 0.96 +## edf F p-value +## s(productivity) 3.6 8.31 5.1e-05 *** +## --- +## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Stan MCMC diagnostics: ## n_eff / iter looks reasonable for all parameters -## Rhats above 1.05 found for 49 parameters +## Rhats above 1.05 found for 28 parameters ## *Diagnose further to investigate why the chains have not mixed -## 3 of 2000 iterations ended with a divergence (0.15%) -## *Try running with larger adapt_delta to remove the divergences +## 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 Fri Jan 12 3:58:16 PM 2024. +## Samples were drawn using NUTS(diag_e) at Fri Jan 19 11:59:46 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) diff --git a/docs/articles/shared_states_files/figure-html/unnamed-chunk-10-1.png b/docs/articles/shared_states_files/figure-html/unnamed-chunk-10-1.png index a7226bac..f154a81b 100644 Binary files a/docs/articles/shared_states_files/figure-html/unnamed-chunk-10-1.png and b/docs/articles/shared_states_files/figure-html/unnamed-chunk-10-1.png differ diff --git a/docs/articles/shared_states_files/figure-html/unnamed-chunk-10-2.png b/docs/articles/shared_states_files/figure-html/unnamed-chunk-10-2.png index 20c81dec..f536f4fc 100644 Binary files a/docs/articles/shared_states_files/figure-html/unnamed-chunk-10-2.png and b/docs/articles/shared_states_files/figure-html/unnamed-chunk-10-2.png differ diff --git a/docs/articles/shared_states_files/figure-html/unnamed-chunk-10-3.png b/docs/articles/shared_states_files/figure-html/unnamed-chunk-10-3.png index f744c96c..3d449058 100644 Binary files a/docs/articles/shared_states_files/figure-html/unnamed-chunk-10-3.png and b/docs/articles/shared_states_files/figure-html/unnamed-chunk-10-3.png differ diff --git a/docs/articles/shared_states_files/figure-html/unnamed-chunk-19-1.png b/docs/articles/shared_states_files/figure-html/unnamed-chunk-19-1.png index ca0cb1e2..70bc52fe 100644 Binary files a/docs/articles/shared_states_files/figure-html/unnamed-chunk-19-1.png and b/docs/articles/shared_states_files/figure-html/unnamed-chunk-19-1.png differ diff --git a/docs/articles/shared_states_files/figure-html/unnamed-chunk-20-1.png b/docs/articles/shared_states_files/figure-html/unnamed-chunk-20-1.png index 282edf16..a94c305e 100644 Binary files a/docs/articles/shared_states_files/figure-html/unnamed-chunk-20-1.png and b/docs/articles/shared_states_files/figure-html/unnamed-chunk-20-1.png differ diff --git a/docs/articles/shared_states_files/figure-html/unnamed-chunk-21-1.png b/docs/articles/shared_states_files/figure-html/unnamed-chunk-21-1.png index c3c9fb90..bfbe1bb0 100644 Binary files a/docs/articles/shared_states_files/figure-html/unnamed-chunk-21-1.png and b/docs/articles/shared_states_files/figure-html/unnamed-chunk-21-1.png differ diff --git a/docs/articles/shared_states_files/figure-html/unnamed-chunk-21-2.png b/docs/articles/shared_states_files/figure-html/unnamed-chunk-21-2.png index 917a7274..8d31c6f4 100644 Binary files a/docs/articles/shared_states_files/figure-html/unnamed-chunk-21-2.png and b/docs/articles/shared_states_files/figure-html/unnamed-chunk-21-2.png differ diff --git a/docs/articles/shared_states_files/figure-html/unnamed-chunk-21-3.png b/docs/articles/shared_states_files/figure-html/unnamed-chunk-21-3.png index 09f7e839..b77896af 100644 Binary files a/docs/articles/shared_states_files/figure-html/unnamed-chunk-21-3.png and b/docs/articles/shared_states_files/figure-html/unnamed-chunk-21-3.png differ diff --git a/docs/articles/shared_states_files/figure-html/unnamed-chunk-21-4.png b/docs/articles/shared_states_files/figure-html/unnamed-chunk-21-4.png index f6acd02e..64256cb0 100644 Binary files a/docs/articles/shared_states_files/figure-html/unnamed-chunk-21-4.png and b/docs/articles/shared_states_files/figure-html/unnamed-chunk-21-4.png differ diff --git a/docs/articles/shared_states_files/figure-html/unnamed-chunk-22-1.png b/docs/articles/shared_states_files/figure-html/unnamed-chunk-22-1.png index 8094168e..0acfe901 100644 Binary files a/docs/articles/shared_states_files/figure-html/unnamed-chunk-22-1.png and b/docs/articles/shared_states_files/figure-html/unnamed-chunk-22-1.png differ diff --git a/docs/articles/shared_states_files/figure-html/unnamed-chunk-23-1.png b/docs/articles/shared_states_files/figure-html/unnamed-chunk-23-1.png index 4b1dcd93..6cfd918f 100644 Binary files a/docs/articles/shared_states_files/figure-html/unnamed-chunk-23-1.png and b/docs/articles/shared_states_files/figure-html/unnamed-chunk-23-1.png differ diff --git a/docs/articles/shared_states_files/figure-html/unnamed-chunk-24-1.png b/docs/articles/shared_states_files/figure-html/unnamed-chunk-24-1.png index 5126e752..1a7a9c69 100644 Binary files a/docs/articles/shared_states_files/figure-html/unnamed-chunk-24-1.png and b/docs/articles/shared_states_files/figure-html/unnamed-chunk-24-1.png differ diff --git a/docs/articles/shared_states_files/figure-html/unnamed-chunk-8-1.png b/docs/articles/shared_states_files/figure-html/unnamed-chunk-8-1.png index 6efebe2e..2228fe12 100644 Binary files a/docs/articles/shared_states_files/figure-html/unnamed-chunk-8-1.png and b/docs/articles/shared_states_files/figure-html/unnamed-chunk-8-1.png differ diff --git a/docs/articles/shared_states_files/figure-html/unnamed-chunk-8-2.png b/docs/articles/shared_states_files/figure-html/unnamed-chunk-8-2.png index c4e75afe..64cb21d4 100644 Binary files a/docs/articles/shared_states_files/figure-html/unnamed-chunk-8-2.png and b/docs/articles/shared_states_files/figure-html/unnamed-chunk-8-2.png differ diff --git a/docs/articles/shared_states_files/figure-html/unnamed-chunk-9-1.png b/docs/articles/shared_states_files/figure-html/unnamed-chunk-9-1.png index 2d821f81..45cfc15c 100644 Binary files a/docs/articles/shared_states_files/figure-html/unnamed-chunk-9-1.png and b/docs/articles/shared_states_files/figure-html/unnamed-chunk-9-1.png differ diff --git a/docs/articles/shared_states_files/figure-html/unnamed-chunk-9-2.png b/docs/articles/shared_states_files/figure-html/unnamed-chunk-9-2.png index 39df12fa..c24b1488 100644 Binary files a/docs/articles/shared_states_files/figure-html/unnamed-chunk-9-2.png and b/docs/articles/shared_states_files/figure-html/unnamed-chunk-9-2.png differ diff --git a/docs/articles/shared_states_files/figure-html/unnamed-chunk-9-3.png b/docs/articles/shared_states_files/figure-html/unnamed-chunk-9-3.png index e1ce6553..2bb9461a 100644 Binary files a/docs/articles/shared_states_files/figure-html/unnamed-chunk-9-3.png and b/docs/articles/shared_states_files/figure-html/unnamed-chunk-9-3.png differ diff --git a/docs/articles/time_varying_effects.html b/docs/articles/time_varying_effects.html index 59914c4c..4f59e462 100644 --- a/docs/articles/time_varying_effects.html +++ b/docs/articles/time_varying_effects.html @@ -26,7 +26,7 @@ mvgam - 1.0.9 + 1.0.91