Estimating diagnostic classification models

With Stan and measr

W. Jake Thompson, Ph.D.

Example data

Data for examples

  • Examination for the Certificate of Proficiency in English (ECPE; Templin & Hoffman, 2013)
    • 28 items measuring 3 total attributes
    • 2,922 respondents
  • 3 attributes
    • Morphosyntactic rules
    • Cohesive rules
    • Lexical rules

ECPE data

library(measr)

ecpe_data
# A tibble: 2,922 × 29
   resp_id    E1    E2    E3    E4    E5    E6    E7    E8    E9   E10   E11
     <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int>
 1       1     1     1     1     0     1     1     1     1     1     1     1
 2       2     1     1     1     1     1     1     1     1     1     1     1
 3       3     1     1     1     1     1     1     0     1     1     1     1
 4       4     1     1     1     1     1     1     1     1     1     1     1
 5       5     1     1     1     1     1     1     1     1     1     1     1
 6       6     1     1     1     1     1     1     1     1     1     1     1
 7       7     1     1     1     1     1     1     1     1     1     1     1
 8       8     0     1     1     1     1     1     0     1     1     1     0
 9       9     1     1     1     1     1     1     1     1     1     1     1
10      10     1     1     1     1     0     0     1     1     1     1     1
# ℹ 2,912 more rows
# ℹ 17 more variables: E12 <int>, E13 <int>, E14 <int>, E15 <int>, E16 <int>,
#   E17 <int>, E18 <int>, E19 <int>, E20 <int>, E21 <int>, E22 <int>,
#   E23 <int>, E24 <int>, E25 <int>, E26 <int>, E27 <int>, E28 <int>

ECPE Q-matrix

ecpe_qmatrix
# A tibble: 28 × 4
   item_id morphosyntactic cohesive lexical
   <chr>             <int>    <int>   <int>
 1 E1                    1        1       0
 2 E2                    0        1       0
 3 E3                    1        0       1
 4 E4                    0        0       1
 5 E5                    0        0       1
 6 E6                    0        0       1
 7 E7                    1        0       1
 8 E8                    0        1       0
 9 E9                    0        0       1
10 E10                   1        0       0
# ℹ 18 more rows

Data for exercises

Exercise 1

  1. Download the exercise files
usethis::use_course("measr-workshop/stancon2023-exercises")
  1. Open estimation.Rmd

  2. Run the setup chunk

  3. Explore mdm_data and mdm_qmatrix

    • How many items are in the data?
    • How many respondents are in the data?
    • How many attributes are measured?
03:00

mdm_data
# A tibble: 142 × 5
   respondent  mdm1  mdm2  mdm3  mdm4
   <chr>      <int> <int> <int> <int>
 1 7gray          1     1     1     1
 2 3dvje          1     1     1     1
 3 sl7fa          1     1     1     1
 4 bo4ps          1     1     1     1
 5 woxc7          1     1     1     1
 6 s9exl          1     1     1     1
 7 6cboa          1     1     1     1
 8 tckp1          1     1     1     1
 9 zfn2i          1     1     1     1
10 pre1b          1     1     1     1
# ℹ 132 more rows
  • 4 items
  • 142 respondents
mdm_qmatrix
# A tibble: 4 × 2
  item  multiplication
  <chr>          <int>
1 mdm1               1
2 mdm2               1
3 mdm3               1
4 mdm4               1
  • 4 items
  • 1 attribute
    • Multiplication skills

Estimation options

Existing software

Software Programs

  • Mplus, flexMIRT, mdltm
  • Limitations
    • Tedious to implement, expensive, limited licenses, etc.

R Packages

  • CDM, GDINA, mirt, blatant
  • Limitations
    • Limited to constrained DCMs, under-documented
    • Different packages have different functionality, and don’t talk to each other

DCMs with Stan

Stan logo.

Stan for DCMs: Structure

Like most Stan programs, there are four main blocks:

  1. data
  2. parameters
  3. transformed parameters
  4. model

Stan for DCMs: data

data {
  int<lower=1> I;           // # of items
  int<lower=1> R;           // # of respondents 
  int<lower=1> A;           // # of attributes
  int<lower=1> C;           // # of attribute profiles (latent classes) 
  matrix[R,I] y;            // response matrix
}

Stan for DCMs: parameters

parameters {
  simplex[C] Vc;            // probability of class membership
  
  // item parameters
  real l1_0;
  real<lower=0> l1_11;
  real<lower=0> l1_12;
  real<lower=-1 * min([l1_11,l1_12])> l1_212;
  
  real l2_0;
  real_l2_12;
  ...
}

Item parameters in Stan

