The model takes a long time to run. I have tried to vectorise it as far as I can, but it's still taking about a day to compute. I'm wondering if it's possible to speed it up.

Maybe the problem isn't even with the code, but with the model being improper/bad priors/...? Or maybe the data set is simply very large and supposed to take hours to compute. Either way, for now, I dread adding another predictor.

# data is in wide format, one row with one outcome and a range of predictors

set_cppo(mode = "fast")

code ="                     // the model itself
data        {
       int<lower=1> nrow;   // 72000 (3 items per W*N*R)
       int W;               // n words; 24
       int N;               // n subjects; 100
       int R;               // n raters; 10
       int<lower=1,upper=nrow> ind[nrow]; // this vector is used to index fields per row
       int<lower=0,upper=1> responses[nrow];  // outcomes
       int subjects[nrow];
       int words[nrow];
       int sex[nrow];
       int languages[nrow];
       int raters[nrow];

parameters  {
real intercept;               // "fixed"
real beta_language;
real beta_sex;
real interaction_s_l;
real beta_rater[R];           // "random"
real beta_word[W];
real beta_subject[N];

model       {
real theta[nrow];

 beta_word ~ normal(0,100);         // priors
 beta_subject ~ normal(0,100);
 beta_language ~ normal(0,100);
 beta_sex ~ normal(0,100);
 interaction_s_l ~ normal(0,100);

for (i in 1:nrow){                    // for each row, construct theta
theta[i] <- intercept
          + beta_sex * sex[ind[i]]
          + beta_language * languages[ind[i]]
          + interaction_s_l * languages[ind[i]] * sex[ind[i]]
          + beta_subject[subjects[ind[i]]]
          + beta_rater[raters[ind[i]]];}

responses ~ bernoulli_logit(theta);   // sample

generated quantities {               // get odds ratios
real odds_intercept;
real odds_sex;
real odds_language;
real odds_interaction;

odds_intercept <- exp(intercept);
odds_sex <- exp(beta_sex);
odds_language <- exp(beta_language);
odds_interaction <- exp(interaction_s_l);

fit <- stan(model_code = code, data = dat,     // run the model
            iter = 4000, chains = 3)
  • \$\begingroup\$ Could you describe the purpose of this code? \$\endgroup\$
    – Jamal
    Commented Jun 29, 2014 at 22:05
  • \$\begingroup\$ It's supposed to estimate the odds ratios for my predictors, with "sex" and "language" as well as their interaction as 'fixed' and "word" and "subject" as 'random' effects. Does it not do that? I have added some comments. \$\endgroup\$
    – jona
    Commented Jun 29, 2014 at 22:26
  • \$\begingroup\$ Oh no, I'm asking because it'll help give a more descriptive title. You may make that change now. \$\endgroup\$
    – Jamal
    Commented Jun 29, 2014 at 22:35
  • \$\begingroup\$ This is a bit late, but I'd suggest taking a look on the Stan-users mailing list where these kinds of questions get asked daily. \$\endgroup\$ Commented Jul 23, 2014 at 17:59
  • \$\begingroup\$ I happened to come across this again and finally got to take a look at it. One easy transformation you can do is use standard normal priors, which you then rescale in the transformed parameters block. e.g. transformed parameters{ beta_word <- beta_word_standardized * 10; } then model { beta_word_standardized ~ normal(0,1); }. \$\endgroup\$ Commented Aug 17, 2014 at 19:16

1 Answer 1


You made a typo. The line

responses ~ bernoulli_logit(theta);

should be outside the for loop in the model block.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.