diff --git a/NAMESPACE b/NAMESPACE index 9cc52ceb..8c08ea78 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -186,6 +186,7 @@ importFrom(stats,cor) importFrom(stats,cov) importFrom(stats,cov2cor) importFrom(stats,dbeta) +importFrom(stats,dbinom) importFrom(stats,density) importFrom(stats,dgamma) importFrom(stats,dlnorm) @@ -208,6 +209,7 @@ importFrom(stats,na.fail) importFrom(stats,na.pass) importFrom(stats,pacf) importFrom(stats,pbeta) +importFrom(stats,pbinom) importFrom(stats,pgamma) importFrom(stats,plnorm) importFrom(stats,plogis) @@ -216,6 +218,7 @@ importFrom(stats,poisson) importFrom(stats,ppois) importFrom(stats,predict) importFrom(stats,printCoefmat) +importFrom(stats,qbinom) importFrom(stats,qcauchy) importFrom(stats,qnorm) importFrom(stats,qqline) diff --git a/R/families.R b/R/families.R index 73904953..7711267d 100644 --- a/R/families.R +++ b/R/families.R @@ -1,5 +1,7 @@ #' Supported mvgam families -#' @importFrom stats make.link dgamma pgamma rgamma qnorm plnorm runif pbeta dlnorm dpois pnorm ppois plogis gaussian poisson Gamma dnbinom rnbinom dnorm dbeta +#' @importFrom stats make.link dgamma pgamma rgamma qnorm plnorm runif pbeta dlnorm dpois +#' @importFrom stats pnorm ppois plogis gaussian poisson Gamma dnbinom rnbinom dnorm dbeta +#' @importFrom stats rbinom pbinom dbinom qbinom qlogis #' @importFrom brms lognormal student rstudent_t qstudent_t dstudent_t pstudent_t #' @param link a specification for the family link function. At present these cannot #' be changed diff --git a/R/globals.R b/R/globals.R index 89eb33e6..0bb42e71 100644 --- a/R/globals.R +++ b/R/globals.R @@ -10,4 +10,4 @@ utils::globalVariables(c("y", "year", "smooth_vals", "smooth_num", "term", "data_test", "object", "row_num", "trends_test", "trend", "trend_series", "trend_y", ".", "gam", "group", "mod", "row_id", "byvar", "direction", - "index..time..index")) + "index..time..index", "trend_test")) diff --git a/README.Rmd b/README.Rmd index e37f6827..df593d54 100644 --- a/README.Rmd +++ b/README.Rmd @@ -26,6 +26,7 @@ The goal of `mvgam` is to fit Bayesian Dynamic Generalized Additive Models to ti ## Resources A series of [vignettes cover data formatting, forecasting and several extended case studies of DGAMs](https://nicholasjclark.github.io/mvgam/){target="_blank"}. A number of other examples have also been compiled: +* [Ecological Forecasting with Dynamic Generalized Additive Models](https://www.youtube.com/watch?v=0zZopLlomsQ){target="_blank"} * [mvgam case study 1: model comparison and data assimilation](https://rpubs.com/NickClark47/mvgam){target="_blank"} * [mvgam case study 2: multivariate models](https://rpubs.com/NickClark47/mvgam2){target="_blank"} * [mvgam case study 3: distributed lag models](https://rpubs.com/NickClark47/mvgam3){target="_blank"} @@ -212,6 +213,7 @@ plot(lynx_mvgam, type = 'residuals') * `poisson()` for count data * `nb()` for overdispersed count data * `tweedie()` for overdispersed count data +* `nmix()` for count data with imperfect detection Note that only `poisson()`, `nb()`, and `tweedie()` are available if using `JAGS`. All families, apart from `tweedie()`, are supported if using `Stan`. See `??mvgam_families` for more information. Below is a simple example for simulating and modelling proportional data with `Beta` observations over a set of seasonal series with independent Gaussian Process dynamic trends: ```{r beta_sim, message=FALSE, warning=FALSE, fig.width=6.5, fig.height=6.5, dpi=160} diff --git a/README.md b/README.md index 6230bd4d..19c7be45 100644 --- a/README.md +++ b/README.md @@ -24,6 +24,9 @@ target="_blank">vignettes cover data formatting, forecasting and several extended case studies of DGAMs. A number of other examples have also been compiled: +- Ecological Forecasting with Dynamic Generalized Additive + Models - mvgam case study 1: model comparison and data assimilation - mvgam @@ -314,29 +317,29 @@ summary(lynx_mvgam) #> #> #> GAM coefficient (beta) estimates: -#> 2.5% 50% 97.5% Rhat n_eff -#> (Intercept) 5.900 6.600 7.00 1.01 302 -#> s(season).1 -0.610 0.035 0.69 1.00 910 -#> s(season).2 -0.180 0.840 1.90 1.00 341 -#> s(season).3 -0.031 1.300 2.50 1.00 313 -#> s(season).4 -0.470 0.420 1.40 1.00 901 -#> s(season).5 -1.300 -0.180 0.91 1.00 385 -#> s(season).6 -1.100 -0.057 1.10 1.00 534 -#> s(season).7 -0.700 0.390 1.40 1.00 631 -#> s(season).8 -0.970 0.340 1.90 1.00 282 -#> s(season).9 -1.100 -0.230 0.75 1.00 392 -#> s(season).10 -1.300 -0.680 -0.01 1.01 548 +#> 2.5% 50% 97.5% Rhat n_eff +#> (Intercept) 6.10 6.600 7.000 1.00 355 +#> s(season).1 -0.59 0.046 0.710 1.00 1022 +#> s(season).2 -0.26 0.770 1.800 1.00 443 +#> s(season).3 -0.12 1.100 2.400 1.00 399 +#> s(season).4 -0.51 0.410 1.300 1.00 890 +#> s(season).5 -1.20 -0.130 0.950 1.01 503 +#> s(season).6 -1.00 -0.011 1.100 1.01 699 +#> s(season).7 -0.71 0.340 1.400 1.00 711 +#> s(season).8 -0.92 0.180 1.800 1.00 371 +#> s(season).9 -1.10 -0.290 0.710 1.00 476 +#> s(season).10 -1.30 -0.660 0.027 1.00 595 #> #> Approximate significance of GAM observation smooths: #> edf Chi.sq p-value -#> s(season) 4.43 19440 0.24 +#> s(season) 5.08 17751 0.25 #> #> Latent trend AR parameter estimates: #> 2.5% 50% 97.5% Rhat n_eff -#> ar1[1] 0.71 1.10 1.400 1 594 -#> ar2[1] -0.81 -0.38 0.074 1 1367 -#> ar3[1] -0.47 -0.11 0.340 1 401 -#> sigma[1] 0.40 0.50 0.640 1 1125 +#> ar1[1] 0.74 1.10 1.400 1 635 +#> ar2[1] -0.84 -0.40 0.062 1 1514 +#> ar3[1] -0.47 -0.13 0.290 1 540 +#> sigma[1] 0.40 0.50 0.640 1 1154 #> #> Stan MCMC diagnostics: #> n_eff / iter looks reasonable for all parameters @@ -345,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 Wed Dec 06 5:18:32 PM 2023. +#> Samples were drawn using NUTS(diag_e) at Mon Jan 15 8:57:57 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) @@ -450,6 +453,7 @@ plotted on the outcome scale, for example: ``` r require(ggplot2) #> Loading required package: ggplot2 +#> Warning: package 'ggplot2' was built under R version 4.3.2 plot_predictions(lynx_mvgam, condition = 'season', points = 0.5) + theme_classic() ``` @@ -466,7 +470,7 @@ plot(lynx_mvgam, type = 'forecast', newdata = lynx_test) #> Out of sample CRPS: - #> [1] 2832.494 + #> [1] 2892.767 And the estimated latent trend component, again using the more flexible `plot_mvgam_...()` option to show first derivatives of the estimated @@ -524,6 +528,7 @@ families. Currently, the package can handle data for the following: - `poisson()` for count data - `nb()` for overdispersed count data - `tweedie()` for overdispersed count data +- `nmix()` for count data with imperfect detection Note that only `poisson()`, `nb()`, and `tweedie()` are available if using `JAGS`. All families, apart from `tweedie()`, are supported if @@ -587,41 +592,41 @@ summary(mod, include_betas = FALSE) #> #> Observation precision parameter estimates: #> 2.5% 50% 97.5% Rhat n_eff -#> phi[1] 5.5 8.3 12 1 1641 -#> phi[2] 5.6 8.6 13 1 1725 -#> phi[3] 5.8 8.5 12 1 2063 +#> 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 #> #> GAM coefficient (beta) estimates: -#> 2.5% 50% 97.5% Rhat n_eff -#> (Intercept) -0.24 0.18 0.45 1 799 +#> 2.5% 50% 97.5% Rhat n_eff +#> (Intercept) -0.15 0.2 0.46 1.01 891 #> #> Approximate significance of GAM observation smooths: #> edf Chi.sq p-value -#> s(season) 5.00 16.30 1.1e-05 *** -#> s(season):seriesseries_1 3.83 0.12 0.97 -#> s(season):seriesseries_2 3.72 0.08 0.99 -#> s(season):seriesseries_3 3.78 0.83 0.57 +#> 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 #> --- #> 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.046 0.42 0.95 1 573 -#> alpha_gp[2] 0.350 0.71 1.30 1 1098 -#> alpha_gp[3] 0.180 0.48 0.97 1 1077 -#> rho_gp[1] 1.200 3.90 13.00 1 1487 -#> rho_gp[2] 1.700 7.30 34.00 1 399 -#> rho_gp[3] 1.400 4.90 20.00 1 834 +#> 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 #> #> Stan MCMC diagnostics: #> n_eff / iter looks reasonable for all parameters #> Rhat looks reasonable for all parameters -#> 4 of 2000 iterations ended with a divergence (0.2%) +#> 3 of 2000 iterations ended with a divergence (0.15%) #> *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 Wed Dec 06 5:19:47 PM 2023. +#> Samples were drawn using NUTS(diag_e) at Mon Jan 15 8:59:01 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) @@ -637,11 +642,11 @@ for(i in 1:3){ ``` #> Out of sample CRPS: - #> [1] 2.09906 + #> [1] 2.079879 #> Out of sample CRPS: - #> [1] 1.843141 + #> [1] 1.843941 #> Out of sample CRPS: - #> [1] 1.754107 + #> [1] 1.728574 diff --git a/docs/articles/mvgam_overview.html b/docs/articles/mvgam_overview.html index c131ab7f..71b41a49 100644 --- a/docs/articles/mvgam_overview.html +++ b/docs/articles/mvgam_overview.html @@ -77,7 +77,7 @@

