
We are very pleased to announce the release of a new package, dcmstan. The goal of dcmstan is to provide users with a friendly interface for creating Stan scripts necessary for estimating diagnostic classification models (DCMs; also called cognitive diagnostic models [CDMs]). dcmstan is primarily intended to serve as a backend for {measr}, which will interface with {rstan} or {cmdstanr} to actually estimate a DCM. However, dcmstan can also be used independently if you want finer control of the estimation process or would like to customize the Stan script used for estimation.
You can install dcmstan from CRAN with:
install.packages("dcmstan")This blog post will highlight the major features and describe how dcmstan fits into the larger universe of r-dcm packages.
library(dcmstan)DCM specification
For this example, we’ll create a specification for a DCM that we want to fit to ECPE data, which is available in the {dcmdata} package.
library(dcmdata)
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
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>We can create a DCM specification using dcm_specify(). We must provide our Q-matrix and, if present, the name of the Q-matrix column that contains the item identifiers. We then must choose a measurement and structural model to be used. By default, dcm_specify() will fit a loglinear cognitive diagnostic model (LCDM; Henson et al., 2009; Henson & Templin, 2019) with an unconstrained structural model. These are the measurement models used by Templin & Hoffman (2013) in their examination of the ECPE data, so we will echo those choices here. However, there are many other measurement and structural models we could choose from. For details on the specification options, see vignette("dcmstan", package = "dcmstan").
ecpe_spec <- dcm_specify(
qmatrix = ecpe_qmatrix,
identifier = "item_id",
measurement_model = lcdm(),
structural_model = unconstrained()
)
ecpe_spec
#> A loglinear cognitive diagnostic model (LCDM) measuring 3 attributes with 28
#> items.
#>
#> ℹ Attributes:
#> • "morphosyntactic" (13 items)
#> • "cohesive" (6 items)
#> • "lexical" (18 items)
#>
#> ℹ Attribute structure:
#> Unconstrained
#>
#> ℹ Prior distributions:
#> intercept ~ normal(0, 2)
#> maineffect ~ lognormal(0, 1)
#> interaction ~ normal(0, 2)
#> `Vc` ~ dirichlet(1, 1, 1)We can also specify prior distributions for the model. Reasonable priors are defined by default, but custom priors can be specified for specific parameters or an entire type of parameters (e.g., applies to all intercept parameters). For a list of possible parameters for given specification, see get_parameters().
get_parameters(ecpe_spec)
#> # A tibble: 75 × 4
#> item_id type attributes coefficient
#> <chr> <chr> <chr> <chr>
#> 1 E1 intercept <NA> l1_0
#> 2 E1 maineffect morphosyntactic l1_11
#> 3 E1 maineffect cohesive l1_12
#> 4 E1 interaction morphosyntactic__cohesive l1_212
#> 5 E2 intercept <NA> l2_0
#> 6 E2 maineffect cohesive l2_12
#> 7 E3 intercept <NA> l3_0
#> 8 E3 maineffect morphosyntactic l3_11
#> 9 E3 maineffect lexical l3_13
#> 10 E3 interaction morphosyntactic__lexical l3_213
#> # ℹ 65 more rowsStan code and data
Once we have our specification, we can create the Stan code for estimating the model with stan_code(). This code can either be passed directly to the model_code argument of rstan::stan(), or written to a .stan file, which can then be passed to the file argument of rstan::stan() or the stan_file argument of cmdstanr::cmdstan_model().
stan_code(ecpe_spec)#> data {
#> int<lower=1> I; // number of items
#> int<lower=1> R; // number of respondents
#> int<lower=1> N; // number of observations
#> int<lower=1> C; // number of classes
#> array[N] int<lower=1,upper=I> ii; // item for observation n
#> array[N] int<lower=1,upper=R> rr; // respondent for observation n
#> array[N] int<lower=0,upper=1> y; // score for observation n
#> array[R] int<lower=1,upper=N> start; // starting row for respondent R
#> array[R] int<lower=1,upper=I> num; // number items for respondent R
#> }
#> parameters {
#> simplex[C] Vc;
#>
#> ////////////////////////////////// item intercepts
#> real l1_0;
#> real l2_0;
#> real l3_0;
#> real l4_0;
#> real l5_0;
#> real l6_0;
#> real l7_0;
#> real l8_0;
#> real l9_0;
#> real l10_0;
#> real l11_0;
#> real l12_0;
#> real l13_0;
#> real l14_0;
#> real l15_0;
#> real l16_0;
#> real l17_0;
#> real l18_0;
#> real l19_0;
#> real l20_0;
#> real l21_0;
#> real l22_0;
#> real l23_0;
#> real l24_0;
#> real l25_0;
#> real l26_0;
#> real l27_0;
#> real l28_0;
#>
#> ////////////////////////////////// item main effects
#> real<lower=0> l1_11;
#> real<lower=0> l1_12;
#> real<lower=0> l2_12;
#> real<lower=0> l3_11;
#> real<lower=0> l3_13;
#> real<lower=0> l4_13;
#> real<lower=0> l5_13;
#> real<lower=0> l6_13;
#> real<lower=0> l7_11;
#> real<lower=0> l7_13;
#> real<lower=0> l8_12;
#> real<lower=0> l9_13;
#> real<lower=0> l10_11;
#> real<lower=0> l11_11;
#> real<lower=0> l11_13;
#> real<lower=0> l12_11;
#> real<lower=0> l12_13;
#> real<lower=0> l13_11;
#> real<lower=0> l14_11;
#> real<lower=0> l15_13;
#> real<lower=0> l16_11;
#> real<lower=0> l16_13;
#> real<lower=0> l17_12;
#> real<lower=0> l17_13;
#> real<lower=0> l18_13;
#> real<lower=0> l19_13;
#> real<lower=0> l20_11;
#> real<lower=0> l20_13;
#> real<lower=0> l21_11;
#> real<lower=0> l21_13;
#> real<lower=0> l22_13;
#> real<lower=0> l23_12;
#> real<lower=0> l24_12;
#> real<lower=0> l25_11;
#> real<lower=0> l26_13;
#> real<lower=0> l27_11;
#> real<lower=0> l28_13;
#>
#> ////////////////////////////////// item interactions
#> real<lower=-1 * min([l1_11,l1_12])> l1_212;
#> real<lower=-1 * min([l3_11,l3_13])> l3_213;
#> real<lower=-1 * min([l7_11,l7_13])> l7_213;
#> real<lower=-1 * min([l11_11,l11_13])> l11_213;
#> real<lower=-1 * min([l12_11,l12_13])> l12_213;
#> real<lower=-1 * min([l16_11,l16_13])> l16_213;
#> real<lower=-1 * min([l17_12,l17_13])> l17_223;
#> real<lower=-1 * min([l20_11,l20_13])> l20_213;
#> real<lower=-1 * min([l21_11,l21_13])> l21_213;
#> }
#> transformed parameters {
#> vector[C] log_Vc = log(Vc);
#> matrix[I,C] pi;
#>
#> ////////////////////////////////// probability of correct response
#> pi[1,1] = inv_logit(l1_0);
#> pi[1,2] = inv_logit(l1_0+l1_11);
#> pi[1,3] = inv_logit(l1_0+l1_12);
#> pi[1,4] = inv_logit(l1_0);
#> pi[1,5] = inv_logit(l1_0+l1_11+l1_12+l1_212);
#> pi[1,6] = inv_logit(l1_0+l1_11);
#> pi[1,7] = inv_logit(l1_0+l1_12);
#> pi[1,8] = inv_logit(l1_0+l1_11+l1_12+l1_212);
#> pi[2,1] = inv_logit(l2_0);
#> pi[2,2] = inv_logit(l2_0);
#> pi[2,3] = inv_logit(l2_0+l2_12);
#> pi[2,4] = inv_logit(l2_0);
#> pi[2,5] = inv_logit(l2_0+l2_12);
#> pi[2,6] = inv_logit(l2_0);
#> pi[2,7] = inv_logit(l2_0+l2_12);
#> pi[2,8] = inv_logit(l2_0+l2_12);
#> pi[3,1] = inv_logit(l3_0);
#> pi[3,2] = inv_logit(l3_0+l3_11);
#> pi[3,3] = inv_logit(l3_0);
#> pi[3,4] = inv_logit(l3_0+l3_13);
#> pi[3,5] = inv_logit(l3_0+l3_11);
#> pi[3,6] = inv_logit(l3_0+l3_11+l3_13+l3_213);
#> pi[3,7] = inv_logit(l3_0+l3_13);
#> pi[3,8] = inv_logit(l3_0+l3_11+l3_13+l3_213);
#> pi[4,1] = inv_logit(l4_0);
#> pi[4,2] = inv_logit(l4_0);
#> pi[4,3] = inv_logit(l4_0);
#> pi[4,4] = inv_logit(l4_0+l4_13);
#> pi[4,5] = inv_logit(l4_0);
#> pi[4,6] = inv_logit(l4_0+l4_13);
#> pi[4,7] = inv_logit(l4_0+l4_13);
#> pi[4,8] = inv_logit(l4_0+l4_13);
#> pi[5,1] = inv_logit(l5_0);
#> pi[5,2] = inv_logit(l5_0);
#> pi[5,3] = inv_logit(l5_0);
#> pi[5,4] = inv_logit(l5_0+l5_13);
#> pi[5,5] = inv_logit(l5_0);
#> pi[5,6] = inv_logit(l5_0+l5_13);
#> pi[5,7] = inv_logit(l5_0+l5_13);
#> pi[5,8] = inv_logit(l5_0+l5_13);
#> pi[6,1] = inv_logit(l6_0);
#> pi[6,2] = inv_logit(l6_0);
#> pi[6,3] = inv_logit(l6_0);
#> pi[6,4] = inv_logit(l6_0+l6_13);
#> pi[6,5] = inv_logit(l6_0);
#> pi[6,6] = inv_logit(l6_0+l6_13);
#> pi[6,7] = inv_logit(l6_0+l6_13);
#> pi[6,8] = inv_logit(l6_0+l6_13);
#> pi[7,1] = inv_logit(l7_0);
#> pi[7,2] = inv_logit(l7_0+l7_11);
#> pi[7,3] = inv_logit(l7_0);
#> pi[7,4] = inv_logit(l7_0+l7_13);
#> pi[7,5] = inv_logit(l7_0+l7_11);
#> pi[7,6] = inv_logit(l7_0+l7_11+l7_13+l7_213);
#> pi[7,7] = inv_logit(l7_0+l7_13);
#> pi[7,8] = inv_logit(l7_0+l7_11+l7_13+l7_213);
#> pi[8,1] = inv_logit(l8_0);
#> pi[8,2] = inv_logit(l8_0);
#> pi[8,3] = inv_logit(l8_0+l8_12);
#> pi[8,4] = inv_logit(l8_0);
#> pi[8,5] = inv_logit(l8_0+l8_12);
#> pi[8,6] = inv_logit(l8_0);
#> pi[8,7] = inv_logit(l8_0+l8_12);
#> pi[8,8] = inv_logit(l8_0+l8_12);
#> pi[9,1] = inv_logit(l9_0);
#> pi[9,2] = inv_logit(l9_0);
#> pi[9,3] = inv_logit(l9_0);
#> pi[9,4] = inv_logit(l9_0+l9_13);
#> pi[9,5] = inv_logit(l9_0);
#> pi[9,6] = inv_logit(l9_0+l9_13);
#> pi[9,7] = inv_logit(l9_0+l9_13);
#> pi[9,8] = inv_logit(l9_0+l9_13);
#> pi[10,1] = inv_logit(l10_0);
#> pi[10,2] = inv_logit(l10_0+l10_11);
#> pi[10,3] = inv_logit(l10_0);
#> pi[10,4] = inv_logit(l10_0);
#> pi[10,5] = inv_logit(l10_0+l10_11);
#> pi[10,6] = inv_logit(l10_0+l10_11);
#> pi[10,7] = inv_logit(l10_0);
#> pi[10,8] = inv_logit(l10_0+l10_11);
#> pi[11,1] = inv_logit(l11_0);
#> pi[11,2] = inv_logit(l11_0+l11_11);
#> pi[11,3] = inv_logit(l11_0);
#> pi[11,4] = inv_logit(l11_0+l11_13);
#> pi[11,5] = inv_logit(l11_0+l11_11);
#> pi[11,6] = inv_logit(l11_0+l11_11+l11_13+l11_213);
#> pi[11,7] = inv_logit(l11_0+l11_13);
#> pi[11,8] = inv_logit(l11_0+l11_11+l11_13+l11_213);
#> pi[12,1] = inv_logit(l12_0);
#> pi[12,2] = inv_logit(l12_0+l12_11);
#> pi[12,3] = inv_logit(l12_0);
#> pi[12,4] = inv_logit(l12_0+l12_13);
#> pi[12,5] = inv_logit(l12_0+l12_11);
#> pi[12,6] = inv_logit(l12_0+l12_11+l12_13+l12_213);
#> pi[12,7] = inv_logit(l12_0+l12_13);
#> pi[12,8] = inv_logit(l12_0+l12_11+l12_13+l12_213);
#> pi[13,1] = inv_logit(l13_0);
#> pi[13,2] = inv_logit(l13_0+l13_11);
#> pi[13,3] = inv_logit(l13_0);
#> pi[13,4] = inv_logit(l13_0);
#> pi[13,5] = inv_logit(l13_0+l13_11);
#> pi[13,6] = inv_logit(l13_0+l13_11);
#> pi[13,7] = inv_logit(l13_0);
#> pi[13,8] = inv_logit(l13_0+l13_11);
#> pi[14,1] = inv_logit(l14_0);
#> pi[14,2] = inv_logit(l14_0+l14_11);
#> pi[14,3] = inv_logit(l14_0);
#> pi[14,4] = inv_logit(l14_0);
#> pi[14,5] = inv_logit(l14_0+l14_11);
#> pi[14,6] = inv_logit(l14_0+l14_11);
#> pi[14,7] = inv_logit(l14_0);
#> pi[14,8] = inv_logit(l14_0+l14_11);
#> pi[15,1] = inv_logit(l15_0);
#> pi[15,2] = inv_logit(l15_0);
#> pi[15,3] = inv_logit(l15_0);
#> pi[15,4] = inv_logit(l15_0+l15_13);
#> pi[15,5] = inv_logit(l15_0);
#> pi[15,6] = inv_logit(l15_0+l15_13);
#> pi[15,7] = inv_logit(l15_0+l15_13);
#> pi[15,8] = inv_logit(l15_0+l15_13);
#> pi[16,1] = inv_logit(l16_0);
#> pi[16,2] = inv_logit(l16_0+l16_11);
#> pi[16,3] = inv_logit(l16_0);
#> pi[16,4] = inv_logit(l16_0+l16_13);
#> pi[16,5] = inv_logit(l16_0+l16_11);
#> pi[16,6] = inv_logit(l16_0+l16_11+l16_13+l16_213);
#> pi[16,7] = inv_logit(l16_0+l16_13);
#> pi[16,8] = inv_logit(l16_0+l16_11+l16_13+l16_213);
#> pi[17,1] = inv_logit(l17_0);
#> pi[17,2] = inv_logit(l17_0);
#> pi[17,3] = inv_logit(l17_0+l17_12);
#> pi[17,4] = inv_logit(l17_0+l17_13);
#> pi[17,5] = inv_logit(l17_0+l17_12);
#> pi[17,6] = inv_logit(l17_0+l17_13);
#> pi[17,7] = inv_logit(l17_0+l17_12+l17_13+l17_223);
#> pi[17,8] = inv_logit(l17_0+l17_12+l17_13+l17_223);
#> pi[18,1] = inv_logit(l18_0);
#> pi[18,2] = inv_logit(l18_0);
#> pi[18,3] = inv_logit(l18_0);
#> pi[18,4] = inv_logit(l18_0+l18_13);
#> pi[18,5] = inv_logit(l18_0);
#> pi[18,6] = inv_logit(l18_0+l18_13);
#> pi[18,7] = inv_logit(l18_0+l18_13);
#> pi[18,8] = inv_logit(l18_0+l18_13);
#> pi[19,1] = inv_logit(l19_0);
#> pi[19,2] = inv_logit(l19_0);
#> pi[19,3] = inv_logit(l19_0);
#> pi[19,4] = inv_logit(l19_0+l19_13);
#> pi[19,5] = inv_logit(l19_0);
#> pi[19,6] = inv_logit(l19_0+l19_13);
#> pi[19,7] = inv_logit(l19_0+l19_13);
#> pi[19,8] = inv_logit(l19_0+l19_13);
#> pi[20,1] = inv_logit(l20_0);
#> pi[20,2] = inv_logit(l20_0+l20_11);
#> pi[20,3] = inv_logit(l20_0);
#> pi[20,4] = inv_logit(l20_0+l20_13);
#> pi[20,5] = inv_logit(l20_0+l20_11);
#> pi[20,6] = inv_logit(l20_0+l20_11+l20_13+l20_213);
#> pi[20,7] = inv_logit(l20_0+l20_13);
#> pi[20,8] = inv_logit(l20_0+l20_11+l20_13+l20_213);
#> pi[21,1] = inv_logit(l21_0);
#> pi[21,2] = inv_logit(l21_0+l21_11);
#> pi[21,3] = inv_logit(l21_0);
#> pi[21,4] = inv_logit(l21_0+l21_13);
#> pi[21,5] = inv_logit(l21_0+l21_11);
#> pi[21,6] = inv_logit(l21_0+l21_11+l21_13+l21_213);
#> pi[21,7] = inv_logit(l21_0+l21_13);
#> pi[21,8] = inv_logit(l21_0+l21_11+l21_13+l21_213);
#> pi[22,1] = inv_logit(l22_0);
#> pi[22,2] = inv_logit(l22_0);
#> pi[22,3] = inv_logit(l22_0);
#> pi[22,4] = inv_logit(l22_0+l22_13);
#> pi[22,5] = inv_logit(l22_0);
#> pi[22,6] = inv_logit(l22_0+l22_13);
#> pi[22,7] = inv_logit(l22_0+l22_13);
#> pi[22,8] = inv_logit(l22_0+l22_13);
#> pi[23,1] = inv_logit(l23_0);
#> pi[23,2] = inv_logit(l23_0);
#> pi[23,3] = inv_logit(l23_0+l23_12);
#> pi[23,4] = inv_logit(l23_0);
#> pi[23,5] = inv_logit(l23_0+l23_12);
#> pi[23,6] = inv_logit(l23_0);
#> pi[23,7] = inv_logit(l23_0+l23_12);
#> pi[23,8] = inv_logit(l23_0+l23_12);
#> pi[24,1] = inv_logit(l24_0);
#> pi[24,2] = inv_logit(l24_0);
#> pi[24,3] = inv_logit(l24_0+l24_12);
#> pi[24,4] = inv_logit(l24_0);
#> pi[24,5] = inv_logit(l24_0+l24_12);
#> pi[24,6] = inv_logit(l24_0);
#> pi[24,7] = inv_logit(l24_0+l24_12);
#> pi[24,8] = inv_logit(l24_0+l24_12);
#> pi[25,1] = inv_logit(l25_0);
#> pi[25,2] = inv_logit(l25_0+l25_11);
#> pi[25,3] = inv_logit(l25_0);
#> pi[25,4] = inv_logit(l25_0);
#> pi[25,5] = inv_logit(l25_0+l25_11);
#> pi[25,6] = inv_logit(l25_0+l25_11);
#> pi[25,7] = inv_logit(l25_0);
#> pi[25,8] = inv_logit(l25_0+l25_11);
#> pi[26,1] = inv_logit(l26_0);
#> pi[26,2] = inv_logit(l26_0);
#> pi[26,3] = inv_logit(l26_0);
#> pi[26,4] = inv_logit(l26_0+l26_13);
#> pi[26,5] = inv_logit(l26_0);
#> pi[26,6] = inv_logit(l26_0+l26_13);
#> pi[26,7] = inv_logit(l26_0+l26_13);
#> pi[26,8] = inv_logit(l26_0+l26_13);
#> pi[27,1] = inv_logit(l27_0);
#> pi[27,2] = inv_logit(l27_0+l27_11);
#> pi[27,3] = inv_logit(l27_0);
#> pi[27,4] = inv_logit(l27_0);
#> pi[27,5] = inv_logit(l27_0+l27_11);
#> pi[27,6] = inv_logit(l27_0+l27_11);
#> pi[27,7] = inv_logit(l27_0);
#> pi[27,8] = inv_logit(l27_0+l27_11);
#> pi[28,1] = inv_logit(l28_0);
#> pi[28,2] = inv_logit(l28_0);
#> pi[28,3] = inv_logit(l28_0);
#> pi[28,4] = inv_logit(l28_0+l28_13);
#> pi[28,5] = inv_logit(l28_0);
#> pi[28,6] = inv_logit(l28_0+l28_13);
#> pi[28,7] = inv_logit(l28_0+l28_13);
#> pi[28,8] = inv_logit(l28_0+l28_13);
#> }
#> model {
#> ////////////////////////////////// priors
#> Vc ~ dirichlet(rep_vector(1, C));
#> l1_0 ~ normal(0, 2);
#> l1_11 ~ lognormal(0, 1);
#> l1_12 ~ lognormal(0, 1);
#> l1_212 ~ normal(0, 2);
#> l2_0 ~ normal(0, 2);
#> l2_12 ~ lognormal(0, 1);
#> l3_0 ~ normal(0, 2);
#> l3_11 ~ lognormal(0, 1);
#> l3_13 ~ lognormal(0, 1);
#> l3_213 ~ normal(0, 2);
#> l4_0 ~ normal(0, 2);
#> l4_13 ~ lognormal(0, 1);
#> l5_0 ~ normal(0, 2);
#> l5_13 ~ lognormal(0, 1);
#> l6_0 ~ normal(0, 2);
#> l6_13 ~ lognormal(0, 1);
#> l7_0 ~ normal(0, 2);
#> l7_11 ~ lognormal(0, 1);
#> l7_13 ~ lognormal(0, 1);
#> l7_213 ~ normal(0, 2);
#> l8_0 ~ normal(0, 2);
#> l8_12 ~ lognormal(0, 1);
#> l9_0 ~ normal(0, 2);
#> l9_13 ~ lognormal(0, 1);
#> l10_0 ~ normal(0, 2);
#> l10_11 ~ lognormal(0, 1);
#> l11_0 ~ normal(0, 2);
#> l11_11 ~ lognormal(0, 1);
#> l11_13 ~ lognormal(0, 1);
#> l11_213 ~ normal(0, 2);
#> l12_0 ~ normal(0, 2);
#> l12_11 ~ lognormal(0, 1);
#> l12_13 ~ lognormal(0, 1);
#> l12_213 ~ normal(0, 2);
#> l13_0 ~ normal(0, 2);
#> l13_11 ~ lognormal(0, 1);
#> l14_0 ~ normal(0, 2);
#> l14_11 ~ lognormal(0, 1);
#> l15_0 ~ normal(0, 2);
#> l15_13 ~ lognormal(0, 1);
#> l16_0 ~ normal(0, 2);
#> l16_11 ~ lognormal(0, 1);
#> l16_13 ~ lognormal(0, 1);
#> l16_213 ~ normal(0, 2);
#> l17_0 ~ normal(0, 2);
#> l17_12 ~ lognormal(0, 1);
#> l17_13 ~ lognormal(0, 1);
#> l17_223 ~ normal(0, 2);
#> l18_0 ~ normal(0, 2);
#> l18_13 ~ lognormal(0, 1);
#> l19_0 ~ normal(0, 2);
#> l19_13 ~ lognormal(0, 1);
#> l20_0 ~ normal(0, 2);
#> l20_11 ~ lognormal(0, 1);
#> l20_13 ~ lognormal(0, 1);
#> l20_213 ~ normal(0, 2);
#> l21_0 ~ normal(0, 2);
#> l21_11 ~ lognormal(0, 1);
#> l21_13 ~ lognormal(0, 1);
#> l21_213 ~ normal(0, 2);
#> l22_0 ~ normal(0, 2);
#> l22_13 ~ lognormal(0, 1);
#> l23_0 ~ normal(0, 2);
#> l23_12 ~ lognormal(0, 1);
#> l24_0 ~ normal(0, 2);
#> l24_12 ~ lognormal(0, 1);
#> l25_0 ~ normal(0, 2);
#> l25_11 ~ lognormal(0, 1);
#> l26_0 ~ normal(0, 2);
#> l26_13 ~ lognormal(0, 1);
#> l27_0 ~ normal(0, 2);
#> l27_11 ~ lognormal(0, 1);
#> l28_0 ~ normal(0, 2);
#> l28_13 ~ lognormal(0, 1);
#>
#> ////////////////////////////////// likelihood
#> for (r in 1:R) {
#> row_vector[C] ps;
#> for (c in 1:C) {
#> array[num[r]] real log_items;
#> for (m in 1:num[r]) {
#> int i = ii[start[r] + m - 1];
#> log_items[m] = y[start[r] + m - 1] * log(pi[i,c]) +
#> (1 - y[start[r] + m - 1]) * log(1 - pi[i,c]);
#> }
#> ps[c] = log_Vc[c] + sum(log_items);
#> }
#> target += log_sum_exp(ps);
#> }
#> }
Both rstan and cmdstanr also require a list of data objects that correspond to the data block of the Stan code. This can be created with stan_data(). Note that if you edit the generated Stan code to customize the estimation and add to the data block, you will need to also add corresponding objects to the data list before estimation.
dat <- stan_data(ecpe_spec, data = ecpe_data, identifier = "resp_id")
str(dat)
#> List of 9
#> $ I : int 28
#> $ R : int 2922
#> $ N : int 81816
#> $ C : int 8
#> $ ii : num [1:81816] 1 2 3 4 5 6 7 8 9 10 ...
#> $ rr : num [1:81816] 1 1 1 1 1 1 1 1 1 1 ...
#> $ y : int [1:81816] 1 1 1 0 1 1 1 1 1 1 ...
#> $ start: int [1:2922] 1 29 57 85 113 141 169 197 225 253 ...
#> $ num : int [1:2922] 28 28 28 28 28 28 28 28 28 28 ...dcmstan + measr
dcmstan is primarily intended to be a backend to {measr}. Many dcmstan functions are reexported by measr so that if you are using measr to estimate your model, you should not need to directly load or interact with dcmstan. For example, dcm_specify() is reexported by measr, so one can simply create a model specification and pass that directly to measr::dcm_estimate().
library(measr)
ecpe_spec <- dcm_specify(
qmatrix = ecpe_qmatrix,
identifier = "item_id",
measurement_model = lcdm(),
structural_model = unconstrained()
)
model <- dcm_estimate(
dcm_spec = ecpe_spec,
data = ecpe_data,
identifier = "resp_id"
)dcm_estimate() calls stan_code() and stan_data() internally to create the necessary Stan script and data list and then estimates the model using your chosen backend (i.e., rstan or cmdstanr). Thus, a direct interface with dcmstan should only be necessary if you want to modify the Stan code that is generated by default.
Acknowledgments
The research reported here was supported by the Institute of Education Sciences, U.S. Department of Education, through Grants R305D210045 and R305D240032 to the University of Kansas Center for Research, Inc., ATLAS. The opinions expressed are those of the authors and do not represent the views of the the Institute or the U.S. Department of Education.
Featured photo by Andreas Haslinger on Unsplash.
References
Citation
@online{thompson2025,
author = {Thompson, W. Jake},
title = {Dcmstan 0.1.0},
date = {2025-11-26},
url = {https://r-dcm.org/blog/2025-11-dcmstan-0.1.0/},
doi = {10.59350/1kz71-zdj91},
langid = {en}
}