diff --git a/docs/articles/mvgam_overview.html b/docs/articles/mvgam_overview.html index 331c738f..8e128e53 100644 --- a/docs/articles/mvgam_overview.html +++ b/docs/articles/mvgam_overview.html @@ -368,8 +368,11 @@

Regression formulae?mvgam_formulae for more guidance.

-

For setting up State-Space models, the ptional process model formula -can be used (see and for guidance on using trend formulae).

+

For setting up State-Space models, the optional process model formula +can be used (see the +State-Space model vignette and the +shared latent states vignette for guidance on using trend +formulae).

Example time series data @@ -381,7 +384,7 @@

Example time series datain this preprint by Ernest et al on the Biorxiv.

+disruptions etc…). You can read about the full sampling protocol in this preprint by Ernest et al on the Biorxiv.

 data("portal_data")

As the data come pre-loaded with the mvgam package, you @@ -433,18 +436,18 @@

Manipulating data for modellingdata <- sim_mvgam(n_series = 4, T = 24) head(data$data_train, 12)

##    y season year   series time
-## 1  3      1    1 series_1    1
-## 2  2      1    1 series_2    1
-## 3  0      1    1 series_3    1
-## 4  1      1    1 series_4    1
-## 5  0      2    1 series_1    2
-## 6  0      2    1 series_2    2
+## 1  0      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  3      2    1 series_1    2
+## 6  1      2    1 series_2    2
 ## 7  1      2    1 series_3    2
-## 8  0      2    1 series_4    2
-## 9  1      3    1 series_1    3
-## 10 1      3    1 series_2    3
-## 11 0      3    1 series_3    3
-## 12 0      3    1 series_4    3
+## 8 1 2 1 series_4 2 +## 9 2 3 1 series_1 3 +## 10 0 3 1 series_2 3 +## 11 2 3 1 series_3 3 +## 12 1 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 @@ -637,8 +640,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.4, 0.98); -## 2 sigma_raw ~ student_t(3, 0, 2.5); sigma_raw ~ exponential(0.79); +## 1 mu_raw ~ std_normal(); mu_raw ~ normal(-0.22, 0.96); +## 2 sigma_raw ~ student_t(3, 0, 2.5); sigma_raw ~ exponential(0.05); ## new_lowerbound new_upperbound ## 1 NA NA ## 2 NA NA @@ -674,45 +677,44 @@

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 2475 -## s(year_fac).2 2.50 2.7 2.8 1 2522 -## s(year_fac).3 3.00 3.1 3.2 1 2959 -## s(year_fac).4 3.10 3.3 3.4 1 2138 -## s(year_fac).5 1.90 2.1 2.3 1 3027 -## s(year_fac).6 1.50 1.7 2.0 1 3344 -## s(year_fac).7 1.80 2.0 2.3 1 2450 -## s(year_fac).8 2.80 3.0 3.1 1 3241 -## s(year_fac).9 3.10 3.3 3.4 1 2615 -## s(year_fac).10 2.60 2.8 2.9 1 2894 -## s(year_fac).11 2.90 3.1 3.2 1 2982 -## s(year_fac).12 3.10 3.2 3.3 1 1720 -## s(year_fac).13 2.00 2.3 2.4 1 2544 -## s(year_fac).14 2.50 2.6 2.8 1 2866 -## s(year_fac).15 1.90 2.2 2.4 1 2872 -## s(year_fac).16 1.90 2.1 2.3 1 3523 -## s(year_fac).17 -0.21 1.0 1.9 1 504 +## 2.5% 50% 97.5% Rhat n_eff +## s(year_fac).1 1.8 2.1 2.3 1.00 2505 +## s(year_fac).2 2.5 2.7 2.8 1.00 3005 +## s(year_fac).3 3.0 3.1 3.2 1.00 2744 +## s(year_fac).4 3.1 3.3 3.4 1.00 2373 +## s(year_fac).5 1.9 2.1 2.3 1.00 2832 +## s(year_fac).6 1.5 1.8 2.0 1.00 2831 +## s(year_fac).7 1.8 2.0 2.3 1.00 2250 +## s(year_fac).8 2.8 3.0 3.1 1.00 3432 +## s(year_fac).9 3.1 3.2 3.4 1.00 2494 +## s(year_fac).10 2.6 2.8 2.9 1.00 3609 +## s(year_fac).11 2.9 3.1 3.2 1.00 2475 +## s(year_fac).12 3.1 3.2 3.3 1.00 2967 +## s(year_fac).13 2.0 2.2 2.5 1.00 2587 +## s(year_fac).14 2.5 2.6 2.8 1.00 2688 +## s(year_fac).15 1.9 2.2 2.4 1.00 2725 +## s(year_fac).16 1.9 2.1 2.3 1.00 2817 +## s(year_fac).17 -0.3 1.1 2.0 1.01 529 ## ## GAM group-level estimates: ## 2.5% 50% 97.5% Rhat n_eff -## mean(s(year_fac)) 2.00 2.40 2.7 1.02 214 -## sd(s(year_fac)) 0.47 0.69 1.2 1.03 199 +## mean(s(year_fac)) 2.00 2.40 2.7 1.02 217 +## sd(s(year_fac)) 0.45 0.67 1.1 1.02 269 ## ## Approximate significance of GAM observation smooths: -## edf Chi.sq p-value -## s(year_fac) 14 24783 <2e-16 *** +## edf Chi.sq p-value +## s(year_fac) 13.6 23562 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Stan MCMC diagnostics: ## n_eff / iter looks reasonable for all parameters ## Rhat looks reasonable for all parameters -## 1 of 2000 iterations ended with a divergence (0.05%) -## *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 Wed Jan 03 9:19:31 AM 2024. +## Samples were drawn using NUTS(diag_e) at Wed Jan 03 12:39:57 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) @@ -729,23 +731,23 @@