Nicholas J Clark

-

2024-01-12

+

2024-01-15

Source: vignettes/mvgam_overview.Rmd
mvgam_overview.Rmd
@@ -93,7 +93,9 @@

Dynamic GAMs\(\tilde{\boldsymbol{y}}_{i,t}\) is the +independently of the observation process. An introduction to the package +and some worked examples are also shown in this seminar: Ecological Forecasting with Dynamic Generalized Additive +Models. Briefly, assume \(\tilde{\boldsymbol{y}}_{i,t}\) is the conditional expectation of response variable \(\boldsymbol{i}\) at time \(\boldsymbol{t}\). Assuming \(\boldsymbol{y_i}\) is drawn from an exponential distribution with an invertible link function, the linear predictor for a multivariate Dynamic GAM can be written as:

@@ -443,17 +445,17 @@

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

Notice how we have four different time series in these simulated data, but we do not spread the outcome values into different columns. @@ -647,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.9, 0.59); -## 2 sigma_raw ~ student_t(3, 0, 2.5); sigma_raw ~ exponential(0.61); +## 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); ## new_lowerbound new_upperbound ## 1 NA NA ## 2 NA NA @@ -684,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.80 2.1 2.3 1 2467 -## s(year_fac).2 2.50 2.7 2.8 1 2762 -## s(year_fac).3 3.00 3.1 3.2 1 3040 -## s(year_fac).4 3.10 3.3 3.4 1 3381 -## s(year_fac).5 1.90 2.1 2.3 1 3455 -## s(year_fac).6 1.50 1.7 2.0 1 3422 -## s(year_fac).7 1.80 2.0 2.3 1 3417 -## s(year_fac).8 2.80 3.0 3.1 1 3207 -## s(year_fac).9 3.10 3.3 3.4 1 3011 -## s(year_fac).10 2.60 2.8 2.9 1 2673 -## s(year_fac).11 2.90 3.1 3.2 1 3767 -## s(year_fac).12 3.10 3.2 3.3 1 2400 -## s(year_fac).13 2.00 2.3 2.5 1 2857 -## s(year_fac).14 2.50 2.6 2.8 1 2272 -## s(year_fac).15 1.90 2.2 2.4 1 2637 -## s(year_fac).16 1.90 2.1 2.3 1 2868 -## s(year_fac).17 -0.31 1.1 1.9 1 583 +## 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 ## ## GAM group-level estimates: ## 2.5% 50% 97.5% Rhat n_eff -## mean(s(year_fac)) 2.00 2.40 2.7 1.02 229 -## sd(s(year_fac)) 0.46 0.68 1.0 1.01 302 +## 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 ## ## Approximate significance of GAM observation smooths: ## edf Chi.sq p-value -## s(year_fac) 13.4 23802 <2e-16 *** +## s(year_fac) 13.2 24457 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## @@ -721,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 Fri Jan 12 3:51:59 PM 2024. +## Samples were drawn using NUTS(diag_e) at Mon Jan 15 9:06:57 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) @@ -738,23 +740,23 @@

