problems estimating scale of multivariate prior for hierarchical regression

216 views
Skip to first unread message

Guido Biele

unread,
Mar 4, 2015, 10:18:51 AM3/4/15
to stan-...@googlegroups.com
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    
HLRegressionP.stan
test1.rdump
accept_treedepth_etc_modified_sampler_paras.png
tau_modified_sampler_paras.png
accept_treedepth_etc_default_sampler_paras.png
tau_default_sampler_paras.png

Daniel Lee

unread,
Mar 4, 2015, 10:34:20 AM3/4/15
to stan-...@googlegroups.com
Hi Guido,

It might help to also look at tau_unif. Could you plot that also?
[edit: I just realized that tau was a transformation of tau_unif]

On Wed, Mar 4, 2015 at 10:18 AM, Guido Biele <guido...@gmail.com> wrote:
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?

The defaults certainly look better than the modified. In the modified one, your marginal posterior estimates for your individual taus will be off.
 
2) Are the modified sampler parameters I used OK, or should I use different parameter values?

I'll like Mike chime in.
 
3) Does the model look all right?

What's the prior you're trying to induce on tau? It looks like that's what you're interested in, but you've set up tau_unif as between 0 and pi/2 and then 'tau <- tan(tau_unif)'. It's a little too early for me to reason through what that means. Perhaps you can explain?
 

--
You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Rob Trangucci

unread,
Mar 4, 2015, 1:29:13 PM3/4/15
to stan-...@googlegroups.com
Guido,

Try tightening the priors on tau. Instead of cauchy(0,1), use normal(0,1). 

Rob
--
You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
<HLRegressionP.stan>
<test1.rdump>
<accept_treedepth_etc_modified_sampler_paras.png>
<tau_modified_sampler_paras.png>
<accept_treedepth_etc_default_sampler_paras.png>
<tau_default_sampler_paras.png>

Guido Biele

unread,
Mar 4, 2015, 2:17:15 PM3/4/15
to stan-...@googlegroups.com
Thanks for your quick response Daniel!
I am attaching the plots for tau_unif and for tau for a new run.

the prior for tau_unif was on purpose a uniform. I omitted the sampling statement because the variable will be samples as a from a uniform with the boundaries specified during variable declaration. For the sake o clarity  I should have included the sampling statement.

following the stan reference i replaced 
chauchy(0,1)
with 
tan_unif<lower=0,uper=pi/2>
tau <- tan(tan_unif)

(the example in stan reference is to substitute beta ~ chauchy(mu,tau) with beta_unif<lower = -pi()/2,upper = pi(2)> and beta <- mu + tau*tan(beta_unif)

cheers - guido

Guido Biele

unread,
Mar 4, 2015, 3:38:53 PM3/4/15
to stan-...@googlegroups.com
Thanks for the pointer Rob!
tau ~ normal(0,1) works much better.
reading about stan sampler parameters today made me overlook this much simpler solution ...
Best - Guido
Reply all
Reply to author
Forward
0 new messages