GLMs with temporal random effectsdplyr::glimpse(beta_post) 
## Rows: 2,000
 ## Columns: 17
-## $ `s(year_fac).1`  <dbl> 2.04273, 1.95493, 2.24298, 2.19469, 2.07041, 2.00658,…
-## $ `s(year_fac).2`  <dbl> 2.81068, 2.53752, 2.80149, 2.80148, 2.60593, 2.81318,…
-## $ `s(year_fac).3`  <dbl> 3.11034, 3.14625, 3.03162, 3.05200, 3.16094, 3.01492,…
-## $ `s(year_fac).4`  <dbl> 3.29030, 3.34690, 3.16851, 3.20425, 3.33975, 3.18636,…
-## $ `s(year_fac).5`  <dbl> 2.10375, 2.05522, 2.12469, 2.05562, 2.15567, 2.10311,…
-## $ `s(year_fac).6`  <dbl> 1.68127, 1.67443, 1.90642, 1.79536, 1.70140, 1.68694,…
-## $ `s(year_fac).7`  <dbl> 1.89145, 2.07592, 1.93246, 1.92188, 2.03806, 2.02097,…
-## $ `s(year_fac).8`  <dbl> 3.03259, 2.86407, 3.01047, 3.05813, 2.81943, 3.03458,…
-## $ `s(year_fac).9`  <dbl> 3.18867, 3.25225, 3.24096, 3.18450, 3.34753, 3.16951,…
-## $ `s(year_fac).10` <dbl> 2.53856, 2.92960, 2.79222, 2.73489, 2.51266, 2.64136,…
-## $ `s(year_fac).11` <dbl> 3.10996, 3.13094, 3.05378, 3.04545, 3.13950, 3.02744,…
-## $ `s(year_fac).12` <dbl> 3.09751, 3.18717, 3.21799, 3.28223, 3.17069, 3.25742,…
-## $ `s(year_fac).13` <dbl> 2.37870, 2.00017, 2.35753, 2.30654, 2.26299, 2.21293,…
-## $ `s(year_fac).14` <dbl> 2.52393, 2.71019, 2.67728, 2.63873, 2.53654, 2.81299,…
-## $ `s(year_fac).15` <dbl> 1.97586, 2.20523, 2.28121, 2.15232, 2.26436, 2.27818,…
-## $ `s(year_fac).16` <dbl> 1.86645, 2.16197, 2.11100, 2.21559, 2.09984, 2.15585,…
-## $ `s(year_fac).17` <dbl> 1.2658000, 1.0175300, 1.5136500, 1.2820000, 0.8126730…
+## $ `s(year_fac).1` <dbl> 1.98962, 2.15138, 1.92684, 1.95376, 2.15078, 1.99937,… +## $ `s(year_fac).2` <dbl> 2.80707, 2.53347, 2.63469, 2.70882, 2.70084, 2.67736,… +## $ `s(year_fac).3` <dbl> 3.06250, 3.13805, 3.16438, 3.05216, 3.01915, 3.12451,… +## $ `s(year_fac).4` <dbl> 3.22576, 3.25272, 3.25449, 3.18286, 3.21455, 3.22349,… +## $ `s(year_fac).5` <dbl> 2.08546, 2.21868, 2.06373, 2.02999, 2.24433, 2.05023,… +## $ `s(year_fac).6` <dbl> 1.65832, 1.90175, 1.67368, 1.54447, 1.87506, 1.67329,… +## $ `s(year_fac).7` <dbl> 2.01131, 2.10245, 1.96819, 2.00301, 2.01704, 2.07451,… +## $ `s(year_fac).8` <dbl> 3.01786, 2.83573, 2.88175, 3.07002, 2.91504, 2.89750,… +## $ `s(year_fac).9` <dbl> 3.24757, 3.23585, 3.23521, 3.29022, 3.35363, 3.35950,… +## $ `s(year_fac).10` <dbl> 2.71082, 2.78666, 2.62978, 2.66646, 2.86251, 2.64835,… +## $ `s(year_fac).11` <dbl> 2.99183, 3.06531, 3.01332, 3.13231, 3.10931, 3.04595,… +## $ `s(year_fac).12` <dbl> 3.25333, 3.19118, 3.18056, 3.21193, 3.23475, 3.24401,… +## $ `s(year_fac).13` <dbl> 2.45195, 2.04804, 2.47966, 2.47744, 2.03681, 2.47318,… +## $ `s(year_fac).14` <dbl> 2.58313, 2.68389, 2.55350, 2.60712, 2.69258, 2.60655,… +## $ `s(year_fac).15` <dbl> 2.09122, 2.22885, 2.15209, 2.21427, 2.13389, 2.24489,… +## $ `s(year_fac).16` <dbl> 2.04307, 2.09236, 2.06631, 2.06121, 2.11456, 2.07167,… +## $ `s(year_fac).17` <dbl> 0.294882, 1.309190, 1.631300, 1.572010, 1.252260, 1.8…

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

