5.14 MCMC for DINA model in STAN

STAN is another very popular MCMC program and one can use rstan in R to run stan code. STAN uses different algorithms from JAGS and nimble. A limitation of STAN is that it cannot handle discrete variable directly.

Recall that in nimble and JAGS, one needs to specify the probability of obtaining a score for each student on each item, which serves as the likelihood.

In STAN, we have to marginalize the discrete variables.

Recall that the log likelihood of item response matrix $\mathbf{Y}$ of all $N$ students can be written by $$\begin{equation} \ell(\mathbf{Y}) = \log \prod\limits_{i = 1}^N L ({\mathbf{Y}_i}) = \sum\limits_{i = 1}^N \log L({\mathbf{Y}_i}) \end{equation}$$ where the marginalized log likelihood observing a response vector $\mathbf{Y}_i$ for student $i$ is $$\begin{equation} \log L({\mathbf{Y}_i})=\log \sum_cL( {\mathbf{Y}_i}|{\alpha_i=\alpha_c})p(\mathbf{\alpha}_c) \end{equation}$$

The code below, which was adopted from here, can be used.

Code

data{
  int<lower=1> J;         // # of items
  int<lower=1> N;         // # of respondents 
  int<lower=1> K;         // # of attributes
  int<lower=1> C;         // # of attribute profiles (latent classes) 
  matrix[N,J] y;          // response matrix
  matrix[C,J] eta;      // eta matrix
  matrix[J,K] Q; 
}

parameters{
  simplex[C] pi;                   // probabilities of latent class membership
  real<lower=0,upper=1> s[J];      // slip parameter   
  real<lower=0,upper=1> g[J];      // guess parameter
}

transformed parameters{
  vector[C] log_pi;
  log_pi = log(pi);
}

model{
  real ps[C];             
  real pj;
  real log_items[J];
  s ~ beta(1,1);
  g ~ beta(1,1);
  for (i in 1:N){
    for (c in 1:C){
      for (j in 1:J){
        pj = g[j] + (1-s[j]-g[j])*eta[c,j];
        log_items[j] = y[i,j] * log(pj) + (1 - y[i,j]) * log(1 - pj);
      }
      ps[c] = log_pi[c] + sum(log_items); 
    }
    target += log_sum_exp(ps);
  }
}

Code

library(GDINA)
score <- sim10GDINA$simdat
Q <- sim10GDINA$simQ
J <- ncol(score)
K=ncol(Q)
N=nrow(score)
all.patterns <- GDINA::attributepattern(K)
C <- nrow(all.patterns)

eta <- matrix(NA,C,J)
for(l in 1:C){
  for(j in 1:J){
    eta[l,j] <- 1 * (drop(all.patterns[l,]%*%Q[j,])>=sum(Q[j,]))
  }
}

dina_data <- list(N = N, J = J, K = K, C = C, y = score, eta = eta)


stan_inits <- list(list(g = runif(J, 0.1, 0.3), s = runif(J, 0.1, 0.3)), 
                   list(g = runif(J, 0.2, 0.4), s = runif(J, 0.2, 0.4)))

# Fit model to simulated data
library(rstan)
dina_stan_mcmc <- stan(file = "dina.stan", data = dina_data, 
                         chains = 2, iter = 5000,warmup =  2500, init = stan_inits)

Code

library(MCMCvis)
MCMCvis::MCMCsummary(dina_stan_mcmc, params = c('s','g'), 
                     Rhat = TRUE,round = 2)

##       mean   sd 2.5%  50% 97.5% Rhat n.eff
## s[1]  0.26 0.03 0.20 0.26  0.31    1  4517
## s[2]  0.28 0.03 0.22 0.28  0.33    1  5522
## s[3]  0.11 0.02 0.07 0.11  0.16    1  5857
## s[4]  0.20 0.03 0.14 0.20  0.25    1  4838
## s[5]  0.29 0.04 0.22 0.29  0.36    1  5319
## s[6]  0.07 0.02 0.04 0.07  0.11    1  6961
## s[7]  0.32 0.03 0.26 0.32  0.38    1  6037
## s[8]  0.30 0.04 0.23 0.30  0.37    1  3899
## s[9]  0.20 0.03 0.15 0.20  0.26    1  8282
## s[10] 0.16 0.03 0.10 0.16  0.22    1  6326
## g[1]  0.07 0.05 0.00 0.07  0.18    1  3524
## g[2]  0.09 0.03 0.04 0.09  0.16    1  4796
## g[3]  0.07 0.03 0.02 0.07  0.13    1  4278
## g[4]  0.22 0.02 0.18 0.22  0.27    1  6345
## g[5]  0.08 0.01 0.06 0.08  0.11    1  8461
## g[6]  0.59 0.03 0.53 0.59  0.64    1  5092
## g[7]  0.26 0.02 0.22 0.26  0.30    1  7905
## g[8]  0.15 0.02 0.11 0.15  0.19    1  6004
## g[9]  0.27 0.02 0.23 0.27  0.31    1  6979
## g[10] 0.34 0.02 0.30 0.33  0.37    1  6575

Code

MCMCvis::MCMCtrace(dina_stan_mcmc,params = c("g[1]","s[1]"),ISB = FALSE,
                   exact = TRUE,
                   ind = TRUE,
                   pdf = FALSE,
                   Rhat = TRUE,
                   n.eff = TRUE)