Getting started with regDIF

Using LASSO regularization to evaluate measurement bias in a 2-parameter logistic IRT model.

In this example, data in ida were generated to mimic an integrative data analysis, where data are pooled across multiple studies and the measurement model is evaluated for between-study and within-study (e.g., gender, age) measurement bias. These data include 6 item responses (binary) and 3 background characteristics – namely, age (continuous, centered), gender (categorical, effect-coded), and study (categorical, effect-coded).

DIF was generated to be on items 2 (age, gender, study), 3 (age, gender, study), 4 (age), and 5 (gender, study), for both intercepts and slopes.

library(regDIF)
head(ida)

##   item1 item2 item3 item4 item5 item6 age gender study
## 1     0     0     0     0     0     0  -2     -1    -1
## 2     0     0     0     0     0     0   0     -1    -1
## 3     0     0     0     0     0     0   3     -1    -1
## 4     0     1     1     1     1     1   1     -1    -1
## 5     0     0     0     0     0     0  -2     -1    -1
## 6     1     0     0     0     0     0   1     -1    -1

The item response data must first be separated from the predictor data (background variables) before running regDIF(). A single value of the tuning parameter, tau = 2, is then fit to the data.

item.data <- regDIF::ida[,1:6]
pred.data <- regDIF::ida[,7:9]
fit <- regDIF(item.data, pred.data, tau = 2)

summary(fit)

## Call:
## regDIF(item.data = item.data, pred.data = pred.data, tau = 2)
## 
## Optimal model (out of 1):
##     tau     bic 
##    2.00 4133.96 
## 
## Non-zero DIF effects:

The summary() function shows that no DIF effects remain in the model. Only the latent variable parameters and base item parameters, which were not penalized at all, are estimated to be non-zero.

Now that the data have been properly specified in regDIF, a more thorough investigation of DIF is warranted. The regDIF() function defaults to estimating 100 values of the tuning parameter, starting with a value large enough to penalize all DIF effects to zero. However, for brevity, only 10 values of tau are specified with the num.tau argument.

fit2 <- regDIF(item.data, pred.data, num.tau = 10)

fit2

## Call:
## regDIF(item.data = item.data, pred.data = pred.data, num.tau = 10)
## 
## regDIF results:
##           tau      bic
## 1  0.65292855 4133.960
## 2  0.45857259 4118.947
## 3  0.30720781 4102.819
## 4  0.19346031 4074.727
## 5  0.11195620 4094.030
## 6  0.05732157 4106.946
## 7          NA       NA
## 8          NA       NA
## 9          NA       NA
## 10         NA       NA

By printing the model object, 10 rows of results appear, one for each value of the tuning parameter. The first thing to notice is that 4 rows are missing. This occurs because regDIF() automatically stops model-fitting when a small value of tau would produce a non-identified model. In focusing attention to the BIC column, it is evident that the smallest value occurs well before the model would be non-identified. This is an encouraging result. The summary() function may be used again, which, with multiple values of the tuning parameter fit to the data, produces non-zero DIF effects corresponding to the model with the minimum value of BIC.

summary(fit2)

## Call:
## regDIF(item.data = item.data, pred.data = pred.data, num.tau = 10)
## 
## Optimal model (out of 10):
##          tau          bic 
##    0.1934603 4074.7270000 
## 
## Non-zero DIF effects:
##    item4.int.age    item5.int.age item5.int.gender  item5.int.study 
##            0.263           -0.247           -0.656            0.349

A plot of the regularization path also shows the remaining DIF effects.

plot(fit2)

To produce other model results, the fit2 object contains lists of the impact (latent variable) parameters, base (intercept and slope) item parameters, and DIF parameters for all values of the tuning parameter. For instance, the impact parameters are printed below.

fit2$impact

