Skip to contents

In this article, we will walk you through the steps for estimating diagnostic classification models (DCMs; also known as cognitive diagnostic models [CDMs]) using measr. We start with the data to analyze, understand how to specify a DCM to estimate, and learn how to customize the model estimation process (e.g., prior distributions).

To use the code in this article, you will need to install and load the measr package.

Examination for the Certificate of Proficiency in English Data

To demonstrate how measr works, we will examine data from the grammar section of the Examination for the Certificate of Proficiency in English (ECPE). This section of the ECPE consists of 28 item multiple-choice items that together measure three attributes: morphosyntactic rules, cohesive rules, and lexical rules. In total, data is available for 2,922 respondents. This data set has been widely studied in the DCM literature (e.g., Chen et al., 2018; Liu & Johnson, 2019; Templin & Bradshaw, 2014). Additionally, parameter estimates from other software (Mplus) are readily available from Templin & Hoffman (2013), so we can compare the parameter estimates from measr to estimates derived from other software.

The ECPE data is included within the measr package. For additional information on these data, see ?ecpe_data.

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_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

Specifying a DCM for Estimation

In measr, DCMs are specified and estimated using the measr_dcm() function. We’ll start by estimating a loglinear cognitive diagnostic model (LCDM). The LCDM is a general DCM that subsumes many other DCM subtypes (Henson et al., 2009).

First, we specify the data (data) and the Q-matrix (qmatrix) that should be used to estimate the model. Note that these are the only required arguments to the measr_dcm() function. If no other arguments are provided, sensible defaults (described below) will take care of the rest of the specification and estimation. Next, we can specify which columns, if any, in our data and qmatrix contain respondent identifiers and item identifiers in, respectively. If one or more of these variables are not present in the data, these arguments can be omitted, and measr will assign identifiers based on the row number (i.e., row 1 in qmatrix becomes item 1). We can then specify the type of DCM we want to estimate. The current options are "lcdm" (the default) "dina", and "dino" (see Estimating Other DCM Sub-Types below).

You also have the option to choose which estimation engine to use, via the "backend" argument. The default backend is backend = "rstan", which will use the rstan package to estimate the model. Alternatively, you can use the cmdstanr package to estimate the model by specifying backend = "cmdstanr". The cmdstanr package works by using a local installation of Stan to estimate the models, rather than the version that is pre-compiled in rstan. Once a backend has been chosen, we can supply additional arguments to those specific estimating functions. In the example below, I specify 500 warm-up iterations per chain, 500 post-warm-up iterations per chain, and 4 cores to run the chains in parallel. The full set up options available for rstan and cmdstanr can be found by looking at the help pages for rstan::sampling() and cmdstanr::`model-method-sample`, respectively.

Finally, because estimating these models can be time intensive, you can specify a file. If a file is specified, an R object of the fitted model will be automatically saved to the specified file. If the specified file already exists, then the fitted model will be read back into R, eliminating the need to re-estimate the model.

ecpe <- measr_dcm(data = ecpe_data, qmatrix = ecpe_qmatrix,
                  resp_id = "resp_id", item_id = "item_id", type = "lcdm",
                  method = "mcmc", backend = "cmdstanr",
                  chains = 4, iter_warmup = 1000, iter_sampling = 500,
                  parallel_chains = 4,
                  file = "fits/ecpe-lcdm")

Examining Parameter Estimates

Now that we’ve estimated a model, let’s compare our parameter estimates to those reported for the ECPE data by Templin & Hoffman (2013). We can start be looking at our estimates using measr_extract(). This function extracts different aspects of a model estimated with measr. Here, the estimate column reports estimated value for each parameter and a measure of the associated error (i.e., the standard deviation of the posterior distribution). For example, item E1 has four parameters, as it measures two attributes:

  1. An intercept, which represents the log-odds of providing a correct response for a respondent who is proficient in neither of the attributes this item measures (i.e., morphosyntactic rules and cohesive rules).
  2. A main effect for morphosyntactic rules, which represents the increase in the log-odds of providing a correct response for a respondent who is proficient in that attribute.
  3. A main effect for cohesive rules, which represents the increase in the log-odds of providing a correct response for a respondent who is proficient in that attribute.
  4. An interaction between morphosyntactic and cohesive rules, which is the change in the log-odds for a respondent who is proficient in both attributes.
