diff --git a/R/jsdgam.R b/R/jsdgam.R index c46f9a22..5b21e693 100644 --- a/R/jsdgam.R +++ b/R/jsdgam.R @@ -329,8 +329,8 @@ jsdgam = function(formula, parallel = TRUE, threads = 1, silent = 1, - max_treedepth = 12, - adapt_delta = 0.85, + max_treedepth = 10, + adapt_delta = 0.8, backend = getOption("brms.backend", "cmdstanr"), algorithm = getOption("brms.algorithm", "sampling"), run_model = TRUE, @@ -384,7 +384,7 @@ jsdgam = function(formula, model_file[grep('int n_lv; // number of dynamic factors', model_file, fixed = TRUE)] <- paste0( 'int n_lv; // number of dynamic factors\n', - 'int M; // number of nonzero factor loadings' + 'int M; // number of nonzero lower-triangular factor loadings' ) model_file <- readLines(textConnection(model_file), n = -1) @@ -396,8 +396,10 @@ jsdgam = function(formula, model_file[grep("matrix[n, n_lv] LV_raw;", model_file, fixed = TRUE)] <- paste0( "matrix[n, n_lv] LV_raw;\n\n", - "// factor lower triangle loading coefficients\n", - "vector[M] L;") + "// factor lower triangle loadings\n", + "vector[M] L_lower;\n", + "// factor diagonal loadings\n", + "vector[n_lv] L_diag;") model_file <- readLines(textConnection(model_file), n = -1) # Update transformed parameters @@ -424,17 +426,19 @@ jsdgam = function(formula, "trend_mus = X_trend * b_trend;\n\n", "// constraints allow identifiability of loadings\n", "{\n", - "int index;\n", - "index = 0;\n", - "for (j in 1 : n_lv) {\n", - "for (i in j : n_series) {\n", - "index = index + 1;\n", - "lv_coefs[i, j] = L[index];\n", + "int idx;\n", + "idx = 0;\n", + "for(j in 1 : n_lv) lv_coefs[j, j] = L_diag[j];\n", + "for(j in 1 : n_lv) {\n", + "for(k in (j + 1) : n_series) {\n", + "idx = idx + 1;\n", + "lv_coefs[k, j] = L_lower[idx];\n", "}\n", "}\n", "}\n\n", "// raw latent factors (with linear predictors)\n", - "for (j in 1 : n_lv) {\n", "for (i in 1 : n) {\n", + "for (j in 1 : n_lv) {\n", + "for (i in 1 : n) {\n", "LV[i, j] = trend_mus[ytimes_trend[i, j]] + LV_raw[i, j];\n", "}\n}\n") @@ -448,8 +452,9 @@ jsdgam = function(formula, model_file <- model_file[-sigma_prior] model_file[grep("// priors for latent state SD parameters", model_file, fixed = TRUE)] <- paste0( - "// priors for factor loading coefficients\n", - "L ~ std_normal();") + "// priors for factors and loading coefficients\n", + "L_lower ~ student_t(3, 0, 1);\n", + "L_diag ~ student_t(3, 0, 1);") model_file <- readLines(textConnection(model_file), n = -1) # Update generated quantities @@ -471,7 +476,7 @@ jsdgam = function(formula, # Add M to model_data n_series <- NCOL(model_data$ytimes) - model_data$M <- n_lv * (n_series - n_lv) + n_lv * (n_lv - 1) / 2 + n_lv + model_data$M <- n_lv * (n_series - n_lv) + n_lv * (n_lv - 1) / 2 #### Autoformat the Stan code #### if(requireNamespace('cmdstanr', quietly = TRUE) & backend == 'cmdstanr'){ diff --git a/R/mvgam.R b/R/mvgam.R index 1305d92d..58766e0f 100644 --- a/R/mvgam.R +++ b/R/mvgam.R @@ -208,7 +208,7 @@ #' variables defined in Stan's \code{parameters} block should be saved #' (default is \code{FALSE}). #'@param max_treedepth positive integer placing a cap on the number of simulation steps evaluated during each iteration when -#'`use_stan == TRUE`. Default is `12`. Increasing this value can sometimes help with exploration of complex +#'`use_stan == TRUE`. Default is `10`. Increasing this value can sometimes help with exploration of complex #'posterior geometries, but it is rarely fruitful to go above a `max_treedepth` of `14` #'@param adapt_delta positive numeric between `0` and `1` defining the target average proposal acceptance probability #'during Stan's adaptation period, if `use_stan == TRUE`. Default is `0.8`. In general you should not need to change adapt_delta @@ -654,8 +654,8 @@ mvgam = function(formula, algorithm = getOption("brms.algorithm", "sampling"), autoformat = TRUE, save_all_pars = FALSE, - max_treedepth = 12, - adapt_delta = 0.85, + max_treedepth = 10, + adapt_delta = 0.8, silent = 1, jags_path, ...){ diff --git a/R/stan_utils.R b/R/stan_utils.R index 06a3ccab..84970746 100644 --- a/R/stan_utils.R +++ b/R/stan_utils.R @@ -3372,6 +3372,14 @@ check_n_eff <- function(fit, quiet=FALSE, fit_summary, ignore_b_trend = FALSE) { fit_summary <- fit_summary[-grep('L[', rownames(fit_summary), fixed = TRUE), ] } + if(any(grepl('L_diag[', rownames(fit_summary), fixed = TRUE))){ + fit_summary <- fit_summary[-grep('L_diag[', + rownames(fit_summary), fixed = TRUE), ] + } + if(any(grepl('L_lower[', rownames(fit_summary), fixed = TRUE))){ + fit_summary <- fit_summary[-grep('L_lower[', + rownames(fit_summary), fixed = TRUE), ] + } if(any(grepl('LV_raw[', rownames(fit_summary), fixed = TRUE))){ fit_summary <- fit_summary[-grep('LV_raw[', rownames(fit_summary), fixed = TRUE), ] @@ -3432,6 +3440,14 @@ check_rhat <- function(fit, quiet = FALSE, fit_summary, ignore_b_trend = FALSE) fit_summary <- fit_summary[-grep('L[', rownames(fit_summary), fixed = TRUE), ] } + if(any(grepl('L_diag[', rownames(fit_summary), fixed = TRUE))){ + fit_summary <- fit_summary[-grep('L_diag[', + rownames(fit_summary), fixed = TRUE), ] + } + if(any(grepl('L_lower[', rownames(fit_summary), fixed = TRUE))){ + fit_summary <- fit_summary[-grep('L_lower[', + rownames(fit_summary), fixed = TRUE), ] + } if(any(grepl('LV_raw[', rownames(fit_summary), fixed = TRUE))){ fit_summary <- fit_summary[-grep('LV_raw[', rownames(fit_summary), fixed = TRUE), ] diff --git a/man/jsdgam.Rd b/man/jsdgam.Rd index aab3fc41..a5344d30 100644 --- a/man/jsdgam.Rd +++ b/man/jsdgam.Rd @@ -24,8 +24,8 @@ jsdgam( parallel = TRUE, threads = 1, silent = 1, - max_treedepth = 12, - adapt_delta = 0.85, + max_treedepth = 10, + adapt_delta = 0.8, backend = getOption("brms.backend", "cmdstanr"), algorithm = getOption("brms.algorithm", "sampling"), run_model = TRUE, @@ -139,7 +139,7 @@ actual sampling progress is still printed. Set \code{refresh = 0} to turn this o progress bars.} \item{max_treedepth}{positive integer placing a cap on the number of simulation steps evaluated during each iteration when -\code{use_stan == TRUE}. Default is \code{12}. Increasing this value can sometimes help with exploration of complex +\code{use_stan == TRUE}. Default is \code{10}. Increasing this value can sometimes help with exploration of complex posterior geometries, but it is rarely fruitful to go above a \code{max_treedepth} of \code{14}} \item{adapt_delta}{positive numeric between \code{0} and \code{1} defining the target average proposal acceptance probability diff --git a/man/mvgam.Rd b/man/mvgam.Rd index 9f1515eb..1f574111 100644 --- a/man/mvgam.Rd +++ b/man/mvgam.Rd @@ -39,8 +39,8 @@ mvgam( algorithm = getOption("brms.algorithm", "sampling"), autoformat = TRUE, save_all_pars = FALSE, - max_treedepth = 12, - adapt_delta = 0.85, + max_treedepth = 10, + adapt_delta = 0.8, silent = 1, jags_path, ... @@ -278,7 +278,7 @@ variables defined in Stan's \code{parameters} block should be saved (default is \code{FALSE}).} \item{max_treedepth}{positive integer placing a cap on the number of simulation steps evaluated during each iteration when -\code{use_stan == TRUE}. Default is \code{12}. Increasing this value can sometimes help with exploration of complex +\code{use_stan == TRUE}. Default is \code{10}. Increasing this value can sometimes help with exploration of complex posterior geometries, but it is rarely fruitful to go above a \code{max_treedepth} of \code{14}} \item{adapt_delta}{positive numeric between \code{0} and \code{1} defining the target average proposal acceptance probability