Changes:

- priors on variances - adding option to monitor user specified pars list in fit_dfa
fate-ewi · Aug 28, 2020 · 2d3dbdd · 2d3dbdd
1 parent efd5241
commit 2d3dbdd
Show file tree

Hide file tree

Showing 36 changed files with 226 additions and 69 deletions.
diff --git a/R/find_dfa_trends.R b/R/find_dfa_trends.R
@@ -22,7 +22,7 @@
 #' s <- sim_dfa(num_trends = 2, num_years = 20, num_ts = 3)
 #' # only 1 chain and 180 iterations used so example runs quickly:
 #' m <- find_dfa_trends(
-#'   y = s$y_sim, iter = 180,
+#'   y = s$y_sim, iter = 50,
 #'   kmin = 1, kmax = 2, chains = 1, compare_normal = FALSE,
 #'   variance = "equal", convergence_threshold = 1.1,
 #'   control = list(adapt_delta = 0.95, max_treedepth = 20))

diff --git a/R/find_regimes.R b/R/find_regimes.R
@@ -15,7 +15,7 @@
 #'
 #' @examples
 #' data(Nile)
-#' find_regimes(log(Nile), iter = 500, chains = 1, max_regimes = 2)
+#' find_regimes(log(Nile), iter = 50, chains = 1, max_regimes = 2)
 
 find_regimes <- function(y,
   sds = NULL,

diff --git a/R/find_swans.R b/R/find_swans.R
@@ -17,7 +17,7 @@
 #' plot(s$y_sim[1,], type = "o")
 #' abline(v = 15, col = "red")
 #' # only 1 chain and 250 iterations used so example runs quickly:
-#' m <- fit_dfa(y = s$y_sim, num_trends = 1, iter = 250, chains = 1, nu_fixed = 2)
+#' m <- fit_dfa(y = s$y_sim, num_trends = 1, iter = 50, chains = 1, nu_fixed = 2)
 #' r <- rotate_trends(m)
 #' p <- plot_trends(r) #+ geom_vline(xintercept = 15, colour = "red")
 #' print(p)

diff --git a/R/fit_dfa.R b/R/fit_dfa.R
@@ -50,6 +50,12 @@
 #' @param z_bound Optional hard constraints for estimated factor loadings -- really only applies to model with 1 trend. Passed in as a 2-element vector representing the lower and upper bound, e.g. (0, 100) to constrain positive
 #' @param z_model Optional argument allowing for elements of Z to be constrained to be proportions (each time series modeled as a mixture of trends). Arguments can be "dfa" (default) or "proportion"
 #' @param ... Any other arguments to pass to [rstan::sampling()].
+#' @param par_list A vector of parameter names of variables to be estimated by Stan. If NULL, this will default to
+#'   c("x", "Z", "sigma", "log_lik", "psi","xstar") for most models -- though if AR / MA, or Student-t models are used
+#'   additional parameters will be monitored. If you want to use diagnostic tools in rstan, including moment_matching,
+#'   you will need to pass in a larger list. Setting this argument to "all" will monitor all parameters, enabling the use
+#'   of diagnostic functions -- but making the models a lot larger for storage. Finally, this argument may be a custom string
+#'   of parameters to monitor, e.g. c("x","sigma")
 #' @details Note that there is nothing restricting the loadings and trends from
 #'   being inverted (i.e. multiplied by `-1`) for a given chain. Therefore, if
 #'   you fit multiple chains, the package will attempt to determine which chains
@@ -67,33 +73,33 @@
 #' set.seed(42)
 #' s <- sim_dfa(num_trends = 1, num_years = 20, num_ts = 3)
 #' # only 1 chain and 250 iterations used so example runs quickly:
-#' m <- fit_dfa(y = s$y_sim, iter = 250, chains = 1)
+#' m <- fit_dfa(y = s$y_sim, iter = 50, chains = 1)
 #'\dontrun{
 #' # example of observation error covariates
 #' obs_covar = expand.grid("time"=1:20,"timeseries"=1:3,"covariate"=1)
 #' obs_covar$value=rnorm(nrow(obs_covar),0,0.1)
-#' m <- fit_dfa(y = s$y_sim, iter = 250, chains = 1, obs_covar=obs_covar)
+#' m <- fit_dfa(y = s$y_sim, iter = 50, chains = 1, obs_covar=obs_covar)
 #'
 #' # example of process error covariates
 #' pro_covar = expand.grid("time"=1:20,"trend"=1:3,"covariate"=1)
 #' pro_covar$value=rnorm(nrow(pro_covar),0,0.1)
-#' m <- fit_dfa(y = s$y_sim, iter = 250, chains = 1, pro_covar=pro_covar)
+#' m <- fit_dfa(y = s$y_sim, iter = 50, chains = 1, pro_covar=pro_covar)
 #'
 #' # example of long format data
 #' s <- sim_dfa(num_trends = 1, num_years = 20, num_ts = 3)
 #' obs <- c(s$y_sim[1,], s$y_sim[2,], s$y_sim[3,])
 #' long = data.frame("obs" = obs, "ts" = sort(rep(1:3,20)), "time" = rep(1:20,3))
-#' m = fit_dfa(y = long, data_shape = "long", iter = 500, chains = 1)
+#' m = fit_dfa(y = long, data_shape = "long", iter = 50, chains = 1)
 #'
 #' # example of model with Z constrained to be proportions and wide format data
 #' s <- sim_dfa(num_trends = 1, num_years = 20, num_ts = 3)
-#' m = fit_dfa(y = s$y_sim, z_model = "proportion", iter = 500, chains = 1)
+#' m = fit_dfa(y = s$y_sim, z_model = "proportion", iter = 50, chains = 1)
 #'
 #' # example of model with Z constrained to be proportions and long format data
 #' s <- sim_dfa(num_trends = 1, num_years = 20, num_ts = 3)
 #' obs <- c(s$y_sim[1,], s$y_sim[2,], s$y_sim[3,])
 #' long = data.frame("obs" = obs, "ts" = sort(rep(1:3,20)), "time" = rep(1:20,3))
-#' m = fit_dfa(y = long, data_shape = "long", z_model = "proportion", iter = 500, chains = 1)
+#' m = fit_dfa(y = long, data_shape = "long", z_model = "proportion", iter = 50, chains = 1)
 #'}
 fit_dfa <- function(y = y,
                     num_trends = 1,
@@ -116,6 +122,7 @@ fit_dfa <- function(y = y,
                     pro_covar = NULL,
                     z_bound = NULL,
                     z_model = c("dfa","proportion"),
+                    par_list = NULL,
                     ...) {
   # check arguments
   data_shape <- match.arg(data_shape, c("wide","long"))
@@ -330,7 +337,13 @@ fit_dfa <- function(y = y,
     n_sigma_process = n_sigma_process
   )
 
-  pars <- c("x", "Z", "sigma", "log_lik", "psi","xstar") # removed pred
+  if(is.null(par_list)) {
+    pars <- c("x", "Z", "sigma", "log_lik", "psi","xstar")
+  } else {
+    if(par_list == "all") pars <- c("x", "Z", "sigma", "log_lik", "psi","xstar",
+      "devs","x0","z","zpos","sigma_process","p_z",
+      "b_obs","b_pro","phi","theta","Lcorr","ymiss","nu") # removed pred
+  }
   if (est_correlation) pars <- c(pars, "Omega", "Sigma") # add correlation matrix
   if (estimate_nu) pars <- c(pars, "nu")
   if (estimate_trend_ar) pars <- c(pars, "phi")

diff --git a/R/fit_regimes.R b/R/fit_regimes.R
@@ -22,7 +22,7 @@
 #'
 #' @examples
 #' data(Nile)
-#' fit_regimes(log(Nile), iter = 1000, n_regimes = 1)
+#' fit_regimes(log(Nile), iter = 50, n_regimes = 1)
 
 fit_regimes <- function(y,
   sds = NULL,

diff --git a/R/invert_chains.R b/R/invert_chains.R
@@ -17,7 +17,7 @@
 #' set.seed(2)
 #' s <- sim_dfa(num_trends = 2)
 #' set.seed(1)
-#' m <- fit_dfa(y = s$y_sim, num_trends = 1, iter = 500, chains = 2)
+#' m <- fit_dfa(y = s$y_sim, num_trends = 1, iter = 30, chains = 2)
 #' # chains were already inverted, but we can redo that, as an example, with:
 #' find_inverted_chains(m$model, plot = TRUE)
 

diff --git a/R/loo.R b/R/loo.R
@@ -19,8 +19,8 @@
 #' @examples
 #' \donttest{
 #' set.seed(1)
-#' s <- sim_dfa(num_trends = 1, num_years = 50, num_ts = 3)
-#' m <- fit_dfa(y = s$y_sim, iter = 300, chains = 1, num_trends = 1)
+#' s <- sim_dfa(num_trends = 1, num_years = 20, num_ts = 3)
+#' m <- fit_dfa(y = s$y_sim, iter = 50, chains = 1, num_trends = 1)
 #' loo(m)
 #' }
 #' @rdname loo

diff --git a/R/plot_fitted.R b/R/plot_fitted.R
@@ -11,7 +11,7 @@
 #' @examples
 #' \donttest{
 #' y <- sim_dfa(num_trends = 2, num_years = 20, num_ts = 4)
-#' m <- fit_dfa(y = y$y_sim, num_trends = 2, iter = 200, chains = 1)
+#' m <- fit_dfa(y = y$y_sim, num_trends = 2, iter = 50, chains = 1)
 #' p <- plot_fitted(m)
 #' print(p)
 #' }

diff --git a/R/plot_loadings.R b/R/plot_loadings.R
@@ -24,7 +24,7 @@
 #' set.seed(42)
 #' s <- sim_dfa(num_trends = 2, num_ts = 4, num_years = 10)
 #' # only 1 chain and 180 iterations used so example runs quickly:
-#' m <- fit_dfa(y = s$y_sim, num_trends = 2, iter = 180, chains = 1)
+#' m <- fit_dfa(y = s$y_sim, num_trends = 2, iter = 50, chains = 1)
 #' r <- rotate_trends(m)
 #' plot_loadings(r, violin = FALSE, facet = TRUE)
 #' plot_loadings(r, violin = FALSE, facet = FALSE)

diff --git a/R/plot_regime_model.R b/R/plot_regime_model.R
@@ -17,7 +17,7 @@
 #' @examples
 #' \donttest{
 #' data(Nile)
-#' m <- fit_regimes(log(Nile), n_regimes = 2, chains = 1, iter = 800)
+#' m <- fit_regimes(log(Nile), n_regimes = 2, chains = 1, iter = 50)
 #' plot_regime_model(m)
 #' plot_regime_model(m, plot_prob_indices=c(2))
 #' plot_regime_model(m, type = "means")

diff --git a/R/plot_trends.R b/R/plot_trends.R
@@ -16,7 +16,7 @@
 #' @examples
 #' set.seed(1)
 #' s <- sim_dfa(num_trends = 1)
-#' m <- fit_dfa(y = s$y_sim, num_trends = 1, iter = 500, chains = 1)
+#' m <- fit_dfa(y = s$y_sim, num_trends = 1, iter = 50, chains = 1)
 #' r <- rotate_trends(m)
 #' p <- plot_trends(r)
 #' print(p)

diff --git a/R/predicted.R b/R/predicted.R
@@ -10,7 +10,7 @@
 #' set.seed(42)
 #' s <- sim_dfa(num_trends = 1, num_years = 20, num_ts = 3)
 #' # only 1 chain and 1000 iterations used so example runs quickly:
-#' m <- fit_dfa(y = s$y_sim, iter = 1000, chains = 1)
+#' m <- fit_dfa(y = s$y_sim, iter = 50, chains = 1)
 #' pred <- predicted(m)
 #'
 predicted <- function(fitted_model) {

diff --git a/R/rotate_trends.R b/R/rotate_trends.R
@@ -12,7 +12,7 @@
 #' set.seed(42)
 #' s <- sim_dfa(num_trends = 1, num_years = 20, num_ts = 3)
 #' # only 1 chain and 800 iterations used so example runs quickly:
-#' m <- fit_dfa(y = s$y_sim, iter = 800, chains = 1)
+#' m <- fit_dfa(y = s$y_sim, iter = 50, chains = 1)
 #' r <- rotate_trends(m)
 #' plot_trends(r)
 

diff --git a/R/trend_cor.R b/R/trend_cor.R
@@ -33,7 +33,7 @@
 #' @examples
 #' set.seed(1)
 #' s <- sim_dfa(num_trends = 1, num_years = 15)
-#' m <- fit_dfa(y = s$y_sim, num_trends = 1, iter = 500, chains = 1)
+#' m <- fit_dfa(y = s$y_sim, num_trends = 1, iter = 50, chains = 1)
 #' r <- rotate_trends(m)
 #' n_years <- ncol(r$trends[,1,])
 #' fake_dat <- rnorm(n_years, 0, 1)

diff --git a/man/find_dfa_trends.Rd b/man/find_dfa_trends.Rd
diff --git a/man/find_inverted_chains.Rd b/man/find_inverted_chains.Rd
diff --git a/man/find_regimes.Rd b/man/find_regimes.Rd
diff --git a/man/find_swans.Rd b/man/find_swans.Rd
diff --git a/man/fit_dfa.Rd b/man/fit_dfa.Rd
diff --git a/man/fit_regimes.Rd b/man/fit_regimes.Rd
diff --git a/man/loo.Rd b/man/loo.Rd
diff --git a/man/plot_fitted.Rd b/man/plot_fitted.Rd
diff --git a/man/plot_loadings.Rd b/man/plot_loadings.Rd
diff --git a/man/plot_regime_model.Rd b/man/plot_regime_model.Rd
diff --git a/man/plot_trends.Rd b/man/plot_trends.Rd
diff --git a/man/predicted.Rd b/man/predicted.Rd
diff --git a/man/rotate_trends.Rd b/man/rotate_trends.Rd
diff --git a/man/trend_cor.Rd b/man/trend_cor.Rd
diff --git a/src/stan_files/dfa.hpp b/src/stan_files/dfa.hpp
@@ -2076,7 +2076,7 @@ class model_dfa
             if (as_bool(est_sigma_process)) {
 
                 current_statement_begin__ = 388;
-                lp_accum__.add(student_t_log<propto__>(sigma_process, 3, 0, 2));
+                lp_accum__.add(student_t_log<propto__>(sigma_process, 3, 0, 5));
             }
             current_statement_begin__ = 391;
             if (as_bool(logical_eq(proportional_model, 0))) {
@@ -2095,7 +2095,7 @@ class model_dfa
                 }
             }
             current_statement_begin__ = 402;
-            lp_accum__.add(student_t_log<propto__>(sigma, 3, 0, 2));
+            lp_accum__.add(student_t_log<propto__>(sigma, 3, 0, 5));
             current_statement_begin__ = 403;
             if (as_bool(logical_eq(est_cor, 1))) {