From 8b11e80a2de671db29ee23f3c0a3689c86be4322 Mon Sep 17 00:00:00 2001 From: ecmerkle Date: Wed, 11 Sep 2024 09:57:33 -0500 Subject: [PATCH 1/8] allow shrink_t through prior check --- R/stanmarg_data.R | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/R/stanmarg_data.R b/R/stanmarg_data.R index 5a729162..2d7344b1 100644 --- a/R/stanmarg_data.R +++ b/R/stanmarg_data.R @@ -184,12 +184,18 @@ format_priors <- function(lavpartable, level = 1L) { # @return nothing check_priors <- function(lavpartable) { right_pris <- sapply(dpriors(target = "stan"), function(x) strsplit(x, "[, ()]+")[[1]][1]) + ## add additional prior options here + new_pri <- rep("shrink_t", 4); names(new_pri) <- c("nu", "alpha", "lambda", "beta") + right_pris <- c(right_pris, new_pri) pt_pris <- sapply(lavpartable$prior[lavpartable$prior != ""], function(x) strsplit(x, "[, ()]+")[[1]][1]) names(pt_pris) <- lavpartable$mat[lavpartable$prior != ""] right_pris <- c(right_pris, lvrho = "lkj_corr") primatch <- match(names(pt_pris), names(right_pris)) - badpris <- which(pt_pris != right_pris[primatch]) - + badpris <- rep(FALSE, length(pt_pris)) + for (i in 1:length(pt_pris)) { + badpris[i] <- !(pt_pris[i] %in% right_pris[names(right_pris) == names(pt_pris)[i]]) + } + badpris <- which(badpris) ## lvrho entries could also receive beta priors okpris <- which(names(pt_pris[badpris]) == "lvrho" & pt_pris[badpris] == "beta") if (length(okpris) > 0) badpris <- badpris[-okpris] From 23c0a30671042968fb99b3441010863c2da662bd Mon Sep 17 00:00:00 2001 From: ecmerkle Date: Wed, 11 Sep 2024 10:47:12 -0500 Subject: [PATCH 2/8] remove unused prior matrices --- R/stanmarg_data.R | 10 ++-------- 1 file changed, 2 insertions(+), 8 deletions(-) diff --git a/R/stanmarg_data.R b/R/stanmarg_data.R index 2d7344b1..c0dba52d 100644 --- a/R/stanmarg_data.R +++ b/R/stanmarg_data.R @@ -74,23 +74,17 @@ format_priors <- function(lavpartable, level = 1L) { } transtab <- list(c('lambda_y_mn', 'lambda_y_sd', 'len_lam_y'), - c('lambda_x_mn', 'lambda_x_sd', 'len_lam_x'), - c('gamma_mn', 'gamma_sd', 'len_gam'), c('b_mn', 'b_sd', 'len_b'), c('theta_sd_shape', 'theta_sd_rate', 'len_thet_sd', 'theta_pow'), - c('theta_x_sd_shape', 'theta_x_sd_rate', 'len_thet_x_sd', 'theta_x_pow'), c('theta_r_alpha', 'theta_r_beta', 'len_thet_r'), - c('theta_x_r_alpha', 'theta_x_r_beta', 'len_thet_x_r'), c('psi_sd_shape', 'psi_sd_rate', 'len_psi_sd', 'psi_pow'), c('psi_r_alpha', 'psi_r_beta', 'len_psi_r'), - c('phi_sd_shape', 'phi_sd_rate', 'len_phi_sd', 'phi_pow'), - c('phi_r_alpha', 'phi_r_beta', 'len_phi_r'), c('nu_mn', 'nu_sd', 'len_nu'), c('alpha_mn', 'alpha_sd', 'len_alph'), c('tau_mn', 'tau_sd', 'len_tau')) - mats <- c('lambda', 'lambda_x', 'gamma', 'beta', 'thetavar', 'cov.xvar', 'thetaoff', - 'cov.xoff', 'psivar', 'psioff', 'phivar', 'phioff', 'nu', 'alpha', 'tau') + mats <- c('lambda', 'beta', 'thetavar', 'thetaoff', + 'psivar', 'psioff', 'nu', 'alpha', 'tau') if (level == 2L) { newmats <- c('lambda', 'beta', 'thetavar', 'thetaoff', 'psivar', 'psioff', 'nu', 'alpha') subloc <- match(newmats, mats) From ee290c3614d5770ef99ef2908483271e7a7fc0cf Mon Sep 17 00:00:00 2001 From: ecmerkle Date: Wed, 11 Sep 2024 11:23:41 -0500 Subject: [PATCH 3/8] format data on shrinkage priors and blocks --- R/stanmarg_data.R | 33 ++++++++++++++++++++++++++++----- 1 file changed, 28 insertions(+), 5 deletions(-) diff --git a/R/stanmarg_data.R b/R/stanmarg_data.R index c0dba52d..f10bf56a 100644 --- a/R/stanmarg_data.R +++ b/R/stanmarg_data.R @@ -73,14 +73,14 @@ format_priors <- function(lavpartable, level = 1L) { lavpartable <- lavpartable[order(lavpartable$col, lavpartable$row),] } - transtab <- list(c('lambda_y_mn', 'lambda_y_sd', 'len_lam_y'), - c('b_mn', 'b_sd', 'len_b'), + transtab <- list(c('lambda_y_mn', 'lambda_y_sd', 'len_lam_y', 'lambda_y_pri', 'lambda_y_blk'), + c('b_mn', 'b_sd', 'len_b', 'b_pri', 'b_blk'), c('theta_sd_shape', 'theta_sd_rate', 'len_thet_sd', 'theta_pow'), c('theta_r_alpha', 'theta_r_beta', 'len_thet_r'), c('psi_sd_shape', 'psi_sd_rate', 'len_psi_sd', 'psi_pow'), c('psi_r_alpha', 'psi_r_beta', 'len_psi_r'), - c('nu_mn', 'nu_sd', 'len_nu'), - c('alpha_mn', 'alpha_sd', 'len_alph'), + c('nu_mn', 'nu_sd', 'len_nu', 'nu_pri', 'nu_blk'), + c('alpha_mn', 'alpha_sd', 'len_alph', 'alpha_pri', 'alpha_blk'), c('tau_mn', 'tau_sd', 'len_tau')) mats <- c('lambda', 'beta', 'thetavar', 'thetaoff', @@ -94,6 +94,14 @@ format_priors <- function(lavpartable, level = 1L) { } out <- list() + + ## if we have prior blocks specified via <.>, number them for the whole partable + blkpris <- grep("", lavpartable$prior) + blknum <- rep(0, length(lavpartable$prior)) + if (length(blkpris) > 0) { + blknum[blkpris] <- as.numeric( as.factor(lavpartable$prior[blkpris]) ) + } + lavpartable$blknum <- blknum for (i in 1:length(mats)) { mat <- origmat <- mats[i] @@ -119,6 +127,8 @@ format_priors <- function(lavpartable, level = 1L) { prisel <- prisel & (lavpartable$free > 0) thepris <- lavpartable$prior[prisel] + priblks <- lavpartable$blknum[prisel] + blkmats <- mat %in% c("nu", "lambda", "beta", "alpha") if (length(thepris) > 0) { textpris <- thepris[thepris != ""] @@ -126,14 +136,21 @@ format_priors <- function(lavpartable, level = 1L) { prisplit <- strsplit(textpris, "[, ()]+") param1 <- sapply(prisplit, function(x) x[2]) + prinms <- sapply(prisplit, function(x) x[1]) - if (!grepl("\\[", prisplit[[1]][3])) { + if (!grepl("\\[", prisplit[[1]][3]) & !blkmats) { param2 <- sapply(prisplit, function(x) x[3]) if (any(is.na(param2)) & mat == "lvrho") { ## omit lkj here param1 <- param1[!is.na(param2)] param2 <- param2[!is.na(param2)] } + } else if (blkmats) { + param2 <- sapply(prisplit, function(x) x[3]) + param2 <- as.numeric(param2) + param2[is.na(param2)] <- 1 ## we don't need it, but we don't want NA to go to Stan + pritype <- array(0, length(param2)) + pritype[prinms == "shrink_t"] <- 1 } else { param2 <- rep(NA, length(param1)) } @@ -157,6 +174,8 @@ format_priors <- function(lavpartable, level = 1L) { param1 <- array(0, 0) param2 <- array(0, 0) powpar <- 1 + pritype <- array(0, 0) + priblks <- array(0, 0) } out[[ transtab[[i]][1] ]] <- param1 @@ -166,6 +185,10 @@ format_priors <- function(lavpartable, level = 1L) { if (origmat %in% c('thetavar', 'cov.xvar', 'psivar', 'phivar')) { out[[ transtab[[i]][4] ]] <- powpar } + if (blkmats) { + out[[ transtab[[i]][4] ]] <- pritype + out[[ transtab[[i]][5] ]] <- priblks + } } # mats return(out) From a4db08348106898dcb7fedb22659c17db14e78b4 Mon Sep 17 00:00:00 2001 From: ecmerkle Date: Wed, 11 Sep 2024 11:35:12 -0500 Subject: [PATCH 4/8] avoid coercion warning --- R/stanmarg_data.R | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/R/stanmarg_data.R b/R/stanmarg_data.R index f10bf56a..ddd64587 100644 --- a/R/stanmarg_data.R +++ b/R/stanmarg_data.R @@ -146,11 +146,11 @@ format_priors <- function(lavpartable, level = 1L) { param2 <- param2[!is.na(param2)] } } else if (blkmats) { + pritype <- array(0, length(param1)) + pritype[prinms == "shrink_t"] <- 1 param2 <- sapply(prisplit, function(x) x[3]) + param2[pritype == 1] <- 1 ## we don't need it, but we don't want NA to go to Stan param2 <- as.numeric(param2) - param2[is.na(param2)] <- 1 ## we don't need it, but we don't want NA to go to Stan - pritype <- array(0, length(param2)) - pritype[prinms == "shrink_t"] <- 1 } else { param2 <- rep(NA, length(param1)) } From a392e8cdcdf606569fb55ded616c9fd58260ac3b Mon Sep 17 00:00:00 2001 From: ecmerkle Date: Wed, 11 Sep 2024 11:40:01 -0500 Subject: [PATCH 5/8] collect info about prior blocks --- R/stanmarg_data.R | 3 +++ 1 file changed, 3 insertions(+) diff --git a/R/stanmarg_data.R b/R/stanmarg_data.R index ddd64587..27ca42a6 100644 --- a/R/stanmarg_data.R +++ b/R/stanmarg_data.R @@ -458,6 +458,9 @@ stanmarg_data <- function(YX = NULL, S = NULL, YXo = NULL, N, Ng, grpnum, # data dat <- c(dat, format_priors(lavpartable[lavpartable$level == levlabs[1],])) dat <- c(dat, format_priors(lavpartable[lavpartable$level == levlabs[2],], level = 2L)) } + priblks <- table(with(dat, c(lambda_y_blk, b_blk, nu_blk, alpha_blk)))[-1] + dat$npriblks <- length(priblks) + dat$priblklen <- max(priblks) return(dat) } From 0a91c0b8c792ead46e06e63be48154bc0a5811f9 Mon Sep 17 00:00:00 2001 From: ecmerkle Date: Sat, 14 Sep 2024 15:09:15 -0500 Subject: [PATCH 6/8] improve computation of prior blocks --- R/stanmarg_data.R | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/R/stanmarg_data.R b/R/stanmarg_data.R index 27ca42a6..4a9796f5 100644 --- a/R/stanmarg_data.R +++ b/R/stanmarg_data.R @@ -460,7 +460,8 @@ stanmarg_data <- function(YX = NULL, S = NULL, YXo = NULL, N, Ng, grpnum, # data } priblks <- table(with(dat, c(lambda_y_blk, b_blk, nu_blk, alpha_blk)))[-1] dat$npriblks <- length(priblks) - dat$priblklen <- max(priblks) + dat$priblklen <- 0 + if (dat$npriblks > 0) dat$priblklen <- max(priblks) return(dat) } From 83d2a1c31cbfe0e5f506a6173324e81a7c74dcdf Mon Sep 17 00:00:00 2001 From: ecmerkle Date: Wed, 18 Sep 2024 12:26:39 -0500 Subject: [PATCH 7/8] additions to stanmarg for shrinkage priors --- R/stanmarg_data.R | 20 ++++++++-- inst/stan/stanmarg.stan | 86 +++++++++++++++++++++++++++++++++++++---- 2 files changed, 94 insertions(+), 12 deletions(-) diff --git a/R/stanmarg_data.R b/R/stanmarg_data.R index 4a9796f5..8045195c 100644 --- a/R/stanmarg_data.R +++ b/R/stanmarg_data.R @@ -95,6 +95,11 @@ format_priors <- function(lavpartable, level = 1L) { out <- list() + ## shrinkage priors without <.> + shrpris <- which(grepl("shrink_t", lavpartable$prior) & !grepl("", lavpartable$prior)) + if (length(shrpris) > 0) { + lavpartable$prior[shrpris] <- paste0(lavpartable$prior[shrpris], "<999>") + } ## if we have prior blocks specified via <.>, number them for the whole partable blkpris <- grep("", lavpartable$prior) blknum <- rep(0, length(lavpartable$prior)) @@ -149,7 +154,6 @@ format_priors <- function(lavpartable, level = 1L) { pritype <- array(0, length(param1)) pritype[prinms == "shrink_t"] <- 1 param2 <- sapply(prisplit, function(x) x[3]) - param2[pritype == 1] <- 1 ## we don't need it, but we don't want NA to go to Stan param2 <- as.numeric(param2) } else { param2 <- rep(NA, length(param1)) @@ -458,11 +462,19 @@ stanmarg_data <- function(YX = NULL, S = NULL, YXo = NULL, N, Ng, grpnum, # data dat <- c(dat, format_priors(lavpartable[lavpartable$level == levlabs[1],])) dat <- c(dat, format_priors(lavpartable[lavpartable$level == levlabs[2],], level = 2L)) } - priblks <- table(with(dat, c(lambda_y_blk, b_blk, nu_blk, alpha_blk)))[-1] + allblks <- with(dat, c(lambda_y_blk, b_blk, nu_blk, alpha_blk)) + priblks <- table(c(0, allblks))[-1] dat$npriblks <- length(priblks) dat$priblklen <- 0 - if (dat$npriblks > 0) dat$priblklen <- max(priblks) - + dat$blkparm1 <- array(0, 0) + dat$blkparm2 <- array(0, 0) + if (dat$npriblks > 0) { + dat$priblklen <- max(priblks) + allparm1 <- with(dat, c(lambda_y_mn, b_mn, nu_mn, alpha_mn)) + dat$blkparm1 <- array(tapply(allparm1[allblks > 0], allblks[allblks > 0], head, 1), dat$npriblks) + allparm2 <- with(dat, c(lambda_y_sd, b_sd, nu_sd, alpha_sd)) + dat$blkparm2 <- array(tapply(allparm2[allblks > 0], allblks[allblks > 0], head, 1), dat$npriblks) + } return(dat) } diff --git a/inst/stan/stanmarg.stan b/inst/stan/stanmarg.stan index 639d1b38..925ae394 100644 --- a/inst/stan/stanmarg.stan +++ b/inst/stan/stanmarg.stan @@ -557,7 +557,68 @@ functions { // you can use these in R following `rstan::expose_stan_functions("f } return out; - } + } + + real eval_priors(vector Lambda_y_free, vector B_free, vector Nu_free, vector Alpha_free, vector sd0, vector lambda_y_primn, array[] real lambda_y_sd, array[] int lambda_y_pri, array[] int lambda_y_blk, vector b_primn, array[] real b_sd, array[] int b_pri, array[] int b_blk, vector nu_primn, array[] real nu_sd, array[] int nu_pri, array[] int nu_blk, vector alpha_primn, array[] real alpha_sd, array[] int alpha_pri, array[] int alpha_blk, array[] int len_free, vector blkparm1, vector blkparm2, int priblklen, int npriblks) { + real out = 0.0; + if (npriblks == 0) { + out += normal_lpdf(Lambda_y_free | lambda_y_primn, lambda_y_sd); + out += normal_lpdf(B_free | b_primn, b_sd); + out += normal_lpdf(Nu_free | nu_primn, nu_sd); + out += normal_lpdf(Alpha_free | alpha_primn, alpha_sd); + } else { + for (i in 1:npriblks) { + vector[priblklen] parvec; + int npars = 1; + if (len_free[1] > 0) { + for (j in 1:len_free[1]) { + if (lambda_y_pri[j] == 0) { + out += normal_lpdf(Lambda_y_free[j] | lambda_y_primn[j], lambda_y_sd[j]); + } else if (lambda_y_blk[j] == i && lambda_y_pri[j] == 1) { + parvec[npars] = Lambda_y_free[j]; + npars += 1; + } + } + } + + if (len_free[4] > 0) { + for (j in 1:len_free[4]) { + if (b_pri[j] == 0) { + out += normal_lpdf(B_free[j] | b_primn[j], b_sd[j]); + } else if (b_blk[j] == i && b_pri[j] == 1) { + parvec[npars] = B_free[j]; + npars += 1; + } + } + } + + if (len_free[13] > 0) { + for (j in 1:len_free[13]) { + if (nu_pri[j] == 0) { + out += normal_lpdf(Nu_free[j] | nu_primn[j], nu_sd[j]); + } else if (nu_blk[j] == i && nu_pri[j] == 1) { + parvec[npars] = Nu_free[j]; + npars += 1; + } + } + } + + if (len_free[14] > 0) { + for (j in 1:len_free[14]) { + if (alpha_pri[j] == 0) { + out += normal_lpdf(Alpha_free[j] | alpha_primn[j], alpha_sd[j]); + } else if (alpha_blk[j] == i && alpha_pri[j] == 1) { + parvec[npars] = Alpha_free[j]; + npars += 1; + } + } + } + out += normal_lpdf(parvec[1:(npars - 1)] | 0, sd0[i]); + out += student_t_lpdf(sd0[i] | blkparm1[i], 0, blkparm2[i]); + } + } + return out; + } } data { // see p. 2 https://books.google.com/books?id=9AC-s50RjacC @@ -602,6 +663,10 @@ data { int do_test; // should we do everything in generated quantities? array[Np] vector[multilev ? p_tilde : p + q - Nord] YXbar; // sample means of continuous manifest variables array[Np] matrix[multilev ? (p_tilde + 1) : (p + q - Nord + 1), multilev ? (p_tilde + 1) : (p + q - Nord + 1)] S; // sample covariance matrix among all continuous manifest variables NB!! multiply by (N-1) to use wishart lpdf!! + int npriblks; // how many blocks of parameters for prior distributions? + int priblklen; // max number of parameters in a block + vector[npriblks] blkparm1; // parameters of block priors + vector[npriblks] blkparm2; array[sum(nclus[,2])] int cluster_size; // number of obs per cluster array[Ng] int ncluster_sizes; // number of unique cluster sizes @@ -631,7 +696,6 @@ data { vector[multilev ? sum(ncluster_sizes) : Ng] log_lik_x; // ll of fixed x variables by unique cluster size vector[multilev ? sum(nclus[,2]) : Ng] log_lik_x_full; // ll of fixed x variables by cluster - /* sparse matrix representations of skeletons of coefficient matrices, which is not that interesting but necessary because you cannot pass missing values into the data block of a Stan program from R */ @@ -645,6 +709,8 @@ data { int len_lam_y; // number of free elements minus equality constraints array[len_lam_y] real lambda_y_mn; // prior array[len_lam_y] real lambda_y_sd; + array[len_lam_y] int lambda_y_pri; + array[len_lam_y] int lambda_y_blk; // same things but for B int len_w4; @@ -657,6 +723,8 @@ data { int len_b; array[len_b] real b_mn; array[len_b] real b_sd; + array[len_b] int b_pri; + array[len_b] int b_blk; // same things but for diag(Theta) int len_w5; @@ -726,6 +794,8 @@ data { int len_nu; array[len_nu] real nu_mn; array[len_nu] real nu_sd; + array[len_nu] int nu_pri; + array[len_nu] int nu_blk; // same things but for Alpha int len_w14; @@ -737,7 +807,9 @@ data { int len_alph; array[len_alph] real alpha_mn; array[len_alph] real alpha_sd; - + array[len_alph] int alpha_pri; + array[len_alph] int alpha_blk; + // same things but for Tau int len_w15; array[Ng] int wg15; @@ -1170,6 +1242,7 @@ parameters { vector[len_free[13]] Nu_free; vector[len_free[14]] Alpha_free; vector[len_free[15]] Tau_ufree; + vector[npriblks] sd0; // shrink_t parameters vector[Noent] z_aug; //augmented ordinal data vector[len_free_c[1]] Lambda_y_free_c; @@ -1550,11 +1623,8 @@ model { // N.B.: things declared in the model block do not get saved in the outp } } - /* prior densities in log-units */ - target += normal_lpdf(Lambda_y_free | lambda_y_primn, lambda_y_sd); - target += normal_lpdf(B_free | b_primn, b_sd); - target += normal_lpdf(Nu_free | nu_primn, nu_sd); - target += normal_lpdf(Alpha_free | alpha_primn, alpha_sd); + /* prior densities in log-units, first for unbounded parameters that could have shrinkage priors */ + target += eval_priors(Lambda_y_free, B_free, Nu_free, Alpha_free, sd0, lambda_y_primn, lambda_y_sd, lambda_y_pri, lambda_y_blk, b_primn, b_sd, b_pri, b_blk, nu_primn, nu_sd, nu_pri, nu_blk, alpha_primn, alpha_sd, alpha_pri, alpha_blk, len_free, blkparm1, blkparm2, priblklen, npriblks); target += normal_lpdf(Tau_ufree | tau_primn, tau_sd); target += normal_lpdf(Lambda_y_free_c | lambda_y_primn_c, lambda_y_sd_c); From 283e804d8c8675290d4350ec6d15cca84d7e6e9d Mon Sep 17 00:00:00 2001 From: ecmerkle Date: Wed, 18 Sep 2024 15:29:58 -0500 Subject: [PATCH 8/8] add normalizing constant for truncation of shrinkage prior --- inst/stan/stanmarg.stan | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/inst/stan/stanmarg.stan b/inst/stan/stanmarg.stan index 925ae394..17254626 100644 --- a/inst/stan/stanmarg.stan +++ b/inst/stan/stanmarg.stan @@ -614,7 +614,7 @@ functions { // you can use these in R following `rstan::expose_stan_functions("f } } out += normal_lpdf(parvec[1:(npars - 1)] | 0, sd0[i]); - out += student_t_lpdf(sd0[i] | blkparm1[i], 0, blkparm2[i]); + out += student_t_lpdf(sd0[i] | blkparm1[i], 0, blkparm2[i]) - log(.5); // left-truncated at 0 } } return out;