GLMs with temporal random effectsdplyr::glimpse(beta_post) 
## Rows: 2,000
 ## Columns: 17
-## $ `s(year_fac).1`  <dbl> 2.06435, 2.18026, 2.10259, 2.15858, 2.08945, 2.03449,…
-## $ `s(year_fac).2`  <dbl> 2.66297, 2.76049, 2.58599, 2.84501, 2.65126, 2.77843,…
-## $ `s(year_fac).3`  <dbl> 2.99574, 3.07175, 3.07301, 3.05926, 3.09694, 2.96840,…
-## $ `s(year_fac).4`  <dbl> 3.24798, 3.28896, 3.30628, 3.25715, 3.22578, 3.27270,…
-## $ `s(year_fac).5`  <dbl> 2.13139, 2.07544, 1.97331, 2.08521, 2.29660, 2.08273,…
-## $ `s(year_fac).6`  <dbl> 1.42739, 2.02722, 1.71548, 1.92531, 1.73130, 1.49182,…
-## $ `s(year_fac).7`  <dbl> 1.93186, 2.21940, 2.11423, 2.18194, 2.03794, 2.12441,…
-## $ `s(year_fac).8`  <dbl> 2.94445, 3.01923, 3.01388, 2.95011, 3.01846, 2.93042,…
-## $ `s(year_fac).9`  <dbl> 3.19531, 3.27323, 3.38481, 3.30585, 3.13491, 3.18199,…
-## $ `s(year_fac).10` <dbl> 2.74548, 2.88965, 2.78905, 2.90565, 2.72434, 2.75425,…
-## $ `s(year_fac).11` <dbl> 3.01341, 3.09803, 3.14241, 3.10061, 3.05058, 3.02510,…
-## $ `s(year_fac).12` <dbl> 3.13730, 3.30920, 3.24322, 3.20736, 3.18905, 3.23454,…
-## $ `s(year_fac).13` <dbl> 2.18607, 2.26127, 2.15356, 2.14097, 2.43450, 2.22882,…
-## $ `s(year_fac).14` <dbl> 2.38431, 2.50610, 2.49737, 2.66538, 2.47568, 2.75505,…
-## $ `s(year_fac).15` <dbl> 2.22367, 2.34576, 2.24644, 2.22469, 2.27891, 2.17492,…
-## $ `s(year_fac).16` <dbl> 2.27256, 2.11722, 1.97474, 2.07540, 2.08201, 1.86904,…
-## $ `s(year_fac).17` <dbl> 1.189810, 1.151540, 0.804486, 1.526610, 0.851250, 1.1…
+## $ `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…

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

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

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

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

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

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

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

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] 180.2378
+## [1] 182.1551

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 @@ -975,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] 10 7 11 7 8 11 5 9 7 11 ... +## ..$ PP: num [1:2000, 1:160] 11 10 4 6 7 3 12 5 4 8 ... ## .. ..- 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 8 6 14 12 11 11 9 6 10 ... +## ..$ PP: num [1:2000, 1:39] 9 9 16 3 9 8 11 20 9 10 ... ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : NULL ## .. .. ..$ : chr [1:39] "ypred[161,1]" "ypred[162,1]" "ypred[163,1]" "ypred[164,1]" ... @@ -1037,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 1941 -## s(year_fac).1 1.10 1.40 1.70 1 2704 -## s(year_fac).2 1.80 2.00 2.20 1 2165 -## s(year_fac).3 2.20 2.40 2.60 1 2034 -## s(year_fac).4 2.30 2.50 2.70 1 2187 -## s(year_fac).5 1.20 1.40 1.60 1 2687 -## s(year_fac).6 1.00 1.30 1.50 1 2437 -## s(year_fac).7 1.10 1.40 1.70 1 3079 -## s(year_fac).8 2.10 2.30 2.50 1 2295 -## s(year_fac).9 2.70 2.90 3.00 1 2301 -## s(year_fac).10 2.00 2.20 2.40 1 2581 -## s(year_fac).11 2.30 2.40 2.60 1 2190 -## s(year_fac).12 2.50 2.70 2.80 1 2097 -## s(year_fac).13 1.40 1.60 1.80 1 3245 -## s(year_fac).14 0.56 2.00 3.40 1 1688 -## s(year_fac).15 0.73 2.00 3.40 1 1698 -## s(year_fac).16 0.52 2.00 3.20 1 1581 -## s(year_fac).17 0.66 2.00 3.30 1 1633 +## 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 +## 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 ## ## GAM group-level estimates: -## 2.5% 50% 97.5% Rhat n_eff -## mean(s(year_fac)) 1.6 2.0 2.3 1.01 394 -## sd(s(year_fac)) 0.4 0.6 1.0 1.00 493 +## 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 ## ## Approximate significance of GAM observation smooths: ## edf Chi.sq p-value -## s(year_fac) 10.8 2810 <2e-16 *** +## s(year_fac) 10.3 2802 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## @@ -1074,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 Fri Jan 12 3:53:06 PM 2024. +## Samples were drawn using NUTS(diag_e) at Mon Jan 15 9:07:58 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) @@ -1085,24 +1087,24 @@