item_parameters <- measr_extract(ecpe, what = "item_param")
item_parameters
#> # 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

We can compare these estimates to those that Templin & Hoffman (2013) reported when using different software to estimate the same model. In the figure below, most parameters fall on or very close to the dashed line, which represents perfect agreement.

Figure shows a strong correlation between item parameters, with only a few discrepancies off of the line of perfect agreement.

There are some parameters that deviate from the line of perfect agreement, but these are expected. For example, take item E7, which measures morphosyntactic and lexical rules. Both measr and Templin & Hoffman (2013) report values of approximately -0.09 for the intercept and 0.93 for the main effect of lexical rules. For the main effect of morphosyntactic rules, measr estimated a value of 1.59, compared to a value of 2.86 reported by Templin & Hoffman (2013), a difference of -1.26. Similarly, the interaction term estimated by measr is 0.36, compared to a value of -0.95 reported by Templin & Hoffman (2013), a difference of 1.31. This indicates that the log-odds of providing a correct response for an individual who has mastered both attributes is approximately the same, regardless of software. That is, for measr, we get a log-odds of -0.07 + 1.59 + 0.91 + 0.36 = 2.8, and from Templin & Hoffman (2013), we get a log-odds of -0.11 + 2.86 + 0.95 + -0.95 = 2.75. This is true for all of the differences in the figure. There is a change to the main effect for morphosyntactic rules and corresponding change to the interaction term that “cancels out” the difference.

So why is this happening? Let’s look at the proportion of respondents in each class. There are very few respondents who are proficient in morphosyntactic rules without also being proficient in the other attributes (classes 2, 5, and 6; less than 4% of all respondents). Therefore, there is less information for estimating the morphosyntactic main effects, which for items that measure multiple attributes, represent the increase in log-odds for proficiency in morphosyntactic rules conditional on not being proficient on the other attribute.

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

This means that the prior will have more influence on these parameters. Note in the above figure that the main effect estimates that are off the diagonal are less extreme when using measr. For example, the triangle at the top right is a main effect that is nearly 3 as reported by Templin & Hoffman (2013), but is just over 1.5 when using with measr. Thus, there is a regularizing effect, where the prior is pulling in extreme values, which is an intended outcome. Prior distributions for measr are discussed more in the following section.

Customizing the Model Estimation Process

Prior Distributions

In the code to estimate the LCDM above, we did not specify any prior distributions in the call to measr_dcm(). By default, measr uses the following prior distributions for the LCDM:

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))

As you can see, main effect parameters get a lognormal(0, 1) prior by default. Different prior distributions can be specified with the prior() function. For example, we can specify a normal(0, 10) prior for the main effects with:

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

By default, the prior is applied to all parameters in the class (i.e., all main effects). However, we can also apply a prior to a specific parameter. For example, here we specify a \(\chi^2\) distribution with 2 degrees of freedom as the default prior for main effects, and an exponential distribution with a rate of 2 for the main effect of attribute 1 (morphosyntactic rules) on item 7. To see all parameters (including class and coef) that can be specified, we can use get_parameters().

c(prior(chi_square(2), class = "maineffect"),
  prior(exponential(2), class = "maineffect", coef = "l7_11"))
#> # A tibble: 2 × 3
#>   class      coef  prior_def     
#>   <chr>      <chr> <chr>         
#> 1 maineffect NA    chi_square(2) 
#> 2 maineffect l7_11 exponential(2)

get_parameters(ecpe_qmatrix, item_id = "item_id", type = "lcdm")
#> # A tibble: 75 × 4
#>    item_id class       attributes                coef  
#>      <int> <chr>       <chr>                     <glue>
#>  1       1 intercept   NA                        l1_0  
#>  2       1 maineffect  morphosyntactic           l1_11 
#>  3       1 maineffect  cohesive                  l1_12 
#>  4       1 interaction morphosyntactic__cohesive l1_212
#>  5       2 intercept   NA                        l2_0  
#>  6       2 maineffect  cohesive                  l2_12 
#>  7       3 intercept   NA                        l3_0  
#>  8       3 maineffect  morphosyntactic           l3_11 
#>  9       3 maineffect  lexical                   l3_13 
#> 10       3 interaction morphosyntactic__lexical  l3_213
#> # ℹ 65 more rows

