rStan "hangs" after completing multivariate probit with 19 items - fine with 17

Modeling
1

Edit: issue solved. I create a vector on the last block containing each possible response pattern and as I include more items that vector expands in length like the grains of rice on the chess board squares in the old story. Oops.

I built a statistical model around the multivariate probit example in the Stan manual.

I have used this to estimate the correlation structure of test items answered by school children and draw posterior predictive samples. The model works as intended in multiple data sets, but it seems to hang after finishing the chains if there are too many items. When I recently tried to estimate on a dataset with all 19 binary test items (in a data file with about 260 resondents) the model finishes each chain but never returns to the R system prompt - even when I left it running overnight.

The weird thing:

  • Using only five items from this datafile takes about 50 seconds to finish and has no hanging issues.
  • Using only 10 of the 19 items take 220 seconds to run - no hanging
  • Only the last 10, estimates in 115 seconds, no issues with hanging.
  • Using the first 15 items takes 330 seconds - there is a perceptible “lag” after finishing - after which it tells me Tail effective sizes and Effective samples are low and then finishes and returns to the R prompt
  • Using the first 17 takes 839 seconds and a perceptible lag of about 10-15 minutes after which it notes the ESS and TSS issues
  • And - as noted initially - using 19 items takes 2 hours to finish and a lag that I have yet to see the end of

My practical solution for now is just to leave out the (least interesting) subset of the items to get the item count reduced, but if anyone has insights into why things seem to change so dramatically between 17 and 19 items please let me know what it might reflect and how to address it.

One other thing: When estimating I leave out most of the variables from the returned file to reduce file size:

stan(
file = here(“stan/binary - multivariate probit.stan”),
data = list(
obs_n = nrow(temp),
items_n = ncol(temp),
responses = temp),
pars = c(“z_pos”,
“z_neg”,
“z”),
include = F,
refresh = 50)

Model code:

functions {
  // A helper function to sum all elements of a 2D int array
  int sum2d(int[,] a) {
    int s = 0;
    for (i in 1:size(a)) {
      s += sum(a[i]);
    }
    return s;
  }
}

data {
  int<lower=0> obs_n;                      // Number of observations
  int<lower=0> items_n;                    // Number of items
  int<lower=0, upper=1> responses[obs_n, items_n];  // 0/1 response matrix
}

transformed data {
  int total_ones = sum2d(responses);
  int total_zeros = obs_n * items_n - total_ones;
  int obs_pos[total_ones];
  int item_pos[total_ones];
  int obs_neg[total_zeros];
  int item_neg[total_zeros];
  int sequence_count = 1;

  {
    int pos_idx = 1;
    int neg_idx = 1;
    for (n in 1:obs_n) {
      for (j in 1:items_n) {
        if (responses[n,j] == 1) {
          obs_pos[pos_idx]   = n;
          item_pos[pos_idx]  = j;
          pos_idx += 1;
        } else {
          obs_neg[neg_idx]   = n;
          item_neg[neg_idx]  = j;
          neg_idx += 1;
        }
      }
    }
  }

  for (i in 1:items_n) {
    sequence_count *= 2;
  }
}

parameters {
  cholesky_factor_corr[items_n] L_Omega;
  vector[items_n] alpha;
  vector<lower=0>[total_ones] z_pos;
  vector<upper=0>[total_zeros] z_neg;
}

transformed parameters {
  vector[items_n] z[obs_n];
  for (p in 1:total_ones) {
    z[obs_pos[p], item_pos[p]] = z_pos[p];
  }
  for (q in 1:total_zeros) {
    z[obs_neg[q], item_neg[q]] = z_neg[q];
  }
}

model {
  L_Omega ~ lkj_corr_cholesky(2);   
  alpha   ~ normal(0, 5);         
  z ~ multi_normal_cholesky(alpha,
  L_Omega);
}

generated quantities {
  corr_matrix[items_n] Omega; 
  int sequences[sequence_count] = rep_array(0, sequence_count);
  int sumscore[items_n + 1] = rep_array(0, items_n + 1);
  Omega = multiply_lower_tri_self_transpose(L_Omega);

  {
    vector[items_n] z_sim[obs_n];
    int temp_array[obs_n, items_n];

    for (n in 1:obs_n) {
      z_sim[n] = multi_normal_cholesky_rng(alpha, L_Omega);
    }

    // Convert those latent z_sim's to 0/1 by thresholding at 0
    for (n in 1:obs_n) {
      for (j in 1:items_n) {
        if (z_sim[n, j] > 0) {
          temp_array[n, j] = 1;
        } else {
          temp_array[n, j] = 0;
        }
      }
    }

    for (obs_i in 1:obs_n) {
      int temp_sumscore = sum(temp_array[obs_i, ]);
      int temp_sequence = 0;
      int factor = 1;

      for (item_i in 1:items_n) {
        if (temp_array[obs_i, item_i] == 1) {
          temp_sequence += factor;
        }
        factor *= 2;
      }

      sumscore[temp_sumscore + 1] += 1;
      sequences[temp_sequence + 1] += 1;
    }
  }
}