Adding predictors as “fixed” eff 
 coef(model2)
##                     2.5%       50%     97.5% Rhat n_eff
-## ndvi           0.3219343 0.3896115 0.4561312    1  1941
-## s(year_fac).1  1.1213810 1.4014750 1.6763905    1  2704
-## s(year_fac).2  1.8000682 2.0033650 2.2080812    1  2165
-## s(year_fac).3  2.1762575 2.3832750 2.5682495    1  2034
-## s(year_fac).4  2.3162657 2.5098600 2.6993935    1  2187
-## s(year_fac).5  1.1848087 1.4253400 1.6454560    1  2687
-## s(year_fac).6  1.0015992 1.2699300 1.5162433    1  2437
-## s(year_fac).7  1.1421785 1.4162100 1.6830155    1  3079
-## s(year_fac).8  2.0869212 2.2735250 2.4629585    1  2295
-## s(year_fac).9  2.7231583 2.8611750 2.9968915    1  2301
-## s(year_fac).10 1.9826495 2.1869100 2.3888603    1  2581
-## s(year_fac).11 2.2645230 2.4408950 2.6059805    1  2190
-## s(year_fac).12 2.5381555 2.6944350 2.8382668    1  2097
-## s(year_fac).13 1.3885932 1.6143500 1.8470717    1  3245
-## s(year_fac).14 0.5566057 1.9835550 3.3532650    1  1688
-## s(year_fac).15 0.7313315 1.9786150 3.3730625    1  1698
-## s(year_fac).16 0.5209962 1.9957750 3.1762922    1  1581
-## s(year_fac).17 0.6614678 1.9759650 3.2837728    1  1633
+## 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

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

@@ -1118,24 +1120,24 @@ 
Adding predictors as “fixed” eff
 dplyr::glimpse(beta_post)
## Rows: 2,000
 ## Columns: 18