@@ -866,7 +868,7 @@

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

 hc <- hindcast(model1, type = 'link')
 range(hc$hindcasts$PP)
-
## [1] -2.87512  3.44778
+
## [1] -1.23492  3.46318

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

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

## Out of sample DRPS:
-## [1] 183.0032
+## [1] 182.79

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] 183.0032
+## [1] 182.79

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 @@ -966,12 +968,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] 8 7 10 9 7 11 6 9 14 7 ... +## ..$ PP: num [1:2000, 1:160] 12 12 8 9 6 8 16 5 7 10 ... ## .. ..- 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] 12 10 8 8 9 9 7 14 11 7 ... +## ..$ PP: num [1:2000, 1:39] 8 11 11 13 12 8 9 5 13 12 ... ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : NULL ## .. .. ..$ : chr [1:39] "ypred[161,1]" "ypred[162,1]" "ypred[163,1]" "ypred[164,1]" ... @@ -1028,33 +1030,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 1741 -## s(year_fac).1 1.10 1.40 1.70 1 2348 -## s(year_fac).2 1.80 2.00 2.20 1 2073 -## s(year_fac).3 2.20 2.40 2.60 1 2093 -## s(year_fac).4 2.30 2.50 2.70 1 1999 -## s(year_fac).5 1.20 1.40 1.70 1 2163 -## s(year_fac).6 1.00 1.30 1.50 1 2451 -## s(year_fac).7 1.10 1.40 1.70 1 2451 -## s(year_fac).8 2.10 2.30 2.50 1 2233 -## s(year_fac).9 2.70 2.90 3.00 1 2115 -## s(year_fac).10 2.00 2.20 2.40 1 2460 -## s(year_fac).11 2.30 2.40 2.60 1 2005 -## s(year_fac).12 2.50 2.70 2.80 1 2112 -## s(year_fac).13 1.40 1.60 1.80 1 2266 -## s(year_fac).14 0.69 2.00 3.30 1 1498 -## s(year_fac).15 0.59 2.00 3.20 1 1485 -## s(year_fac).16 0.63 2.00 3.30 1 1740 -## s(year_fac).17 0.65 2.00 3.30 1 1478 +## ndvi 0.33 0.39 0.46 1 1670 +## s(year_fac).1 1.10 1.40 1.70 1 2608 +## s(year_fac).2 1.80 2.00 2.20 1 2084 +## s(year_fac).3 2.20 2.40 2.60 1 2097 +## s(year_fac).4 2.30 2.50 2.70 1 2104 +## s(year_fac).5 1.20 1.40 1.60 1 2221 +## s(year_fac).6 1.00 1.30 1.50 1 2309 +## s(year_fac).7 1.10 1.40 1.70 1 2489 +## s(year_fac).8 2.10 2.30 2.50 1 1937 +## s(year_fac).9 2.70 2.90 3.00 1 2036 +## s(year_fac).10 2.00 2.20 2.40 1 2471 +## s(year_fac).11 2.30 2.40 2.60 1 2179 +## s(year_fac).12 2.50 2.70 2.80 1 2070 +## s(year_fac).13 1.40 1.60 1.80 1 2645 +## s(year_fac).14 0.66 1.90 3.20 1 1505 +## s(year_fac).15 0.57 2.00 3.20 1 1333 +## s(year_fac).16 0.56 2.00 3.40 1 1601 +## s(year_fac).17 0.53 2.00 3.30 1 1215 ## ## GAM group-level estimates: ## 2.5% 50% 97.5% Rhat n_eff -## mean(s(year_fac)) 1.60 2.0 2.40 1.01 286 -## sd(s(year_fac)) 0.41 0.6 0.97 1.01 380 +## mean(s(year_fac)) 1.60 2.0 2.30 1.01 363 +## sd(s(year_fac)) 0.42 0.6 0.98 1.01 305 ## ## Approximate significance of GAM observation smooths: -## edf Chi.sq p-value -## s(year_fac) 11 2892 <2e-16 *** +## edf Chi.sq p-value +## s(year_fac) 10.7 2788 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## @@ -1065,7 +1067,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 Wed Jan 03 9:21:08 AM 2024. +## Samples were drawn using NUTS(diag_e) at Wed Jan 03 12:41:36 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) @@ -1075,25 +1077,25 @@