parameters {
  simplex[C] Vc;
  
  // item parameters
  real l1_0;
  real<lower=0> l1_11;
  real<lower=0> l1_12;
  real<lower=-1 * min([l1_11,l1_12])> l1_212;
  
  real l2_0;
  real_l2_12;
  ...
}
# A tibble: 2 × 4
  item_id morphosyntactic cohesive lexical
  <chr>             <int>    <int>   <int>
1 E1                    1        1       0
2 E2                    0        1       0

\(\lambda_{1,0} + \lambda_{1,1(1)} + \lambda_{1,1(2)} + \lambda_{1,2(1,2)}\)

  • Item 1 measures attributes 1 and 2
    • Intercept (l1_0)
    • Two main effects (l1_11 and l1_12)
    • One two-way interaction (l1_212)

\(\lambda_{2,0} + \lambda_{2,1(2)}\)

  • Item 2 measures only attribute 2
    • Intercept (l2_0)
    • One main effect (l2_12)

Repeat for all remaining items.

Stan for DCMs: transformed parameters

transformed parameters {
  matrix[I,C] PImat;
  
  PImat[1,1] = inv_logit(l1_0);
  PImat[1,2] = inv_logit(l1_0 + l1_1);
  PImat[1,3] = inv_logit(l1_0 + l1_2);
  PImat[1,4] = inv_logit(l1_0);
  PImat[1,5] = inv_logit(l1_0 + l1_11 + l1_12 + l1_212);
  PImat[1,6] = inv_logit(l1_0 + l1_1);
  PImat[1,7] = inv_logit(l1_0 + l1_2);
  PImat[1,8] = inv_logit(l1_0 + l1_11 + l1_12 + l1_212);
  ...
}

DCM classes

  • When using binary attributes, there are 2A possible profiles
create_profiles(3)
# A tibble: 8 × 3
   att1  att2  att3
  <int> <int> <int>
1     0     0     0
2     1     0     0
3     0     1     0
4     0     0     1
5     1     1     0
6     1     0     1
7     0     1     1
8     1     1     1

Mapping classes to PImat

transformed parameters {
  matrix[I,C] PImat;
  
  PImat[1,1] = inv_logit(l1_0);
  PImat[1,2] = inv_logit(l1_0 + l1_1);
  PImat[1,3] = inv_logit(l1_0 + l1_2);
  PImat[1,4] = inv_logit(l1_0);
  PImat[1,5] = inv_logit(l1_0 + l1_11 + l1_12 + l1_212);
  PImat[1,6] = inv_logit(l1_0 + l1_1);
  PImat[1,7] = inv_logit(l1_0 + l1_2);
  PImat[1,8] = inv_logit(l1_0 + l1_11 + l1_12 + l1_212);
  ...
}
create_profiles(3)
# A tibble: 8 × 3
   att1  att2  att3
  <int> <int> <int>
1     0     0     0
2     1     0     0
3     0     1     0
4     0     0     1
5     1     1     0
6     1     0     1
7     0     1     1
8     1     1     1

Stan for DCMs: model

\(P(X_r=x_r) = \sum_{c=1}^C\nu_c \prod_{i=1}^I\pi_{ic}^{x_{ir}}(1-\pi_{ic})^{1 - x_{ir}}\)

model {
  for (r in 1:R) {
    real ps[C];
    for (c in 1:C) {
      real log_items[I];
      for (i in 1:I) {
        log_items[i] = y[r,i] * log(PImat[i,c]) + (1 - y[r,i]) * log(1 - pi[i,c]);
      }
      ps[c] = log(Vc[c]) + sum(log_items);
    }
    target += log_sum_exp(ps);
  }
}

Exercise 2

Complete the parameters and transformed parameters blocks for the MDM data.

05:00

parameters {
  simplex[C] Vc;
  
  // item parameters
  real l1_0;
  real<lower=0> l1_11;
  real l2_0;
  real<lower=0> l2_11;
  real l3_0;
  real<lower=0> l3_11;
  real l4_0;
  real<lower=0> l4_11;
}
transformed parameters {
  matrix[I,C] PImat;
  
  PImat[1,1] = inv_logit(l1_0);
  PImat[1,2] = inv_logit(l1_0 + l1_11);
  PImat[2,1] = inv_logit(l2_0);
  PImat[2,2] = inv_logit(l2_0 + l2_11);
  PImat[3,1] = inv_logit(l3_0);
  PImat[3,2] = inv_logit(l3_0 + l3_11);
  PImat[4,1] = inv_logit(l4_0);
  PImat[4,2] = inv_logit(l4_0 + l4_11);
}

Limitations of Stan for DCMs

  • Very tedious—prone to typos
  • Complexity increases with the number of attributes and item structure
  • parameters and transformed parameters blocks have to be customized to each particular Q-matrix

