Good afternoon dear Community,
I have been learning how to use Stan for bayesian statistics for almost 1 month. In this occasion, I am trying to do a Ordered Logistic Model of Item Response Model for a test with different number of categories per item. When I run my code I have these warning:
fitasime<- stan(file = “asime.stan”, data = dat,iter = 10, chains = 1)
SAMPLING FOR MODEL ‘asime’ NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.004 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 40 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: WARNING: No variance estimation is
Chain 1: performed for num_warmup < 20
Chain 1:
Chain 1: Iteration: 1 / 10 [ 10%] (Warmup)
Chain 1: Iteration: 2 / 10 [ 20%] (Warmup)
Chain 1: Iteration: 3 / 10 [ 30%] (Warmup)
Chain 1: Iteration: 4 / 10 [ 40%] (Warmup)
Chain 1: Iteration: 5 / 10 [ 50%] (Warmup)
Chain 1: Iteration: 6 / 10 [ 60%] (Sampling)
Chain 1: Iteration: 7 / 10 [ 70%] (Sampling)
Chain 1: Iteration: 8 / 10 [ 80%] (Sampling)
Chain 1: Iteration: 9 / 10 [ 90%] (Sampling)
Chain 1: Iteration: 10 / 10 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 0.114 seconds (Warm-up)
Chain 1: 0.099 seconds (Sampling)
Chain 1: 0.213 seconds (Total)
Chain 1:
Warning messages:
1: There were 5 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
2: Examine the pairs() plot to diagnose sampling problems
fitasime2<- stan(fit =fitasime , data =dat, iter = 4000, chains = 4, control = list(max_treedepth = 12))
Error in unserialize(socklist[[n]]) : error al leer desde conexión
Error in serialize(data, node$con, xdr = FALSE) :
error al escribir en una conexión
I tried to increase the adapt_delta to 0.99 but does not improve. More over, several time after this messages, appear “R session is aborted. R encoutered a fatal error” and then the R program closes .
I attache my data and my code and will be so gratefull your help.
[datosoriginal2.csv (7.8 KB) ](http://)
data{
int<lower=1> N; // number of answers
int<lower=10> N_fami; // number of entrepreneur
int<lower=2> N_item; // number of items
int<lower=2, upper=5>N_cat[N];// cut points
int<lower=1, upper=4>N_ord_cat[N];// order question with the same cut points
int<lower=1> y[N]; // respuestas; y[N], N-th answers
int<lower=1,upper=N_fami> famiempresas[N];// respuestas del empresario[N]
int<lower=1,upper=N_item> item[N]; // número de items [N]
}//end data
parameters{
vector[N_fami] theta; // latent traits of microempresarios
ordered[2] beta_1[2];// coefficients of predictors (cut p)
ordered[3] beta_2[4];// coefficients of predictors (cut p)
ordered[4] beta_3;// coefficients of predictors (cut p)
ordered[5] beta_4;// coefficients of predictors (cut p)
real mu_beta; // mean of the beta-parameters: difficulty of the item
real<lower=0> sigma_beta; //sd of the prior distributions of category difficulty
vector<lower=0>[N_item] alpha; // parámetro de discriminación
real pho_theta; // parámetro de asimétria del trazo latente
}//end parameters
model{
theta ~ skew_normal(0,1, pho_theta);// latent traits
for(j in 1:N_item){
alpha[j] ~ normal(0,1); // prior discriminación
}
for(j in 1:2){
beta_1[j]~ normal(mu_beta,sigma_beta);// prior of the item dificultad
}
for(j in 1:4){
beta_2[j]~ normal(mu_beta,sigma_beta);// prior of the item dificultad
}
beta_3~ normal(mu_beta,sigma_beta);// prior of the item dificultad
beta_4~ normal(mu_beta,sigma_beta);// prior of the item dificultad
pho_theta ~ normal(0,1); // hyper prior for theta, parámetro de asimétria.
mu_beta ~ normal(0,2); // hyper prior for mu_beta
sigma_beta ~ cauchy(0,2); // hyper prior for sigma_beta
//likelihood
for(n in 1:N){
if (N_cat[n]==2){
y[n] ~ ordered_logistic(alpha[item[n]]*theta[famiempresas[n]],beta_1[N_ord_cat[n]]);
}
if (N_cat[n]==3){
y[n] ~ ordered_logistic(alpha[item[n]]*theta[famiempresas[n]],beta_2[N_ord_cat[n]]);
}
if (N_cat[n]==4){
y[n] ~ ordered_logistic(alpha[item[n]]*theta[famiempresas[n]],beta_3);
}
if (N_cat[n]==5){
y[n] ~ ordered_logistic(alpha[item[n]]*theta[famiempresas[n]],beta_4);
}
}
}//end model
R code
rm(list=ls(all=TRUE))
datos=read.table(“datosoriginal2.csv”, header=T, sep=“;”)
attach(datos)
names(datos)
#install.packages(‘tidyr’) datos en forma long
library(tidyr)
df<- gather(datos, Item, Valor, 2:9)
attach(df)
names(df)
head(df)
dim(df)
dim(datos)
rstan
library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
parameters
NN = nrow(df) #3072; num. registros
NP = nrow (datos) #384; num. famiempresas
NI = ncol(datos)-1 #8; num. items
cut point
datos1=datos[,-1]
NC=rep(NA, NI)
for (i in 1:NI){
NC[i]=(dim(table(datos1[[i]]))-1)
}
n_cat_fami=rep(c(2,3,4,5),c(768,1536,384,384))
n_orden_cat=rep(c(1,2,1,2,3,4,1,1),c(384, 384,384, 384,384, 384,384, 384))
library(plyr)
df$Item ← revalue(df$Item, c(“X1”=1))
df$Item ← revalue(df$Item, c(“X2”=2))
df$Item ← revalue(df$Item, c(“X3”=3))
df$Item ← revalue(df$Item, c(“X4”=4))
df$Item ← revalue(df$Item, c(“X5”=5))
df$Item ← revalue(df$Item, c(“X6”=6))
df$Item ← revalue(df$Item, c(“X7”=7))
df$Item ← revalue(df$Item, c(“X8”=8))
i=as.integer(df[,2])
For stan
dat = list( famiempresas=df[,1], item=i, y= df[,3], N=NN, N_fami=NP,
N_item=NI, N_cat=n_cat_fami, N_ord_cat=n_orden_cat)
attach(dat)
str(dat)
head(dat)
tail(item)
fitasime<- stan(file = “asime.stan”, data = dat,iter = 10, chains = 1)
fitasime2<- stan(fit =fitasime , data =dat, iter = 4000, chains = 4, control = list(max_treedepth = 12))