Adding predictors as “fixed” eff coef:

 coef(model2)
-
##                     2.5%       50%     97.5% Rhat n_eff
-## ndvi           0.3179267 0.3884385 0.4599065    1  1741
-## s(year_fac).1  1.1381883 1.4073350 1.6664487    1  2348
-## s(year_fac).2  1.8113835 2.0010700 2.2082388    1  2073
-## s(year_fac).3  2.1891040 2.3816450 2.5858488    1  2093
-## s(year_fac).4  2.3338242 2.5116550 2.6970703    1  1999
-## s(year_fac).5  1.2016132 1.4208700 1.6612615    1  2163
-## s(year_fac).6  1.0220960 1.2824050 1.5243717    1  2451
-## s(year_fac).7  1.1214145 1.4186650 1.6875790    1  2451
-## s(year_fac).8  2.0748402 2.2711050 2.4729067    1  2233
-## s(year_fac).9  2.7158482 2.8566350 2.9914740    1  2115
-## s(year_fac).10 1.9840612 2.1945950 2.3889882    1  2460
-## s(year_fac).11 2.2720697 2.4390350 2.6083575    1  2005
-## s(year_fac).12 2.5372775 2.6940600 2.8448787    1  2112
-## s(year_fac).13 1.3821790 1.6218200 1.8466833    1  2266
-## s(year_fac).14 0.6922397 1.9795850 3.2868190    1  1498
-## s(year_fac).15 0.5921623 2.0114550 3.2356638    1  1485
-## s(year_fac).16 0.6340315 1.9646000 3.2555120    1  1740
-## s(year_fac).17 0.6489102 2.0133450 3.2578187    1  1478
+
##                     2.5%      50%     97.5% Rhat n_eff
+## ndvi           0.3269290 0.389354 0.4596232    1  1670
+## s(year_fac).1  1.1229422 1.400015 1.6541050    1  2608
+## s(year_fac).2  1.7983470 2.000555 2.1934580    1  2084
+## s(year_fac).3  2.1852790 2.381535 2.5645077    1  2097
+## s(year_fac).4  2.3394668 2.505765 2.6713190    1  2104
+## s(year_fac).5  1.1993645 1.425985 1.6406352    1  2221
+## s(year_fac).6  1.0160217 1.270055 1.5001320    1  2309
+## s(year_fac).7  1.1463117 1.415585 1.6716445    1  2489
+## s(year_fac).8  2.0796585 2.271825 2.4550455    1  1937
+## s(year_fac).9  2.7243590 2.853775 2.9882708    1  2036
+## s(year_fac).10 1.9837650 2.186395 2.3793178    1  2471
+## s(year_fac).11 2.2686295 2.436015 2.6015627    1  2179
+## s(year_fac).12 2.5308762 2.694650 2.8376307    1  2070
+## s(year_fac).13 1.3632060 1.610060 1.8466905    1  2645
+## s(year_fac).14 0.6576451 1.947180 3.2149295    1  1505
+## s(year_fac).15 0.5655013 1.964790 3.1689792    1  1333
+## s(year_fac).16 0.5624482 1.968640 3.4042000    1  1601
+## s(year_fac).17 0.5301895 2.005110 3.3217602    1  1215

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

