Author: Tony

I wrote last time about the difR package (Magis, Beland, Tuerlinckx, & De Boeck, 2010) and how it doesn’t account for missing data in Mantel-Haenszel DIF analysis. I’ve noticed two more issues as I’ve continued testing the package (version 5.1).

1. The problem with Mantel-Haenszel also appears in the code for the standardization method, accessed via difR:::difStd, which calls difR:::stdPDIF. Look there and you’ll see base:::length used to obtain counts (e.g., number of correct/incorrect for focal and reference groups at a given score level). Missing data will throw off these counts. So, difR standardization and MH are only recommended with complete data.
2. In the likelihood ratio method, code for pseudo $R^2$ (used as a measure of DIF effect size) can lead to errors for some models. The code also seems to assume no missing data. More on these issues below.

DIF with the likelihood ratio method is performed using the difR:::difLogistic function, which ultimately calls difR:::Logistik to do the modeling (via glm) and calculate the $R^2$. The functions for calculating $R^2$ are embedded within the difR:::Logistik function.

R2 <- function(m, n) {
  1 - (exp(-m$null.deviance / 2 + m$deviance / 2))^(2 / n)
}
R2max <- function(m, n) {
  1 - (exp(-m$null.deviance / 2))^(2 / n)
}
R2DIF <- function(m, n) {
  R2(m, n) / R2max(m, n)
}

These functions capture $R^2$ as defined by Nagelkerke (1991), which is a modification to Cox and Snell (1989). When these are run via difR:::Logistik, the sample size n argument is set to the number of rows in the data set, which ignores missing data on a particular item. So, n will be inflated for items with missing data, and $R^2$ will be reduced (assuming a constant deviance).

In addition to the missing data issue, because of the way they’re written, these functions stretch the precision limits of R. In the R2max function specifically, the model deviance is first converted to a log-likelihood, and then a likelihood, before raising to 2/n. The problem is, large deviances correspond to very small likelihoods. A deviance of 500 gives us a likelihood of 7.175096e-66, which R can manage. But a deviance of 1500 gives us a likelihood of 0, which produces $R^2 = 1$.

The workaround is simple – avoid calculating likelihoods by rearranging terms. Here’s how I’ve written them in the epmr package.

r2_cox <- function(object, n = length(object$y)) {
  1 - exp((object\$deviance - object\$null.deviance) / n)
}
r2_nag <- function(object, n = length(object$y)) {
  r2_cox(object, n) / (1 - exp(-object$null.deviance / n))
}

And here are two examples that compare results from difR with epmr and DescTools. The first example shows how roughly 10% missing data reduces $R^2$ by as much as 0.02 when using difR. Data come from the verbal data set, included in difR.

# Load example data from the difR package
# See ?difR:::verbal for details
data("verbal", package = "difR")

# Insert missing data on first half of items
set.seed(42)
np <- nrow(verbal)
ni <- 24
na_index <- matrix(
  sample(c(TRUE, FALSE), size = np * ni / 2,
    prob = c(.1, .9), replace = TRUE),
  nrow = np, ncol = ni / 2)
verbal[, 1:(ni / 2)][na_index] <- NA

# Get R2 from difR
# verbal[, 26] is the grouping variable gender
verb_total <- rowSums(verbal[, 1:ni], na.rm = TRUE)
verb_difr <- difR:::Logistik(verbal[, 1:ni],
  match = verb_total, member = verbal[, 26],
  type = "udif")

# Fit the uniform DIF models by hand
# To test for DIF, we would compare these with base
# models, not fit here
verb_glm <- vector("list", ni)
for (i in 1:ni) {
  verbal_sub <- data.frame(y = verbal[, i],
    t = verb_total, g = verbal[, 26])
  verb_glm[[i]] <- glm(y ~ t + g, family = "binomial",
    data = verbal_sub)
}

# Get R2 from epmr and DescTools packages
verb_epmr <- sapply(verb_glm, epmr:::r2_nag)
verb_desc <- sapply(verb_glm, DescTools:::PseudoR2,
  which = "Nag")

# Compare
# epmr and DescTools match for all items
# difR matches for the last 12 items, but R2 on the
# first 12 are depressed because of missing data
verb_tab <- data.frame(item = 1:24,
  pct_na = apply(verbal[, 1:ni], 2, epmr:::summiss) / np,
  difR = verb_difr$R2M0, epmr = verb_epmr,
  DescTools = verb_desc)

This table shows results for items 9 through 16, the last four items with missing data and the first four with complete data.

The second example shows a situation where $R^2$ in the difR package comes to 1. Data are from the 2009 administration of PISA, included in epmr.

# Prep data from epmr::PISA09
# Vector of item names
rsitems <- c("r414q02s", "r414q11s", "r414q06s",
  "r414q09s", "r452q03s", "r452q04s", "r452q06s",
  "r452q07s", "r458q01s", "r458q07s", "r458q04s")

# Subset to USA and Canada
pisa <- subset(PISA09, cnt %in% c("USA", "CAN"))

# Get R2 from difR
pisa_total <- rowSums(pisa[, rsitems],
  na.rm = TRUE)
pisa_difr <- difR:::Logistik(pisa[, rsitems],
  match = pisa_total, member = pisa$cnt,
  type = "udif")

# Fit the uniform DIF models by hand
pisa_glm <- vector("list", length(rsitems))
for (i in seq_along(rsitems)) {
  pisa_sub <- data.frame(y = pisa[, rsitems[i]],
    t = pisa_total, g = pisa$cnt)
  pisa_glm[[i]] <- glm(y ~ t + g, family = "binomial",
    data = pisa_sub)
}

# Get R2 from epmr and DescTools packages
pisa_epmr <- sapply(pisa_glm, epmr:::r2_nag)
pisa_desc <- sapply(pisa_glm, DescTools:::PseudoR2,
  which = "Nag")

# Compare
pisa_tab <- data.frame(item = seq_along(rsitems),
  difR = pisa_difr$R2M0, epmr = pisa_epmr,
  DescTools = pisa_desc)

Here are the resulting $R^2$ for each package, across all items.

References

Cox, D. R. & Snell, E. J. (1989). The analysis of binary data. London: Chapman and Hall.

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.

Nagelkerke, N. J. D. (1991). A note on a general definition of the coefficient of determination. Biometrika, 78, 691-692.

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