-## $ ndvi             <dbl> 0.355167, 0.354315, 0.400233, 0.388748, 0.341074, 0.4…
-## $ `s(year_fac).1`  <dbl> 1.45330, 1.39405, 1.45254, 1.11547, 1.66909, 1.17323,…
-## $ `s(year_fac).2`  <dbl> 2.05506, 2.07679, 1.97119, 1.99876, 2.01560, 2.07635,…
-## $ `s(year_fac).3`  <dbl> 2.59536, 2.35567, 2.37899, 2.32072, 2.48880, 2.40466,…
-## $ `s(year_fac).4`  <dbl> 2.63194, 2.56439, 2.55642, 2.43839, 2.67447, 2.47027,…
-## $ `s(year_fac).5`  <dbl> 1.54866, 1.44661, 1.48310, 1.35734, 1.53211, 1.42122,…
-## $ `s(year_fac).6`  <dbl> 1.28962, 1.44256, 1.13824, 1.43169, 1.27521, 1.37027,…
-## $ `s(year_fac).7`  <dbl> 1.60251, 1.39342, 1.47103, 1.40323, 1.46742, 1.40967,…
-## $ `s(year_fac).8`  <dbl> 2.33166, 2.34704, 2.19045, 2.30872, 2.34791, 2.26251,…
-## $ `s(year_fac).9`  <dbl> 2.96146, 2.87900, 2.85271, 2.76788, 2.92536, 2.87675,…
-## $ `s(year_fac).10` <dbl> 2.16140, 2.41945, 1.98341, 2.39724, 2.20351, 2.18907,…
-## $ `s(year_fac).11` <dbl> 2.51912, 2.52120, 2.45377, 2.36322, 2.59581, 2.37504,…
-## $ `s(year_fac).12` <dbl> 2.75224, 2.76068, 2.70921, 2.68268, 2.78127, 2.68005,…
-## $ `s(year_fac).13` <dbl> 1.63046, 1.77451, 1.48580, 1.66704, 1.72238, 1.52178,…
-## $ `s(year_fac).14` <dbl> 1.496390, 1.016750, 1.710010, 2.178540, 2.644870, 2.6…
-## $ `s(year_fac).15` <dbl> 1.871760, 2.004820, 1.210480, 2.188060, 2.624220, 2.5…
-## $ `s(year_fac).16` <dbl> 2.634900, 2.017970, 2.073360, 2.412940, 2.959150, 2.1…
-## $ `s(year_fac).17` <dbl> 1.91701, 1.50308, 1.17690, 1.54644, 1.92542, 2.15811,…
+## $ 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…

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 @@ -1272,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 976 -## ndvi 0.26 0.33 0.400 1.00 1005 -## s(time).1 -2.20 -1.10 -0.066 1.01 481 -## s(time).2 0.45 1.30 2.200 1.00 433 -## s(time).3 -0.49 0.46 1.400 1.01 377 -## s(time).4 1.60 2.50 3.400 1.01 376 -## s(time).5 -1.20 -0.20 0.730 1.01 373 -## s(time).6 -0.60 0.36 1.400 1.01 396 -## s(time).7 -1.50 -0.50 0.480 1.00 436 -## s(time).8 0.59 1.50 2.400 1.01 377 -## s(time).9 1.10 2.10 3.000 1.01 368 -## s(time).10 -0.37 0.55 1.500 1.01 370 -## s(time).11 0.80 1.70 2.700 1.01 372 -## s(time).12 0.65 1.50 2.400 1.01 365 -## s(time).13 -1.20 -0.30 0.560 1.00 601 -## s(time).14 -7.20 -4.20 -1.000 1.00 509 +## (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 ## ## Approximate significance of GAM observation smooths: ## edf Chi.sq p-value -## s(time) 9.76 679 <2e-16 *** +## s(time) 8.57 797 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## @@ -1302,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 Fri Jan 12 3:53:45 PM 2024. +## Samples were drawn using NUTS(diag_e) at Mon Jan 15 9:08:44 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) @@ -1429,7 +1431,7 @@

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

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

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 @@ -1512,33 +1514,32 @@

Latent dynamics in mvgam## ## GAM coefficient (beta) estimates: ## 2.5% 50% 97.5% Rhat n_eff -## (Intercept) 0.580 1.9000 2.700 1.17 36 -## s(ndvi).1 -0.078 0.0069 0.120 1.01 1388 -## s(ndvi).2 -0.120 0.0150 0.210 1.00 915 -## s(ndvi).3 -0.038 -0.0012 0.038 1.00 2092 -## s(ndvi).4 -0.200 0.1000 0.900 1.01 561 -## s(ndvi).5 -0.049 0.1600 0.350 1.00 489 +## (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 ## ## Approximate significance of GAM observation smooths: -## edf Chi.sq p-value -## s(ndvi) 1.05 89.1 0.055 . +## edf Chi.sq p-value +## s(ndvi) 0.978 95.5 0.088 . ## --- ## 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.70 0.81 0.93 1.05 76 -## sigma[1] 0.67 0.80 0.95 1.00 555 +## ar1[1] 0.69 0.81 0.94 1.01 173 +## sigma[1] 0.67 0.80 0.96 1.01 518 ## ## Stan MCMC diagnostics: ## n_eff / iter looks reasonable for all parameters -## Rhats above 1.05 found for 131 parameters -## *Diagnose further to investigate why the chains have not mixed +## Rhat looks reasonable for all parameters ## 0 of 2000 iterations ended with a divergence (0%) ## 0 of 2000 iterations saturated the maximum tree depth of 12 (0%) ## E-FMI indicated no pathological behavior ## -## Samples were drawn using NUTS(diag_e) at Fri Jan 12 3:54:58 PM 2024. +## Samples were drawn using NUTS(diag_e) at Mon Jan 15 9:10:06 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) @@ -1553,8 +1554,8 @@

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

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

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

 plot(model4, type = 'trend', newdata = data_test)
@@ -1569,7 +1570,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.0      66.5
+## model3 -560.2 66.5

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

@@ -1584,7 +1585,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] -136.3461
+
## [1] -135.6337

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 f099c576..56b72cb6 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 f7c968dd..bcdffa60 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 3d4ab4ea..1109c200 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 f69db9ae..a001b44e 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 7300d72c..a0d871f8 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 c9993f7a..08683df5 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 26da8d53..d794301d 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 7a93a32c..7c90e2d2 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 c5810b2a..3a4bedc4 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 b6536a5e..d23650fc 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 4bed7b14..90fc8860 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 c5810b2a..3a4bedc4 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 4a50ce1b..b8408084 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 beab5884..52974248 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 7dbf0d61..ae02aec2 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 e1c1981d..036dd03e 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 f69b9b64..a244e649 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 2a52a1a3..66a4b17b 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 3dc692a2..84cd11b3 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 50f0a428..b289cefe 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 bad44e54..f5219b5c 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 2eb5a66f..275429b3 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 7b9f36f9..29cbf193 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 d0b5a580..1cb85882 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 eb6ffc87..c03916e2 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/index.html b/docs/index.html index 575f98d4..f1389398 100644 --- a/docs/index.html +++ b/docs/index.html @@ -87,7 +87,7 @@