@@ -1109,24 +1111,24 @@ 
Adding predictors as “fixed” eff
 dplyr::glimpse(beta_post)
## Rows: 2,000
 ## Columns: 18
-## $ ndvi             <dbl> 0.417499, 0.470458, 0.312997, 0.410634, 0.403887, 0.3…
-## $ `s(year_fac).1`  <dbl> 1.29733, 1.38725, 1.49245, 1.43749, 1.40625, 1.36877,…
-## $ `s(year_fac).2`  <dbl> 1.85292, 1.94953, 2.01111, 1.86587, 1.90230, 2.01867,…
-## $ `s(year_fac).3`  <dbl> 2.25025, 2.21731, 2.51190, 2.34818, 2.37044, 2.25006,…
-## $ `s(year_fac).4`  <dbl> 2.45398, 2.38964, 2.52105, 2.45516, 2.47432, 2.50800,…
-## $ `s(year_fac).5`  <dbl> 1.46950, 1.22011, 1.63476, 1.51507, 1.41279, 1.60009,…
-## $ `s(year_fac).6`  <dbl> 1.330690, 1.158300, 1.430470, 1.059190, 1.237660, 1.2…
-## $ `s(year_fac).7`  <dbl> 1.29866, 1.39150, 1.43582, 1.25721, 1.34511, 1.42439,…
-## $ `s(year_fac).8`  <dbl> 2.26303, 2.04204, 2.42008, 2.29731, 2.31482, 2.32182,…
-## $ `s(year_fac).9`  <dbl> 2.73093, 2.68393, 2.86665, 2.80806, 2.83570, 2.93193,…
-## $ `s(year_fac).10` <dbl> 2.25636, 1.93898, 2.38226, 2.11624, 2.16484, 2.14604,…
-## $ `s(year_fac).11` <dbl> 2.33500, 2.35651, 2.50488, 2.38488, 2.40927, 2.48524,…
-## $ `s(year_fac).12` <dbl> 2.71398, 2.50772, 2.73023, 2.65270, 2.68361, 2.64094,…
-## $ `s(year_fac).13` <dbl> 1.55455, 1.43178, 1.80714, 1.43513, 1.61287, 1.54666,…
-## $ `s(year_fac).14` <dbl> 0.300022, 1.241960, 1.963160, 1.077230, 1.878700, 2.0…
-## $ `s(year_fac).15` <dbl> 1.545490, 1.669830, 2.477080, 1.984560, 2.401220, 2.1…
-## $ `s(year_fac).16` <dbl> 1.920740, 1.570980, 1.930320, 1.206490, 0.971772, 3.1…
-## $ `s(year_fac).17` <dbl> 1.40739, 3.28222, 1.52013, 1.78129, 1.10294, 2.41152,…
+## $ ndvi <dbl> 0.458918, 0.415930, 0.424723, 0.371165, 0.418482, 0.4… +## $ `s(year_fac).1` <dbl> 1.35096, 1.39754, 1.50055, 1.38657, 1.36258, 1.30886,… +## $ `s(year_fac).2` <dbl> 1.69892, 2.15042, 1.81957, 1.93192, 2.03218, 1.81760,… +## $ `s(year_fac).3` <dbl> 2.17382, 2.42260, 2.20046, 2.32631, 2.27012, 2.33111,… +## $ `s(year_fac).4` <dbl> 2.36153, 2.35569, 2.44385, 2.55106, 2.45282, 2.38898,… +## $ `s(year_fac).5` <dbl> 1.33112, 1.37387, 1.51102, 1.54072, 1.46647, 1.15258,… +## $ `s(year_fac).6` <dbl> 1.386510, 1.151570, 1.449020, 1.268060, 1.364190, 1.1… +## $ `s(year_fac).7` <dbl> 1.20461, 1.41491, 1.29277, 1.62229, 1.06709, 1.46546,… +## $ `s(year_fac).8` <dbl> 2.09838, 2.25496, 2.18309, 2.31377, 2.21560, 2.18134,… +## $ `s(year_fac).9` <dbl> 2.90369, 2.82270, 2.80546, 2.85063, 2.85093, 2.82343,… +## $ `s(year_fac).10` <dbl> 1.97290, 2.27780, 2.03467, 2.18063, 2.09064, 2.11179,… +## $ `s(year_fac).11` <dbl> 2.29020, 2.47951, 2.26861, 2.35658, 2.38699, 2.40411,… +## $ `s(year_fac).12` <dbl> 2.60218, 2.59223, 2.64087, 2.79068, 2.76861, 2.57621,… +## $ `s(year_fac).13` <dbl> 1.48318, 1.66002, 1.57568, 1.71900, 1.45772, 1.46513,… +## $ `s(year_fac).14` <dbl> 1.19630, 2.63723, 2.01591, 2.19931, 2.09856, 1.43919,… +## $ `s(year_fac).15` <dbl> 1.599080, 1.969150, 1.772600, 1.853290, 2.423470, 1.1… +## $ `s(year_fac).16` <dbl> 0.721965, 2.355850, 1.942030, 1.961890, 2.169870, 2.5… +## $ `s(year_fac).17` <dbl> 2.014870, 2.480350, 2.472000, 1.531710, 2.854490, 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 @@ -1157,10 +1159,10 @@

