R constantly crashed after a few iterations


#1

Hi all,

I am totally new to Stan. So please forgive me if I ask some “dumb questions”.

I am recently trying to estimate a complex Item Response Theory (IRT) model. However, R just constantly crashed after several iterations (usually less than 10). I am not sure what’s going on. Here blow is my code.

BRES.Ordinal<-'

functions {
  
real RES(int y, int K, real theta_R, real theta_E, real theta_S, real alpha_R, real alpha_E, real alpha_S, vector beta_R, real beta_E, vector beta_S) {

////////////////////////////////////////////
///////////////////////////////////////////

    row_vector[K] prob_R;
    real prob_E;  
    matrix[K,K] prob_pi;  
    vector[K] prob_obs;     
    vector[K-1] omegaA;     

////////////////////////////////////////////
///////////////////////////////////////////

    prob_R[1] = 1-inv_logit(alpha_R*theta_R-beta_R[1]);
 
for (e in 2:(K-1)){
    prob_R[e]=inv_logit(alpha_R*theta_R-beta_R[e-1])-inv_logit(alpha_R*theta_R-beta_R[e]);
}
     prob_R[K]=inv_logit(alpha_R*theta_R-beta_R[K-1]);


////////////////////////////////////////////
///////////////////////////////////////////

    prob_E=inv_logit(alpha_E*theta_E-beta_E);

////////////////////////////////////////////
///////////////////////////////////////////

for (c in 1:(K-2)){
    omegaA[c]=exp(alpha_S*theta_S - beta_S[c]);  
}
    omegaA[K-1]=1;

for (h in 1:(K-1)){
    prob_pi[(K-1),h]=0;
    prob_pi[K,h]=0;
}
    prob_pi[(K-1),K]=1;
    prob_pi[K,K]=1;

for (X in 1:(K-2)){
  for (Y in 1:X){
    prob_pi[X,Y]=0;
  }
}

for (m in 1:(K-2)){
  for (n in (m+1):K){
    prob_pi[m,n]=omegaA[n-1]/sum(omegaA[m:(K-1)]);
  }
}
////////////////////////////////////////////
///////////////////////////////////////////

    prob_obs[1]=prob_R[1]*(1-prob_E);


for (f in 2:(K-1)){
    prob_obs[f]=prob_R[f]*(1-prob_E)+prob_E*(prob_R[1:(f-1)]*prob_pi[1:(f-1),f]);

}
    prob_obs[K]=prob_R[K]+prob_E*(prob_R[1:(K-1)]*prob_pi[1:(K-1),K]);


////////////////////////////////////////////
///////////////////////////////////////////

    return categorical_lpmf(y|prob_obs);
}
}

data {
  int<lower=1> K;               
  int<lower=1> I;               
  int<lower=1> J;               
  int<lower=1> N;               
  int<lower=0> N_mis;          
  int<lower=1, upper=I> II[N];  
  int<lower=1, upper=J> JJ[N];  
  int<lower=0, upper=6> y[N];   
  real alpha_S[I]; 
  vector[K-2] beta_S;
  vector[3] sigma;
  vector[3] theta_mu;
}

parameters {

  real<lower=0> alpha_R[I];
  real<lower=0> alpha_E[I];
  real beta_E[I];
  ordered[K-1] beta_R[I];
  matrix[J,3] theta;                       
  cholesky_factor_corr[3] theta_cor;                
}
transformed parameters{

  vector[J] theta_R;
  vector[J] theta_E;
  vector[J] theta_S;

  theta_R=theta[,1];
  theta_E=theta[,2];
  theta_S=theta[,3];
}

model {
    theta_cor ~ lkj_corr_cholesky(1);
    alpha_R[I] ~ lognormal(0,0.3);
    alpha_E[I] ~ lognormal(0,0.3);
    beta_E[I] ~ normal(0, 10);
    
    for (i in 1:I){
    beta_R[i,] ~ normal(0,10);
    }

for (q in 1:J){

    theta[q,] ~ multi_normal_cholesky(theta_mu,diag_pre_multiply(sigma,theta_cor));
}  

for (n in 1:N)
  target += RES(y[n],K,theta[JJ[n],1],theta[JJ[n],1],theta[JJ[n],3],alpha_R[II[n]],alpha_E[II[n]],alpha_S[II[n]],beta_R[II[n],],beta_E[II[n]],beta_S);
}
'

I also attached my code for generating data and fitting the model.

library(rstan)
library(edstan)
library(MASS)
library(psych)

K=5       # Number of response category 
JJ=100    # Number of respondent
II=10      # Number of items 

# covariance matrix among theta.R, theta.E and theta.S

cov<-matrix(c(1,0.5,0.5,
              0.5,1,0.5,
              0.5,0.5,1),byrow = T,nrow = 3)