mvgam

MultiVariate (Dynamic) Generalized Addivite Models

-

The goal of mvgam is to use a Bayesian framework to estimate parameters of Dynamic Generalized Additive Models (DGAMs) for time series with dynamic trend components. The package provides an interface to fit Bayesian DGAMs using either JAGS or Stan as the backend, but note that users are strongly encouraged to opt for Stan over JAGS. The formula syntax is based on that of the package mgcv to provide a familiar GAM modelling interface. The motivation for the package and some of its primary objectives are described in detail by Clark & Wells 2022 (published in Methods in Ecology and Evolution).

+

The goal of mvgam is to use a Bayesian framework to estimate parameters of Dynamic Generalized Additive Models (DGAMs) for time series with dynamic trend components. The package provides an interface to fit Bayesian DGAMs using either JAGS or Stan as the backend, but note that users are strongly encouraged to opt for Stan over JAGS. The formula syntax is based on that of the package mgcv to provide a familiar GAM modelling interface. The motivation for the package and some of its primary objectives are described in detail by Clark & Wells 2022 (published in Methods in Ecology and Evolution). An introduction to the package and some worked examples are also shown in this seminar: Ecological Forecasting with Dynamic Generalized Additive Models.

Installation @@ -160,8 +160,9 @@

Vignettes

Other resources

-

A number of case studies have been compiled to highlight how DGAMs can be estimated using MCMC sampling. These are hosted currently on RPubs at the following links:

+