Benefits of Stan for DCMs

  • Powerful and flexible

  • Access to other packages in the Stan ecosystem

  • Just need to automate the creation of the Stan scripts…

Hex logo for the measr R package.

What is measr?

  • R package that automates the creation of Stan scripts for DCMs
  • Wraps rstan or cmdstanr to estimate the models
  • Provides additional functions to automate the evaluation of DCMs
    • Model fit
    • Classification accuracy and consistency

measr_dcm()

Estimate a DCM with Stan

ecpe <- measr_dcm(
  data = ecpe_data, qmatrix = ecpe_qmatrix,
  resp_id = "resp_id", item_id = "item_id",
  type = "lcdm",
  method = "mcmc", backend = "cmdstanr",
  iter_warmup = 1000, iter_sampling = 500,
  chains = 4, parallel_chains = 4,
  file = "fits/ecpe-lcdm"
)
1
Specify your data, Q-matrix, and ID columns
2
Choose the DCM to estimate (e.g., LCDM, DINA, etc.)
3
Choose the estimation engine
4
Pass additional arguments to rstan or cmdstanr
5
Save the model to save time in the future

measr_dcm() options

  • type: Declare the type of DCM to estimate. Currently support LCDM, DINA, DINO, and C-RUM

  • method: How to estimate the model. To sample, use “mcmc”. To use Stan’s optimizer, use “optim”

  • backend: Which engine to use, either “rstan” or “cmdstanr”

  • ...: Additional arguments that are passed to, depending on the method and backend:

Exercise 3

Estimate and LCDM on the MDM data. Your model should have:

  • 2 chains
  • 1000 warmup and 500 sampling iterations
  • Use whichever backend you prefer
05:00

mdm_lcdm <- measr_dcm(
  data = mdm_data, qmatrix = mdm_qmatrix,
  resp_id = "respondent", item_id = "item",
  type = "lcdm",
  method = "mcmc", backend = "rstan",
  warmup = 1000, iter = 1500, chains = 2,
  file = "fits/mdm-lcdm"
)

Extracting item parameters

measr_extract() can be used to pull out components of an estimated model.

measr_extract(ecpe, "item_param")
# A tibble: 74 × 5
   item_id class       attributes                coef         estimate
   <fct>   <chr>       <chr>                     <glue>     <rvar[1d]>
 1 E1      intercept   <NA>                      l1_0     0.82 ± 0.074
 2 E1      maineffect  morphosyntactic           l1_11    0.59 ± 0.356
 3 E1      maineffect  cohesive                  l1_12    0.64 ± 0.217
 4 E1      interaction morphosyntactic__cohesive l1_212   0.54 ± 0.481
 5 E2      intercept   <NA>                      l2_0     1.04 ± 0.077
 6 E2      maineffect  cohesive                  l2_12    1.24 ± 0.147
 7 E3      intercept   <NA>                      l3_0    -0.35 ± 0.073
 8 E3      maineffect  morphosyntactic           l3_11    0.74 ± 0.379
 9 E3      maineffect  lexical                   l3_13    0.36 ± 0.115
10 E3      interaction morphosyntactic__lexical  l3_213   0.53 ± 0.395
# ℹ 64 more rows

Extracting structural parameters

measr_extract(ecpe, "strc_param")
# A tibble: 8 × 2
  class           estimate
  <chr>         <rvar[1d]>
1 [0,0,0]  0.2981 ± 0.0165
2 [1,0,0]  0.0118 ± 0.0063
3 [0,1,0]  0.0164 ± 0.0109
4 [0,0,1]  0.1279 ± 0.0198
5 [1,1,0]  0.0093 ± 0.0057
6 [1,0,1]  0.0184 ± 0.0100
7 [0,1,1]  0.1730 ± 0.0203
8 [1,1,1]  0.3451 ± 0.0167

Extracting respondent probabilities

ecpe <- add_respondent_estimates(ecpe)
measr_extract(ecpe, "attribute_prob")
# A tibble: 2,922 × 4
   resp_id morphosyntactic cohesive lexical
   <fct>             <dbl>    <dbl>   <dbl>
 1 1               0.997      0.955   1.00 
 2 2               0.995      0.899   1.00 
 3 3               0.983      0.988   1.00 
 4 4               0.998      0.990   1.00 
 5 5               0.988      0.980   0.955
 6 6               0.992      0.989   1.00 
 7 7               0.992      0.989   1.00 
 8 8               0.00436    0.444   0.961
 9 9               0.945      0.984   0.999
10 10              0.545      0.123   0.115
# ℹ 2,912 more rows

Exercise 4

