Please share your Stan program and accompanying data if possible.
When including Stan code in your post it really helps if you make it as readable as possible by using Stan code chunks (```stan) with clear spacing and indentation. For example, use
data {
// ----------------------------
// Basic dimensions
// ----------------------------
int<lower=1> L; // number of areas
vector[L] c; // lower-bound c_i
real<lower=0> a; // sum(theta_i) < a
// Sample sizes
int<lower=2> n[L];
// Direct estimates
real theta_hat[L]; // sample mean in each county
real S2[L]; // sample variance in each county
// Design for mean model
int<lower=1> pBeta;
matrix[L, pBeta] X_beta;
// Design for variance model
int<lower=1> pGamma;
matrix[L, pGamma] X_gamma;
}
parameters {
// ----------------------------
// Hyperparameters
// ----------------------------
vector[pBeta] beta; // regression coeffs for mean
real<lower=0> delta2; // variance for mean prior
vector[pGamma] gamma; // regression coeffs for log-precision
real<lower=0> alpha; // shape param for tau
// ----------------------------
// Area-level parameters
// ----------------------------
// tau_i = 1/sigma_i^2 > 0
vector<lower=0>[L] tau;
// theta_i is unconstrained here; I will do the constraints in model block
vector[L] theta;
}
transformed parameters {
// sigma_i^2 = 1 / tau_i
vector[L] sigma2;
for (i in 1:L) {
sigma2[i] = 1 / tau[i];
}
}
model {
// =================================================================
// 0) Enforce constraints: theta_i >= c_i and sum(theta_i) < a
// =================================================================
for (i in 1:L) {
if (theta[i] < c[i]) {
target += negative_infinity();
}
}
if (sum(theta) >= a) {
target += negative_infinity();
}
// =================================================================
// 1) Hyperpriors (weak or vague) on beta, delta^2, gamma, alpha
// =================================================================
// e.g. wide normal for beta
target += normal_lpdf(beta | 0, 100);
// delta^2 ~ InvGamma(0.001, 0.001)
target += inv_gamma_lpdf(delta2 | 0.001, 0.001);
// wide normal for gamma
target += normal_lpdf(gamma | 0, 100);
// alpha ~ Gamma(0.5, 0.5)
target += gamma_lpdf(alpha | 0.5, 0.5);
// =================================================================
// 2) Priors: tau_i = sigma_i^-2 ~ Gamma(alpha/2, [alpha * exp(-x_i gamma)]/2)
// =================================================================
for (i in 1:L) {
real rate_i = 0.5 * alpha * exp(-dot_product(X_gamma[i], gamma));
target += gamma_lpdf(tau[i] | alpha/2, rate_i);
}
// =================================================================
// 3) Priors: theta_i ~ Normal( X_beta[i]*beta, delta^2 )
// =================================================================
for (i in 1:L) {
target += normal_lpdf(theta[i] | dot_product(X_beta[i], beta), sqrt(delta2));
}
// =================================================================
// 4) Likelihood: Double Shrinkage
// (a) theta_hat[i] ~ Normal(theta[i], sigma2[i])
// (b) (n_i - 1)*S2[i]/sigma2[i] ~ Gamma( shape=(n_i-1)/2, rate=1/2 )
// =================================================================
for (i in 1:L) {
// (a) direct estimate
target += normal_lpdf(theta_hat[i] | theta[i], sqrt(sigma2[i]));
// (b) sample variance
// (n_i-1)*S2[i]/sigma2[i] ~ Chi^2_{(n_i-1)} => Gamma( (n_i-1)/2, 1/2 )
real shape_i = 0.5 * (n[i] - 1);
real rate_i = 0.5;
target += gamma_lpdf((n[i] - 1)*S2[i]/sigma2[i] | shape_i, rate_i);
}
}
2) FITTING THE MODEL IN STAN
=====================================================================
pBeta ← 1
X_beta ← matrix(1, nrow=L, ncol=1)
pGamma ← 1
X_gamma ← matrix(1, nrow=L, ncol=1)
stanData ← list(
L = L,
c = simData$c_i,
a = a,
n = simData$n_i,
theta_hat= simData$theta_hat,
S2 = simData$S2,
pBeta = pBeta,
X_beta = X_beta,
pGamma = pGamma,
X_gamma = X_gamma
)
init_fun ← function(chain_id = 1) {
leftover ← stanData$a - sum(stanData$c)
leftover_portion ← leftover * 0.5 / stanData$L
init_theta ← stanData$c + leftover_portion
list(
theta = init_theta,
tau = rep(1, stanData$L),
beta = as.array(rep(0, stanData$pBeta)), # force dimension [1]
gamma = as.array(rep(0, stanData$pGamma)), # likewise
alpha = 1,
delta2 = 1
)
}
fit ← stan(
file = “nandram_double_shrinkage.stan”,
data = stanData,
chains = 2,
iter = 10000,
warmup = 5000,
seed = 123,
control = list(adapt_delta = 0.999, max_treedepth = 20),
init = init_fun # pass the function here
)
I am getting these warning after running and I don’t know why:
Warning messages:
1: There were 10000 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.
2: Examine the pairs() plot to diagnose sampling problems
3: The largest R-hat is 2.64, indicating chains have not mixed.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#r-hat
4: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#bulk-ess
5: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#tail-ess