diff --git a/README.Rmd b/README.Rmd index f95be11c..19d55533 100644 --- a/README.Rmd +++ b/README.Rmd @@ -18,6 +18,8 @@ knitr::opts_chunk$set( ) ``` +brms Logo[Stan Logo](https://mc-stan.org/) + *mvgam* ================ diff --git a/README.md b/README.md index 986adb5e..141647e4 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,8 @@ +brms Logo[Stan Logo](https://mc-stan.org/) + # *mvgam* The goal of `mvgam` is to fit Bayesian Dynamic Generalized Additive @@ -317,29 +319,29 @@ summary(lynx_mvgam) #> #> #> GAM coefficient (beta) estimates: -#> 2.5% 50% 97.5% Rhat n_eff -#> (Intercept) 6.100 6.60000 7.0000 1.00 659 -#> s(season).1 -0.630 0.02900 0.7500 1.00 1162 -#> s(season).2 -0.250 0.81000 1.8000 1.01 422 -#> s(season).3 -0.039 1.20000 2.5000 1.01 371 -#> s(season).4 -0.560 0.43000 1.4000 1.01 902 -#> s(season).5 -1.200 -0.14000 1.0000 1.01 565 -#> s(season).6 -1.100 0.00089 1.1000 1.01 703 -#> s(season).7 -0.810 0.37000 1.5000 1.00 873 -#> s(season).8 -0.980 0.26000 1.8000 1.01 366 -#> s(season).9 -1.100 -0.29000 0.6600 1.00 474 -#> s(season).10 -1.400 -0.67000 0.0049 1.00 673 +#> 2.5% 50% 97.5% Rhat n_eff +#> (Intercept) 6.10 6.6000 7.0000 1.00 771 +#> s(season).1 -0.59 0.0480 0.7200 1.00 973 +#> s(season).2 -0.25 0.8000 1.8000 1.01 389 +#> s(season).3 -0.12 1.2000 2.5000 1.01 321 +#> s(season).4 -0.55 0.4100 1.3000 1.00 892 +#> s(season).5 -1.20 -0.1000 0.9300 1.00 479 +#> s(season).6 -1.00 -0.0089 1.1000 1.00 577 +#> s(season).7 -0.72 0.3400 1.4000 1.01 659 +#> s(season).8 -0.95 0.2200 1.8000 1.01 328 +#> s(season).9 -1.10 -0.3100 0.6700 1.00 398 +#> s(season).10 -1.30 -0.6700 -0.0054 1.01 595 #> #> Approximate significance of GAM observation smooths: #> edf Chi.sq p-value -#> s(season) 3.71 18666 0.28 +#> s(season) 4.44 17999 0.27 #> #> Latent trend AR parameter estimates: #> 2.5% 50% 97.5% Rhat n_eff -#> ar1[1] 0.75 1.10 1.400 1.00 614 -#> ar2[1] -0.84 -0.41 0.042 1.00 1688 -#> ar3[1] -0.48 -0.12 0.300 1.01 571 -#> sigma[1] 0.40 0.50 0.640 1.00 868 +#> ar1[1] 0.73 1.10 1.400 1.00 609 +#> ar2[1] -0.82 -0.41 0.043 1.00 1562 +#> ar3[1] -0.46 -0.12 0.300 1.01 510 +#> sigma[1] 0.40 0.50 0.630 1.00 1051 #> #> Stan MCMC diagnostics: #> n_eff / iter looks reasonable for all parameters @@ -348,7 +350,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 Fri Mar 01 9:45:36 AM 2024. +#> Samples were drawn using NUTS(diag_e) at Tue Mar 12 1:50: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) @@ -453,7 +455,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 +#> Warning: package 'ggplot2' was built under R version 4.3.3 plot_predictions(lynx_mvgam, condition = 'season', points = 0.5) + theme_classic() ``` @@ -470,7 +472,7 @@ plot(lynx_mvgam, type = 'forecast', newdata = lynx_test) #> Out of sample CRPS: - #> [1] 2901.021 + #> [1] 2852.181 And the estimated latent trend component, again using the more flexible `plot_mvgam_...()` option to show first derivatives of the estimated @@ -626,7 +628,7 @@ summary(mod, include_betas = FALSE) #> 0 of 2000 iterations saturated the maximum tree depth of 12 (0%) #> E-FMI indicated no pathological behavior #> -#> Samples were drawn using NUTS(diag_e) at Fri Mar 01 9:46:54 AM 2024. +#> Samples were drawn using NUTS(diag_e) at Tue Mar 12 1:51:32 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) diff --git a/docs/articles/data_in_mvgam.html b/docs/articles/data_in_mvgam.html index 6282b2f6..00004786 100644 --- a/docs/articles/data_in_mvgam.html +++ b/docs/articles/data_in_mvgam.html @@ -78,7 +78,7 @@

Nicholas J Clark

-

2024-03-01

+

2024-03-12

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

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

series as a factor variable @@ -174,24 +174,24 @@

