I’ve been using the difR package (Magis, Beland, Tuerlinckx, & De Boeck, 2010) to run differential item functioning (DIF) analysis in R. Here’s the package on CRAN.
https://cran.r-project.org/package=difR
I couldn’t get my own code to match the Mantel-Haenszel (MH) results from the difR package and it looks like it’s because there’s an issue in how the difR:::difMH function handles missing data. My code is on GitHub.
https://github.com/talbano/epmr/blob/master/R/difstudy.R
The MH DIF method is based on counts for correct vs incorrect responses in focal vs reference groups of test takers across levels of the construct (usually total scores). The code for difR:::difMH uses the length of a vector that is subset with logical indices to get the counts of test takers in each group. But missing data here will return NA in the logical comparisons, and NA isn’t omitted from length.
I’m pasting below the code from difR:::mantelHaenszel, which is called by difR:::difMH to run the MH analysis. Lines 19 to 33 all use length to find counts. This works fine with complete data, but as soon as someone has NA for an item score, captured in data[, item], they’ll figure into the count regardless of the comparisons being examined.
function (data, member, match = "score", correct = TRUE, exact = FALSE,
anchor = 1:ncol(data))
{
res <- resAlpha <- varLambda <- RES <- NULL
for (item in 1:ncol(data)) {
data2 <- data[, anchor]
if (sum(anchor == item) == 0)
data2 <- cbind(data2, data[, item])
if (!is.matrix(data2))
data2 <- cbind(data2)
if (match[1] == "score")
xj <- rowSums(data2, na.rm = TRUE)
else xj <- match
scores <- sort(unique(xj))
prov <- NULL
ind <- 1:nrow(data)
for (j in 1:length(scores)) {
Aj <- length(ind[xj == scores[j] & member == 0 &
data[, item] == 1])
Bj <- length(ind[xj == scores[j] & member == 0 &
data[, item] == 0])
Cj <- length(ind[xj == scores[j] & member == 1 &
data[, item] == 1])
Dj <- length(ind[xj == scores[j] & member == 1 &
data[, item] == 0])
nrj <- length(ind[xj == scores[j] & member == 0])
nfj <- length(ind[xj == scores[j] & member == 1])
m1j <- length(ind[xj == scores[j] & data[, item] ==
1])
m0j <- length(ind[xj == scores[j] & data[, item] ==
0])
Tj <- length(ind[xj == scores[j]])
if (exact) {
if (Tj > 1)
prov <- c(prov, c(Aj, Bj, Cj, Dj))
}
else {
if (Tj > 1)
prov <- rbind(prov, c(Aj, nrj * m1j/Tj, (((nrj *
nfj)/Tj) * (m1j/Tj) * (m0j/(Tj - 1))), scores[j],
Bj, Cj, Dj, Tj))
}
}
if (exact) {
tab <- array(prov, c(2, 2, length(prov)/4))
pr <- mantelhaen.test(tab, exact = TRUE)
RES <- rbind(RES, c(item, pr$statistic, pr$p.value))
}
else {
if (correct)
res[item] <- (abs(sum(prov[, 1] - prov[, 2])) -
0.5)^2/sum(prov[, 3])
else res[item] <- (abs(sum(prov[, 1] - prov[, 2])))^2/sum(prov[,
3])
resAlpha[item] <- sum(prov[, 1] * prov[, 7]/prov[,
8])/sum(prov[, 5] * prov[, 6]/prov[, 8])
varLambda[item] <- sum((prov[, 1] * prov[, 7] + resAlpha[item] *
prov[, 5] * prov[, 6]) * (prov[, 1] + prov[,
7] + resAlpha[item] * (prov[, 5] + prov[, 6]))/prov[,
8]^2)/(2 * (sum(prov[, 1] * prov[, 7]/prov[,
8]))^2)
}
}
if (match[1] != "score")
mess <- "matching variable"
else mess <- "score"
if (exact)
return(list(resMH = RES[, 2], Pval = RES[, 3], match = mess))
else return(list(resMH = res, resAlpha = resAlpha, varLambda = varLambda,
match = mess))
}
Here’s a very simplified example of the issue. The vector 1:4 is in place of the ind object in the mantelHaenszel function (created on line 17). The vector c(1, 1, NA, 0) is in place of data[, item] (e.g., on line 20). One person has a score of 0 on this item, and two have scores of 1, but length returns count 2 for item score 0 and 3 for item score 1 because the NA is not removed by default.
length((1:4)[c(1, 1, NA, 0) == 0]) ## [1] 2 length((1:4)[c(1, 1, NA, 0) == 1]) ## [1] 3
With missing data, the MH counts from difR:::mantelHaenszel will all be padded by the number of people with NA for their item score. It could be that the authors are accounting for this somewhere else in the code, but I couldn’t find it.
Here’s what happens to the MH results with some made up testing data. For 200 people taking a test with five items, I’m giving a boost on two items to 20 of the reference group test takers (to generate DIF), and then inserting NA for 20 people on one of those items. MH stats are consistent across packages for the first DIF item (item 4) but not the second (item 5).
# Number of items and people
ni <- 5
np <- 200
# Create focal and reference groups
groups <- rep(c("foc", "ref"), each = np / 2)
# Generate scores
set.seed(220821)
item_scores <- matrix(sample(0:1, size = ni * np,
replace = T), nrow = np, ncol = ni)
# Give 20 people from the reference group a boost on
# items 4 and 5
boost_ref_index <- sample((1:np)[groups == "ref"], 20)
item_scores[boost_ref_index, 4:5] <- 1
# Fix 20 scores on item 5 to be NA
item_scores[sample(1:np, 20), 5] <- NA
# Find total scores on the first three items,
# treated as anchor
total_scores <- rowSums(item_scores[, 1:3])
# Comparing MH stats, chi square matches for item 4
# with no NA but differs for item 5
epmr:::difstudy(item_scores, groups = groups,
focal = "foc", scores = total_scores, anchor_items = 1:3,
dif_items = 4:5, complete = FALSE)
##
## Differential Item Functioning Study
##
## item rn fn r1 f1 r0 f0 mh delta delta_abs chisq chisq_p ets_level
## 1 4 100 100 61 52 39 48 1.50 -0.946 0.946 1.58 0.2083 a
## 2 5 88 92 55 40 33 52 2.06 -1.701 1.701 4.84 0.0278 c
difR:::difMH(data.frame(item_scores), group = groups,
focal.name = "foc", anchor = 1:3, match = total_scores)
##
## Detection of Differential Item Functioning using Mantel-Haenszel method
## with continuity correction and without item purification
##
## Results based on asymptotic inference
##
## Matching variable: specified matching variable
##
## Anchor items (provided by the user):
##
## X1
## X2
## X3
##
##
## No p-value adjustment for multiple comparisons
##
## Mantel-Haenszel Chi-square statistic:
##
## Stat. P-value
## X4 1.5834 0.2083
## X5 4.8568 0.0275 *
##
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Detection threshold: 3.8415 (significance level: 0.05)
##
## Items detected as DIF items:
##
## X5
##
##
## Effect size (ETS Delta scale):
##
## Effect size code:
## 'A': negligible effect
## 'B': moderate effect
## 'C': large effect
##
## alphaMH deltaMH
## X4 1.4955 -0.9457 A
## X5 1.8176 -1.4041 B
##
## Effect size codes: 0 'A' 1.0 'B' 1.5 'C'
## (for absolute values of 'deltaMH')
##
## Output was not captured!
One more note, when reporting MH results, the difR package only uses the absolute delta values to assign ETS significance levels (A, B, C). You can see this in the difR:::print.MH function (not shown here). Usually, the MH approach also incorporates the p-value for the chi square (Zwick, 2012).
References
Magis, D., Beland, S, Tuerlinckx, F, & De Boeck, P. (2010). A general framework and an R package for the detection of dichotomous differential item functioning. Behavior Research Methods, 42, 847–862.
Zwick, R. (2012). A review of ETS differential item functioning assessment procedures: Flagging rules, minimum sample size requirements, and criterion refinement. Princeton, NJ: Educational Testing Service. https://files.eric.ed.gov/fulltext/EJ1109842.pdf