hierarchical logistic low speed

[Roula Nezi]

Dear Stan Users, 

I have another Stan newbies question. I am trying to fit a hierarchical logistic regression,  following the "election88_full.stan" example, but it takes more than 2 hours for iter=100 and chains=2. Most probably something is wrong with my syntax, but i cant figure out what is it. Any ideas what I might be doing wrong?  My code is below. 

Thank you!  

micro_code <- '

data {

int<lower=0> N; 

int<lower=0> n_cnt; 

int<lower=3> n_g;

int<lower=5> n_lag;

int<lower=0, upper=n_lag> lag [N];

int<lower= 1 > n_gdp;

int<lower=0,upper=n_cnt> cnt[N];

int<lower=0, upper=n_g>r0g [N];

int<lower=1, upper=34>rgdp[N];

int<lower=0,upper=97> age [N];

vector<lower=0,upper=1> [N] gender;

int<lower=1,upper=10>  income [N];

int<lower=1,upper=3>  class [N];

int<lower=1,upper=7>  lr [N];

int<lower=0,upper=1> y[N];

} 

parameters {

vector[n_lag] a;

vector[n_g] b;

vector[n_gdp] c;

vector[n_cnt] d;

vector[7] beta;

real<lower=0,upper=100> sigma_a;

  real<lower=0,upper=100> sigma_b;

  real<lower=0,upper=100> sigma_c;

 real<lower=0,upper=100> sigma_d;

}

transformed parameters {

vector[N] y_hat;

for (i in 1:N)

y_hat[i] <- beta[1] + beta[2] * class[i] + beta[3] * gender[i]

+ beta[4] * income[i] + beta[5] * lr[i] 

+ beta[6] * age[i] + a[lag[i]]+ b[r0g[i]]+c[rgdp[i]]+d[cnt[i]];

} 

model {

a ~ normal (0, sigma_a);

b ~ normal (0, sigma_b);

c ~ normal (0, sigma_c);

d ~ normal (0, sigma_d);

beta ~ normal(0, 100);

y ~ bernoulli_logit(y_hat);

}

'

micro.result <- stan(model_code=micro_code, data=dataList.1, iter=100, chains=2)

[Ben Goodrich]

It looks mostly fine. The normal(0,100) priors are probably overly diffuse. You could use matrix algebra instead of that scalar thing for yhat. But the main thing I think is to declare and define yhat inside the model block rather than the transformed parameters block. All those yhats get written to the disk after the warmup period which involves a lot of I/O for something that you probably don't care about and could recreate in R if you did care about them.

Ben

[Michael Betancourt]

Putting real priors on the standard deviations and employing a

non-centered parameterization of y_hat will also help.

--
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.

[Lisa Mu]

For a random intercept model  : logit (Y.ij) = beta %*% X.ij +  gamma.i , you may try  :   

1)   In data statements, a matrix to store individual level predictor variables  e.g age, plus a column vector '1' for intercept. Then you could use MVN as prior distribution for beta. To make coding easy,  data  had better be sorted by the order of groups ( e.g  all individuals from group 1 were followed by all individuals from group 2 ...)

2)  In data statements, declare a vector with length = N, to store group index for each individual , following the same order in 1) , e.g vector (1,2, ....I)

3) In model statement, a) sample gamma.i for each group ; b) sample beta ; c) sample bernoulli for each individual based on a) and b)

You could refer to Stan manual about hierarchical  model where Y.ij is normal.