scenario-based predictions, conditional and marginal effects, all on the outcome scale. Here we will use the plot_predictions function from marginaleffects to inspect the conditional -effect of ndvi (use ?plot_predictions for +effect of ndvi (use ?plot_predictions for guidance on how to modify these plots):

-plot_predictions(model2, 
+plot_predictions(model2, 
                  condition = "ndvi",
                  # include the observed count values
                  # as points, and show rugs for the observed
@@ -1262,27 +1264,27 @@ 
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   948
-## ndvi         0.26  0.33  0.400 1.00  1074
-## s(time).1   -2.20 -1.10  0.062 1.01   465
-## s(time).2    0.41  1.30  2.300 1.00   362
-## s(time).3   -0.52  0.48  1.500 1.01   340
-## s(time).4    1.50  2.50  3.600 1.01   326
-## s(time).5   -1.20 -0.17  0.840 1.01   361
-## s(time).6   -0.65  0.40  1.500 1.01   341
-## s(time).7   -1.50 -0.50  0.570 1.00   355
-## s(time).8    0.48  1.50  2.600 1.01   341
-## s(time).9    1.10  2.10  3.100 1.01   317
-## s(time).10  -0.41  0.57  1.600 1.01   332
-## s(time).11   0.78  1.80  2.800 1.01   320
-## s(time).12   0.62  1.50  2.500 1.01   341
-## s(time).13  -1.20 -0.29  0.710 1.01   400
-## s(time).14  -7.60 -4.30 -0.950 1.01   366
+##              2.5%   50%   97.5% Rhat n_eff
+## (Intercept)  2.00  2.10  2.2000 1.00   864
+## ndvi         0.26  0.33  0.3900 1.00   826
+## s(time).1   -2.20 -1.10 -0.0012 1.01   567
+## s(time).2    0.45  1.30  2.4000 1.01   420
+## s(time).3   -0.56  0.45  1.5000 1.01   435
+## s(time).4    1.60  2.40  3.5000 1.01   361
+## s(time).5   -1.20 -0.20  0.8800 1.01   440
+## s(time).6   -0.58  0.36  1.5000 1.01   416
+## s(time).7   -1.60 -0.54  0.5900 1.01   461
+## s(time).8    0.53  1.50  2.6000 1.01   361
+## s(time).9    1.20  2.00  3.1000 1.01   409
+## s(time).10  -0.36  0.54  1.6000 1.01   411
+## s(time).11   0.80  1.70  2.8000 1.01   402
+## s(time).12   0.66  1.50  2.5000 1.01   404
+## s(time).13  -1.20 -0.32  0.6800 1.01   535
+## s(time).14  -7.60 -4.10 -1.0000 1.01   544
 ## 
 ## Approximate significance of GAM observation smooths:
 ##          edf Chi.sq p-value    
-## s(time) 8.55    684  <2e-16 ***
+## s(time) 8.91    673  <2e-16 ***
 ## ---
 ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
 ## 
@@ -1293,7 +1295,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 Wed Jan 03 9:22:07 AM 2024.
+## Samples were drawn using NUTS(diag_e) at Wed Jan 03 12:42:35 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)
@@ -1420,7 +1422,7 @@

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

## Out of sample DRPS:
-## [1] 287.4924
+## [1] 286.2666

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 @@ -1502,40 +1504,40 @@

Latent dynamics in mvgam## ## ## GAM coefficient (beta) estimates: -## 2.5% 50% 97.5% Rhat n_eff -## (Intercept) 0.970 1.9000 2.800 1.19 22 -## s(ndvi).1 -0.160 -0.0084 0.072 1.01 573 -## s(ndvi).2 -0.160 0.0150 0.270 1.01 445 -## s(ndvi).3 -0.049 -0.0018 0.045 1.00 965 -## s(ndvi).4 -0.270 0.1100 1.200 1.01 318 -## s(ndvi).5 -0.066 0.1500 0.360 1.00 602 +## 2.5% 50% 97.5% Rhat n_eff +## (Intercept) 0.790 2.000 2.600 1.08 44 +## s(ndvi).1 -0.150 -0.010 0.077 1.00 674 +## s(ndvi).2 -0.130 0.016 0.310 1.01 410 +## s(ndvi).3 -0.048 -0.001 0.042 1.00 927 +## s(ndvi).4 -0.240 0.120 1.200 1.01 272 +## s(ndvi).5 -0.062 0.150 0.370 1.00 438 ## ## Approximate significance of GAM observation smooths: ## edf Chi.sq p-value -## s(ndvi) 1.64 78.9 0.07 . +## s(ndvi) 1.13 87.5 0.079 . ## --- ## 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.92 1.02 315 -## sigma[1] 0.68 0.81 0.97 1.02 417 +## ar1[1] 0.69 0.81 0.93 1.02 155 +## sigma[1] 0.67 0.80 0.95 1.01 386 ## ## Stan MCMC diagnostics: ## n_eff / iter looks reasonable for all parameters -## Rhats above 1.05 found for 124 parameters +## Rhats above 1.05 found for 64 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 Wed Jan 03 9:23:33 AM 2024. +## Samples were drawn using NUTS(diag_e) at Wed Jan 03 12:44:01 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)

View conditional smooths for the ndvi effect:

-plot_predictions(model4, 
+plot_predictions(model4, 
                  condition = "ndvi",
                  points = 0.5, rug = TRUE)

@@ -1544,8 +1546,8 @@

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

-
## Out of sample DRPS:
-## [1] 153.0846
+
## Out of sample CRPS:
+## [1] 150.1441

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

 plot(model4, type = 'trend', newdata = data_test)
@@ -1560,7 +1562,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 -562.1      66.8
+## model3 -558.4 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.

@@ -1575,7 +1577,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] -134.4078
+
## [1] -136.1226

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 4643415f..e5ae5807 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 7e942103..99810ad4 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 6820b028..19815a75 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 647e7a1a..79ec8fb4 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 6850ec54..c72a713e 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 3d6034bf..0e28c3fd 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 8648b92d..b49d4b5c 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 1fc049a9..830828d7 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 eb1f24e6..ef3af78d 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 d135d82b..26013ece 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 6c4b58e3..29738616 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 eb1f24e6..ef3af78d 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 a2665159..594df361 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 0994e12b..c389e369 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 9108dd28..c69feb58 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 3cfae64f..4f6bc71a 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 17c5cb99..fdf53631 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 a0dcec43..6443254a 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 cfc2547d..73e30b96 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 192dec4a..f77e235c 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 52dd51d1..4254f45d 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 5d7cd8cb..cd58fdb0 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 fda2ea4b..30ed8476 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 ce3d84fe..fdd8f4ea 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 7849c09f..90247132 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/src/mvgam.dll b/src/mvgam.dll index 4989a96e..b5b832a8 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 252fccb5..f2edf497 100644 --- a/vignettes/mvgam_overview.Rmd +++ b/vignettes/mvgam_overview.Rmd @@ -126,10 +126,10 @@ Sean J. Taylor and Benjamin Letham. "[Forecasting at scale.](https://www.tandfon `mvgam` supports an observation model regression formula, built off the `mvgcv` package, as well as an optional process model regression formula. The formulae supplied to \code{\link{mvgam}} are exactly like those supplied to `glm()` except that smooth terms, `s()`, `te()`, `ti()` and `t2()`, time-varying effects using `dynamic()`, monotonically increasing (using `s(x, bs = 'moi')`) or decreasing splines (using `s(x, bs = 'mod')`; see `?smooth.construct.moi.smooth.spec` for details), as well as Gaussian Process functions using `gp()`, can be added to the right hand side (and `.` is not supported in `mvgam` formulae). See `?mvgam_formulae` for more guidance. -For setting up State-Space models, the ptional process model formula can be used (see \href{https://nicholasjclark.github.io/mvgam/articles/trend_formulas.html}{the State-Space model vignette} and \href{https://nicholasjclark.github.io/mvgam/articles/trend_formulas.html}{the shared latent states vignette} for guidance on using trend formulae). +For setting up State-Space models, the optional process model formula can be used (see [the State-Space model vignette](https://nicholasjclark.github.io/mvgam/articles/trend_formulas.html) and [the shared latent states vignette](https://nicholasjclark.github.io/mvgam/articles/trend_formulas.html) for guidance on using trend formulae). ## Example time series data -The 'portal_data' object contains time series of rodent captures from the Portal Project, [a long-term monitoring study based near the town of Portal, Arizona](https://portal.weecology.org/){target="_blank"}. Researchers have been operating a standardized set of baited traps within 24 experimental plots at this site since the 1970's. Sampling follows the lunar monthly cycle, with observations occurring on average about 28 days apart. However, missing observations do occur due to difficulties accessing the site (weather events, COVID disruptions etc..). You can read about the full sampling protocol [in this preprint by Ernest et al on the Biorxiv](https://www.biorxiv.org/content/10.1101/332783v3.full){target="_blank"}. +The 'portal_data' object contains time series of rodent captures from the Portal Project, [a long-term monitoring study based near the town of Portal, Arizona](https://portal.weecology.org/){target="_blank"}. Researchers have been operating a standardized set of baited traps within 24 experimental plots at this site since the 1970's. Sampling follows the lunar monthly cycle, with observations occurring on average about 28 days apart. However, missing observations do occur due to difficulties accessing the site (weather events, COVID disruptions etc...). You can read about the full sampling protocol [in this preprint by Ernest et al on the Biorxiv](https://www.biorxiv.org/content/10.1101/332783v3.full){target="_blank"}. ```{r Access time series data} data("portal_data") ```