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