Hi avehtari,

Thanks for your reply.

The code below is for the mixture 2PL IRT model, for the 2-group solution, where “K = 2” in the code means fit the 2-group solution.

y = resp

K = 2

N = nrow(y)

J = ncol(y)

dir_alpha = c(1 , 1)

modelString="data {

int<lower=1> K ;

int<lower=1> N ;

int<lower=1> J ;

int<lower=0,upper=1> y[N,J] ;

vector<lower=0>[K] dir_alpha ;

}

parameters {

simplex[K] pi ;

vector[N] theta ;

ordered[K] mu ;

ordered[K] beta[J] ;

vector<lower=0>[K] alpha[J] ;

}

model {

real lps[K] ;

real p[K] ;

real lpth[K] ;

for (k in 1:K){

mu[k] ~ normal(0, 1) ;

}

pi ~ dirichlet(dir_alpha) ;

for (j in 1:J){

for (k in 1:K){

beta[j,k] ~ normal (0,1) ;

alpha[j,k] ~ normal(0,1) ;

}}

for (i in 1:N){

for (k in 1:K){

lpth[k]= log(pi[k]) + normal_lpdf(theta[i] | mu[k], 1) ;

}

target += log_sum_exp(lpth ) ;

}

for (i in 1:N){

for (j in 1:J){

for (k in 1:K){

p[k] = inv_logit( alpha[j,k]*(theta[i]-beta[j,k]) ) ;

lps[k]= log(pi[k]) + bernoulli_lpmf(y[i,j] | p[k]) ;

}

target += log_sum_exp(lps) ;

}}

}

generated quantities{

real p[K] ;

vector[J]log_lik[N] ;

for (i in 1:N){

for (j in 1:J){

for (k in 1:K){

p[k] = inv_logit( alpha[j,k]*(theta[i]-beta[j,k]) ) ;

log_lik[i,j] = log( exp(log(pi[k]) + bernoulli_lpmf(y[i,j] | p[k])) ) ;

}}}

}

"

stanDso <- stan_model( model_code=modelString )

Mix_2C = sampling( object=stanDso, data=list(K,N,J, y, dir_alpha ), warmup = 3000, iter = 5000, chains = 3 )

The results is as follows:

- one-group solution, K =1 (which means fitting the conventional 2PL IRT)

loo1

Computed from 10000 by 30000 log-likelihood matrix

```
Estimate SE
```

## elpd_loo -15151.2 90.5

p_loo 917.0 11.8

looic 30302.4 181.1

Monte Carlo SE of elpd_loo is 0.4.

Pareto k diagnostic values:

Count Pct. Min. n_eff

(-Inf, 0.5] (good) 29991 100.0% 1875

(0.5, 0.7] (ok) 9 0.0% 1155

(0.7, 1] (bad) 0 0.0%

(1, Inf) (very bad) 0 0.0%

All Pareto k estimates are ok (k < 0.7).

See help(‘pareto-k-diagnostic’) for details.

- The two-group solution (K = 2) (This is the correct solution)

loo2

Computed from 6000 by 30000 log-likelihood matrix

```
Estimate SE
```

## elpd_loo -43311.6 235.7

p_loo 7969.5 141.4

looic 86623.3 471.4

Monte Carlo SE of elpd_loo is NA.

Pareto k diagnostic values:

Count Pct. Min. n_eff

(-Inf, 0.5] (good) 24692 82.3% 296

(0.5, 0.7] (ok) 2538 8.5% 62

(0.7, 1] (bad) 1831 6.1% 11

(1, Inf) (very bad) 939 3.1% 2

See help(‘pareto-k-diagnostic’) for details.

- The 3-group solution (K = 3)

loo3

Computed from 6000 by 30000 log-likelihood matrix

```
Estimate SE
```

## elpd_loo -50537.5 301.5

p_loo 11095.7 192.6

looic 101074.9 603.1

Monte Carlo SE of elpd_loo is NA.

Pareto k diagnostic values:

Count Pct. Min. n_eff

(-Inf, 0.5] (good) 23163 77.2% 465

(0.5, 0.7] (ok) 2615 8.7% 102

(0.7, 1] (bad) 2211 7.4% 13

(1, Inf) (very bad) 2011 6.7% 1

See help(‘pareto-k-diagnostic’) for details.