Data model code running too slowly

Question

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)

Answer 1

You made a typo. The line

responses ~ bernoulli_logit(theta);

should be outside the for loop in the model block.