##               [,1]   [,2]   [,3]   [,4]   [,5]   [,6] [,7] [,8] [,9] [,10]
## mean.age     0.706  0.725  0.742  0.711  0.724  0.884   NA   NA   NA    NA
## mean.gender -0.174 -0.134 -0.107 -0.072 -0.075 -0.107   NA   NA   NA    NA
## mean.study   0.961  0.957  0.957  0.895  0.890  0.943   NA   NA   NA    NA
## var.age      0.436  0.435  0.419  0.406  0.434  0.412   NA   NA   NA    NA
## var.gender  -0.102 -0.098 -0.102 -0.103 -0.045 -0.019   NA   NA   NA    NA
## var.study   -0.122 -0.126 -0.100 -0.091 -0.056  0.172   NA   NA   NA    NA

EAP scores and standard deviations may also be produced.

lapply(fit2$eap, head)

## $scores
##            [,1]       [,2]       [,3]       [,4]       [,5]       [,6] [,7]
## [1,] -1.8154255 -1.8707496 -1.9152125 -1.8490203 -1.8191629 -1.9636622   NA
## [2,] -1.3354071 -1.3647273 -1.3736359 -1.3435034 -1.2882957 -1.1873038   NA
## [3,] -1.2323978 -1.2210124 -1.1644336 -1.1753163 -1.1495209 -0.7501568   NA
## [4,]  2.0818664  2.0633270  2.0218976  1.9922185  1.9668574  1.7628487   NA
## [5,] -1.8154255 -1.8707496 -1.9152125 -1.8490203 -1.8191629 -1.9636622   NA
## [6,] -0.3815382 -0.3929121 -0.4053502 -0.3809949 -0.3741594 -0.2721668   NA
##      [,8] [,9] [,10]
## [1,]   NA   NA    NA
## [2,]   NA   NA    NA
## [3,]   NA   NA    NA
## [4,]   NA   NA    NA
## [5,]   NA   NA    NA
## [6,]   NA   NA    NA
## 
## $sd
##           [,1]      [,2]      [,3]      [,4]      [,5]      [,6] [,7] [,8] [,9]
## [1,] 0.7948612 0.7949724 0.7915591 0.7934734 0.7556749 0.6873146   NA   NA   NA
## [2,] 0.9321327 0.9315537 0.9199510 0.9140158 0.8823574 0.7836142   NA   NA   NA
## [3,] 1.1492558 1.1477015 1.1174331 1.0943602 1.0720908 0.9062620   NA   NA   NA
## [4,] 0.6987864 0.6964314 0.6891664 0.6781034 0.6556354 0.5806772   NA   NA   NA
## [5,] 0.7948612 0.7949724 0.7915591 0.7934734 0.7556749 0.6873146   NA   NA   NA
## [6,] 0.8568880 0.8546078 0.8430549 0.8292039 0.8073483 0.7058149   NA   NA   NA
##      [,10]
## [1,]    NA
## [2,]    NA
## [3,]    NA
## [4,]    NA
## [5,]    NA
## [6,]    NA

Finally, when data include a large number of items, observations, and predictors, regDIF() can run relatively slowly. An alternative approach, which yields much faster results, is to provide an observed proxy for the latent scores. In the case of binary data, this might be sum scores. Note that using observed proxy scores is identical to performing a multivariate regression, where the item responses are regressed on the proxy scores and background variables. (The proxy scores are simultaneously regressed on the background variables as well.)

fit3 <- regDIF(item.data, pred.data, prox.data = rowSums(item.data))

summary(fit3)

## Call:
## regDIF(item.data = item.data, pred.data = pred.data, prox.data = rowSums(item.data))
## 
## Optimal model (out of 100):
##          tau          bic 
##    0.2241118 3550.1030000 
## 
## Non-zero DIF effects:
## item3.int.gender    item4.int.age    item5.int.age item5.int.gender 
##            0.121            0.339           -0.192           -0.589 
##  item5.int.study item6.int.gender item2.slp.gender  item3.slp.study 
##            0.519            0.035            0.144           -0.062 
## item5.slp.gender 
##           -0.126

The results show more DIF on both the intercepts (Items 3 and 6) and slopes (Items 2, 3, and 5).

More modeling possibilities with regDIF.