Any distribution that is supported by the Stan language can be used as a prior. A list of all distributions is available in the Stan documentation, and are linked to from the ?prior() help page.

Priors can be defined before estimating the function, or created at the same time the model is estimated. For example, both of the following are equivalent. Here we set the prior for main effects to be a truncated normal distribution with a lower bound of 0. This is done because the main effects in the LCDM are constrained to be positive to ensure monotonicity in the model. Additionally note that I’ve set method = "optim". This means that we will estimate the model using Stan’s optimizer, rather than using full Markov Chain Monte Carlo. Note that the prior still influences the model when using method = "optim", just as they do when using method = "mcmc" (the default).

new_prior <- prior(normal(0, 15), class = "maineffect", lb = 0)
new_lcdm <- measr_dcm(data = ecpe_data, qmatrix = ecpe_qmatrix,
                      resp_id = "resp_id", item_id = "item_id",
                      type = "lcdm", backend = "cmdstanr",
                      method = "optim",
                      prior = new_prior,
                      file = "fits/ecpe-optim-lcdm")

new_lcdm <- measr_dcm(data = ecpe_data, qmatrix = ecpe_qmatrix,
                      resp_id = "resp_id", item_id = "item_id",
                      type = "lcdm", backend = "cmdstanr",
                      method = "optim",
                      prior = c(prior(normal(0, 15), class = "maineffect",
                                      lb = 0)),
                      file = "fits/ecpe-optim-lcdm")

The priors used to estimate the model are saved in the returned model object, so we can always go back and see which priors were used if we are unsure. We can see that for the new_lcdm model, our specified normal prior was used for the main effects, but the default priors were still applied to the parameters for which we did not explicitly state a prior distribution.

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

Other DCM Sub-Types

Although a primary motivation for measr is to provide researchers with software that makes the LCDM readily accessible, a few other popular DCM subtypes are also supported. For example, we can estimate the deterministic inputs, noisy “and” gate (DINA, Junker & Sijtsma, 2001) or the deterministic inputs, noisy “or” gate (DINO, Templin & Henson, 2006) models by specifying a different type in the measr_dcm() function.

Future development work will continue to add functionality for more DCM subtypes. If there is a specific subtype you are interested in, or would like to see supported, please open an issue on the GitHub repository.

References

Chen, F., Liu, Y., Xin, T., & Cui, Y. (2018). Applying the \({M}_2\) statistic to evaluate the fit of diagnostic classification models in the presence of attribute hierarchies. Frontiers in Psychology, 9. https://doi.org/10.3389/fpsyg.2018.01875
Henson, R., Templin, J. L., & Willse, J. T. (2009). Defining a family of cognitive diagnosis models using log-linear models with latent variables. Psychometrika, 74(2), 191–210. https://doi.org/10.1007/s11336-008-9089-5
Junker, B. W., & Sijtsma, K. (2001). Cognitive assessment models with few assumptions, and connections with nonparametric item response theory. Applied Psychological Measurement, 25(3), 258–272. https://doi.org/10.1177/01466210122032064
Liu, X., & Johnson, M. S. (2019). Estimating CDMs using MCMC. In M. von Davier & Y.-S. Lee (Eds.), Handbook of Diagnostic Classification Models (pp. 629–646). Springer International Publishing. https://doi.org/10.1007/978-3-030-05584-4_31
Templin, J., & Bradshaw, L. (2014). Hierarchical diagnostic classification models: A family of models for estimating and testing attribute hierarchies. Psychometrika, 79(2), 317–339. https://doi.org/10.1007/s11336-013-9362-0
Templin, J., & Henson, R. A. (2006). Measurement of psychological disorders using cognitive diagnosis models. Psychological Methods, 11(3), 287–305. https://doi.org/10.1037/1082-989X.11.3.287
Templin, J., & Hoffman, L. (2013). Obtaining diagnostic classification model estimates using Mplus. Educational Measurement: Issues and Practice, 32(2), 37–50. https://doi.org/10.1111/emip.12010