Skip to content

Commit 8d839c8

Browse files
Merge branch 'main' into doc-cleanup
2 parents 6961f2a + 4304305 commit 8d839c8

File tree

15 files changed

+370
-17
lines changed

15 files changed

+370
-17
lines changed

.github/workflows/lint.yaml

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -22,7 +22,7 @@ jobs:
2222

2323
- uses: r-lib/actions/setup-r-dependencies@v2
2424
with:
25-
extra-packages: any::lintr, local::.
25+
extra-packages: r-lib/lintr, local::.
2626
needs: lint
2727

2828
- name: Lint

.lintr

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,6 @@
11
linters: linters_with_defaults(
2-
indentation_linter = NULL
2+
indentation_linter = NULL,
3+
return_linter = NULL
34
)
45
exclusions: list(
56
"R/stanmodels.R",

NAMESPACE

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -24,6 +24,7 @@ export(as_draws)
2424
export(as_label)
2525
export(as_measrfit)
2626
export(as_name)
27+
export(cdi)
2728
export(create_profiles)
2829
export(default_dcm_priors)
2930
export(enquo)

R/discrimination.R

Lines changed: 174 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,174 @@
1+
#' Item, attribute, and test-level discrimination indices
2+
#'
3+
#' The cognitive diagnostic index (CDI) is a measure of how well an assessment
4+
#' is able to distinguish between attribute profiles. The index was originally
5+
#' proposed by Henson & Douglas (2005) for item- and test-level discrimination,
6+
#' and then expanded by Henson et al. (2008) to include attribute-level
7+
#' discrimination indices.
8+
#'
9+
#' @param model The estimated model to be evaluated.
10+
#' @param weight_prevalence Logical indicating whether the discrimination
11+
#' indices should be weighted by the prevalence of the attribute profiles. See
12+
#' details for additional information.
13+
#'
14+
#' @details
15+
#' Henson et al. (2008) described two attribute-level discrimination indices,
16+
#' \eqn{\mathbf{d}_{(A)\mathbf{\cdot}}} (Equation 8) and
17+
#' \eqn{\mathbf{d}_{(B)\mathbf{\cdot}}} (Equation 13), which are similar in that
18+
#' both are the sum of item-level discrimination indices.
19+
#' In both cases, item-level discrimination indices are calculated as the
20+
#' average of Kullback-Leibler information for all pairs of attributes profiles
21+
#' for the item.
22+
#' The item-level indices are then summed to achieve the test-level
23+
#' discrimination index for each attribute, or the test overall.
24+
#' However, whereas \eqn{\mathbf{d}_{(A)\mathbf{\cdot}}} is an unweighted
25+
#' average of the Kullback-Leibler information,
26+
#' \eqn{\mathbf{d}_{(B)\mathbf{\cdot}}} is a weighted average, where the weight
27+
#' is defined by the prevalence of each profile (i.e.,
28+
#' [`measr_extract(model, what = "strc_param")`][measr_extract()]).
29+
#'
30+
#' @return A list with two elements:
31+
#' * `item_discrimination`: A [tibble][tibble::tibble-package] with one row
32+
#' per item containing the CDI for the item and any relevant attributes.
33+
#' * `test_discrimination`: A [tibble][tibble::tibble-package] with one row
34+
#' containing the total CDI for the assessment and for each attribute.
35+
#' @export
36+
#'
37+
#' @references Henson, R., & Douglas, J. (2005). Test construction for cognitive
38+
#' diagnosis. *Applied Psychological Measurement, 29*(4), 262-277.
39+
#' \doi{10.1177/0146621604272623}
40+
#' @references Henson, R., Roussos, L., Douglas, J., & Xuming, H. (2008).
41+
#' Cognitive diagnostic attribute-level discrimination indices.
42+
#' *Applied Psychological Measurement, 32*(4), 275-288.
43+
#' \doi{10.1177/0146621607302478}
44+
#' @examplesIf measr_examples()
45+
#' rstn_ecpe_lcdm <- measr_dcm(
46+
#' data = ecpe_data, missing = NA, qmatrix = ecpe_qmatrix,
47+
#' resp_id = "resp_id", item_id = "item_id", type = "lcdm",
48+
#' method = "optim", seed = 63277, backend = "rstan"
49+
#' )
50+
#'
51+
#' cdi(rstn_ecpe_lcdm)
52+
cdi <- function(model, weight_prevalence = TRUE) {
53+
model <- check_model(model, required_class = "measrfit", name = "model")
54+
weight_prevalence <- check_logical(weight_prevalence,
55+
name = "weight_prevalence")
56+
57+
stan_draws <- switch(model$method,
58+
"mcmc" = get_mcmc_draws(model),
59+
"optim" = get_optim_draws(model))
60+
61+
pi_matrix <- stan_draws %>%
62+
posterior::subset_draws(variable = "pi") %>%
63+
posterior::as_draws_df() %>%
64+
tibble::as_tibble() %>%
65+
tidyr::pivot_longer(cols = -c(".chain", ".iteration", ".draw")) %>%
66+
dplyr::summarize(value = mean(.data$value), .by = "name") %>%
67+
tidyr::separate_wider_regex(
68+
cols = "name",
69+
patterns = c("pi\\[", item = "[0-9]*", ",", class = "[0-9]*", "\\]")
70+
) %>%
71+
dplyr::mutate(item = as.integer(.data$item),
72+
class = as.integer(.data$class))
73+
74+
hamming <- profile_hamming(
75+
dplyr::select(measr_extract(model, "classes"), -"class")
76+
)
77+
att_names <- hamming %>%
78+
dplyr::select(-c("profile_1", "profile_2", "hamming")) %>%
79+
colnames()
80+
81+
item_discrim <- tidyr::crossing(item = unique(pi_matrix$item),
82+
profile_1 = unique(pi_matrix$class),
83+
profile_2 = unique(pi_matrix$class)) %>%
84+
dplyr::left_join(pi_matrix, by = c("item", "profile_1" = "class"),
85+
relationship = "many-to-one") %>%
86+
dplyr::rename("prob_1" = "value") %>%
87+
dplyr::left_join(pi_matrix, by = c("item", "profile_2" = "class"),
88+
relationship = "many-to-one") %>%
89+
dplyr::rename("prob_2" = "value") %>%
90+
dplyr::mutate(kli = (.data$prob_1 * log(.data$prob_1 / .data$prob_2)) +
91+
((1 - .data$prob_1) *
92+
log((1 - .data$prob_1) / (1 - .data$prob_2)))) %>%
93+
dplyr::left_join(hamming, by = c("profile_1", "profile_2"),
94+
relationship = "many-to-one") %>%
95+
dplyr::mutate(dplyr::across(dplyr::where(is.logical),
96+
\(x) {
97+
dplyr::case_when(
98+
x & .data$hamming == 1L ~ TRUE,
99+
.default = NA
100+
)
101+
}),
102+
dplyr::across(dplyr::where(is.logical),
103+
\(x) as.integer(x) * .data$kli)) %>%
104+
dplyr::filter(.data$hamming > 0) %>%
105+
dplyr::mutate(weight = 1 / .data$hamming)
106+
107+
if (weight_prevalence) {
108+
vc <- stan_draws %>%
109+
posterior::subset_draws(variable = "log_Vc") %>%
110+
posterior::as_draws_df() %>%
111+
tibble::as_tibble() %>%
112+
tidyr::pivot_longer(cols = -c(".chain", ".iteration", ".draw")) %>%
113+
dplyr::summarize(value = mean(.data$value), .by = "name") %>%
114+
dplyr::mutate(value = exp(.data$value)) %>%
115+
tidyr::separate_wider_regex(
116+
cols = "name",
117+
patterns = c("log_Vc\\[", class = "[0-9]*", "\\]")
118+
) %>%
119+
dplyr::mutate(class = as.integer(.data$class))
120+
121+
item_discrim <- item_discrim %>%
122+
dplyr::left_join(vc, by = c("profile_1" = "class")) %>%
123+
dplyr::mutate(weight = .data$weight * .data$value) %>%
124+
dplyr::select(-"value")
125+
}
126+
127+
item_discrim <- item_discrim %>%
128+
dplyr::summarize(
129+
overall = stats::weighted.mean(.data$kli, w = .data$weight),
130+
dplyr::across(
131+
dplyr::all_of(att_names),
132+
\(x) stats::weighted.mean(x, w = .data$weight, na.rm = TRUE)
133+
),
134+
.by = "item"
135+
)
136+
137+
test_discrim <- item_discrim %>%
138+
dplyr::summarize(dplyr::across(-"item", sum))
139+
140+
return(
141+
list(item_discrimination = item_discrim,
142+
test_discrimination = test_discrim)
143+
)
144+
}
145+
146+
profile_hamming <- function(profiles) {
147+
profile_combos <- tidyr::crossing(profile_1 = seq_len(nrow(profiles)),
148+
profile_2 = seq_len(nrow(profiles)))
149+
150+
151+
hamming <- mapply(hamming_distance, profile_combos$profile_1,
152+
profile_combos$profile_2,
153+
MoreArgs = list(profiles = profiles),
154+
SIMPLIFY = FALSE, USE.NAMES = FALSE) %>%
155+
dplyr::bind_rows()
156+
157+
dplyr::bind_cols(profile_combos, hamming)
158+
}
159+
160+
hamming_distance <- function(prof1, prof2, profiles) {
161+
pattern1 <- profiles[prof1, ]
162+
pattern2 <- profiles[prof2, ]
163+
164+
pattern1 %>%
165+
tidyr::pivot_longer(cols = dplyr::everything(),
166+
names_to = "att", values_to = "patt1") %>%
167+
dplyr::left_join(tidyr::pivot_longer(pattern2, cols = dplyr::everything(),
168+
names_to = "att", values_to = "patt2"),
169+
by = "att", relationship = "one-to-one") %>%
170+
dplyr::mutate(mismatch = .data$patt1 != .data$patt2,
171+
hamming = sum(.data$mismatch)) %>%
172+
dplyr::select("att", "mismatch", "hamming") %>%
173+
tidyr::pivot_wider(names_from = "att", values_from = "mismatch")
174+
}

R/model-evaluation.R

Lines changed: 4 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -50,7 +50,8 @@
5050
#' `$reliability` element of the fitted model. Pattern level reliability is
5151
#' described by Cui et al. (2012). Classification reliability and posterior
5252
#' probability reliability are described by Johnson & Sinharay (2018, 2020),
53-
#' respectively. This function wraps [reliability()].
53+
#' respectively. This function wraps [reliability()]. Arguments supplied to
54+
#' `...` are passed to [reliability()].
5455
#'
5556
#' @return A modified [measrfit] object with the corresponding slot populated
5657
#' with the specified information.
@@ -160,7 +161,7 @@ add_criterion <- function(x, criterion = c("loo", "waic"), overwrite = FALSE,
160161

161162
#' @export
162163
#' @rdname model_evaluation
163-
add_reliability <- function(x, overwrite = FALSE, save = TRUE) {
164+
add_reliability <- function(x, overwrite = FALSE, save = TRUE, ...) {
164165
model <- check_model(x, required_class = "measrfit", name = "x")
165166
overwrite <- check_logical(overwrite, name = "overwrite")
166167
save <- check_logical(save, name = "force_save")
@@ -169,7 +170,7 @@ add_reliability <- function(x, overwrite = FALSE, save = TRUE) {
169170
run_reli <- length(model$reliability) == 0 || overwrite
170171

171172
if (run_reli) {
172-
model$reliability <- reliability(model)
173+
model$reliability <- reliability(model, force = TRUE, ...)
173174
}
174175

175176
# re-save model object (if applicable)

R/reliability.R

Lines changed: 30 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -28,7 +28,11 @@ reliability <- function(model, ...) {
2828
#'
2929
#' @param threshold For `map_reliability`, the threshold applied to the
3030
#' attribute-level probabilities for determining the binary attribute
31-
#' classifications.
31+
#' classifications. Should be a numeric vector of length 1 (the same threshold
32+
#' is applied to all attributes), or length equal to the number of attributes.
33+
#' If a named vector is supplied, names should match the attribute names in the
34+
#' Q-matrix used to estimate the model. If unnamed, thresholds should be in the
35+
#' order the attributes were defined in the Q-matrix.
3236
#'
3337
#' @details The pattern-level reliability (`pattern_reliability`) statistics are
3438
#' described in Cui et al. (2012). Attribute-level classification reliability
@@ -71,16 +75,38 @@ reliability <- function(model, ...) {
7175
#'
7276
#' reliability(rstn_mdm_lcdm)
7377
reliability.measrdcm <- function(model, ..., threshold = 0.5, force = FALSE) {
78+
threshold <- check_double(threshold, lb = 0, ub = 1, inclusive = FALSE,
79+
name = "threshold")
80+
force <- check_logical(force, name = "force")
81+
82+
att_names <- colnames(dplyr::select(model$data$qmatrix, -"item_id"))
83+
if (length(threshold) == 1) {
84+
threshold <- rep(threshold, times = length(att_names)) %>%
85+
rlang::set_names(att_names)
86+
} else if (length(threshold) == length(att_names)) {
87+
if (is.null(names(threshold))) {
88+
threshold <- rlang::set_names(threshold, att_names)
89+
} else if (!all(names(threshold) %in% att_names)) {
90+
bad_names <- setdiff(names(threshold), att_names)
91+
rlang::abort(
92+
message = glue::glue("Unknown attribute names provided: ",
93+
"{paste(bad_names, collapse = ', ')}")
94+
)
95+
}
96+
} else {
97+
rlang::abort(
98+
message = glue::glue("`threshold` must be of length 1 or length ",
99+
"{length(att_names)} (the number of attributes).")
100+
)
101+
}
102+
74103
if ((!is.null(model$reliability) && length(model$reliability) > 0) &&
75104
!force) {
76105
return(model$reliability)
77106
}
78107

79108
# coerce model into a list of values required for reliability
80109
obj <- reli_list(model, threshold = threshold)
81-
att_names <- model$data$qmatrix %>%
82-
dplyr::select(-"item_id") %>%
83-
colnames()
84110

85111
tbl <- obj$acc
86112
p <- obj$prev

R/utils-reliability.R

Lines changed: 13 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -159,9 +159,19 @@ reli_list <- function(model, threshold) {
159159

160160
# map estimates
161161
binary_att <- attr_probs %>%
162-
dplyr::mutate(dplyr::across(dplyr::everything(),
163-
~dplyr::case_when(.x >= threshold ~ 1L,
164-
TRUE ~ 0L))) %>%
162+
tibble::rowid_to_column(var = "resp_id") %>%
163+
tidyr::pivot_longer(cols = -"resp_id",
164+
names_to = "attribute", values_to = "probability") %>%
165+
dplyr::left_join(tibble::enframe(threshold, name = "attribute",
166+
value = "threshold"),
167+
by = "attribute",
168+
relationship = "many-to-one") %>%
169+
dplyr::mutate(class = dplyr::case_when(.data$probability >=
170+
.data$threshold ~ 1L,
171+
.default = 0L)) %>%
172+
dplyr::select("resp_id", "attribute", "class") %>%
173+
tidyr::pivot_wider(names_from = "attribute", values_from = "class") %>%
174+
dplyr::select(dplyr::all_of(names(threshold))) %>%
165175
as.matrix() %>%
166176
unname() %>%
167177
as.vector()

inst/WORDLIST

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,4 @@
1+
CDI
12
CDM
23
CDMs
34
CMD
@@ -17,7 +18,9 @@ Gelman
1718
Gierl
1819
HDCM
1920
Kruskal's
21+
Kullback
2022
LCDM
23+
Leibler
2124
Liu
2225
MDM
2326
MacReady
@@ -33,6 +36,7 @@ PPMCs
3336
Psychometrika
3437
RMSEA
3538
RStan
39+
Roussos
3640
Rtools
3741
Rupp
3842
SRMSR
@@ -45,6 +49,7 @@ WAIC
4549
Watanabe
4650
Xcode
4751
Xin
52+
Xuming
4853
Youden's
4954
al
5055
att

0 commit comments

Comments
 (0)