Estimating diagnostic classification models
Source:vignettes/articles/model-estimation.Rmd
model-estimation.RmdIn 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 rowsSpecifying 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:
- 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).
- 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.
- 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.
- 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 rowsWe 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.
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.0167This 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 rowsAny 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.