diff --git a/README.md b/README.md index 46c86243..26e4a39a 100644 --- a/README.md +++ b/README.md @@ -234,29 +234,29 @@ summary(lynx_mvgam) #> #> #> GAM coefficient (beta) estimates: -#> 2.5% 50% 97.5% Rhat n_eff -#> (Intercept) 6.4000 6.60 6.900 1 718 -#> s(season).1 -0.6500 -0.13 0.420 1 987 -#> s(season).2 0.7600 1.30 1.900 1 865 -#> s(season).3 1.3000 1.90 2.600 1 816 -#> s(season).4 -0.0360 0.53 1.100 1 909 -#> s(season).5 -1.4000 -0.69 -0.096 1 759 -#> s(season).6 -1.3000 -0.56 0.150 1 869 -#> s(season).7 0.0092 0.71 1.400 1 970 -#> s(season).8 0.6100 1.40 2.000 1 802 -#> s(season).9 -0.3500 0.21 0.840 1 897 -#> s(season).10 -1.4000 -0.86 -0.370 1 1366 +#> 2.5% 50% 97.5% Rhat n_eff +#> (Intercept) 6.400 6.60 6.900 1.00 961 +#> s(season).1 -0.650 -0.14 0.380 1.00 984 +#> s(season).2 0.730 1.30 1.900 1.00 1113 +#> s(season).3 1.300 1.90 2.500 1.00 727 +#> s(season).4 -0.054 0.53 1.100 1.00 727 +#> s(season).5 -1.300 -0.70 -0.069 1.00 919 +#> s(season).6 -1.200 -0.55 0.170 1.00 1012 +#> s(season).7 0.012 0.73 1.400 1.00 1134 +#> s(season).8 0.600 1.40 2.100 1.01 756 +#> s(season).9 -0.410 0.22 0.810 1.00 710 +#> s(season).10 -1.400 -0.87 -0.350 1.00 1164 #> #> Approximate significance of GAM smooths: #> edf Ref.df Chi.sq p-value -#> s(season) 9.96 10 49.6 <2e-16 *** +#> s(season) 9.95 10 53.4 <2e-16 *** #> --- #> 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.57 0.83 0.98 1.01 659 -#> sigma[1] 0.38 0.48 0.60 1.01 745 +#> ar1[1] 0.59 0.83 0.99 1.00 519 +#> sigma[1] 0.38 0.48 0.60 1.01 485 #> #> Stan MCMC diagnostics: #> n_eff / iter looks reasonable for all parameters @@ -265,7 +265,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 Thu Jul 25 3:12:50 PM 2024. +#> Samples were drawn using NUTS(diag_e) at Tue Sep 03 1:56: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) @@ -370,14 +370,12 @@ partial effects of smooths on the link scale ``` r require(gratia) #> Loading required package: gratia +#> Warning: package 'gratia' was built under R version 4.2.3 #> #> Attaching package: 'gratia' #> The following object is masked from 'package:mvgam': #> #> add_residuals -#> The following object is masked from 'package:brms': -#> -#> posterior_samples draw(lynx_mvgam) ``` @@ -405,7 +403,7 @@ series (testing and training) ``` r plot(lynx_mvgam, type = 'forecast', newdata = lynx_test) #> Out of sample DRPS: -#> 2380.85453525 +#> 2356.32008125 ``` Plotting forecast distributions using mvgam in R @@ -535,41 +533,41 @@ summary(mod, include_betas = FALSE) #> #> Observation precision parameter estimates: #> 2.5% 50% 97.5% Rhat n_eff -#> phi[1] 5.4 8.3 12 1 1248 -#> phi[2] 5.7 8.6 13 1 1312 -#> phi[3] 5.6 8.5 12 1 1724 +#> phi[1] 5.4 8.3 12 1 2019 +#> phi[2] 5.7 8.7 13 1 1795 +#> phi[3] 5.5 8.4 12 1 1680 #> #> GAM coefficient (beta) estimates: -#> 2.5% 50% 97.5% Rhat n_eff -#> (Intercept) -0.2 0.19 0.46 1.01 566 +#> 2.5% 50% 97.5% Rhat n_eff +#> (Intercept) -0.18 0.19 0.45 1.01 757 #> #> Approximate significance of GAM smooths: #> edf Ref.df Chi.sq p-value -#> s(season) 3.872 5 29.63 1.6e-05 *** -#> s(season):seriesseries_1 0.615 4 0.77 0.98 -#> s(season):seriesseries_2 1.012 4 0.30 0.99 -#> s(season):seriesseries_3 1.106 4 1.54 0.81 +#> s(season) 4.231 5 28.89 2.6e-06 *** +#> s(season):seriesseries_1 0.687 4 0.69 0.98 +#> s(season):seriesseries_2 0.664 4 0.60 0.99 +#> s(season):seriesseries_3 1.286 4 1.39 0.82 #> --- #> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 #> #> Latent trend marginal deviation (alpha) and length scale (rho) estimates: -#> 2.5% 50% 97.5% Rhat n_eff -#> alpha_gp[1] 0.051 0.41 0.92 1.01 525 -#> alpha_gp[2] 0.360 0.72 1.20 1.00 946 -#> alpha_gp[3] 0.150 0.46 1.00 1.00 659 -#> rho_gp[1] 1.100 3.80 15.00 1.01 370 -#> rho_gp[2] 1.900 7.80 37.00 1.01 365 -#> rho_gp[3] 1.400 5.10 21.00 1.00 645 +#> 2.5% 50% 97.5% Rhat n_eff +#> alpha_gp[1] 0.06 0.41 0.96 1.00 823 +#> alpha_gp[2] 0.38 0.73 1.30 1.00 1073 +#> alpha_gp[3] 0.15 0.45 0.98 1.00 903 +#> rho_gp[1] 1.10 3.90 14.00 1.01 302 +#> rho_gp[2] 1.70 7.00 32.00 1.01 364 +#> rho_gp[3] 1.30 4.80 22.00 1.01 373 #> #> Stan MCMC diagnostics: #> n_eff / iter looks reasonable for all parameters #> Rhat looks reasonable for all parameters -#> 12 of 2000 iterations ended with a divergence (0.6%) +#> 18 of 2000 iterations ended with a divergence (0.9%) #> *Try running with larger adapt_delta to remove the divergences #> 0 of 2000 iterations saturated the maximum tree depth of 12 (0%) #> E-FMI indicated no pathological behavior #> -#> Samples were drawn using NUTS(diag_e) at Thu Jul 25 3:13:40 PM 2024. +#> Samples were drawn using NUTS(diag_e) at Tue Sep 03 1:57: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) diff --git a/docs/articles/data_in_mvgam.html b/docs/articles/data_in_mvgam.html index ff06e687..17989235 100644 --- a/docs/articles/data_in_mvgam.html +++ b/docs/articles/data_in_mvgam.html @@ -81,7 +81,7 @@

Nicholas J Clark

-

2024-08-24

+

2024-09-03

Source: vignettes/data_in_mvgam.Rmd
data_in_mvgam.Rmd
@@ -111,21 +111,21 @@

Required long data formathead(simdat$data_train, 16) #> y season year series time #> 1 NA 1 1 series_1 1 -#> 2 0 1 1 series_2 1 +#> 2 1 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 2 2 1 series_2 2 +#> 4 NA 1 1 series_4 1 +#> 5 1 2 1 series_1 2 +#> 6 NA 2 1 series_2 2 #> 7 1 2 1 series_3 2 #> 8 NA 2 1 series_4 2 #> 9 0 3 1 series_1 3 -#> 10 NA 3 1 series_2 3 +#> 10 3 3 1 series_2 3 #> 11 1 3 1 series_3 3 -#> 12 2 3 1 series_4 3 -#> 13 0 4 1 series_1 4 +#> 12 0 3 1 series_4 3 +#> 13 NA 4 1 series_1 4 #> 14 2 4 1 series_2 4 -#> 15 0 4 1 series_3 4 -#> 16 3 4 1 series_4 4 +#> 15 1 4 1 series_3 4 +#> 16 1 4 1 series_4 4

series as a factor variable @@ -176,28 +176,28 @@

A single outcome variable#> #> Deviance Residuals: #> Min 1Q Median 3Q Max -#> -2.3329 -1.2350 -0.6016 0.3510 2.5133 +#> -1.7390 -0.9456 -0.3864 0.7284 2.8061 #> #> Coefficients: -#> Estimate Std. Error z value Pr(>|z|) -#> (Intercept) -1.24362 0.61486 -2.023 0.043114 * -#> seriesseries_2 2.27379 0.60407 3.764 0.000167 *** -#> seriesseries_3 1.36576 0.63696 2.144 0.032018 * -#> seriesseries_4 2.32647 0.60468 3.847 0.000119 *** -#> time -0.02912 0.02156 -1.351 0.176694 +#> Estimate Std. Error z value Pr(>|z|) +#> (Intercept) 0.25386 0.36697 0.692 0.4891 +#> seriesseries_2 0.60796 0.34442 1.765 0.0775 . +#> seriesseries_3 -1.11102 0.52657 -2.110 0.0349 * +#> seriesseries_4 0.01518 0.39240 0.039 0.9691 +#> time -0.03449 0.02642 -1.305 0.1918 #> --- #> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 #> #> (Dispersion parameter for poisson family taken to be 1) #> -#> Null deviance: 134.349 on 58 degrees of freedom -#> Residual deviance: 96.445 on 54 degrees of freedom -#> (13 observations deleted due to missingness) -#> AIC: 190.39 +#> Null deviance: 79.481 on 57 degrees of freedom +#> Residual deviance: 61.529 on 53 degrees of freedom +#> (14 observations deleted due to missingness) +#> AIC: 148.43 #> -#> Number of Fisher Scoring iterations: 6