In addition to LASSO, other penalty functions are possible with regDIF(). For instance, the elastic net penalty combines LASSO and ridge functions, which is useful when many correlated predictors are evaluated for DIF. The elastic net is controlled by a second tuning parameter, alpha, and defaults to alpha = 1, corresponding to the LASSO penalty. In contrast, alpha = 0 corresponds to the ridge penalty. When alpha is between 0 and 1, however, the elastic net is used to perform DIF selection. For brevity, observed proxy scores are used in all model fitting below.

fit_net <- regDIF(item.data, pred.data, prox.data = rowSums(item.data), alpha = .5)

summary(fit_net)

## Call:
## regDIF(item.data = item.data, pred.data = pred.data, prox.data = rowSums(item.data), 
##     alpha = 0.5)
## 
## Optimal model (out of 100):
##          tau          bic 
##    0.4377184 3562.5710000 
## 
## Non-zero DIF effects:
## item3.int.gender    item4.int.age    item5.int.age item5.int.gender 
##            0.104            0.277           -0.209           -0.466 
##  item5.int.study item6.int.gender item2.slp.gender  item3.slp.study 
##            0.369            0.033            0.119           -0.054 
## item5.slp.gender 
##           -0.099

The final elastic net results yield the same DIF effects as the LASSO results, although the amount of penalization is greater for elastic net (i.e., larger tau).

Because there are two tuning parameters – namely, tau and alpha – a grid search may be performed to identify the optimal degree of penalization.

alpha_vec <- seq(.1, 1, .1)
fit_net2 <- replicate(length(alpha_vec), NA, simplify = FALSE)

for(a in seq_along(alpha_vec)) {
  fit_net2[[a]] <- regDIF(item.data, pred.data, prox.data = rowSums(item.data), alpha = alpha_vec[a])
}

## Warning: Large increase in the number of DIF parameters from iteration 1 to 2.
##   Two Options:
##   1. Provide smaller differences between tau values.
##   2. Provide anchor item(s).

which.min(
  sapply(fit_net2, function(fit) min(fit$bic, na.rm = T))
)

## [1] 10

The minimum BIC across all models corresponds to the LASSO penalty model, i.e., alpha = 1.

Other penalty functions include the minimax concave penalty (MCP) and the group extensions of LASSO and MCP, which penalize the intercept and slope DIF effects in tandem. The group LASSO function is shown below.

fit_grp_mcp <- regDIF(item.data, pred.data, prox.data = rowSums(item.data), pen.type = "grp.mcp")

summary(fit_grp_mcp)

## Call:
## regDIF(item.data = item.data, pred.data = pred.data, prox.data = rowSums(item.data), 
##     pen.type = "grp.mcp")
## 
## Optimal model (out of 100):
##          tau          bic 
##    0.2960388 3543.9880000 
## 
## Non-zero DIF effects:
## item2.int.gender item3.int.gender    item4.int.age item5.int.gender 
##            0.030            0.099            0.344           -0.793 
##  item5.int.study item2.slp.gender item3.slp.gender    item4.slp.age 
##            0.950            0.112           -0.025            0.095 
## item5.slp.gender  item5.slp.study 
##           -0.499            0.170

Although the MCP results appear largely the same as the LASSO results, the group MCP function included both the intercept and slope for each background variable remaining in the final model.

In summary, the regDIF R package provides a flexible implementation of using regularization to identify DIF across multiple background characteristics.

Please reach out to wbelzak@gmail.com for any questions, and remember to cite regDIF in your work. Thank you kindly!

## 
## To cite regDIF in publications use:
## 
##   Belzak, W. C. M (2021). regDIF: Regularized Differential Item
##   Functioning. version 1.0.0. Duolingo. Pittsburgh, PA.
##   https://github.com/wbelzak/regDIF/
## 
## A BibTeX entry for LaTeX users is
## 
##   @Manual{,
##     title = {regDIF: Regularized Differential Item Functioning},
##     author = {William C. M. Belzak},
##     organization = {Duolingo},
##     address = {Pittsburgh, PA},
##     note = {version 1.0.0},
##     year = {2021},
##     url = {https://github.com/wbelzak/regDIF/},
##   }