Hello,
I am trying to estimate an I believe straight forward hierarchical logistic regression, but the traces for some parameters (mainly tau, the coefficient scales of the multivariate hierarchical prior ) look “spiky” (see attached figure ) and I can’t seem to get the behave better. At times, these parameters also get stuck.
The essential parts of the model are below.
Here is what I have done so far to get nicer sampling behavior
- Reparameterized the chauchy distribution for tau (the scale of the multivariate priors) as described on pages 198-190 in the stan reference
- Increased the sampler parameter max_treedepth to 15
- Increased target acceptance rate to .9 and jittered stepsize (.5)
All this has all little effect (If anything, the results with the modified parameters are worse).
I am attaching traceplots and histograms of acceptance probability, step size, treedepth.
So I basically have three questions now:
1) Does the sample trace look OK after all, and I am worrying for nothing?
2) Are the modified sampler parameters I used OK, or should I use different parameter values?
3) Does the model look all right?
The full model and a small data set are attached. (Full model and data are attached. Note that the full model has lots of code to scale variables and to prepare and do average predictive comparisons)
Cheers - Guido
############# essential sections of the model #################
parameters {
// covariances for regression model parameters for mean of negative binomial
// (from stan language reference)
matrix[K,M] z;
cholesky_factor_corr[K] L_Omega;
//vector<lower=0>[K] tau; // prior scale
vector<lower=0, upper=pi()/2>[K] tau_unif;
matrix[L,K] gamma; // group coeffs
}
transformed parameters {
// regression weights for mean of negative binomial model
vector<lower=0>[K] tau;
matrix[M,K] beta;
for (k in 1:K) tau[k] <- tan(tau_unif[k]);
beta <- gx * gamma + (diag_pre_multiply(tau,L_Omega) * z)';
}
model {
// log(mean) of negative binomial
vector[N] alpha;
// sample parameters
//tau ~ cauchy(0,2.5); # orig scale = 2.5
to_vector(z) ~ normal(0,1);
L_Omega ~ lkj_corr_cholesky(2);
to_vector(gamma) ~ normal(0,5); # orig scale = 5
alpha <- to_vector(sx * beta');
// increment binomial
y ~ bernoulli_logit(alpha);
}
########### model calls ################
called the model with
stanfit = stan(file = "stanmodels/HLRegP/HLRegressionP.stan",data = standata,chains = 1, iter = 500, control=list(refresh = 25))
and
stanfit_plus = stan(fit = stanfit,data = standata,chains = 4, iter = 500,control = list(max_treedepth = 15, adapt_delta = .85, stepsize_jitter = .25))
(I a longer warmup, i tried 750, does not make a difference)
data is in test1.rdump
########### R session info ################
> sessionInfo()
R version 3.1.2 (2014-10-31)
Platform: x86_64-w64-mingw32/x64 (64-bit)
locale:
[1] LC_COLLATE=Norwegian (Bokmål)_Norway.1252 LC_CTYPE=Norwegian (Bokmål)_Norway.1252
[3] LC_MONETARY=Norwegian (Bokmål)_Norway.1252 LC_NUMERIC=C
[5] LC_TIME=Norwegian (Bokmål)_Norway.1252
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] rstan_2.6.0 inline_0.3.13 Rcpp_0.11.4
loaded via a namespace (and not attached):
[1] codetools_0.2-9 stats4_3.1.2 tools_3.1.2