add vignette 1

.Rbuildignore

@@ -1,3 +1,5 @@
 ^mvgam\.Rproj$
 ^\.Rproj\.user$
 ^LICENSE\.md$
+^doc$
+^Meta$

.gitignore

@@ -3,3 +3,5 @@
 .RData
 .Ruserdata
 NEON_manuscript/Results/
+doc
+Meta

DESCRIPTION

@@ -33,3 +33,5 @@ Encoding: UTF-8
 LazyData: true
 Roxygen: list(markdown = TRUE)
 RoxygenNote: 7.1.1
+Suggests: knitr, rmarkdown
+VignetteBuilder: knitr

NEON_manuscript/amb_hyp_tests.R

@@ -47,9 +47,9 @@ hyp3 = y ~
   s(season, by = series, m = 1, k = 8) - 1
 
 # Fit each hypothesis
-n.burnin = 1000
-n.iter = 1000
-thin = 1
+n.burnin = 100000
+n.iter = 5000
+thin = 5
 
 # Use default priors for latent drift terms and for smooth penalties
 fit_null <- fit_mvgam(data_train = all_data$data_train,

NEON_manuscript/ixodes_hyp_tests.R

@@ -10,7 +10,7 @@ source('NEON_manuscript/neon_utility_functions.R')
 # Prep data
 all_data <- prep_neon_data(species = 'Ixodes_scapularis', split_prop = 0.8)
 
-#### Set hypothtesis formulae ####
+#### Set hypothesis formulae ####
 # NULL. There is no seasonal pattern to be estimated, and we simply let the latent
 # factors and site-level effects of growing days influence the series dynamics
 null_hyp = y ~ siteID + s(cum_gdd, by = siteID, k = 3) - 1
@@ -47,9 +47,9 @@ hyp3 = y ~
   s(season, by = series, m = 1, k = 8) - 1
 
 # Fit each hypothesis
-n.burnin = 2000
-n.iter = 2000
-thin = 1
+n.burnin = 1000
+n.iter = 5000
+thin = 5
 
 fit_null <- fit_mvgam(data_train = all_data$data_train,
                   data_test = all_data$data_test,
@@ -75,13 +75,6 @@ fit_hyp1 <- fit_mvgam(data_train = all_data$data_train,
                       auto_update = F,
                       interval_width = 0.9)
 
-par(mfrow = c(1,1))
-compare_mvgams(model1 = fit_null$out_gam_mod,
-               model2 = fit_hyp1$out_gam_mod,
-               fc_horizon = 6,
-               n_evaluations = 20, n_samples = 1000,
-               n_cores = 2)
-
 fit_hyp2 <- fit_mvgam(data_train = all_data$data_train,
                       data_test = all_data$data_test,
                       formula = hyp2,

NEON_manuscript/simulation_component.R

@@ -6,7 +6,6 @@ source('neon_utility_functions.R')
 # Use hierarchical seasonality for all;
 # trend components: 0.3 and 0.7
 # N_series: 4, 12
-# Overdispersion: 0.5, 100
 # Length: 6 years
 # Missingness: none, 10%, 50%
 set.seed(110)
@@ -22,18 +21,18 @@ sim_results <- lapply(seq_len(nrow(run_parameters)), function(x){
   sim_data <- sim_mvgam(T = run_parameters$T[x],
                         n_series = run_parameters$n_series[x],
                         prop_missing = run_parameters$prop_missing[x],
-                        train_prop = 0.85,
+                        train_prop = 0.833,
                         trend_rel = run_parameters$trend_rel[x],
                         seasonality = 'hierarchical')
 
   # Fit a multivariate gam with no seasonality (null model)
   t <- Sys.time()
   nullmod <- mvjagam(data_train = sim_data$data_train,
                   data_test = sim_data$data_test,
-                  formula = y ~ s(series, bs = 're') - 1,
-                  use_nb = TRUE,
+                  formula = y ~ s(series, bs = 're'),
+                  family = 'nb',
                   use_lv = F,
-                  n.burnin = 60000,
+                  n.burnin = 25000,
                   n.iter = 5000,
                   thin = 5,
                   auto_update = T)
@@ -103,12 +102,13 @@ sim_results <- lapply(seq_len(nrow(run_parameters)), function(x){
   # Fit a multivariate gam with hierarchical seasonality
   t <- Sys.time()
   hier_mod <- mvjagam(data_train = sim_data$data_train,
-                        data_test = sim_data$data_test,
-                        formula = y ~ s(season, k = 4, m = 2, bs = 'cc') +
-                                        s(season, by = series, m = 1, k = 8) - 1,
-                        use_nb = TRUE,
-                        use_lv = T,
-                      n.burnin = 90000,
+                      data_test = sim_data$data_test,
+                      formula = y ~ s(season, k = 4, m = 2, bs = 'cc') +
+                                        s(season, by = series, m = 1, k = 8),
+                      family = 'nb',
+                      trend_model = 'AR3',
+                      use_lv = T,
+                      n.burnin = 25000,
                       n.iter = 5000,
                       thin = 5,
                       auto_update = T)
@@ -132,10 +132,10 @@ sim_results <- lapply(seq_len(nrow(run_parameters)), function(x){
                                                              correlations$mean_correlations))
 
   # Fit a standard mgcv version of the hierarchical model as a final comparison, including
-  # non-wiggly series-specific year smooths
+  # non-wiggly series-specific year smooths with penalty on the first derivative
   mgcv_hier_mod <- gam(y ~ s(season, k = 4, m = 2, bs = 'cc') +
                          s(season, by = series, m = 1, k = 8) +
-                         s(year, by = series, k = 3, bs = 'gp'),
+                         s(year, by = series, k = 3, bs = 'gp', m = 1),
                        family = nb(), data = sim_data$data_train)
 
   # Extract posterior predictions from the mgcv model and calculate out of sample DRPS

R/mvjagam.R

@@ -85,7 +85,7 @@ mvjagam = function(formula,
                     prior_simulation = FALSE,
                     family = 'nb',
                     use_lv = FALSE,
-                    n_lv = 5,
+                    n_lv,
                     trend_model = 'RW',
                     n.chains = 2,
                     n.burnin = 5000,
@@ -111,18 +111,20 @@ mvjagam = function(formula,
   # Estimate the GAM model using mgcv so that the linear predictor matrix can be easily calculated
   # when simulating from the JAGS model later on
   if(!missing(knots)){
-    ss_gam <- mgcv::gam(formula(formula),
+    ss_gam <- mgcv::bam(formula(formula),
                         data = data_train,
-                        method = "REML",
+                        method = "fREML",
                         family = poisson(),
                         drop.unused.levels = FALSE,
-                        knots = knots)
+                        knots = knots,
+                        nthreads = parallel::detectCores()-1)
   } else {
-    ss_gam <- mgcv::gam(formula(formula),
+    ss_gam <- mgcv::bam(formula(formula),
                         data = data_train,
-                        method = "REML",
+                        method = "fREML",
                         family = poisson(),
-                        drop.unused.levels = FALSE)
+                        drop.unused.levels = FALSE,
+                        nthreads = parallel::detectCores()-1)
   }
 
   # Make jags file and appropriate data structures