+#> Number of Fisher Scoring iterations: 5
-summary(mgcv::gam(y ~ series + s(time, by = series),
+summary(mgcv::gam(y ~ series + s(time, by = series),
             data = simdat$data_train,
             family = poisson()))
 #> 
@@ -208,25 +208,23 @@ 

A single outcome variable#> y ~ series + s(time, by = series) #> #> Parametric coefficients: -#> Estimate Std. Error z value Pr(>|z|) -#> (Intercept) -1.6035 0.6005 -2.671 0.00757 ** -#> seriesseries_2 1.9155 0.6601 2.902 0.00371 ** -#> seriesseries_3 1.2765 0.6790 1.880 0.06012 . -#> seriesseries_4 2.0982 0.6508 3.224 0.00126 ** +#> Estimate Std. Error z value Pr(>|z|) +#> (Intercept) -0.1536 0.2918 -0.526 0.5987 +#> seriesseries_2 0.6781 0.3582 1.893 0.0584 . +#> seriesseries_3 -2.2093 1.1016 -2.005 0.0449 * +#> seriesseries_4 -0.3154 0.4940 -0.639 0.5231 #> --- #> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 #> #> Approximate significance of smooth terms: -#> edf Ref.df Chi.sq p-value -#> s(time):seriesseries_1 1.494 1.833 0.575 0.73049 -#> s(time):seriesseries_2 2.591 3.278 13.476 0.00574 ** -#> s(time):seriesseries_3 4.750 5.799 9.781 0.12285 -#> s(time):seriesseries_4 4.592 5.621 13.802 0.03315 * -#> --- -#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 +#> edf Ref.df Chi.sq p-value +#> s(time):seriesseries_1 2.752 3.397 1.887 0.608 +#> s(time):seriesseries_2 1.236 1.436 1.218 0.507 +#> s(time):seriesseries_3 4.664 5.653 4.517 0.580 +#> s(time):seriesseries_4 3.461 4.274 6.525 0.194 #> -#> R-sq.(adj) = 0.628 Deviance explained = 68.3% -#> UBRE = 0.31332 Scale est. = 1 n = 59

+#> R-sq.(adj) = 0.248 Deviance explained = 50.5% +#> UBRE = 0.23426 Scale est. = 1 n = 58

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 @@ -248,17 +246,17 @@

A single outcome variable levels = 'series1'), time = 1:10) gauss_dat -#> outcome series time -#> 1 0.1754683 series1 1 -#> 2 0.2411510 series1 2 -#> 3 1.0221725 series1 3 -#> 4 0.7186505 series1 4 -#> 5 0.4864740 series1 5 -#> 6 1.6725269 series1 6 -#> 7 -0.2081990 series1 7 -#> 8 0.2345187 series1 8 -#> 9 -0.8994249 series1 9 -#> 10 -0.8554195 series1 10 +#> outcome series time +#> 1 -1.105964948 series1 1 +#> 2 0.299966201 series1 2 +#> 3 -0.942046180 series1 3 +#> 4 -0.001887559 series1 4 +#> 5 0.339117564 series1 5 +#> 6 1.458563643 series1 6 +#> 7 -0.024058830 series1 7 +#> 8 0.646390934 series1 8 +#> 9 -0.686607550 series1 9 +#> 10 0.126273947 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):

@@ -269,14 +267,14 @@

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

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

@@ -369,15 +367,15 @@ 

Checking data with get_mvgam_ series = factor('series_1'), outcome = rnorm(8)) bad_times -#> time series outcome -#> 1 1 series_1 -0.8227444 -#> 2 3 series_1 -0.6851957 -#> 3 5 series_1 1.5221436 -#> 4 7 series_1 -0.8227729 -#> 5 9 series_1 1.0787734 -#> 6 11 series_1 -0.7118972 -#> 7 13 series_1 -1.2150116 -#> 8 15 series_1 1.5165404

+#> time series outcome +#> 1 1 series_1 1.48518852 +#> 2 3 series_1 -0.56975888 +#> 3 5 series_1 -0.63240452 +#> 4 7 series_1 -2.30914822 +#> 5 9 series_1 -0.03495886 +#> 6 11 series_1 0.54345770 +#> 7 13 series_1 0.99494365 +#> 8 15 series_1 -0.35716657

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

@@ -398,22 +396,22 @@ 

Checking data with get_mvgam_ dplyr::arrange(time) -> good_times #> Joining with `by = join_by(time, series)` good_times -#> time series outcome -#> 1 1 series_1 -0.8227444 -#> 2 2 series_1 NA -#> 3 3 series_1 -0.6851957 -#> 4 4 series_1 NA -#> 5 5 series_1 1.5221436 -#> 6 6 series_1 NA -#> 7 7 series_1 -0.8227729 -#> 8 8 series_1 NA -#> 9 9 series_1 1.0787734 -#> 10 10 series_1 NA -#> 11 11 series_1 -0.7118972 -#> 12 12 series_1 NA -#> 13 13 series_1 -1.2150116 -#> 14 14 series_1 NA -#> 15 15 series_1 1.5165404

+#> time series outcome +#> 1 1 series_1 1.48518852 +#> 2 2 series_1 NA +#> 3 3 series_1 -0.56975888 +#> 4 4 series_1 NA +#> 5 5 series_1 -0.63240452 +#> 6 6 series_1 NA +#> 7 7 series_1 -2.30914822 +#> 8 8 series_1 NA +#> 9 9 series_1 -0.03495886 +#> 10 10 series_1 NA +#> 11 11 series_1 0.54345770 +#> 12 12 series_1 NA +#> 13 13 series_1 0.99494365 +#> 14 14 series_1 NA +#> 15 15 series_1 -0.35716657

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

@@ -423,9 +421,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.91, 0.9); +#> 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.47, 0.56); #> new_lowerbound new_upperbound #> 1 NA NA #> 2 NA NA

@@ -473,9 +471,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.5, 2.5); (Intercept) ~ normal(0, 1); -#> 2 sigma_obs ~ student_t(3, 0, 2.5); sigma_obs ~ normal(-0.01, 0.79); +#> prior example_change +#> 1 (Intercept) ~ student_t(3, -0.4, 2.5); (Intercept) ~ normal(0, 1); +#> 2 sigma_obs ~ student_t(3, 0, 2.5); sigma_obs ~ normal(-0.62, 1); #> new_lowerbound new_upperbound #> 1 NA NA #> 2 NA NA @@ -498,17 +496,17 @@