What proportion of respondents have mastered multiplication in the MDM data?

What is the probability that respondent zfn2i has mastered multiplication?

02:00

measr_extract(mdm_lcdm, "strc_param")
# A tibble: 2 × 2
  class     estimate
  <chr>   <rvar[1d]>
1 [0]    0.5 ± 0.074
2 [1]    0.5 ± 0.074
mdm_lcdm <- add_respondent_estimates(mdm_lcdm)
measr_extract(mdm_lcdm, "attribute_prob")
# A tibble: 142 × 2
   respondent multiplication
   <fct>               <dbl>
 1 7gray               0.997
 2 3dvje               0.997
 3 sl7fa               0.997
 4 bo4ps               0.997
 5 woxc7               0.997
 6 s9exl               0.997
 7 6cboa               0.997
 8 tckp1               0.997
 9 zfn2i               0.997
10 pre1b               0.997
# ℹ 132 more rows

Working with Stan objects

Extract Stan code with:

ecpe$stancode

Extract the Stan object:

ecpe$model

Priors

Default priors

default_dcm_priors(type = "lcdm")
# A tibble: 4 × 3
  class       coef  prior_def                  
  <chr>       <chr> <chr>                      
1 intercept   <NA>  normal(0, 2)               
2 maineffect  <NA>  lognormal(0, 1)            
3 interaction <NA>  normal(0, 2)               
4 structural  Vc    dirichlet(rep_vector(1, C))

Weakly informative priors

The distribution of expected probabilities of providing a correct response for each class, based on the default priors.

Creating custom priors

prior(normal(0, 10), class = "maineffect")
# A tibble: 1 × 3
  class      coef  prior_def    
  <chr>      <chr> <chr>        
1 maineffect <NA>  normal(0, 10)

prior(normal(0, 1), class = "intercept", coef = "l3_0")
# A tibble: 1 × 3
  class     coef  prior_def   
  <chr>     <chr> <chr>       
1 intercept l3_0  normal(0, 1)

Specifying custom priors

my_prior <- prior(normal(-1, 1), class = "intercept")
new_ecpe <- measr_dcm(
  data = ecpe_data, qmatrix = ecpe_qmatrix,
  resp_id = "resp_id", item_id = "item_id", 
  type = "lcdm",  
  prior = my_prior,
  method = "optim", backend = "cmdstanr",
  file = "fits/ecpe-new-prior"
)
new_ecpe <- measr_dcm(
  data = ecpe_data, qmatrix = ecpe_qmatrix,
  resp_id = "resp_id", item_id = "item_id", 
  type = "lcdm",  
  prior = prior(normal(-1, 1), class = "intercept"),
  method = "optim", backend = "cmdstanr",
  file = "fits/ecpe-new-prior"
)

Extract priors

measr_extract(ecpe, "prior")
# A tibble: 4 × 3
  class       coef  prior_def                  
  <chr>       <chr> <chr>                      
1 intercept   <NA>  normal(0, 2)               
2 maineffect  <NA>  lognormal(0, 1)            
3 interaction <NA>  normal(0, 2)               
4 structural  Vc    dirichlet(rep_vector(1, C))

measr_extract(new_ecpe, "prior")
# A tibble: 4 × 3
  class       coef  prior_def                  
  <chr>       <chr> <chr>                      
1 intercept   <NA>  normal(-1, 1)              
2 maineffect  <NA>  lognormal(0, 1)            
3 interaction <NA>  normal(0, 2)               
4 structural  Vc    dirichlet(rep_vector(1, C))

Exercise 5

Estimate a new LCDM model for the MDM data with the following priors:

  • Intercepts: normal(-1, 2)
    • Except item 3, which should use a normal(0, 1) prior
  • Main effects: lognormal(0, 1)

Extract the prior to see that the specifications were applied

05:00

new_mdm <- measr_dcm(
  data = mdm_data, qmatrix = mdm_qmatrix,
  resp_id = "respondent", item_id = "item",
  type = "lcdm",
  prior = c(prior(normal(-1, 2), class = "intercept"),
            prior(normal(0, 1), class = "intercept", coef = "l3_0")),
  method = "mcmc", backend = "rstan",
  warmup = 1000, iter = 1500, chains = 2,
  file = "fits/mdm-new-prior"
)

measr_extract(new_mdm, "prior")
# A tibble: 4 × 3
  class      coef  prior_def                  
  <chr>      <chr> <chr>                      
1 intercept  <NA>  normal(-1, 2)              
2 intercept  l3_0  normal(0, 1)               
3 maineffect <NA>  lognormal(0, 1)            
4 structural Vc    dirichlet(rep_vector(1, C))

Estimating diagnostic classification models

With Stan and measr