A number of case studies have been compiled to highlight how DGAMs can be estimated using MCMC sampling:

    +
  • Ecological Forecasting with Dynamic Generalized Additive Models
  • mvgam case study 1: model comparison and data assimilation
  • mvgam case study 2: multivariate models
  • mvgam case study 3: distributed lag models
    • diff --git a/docs/reference/figures/README-beta_fc-1.png b/docs/reference/figures/README-beta_fc-1.png index fb343347..2fda386d 100644 Binary files a/docs/reference/figures/README-beta_fc-1.png and b/docs/reference/figures/README-beta_fc-1.png differ diff --git a/docs/reference/figures/README-unnamed-chunk-13-1.png b/docs/reference/figures/README-unnamed-chunk-13-1.png index a2c84c34..8e55a768 100644 Binary files a/docs/reference/figures/README-unnamed-chunk-13-1.png and b/docs/reference/figures/README-unnamed-chunk-13-1.png differ diff --git a/docs/reference/figures/README-unnamed-chunk-14-1.png b/docs/reference/figures/README-unnamed-chunk-14-1.png index 0fec3676..60f59430 100644 Binary files a/docs/reference/figures/README-unnamed-chunk-14-1.png and b/docs/reference/figures/README-unnamed-chunk-14-1.png differ diff --git a/docs/reference/figures/README-unnamed-chunk-15-1.png b/docs/reference/figures/README-unnamed-chunk-15-1.png index 35f3787c..0b1e2b90 100644 Binary files a/docs/reference/figures/README-unnamed-chunk-15-1.png and b/docs/reference/figures/README-unnamed-chunk-15-1.png differ diff --git a/docs/reference/figures/README-unnamed-chunk-16-1.png b/docs/reference/figures/README-unnamed-chunk-16-1.png index 7c34523c..6fea81f3 100644 Binary files a/docs/reference/figures/README-unnamed-chunk-16-1.png and b/docs/reference/figures/README-unnamed-chunk-16-1.png differ diff --git a/docs/reference/figures/README-unnamed-chunk-17-1.png b/docs/reference/figures/README-unnamed-chunk-17-1.png index a6688002..5b6a05c1 100644 Binary files a/docs/reference/figures/README-unnamed-chunk-17-1.png and b/docs/reference/figures/README-unnamed-chunk-17-1.png differ diff --git a/docs/reference/figures/README-unnamed-chunk-18-1.png b/docs/reference/figures/README-unnamed-chunk-18-1.png index 52b8c609..87c39ae7 100644 Binary files a/docs/reference/figures/README-unnamed-chunk-18-1.png and b/docs/reference/figures/README-unnamed-chunk-18-1.png differ diff --git a/docs/reference/figures/README-unnamed-chunk-19-1.png b/docs/reference/figures/README-unnamed-chunk-19-1.png index 9ec2a55f..aacb08f5 100644 Binary files a/docs/reference/figures/README-unnamed-chunk-19-1.png and b/docs/reference/figures/README-unnamed-chunk-19-1.png differ diff --git a/docs/reference/figures/README-unnamed-chunk-20-1.png b/docs/reference/figures/README-unnamed-chunk-20-1.png index 24ac38b0..5a6df81b 100644 Binary files a/docs/reference/figures/README-unnamed-chunk-20-1.png and b/docs/reference/figures/README-unnamed-chunk-20-1.png differ diff --git a/docs/reference/figures/README-unnamed-chunk-21-1.png b/docs/reference/figures/README-unnamed-chunk-21-1.png index 7239bca0..65921406 100644 Binary files a/docs/reference/figures/README-unnamed-chunk-21-1.png and b/docs/reference/figures/README-unnamed-chunk-21-1.png differ diff --git a/docs/reference/figures/README-unnamed-chunk-22-1.png b/docs/reference/figures/README-unnamed-chunk-22-1.png index fe812642..ede996b0 100644 Binary files a/docs/reference/figures/README-unnamed-chunk-22-1.png and b/docs/reference/figures/README-unnamed-chunk-22-1.png differ diff --git a/docs/reference/figures/README-unnamed-chunk-23-1.png b/docs/reference/figures/README-unnamed-chunk-23-1.png index 9abc7f94..83596c37 100644 Binary files a/docs/reference/figures/README-unnamed-chunk-23-1.png and b/docs/reference/figures/README-unnamed-chunk-23-1.png differ diff --git a/docs/reference/figures/README-unnamed-chunk-24-1.png b/docs/reference/figures/README-unnamed-chunk-24-1.png index 2a03f3fd..d487efc7 100644 Binary files a/docs/reference/figures/README-unnamed-chunk-24-1.png and b/docs/reference/figures/README-unnamed-chunk-24-1.png differ diff --git a/index.Rmd b/index.Rmd index 82185a05..8b9f990b 100644 --- a/index.Rmd +++ b/index.Rmd @@ -7,7 +7,7 @@ output: github_document ## mvgam **M**ulti**V**ariate (Dynamic) **G**eneralized **A**ddivite **M**odels -The goal of `mvgam` is to use a Bayesian framework to estimate parameters of Dynamic Generalized Additive Models (DGAMs) for time series with dynamic trend components. The package provides an interface to fit Bayesian DGAMs using either `JAGS` or `Stan` as the backend, but note that users are strongly encouraged to opt for `Stan` over `JAGS`. The formula syntax is based on that of the package `mgcv` to provide a familiar GAM modelling interface. The motivation for the package and some of its primary objectives are described in detail by [Clark & Wells 2022](https://besjournals.onlinelibrary.wiley.com/doi/10.1111/2041-210X.13974) (published in *Methods in Ecology and Evolution*). +The goal of `mvgam` is to use a Bayesian framework to estimate parameters of Dynamic Generalized Additive Models (DGAMs) for time series with dynamic trend components. The package provides an interface to fit Bayesian DGAMs using either `JAGS` or `Stan` as the backend, but note that users are strongly encouraged to opt for `Stan` over `JAGS`. The formula syntax is based on that of the package `mgcv` to provide a familiar GAM modelling interface. The motivation for the package and some of its primary objectives are described in detail by [Clark & Wells 2022](https://besjournals.onlinelibrary.wiley.com/doi/10.1111/2041-210X.13974) (published in *Methods in Ecology and Evolution*). An introduction to the package and some worked examples are also shown in this seminar: [Ecological Forecasting with Dynamic Generalized Additive Models](https://www.youtube.com/watch?v=0zZopLlomsQ){target="_blank"}. ## Installation Install the development version from `GitHub` using: @@ -74,8 +74,9 @@ Various `S3` functions can be used to inspect parameter estimates, plot smooth f You can set `build_vignettes = TRUE` when installing with either `devtools::install_github` or `remotes::install_github`, but be aware this will slow down the installation drastically. Instead, you can always access the vignette htmls online at [https://nicholasjclark.github.io/mvgam/articles/](https://nicholasjclark.github.io/mvgam/articles/) ## Other resources -A number of case studies have been compiled to highlight how DGAMs can be estimated using MCMC sampling. These are hosted currently on `RPubs` at the following links: +A number of case studies have been compiled to highlight how DGAMs can be estimated using MCMC sampling: +* [Ecological Forecasting with Dynamic Generalized Additive Models](https://www.youtube.com/watch?v=0zZopLlomsQ){target="_blank"} * [mvgam case study 1: model comparison and data assimilation](https://rpubs.com/NickClark47/mvgam) * [mvgam case study 2: multivariate models](https://rpubs.com/NickClark47/mvgam2) * [mvgam case study 3: distributed lag models](https://rpubs.com/NickClark47/mvgam3) diff --git a/index.md b/index.md index 4b0a9d72..a087194b 100644 --- a/index.md +++ b/index.md @@ -14,7 +14,11 @@ The formula syntax is based on that of the package `mgcv` to provide a familiar GAM modelling interface. The motivation for the package and some of its primary objectives are described in detail by [Clark & Wells 2022](https://besjournals.onlinelibrary.wiley.com/doi/10.1111/2041-210X.13974) -(published in *Methods in Ecology and Evolution*). +(published in *Methods in Ecology and Evolution*). An introduction to +the package and some worked examples are also shown in this seminar: +Ecological Forecasting with Dynamic Generalized Additive +Models. ## Installation @@ -128,9 +132,11 @@ always access the vignette htmls online at ## Other resources A number of case studies have been compiled to highlight how DGAMs can -be estimated using MCMC sampling. These are hosted currently on `RPubs` -at the following links: +be estimated using MCMC sampling: +- Ecological Forecasting with Dynamic Generalized Additive + Models - [mvgam case study 1: model comparison and data assimilation](https://rpubs.com/NickClark47/mvgam) - [mvgam case study 2: multivariate diff --git a/man/figures/README-beta_fc-1.png b/man/figures/README-beta_fc-1.png index fb343347..2fda386d 100644 Binary files a/man/figures/README-beta_fc-1.png and b/man/figures/README-beta_fc-1.png differ diff --git a/man/figures/README-unnamed-chunk-13-1.png b/man/figures/README-unnamed-chunk-13-1.png index a2c84c34..8e55a768 100644 Binary files a/man/figures/README-unnamed-chunk-13-1.png and b/man/figures/README-unnamed-chunk-13-1.png differ diff --git a/man/figures/README-unnamed-chunk-14-1.png b/man/figures/README-unnamed-chunk-14-1.png index 0fec3676..60f59430 100644 Binary files a/man/figures/README-unnamed-chunk-14-1.png and b/man/figures/README-unnamed-chunk-14-1.png differ diff --git a/man/figures/README-unnamed-chunk-15-1.png b/man/figures/README-unnamed-chunk-15-1.png index 35f3787c..0b1e2b90 100644 Binary files a/man/figures/README-unnamed-chunk-15-1.png and b/man/figures/README-unnamed-chunk-15-1.png differ diff --git a/man/figures/README-unnamed-chunk-16-1.png b/man/figures/README-unnamed-chunk-16-1.png index 7c34523c..6fea81f3 100644 Binary files a/man/figures/README-unnamed-chunk-16-1.png and b/man/figures/README-unnamed-chunk-16-1.png differ diff --git a/man/figures/README-unnamed-chunk-17-1.png b/man/figures/README-unnamed-chunk-17-1.png index a6688002..5b6a05c1 100644 Binary files a/man/figures/README-unnamed-chunk-17-1.png and b/man/figures/README-unnamed-chunk-17-1.png differ diff --git a/man/figures/README-unnamed-chunk-18-1.png b/man/figures/README-unnamed-chunk-18-1.png index 52b8c609..87c39ae7 100644 Binary files a/man/figures/README-unnamed-chunk-18-1.png and b/man/figures/README-unnamed-chunk-18-1.png differ diff --git a/man/figures/README-unnamed-chunk-19-1.png b/man/figures/README-unnamed-chunk-19-1.png index 9ec2a55f..aacb08f5 100644 Binary files a/man/figures/README-unnamed-chunk-19-1.png and b/man/figures/README-unnamed-chunk-19-1.png differ diff --git a/man/figures/README-unnamed-chunk-20-1.png b/man/figures/README-unnamed-chunk-20-1.png index 24ac38b0..5a6df81b 100644 Binary files a/man/figures/README-unnamed-chunk-20-1.png and b/man/figures/README-unnamed-chunk-20-1.png differ diff --git a/man/figures/README-unnamed-chunk-21-1.png b/man/figures/README-unnamed-chunk-21-1.png index 7239bca0..65921406 100644 Binary files a/man/figures/README-unnamed-chunk-21-1.png and b/man/figures/README-unnamed-chunk-21-1.png differ diff --git a/man/figures/README-unnamed-chunk-22-1.png b/man/figures/README-unnamed-chunk-22-1.png index fe812642..ede996b0 100644 Binary files a/man/figures/README-unnamed-chunk-22-1.png and b/man/figures/README-unnamed-chunk-22-1.png differ diff --git a/man/figures/README-unnamed-chunk-23-1.png b/man/figures/README-unnamed-chunk-23-1.png index 9abc7f94..83596c37 100644 Binary files a/man/figures/README-unnamed-chunk-23-1.png and b/man/figures/README-unnamed-chunk-23-1.png differ diff --git a/man/figures/README-unnamed-chunk-24-1.png b/man/figures/README-unnamed-chunk-24-1.png index 2a03f3fd..d487efc7 100644 Binary files a/man/figures/README-unnamed-chunk-24-1.png and b/man/figures/README-unnamed-chunk-24-1.png differ diff --git a/src/mvgam.dll b/src/mvgam.dll index 81b15d62..8e6d2a27 100644 Binary files a/src/mvgam.dll and b/src/mvgam.dll differ diff --git a/vignettes/mvgam_overview.Rmd b/vignettes/mvgam_overview.Rmd index f657378b..43f2758e 100644 --- a/vignettes/mvgam_overview.Rmd +++ b/vignettes/mvgam_overview.Rmd @@ -27,7 +27,7 @@ theme_set(theme_bw(base_size = 12, base_family = 'serif')) The purpose of this vignette is to give a general overview of the `mvgam` package and its primary functions. ## Dynamic GAMs -`mvgam` is designed to propagate unobserved temporal processes to capture latent dynamics in the observed time series. This works in a state-space format, with the temporal *trend* evolving independently of the observation process. Briefly, assume $\tilde{\boldsymbol{y}}_{i,t}$ is the conditional expectation of response variable $\boldsymbol{i}$ at time $\boldsymbol{t}$. Assuming $\boldsymbol{y_i}$ is drawn from an exponential distribution with an invertible link function, the linear predictor for a multivariate Dynamic GAM can be written as: +`mvgam` is designed to propagate unobserved temporal processes to capture latent dynamics in the observed time series. This works in a state-space format, with the temporal *trend* evolving independently of the observation process. An introduction to the package and some worked examples are also shown in this seminar: [Ecological Forecasting with Dynamic Generalized Additive Models](https://www.youtube.com/watch?v=0zZopLlomsQ){target="_blank"}. Briefly, assume $\tilde{\boldsymbol{y}}_{i,t}$ is the conditional expectation of response variable $\boldsymbol{i}$ at time $\boldsymbol{t}$. Assuming $\boldsymbol{y_i}$ is drawn from an exponential distribution with an invertible link function, the linear predictor for a multivariate Dynamic GAM can be written as: $$for~i~in~1:N_{series}~...$$ $$for~t~in~1:N_{timepoints}~...$$