cov2cor(cov)
theta<-mvrnorm(n=JJ,mu=c(0,0,0),Sigma=cov)
theta.R<-theta[,1]
theta.E<-theta[,2]
theta.S<-theta[,3]


# Simulated item parameters 

alpha.R<-matrix(rep(c(1,1.25,1.5,1.75,2),4))[1:II,]
alpha.E<-matrix(rep(c(0.5,1,1.5,1,0.5),4))[1:II,]
alpha.S<-matrix(rep(c(1,1,1,1,1,1,1,1,1,1),2))[1:II,]

beta.R<-matrix(rep(c(-1,-0.5,0.5,1),II),byrow = T,ncol = K-1)[1:II,]
beta.E<-matrix(rep(c(0.5,0.2,-0.2,-0.5),5))[1:II,]
beta.S<-matrix(rep(c(-1,-0.5,0),II),byrow = T,ncol = K-2)[1:II,]


# Data generation function 

GenerateRES<-function(K,JJ,II,theta.R,theta.E,theta.S,alpha.R,alpha.E,alpha.S,beta.R,beta.E,beta.S){
  
  ObservedResponse<-matrix(NA,nrow =JJ,ncol = II)  # observed response matrix 
  
  for (a in 1: JJ){
    
    for(b in 1:II){

      prob.obs<-numeric(K)     # probability of observed response       
      prob.R<-numeric(K)       # probability at the retrieval stage   
      prob.E<-numeric(K)       # probability at the editing stage   
      prob.pi<-matrix(NA,K,K)  # transition matrix
      omegaA<-numeric(K-1)       # overall propensity to select the Kth category 
      
      for (c in 1:(K-2)){
        omegaA[c]=exp(alpha.S[b]*theta.S[a]-beta.S[b,(c)])  # define the overall propensity to select the Kth category  
      }
        omegaA[K-1]=1

      
      # calculate probability at the retrieval stage 
      prob.R[1]=1/(1+exp(alpha.R[b]*theta.R[a]-beta.R[b,1]))
      for (e in 2:(K-1)){
        prob.R[e]=exp(alpha.R[b]*theta.R[a]-beta.R[b,e-1])/(1+exp(alpha.R[b]*theta.R[a]-beta.R[b,e-1]))-exp(alpha.R[b]*theta.R[a]-beta.R[b,e])/(1+exp(alpha.R[b]*theta.R[a]-beta.R[b,e]))
      }
      prob.R[K]=exp(alpha.R[b]*theta.R[a]-beta.R[b,K-1])/(1+exp(alpha.R[b]*theta.R[a]-beta.R[b,K-1]))
      
      # calculate probability at the editing stage
      
        prob.E=exp(alpha.E[b]*theta.E[a]-beta.E[b])/(1+exp(alpha.E[b]*theta.E[a]-beta.E[b]))
      
      # Define the transition matrix (category attractiveness)
      prob.pi[(K-1):K,1:(K-1)]=0
      prob.pi[(K-1):K,K]=1

      for (d in 1:(K-2)){
        prob.pi[d,1:d]=0
      }
      for (m in 1:(K-2)){
        for (n in (m+1):K)  {
          prob.pi[m,n]=omegaA[n-1]/sum(omegaA[m:(K-1)])
        }
      }
  
      prob.obs[1]=prob.R[1]*(1-prob.E)
      
      for(f in 2:(K-1)){
        
        prob.obs[f]<-prob.R[f]*(1-prob.E)+prob.E*sum(prob.R[1:(f-1)]*prob.pi[1:(f-1),f])
      }
      
      prob.obs[K]<-prob.R[K]+prob.E*sum(prob.R[1:(K-1)]*prob.pi[1:(K-1),K])
      
      
      ObservedResponse[a,b]<-sample(K,size = 1,replace = T,prob = prob.obs)
    }
  }
  return(ObservedResponse)

}

y.data<-GenerateRES(K,JJ,II,theta.R,theta.E,theta.S,alpha.R,alpha.E,alpha.S,beta.R,beta.E,beta.S)


Data<-irt_data(response_matrix =y.data)

I<-dim(y.data)[1]
J<-dim(y.data)[2]
data_list<-list(
  I=Data$I,
  J=Data$J,
  N=Data$N,
  II=Data$ii,
  JJ=Data$jj,
  K=max(Data$y),
  N_mis=(Data$I*Data$J-Data$N),
  y=Data$y,
  sigma=c(1,1,1),
  theta_mu=c(0,0,0),
  alpha_S=rep(1,II),
  beta_S=c(-0.5,0,0.5))


## Estimate the model with simulated data 

rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores() - 1)


BRES<-stan(model_code=BRES.Ordinal,data=data_list,iter=10,chains = 1,cores = 1,control = list(adapt_delta=0.90,max_treedepth=15))

I greatly appreciate any suggestion.

Best,
Bobby