A single outcome variable## ## Coefficients: ## Estimate Std. Error z value Pr(>|z|) -## (Intercept) 1.02109 0.22851 4.468 7.88e-06 *** -## seriesseries_2 -1.27226 0.30950 -4.111 3.95e-05 *** -## seriesseries_3 -1.13938 0.29265 -3.893 9.89e-05 *** -## seriesseries_4 -1.28613 0.30985 -4.151 3.31e-05 *** -## time 0.01919 0.01929 0.995 0.32 +## (Intercept) -0.50961 0.36743 -1.387 0.165449 +## seriesseries_2 0.88571 0.35268 2.511 0.012025 * +## seriesseries_3 0.53206 0.37538 1.417 0.156370 +## seriesseries_4 1.22735 0.34590 3.548 0.000388 *** +## time 0.02585 0.01963 1.317 0.187870 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## (Dispersion parameter for poisson family taken to be 1) ## -## Null deviance: 120.57 on 57 degrees of freedom -## Residual deviance: 85.19 on 53 degrees of freedom -## (14 observations deleted due to missingness) -## AIC: 196.15 +## Null deviance: 191.22 on 58 degrees of freedom +## Residual deviance: 174.01 on 54 degrees of freedom +## (13 observations deleted due to missingness) +## AIC: 271.15 ## -## Number of Fisher Scoring iterations: 5 +## Number of Fisher Scoring iterations: 6
-summary(gam(y ~ series + s(time, by = series),
+summary(gam(y ~ series + s(time, by = series),
             data = simdat$data_train,
             family = poisson()))
## 
@@ -202,25 +202,23 @@ 

A single outcome variable## y ~ series + s(time, by = series) ## ## Parametric coefficients: -## Estimate Std. Error z value Pr(>|z|) -## (Intercept) 1.1032 0.1497 7.370 1.71e-13 *** -## seriesseries_2 -1.1967 0.3203 -3.737 0.000187 *** -## seriesseries_3 -1.1768 0.3324 -3.540 0.000400 *** -## seriesseries_4 -1.6351 0.4268 -3.831 0.000128 *** -## --- -## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 +## Estimate Std. Error z value Pr(>|z|) +## (Intercept) -3.055 2.893 -1.056 0.291 +## seriesseries_2 1.504 3.671 0.410 0.682 +## seriesseries_3 2.907 2.924 0.994 0.320 +## seriesseries_4 3.277 2.917 1.124 0.261 ## ## Approximate significance of smooth terms: -## edf Ref.df Chi.sq p-value -## s(time):seriesseries_1 3.935 4.864 10.118 0.0592 . -## s(time):seriesseries_2 1.000 1.000 0.886 0.3466 -## s(time):seriesseries_3 3.048 3.773 7.395 0.1220 -## s(time):seriesseries_4 3.579 4.391 11.105 0.0358 * +## edf Ref.df Chi.sq p-value +## s(time):seriesseries_1 3.539 4.329 6.824 0.226478 +## s(time):seriesseries_2 7.600 8.073 17.087 0.045922 * +## s(time):seriesseries_3 7.828 8.555 18.214 0.018101 * +## s(time):seriesseries_4 3.698 4.605 24.162 0.000195 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## -## R-sq.(adj) = 0.5 Deviance explained = 59.5% -## UBRE = 0.37832 Scale est. = 1 n = 58

+## R-sq.(adj) = 0.453 Deviance explained = 72.3% +## UBRE = 0.80158 Scale est. = 1 n = 59

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

A single outcome variable time = 1:10) gauss_dat

##       outcome  series time
-## 1  -0.5999200 series1    1
-## 2  -1.3346981 series1    2
-## 3   1.1674495 series1    3
-## 4   0.2864362 series1    4
-## 5   0.4372319 series1    5
-## 6  -0.4989263 series1    6
-## 7  -0.8585886 series1    7
-## 8   0.0579170 series1    8
-## 9   1.0562994 series1    9
-## 10  1.0117023 series1   10
+## 1 0.1081561 series1 1 +## 2 -0.5962646 series1 2 +## 3 -0.6245927 series1 3 +## 4 -0.6303395 series1 4 +## 5 -0.9292285 series1 5 +## 6 -0.7462733 series1 6 +## 7 1.3029532 series1 7 +## 8 1.7143050 series1 8 +## 9 0.2084016 series1 9 +## 10 -0.8205039 series1 10

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

-gam(outcome ~ time,
-    family = betar(),
+gam(outcome ~ time,
+    family = betar(),
     data = gauss_dat)
## Warning in family$saturated.ll(y, prior.weights, theta): saturated likelihood
 ## may be inaccurate
## 
-## Family: Beta regression(0.104) 
+## Family: Beta regression(0.095) 
 ## Link function: logit 
 ## 
 ## Formula:
 ## outcome ~ time
 ## Total model degrees of freedom 2 
 ## 
-## REML score: -177.7925
+## REML score: -206.4319

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

 mvgam(outcome ~ time,
-      family = betar(),
+      family = betar(),
       data = gauss_dat)
## Error: Values <= 0 not allowed for beta responses

Please see ?mvgam_families for more information on the @@ -348,14 +346,14 @@

Checking data with get_mvgam_ outcome = rnorm(8)) bad_times
##   time   series    outcome
-## 1    1 series_1 -1.8522467
-## 2    3 series_1 -0.5934998
-## 3    5 series_1 -0.7857181
-## 4    7 series_1 -0.5972708
-## 5    9 series_1 -0.8471251
-## 6   11 series_1  1.3899552
-## 7   13 series_1 -1.1900301
-## 8   15 series_1  1.6603869
+## 1 1 series_1 -1.0151796 +## 2 3 series_1 -0.4903414 +## 3 5 series_1 1.3391872 +## 4 7 series_1 1.2580140 +## 5 9 series_1 -1.0748462 +## 6 11 series_1 0.8043384 +## 7 13 series_1 -0.9175521 +## 8 15 series_1 -0.8733635

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

@@ -378,21 +376,21 @@ 

Checking data with get_mvgam_
 good_times
##    time   series    outcome
-## 1     1 series_1 -1.8522467
+## 1     1 series_1 -1.0151796
 ## 2     2 series_1         NA
-## 3     3 series_1 -0.5934998
+## 3     3 series_1 -0.4903414
 ## 4     4 series_1         NA
-## 5     5 series_1 -0.7857181
+## 5     5 series_1  1.3391872
 ## 6     6 series_1         NA
-## 7     7 series_1 -0.5972708
+## 7     7 series_1  1.2580140
 ## 8     8 series_1         NA
-## 9     9 series_1 -0.8471251
+## 9     9 series_1 -1.0748462
 ## 10   10 series_1         NA
-## 11   11 series_1  1.3899552
+## 11   11 series_1  0.8043384
 ## 12   12 series_1         NA
-## 13   13 series_1 -1.1900301
+## 13   13 series_1 -0.9175521
 ## 14   14 series_1         NA
-## 15   15 series_1  1.6603869
+## 15 15 series_1 -0.8733635

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

@@ -402,9 +400,9 @@ 

Checking data with get_mvgam_
##                             param_name param_length           param_info
 ## 1                          (Intercept)            1          (Intercept)
 ## 2 vector<lower=0>[n_series] sigma_obs;            1 observation error sd
-##                                    prior                   example_change
-## 1 (Intercept) ~ student_t(3, -0.7, 2.5);      (Intercept) ~ normal(0, 1);
-## 2      sigma_obs ~ student_t(3, 0, 2.5); sigma_obs ~ normal(-0.24, 0.14);
+##                                    prior                  example_change
+## 1 (Intercept) ~ student_t(3, -0.7, 2.5);     (Intercept) ~ normal(0, 1);
+## 2      sigma_obs ~ student_t(3, 0, 2.5); sigma_obs ~ normal(0.56, 0.24);
 ##   new_lowerbound new_upperbound
 ## 1             NA             NA
 ## 2             NA             NA
@@ -452,9 +450,9 @@

Checking data with get_mvgam_
##                             param_name param_length           param_info
 ## 1                          (Intercept)            1          (Intercept)
 ## 2 vector<lower=0>[n_series] sigma_obs;            1 observation error sd
-##                                   prior                  example_change
-## 1 (Intercept) ~ student_t(3, 0.2, 2.5);     (Intercept) ~ normal(0, 1);
-## 2     sigma_obs ~ student_t(3, 0, 2.5); sigma_obs ~ normal(0.79, 0.33);
+##                                 prior                  example_change
+## 1 (Intercept) ~ student_t(3, 0, 2.5);     (Intercept) ~ normal(0, 1);
+## 2   sigma_obs ~ student_t(3, 0, 2.5); sigma_obs ~ normal(0.45, 0.98);
 ##   new_lowerbound new_upperbound
 ## 1             NA             NA
 ## 2             NA             NA
@@ -462,7 +460,7 @@

Checking data with get_mvgam_

Covariates with no NAs

Covariates can be used in models just as you would when using -mgcv (see ?formula.gam for details of the +mgcv (see ?formula.gam for details of the formula syntax). But although the outcome variable can have NAs, covariates cannot. Most regression software will silently drop any raws in the model matrix that have NAs, @@ -478,22 +476,22 @@

Covariates with no NAs= 1:10) miss_dat

##       outcome         cov  series time
-## 1  -0.6560146          NA series1    1
-## 2  -0.6507284  1.35999475 series1    2
-## 3   0.2211108 -0.63391398 series1    3
-## 4   2.0087873 -0.17705523 series1    4
-## 5   0.2995543 -1.05899110 series1    5
-## 6  -0.3439845 -1.12735849 series1    6
-## 7   1.4274514  1.71495588 series1    7
-## 8   0.3591625 -0.01438412 series1    8
-## 9   0.6607156 -0.34874780 series1    9
-## 10  0.7234716  0.01834838 series1   10
+## 1 0.1793715 NA series1 1 +## 2 -0.2433165 1.86890517 series1 2 +## 3 -1.4273455 -0.05409832 series1 3 +## 4 -1.4002809 -0.65979788 series1 4 +## 5 -1.3791300 0.33206293 series1 5 +## 6 -0.6647606 -0.54721147 series1 6 +## 7 -1.3586209 -0.57260756 series1 7 +## 8 -0.3662525 0.76195325 series1 8 +## 9 -1.5266103 -0.75760094 series1 9 +## 10 1.8979686 0.43271133 series1 10
 get_mvgam_priors(outcome ~ cov,
                  data = miss_dat,
                  family = gaussian())
## Error: Missing values found in data predictors:
-##  Error in na.fail.default(structure(list(outcome = c(-0.656014605060502, : missing values in object
+## Error in na.fail.default(structure(list(outcome = c(0.179371484078396, : missing values in object

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

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

## Error: Missing values found in data predictors:
-##  Error in na.fail.default(structure(list(outcome = c(-0.0751319618758875, : missing values in object
+## Error in na.fail.default(structure(list(outcome = c(-1.79094121688544, : missing values in object
@@ -687,8 +685,8 @@

Example with NEON tick data
-testmod <- mvgam(y ~ s(epiWeek, by = series, bs = 'cc') +
-                   s(series, bs = 're'),
+testmod <- mvgam(y ~ s(epiWeek, by = series, bs = 'cc') +
+                   s(series, bs = 're'),
                  trend_model = 'AR1',
                  data = model_dat,
                  backend = 'cmdstanr',
@@ -714,12 +712,12 @@ 

Example with NEON tick data## $ S6 : num [1:8, 1:8] 1.037 -0.416 0.419 0.117 0.188 ... ## $ S7 : num [1:8, 1:8] 1.037 -0.416 0.419 0.117 0.188 ... ## $ S8 : num [1:8, 1:8] 1.037 -0.416 0.419 0.117 0.188 ... -## $ p_coefs : Named num 0.7 +## $ p_coefs : Named num 0.726 ## ..- attr(*, "names")= chr "(Intercept)" -## $ p_taus : num 1316 +## $ p_taus : num 314 ## $ ytimes : int [1:416, 1:8] 1 9 17 25 33 41 49 57 65 73 ... ## $ n_series : int 8 -## $ sp : Named num [1:9] 1.86 6.34 7.98 23.32 11.21 ... +## $ sp : Named num [1:9] 1.14 1.55 1.08 4.68 3 ... ## ..- attr(*, "names")= chr [1:9] "s(epiWeek):seriesBLAN_005" "s(epiWeek):seriesBLAN_012" "s(epiWeek):seriesSCBI_002" "s(epiWeek):seriesSCBI_013" ... ## $ y_observed : num [1:416, 1:8] 0 0 0 0 0 0 0 0 0 0 ... ## $ total_obs : int 3328 diff --git a/docs/articles/data_in_mvgam_files/figure-html/unnamed-chunk-23-1.png b/docs/articles/data_in_mvgam_files/figure-html/unnamed-chunk-23-1.png index 4dbe1238..4ddbf75a 100644 Binary files a/docs/articles/data_in_mvgam_files/figure-html/unnamed-chunk-23-1.png and b/docs/articles/data_in_mvgam_files/figure-html/unnamed-chunk-23-1.png differ diff --git a/docs/articles/data_in_mvgam_files/figure-html/unnamed-chunk-24-1.png b/docs/articles/data_in_mvgam_files/figure-html/unnamed-chunk-24-1.png index 4a47ca94..292d3229 100644 Binary files a/docs/articles/data_in_mvgam_files/figure-html/unnamed-chunk-24-1.png and b/docs/articles/data_in_mvgam_files/figure-html/unnamed-chunk-24-1.png differ diff --git a/docs/articles/data_in_mvgam_files/figure-html/unnamed-chunk-25-1.png b/docs/articles/data_in_mvgam_files/figure-html/unnamed-chunk-25-1.png index af7f693a..a00e2f27 100644 Binary files a/docs/articles/data_in_mvgam_files/figure-html/unnamed-chunk-25-1.png and b/docs/articles/data_in_mvgam_files/figure-html/unnamed-chunk-25-1.png differ diff --git a/docs/articles/forecast_evaluation.html b/docs/articles/forecast_evaluation.html index 590d6ecd..3912f228 100644 --- a/docs/articles/forecast_evaluation.html +++ b/docs/articles/forecast_evaluation.html @@ -78,7 +78,7 @@

Nicholas J Clark

-

2024-03-01

+

2024-03-12

Source:
vignettes/forecast_evaluation.Rmd
forecast_evaluation.Rmd
@@ -178,8 +178,8 @@

Modelling dynamics with splines
-mod1 <- mvgam(y ~ s(season, bs = 'cc', k = 8) + 
-                s(time, by = series, bs = 'cr', k = 20),
+mod1 <- mvgam(y ~ s(season, bs = 'cc', k = 8) + 
+                s(time, by = series, bs = 'cr', k = 20),
               knots = list(season = c(0.5, 12.5)),
               trend_model = 'None',
               data = simdat$data_train)
@@ -231,7 +231,7 @@

Modelling dynamics with splines## 0 of 2000 iterations saturated the maximum tree depth of 12 (0%) ## E-FMI indicated no pathological behavior ## -## Samples were drawn using NUTS(diag_e) at Fri Mar 01 10:09:14 AM 2024. +## Samples were drawn using NUTS(diag_e) at Tue Mar 12 1:56:47 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)

@@ -247,12 +247,12 @@

Modelling dynamics with GPsgp() function from brms, which can fit Hilbert space approximate GPs. See ?brms::gp for more details.

-mod2 <- mvgam(y ~ s(season, bs = 'cc', k = 8) + 
-                gp(time, by = series, c = 5/4, k = 20,
+mod2 <- mvgam(y ~ s(season, bs = 'cc', k = 8) + 
+                gp(time, by = series, c = 5/4, k = 20,
                    scale = FALSE),
               knots = list(season = c(0.5, 12.5)),
               trend_model = 'None',
@@ -313,7 +313,7 @@ 

Modelling dynamics with GPs## 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 Mar 01 10:10:16 AM 2024. +## Samples were drawn using NUTS(diag_e) at Tue Mar 12 1:57:51 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)

@@ -322,12 +322,12 @@

Modelling dynamics with GPs\(\alpha\)) parameters:

-mcmc_plot(mod2, variable = c('alpha_gp'), regex = TRUE, type = 'areas')
+mcmc_plot(mod2, variable = c('alpha_gp'), regex = TRUE, type = 'areas')

And now the length scale (\(\rho\)) parameters:

-mcmc_plot(mod2, variable = c('rho_gp'), regex = TRUE, type = 'areas')
+mcmc_plot(mod2, variable = c('rho_gp'), regex = TRUE, type = 'areas')

We can also plot the nonlinear effects as before:

@@ -337,7 +337,7 @@ 

Modelling dynamics with GPs
 require('ggplot2')
-plot_predictions(mod2, 
+plot_predictions(mod2, 
                  condition = c('time', 'series', 'series'),
                  type = 'link') +
   theme(legend.position = 'none')

@@ -460,8 +460,8 @@

Forecasting with newdata mod2 but include the testing data for automatic forecasts:

-mod2 <- mvgam(y ~ s(season, bs = 'cc', k = 8) + 
-                gp(time, by = series, c = 5/4, k = 20,
+mod2 <- mvgam(y ~ s(season, bs = 'cc', k = 8) + 
+                gp(time, by = series, c = 5/4, k = 20,
                    scale = FALSE),
               knots = list(season = c(0.5, 12.5)),
               trend_model = 'None',
diff --git a/docs/articles/mvgam_overview.html b/docs/articles/mvgam_overview.html
index ac0ca32f..0a44a5dd 100644
--- a/docs/articles/mvgam_overview.html
+++ b/docs/articles/mvgam_overview.html
@@ -78,7 +78,7 @@
                         

Nicholas J Clark

-

2024-03-01

+

2024-03-12

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

Supported observation families LogNormal (identity link) -lognormal() +lognormal() Positive real values in \([0, \infty)\) \(\sigma\) @@ -169,7 +169,7 @@

Supported observation families Beta (logit link) -betar() +betar() Real values (proportional) in \([0,1]\) \(\phi\) @@ -183,7 +183,7 @@

Supported observation families Negative Binomial2 (log link) -nb() +nb() Non-negative integers in \((0,1,2,...)\) \(\phi\) @@ -369,12 +369,12 @@

Regression formulaeglm() except that smooth terms, -s(), te(), ti() and -t2(), time-varying effects using dynamic(), +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 +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.

@@ -446,14 +446,14 @@

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  1      1    1 series_2    1
-## 3  4      1    1 series_3    1
-## 4  2      1    1 series_4    1
-## 5  0      2    1 series_1    2
+## 1  0      1    1 series_1    1
+## 2  0      1    1 series_2    1
+## 3  0      1    1 series_3    1
+## 4  0      1    1 series_4    1
+## 5  1      2    1 series_1    2
 ## 6  0      2    1 series_2    2
-## 7  1      2    1 series_3    2
-## 8  0      2    1 series_4    2
+## 7  0      2    1 series_3    2
+## 8  1      2    1 series_4    2
 ## 9  2      3    1 series_1    3
 ## 10 0      3    1 series_2    3
 ## 11 0      3    1 series_3    3
@@ -578,8 +578,8 @@ 

GLMs with temporal random effectsyear to a factor so that we can use a random effect -basis in mvgam. See ?smooth.terms and -?smooth.construct.re.smooth.spec for details about the +basis in mvgam. See ?smooth.terms and +?smooth.construct.re.smooth.spec for details about the re basis construction that is used by both mvgam and mgcv

@@ -606,8 +606,8 @@ 

GLMs with temporal random effects## [11] "2014" "2015" "2016" "2017" "2018" "2019" "2020"

We are now ready for our first mvgam model. The syntax will be familiar to users who have previously built models with -mgcv. But for a refresher, see ?formula.gam -and the examples in ?gam. Random effects can be specified +mgcv. But for a refresher, see ?formula.gam +and the examples in ?gam. Random effects can be specified using the s wrapper with the re basis. Note that we can also suppress the primary intercept using the usual R formula syntax - 1. mvgam has a @@ -626,7 +626,7 @@

GLMs with temporal random effects
-model1 <- mvgam(count ~ s(year_fac, bs = 're') - 1,
+model1 <- mvgam(count ~ s(year_fac, bs = 're') - 1,
                 family = poisson(),
                 data = model_data)

The model can be described mathematically for each timepoint \(t\) as follows: \[\begin{align*} @@ -643,15 +643,15 @@

GLMs with temporal random effects\((\mu_{year})\) and \((\sigma_{year})\) can be viewed using the following code:

-get_mvgam_priors(count ~ s(year_fac, bs = 're') - 1,
+get_mvgam_priors(count ~ s(year_fac, bs = 're') - 1,
                  family = poisson(),
                  data = model_data)
##                      param_name param_length           param_info
 ## 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.96, 0.78);
-## 2 sigma_raw ~ student_t(3, 0, 2.5); sigma_raw ~ exponential(0.58);
+## 1            mu_raw ~ std_normal();  mu_raw ~ normal(-0.32, 0.51);
+## 2 sigma_raw ~ student_t(3, 0, 2.5); sigma_raw ~ exponential(0.81);
 ##   new_lowerbound new_upperbound
 ## 1             NA             NA
 ## 2             NA             NA
@@ -688,32 +688,32 @@

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 2213 -## s(year_fac).2 2.50 2.7 2.8 1 2846 -## s(year_fac).3 3.00 3.1 3.2 1 3350 -## s(year_fac).4 3.10 3.3 3.4 1 2494 -## s(year_fac).5 1.90 2.1 2.3 1 3043 -## s(year_fac).6 1.50 1.8 2.0 1 2577 -## s(year_fac).7 1.80 2.0 2.3 1 2386 -## s(year_fac).8 2.80 3.0 3.1 1 3166 -## s(year_fac).9 3.10 3.3 3.4 1 2760 -## s(year_fac).10 2.60 2.8 2.9 1 2872 -## s(year_fac).11 2.90 3.1 3.2 1 2973 -## s(year_fac).12 3.10 3.2 3.3 1 2742 -## s(year_fac).13 2.00 2.2 2.4 1 2713 -## s(year_fac).14 2.50 2.6 2.8 1 2425 -## s(year_fac).15 1.90 2.2 2.4 1 2522 -## s(year_fac).16 1.90 2.1 2.3 1 2808 -## s(year_fac).17 -0.26 1.1 1.9 1 455 +## s(year_fac).1 1.80 2.1 2.3 1 3394 +## s(year_fac).2 2.50 2.7 2.8 1 2659 +## s(year_fac).3 3.00 3.1 3.2 1 3070 +## s(year_fac).4 3.10 3.3 3.4 1 2710 +## s(year_fac).5 1.90 2.1 2.3 1 3262 +## s(year_fac).6 1.50 1.8 2.0 1 3151 +## s(year_fac).7 1.80 2.0 2.3 1 3153 +## s(year_fac).8 2.80 3.0 3.1 1 3226 +## s(year_fac).9 3.10 3.2 3.4 1 2604 +## s(year_fac).10 2.60 2.8 2.9 1 3122 +## s(year_fac).11 2.90 3.1 3.2 1 3207 +## s(year_fac).12 3.10 3.2 3.3 1 2495 +## s(year_fac).13 2.00 2.2 2.5 1 2844 +## s(year_fac).14 2.50 2.6 2.8 1 2654 +## s(year_fac).15 1.90 2.2 2.4 1 2648 +## s(year_fac).16 1.90 2.1 2.3 1 3048 +## s(year_fac).17 -0.22 1.1 1.9 1 608 ## ## GAM group-level estimates: ## 2.5% 50% 97.5% Rhat n_eff -## mean(s(year_fac)) 2.00 2.40 2.7 1.02 161 -## sd(s(year_fac)) 0.48 0.68 1.2 1.01 221 +## mean(s(year_fac)) 2.10 2.40 2.7 1.02 237 +## sd(s(year_fac)) 0.44 0.67 1.0 1.01 284 ## ## Approximate significance of GAM observation smooths: -## edf Chi.sq p-value -## s(year_fac) 13.7 24841 <2e-16 *** +## edf Chi.sq p-value +## s(year_fac) 14 24250 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## @@ -724,7 +724,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 Mar 01 10:11:47 AM 2024. +## Samples were drawn using NUTS(diag_e) at Tue Mar 12 1:59:29 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)

@@ -741,23 +741,23 @@

GLMs with temporal random effectsdplyr::glimpse(beta_post)
## Rows: 2,000
 ## Columns: 17
-## $ `s(year_fac).1`  <dbl> 1.76044, 2.25401, 1.77490, 2.15955, 1.86358, 2.06236,…
-## $ `s(year_fac).2`  <dbl> 2.81054, 2.60305, 2.76941, 2.61998, 2.67258, 2.72234,…
-## $ `s(year_fac).3`  <dbl> 3.15449, 3.11106, 3.14198, 3.01846, 3.17096, 3.07397,…
-## $ `s(year_fac).4`  <dbl> 3.26832, 3.30034, 3.24835, 3.24246, 3.30360, 3.16272,…
-## $ `s(year_fac).5`  <dbl> 2.12277, 2.08530, 2.19325, 2.08468, 2.25269, 2.05183,…
-## $ `s(year_fac).6`  <dbl> 1.85119, 1.71227, 1.74788, 1.73515, 1.70462, 1.82687,…
-## $ `s(year_fac).7`  <dbl> 2.13533, 2.06434, 1.89028, 2.03741, 1.99278, 2.24407,…
-## $ `s(year_fac).8`  <dbl> 2.86929, 3.09250, 2.86992, 3.07441, 2.92868, 3.02844,…
-## $ `s(year_fac).9`  <dbl> 3.22768, 3.28221, 3.24576, 3.19446, 3.31660, 3.23387,…
-## $ `s(year_fac).10` <dbl> 2.67766, 2.78781, 2.77803, 2.71490, 2.93688, 2.51834,…
-## $ `s(year_fac).11` <dbl> 2.96836, 3.08682, 3.10038, 3.11203, 3.04607, 3.12038,…
-## $ `s(year_fac).12` <dbl> 3.21117, 3.27094, 3.26856, 3.23706, 3.20997, 3.19841,…
-## $ `s(year_fac).13` <dbl> 2.28344, 2.24730, 2.46090, 2.14752, 2.34310, 2.28027,…
-## $ `s(year_fac).14` <dbl> 2.61860, 2.69168, 2.65335, 2.66072, 2.67606, 2.61563,…
-## $ `s(year_fac).15` <dbl> 2.30437, 2.08332, 2.27498, 2.14385, 2.08267, 2.26072,…
-## $ `s(year_fac).16` <dbl> 2.17932, 1.92979, 2.09785, 1.98914, 1.91563, 2.05489,…
-## $ `s(year_fac).17` <dbl> 0.535419, 2.245630, 1.379210, 1.009720, 1.429590, 1.0…
+## $ `s(year_fac).1` <dbl> 2.02996, 2.00453, 2.16828, 1.87209, 2.05206, 2.13466,… +## $ `s(year_fac).2` <dbl> 2.65617, 2.64468, 2.63625, 2.73117, 2.68130, 2.60230,… +## $ `s(year_fac).3` <dbl> 3.07788, 3.10003, 3.07758, 3.11288, 3.08044, 3.20094,… +## $ `s(year_fac).4` <dbl> 3.21960, 3.22782, 3.26233, 3.30094, 3.27737, 3.28554,… +## $ `s(year_fac).5` <dbl> 2.11108, 2.07728, 1.98290, 2.20085, 2.00719, 2.09721,… +## $ `s(year_fac).6` <dbl> 1.72869, 1.69724, 1.63484, 1.85265, 1.60296, 1.99641,… +## $ `s(year_fac).7` <dbl> 2.20549, 2.20062, 1.90837, 2.18498, 2.06963, 1.89268,… +## $ `s(year_fac).8` <dbl> 2.91909, 2.86867, 2.94035, 2.90241, 3.03399, 2.94243,… +## $ `s(year_fac).9` <dbl> 3.25630, 3.18218, 3.14735, 3.15474, 3.31125, 3.21197,… +## $ `s(year_fac).10` <dbl> 2.70871, 2.69267, 2.83379, 2.65286, 2.79554, 2.73977,… +## $ `s(year_fac).11` <dbl> 3.01297, 2.99048, 2.99582, 2.96736, 3.21242, 3.01926,… +## $ `s(year_fac).12` <dbl> 3.29335, 3.28567, 3.27609, 3.31069, 3.13586, 3.23824,… +## $ `s(year_fac).13` <dbl> 2.29235, 2.27421, 2.14434, 2.33456, 2.25957, 2.43291,… +## $ `s(year_fac).14` <dbl> 2.64613, 2.66115, 2.56605, 2.69002, 2.60428, 2.62748,… +## $ `s(year_fac).15` <dbl> 2.15095, 2.11449, 2.22458, 2.08566, 2.09324, 2.08875,… +## $ `s(year_fac).16` <dbl> 1.94051, 1.96340, 2.27971, 1.81916, 2.25314, 2.06480,… +## $ `s(year_fac).17` <dbl> 1.175650, 1.142240, 0.511516, 0.132152, 0.726240, 0.9…

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

@@ -840,7 +840,7 @@

distributions and diagnostics (see ?mvgam::mcmc_plot.mvgam for details):

-mcmc_plot(object = model1,
+mcmc_plot(object = model1,
           variable = 'betas',
           type = 'areas')

@@ -878,7 +878,7 @@

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

 hc <- hindcast(model1, type = 'link')
 range(hc$hindcasts$PP)
-
## [1] -3.01087  3.46325
+
## [1] -2.04395  3.45997

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

@@ -929,7 +929,7 @@ 

Automatic forecasting for new datamodel_data %>% dplyr::filter(time > 160) -> data_test

-model1b <- mvgam(count ~ s(year_fac, bs = 're') - 1,
+model1b <- mvgam(count ~ s(year_fac, bs = 're') - 1,
                 family = poisson(),
                 data = data_train,
                 newdata = data_test)
@@ -943,14 +943,14 @@

Automatic forecasting for new dataplot(model1b, type = 'forecast')

## Out of sample DRPS:
-## [1] 180.5596
+## [1] 182.1879

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.5596
+## [1] 182.1879

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 @@ -978,12 +978,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] 13 8 5 12 8 3 5 12 19 8 ... +## ..$ PP: num [1:2000, 1:160] 12 8 8 13 4 7 3 6 11 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] 9 11 2 9 9 11 6 6 16 11 ... +## ..$ PP: num [1:2000, 1:39] 5 3 8 8 12 5 11 9 4 9 ... ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : NULL ## .. .. ..$ : chr [1:39] "ypred[161,1]" "ypred[162,1]" "ypred[163,1]" "ypred[164,1]" ... @@ -998,7 +998,7 @@

Adding predictors as “fixed” eff Here, we will update the model from above by including a parametric (fixed) effect of ndvi as a linear predictor:

-model2 <- mvgam(count ~ s(year_fac, bs = 're') + 
+model2 <- mvgam(count ~ s(year_fac, bs = 're') + 
                   ndvi - 1,
                 family = poisson(),
                 data = data_train,
@@ -1040,33 +1040,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 1915 -## s(year_fac).1 1.10 1.40 1.70 1 3053 -## s(year_fac).2 1.80 2.00 2.20 1 2440 -## s(year_fac).3 2.20 2.40 2.60 1 2190 -## s(year_fac).4 2.30 2.50 2.70 1 2038 -## s(year_fac).5 1.20 1.40 1.60 1 2336 -## s(year_fac).6 1.00 1.30 1.50 1 2751 -## s(year_fac).7 1.10 1.40 1.70 1 2718 -## s(year_fac).8 2.10 2.30 2.50 1 2361 -## s(year_fac).9 2.70 2.90 3.00 1 2002 -## s(year_fac).10 2.00 2.20 2.40 1 2607 -## s(year_fac).11 2.30 2.40 2.60 1 2072 -## s(year_fac).12 2.50 2.70 2.80 1 2256 -## s(year_fac).13 1.40 1.60 1.90 1 2658 -## s(year_fac).14 0.61 2.00 3.30 1 1624 -## s(year_fac).15 0.55 2.00 3.40 1 1411 -## s(year_fac).16 0.63 2.00 3.40 1 1199 -## s(year_fac).17 0.62 2.00 3.30 1 1694 +## ndvi 0.32 0.39 0.46 1 1973 +## s(year_fac).1 1.10 1.40 1.70 1 2799 +## s(year_fac).2 1.80 2.00 2.20 1 2501 +## s(year_fac).3 2.20 2.40 2.60 1 2477 +## s(year_fac).4 2.30 2.50 2.70 1 2200 +## s(year_fac).5 1.20 1.40 1.60 1 2541 +## s(year_fac).6 1.00 1.30 1.50 1 2642 +## s(year_fac).7 1.10 1.40 1.70 1 2736 +## s(year_fac).8 2.10 2.30 2.50 1 2420 +## s(year_fac).9 2.70 2.90 3.00 1 2240 +## s(year_fac).10 2.00 2.20 2.40 1 2436 +## s(year_fac).11 2.30 2.40 2.60 1 2343 +## s(year_fac).12 2.50 2.70 2.80 1 2105 +## s(year_fac).13 1.40 1.60 1.90 1 2678 +## s(year_fac).14 0.68 1.90 3.20 1 1336 +## s(year_fac).15 0.54 2.00 3.20 1 1768 +## s(year_fac).16 0.67 2.00 3.20 1 1611 +## s(year_fac).17 0.68 2.00 3.20 1 1529 ## ## GAM group-level estimates: -## 2.5% 50% 97.5% Rhat n_eff -## mean(s(year_fac)) 1.60 2.00 2.30 1.01 390 -## sd(s(year_fac)) 0.42 0.61 0.98 1.00 500 +## 2.5% 50% 97.5% Rhat n_eff +## mean(s(year_fac)) 1.60 2.0 2.30 1.01 318 +## sd(s(year_fac)) 0.41 0.6 0.94 1.02 343 ## ## Approximate significance of GAM observation smooths: ## edf Chi.sq p-value -## s(year_fac) 10.8 2780 <2e-16 *** +## s(year_fac) 10.8 2921 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## @@ -1077,7 +1077,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 Mar 01 10:13:01 AM 2024. +## Samples were drawn using NUTS(diag_e) at Tue Mar 12 2:00:50 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)

@@ -1088,24 +1088,24 @@

Adding predictors as “fixed” eff
 coef(model2)
##                     2.5%      50%     97.5% Rhat n_eff
-## ndvi           0.3200159 0.389390 0.4594881    1  1915
-## s(year_fac).1  1.1233240 1.407850 1.6614208    1  3053
-## s(year_fac).2  1.7876195 1.998150 2.2032175    1  2440
-## s(year_fac).3  2.1899153 2.375350 2.5708740    1  2190
-## s(year_fac).4  2.3204982 2.503750 2.6854085    1  2038
-## s(year_fac).5  1.1790965 1.425910 1.6496230    1  2336
-## s(year_fac).6  1.0246777 1.270970 1.5133527    1  2751
-## s(year_fac).7  1.1365035 1.417170 1.6649067    1  2718
-## s(year_fac).8  2.0799478 2.272675 2.4598222    1  2361
-## s(year_fac).9  2.7261930 2.855000 2.9889848    1  2002
-## s(year_fac).10 1.9803723 2.180850 2.3823950    1  2607
-## s(year_fac).11 2.2559958 2.439730 2.6102330    1  2072
-## s(year_fac).12 2.5472203 2.691040 2.8458200    1  2256
-## s(year_fac).13 1.3686538 1.617025 1.8608865    1  2658
-## s(year_fac).14 0.6141420 1.992290 3.2861737    1  1624
-## s(year_fac).15 0.5509917 1.972045 3.3651875    1  1411
-## s(year_fac).16 0.6346040 1.980110 3.3646560    1  1199
-## s(year_fac).17 0.6161202 1.981280 3.2613995    1  1694
+## ndvi 0.3192765 0.390068 0.4565132 1 1973 +## s(year_fac).1 1.1297857 1.401075 1.6601210 1 2799 +## s(year_fac).2 1.8005140 2.002200 2.2031165 1 2501 +## s(year_fac).3 2.1920865 2.377540 2.5795502 1 2477 +## s(year_fac).4 2.3215133 2.507255 2.6850642 1 2200 +## s(year_fac).5 1.1973650 1.424105 1.6422465 1 2541 +## s(year_fac).6 1.0088363 1.273885 1.5184867 1 2642 +## s(year_fac).7 1.1395780 1.413890 1.6757230 1 2736 +## s(year_fac).8 2.0815710 2.273310 2.4560483 1 2420 +## s(year_fac).9 2.7174177 2.856920 2.9835225 1 2240 +## s(year_fac).10 1.9838457 2.184940 2.3871055 1 2436 +## s(year_fac).11 2.2670628 2.437260 2.6049930 1 2343 +## s(year_fac).12 2.5429085 2.690940 2.8465860 1 2105 +## s(year_fac).13 1.3885838 1.618830 1.8522713 1 2678 +## s(year_fac).14 0.6765426 1.941810 3.2150647 1 1336 +## s(year_fac).15 0.5352354 2.000345 3.1917455 1 1768 +## s(year_fac).16 0.6744930 1.988625 3.2106910 1 1611 +## s(year_fac).17 0.6755871 1.958040 3.1927537 1 1529

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

@@ -1121,24 +1121,24 @@ 

Adding predictors as “fixed” eff dplyr::glimpse(beta_post)

## Rows: 2,000
 ## Columns: 18
-## $ ndvi             <dbl> 0.410494, 0.413662, 0.442400, 0.402106, 0.384262, 0.3…
-## $ `s(year_fac).1`  <dbl> 1.13651, 1.35671, 1.35124, 1.28985, 1.55463, 1.34357,…
-## $ `s(year_fac).2`  <dbl> 1.97558, 1.96477, 1.93552, 1.95374, 2.02403, 1.96302,…
-## $ `s(year_fac).3`  <dbl> 2.35366, 2.32803, 2.27197, 2.32897, 2.44064, 2.39664,…
-## $ `s(year_fac).4`  <dbl> 2.52377, 2.37831, 2.35059, 2.55368, 2.46688, 2.61205,…
-## $ `s(year_fac).5`  <dbl> 1.45789, 1.49795, 1.25727, 1.28905, 1.54285, 1.46934,…
-## $ `s(year_fac).6`  <dbl> 1.19230, 1.25613, 1.20471, 1.43570, 1.22353, 1.29816,…
-## $ `s(year_fac).7`  <dbl> 1.26501, 1.22271, 1.51465, 1.43871, 1.32029, 1.40930,…
-## $ `s(year_fac).8`  <dbl> 2.38418, 2.22531, 2.21146, 2.20142, 2.36365, 2.19256,…
-## $ `s(year_fac).9`  <dbl> 2.77094, 2.88541, 2.92821, 2.75575, 2.88332, 2.90387,…
-## $ `s(year_fac).10` <dbl> 2.12203, 2.12384, 2.12692, 2.06942, 2.21348, 2.25463,…
-## $ `s(year_fac).11` <dbl> 2.33057, 2.39559, 2.40417, 2.38458, 2.49063, 2.45352,…
-## $ `s(year_fac).12` <dbl> 2.67567, 2.58010, 2.64711, 2.70238, 2.66524, 2.78719,…
-## $ `s(year_fac).13` <dbl> 1.56129, 1.52968, 1.62966, 1.62109, 1.64491, 1.68440,…
-## $ `s(year_fac).14` <dbl> 2.54425, 2.16076, 2.49400, 2.24232, 2.41177, 2.77880,…
-## $ `s(year_fac).15` <dbl> 2.796810, 2.845590, 1.648850, 2.344600, 2.372690, 2.1…
-## $ `s(year_fac).16` <dbl> 3.14343, 3.05613, 2.80045, 3.09318, 3.41703, 3.54078,…
-## $ `s(year_fac).17` <dbl> 1.61158, 1.77748, 1.56756, 1.69925, 1.89734, 1.05068,…
+## $ ndvi <dbl> 0.384510, 0.399786, 0.392306, 0.349502, 0.422103, 0.4… +## $ `s(year_fac).1` <dbl> 1.28841, 1.51449, 1.25610, 1.48108, 1.36640, 1.28543,… +## $ `s(year_fac).2` <dbl> 1.88972, 2.10276, 2.01286, 1.97798, 1.86941, 1.92906,… +## $ `s(year_fac).3` <dbl> 2.31917, 2.44870, 2.31533, 2.47928, 2.35820, 2.24181,… +## $ `s(year_fac).4` <dbl> 2.52986, 2.49321, 2.54039, 2.64189, 2.35362, 2.42710,… +## $ `s(year_fac).5` <dbl> 1.35578, 1.51313, 1.33385, 1.39870, 1.35560, 1.29213,… +## $ `s(year_fac).6` <dbl> 1.311920, 1.386790, 1.239250, 1.148290, 1.233900, 1.1… +## $ `s(year_fac).7` <dbl> 1.52013, 1.09664, 1.13217, 1.59548, 1.43845, 1.20762,… +## $ `s(year_fac).8` <dbl> 2.27497, 2.27513, 2.18857, 2.36107, 2.20560, 2.03529,… +## $ `s(year_fac).9` <dbl> 2.86233, 2.82828, 2.96744, 3.00598, 2.93260, 2.84276,… +## $ `s(year_fac).10` <dbl> 2.06303, 2.21389, 2.17205, 2.20981, 2.09964, 2.09335,… +## $ `s(year_fac).11` <dbl> 2.44185, 2.45756, 2.40478, 2.59571, 2.31021, 2.34485,… +## $ `s(year_fac).12` <dbl> 2.62680, 2.72703, 2.54081, 2.63102, 2.65908, 2.52744,… +## $ `s(year_fac).13` <dbl> 1.76098, 1.67823, 1.48126, 1.67847, 1.52566, 1.51016,… +## $ `s(year_fac).14` <dbl> 3.678590, 2.295410, 1.691400, 1.657580, 2.371780, 1.0… +## $ `s(year_fac).15` <dbl> 1.64002, 1.65026, 1.44489, 0.89341, 2.02970, 2.42390,… +## $ `s(year_fac).16` <dbl> 1.707470, 2.109060, 2.077640, 2.213590, 3.446480, 2.3… +## $ `s(year_fac).17` <dbl> 2.13327, 1.18055, 1.40620, 1.41307, 2.22927, 1.95107,…

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 @@ -1169,10 +1169,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
@@ -1186,7 +1186,7 @@ 

plots for main effects. This will likely be your go-to function for quickly understanding patterns from fitted mvgam models

-plot(conditional_effects(model2), ask = FALSE)
+plot(conditional_effects(model2), ask = FALSE)

@@ -1224,7 +1224,7 @@

Adding predictors as smooths
-model3 <- mvgam(count ~ s(time, bs = 'bs', k = 15) + 
+model3 <- mvgam(count ~ s(time, bs = 'bs', k = 15) + 
                   ndvi,
                 family = poisson(),
                 data = data_train,
@@ -1274,27 +1274,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 977 -## ndvi 0.26 0.33 0.400 1.00 964 -## s(time).1 -2.20 -1.10 -0.003 1.00 530 -## s(time).2 0.43 1.30 2.300 1.01 433 -## s(time).3 -0.50 0.49 1.400 1.01 383 -## s(time).4 1.60 2.50 3.400 1.01 384 -## s(time).5 -1.20 -0.18 0.800 1.01 389 -## s(time).6 -0.61 0.40 1.500 1.01 409 -## s(time).7 -1.50 -0.50 0.500 1.01 414 -## s(time).8 0.58 1.50 2.500 1.01 393 -## s(time).9 1.20 2.10 3.000 1.01 349 -## s(time).10 -0.42 0.58 1.600 1.01 382 -## s(time).11 0.81 1.80 2.700 1.01 367 -## s(time).12 0.66 1.50 2.400 1.01 391 -## s(time).13 -1.20 -0.29 0.650 1.00 530 -## s(time).14 -7.30 -4.30 -1.100 1.01 520 +## 2.5% 50% 97.5% Rhat n_eff +## (Intercept) 2.00 2.10 2.2000 1.00 1045 +## ndvi 0.26 0.33 0.4000 1.00 1056 +## s(time).1 -2.20 -1.10 -0.0021 1.00 475 +## s(time).2 0.39 1.30 2.3000 1.00 414 +## s(time).3 -0.62 0.45 1.5000 1.00 362 +## s(time).4 1.50 2.50 3.5000 1.00 364 +## s(time).5 -1.30 -0.24 0.7600 1.00 375 +## s(time).6 -0.65 0.38 1.4000 1.00 403 +## s(time).7 -1.60 -0.53 0.4700 1.01 405 +## s(time).8 0.50 1.50 2.5000 1.00 356 +## s(time).9 1.10 2.10 3.1000 1.01 358 +## s(time).10 -0.39 0.55 1.6000 1.00 363 +## s(time).11 0.75 1.70 2.8000 1.00 355 +## s(time).12 0.59 1.50 2.4000 1.00 379 +## s(time).13 -1.20 -0.32 0.6300 1.00 441 +## s(time).14 -7.40 -4.10 -1.0000 1.01 398 ## ## Approximate significance of GAM observation smooths: ## edf Chi.sq p-value -## s(time) 8.23 720 <2e-16 *** +## s(time) 9.54 771 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## @@ -1305,7 +1305,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 Mar 01 10:13:55 AM 2024. +## Samples were drawn using NUTS(diag_e) at Tue Mar 12 2:01:29 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)

@@ -1346,11 +1346,11 @@

Conditional effects
-plot(conditional_effects(model3), ask = FALSE)
+plot(conditional_effects(model3), ask = FALSE)

Or on the link scale:

-plot(conditional_effects(model3, type = 'link'), ask = FALSE)
+plot(conditional_effects(model3, type = 'link'), ask = FALSE)

Inspect the underlying Stan code to gain some idea of how the spline is being penalized:

@@ -1432,7 +1432,7 @@

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

## Out of sample DRPS:
-## [1] 287.7197
+## [1] 286.148

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 @@ -1465,7 +1465,7 @@

Latent dynamics in mvgammvgam can include combinations of smooths and dynamic components:

-model4 <- mvgam(count ~ s(ndvi, k = 6),
+model4 <- mvgam(count ~ s(ndvi, k = 6),
                 family = poisson(),
                 data = data_train,
                 newdata = data_test,
@@ -1514,40 +1514,40 @@ 

Latent dynamics in mvgam## ## ## GAM coefficient (beta) estimates: -## 2.5% 50% 97.5% Rhat n_eff -## (Intercept) 0.0066 2.0000 2.600 1.39 11 -## s(ndvi).1 -0.0790 0.0110 0.180 1.01 536 -## s(ndvi).2 -0.1600 0.0150 0.280 1.01 277 -## s(ndvi).3 -0.0550 -0.0009 0.048 1.00 888 -## s(ndvi).4 -0.2600 0.1100 1.300 1.01 197 -## s(ndvi).5 -0.0700 0.1500 0.360 1.01 599 +## 2.5% 50% 97.5% Rhat n_eff +## (Intercept) 1.200 1.90000 2.500 1.07 74 +## s(ndvi).1 -0.087 0.01300 0.200 1.00 400 +## s(ndvi).2 -0.140 0.02000 0.330 1.01 358 +## s(ndvi).3 -0.056 -0.00094 0.049 1.02 469 +## s(ndvi).4 -0.250 0.13000 1.300 1.01 252 +## s(ndvi).5 -0.092 0.15000 0.360 1.00 467 ## ## Approximate significance of GAM observation smooths: -## edf Chi.sq p-value -## s(ndvi) 0.981 87.3 0.079 . +## edf Chi.sq p-value +## s(ndvi) 0.88 95.3 0.065 . ## --- ## 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.82 0.95 1.13 29 -## sigma[1] 0.68 0.80 0.96 1.00 431 +## ar1[1] 0.69 0.81 0.92 1.01 340 +## sigma[1] 0.68 0.80 0.96 1.01 434 ## ## Stan MCMC diagnostics: ## n_eff / iter looks reasonable for all parameters -## Rhats above 1.05 found for 163 parameters +## Rhats above 1.05 found for 33 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 Fri Mar 01 10:15:13 AM 2024. +## Samples were drawn using NUTS(diag_e) at Tue Mar 12 2:02:41 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)

@@ -1556,8 +1556,8 @@

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

-
## Out of sample DRPS:
-## [1] 154.5616
+
## Out of sample CRPS:
+## [1] 150.9749

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

 plot(model4, type = 'trend', newdata = data_test)
@@ -1566,13 +1566,13 @@

Latent dynamics in mvgamloo package (a higher value is preferred for this metric):

-loo_compare(model3, model4)
+
loo_compare(model3, model4)
## Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.
 
 ## Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.
##        elpd_diff se_diff
 ## model4    0.0       0.0 
-## model3 -560.5      66.8
+## model3 -562.9 67.3

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

@@ -1587,7 +1587,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] -133.158
+
## [1] -135.1958

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 61451058..e49ba75c 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 0cd0d9ac..80917602 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 6848ca0c..9190631c 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 f4d05def..61c9225c 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 cb23b7de..58160331 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 ad7237d1..f615aa95 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 9026c4b9..bd8d7ab2 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 980d720b..71b5967c 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 6eab1f16..e5d295a6 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 f829d810..ebd41a7e 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 e461fccc..9c1f965a 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 6eab1f16..e5d295a6 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 f21b5c95..149d2b0e 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 86928d00..497c62be 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 e4c8a89f..d2cde447 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 40cd51b8..8b42d361 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 8f9893a6..730cb280 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 ef96065e..ee186f33 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 6df735dc..37919630 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 9b4ff3e1..7a89a976 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 18e4e551..a7a03e62 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 cadef279..e9913189 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 e2ab2408..c7f25966 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 9875807e..beea1bc7 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 85d59c7c..ccf4a3af 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 156e5388..2c373404 100644 --- a/docs/index.html +++ b/docs/index.html @@ -84,6 +84,7 @@



+

brms LogoStan Logo

mvgam

diff --git a/docs/pkgdown.yml b/docs/pkgdown.yml index 5e2a8923..98d98d3c 100644 --- a/docs/pkgdown.yml +++ b/docs/pkgdown.yml @@ -9,7 +9,7 @@ articles: shared_states: shared_states.html time_varying_effects: time_varying_effects.html trend_formulas: trend_formulas.html -last_built: 2024-01-29T05:37Z +last_built: 2024-03-12T03:54Z urls: reference: https://nicholasjclark.github.io/mvgam/reference article: https://nicholasjclark.github.io/mvgam/articles diff --git a/docs/reference/figures/README-unnamed-chunk-13-1.png b/docs/reference/figures/README-unnamed-chunk-13-1.png index 97f08209..63559a6a 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 8e34a550..1dca9500 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 f801d30b..8b1142de 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 94d4bf12..ba4439e4 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 85e0db39..d722ba66 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 fbb3a6b1..c72e0706 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 7f2f84c9..87c0bc2e 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 84f2b09d..d39944b0 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 246e3eef..eb759d89 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 e1b4751a..a748b548 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 00f1d110..45ab957c 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 24d051b6..40e68091 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/docs/reference/figures/mvgam_logo.png b/docs/reference/figures/mvgam_logo.png new file mode 100644 index 00000000..ce0ace9b Binary files /dev/null and b/docs/reference/figures/mvgam_logo.png differ diff --git a/docs/reference/get_mvgam_priors.html b/docs/reference/get_mvgam_priors.html index 273505f3..a514ddc0 100644 --- a/docs/reference/get_mvgam_priors.html +++ b/docs/reference/get_mvgam_priors.html @@ -402,13 +402,13 @@

Examples#> 6 trend sd sigma ~ student_t(3, 0, 2.5); #> 7 inverse of NB dispsersion phi_inv ~ student_t(3, 0, 0.1); #> example_change new_lowerbound new_upperbound -#> 1 lambda ~ exponential(0.03); NA NA -#> 2 mu_raw ~ normal(0.56, 0.85); NA NA -#> 3 sigma_raw ~ exponential(0.28); NA NA -#> 4 ar1 ~ normal(0.03, 0.94); NA NA -#> 5 ar2 ~ normal(0.08, 0.56); NA NA -#> 6 sigma ~ exponential(0.97); NA NA -#> 7 phi_inv ~ normal(0.75, 0.12); NA NA +#> 1 lambda ~ exponential(0.18); NA NA +#> 2 mu_raw ~ normal(0.23, 0.71); NA NA +#> 3 sigma_raw ~ exponential(0.79); NA NA +#> 4 ar1 ~ normal(-0.57, 0.81); NA NA +#> 5 ar2 ~ normal(0.19, 0.65); NA NA +#> 6 sigma ~ exponential(0.12); NA NA +#> 7 phi_inv ~ normal(-0.61, 0.11); NA NA # Make a few changes; first, change the population mean for the series-level # random intercepts diff --git a/docs/reference/monotonic.html b/docs/reference/monotonic.html index c243b75d..368ad6d9 100644 --- a/docs/reference/monotonic.html +++ b/docs/reference/monotonic.html @@ -64,16 +64,16 @@

Usage

# S3 method for moi.smooth.spec
-smooth.construct(object, data, knots)
+smooth.construct(object, data, knots)
 
 # S3 method for mod.smooth.spec
-smooth.construct(object, data, knots)
+smooth.construct(object, data, knots)
 
 # S3 method for moi.smooth
-Predict.matrix(object, data)
+Predict.matrix(object, data)
 
 # S3 method for mod.smooth
-Predict.matrix(object, data)
+Predict.matrix(object, data)

@@ -152,7 +152,7 @@

Examples # A standard TRPS smooth doesn't capture monotonicity mod_data <- data.frame(y = y, x = x) -mod <- gam(y ~ s(x, k = 16), +mod <- gam(y ~ s(x, k = 16), data = mod_data, family = gaussian()) @@ -166,7 +166,7 @@

Examples # Using the 'moi' basis in mvgam rectifies this mod_data$time <- 1:NROW(mod_data) -mod2 <- mvgam(y ~ s(x, bs = 'moi', k = 18), +mod2 <- mvgam(y ~ s(x, bs = 'moi', k = 18), data = mod_data, family = gaussian()) @@ -199,7 +199,7 @@

Examplesmod_data$time <- 1:NROW(mod_data) # Fit a model with different smooths per factor level -mod <- mvgam(y ~ s(x, bs = 'moi', by = fac, k = 8), +mod <- mvgam(y ~ s(x, bs = 'moi', by = fac, k = 8), data = mod_data, family = gaussian()) diff --git a/docs/reference/mvgam.html b/docs/reference/mvgam.html index a2c50701..cdc9d72f 100644 --- a/docs/reference/mvgam.html +++ b/docs/reference/mvgam.html @@ -91,6 +91,7 @@

Usage prior_simulation = FALSE, return_model_data = FALSE, family = "poisson", + share_obs_params = FALSE, use_lv = FALSE, n_lv, trend_map, @@ -122,7 +123,7 @@

Usage

Arguments

formula

A character string specifying the GAM observation model formula. These are exactly like the formula -for a GLM except that smooth terms, s(), te(), ti(), t2(), as well as time-varying +for a GLM except that smooth terms, s(), te(), ti(), t2(), as well as time-varying dynamic() terms, can be added to the right hand side to specify that the linear predictor depends on smooth functions of predictors (or linear functionals of these). In nmix() family models, the formula is used to @@ -201,22 +202,31 @@

Argumentsnb() for count data

  • poisson() for count data

  • gaussian() for real-valued data

  • -
  • betar() for proportional data on (0,1)

  • -
  • lognormal() for non-negative real-valued data

  • +
  • betar() for proportional data on (0,1)

  • +
  • lognormal() for non-negative real-valued data

  • student_t() for real-valued data

  • Gamma() for non-negative real-valued data

  • nmix() for count data with imperfect detection modeled via a State-Space N-Mixture model. The latent states are Poisson, capturing the 'true' latent abundance, while the observation process is Binomial to account for imperfect detection. See mvgam_families for an example of how to use this family

  • -

    Note that only nb() and poisson() are available if using JAGS as the backend. +

    Note that only nb() and poisson() are available if using JAGS as the backend. Default is poisson(). See mvgam_families for more details

    +
    share_obs_params
    +

    logical. If TRUE and the family +has additional family-specific observation parameters (e.g. variance components in +student_t() or gaussian(), or dispersion parameters in nb() or betar()), +these parameters will be shared across all series. This is handy if you have multiple +time series that you believe share some properties, such as being from the same +species over different spatial units. Default is FALSE.

    + +
    use_lv

    logical. If TRUE, use dynamic factors to estimate series' latent trends in a reduced dimension format. Only available for @@ -453,8 +463,9 @@

    Details
    Observation level parameters: When more than one series is included in data and an observation family that contains more than one parameter is used, additional observation family parameters -(i.e. phi for nb() or sigma for
    gaussian()) are -estimated independently for each series. +(i.e. phi for nb() or sigma for gaussian()) are +by default estimated independently for each series. But if you wish for the series to share +the same observation parameters, set share_obs_params = TRUE

    Factor regularisation: When using a dynamic factor model for the trends with JAGS factor precisions are given regularized penalty priors to theoretically allow some factors to be dropped from the model by squeezing increasing factors' variances to zero. This is done to help protect against selecting too many latent factors than are needed to @@ -507,7 +518,7 @@

    Examples# Formulate a model using Stan where series share a cyclic smooth for # seasonality and each series has an independent random walk temporal process; # Set run_model = FALSE to inspect the returned objects -mod1 <- mvgam(formula = y ~ s(season, bs = 'cc'), +mod1 <- mvgam(formula = y ~ s(season, bs = 'cc'), data = dat$data_train, trend_model = 'RW', family = 'poisson', @@ -538,7 +549,7 @@

    Examples refresh = 100) # Now fit the model using mvgam with the Stan backend -mod1 <- mvgam(formula = y ~ s(season, bs = 'cc'), +mod1 <- mvgam(formula = y ~ s(season, bs = 'cc'), data = dat$data_train, trend_model = 'RW', family = poisson(), @@ -574,7 +585,7 @@

    Examplesplot(mod1, type = 'smooths', realisations = TRUE) # Plot conditional response predictions using marginaleffects -plot(conditional_effects(mod1), ask = FALSE) +plot(conditional_effects(mod1), ask = FALSE) plot_predictions(mod1, condition = 'season', points = 0.5) # Extract observation model beta coefficient draws as a data.frame @@ -597,7 +608,7 @@

    Examples trend = c(1,1,2)) # Fit the model using AR1 trends -mod <- mvgam(y ~ s(season, bs = 'cc'), +mod <- mvgam(y ~ s(season, bs = 'cc'), trend_map = trend_map, trend_model = 'AR1', data = mod_data, @@ -662,9 +673,9 @@

    Examples# Fit a model that includes the offset in the linear predictor as well as # hierarchical seasonal smooths mod <- mvgam(formula = y ~ offset(offset) + - s(series, bs = 're') + - s(season, bs = 'cc') + - s(season, by = series, m = 1, k = 5), + s(series, bs = 're') + + s(season, bs = 'cc') + + s(season, by = series, m = 1, k = 5), data = dat$data_train, trend_model = 'None', use_stan = TRUE) @@ -765,7 +776,7 @@

    Examples # View the changepoints with ggplot2 utilities library(ggplot2) -mcmc_plot(mod, variable = 'delta_trend', +mcmc_plot(mod, variable = 'delta_trend', regex = TRUE) + scale_y_discrete(labels = mod$trend_model$changepoints) + labs(y = 'Potential changepoint', diff --git a/docs/reference/update.mvgam.html b/docs/reference/update.mvgam.html index c3b8b5a0..823f41e0 100644 --- a/docs/reference/update.mvgam.html +++ b/docs/reference/update.mvgam.html @@ -72,6 +72,7 @@

    Usage use_lv, n_lv, family, + share_obs_params, priors, lfo = FALSE, ... @@ -175,6 +176,15 @@

    Argumentsmvgam_families for more details

    +
    share_obs_params
    +

    logical. If TRUE and the family +has additional family-specific observation parameters (e.g. variance components in +student_t() or gaussian(), or dispersion parameters in nb() or betar()), +these parameters will be shared across all series. This is handy if you have multiple +time series that you believe share some properties, such as being from the same +species over different spatial units. Default is FALSE.

    + +
    priors

    An optional data.frame with prior definitions (in JAGS or Stan syntax). if using Stan, this can also be an object of diff --git a/index.Rmd b/index.Rmd index e7112bb5..1e409a8d 100644 --- a/index.Rmd +++ b/index.Rmd @@ -5,6 +5,8 @@ always_allow_html: true

    +brms Logo[Stan Logo](https://mc-stan.org/) + ## mvgam **M**ulti**V**ariate (Dynamic) **G**eneralized **A**ddivite **M**odels diff --git a/index.md b/index.md index 6441ce50..2fbd2c11 100644 --- a/index.md +++ b/index.md @@ -1,6 +1,8 @@

    +brms Logo[Stan Logo](https://mc-stan.org/) + ## mvgam **M**ulti**V**ariate (Dynamic) **G**eneralized **A**ddivite **M**odels diff --git a/man/figures/README-unnamed-chunk-13-1.png b/man/figures/README-unnamed-chunk-13-1.png index 97f08209..63559a6a 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 8e34a550..1dca9500 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 f801d30b..8b1142de 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 94d4bf12..ba4439e4 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 85e0db39..d722ba66 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 fbb3a6b1..c72e0706 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 7f2f84c9..87c0bc2e 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 84f2b09d..d39944b0 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 246e3eef..eb759d89 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 e1b4751a..a748b548 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 00f1d110..45ab957c 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 24d051b6..40e68091 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/man/figures/mvgam_logo.png b/man/figures/mvgam_logo.png new file mode 100644 index 00000000..ce0ace9b Binary files /dev/null and b/man/figures/mvgam_logo.png differ