Covariates with no NAs= 'series1'), time = 1:10) miss_dat -#> outcome cov series time -#> 1 -1.5062280 NA series1 1 -#> 2 0.8873110 -1.419737577 series1 2 -#> 3 -0.5854087 0.006758633 series1 3 -#> 4 -0.8946378 -0.989167769 series1 4 -#> 5 -0.4520713 2.359401432 series1 5 -#> 6 0.5739319 -0.802900274 series1 6 -#> 7 -0.5756935 -0.735234481 series1 7 -#> 8 0.4873677 -1.386964284 series1 8 -#> 9 1.2759287 -0.474650053 series1 9 -#> 10 -1.8709314 -1.539955160 series1 10 +#> outcome cov series time +#> 1 0.1561020 NA series1 1 +#> 2 0.1084784 0.98322738 series1 2 +#> 3 0.0330442 0.28407529 series1 3 +#> 4 0.6572311 -1.45353685 series1 4 +#> 5 0.8854784 0.02783985 series1 5 +#> 6 1.6044159 -1.17292983 series1 6 +#> 7 -0.9053955 0.88927692 series1 7 +#> 8 2.3019384 0.36575792 series1 8 +#> 9 0.1608429 -0.58446531 series1 9 +#> 10 0.9125095 3.75742498 series1 10
 get_mvgam_priors(outcome ~ cov,
                  data = miss_dat,
@@ -517,10 +515,10 @@ 

Covariates with no NAs#> 1 (Intercept) 1 (Intercept) #> 2 cov 1 cov fixed effect #> 3 vector<lower=0>[n_series] sigma_obs; 1 observation error sd -#> prior example_change -#> 1 (Intercept) ~ student_t(3, -0.5, 2.5); (Intercept) ~ normal(0, 1); -#> 2 cov ~ student_t(3, 0, 2); cov ~ normal(0, 1); -#> 3 sigma_obs ~ student_t(3, 0, 2.5); sigma_obs ~ normal(0.47, 0.76); +#> prior example_change +#> 1 (Intercept) ~ student_t(3, 0.4, 2.5); (Intercept) ~ normal(0, 1); +#> 2 cov ~ student_t(3, 0, 2); cov ~ normal(0, 1); +#> 3 sigma_obs ~ student_t(3, 0, 2.5); sigma_obs ~ normal(-0.59, 0.3); #> new_lowerbound new_upperbound #> 1 NA NA #> 2 NA NA @@ -551,14 +549,14 @@

Covariates with no NAs#> 5 cov4 1 cov4 fixed effect #> 6 cov5 1 cov5 fixed effect #> 7 vector<lower=0>[n_series] sigma_obs; 1 observation error sd -#> prior example_change -#> 1 (Intercept) ~ student_t(3, -0.4, 2.5); (Intercept) ~ normal(0, 1); -#> 2 cov1 ~ student_t(3, 0, 2); cov1 ~ normal(0, 1); -#> 3 cov2 ~ student_t(3, 0, 2); cov2 ~ normal(0, 1); -#> 4 cov3 ~ student_t(3, 0, 2); cov3 ~ normal(0, 1); -#> 5 cov4 ~ student_t(3, 0, 2); cov4 ~ normal(0, 1); -#> 6 cov5 ~ student_t(3, 0, 2); cov5 ~ normal(0, 1); -#> 7 sigma_obs ~ student_t(3, 0, 2.5); sigma_obs ~ normal(-0.4, 0.32); +#> prior example_change +#> 1 (Intercept) ~ student_t(3, -0.2, 2.5); (Intercept) ~ normal(0, 1); +#> 2 cov1 ~ student_t(3, 0, 2); cov1 ~ normal(0, 1); +#> 3 cov2 ~ student_t(3, 0, 2); cov2 ~ normal(0, 1); +#> 4 cov3 ~ student_t(3, 0, 2); cov3 ~ normal(0, 1); +#> 5 cov4 ~ student_t(3, 0, 2); cov4 ~ normal(0, 1); +#> 6 cov5 ~ student_t(3, 0, 2); cov5 ~ normal(0, 1); +#> 7 sigma_obs ~ student_t(3, 0, 2.5); sigma_obs ~ normal(-0.32, 0.68); #> new_lowerbound new_upperbound #> 1 NA NA #> 2 NA NA @@ -740,8 +738,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',
@@ -770,7 +768,7 @@ 

Example with NEON tick data#> $ S8 : num [1:8, 1:8] 1.037 -0.416 0.419 0.117 0.188 ... #> $ p_coefs : Named num 0 #> ..- attr(*, "names")= chr "(Intercept)" -#> $ p_taus : num 1.14 +#> $ p_taus : num 0.945 #> $ 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] 0.368 0.368 0.368 0.368 0.368 ... 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 327ff76c..a3825320 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 d8288e81..0c08f9cd 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/data_in_mvgam_files/figure-html/unnamed-chunk-26-1.png b/docs/articles/data_in_mvgam_files/figure-html/unnamed-chunk-26-1.png index 3c39ede9..cb03d0e4 100644 Binary files a/docs/articles/data_in_mvgam_files/figure-html/unnamed-chunk-26-1.png and b/docs/articles/data_in_mvgam_files/figure-html/unnamed-chunk-26-1.png differ diff --git a/docs/articles/forecast_evaluation.html b/docs/articles/forecast_evaluation.html index 804ac365..be3096b2 100644 --- a/docs/articles/forecast_evaluation.html +++ b/docs/articles/forecast_evaluation.html @@ -26,7 +26,7 @@ mvgam - 1.1.2 + 1.1.3

Plotting GAM smooth functions using mvgamPlotting GAM smooth functions using mvgam

@@ -233,12 +233,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',
@@ -275,32 +275,32 @@ 

Modelling dynamics with GPs#> #> GAM coefficient (beta) estimates: #> 2.5% 50% 97.5% Rhat n_eff -#> (Intercept) -1.1 -0.51 0.34 1 694 +#> (Intercept) -1.1 -0.51 0.25 1 731 #> #> GAM gp term marginal deviation (alpha) and length scale (rho) estimates: #> 2.5% 50% 97.5% Rhat n_eff -#> alpha_gp(time):seriesseries_1 0.21 0.78 2.2 1 819 -#> alpha_gp(time):seriesseries_2 0.73 1.30 2.9 1 761 -#> alpha_gp(time):seriesseries_3 0.46 1.10 2.8 1 1262 -#> rho_gp(time):seriesseries_1 1.30 5.40 22.0 1 682 -#> rho_gp(time):seriesseries_2 2.10 10.00 17.0 1 450 -#> rho_gp(time):seriesseries_3 1.50 8.80 24.0 1 769 +#> alpha_gp(time):seriesseries_1 0.19 0.75 2 1.00 932 +#> alpha_gp(time):seriesseries_2 0.76 1.40 3 1.00 845 +#> alpha_gp(time):seriesseries_3 0.48 1.10 3 1.00 968 +#> rho_gp(time):seriesseries_1 1.10 4.90 25 1.01 678 +#> rho_gp(time):seriesseries_2 2.20 10.00 17 1.00 645 +#> rho_gp(time):seriesseries_3 1.50 9.20 24 1.00 908 #> #> Approximate significance of GAM smooths: -#> edf Ref.df Chi.sq p-value -#> s(season) 3.36 6 21.1 0.0093 ** +#> edf Ref.df Chi.sq p-value +#> s(season) 3.4 6 21.1 0.011 * #> --- #> 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%) +#> 7 of 2000 iterations ended with a divergence (0.35%) #> *Try running with larger adapt_delta to remove the divergences #> 0 of 2000 iterations saturated the maximum tree depth of 12 (0%) #> E-FMI indicated no pathological behavior #> -#> Samples were drawn using NUTS(diag_e) at Sat Jun 29 12:35:39 PM 2024. +#> Samples were drawn using NUTS(diag_e) at Tue Sep 03 2:33:39 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)

@@ -309,16 +309,16 @@

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')

Summarising latent Gaussian Process parameters in mvgam

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')

Summarising latent Gaussian Process parameters in mvgam

We can again plot the nonlinear effects:

-conditional_effects(mod2, type = 'link')
+conditional_effects(mod2, type = 'link')

Plotting latent Gaussian Process effects in mvgam and marginaleffectsPlotting latent Gaussian Process effects in mvgam and marginaleffects

The estimates for the temporal trends are fairly similar for the two models, but below we will see if they produce similar forecasts

@@ -374,45 +374,45 @@

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

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

 plot(fc_mod1, series = 1)
 #> Out of sample CRPS:
-#> 14.89051875
+#> 14.58308175

 plot(fc_mod2, series = 1)
 #> Out of sample DRPS:
-#> 10.84228725
+#> 10.7328925

 
 plot(fc_mod1, series = 2)
 #> Out of sample CRPS:
-#> 495050222726067
+#> 1.74633239486101e+20

 plot(fc_mod2, series = 2)
-#> Out of sample CRPS:
-#> 14.7121945
+#> Out of sample DRPS: +#> 14.23353775

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

@@ -429,8 +429,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',
@@ -447,7 +447,7 @@ 

Forecasting with newdata
 plot(fc_mod2, series = 1)
 #> Out of sample DRPS:
-#> 10.78167525
+#> 10.8005245

Plotting posterior forecast distributions using mvgam and R

@@ -465,54 +465,54 @@

Scoring forecast distributionsstr(crps_mod1) #> List of 4 #> $ series_1 :'data.frame': 25 obs. of 5 variables: -#> ..$ score : num [1:25] 0.186 0.129 1.372 NA 0.037 ... +#> ..$ score : num [1:25] 0.1797 0.1341 1.3675 NA 0.0386 ... #> ..$ in_interval : num [1:25] 1 1 1 NA 1 1 1 1 1 1 ... #> ..$ interval_width: num [1:25] 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 ... #> ..$ eval_horizon : int [1:25] 1 2 3 4 5 6 7 8 9 10 ... #> ..$ score_type : chr [1:25] "crps" "crps" "crps" "crps" ... #> $ series_2 :'data.frame': 25 obs. of 5 variables: -#> ..$ score : num [1:25] 0.354 0.334 0.947 0.492 0.542 ... +#> ..$ score : num [1:25] 0.375 0.283 1.003 0.516 0.649 ... #> ..$ in_interval : num [1:25] 1 1 1 1 1 1 1 1 1 NA ... #> ..$ interval_width: num [1:25] 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 ... #> ..$ eval_horizon : int [1:25] 1 2 3 4 5 6 7 8 9 10 ... #> ..$ score_type : chr [1:25] "crps" "crps" "crps" "crps" ... #> $ series_3 :'data.frame': 25 obs. of 5 variables: -#> ..$ score : num [1:25] 0.31 0.616 0.4 0.349 0.215 ... +#> ..$ score : num [1:25] 0.318 0.604 0.4 0.353 0.212 ... #> ..$ in_interval : num [1:25] 1 1 1 1 1 1 1 1 1 1 ... #> ..$ interval_width: num [1:25] 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 ... #> ..$ eval_horizon : int [1:25] 1 2 3 4 5 6 7 8 9 10 ... #> ..$ score_type : chr [1:25] "crps" "crps" "crps" "crps" ... #> $ all_series:'data.frame': 25 obs. of 3 variables: -#> ..$ score : num [1:25] 0.85 1.079 2.719 NA 0.794 ... +#> ..$ score : num [1:25] 0.873 1.021 2.77 NA 0.9 ... #> ..$ eval_horizon: int [1:25] 1 2 3 4 5 6 7 8 9 10 ... #> ..$ score_type : chr [1:25] "sum_crps" "sum_crps" "sum_crps" "sum_crps" ... crps_mod1$series_1 #> score in_interval interval_width eval_horizon score_type -#> 1 0.18582425 1 0.9 1 crps -#> 2 0.12933350 1 0.9 2 crps -#> 3 1.37181050 1 0.9 3 crps +#> 1 0.17965150 1 0.9 1 crps +#> 2 0.13411725 1 0.9 2 crps +#> 3 1.36749550 1 0.9 3 crps #> 4 NA NA 0.9 4 crps -#> 5 0.03698600 1 0.9 5 crps -#> 6 1.53997900 1 0.9 6 crps -#> 7 1.50467675 1 0.9 7 crps -#> 8 0.63460725 1 0.9 8 crps -#> 9 0.61682725 1 0.9 9 crps -#> 10 0.62428875 1 0.9 10 crps -#> 11 1.33824700 1 0.9 11 crps -#> 12 2.06378300 1 0.9 12 crps -#> 13 0.59247200 1 0.9 13 crps -#> 14 0.13560025 1 0.9 14 crps -#> 15 0.66512975 1 0.9 15 crps -#> 16 0.08238525 1 0.9 16 crps -#> 17 0.08152900 1 0.9 17 crps -#> 18 0.09446425 1 0.9 18 crps -#> 19 0.12084700 1 0.9 19 crps +#> 5 0.03860225 1 0.9 5 crps +#> 6 1.55334350 1 0.9 6 crps +#> 7 1.49198325 1 0.9 7 crps +#> 8 0.64088650 1 0.9 8 crps +#> 9 0.61613650 1 0.9 9 crps +#> 10 0.60528025 1 0.9 10 crps +#> 11 1.30624075 1 0.9 11 crps +#> 12 2.06809025 1 0.9 12 crps +#> 13 0.61887550 1 0.9 13 crps +#> 14 0.13920225 1 0.9 14 crps +#> 15 0.67819725 1 0.9 15 crps +#> 16 0.07817500 1 0.9 16 crps +#> 17 0.07567500 1 0.9 17 crps +#> 18 0.09510025 1 0.9 18 crps +#> 19 0.12604375 1 0.9 19 crps #> 20 NA NA 0.9 20 crps -#> 21 0.21286925 1 0.9 21 crps -#> 22 0.85799700 1 0.9 22 crps +#> 21 0.20347825 1 0.9 21 crps +#> 22 0.82202975 1 0.9 22 crps #> 23 NA NA 0.9 23 crps -#> 24 1.14954750 1 0.9 24 crps -#> 25 0.85131425 1 0.9 25 crps

+#> 24 1.06660275 1 0.9 24 crps +#> 25 0.67787450 1 0.9 25 crps

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

Scoring forecast distributionscrps_mod1 <- score(fc_mod1, score = 'crps', interval_width = 0.6) crps_mod1$series_1 #> score in_interval interval_width eval_horizon score_type -#> 1 0.18582425 1 0.6 1 crps -#> 2 0.12933350 1 0.6 2 crps -#> 3 1.37181050 0 0.6 3 crps +#> 1 0.17965150 1 0.6 1 crps +#> 2 0.13411725 1 0.6 2 crps +#> 3 1.36749550 0 0.6 3 crps #> 4 NA NA 0.6 4 crps -#> 5 0.03698600 1 0.6 5 crps -#> 6 1.53997900 0 0.6 6 crps -#> 7 1.50467675 0 0.6 7 crps -#> 8 0.63460725 1 0.6 8 crps -#> 9 0.61682725 1 0.6 9 crps -#> 10 0.62428875 1 0.6 10 crps -#> 11 1.33824700 0 0.6 11 crps -#> 12 2.06378300 0 0.6 12 crps -#> 13 0.59247200 1 0.6 13 crps -#> 14 0.13560025 1 0.6 14 crps -#> 15 0.66512975 1 0.6 15 crps -#> 16 0.08238525 1 0.6 16 crps -#> 17 0.08152900 1 0.6 17 crps -#> 18 0.09446425 1 0.6 18 crps -#> 19 0.12084700 1 0.6 19 crps +#> 5 0.03860225 1 0.6 5 crps +#> 6 1.55334350 0 0.6 6 crps +#> 7 1.49198325 0 0.6 7 crps +#> 8 0.64088650 1 0.6 8 crps +#> 9 0.61613650 1 0.6 9 crps +#> 10 0.60528025 1 0.6 10 crps +#> 11 1.30624075 0 0.6 11 crps +#> 12 2.06809025 0 0.6 12 crps +#> 13 0.61887550 1 0.6 13 crps +#> 14 0.13920225 1 0.6 14 crps +#> 15 0.67819725 1 0.6 15 crps +#> 16 0.07817500 1 0.6 16 crps +#> 17 0.07567500 1 0.6 17 crps +#> 18 0.09510025 1 0.6 18 crps +#> 19 0.12604375 1 0.6 19 crps #> 20 NA NA 0.6 20 crps -#> 21 0.21286925 1 0.6 21 crps -#> 22 0.85799700 1 0.6 22 crps +#> 21 0.20347825 1 0.6 21 crps +#> 22 0.82202975 1 0.6 22 crps #> 23 NA NA 0.6 23 crps -#> 24 1.14954750 1 0.6 24 crps -#> 25 0.85131425 1 0.6 25 crps +#> 24 1.06660275 1 0.6 24 crps +#> 25 0.67787450 1 0.6 25 crps

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

Scoring forecast distributionslink_mod1 <- forecast(mod1, newdata = simdat$data_test, type = 'link') score(link_mod1, score = 'elpd')$series_1 #> score eval_horizon score_type -#> 1 -0.5343784 1 elpd -#> 2 -0.4326190 2 elpd -#> 3 -2.9699450 3 elpd +#> 1 -0.5285206 1 elpd +#> 2 -0.4286994 2 elpd +#> 3 -2.9660940 3 elpd #> 4 NA 4 elpd -#> 5 -0.1998425 5 elpd -#> 6 -3.3976729 6 elpd -#> 7 -3.2989297 7 elpd -#> 8 -2.0490633 8 elpd -#> 9 -2.0690163 9 elpd -#> 10 -2.0822051 10 elpd -#> 11 -3.1101639 11 elpd -#> 12 -3.7240924 12 elpd -#> 13 -2.1578701 13 elpd -#> 14 -0.2899481 14 elpd -#> 15 -2.3811862 15 elpd -#> 16 -0.2085375 16 elpd -#> 17 -0.1960501 17 elpd -#> 18 -0.2036978 18 elpd -#> 19 -0.2154374 19 elpd +#> 5 -0.1988847 5 elpd +#> 6 -3.3821055 6 elpd +#> 7 -3.2797236 7 elpd +#> 8 -2.0571076 8 elpd +#> 9 -2.0794559 9 elpd +#> 10 -2.0882202 10 elpd +#> 11 -3.0870256 11 elpd +#> 12 -3.7065927 12 elpd +#> 13 -2.1601960 13 elpd +#> 14 -0.2931143 14 elpd +#> 15 -2.3694878 15 elpd +#> 16 -0.2110566 16 elpd +#> 17 -0.1986209 17 elpd +#> 18 -0.2069720 18 elpd +#> 19 -0.2193413 19 elpd #> 20 NA 20 elpd -#> 21 -0.2341597 21 elpd -#> 22 -2.6552948 22 elpd +#> 21 -0.2379762 21 elpd +#> 22 -2.6261457 22 elpd #> 23 NA 23 elpd -#> 24 -2.6652717 24 elpd -#> 25 -0.2759126 25 elpd +#> 24 -2.6309438 24 elpd +#> 25 -0.2817444 25 elpd

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

Scoring forecast distributions#> ..$ interval_width: num [1:25] 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 ... #> ..$ eval_horizon : int [1:25] 1 2 3 4 5 6 7 8 9 10 ... #> $ all_series:'data.frame': 25 obs. of 3 variables: -#> ..$ score : num [1:25] 0.771 1.133 1.26 NA 0.443 ... +#> ..$ score : num [1:25] 0.755 1.12 1.245 NA 0.447 ... #> ..$ eval_horizon: int [1:25] 1 2 3 4 5 6 7 8 9 10 ... #> ..$ score_type : chr [1:25] "energy" "energy" "energy" "energy" ...

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

Scoring forecast distributions
 energy_mod2$all_series
 #>        score eval_horizon score_type
-#> 1  0.7705198            1     energy
-#> 2  1.1330328            2     energy
-#> 3  1.2600785            3     energy
+#> 1  0.7546579            1     energy
+#> 2  1.1200630            2     energy
+#> 3  1.2447843            3     energy
 #> 4         NA            4     energy
-#> 5  0.4427578            5     energy
-#> 6  1.8848308            6     energy
-#> 7  1.4186997            7     energy
-#> 8  0.7280518            8     energy
-#> 9  1.0467755            9     energy
+#> 5  0.4465348            5     energy
+#> 6  1.8231460            6     energy
+#> 7  1.4418019            7     energy
+#> 8  0.7172890            8     energy
+#> 9  1.0762943            9     energy
 #> 10        NA           10     energy
-#> 11 1.4172423           11     energy
-#> 12 3.2326925           12     energy
-#> 13 1.5987732           13     energy
-#> 14 1.1798872           14     energy
-#> 15 1.0311968           15     energy
-#> 16 1.8261356           16     energy
+#> 11 1.4112423           11     energy
+#> 12 3.2385416           12     energy
+#> 13 1.5836460           13     energy
+#> 14 1.1953349           14     energy
+#> 15 1.0412578           15     energy
+#> 16 1.8348615           16     energy
 #> 17        NA           17     energy
-#> 18 0.7170961           18     energy
-#> 19 0.8927311           19     energy
+#> 18 0.7142977           18     energy
+#> 19 0.9059773           19     energy
 #> 20        NA           20     energy
-#> 21 1.0544501           21     energy
-#> 22 1.3280321           22     energy
+#> 21 1.1043397           21     energy
+#> 22 1.3292391           22     energy
 #> 23        NA           23     energy
-#> 24 2.1843621           24     energy
-#> 25 1.2352041           25     energy
+#> 24 2.1419570 24 energy +#> 25 1.2610880 25 energy

You can use your score(s) of choice to compare different models. For example, we can compute and plot the difference in CRPS scores for each series in data. Here, a negative value means the Gaussian Process model diff --git a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-13-1.png b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-13-1.png index a928b1f0..dba7cdfb 100644 Binary files a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-13-1.png and b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-13-1.png differ diff --git a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-14-1.png b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-14-1.png index 91929bec..c4c70376 100644 Binary files a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-14-1.png and b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-14-1.png differ diff --git a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-15-1.png b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-15-1.png index c68f404f..f81bd3e5 100644 Binary files a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-15-1.png and b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-15-1.png differ diff --git a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-15-2.png b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-15-2.png index 1da70679..bf48bc83 100644 Binary files a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-15-2.png and b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-15-2.png differ diff --git a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-18-1.png b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-18-1.png index 0199fa02..bce5f434 100644 Binary files a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-18-1.png and b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-18-1.png differ diff --git a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-18-2.png b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-18-2.png index 29b6a888..f08e7dcb 100644 Binary files a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-18-2.png and b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-18-2.png differ diff --git a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-18-3.png b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-18-3.png index ad12aeb3..dcc6592a 100644 Binary files a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-18-3.png and b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-18-3.png differ diff --git a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-18-4.png b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-18-4.png index 2ae4de6d..61ac74bf 100644 Binary files a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-18-4.png and b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-18-4.png differ diff --git a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-22-1.png b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-22-1.png index f5f5a83c..8347eca5 100644 Binary files a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-22-1.png and b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-22-1.png differ diff --git a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-28-1.png b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-28-1.png index 5643424b..46bb5f6e 100644 Binary files a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-28-1.png and b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-28-1.png differ diff --git a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-28-2.png b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-28-2.png index d73a0f58..51002835 100644 Binary files a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-28-2.png and b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-28-2.png differ diff --git a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-28-3.png b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-28-3.png index 365642af..81d4de53 100644 Binary files a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-28-3.png and b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-28-3.png differ diff --git a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-9-1.png b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-9-1.png index 146519aa..3e22fb75 100644 Binary files a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-9-1.png and b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-9-1.png differ diff --git a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-9-2.png b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-9-2.png index 822eea46..275ed111 100644 Binary files a/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-9-2.png and b/docs/articles/forecast_evaluation_files/figure-html/unnamed-chunk-9-2.png differ diff --git a/docs/articles/index.html b/docs/articles/index.html index 9bfbcb90..7de0d62e 100644 --- a/docs/articles/index.html +++ b/docs/articles/index.html @@ -10,7 +10,7 @@ mvgam - 1.0.91 + 1.1.3 @@ -62,13 +62,13 @@

- +
@@ -80,15 +80,15 @@

State-Space models in mvgam

Nicholas J Clark

- -

2024-07-05

- + +

2024-09-03

+ Source: vignettes/trend_formulas.Rmd
trend_formulas.Rmd
- - + +

The purpose of this vignette is to show how the mvgam package can be used to fit and interrogate State-Space models with nonlinear effects.

@@ -96,7 +96,7 @@

2024-07-05

State-Space Models

-Illustration of a basic State-Space model, which assumes that a latent dynamic process (X) can evolve independently from the way we take observations (Y) of that process
Illustration of a basic State-Space model, which +Illustration of a basic State-Space model, which assumes that a latent dynamic process (X) can evolve independently from the way we take observations (Y) of that process
Illustration of a basic State-Space model, which assumes that a latent dynamic process (X) can evolve independently from the way we take observations (Y) of that process
@@ -198,15 +198,15 @@

Lake Washington plankton data
 plankton_data %>%
   dplyr::filter(series == 'Other.algae') %>%
-  ggplot(aes(x = time, y = temp)) +
-  geom_line(size = 1.1) +
-  geom_line(aes(y = y), col = 'white',
+  ggplot(aes(x = time, y = temp)) +
+  geom_line(size = 1.1) +
+  geom_line(aes(y = y), col = 'white',
             size = 1.3) +
-  geom_line(aes(y = y), col = 'darkred',
+  geom_line(aes(y = y), col = 'darkred',
             size = 1.1) +
-  ylab('z-score') +
-  xlab('Time') +
-  ggtitle('Temperature (black) vs Other algae (red)')
+  ylab('z-score') +
+  xlab('Time') +
+  ggtitle('Temperature (black) vs Other algae (red)')
 #> Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
 #>  Please use `linewidth` instead.
 #> This warning is displayed once every 8 hours.
@@ -216,15 +216,15 @@ 

Lake Washington plankton data
 plankton_data %>%
   dplyr::filter(series == 'Diatoms') %>%
-  ggplot(aes(x = time, y = temp)) +
-  geom_line(size = 1.1) +
-  geom_line(aes(y = y), col = 'white',
+  ggplot(aes(x = time, y = temp)) +
+  geom_line(size = 1.1) +
+  geom_line(aes(y = y), col = 'white',
             size = 1.3) +
-  geom_line(aes(y = y), col = 'darkred',
+  geom_line(aes(y = y), col = 'darkred',
             size = 1.1) +
-  ylab('z-score') +
-  xlab('Time') +
-  ggtitle('Temperature (black) vs Diatoms (red)')

+ ylab('z-score') + + xlab('Time') + + ggtitle('Temperature (black) vs Diatoms (red)')

We will have to try and capture this seasonality in our process model, which should be easy to do given the flexibility of GAMs. Next we @@ -266,10 +266,10 @@

Capturing seasonalitynotrend_mod <- mvgam(y ~ # tensor of temp and month to capture # "global" seasonality - te(temp, month, k = c(4, 4)) + + te(temp, month, k = c(4, 4)) + # series-specific deviation tensor products - te(temp, month, k = c(4, 4), by = series) - 1, + te(temp, month, k = c(4, 4), by = series) - 1, family = gaussian(), data = plankton_train, newdata = plankton_test, @@ -293,17 +293,17 @@

Capturing seasonality
 plot(notrend_mod, type = 'forecast', series = 1)
 #> Out of sample CRPS:
-#> 7.07436462451194
+#> 7.01328371945431

 plot(notrend_mod, type = 'forecast', series = 2)
 #> Out of sample CRPS:
-#> 6.60945059665685
+#> 6.54071425133111

 plot(notrend_mod, type = 'forecast', series = 3)
 #> Out of sample CRPS:
-#> 4.05305783104855
+#> 4.06261023205356

This basic model gives us confidence that we can capture the seasonal variation in the observations. But the model has not captured the @@ -387,8 +387,8 @@

Multiseries dynamicsy ~ -1, # process model formula, which includes the smooth functions - trend_formula = ~ te(temp, month, k = c(4, 4)) + - te(temp, month, k = c(4, 4), by = trend) - 1, + trend_formula = ~ te(temp, month, k = c(4, 4)) + + te(temp, month, k = c(4, 4), by = trend) - 1, # VAR1 model with uncorrelated process errors trend_model = VAR(), @@ -429,8 +429,8 @@

Multiseries dynamics
-priors <- c(prior(normal(0.5, 0.1), class = sigma_obs, lb = 0.2),
-            prior(normal(0.5, 0.25), class = sigma))
+priors <- c(prior(normal(0.5, 0.1), class = sigma_obs, lb = 0.2), + prior(normal(0.5, 0.25), class = sigma))

You may have noticed something else unique about this model: there is no intercept term in the observation formula. This is because a shared intercept parameter can sometimes be unidentifiable with respect to the @@ -446,8 +446,8 @@

Multiseries dynamicsy ~ -1, # process model formula, which includes the smooth functions - trend_formula = ~ te(temp, month, k = c(4, 4)) + - te(temp, month, k = c(4, 4), by = trend) - 1, + trend_formula = ~ te(temp, month, k = c(4, 4)) + + te(temp, month, k = c(4, 4), by = trend) - 1, # VAR1 model with uncorrelated process errors trend_model = VAR(), @@ -508,90 +508,90 @@

Inspecting SS models#> #> Observation error parameter estimates: #> 2.5% 50% 97.5% Rhat n_eff -#> sigma_obs[1] 0.20 0.26 0.34 1.00 550 -#> sigma_obs[2] 0.26 0.40 0.54 1.01 198 -#> sigma_obs[3] 0.41 0.63 0.81 1.07 37 -#> sigma_obs[4] 0.25 0.37 0.50 1.01 310 -#> sigma_obs[5] 0.30 0.43 0.54 1.02 309 +#> sigma_obs[1] 0.20 0.26 0.34 1.01 381 +#> sigma_obs[2] 0.26 0.40 0.54 1.05 143 +#> sigma_obs[3] 0.42 0.65 0.82 1.13 41 +#> sigma_obs[4] 0.25 0.38 0.49 1.01 248 +#> sigma_obs[5] 0.31 0.43 0.55 1.02 265 #> #> GAM observation model coefficient (beta) estimates: #> 2.5% 50% 97.5% Rhat n_eff #> (Intercept) 0 0 0 NaN NaN #> #> Process model VAR parameter estimates: -#> 2.5% 50% 97.5% Rhat n_eff -#> A[1,1] -0.012 0.480 0.850 1.07 40 -#> A[1,2] -0.350 -0.025 0.190 1.01 373 -#> A[1,3] -0.460 -0.036 0.360 1.01 446 -#> A[1,4] -0.270 0.031 0.410 1.01 452 -#> A[1,5] -0.099 0.140 0.550 1.03 107 -#> A[2,1] -0.140 0.017 0.200 1.00 825 -#> A[2,2] 0.630 0.800 0.930 1.01 254 -#> A[2,3] -0.380 -0.120 0.044 1.02 346 -#> A[2,4] -0.041 0.099 0.320 1.01 328 -#> A[2,5] -0.058 0.056 0.200 1.00 694 -#> A[3,1] -0.260 0.022 0.470 1.10 32 -#> A[3,2] -0.480 -0.170 0.036 1.01 178 -#> A[3,3] 0.100 0.450 0.760 1.02 181 -#> A[3,4] -0.040 0.200 0.590 1.02 185 -#> A[3,5] -0.091 0.110 0.370 1.04 110 -#> A[4,1] -0.150 0.056 0.340 1.03 257 -#> A[4,2] -0.099 0.057 0.250 1.01 558 -#> A[4,3] -0.400 -0.110 0.120 1.01 397 -#> A[4,4] 0.490 0.740 0.940 1.02 399 -#> A[4,5] -0.190 -0.033 0.120 1.00 821 -#> A[5,1] -0.200 0.077 0.640 1.05 63 -#> A[5,2] -0.380 -0.110 0.078 1.02 239 -#> A[5,3] -0.550 -0.160 0.140 1.01 268 -#> A[5,4] -0.065 0.170 0.540 1.01 270 -#> A[5,5] 0.500 0.720 0.920 1.01 350 +#> 2.5% 50% 97.5% Rhat n_eff +#> A[1,1] -0.0022 0.520 0.850 1.10 48 +#> A[1,2] -0.3700 -0.037 0.200 1.00 470 +#> A[1,3] -0.4900 -0.052 0.310 1.01 590 +#> A[1,4] -0.2600 0.037 0.430 1.00 436 +#> A[1,5] -0.0980 0.130 0.510 1.03 183 +#> A[2,1] -0.1600 0.015 0.220 1.00 718 +#> A[2,2] 0.6100 0.780 0.910 1.01 304 +#> A[2,3] -0.4200 -0.150 0.030 1.01 253 +#> A[2,4] -0.0340 0.120 0.390 1.02 201 +#> A[2,5] -0.0560 0.066 0.210 1.00 749 +#> A[3,1] -0.3000 0.021 0.490 1.06 81 +#> A[3,2] -0.5400 -0.220 0.016 1.01 190 +#> A[3,3] 0.0660 0.410 0.700 1.02 305 +#> A[3,4] -0.0120 0.250 0.670 1.02 177 +#> A[3,5] -0.0680 0.130 0.400 1.03 209 +#> A[4,1] -0.1600 0.059 0.340 1.04 152 +#> A[4,2] -0.1100 0.057 0.250 1.01 384 +#> A[4,3] -0.4400 -0.120 0.120 1.01 334 +#> A[4,4] 0.4800 0.740 0.970 1.02 212 +#> A[4,5] -0.2000 -0.037 0.120 1.00 861 +#> A[5,1] -0.2000 0.085 0.600 1.05 88 +#> A[5,2] -0.4800 -0.150 0.061 1.02 137 +#> A[5,3] -0.6900 -0.210 0.098 1.01 168 +#> A[5,4] -0.0460 0.210 0.730 1.02 142 +#> A[5,5] 0.5300 0.750 0.960 1.00 449 #> #> Process error parameter estimates: #> 2.5% 50% 97.5% Rhat n_eff -#> Sigma[1,1] 0.053 0.30 0.69 1.10 28 +#> Sigma[1,1] 0.038 0.25 0.65 1.18 31 #> Sigma[1,2] 0.000 0.00 0.00 NaN NaN #> Sigma[1,3] 0.000 0.00 0.00 NaN NaN #> Sigma[1,4] 0.000 0.00 0.00 NaN NaN #> Sigma[1,5] 0.000 0.00 0.00 NaN NaN #> Sigma[2,1] 0.000 0.00 0.00 NaN NaN -#> Sigma[2,2] 0.066 0.11 0.18 1.01 601 +#> Sigma[2,2] 0.065 0.11 0.18 1.02 419 #> Sigma[2,3] 0.000 0.00 0.00 NaN NaN #> Sigma[2,4] 0.000 0.00 0.00 NaN NaN #> Sigma[2,5] 0.000 0.00 0.00 NaN NaN #> Sigma[3,1] 0.000 0.00 0.00 NaN NaN #> Sigma[3,2] 0.000 0.00 0.00 NaN NaN -#> Sigma[3,3] 0.059 0.16 0.29 1.03 161 +#> Sigma[3,3] 0.054 0.15 0.30 1.07 74 #> Sigma[3,4] 0.000 0.00 0.00 NaN NaN #> Sigma[3,5] 0.000 0.00 0.00 NaN NaN #> Sigma[4,1] 0.000 0.00 0.00 NaN NaN #> Sigma[4,2] 0.000 0.00 0.00 NaN NaN #> Sigma[4,3] 0.000 0.00 0.00 NaN NaN -#> Sigma[4,4] 0.056 0.13 0.27 1.05 185 +#> Sigma[4,4] 0.043 0.13 0.25 1.04 115 #> Sigma[4,5] 0.000 0.00 0.00 NaN NaN #> Sigma[5,1] 0.000 0.00 0.00 NaN NaN #> Sigma[5,2] 0.000 0.00 0.00 NaN NaN #> Sigma[5,3] 0.000 0.00 0.00 NaN NaN #> Sigma[5,4] 0.000 0.00 0.00 NaN NaN -#> Sigma[5,5] 0.110 0.21 0.35 1.00 343 +#> Sigma[5,5] 0.098 0.20 0.34 1.03 191 #> #> Approximate significance of GAM process smooths: -#> edf Ref.df Chi.sq p-value -#> te(temp,month) 3.41 15 32.00 0.29 -#> te(temp,month):seriestrend1 1.88 15 1.58 1.00 -#> te(temp,month):seriestrend2 1.35 15 8.75 1.00 -#> te(temp,month):seriestrend3 4.50 15 46.89 0.56 -#> te(temp,month):seriestrend4 3.68 15 6.32 0.87 -#> te(temp,month):seriestrend5 2.43 15 3.59 1.00 +#> edf Ref.df Chi.sq p-value +#> te(temp,month) 4.321 15 35.43 0.11 +#> te(temp,month):seriestrend1 0.929 15 2.18 1.00 +#> te(temp,month):seriestrend2 1.113 15 6.96 1.00 +#> te(temp,month):seriestrend3 4.106 15 48.72 0.41 +#> te(temp,month):seriestrend4 2.831 15 7.54 0.94 +#> te(temp,month):seriestrend5 1.542 15 6.76 1.00 #> #> Stan MCMC diagnostics: #> n_eff / iter looks reasonable for all parameters -#> Rhats above 1.05 found for 14 parameters +#> Rhats above 1.05 found for 38 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 Jul 05 1:09:37 PM 2024. +#> Samples were drawn using NUTS(diag_e) at Tue Sep 03 3:05:05 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) @@ -615,7 +615,7 @@

Inspecting SS modelsA_pars[i, j] <- paste0('A[', i, ',', j, ']') } } -mcmc_plot(var_mod, +mcmc_plot(var_mod, variable = as.vector(t(A_pars)), type = 'hist')

@@ -637,7 +637,7 @@

Inspecting SS modelsSigma_pars[i, j] <- paste0('Sigma[', i, ',', j, ']') } } -mcmc_plot(var_mod, +mcmc_plot(var_mod, variable = as.vector(t(Sigma_pars)), type = 'hist')

@@ -645,7 +645,7 @@

Inspecting SS models
-mcmc_plot(var_mod, variable = 'sigma_obs', regex = TRUE, type = 'hist')
+
mcmc_plot(var_mod, variable = 'sigma_obs', regex = TRUE, type = 'hist')

These are still a bit hard to identify overall, especially when trying to estimate both process and observation error. Often we need to @@ -659,8 +659,8 @@

Correlated process errors
-priors <- c(prior(normal(0.5, 0.1), class = sigma_obs, lb = 0.2),
-            prior(normal(0.5, 0.25), class = sigma))
+priors <- c(prior(normal(0.5, 0.1), class = sigma_obs, lb = 0.2), + prior(normal(0.5, 0.25), class = sigma))

And now we can fit the correlated process error model

 varcor_mod <- mvgam(  
@@ -668,8 +668,8 @@ 

Correlated process errors y ~ -1, # process model formula, which includes the smooth functions - trend_formula = ~ te(temp, month, k = c(4, 4)) + - te(temp, month, k = c(4, 4), by = trend) - 1, + trend_formula = ~ te(temp, month, k = c(4, 4)) + + te(temp, month, k = c(4, 4), by = trend) - 1, # VAR1 model with correlated process errors trend_model = VAR(cor = TRUE), @@ -689,7 +689,7 @@

Correlated process errors Sigma_pars[i, j] <- paste0('Sigma[', i, ',', j, ']') } } -mcmc_plot(varcor_mod, +mcmc_plot(varcor_mod, variable = as.vector(t(Sigma_pars)), type = 'hist')

@@ -706,11 +706,11 @@

Correlated process errors round(median_correlations, 2) #> Bluegreens Diatoms Greens Other.algae Unicells -#> Bluegreens 1.00 -0.04 0.16 -0.04 0.34 -#> Diatoms -0.04 1.00 -0.21 0.48 0.17 -#> Greens 0.16 -0.21 1.00 0.18 0.47 -#> Other.algae -0.04 0.48 0.18 1.00 0.28 -#> Unicells 0.34 0.17 0.47 0.28 1.00 +#> Bluegreens 1.00 -0.04 0.16 -0.06 0.29 +#> Diatoms -0.04 1.00 -0.20 0.50 0.17 +#> Greens 0.16 -0.20 1.00 0.16 0.47 +#> Other.algae -0.06 0.50 0.16 1.00 0.28 +#> Unicells 0.29 0.17 0.47 0.28 1.00

But which model is better? We can compute the variogram score for out of sample forecasts to get a sense of which model does a better job of capturing the dependence structure in the true evaluation set:

@@ -806,9 +806,9 @@

Interested in contributing? + - - + diff --git a/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-12-1.png b/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-12-1.png index 283c3cbf..07de7948 100644 Binary files a/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-12-1.png and b/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-12-1.png differ diff --git a/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-13-1.png b/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-13-1.png index 8a04eaa4..2763603f 100644 Binary files a/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-13-1.png and b/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-13-1.png differ diff --git a/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-14-1.png b/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-14-1.png index 683be347..4eb06df1 100644 Binary files a/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-14-1.png and b/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-14-1.png differ diff --git a/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-15-1.png b/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-15-1.png index 82ff8215..f9119163 100644 Binary files a/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-15-1.png and b/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-15-1.png differ diff --git a/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-16-1.png b/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-16-1.png index 7c5bab02..12ceb628 100644 Binary files a/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-16-1.png and b/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-16-1.png differ diff --git a/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-17-1.png b/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-17-1.png index 8cd7a0b4..863969d1 100644 Binary files a/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-17-1.png and b/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-17-1.png differ diff --git a/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-18-1.png b/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-18-1.png index 22cad379..4a0184ee 100644 Binary files a/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-18-1.png and b/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-18-1.png differ diff --git a/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-19-1.png b/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-19-1.png index e93fee4b..10c823db 100644 Binary files a/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-19-1.png and b/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-19-1.png differ diff --git a/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-20-1.png b/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-20-1.png index debfd45e..d8299458 100644 Binary files a/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-20-1.png and b/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-20-1.png differ diff --git a/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-27-1.png b/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-27-1.png index 915391be..c73b0812 100644 Binary files a/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-27-1.png and b/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-27-1.png differ diff --git a/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-27-2.png b/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-27-2.png index 749ca391..2987feaf 100644 Binary files a/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-27-2.png and b/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-27-2.png differ diff --git a/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-28-1.png b/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-28-1.png index fddaf879..8e5737df 100644 Binary files a/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-28-1.png and b/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-28-1.png differ diff --git a/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-29-1.png b/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-29-1.png index 4abcce5b..560c4dce 100644 Binary files a/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-29-1.png and b/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-29-1.png differ diff --git a/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-30-1.png b/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-30-1.png index ffc05795..dd7e963c 100644 Binary files a/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-30-1.png and b/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-30-1.png differ diff --git a/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-33-1.png b/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-33-1.png index f65d766b..c314f2b2 100644 Binary files a/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-33-1.png and b/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-33-1.png differ diff --git a/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-35-1.png b/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-35-1.png index 2618498a..4980f9a5 100644 Binary files a/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-35-1.png and b/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-35-1.png differ diff --git a/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-36-1.png b/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-36-1.png index beca94c6..3ff7e68f 100644 Binary files a/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-36-1.png and b/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-36-1.png differ diff --git a/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-7-1.png b/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-7-1.png index 6e01007a..7790171c 100644 Binary files a/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-7-1.png and b/docs/articles/trend_formulas_files/figure-html/unnamed-chunk-7-1.png differ diff --git a/docs/reference/Rplot001.png b/docs/reference/Rplot001.png index 79bd0396..17a35806 100644 Binary files a/docs/reference/Rplot001.png and b/docs/reference/Rplot001.png differ diff --git a/docs/reference/Rplot002.png b/docs/reference/Rplot002.png index 16da55c0..d14d14bd 100644 Binary files a/docs/reference/Rplot002.png and b/docs/reference/Rplot002.png differ diff --git a/docs/reference/Rplot003.png b/docs/reference/Rplot003.png index 633fb5a2..81095b14 100644 Binary files a/docs/reference/Rplot003.png and b/docs/reference/Rplot003.png differ diff --git a/docs/reference/Rplot004.png b/docs/reference/Rplot004.png index ba602c61..404a2a85 100644 Binary files a/docs/reference/Rplot004.png and b/docs/reference/Rplot004.png differ diff --git a/docs/reference/figures/README-beta_fc-1.png b/docs/reference/figures/README-beta_fc-1.png index ce3cd414..981099a5 100644 Binary files a/docs/reference/figures/README-beta_fc-1.png and b/docs/reference/figures/README-beta_fc-1.png differ diff --git a/docs/reference/figures/README-unnamed-chunk-12-1.png b/docs/reference/figures/README-unnamed-chunk-12-1.png index b69de8aa..e3a0b6ab 100644 Binary files a/docs/reference/figures/README-unnamed-chunk-12-1.png and b/docs/reference/figures/README-unnamed-chunk-12-1.png differ diff --git a/docs/reference/figures/README-unnamed-chunk-13-1.png b/docs/reference/figures/README-unnamed-chunk-13-1.png index 6b45fc1c..d046d092 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 9b4fb100..69060a98 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 2c1f03cc..046c21bf 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 dca85b6f..bb33b48b 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 ee90a749..2c727472 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 54d461b0..20f53733 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 0ebd4198..0ffe82a9 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 7891d3a3..5dca98be 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 0d8c6e57..0147b2f1 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 db0e6a6d..23160dff 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 55dec179..1f660ff7 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 6b1459fc..b9708f8a 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/index.html b/docs/reference/index.html index 50fbe15e..7cc0b23f 100644 --- a/docs/reference/index.html +++ b/docs/reference/index.html @@ -143,6 +143,11 @@

All functionsirf() + +
Calculate latent VAR impulse response functions
+
+ lfo_cv()
Approximate leave-future-out cross-validation of fitted mvgam objects
@@ -188,7 +193,7 @@

All functionsnuts_params(<mvgam>) log_posterior(<mvgam>) rhat(<mvgam>) neff_ratio(<mvgam>) neff_ratio(<mvgam>) + nuts_params(<mvgam>) log_posterior(<mvgam>) rhat(<mvgam>) neff_ratio(<mvgam>)
Extract diagnostic quantities of mvgam models

@@ -198,7 +203,7 @@

All functionstweedie() student_t() betar() nb() nmix() + tweedie() student_t() betar() nb() lognormal() student() bernoulli() beta_binomial() nmix()

Supported mvgam families
@@ -213,6 +218,11 @@

All functionsmvgam_irf-class +

+
mvgam_irf object description
+
+ get_coef(<mvgam>) set_coef(<mvgam>) get_vcov(<mvgam>) get_predict(<mvgam>) get_data(<mvgam>) get_data(<mvgam_prefit>) find_predictors(<mvgam>) find_predictors(<mvgam_prefit>)
Helper functions for mvgam marginaleffects calculations
@@ -238,6 +248,12 @@

All functionsplot(<mvgam_irf>) + +
Plot impulse responses from an mvgam_irf object +This function takes an mvgam_irf object and produces plots of Impulse Response Functions
+

+ plot(<mvgam_lfo>)
Plot Pareto-k and ELPD values from a leave-future-out object
@@ -359,11 +375,6 @@

All functionsti() te() - -
Defining smooths in mvgam formulae
-

- update(<mvgam>)
Update an existing mvgam object
diff --git a/docs/reference/irf.mvgam-1.png b/docs/reference/irf.mvgam-1.png new file mode 100644 index 00000000..0f799a65 Binary files /dev/null and b/docs/reference/irf.mvgam-1.png differ diff --git a/docs/reference/irf.mvgam-2.png b/docs/reference/irf.mvgam-2.png new file mode 100644 index 00000000..c224e93b Binary files /dev/null and b/docs/reference/irf.mvgam-2.png differ diff --git a/docs/reference/irf.mvgam-3.png b/docs/reference/irf.mvgam-3.png new file mode 100644 index 00000000..6f973ed7 Binary files /dev/null and b/docs/reference/irf.mvgam-3.png differ diff --git a/docs/reference/irf.mvgam-4.png b/docs/reference/irf.mvgam-4.png new file mode 100644 index 00000000..1828da25 Binary files /dev/null and b/docs/reference/irf.mvgam-4.png differ diff --git a/docs/reference/irf.mvgam.html b/docs/reference/irf.mvgam.html new file mode 100644 index 00000000..7be1ea04 --- /dev/null +++ b/docs/reference/irf.mvgam.html @@ -0,0 +1,180 @@ + +Calculate latent VAR impulse response functions — irf.mvgam • mvgam + Skip to contents + + +
+
+
+ +
+

Compute Generalized or Orthogonalized Impulse Response Functions (IRFs) from +mvgam models with Vector Autoregressive dynamics

+
+ +
+

Usage

+
irf(object, ...)
+
+# S3 method for mvgam
+irf(object, h = 1, cumulative = FALSE, orthogonal = FALSE, ...)
+
+ +
+

Arguments

+
object
+

list object of class mvgam resulting from a call to mvgam() +that used a Vector Autoregressive latent process model (either as VAR(cor = FALSE) or +VAR(cor = TRUE))

+ + +
...
+

ignored

+ + +
h
+

Positive integer specifying the forecast horizon over which to calculate +the IRF

+ + +
cumulative
+

Logical flag indicating whether the IRF should be cumulative

+ + +
orthogonal
+

Logical flag indicating whether orthogonalized IRFs should be +calculated. Note that the order of the variables matters when calculating these

+ +
+
+

Value

+ + +

An object of class mvgam_irf containing the posterior IRFs. This +object can be used with the supplied S3 functions plot

+ + +
+
+

Details

+

Generalized or Orthogonalized Impulse Response Functions can be computed +using the posterior estimates of Vector Autoregressive parameters. This function +generates a positive "shock" for a target process at time t = 0 and then +calculates how each of the remaining processes in the latent VAR are expected +to respond over the forecast horizon h. The function computes IRFs for all +processes in the object and returns them in an array that can be plotted using +the S3 plot function

+
+
+

See also

+ +
+
+

Author

+

Nicholas J Clark

+
+ +
+

Examples

+
# \donttest{
+# Simulate some time series that follow a latent VAR(1) process
+simdat <- sim_mvgam(family = gaussian(),
+                    n_series = 4,
+                    trend_model = VAR(cor = TRUE),
+                    prop_trend = 1)
+plot_mvgam_series(data = simdat$data_train, series = 'all')
+
+
+# Fit a model that uses a latent VAR(1)
+mod <- mvgam(y ~ -1,
+             trend_formula = ~ 1,
+             trend_model = VAR(cor = TRUE),
+             family = gaussian(),
+             data = simdat$data_train,
+             silent = 2)
+
+# Calulate Generalized IRFs for each series
+irfs <- irf(mod, h = 12, cumulative = FALSE)
+
+# Plot them
+plot(irfs, series = 1)
+
+plot(irfs, series = 2)
+
+plot(irfs, series = 3)
+
+# }
+
+
+
+ + +
+ + + + + + + diff --git a/docs/reference/mvgam_diagnostics.html b/docs/reference/mvgam_diagnostics.html index 983a5bd4..3b631b5f 100644 --- a/docs/reference/mvgam_diagnostics.html +++ b/docs/reference/mvgam_diagnostics.html @@ -76,9 +76,6 @@

Usage rhat(x, pars = NULL, ...) # S3 method for mvgam -neff_ratio(object, pars = NULL, ...) - -# S3 method for mvgam neff_ratio(object, pars = NULL, ...) @@ -129,42 +126,42 @@

Examples#> Chain 2 Iteration: 1 / 1000 [ 0%] (Warmup) #> Chain 1 Iteration: 100 / 1000 [ 10%] (Warmup) #> Chain 1 Iteration: 200 / 1000 [ 20%] (Warmup) -#> Chain 1 Iteration: 300 / 1000 [ 30%] (Warmup) -#> Chain 1 Iteration: 400 / 1000 [ 40%] (Warmup) #> Chain 2 Iteration: 100 / 1000 [ 10%] (Warmup) #> Chain 2 Iteration: 200 / 1000 [ 20%] (Warmup) +#> Chain 1 Iteration: 300 / 1000 [ 30%] (Warmup) +#> Chain 1 Iteration: 400 / 1000 [ 40%] (Warmup) #> Chain 2 Iteration: 300 / 1000 [ 30%] (Warmup) +#> Chain 2 Iteration: 400 / 1000 [ 40%] (Warmup) #> Chain 1 Iteration: 500 / 1000 [ 50%] (Warmup) #> Chain 1 Iteration: 501 / 1000 [ 50%] (Sampling) #> Chain 1 Iteration: 600 / 1000 [ 60%] (Sampling) -#> Chain 2 Iteration: 400 / 1000 [ 40%] (Warmup) #> Chain 2 Iteration: 500 / 1000 [ 50%] (Warmup) #> Chain 2 Iteration: 501 / 1000 [ 50%] (Sampling) #> Chain 1 Iteration: 700 / 1000 [ 70%] (Sampling) #> Chain 2 Iteration: 600 / 1000 [ 60%] (Sampling) -#> Chain 2 Iteration: 700 / 1000 [ 70%] (Sampling) #> Chain 1 Iteration: 800 / 1000 [ 80%] (Sampling) +#> Chain 2 Iteration: 700 / 1000 [ 70%] (Sampling) +#> Chain 1 Iteration: 900 / 1000 [ 90%] (Sampling) #> Chain 2 Iteration: 800 / 1000 [ 80%] (Sampling) +#> Chain 1 Iteration: 1000 / 1000 [100%] (Sampling) #> Chain 2 Iteration: 900 / 1000 [ 90%] (Sampling) -#> Chain 1 Iteration: 900 / 1000 [ 90%] (Sampling) +#> Chain 1 finished in 0.9 seconds. #> Chain 2 Iteration: 1000 / 1000 [100%] (Sampling) -#> Chain 2 finished in 0.7 seconds. -#> Chain 1 Iteration: 1000 / 1000 [100%] (Sampling) -#> Chain 1 finished in 0.8 seconds. +#> Chain 2 finished in 1.0 seconds. #> #> Both chains finished successfully. -#> Mean chain execution time: 0.7 seconds. -#> Total execution time: 0.9 seconds. +#> Mean chain execution time: 0.9 seconds. +#> Total execution time: 1.1 seconds. #> np <- nuts_params(mod) head(np) #> Chain Iteration Parameter Value -#> 1 1 1 accept_stat__ 0.987513 -#> 2 1 2 accept_stat__ 1.000000 -#> 3 1 3 accept_stat__ 0.993570 -#> 4 1 4 accept_stat__ 0.895941 -#> 5 1 5 accept_stat__ 0.764561 -#> 6 1 6 accept_stat__ 0.925034 +#> 1 1 1 accept_stat__ 0.946851 +#> 2 1 2 accept_stat__ 0.971100 +#> 3 1 3 accept_stat__ 0.980961 +#> 4 1 4 accept_stat__ 0.911601 +#> 5 1 5 accept_stat__ 0.926669 +#> 6 1 6 accept_stat__ 0.837443 # extract the number of divergence transitions sum(subset(np, Parameter == "divergent__")$Value) @@ -172,7 +169,7 @@

Examples head(neff_ratio(mod)) #> mus[1,1] mus[2,1] mus[3,1] mus[4,1] mus[5,1] mus[6,1] -#> 0.7328816 0.8904311 0.7912656 0.6653842 0.6170843 0.5234097 +#> 0.6358219 0.3672060 0.7036268 0.8197350 0.6501380 0.4767678 # } diff --git a/docs/reference/mvgam_families.html b/docs/reference/mvgam_families.html index 6616f662..082a0610 100644 --- a/docs/reference/mvgam_families.html +++ b/docs/reference/mvgam_families.html @@ -71,6 +71,14 @@

Usage nb(...) +lognormal(...) + +student(...) + +bernoulli(...) + +beta_binomial(...) + nmix(link = "log") @@ -230,8 +238,8 @@

Examples data = testdat, # priors can be set in the usual way - priors = c(prior(std_normal(), class = b), - prior(normal(1, 1.5), class = Intercept_trend)), + priors = c(prior(std_normal(), class = b), + prior(normal(1, 1.5), class = Intercept_trend)), chains = 2) #> Your model may benefit from using "noncentred = TRUE" #> Compiling Stan program using cmdstanr @@ -262,24 +270,24 @@

Examples#> Chain 2 Iteration: 900 / 1000 [ 90%] (Sampling) #> Chain 1 Iteration: 900 / 1000 [ 90%] (Sampling) #> Chain 2 Iteration: 1000 / 1000 [100%] (Sampling) -#> Chain 2 finished in 45.2 seconds. +#> Chain 2 finished in 45.9 seconds. #> Chain 1 Iteration: 1000 / 1000 [100%] (Sampling) -#> Chain 1 finished in 45.8 seconds. +#> Chain 1 finished in 46.5 seconds. #> #> Both chains finished successfully. -#> Mean chain execution time: 45.5 seconds. -#> Total execution time: 46.0 seconds. +#> Mean chain execution time: 46.2 seconds. +#> Total execution time: 46.6 seconds. #> # The usual diagnostics summary(mod) #> GAM observation formula: #> obs ~ species - 1 -#> <environment: 0x000001341d5cc658> +#> <environment: 0x000001e60fc56eb8> #> #> GAM process formula: #> ~s(time, by = trend, k = 4) + species -#> <environment: 0x000001341d5cc658> +#> <environment: 0x000001e60fc56eb8> #> #> Family: #> nmix @@ -333,7 +341,7 @@

Examples#> 0 of 1000 iterations saturated the maximum tree depth of 12 (0%) #> E-FMI indicated no pathological behavior #> -#> Samples were drawn using NUTS(diag_e) at Thu Aug 29 11:04:03 AM 2024. +#> Samples were drawn using NUTS(diag_e) at Tue Sep 03 2:04:02 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) @@ -375,6 +383,9 @@

Examplesmod <- mvgam(cbind(y, ntrials) ~ series + s(x, by = series), family = binomial(), data = dat) +#> Warning: Binomial and Beta-binomial families require cbind(n_successes, n_trials) +#> in the formula left-hand side. Do not use cbind(n_successes, n_failures)! +#> This warning is displayed once per session. #> Compiling Stan program using cmdstanr #> #> Start sampling @@ -384,63 +395,63 @@

Examples#> Chain 2 Iteration: 1 / 1000 [ 0%] (Warmup) #> Chain 3 Iteration: 1 / 1000 [ 0%] (Warmup) #> Chain 4 Iteration: 1 / 1000 [ 0%] (Warmup) -#> Chain 4 Iteration: 100 / 1000 [ 10%] (Warmup) -#> Chain 3 Iteration: 100 / 1000 [ 10%] (Warmup) #> Chain 1 Iteration: 100 / 1000 [ 10%] (Warmup) #> Chain 2 Iteration: 100 / 1000 [ 10%] (Warmup) -#> Chain 4 Iteration: 200 / 1000 [ 20%] (Warmup) +#> Chain 3 Iteration: 100 / 1000 [ 10%] (Warmup) #> Chain 1 Iteration: 200 / 1000 [ 20%] (Warmup) -#> Chain 3 Iteration: 200 / 1000 [ 20%] (Warmup) +#> Chain 4 Iteration: 100 / 1000 [ 10%] (Warmup) #> Chain 2 Iteration: 200 / 1000 [ 20%] (Warmup) -#> Chain 3 Iteration: 300 / 1000 [ 30%] (Warmup) -#> Chain 4 Iteration: 300 / 1000 [ 30%] (Warmup) #> Chain 1 Iteration: 300 / 1000 [ 30%] (Warmup) -#> Chain 4 Iteration: 400 / 1000 [ 40%] (Warmup) -#> Chain 3 Iteration: 400 / 1000 [ 40%] (Warmup) -#> Chain 1 Iteration: 400 / 1000 [ 40%] (Warmup) +#> Chain 3 Iteration: 200 / 1000 [ 20%] (Warmup) #> Chain 2 Iteration: 300 / 1000 [ 30%] (Warmup) -#> Chain 4 Iteration: 500 / 1000 [ 50%] (Warmup) -#> Chain 4 Iteration: 501 / 1000 [ 50%] (Sampling) -#> Chain 1 Iteration: 500 / 1000 [ 50%] (Warmup) -#> Chain 3 Iteration: 500 / 1000 [ 50%] (Warmup) -#> Chain 4 Iteration: 600 / 1000 [ 60%] (Sampling) -#> Chain 1 Iteration: 501 / 1000 [ 50%] (Sampling) +#> Chain 1 Iteration: 400 / 1000 [ 40%] (Warmup) #> Chain 2 Iteration: 400 / 1000 [ 40%] (Warmup) -#> Chain 3 Iteration: 501 / 1000 [ 50%] (Sampling) -#> Chain 3 Iteration: 600 / 1000 [ 60%] (Sampling) -#> Chain 1 Iteration: 600 / 1000 [ 60%] (Sampling) +#> Chain 4 Iteration: 200 / 1000 [ 20%] (Warmup) +#> Chain 3 Iteration: 300 / 1000 [ 30%] (Warmup) +#> Chain 1 Iteration: 500 / 1000 [ 50%] (Warmup) #> Chain 2 Iteration: 500 / 1000 [ 50%] (Warmup) -#> Chain 4 Iteration: 700 / 1000 [ 70%] (Sampling) +#> Chain 1 Iteration: 501 / 1000 [ 50%] (Sampling) #> Chain 2 Iteration: 501 / 1000 [ 50%] (Sampling) -#> Chain 3 Iteration: 700 / 1000 [ 70%] (Sampling) +#> Chain 3 Iteration: 400 / 1000 [ 40%] (Warmup) +#> Chain 4 Iteration: 300 / 1000 [ 30%] (Warmup) #> Chain 2 Iteration: 600 / 1000 [ 60%] (Sampling) -#> Chain 4 Iteration: 800 / 1000 [ 80%] (Sampling) -#> Chain 1 Iteration: 700 / 1000 [ 70%] (Sampling) -#> Chain 3 Iteration: 800 / 1000 [ 80%] (Sampling) +#> Chain 3 Iteration: 500 / 1000 [ 50%] (Warmup) +#> Chain 4 Iteration: 400 / 1000 [ 40%] (Warmup) +#> Chain 1 Iteration: 600 / 1000 [ 60%] (Sampling) #> Chain 2 Iteration: 700 / 1000 [ 70%] (Sampling) -#> Chain 4 Iteration: 900 / 1000 [ 90%] (Sampling) +#> Chain 3 Iteration: 501 / 1000 [ 50%] (Sampling) +#> Chain 3 Iteration: 600 / 1000 [ 60%] (Sampling) +#> Chain 4 Iteration: 500 / 1000 [ 50%] (Warmup) +#> Chain 4 Iteration: 501 / 1000 [ 50%] (Sampling) +#> Chain 1 Iteration: 700 / 1000 [ 70%] (Sampling) #> Chain 2 Iteration: 800 / 1000 [ 80%] (Sampling) -#> Chain 3 Iteration: 900 / 1000 [ 90%] (Sampling) +#> Chain 3 Iteration: 700 / 1000 [ 70%] (Sampling) +#> Chain 4 Iteration: 600 / 1000 [ 60%] (Sampling) #> Chain 1 Iteration: 800 / 1000 [ 80%] (Sampling) #> Chain 2 Iteration: 900 / 1000 [ 90%] (Sampling) -#> Chain 3 Iteration: 1000 / 1000 [100%] (Sampling) -#> Chain 4 Iteration: 1000 / 1000 [100%] (Sampling) -#> Chain 3 finished in 3.3 seconds. -#> Chain 4 finished in 3.4 seconds. +#> Chain 4 Iteration: 700 / 1000 [ 70%] (Sampling) +#> Chain 3 Iteration: 800 / 1000 [ 80%] (Sampling) #> Chain 1 Iteration: 900 / 1000 [ 90%] (Sampling) +#> Chain 4 Iteration: 800 / 1000 [ 80%] (Sampling) +#> Chain 3 Iteration: 900 / 1000 [ 90%] (Sampling) #> Chain 2 Iteration: 1000 / 1000 [100%] (Sampling) -#> Chain 2 finished in 3.5 seconds. +#> Chain 2 finished in 3.3 seconds. #> Chain 1 Iteration: 1000 / 1000 [100%] (Sampling) -#> Chain 1 finished in 3.9 seconds. +#> Chain 3 Iteration: 1000 / 1000 [100%] (Sampling) +#> Chain 4 Iteration: 900 / 1000 [ 90%] (Sampling) +#> Chain 1 finished in 3.5 seconds. +#> Chain 3 finished in 3.5 seconds. +#> Chain 4 Iteration: 1000 / 1000 [100%] (Sampling) +#> Chain 4 finished in 3.6 seconds. #> #> All 4 chains finished successfully. #> Mean chain execution time: 3.5 seconds. -#> Total execution time: 4.1 seconds. +#> Total execution time: 3.8 seconds. #> summary(mod) #> GAM formula: #> cbind(y, ntrials) ~ series + s(x, by = series) -#> <environment: 0x000001341d5cc658> +#> <environment: 0x000001e60fc56eb8> #> #> Family: #> binomial @@ -465,42 +476,43 @@

Examples#> #> GAM coefficient (beta) estimates: #> 2.5% 50% 97.5% Rhat n_eff -#> (Intercept) -0.550 -0.41000 -0.280 1.01 1053 -#> seriesseries2 -0.170 0.03100 0.220 1.00 1032 -#> s(x):seriesseries1.1 -0.200 -0.00300 0.180 1.00 947 -#> s(x):seriesseries1.2 -0.250 0.00120 0.230 1.00 855 -#> s(x):seriesseries1.3 -0.056 0.00550 0.100 1.00 665 -#> s(x):seriesseries1.4 -0.130 -0.00390 0.120 1.00 922 -#> s(x):seriesseries1.5 -0.037 0.00180 0.047 1.00 1168 -#> s(x):seriesseries1.6 -0.120 -0.00290 0.110 1.00 823 -#> s(x):seriesseries1.7 -0.021 0.00074 0.028 1.00 708 -#> s(x):seriesseries1.8 -0.500 0.00830 0.560 1.00 794 -#> s(x):seriesseries1.9 0.870 1.20000 1.400 1.00 835 -#> s(x):seriesseries2.1 -0.120 0.04300 0.600 1.02 97 -#> s(x):seriesseries2.2 -0.730 -0.02000 0.220 1.02 123 -#> s(x):seriesseries2.3 -0.200 -0.00400 0.076 1.01 235 -#> s(x):seriesseries2.4 -0.140 -0.00027 0.280 1.01 216 -#> s(x):seriesseries2.5 -0.150 -0.00310 0.041 1.02 147 -#> s(x):seriesseries2.6 -0.120 0.00160 0.300 1.02 131 -#> s(x):seriesseries2.7 -0.047 -0.00047 0.028 1.00 502 -#> s(x):seriesseries2.8 -1.300 0.00750 0.560 1.01 141 -#> s(x):seriesseries2.9 -1.200 -0.95000 -0.160 1.02 106 +#> (Intercept) -0.550 -0.41000 -0.260 1.00 1198 +#> seriesseries2 -0.170 0.02500 0.230 1.00 1140 +#> s(x):seriesseries1.1 -0.220 -0.00640 0.190 1.00 765 +#> s(x):seriesseries1.2 -0.240 0.00210 0.290 1.01 603 +#> s(x):seriesseries1.3 -0.057 0.00640 0.110 1.00 448 +#> s(x):seriesseries1.4 -0.160 -0.00510 0.130 1.00 693 +#> s(x):seriesseries1.5 -0.040 0.00200 0.051 1.00 815 +#> s(x):seriesseries1.6 -0.140 -0.00340 0.120 1.01 620 +#> s(x):seriesseries1.7 -0.022 0.00084 0.028 1.00 951 +#> s(x):seriesseries1.8 -0.550 0.01800 0.650 1.00 654 +#> s(x):seriesseries1.9 0.810 1.20000 1.400 1.00 617 +#> s(x):seriesseries2.1 -0.120 0.05200 0.700 1.05 66 +#> s(x):seriesseries2.2 -0.930 -0.02600 0.260 1.03 122 +#> s(x):seriesseries2.3 -0.260 -0.00420 0.091 1.04 97 +#> s(x):seriesseries2.4 -0.200 0.00014 0.390 1.02 341 +#> s(x):seriesseries2.5 -0.190 -0.00480 0.041 1.03 109 +#> s(x):seriesseries2.6 -0.160 0.00360 0.410 1.03 143 +#> s(x):seriesseries2.7 -0.056 -0.00025 0.039 1.00 691 +#> s(x):seriesseries2.8 -1.700 -0.00180 0.700 1.03 146 +#> s(x):seriesseries2.9 -1.200 -0.94000 -0.038 1.05 65 #> #> Approximate significance of GAM smooths: #> edf Ref.df Chi.sq p-value -#> s(x):seriesseries1 1.15 9 71.6 <2e-16 *** -#> s(x):seriesseries2 1.05 9 54.1 0.13 +#> s(x):seriesseries1 1.06 9 71.5 0.00081 *** +#> s(x):seriesseries2 1.65 9 55.5 < 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 +#> Rhats above 1.05 found for 3 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 Thu Aug 29 11:04:14 AM 2024. +#> Samples were drawn using NUTS(diag_e) at Tue Sep 03 2:04:52 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/reference/mvgam_formulae.html b/docs/reference/mvgam_formulae.html index 37986ada..715d572c 100644 --- a/docs/reference/mvgam_formulae.html +++ b/docs/reference/mvgam_formulae.html @@ -101,6 +101,7 @@

See alsojagam, gam, s, +gp, formula

diff --git a/docs/reference/mvgam_irf-class.html b/docs/reference/mvgam_irf-class.html new file mode 100644 index 00000000..7e8f391e --- /dev/null +++ b/docs/reference/mvgam_irf-class.html @@ -0,0 +1,101 @@ + +mvgam_irf object description — mvgam_irf-class • mvgam + Skip to contents + + +
+
+
+ +
+

A mvgam_irf object returned by function irf. +Run methods(class = "mvgam_irf") to see an overview of available methods.

+
+ + +
+

Details

+

A mvgam_irf object contains a list of posterior IRFs, each stored as +its own list

+
+
+

See also

+ +
+
+

Author

+

Nicholas J Clark

+
+ +
+ + +
+ + + + + + + diff --git a/docs/reference/mvgam_marginaleffects.html b/docs/reference/mvgam_marginaleffects.html index eee73390..bac8d22c 100644 --- a/docs/reference/mvgam_marginaleffects.html +++ b/docs/reference/mvgam_marginaleffects.html @@ -156,8 +156,6 @@

Argumentsdatagrid() call to specify a custom grid of regressors. For example:

  • newdata = datagrid(cyl = c(4, 6)): cyl variable equal to 4 and 6 and other regressors fixed at their means or modes.

  • See the Examples section and the datagrid() documentation.

-
  • subset() call with a single argument to select a subset of the dataset used to fit the model, ex: newdata = subset(treatment == 1)

  • -
  • dplyr::filter() call with a single argument to select a subset of the dataset used to fit the model, ex: newdata = filter(treatment == 1)

  • string: