difNLR: Generalized Logistic Regression Models for DIF and DDF Detection

Differential item functioning (DIF) and differential distractor functioning (DDF) are important topics in psychometrics, pointing to potential unfairness in items with respect to minorities or different social groups. Various methods have been proposed to detect these issues. The difNLR R package extends DIF methods currently provided in other packages by offering approaches based on generalized logistic regression models that account for possible guessing or inattention, and by providing methods to detect DIF and DDF among ordinal and nominal data. In the current paper, we describe implementation of the main functions of the difNLR package, from data generation, through the model fitting and hypothesis testing, to graphical representation of the results. Finally, we provide a real data example to bring the concepts together.


Introduction
Differential item functioning (DIF) is a phenomenon studied in the field of psychometrics which occurs when two subjects with the same ability or knowledge but from different social groups have a different probability of answering a specific test item correctly. DIF may signal bias and therefore potential unfairness of the measurement, so, DIF analysis should be used routinely in development and validation of educational and psychological tests .
In recent decades, various methods for DIF detection have been proposed by many different authors and examples of their (non-exhaustive) overviews can be found in Penfield and Camilli (2006) or Magis et al. (2010). Generally, DIF detection methods can be classified into two groups -those based on item response theory (IRT) models and those based on other approaches (referenced here as non-IRT). IRT models assume that ability is latent and needs to be estimated, whereas non-IRT methods use the total test score (or its standardization) for an approximation of ability level. IRT models can be found more conceptually and computationally challenging, while non-IRT methods are in general easier to implement and interpret.
Another distinction of DIF detection methods can be made with respect to the nature of the DIF that can be detected. Some methods, such as the Mantel-Haenszel test (Mantel and Haenszel, 1959) or IRT methods using the Rasch model (Rasch, 1993), can only detect so-called uniform DIF, a situation whereby an item is consistently advantageous to one group across all levels of ability. Other methods, for instance logistic regression (Swaminathan and Rogers, 1990) or methods using a 2 parameter logistic (2PL) IRT model, can also detect a non-uniform DIF, a situation when an item is advantageous to one group for some parts of ability levels while for other parts of ability levels the other group is in a more favourable position. DIF can be described as a difference in item parameters between the two groups. Uniform DIF can be then viewed as the difference in item difficulties (location of the inflexion point in the item characteristic curve) whereas non-uniform DIF also comprises a difference in item discrimination (slope at the inflexion point in the item characteristic curve). Most methods for DIF detection, with the important exception of methods based on 3PL and 4PL IRT models, are limited to testing differences in difficulty and discrimination and do not account for the possibility that an item can be guessed without necessary knowledge or incorrectly answered due to inattention, let alone the possibility to test for differences in related parameters.
There are several packages on the Comprehensive R Archive Network (CRAN, https://CRAN. Rproject.org/) which implement various methods for DIF detection. The difR package (Magis et al., 2010) includes a number of DIF detection methods among dichotomous items, inclusive of non-IRT methods as well as those based on IRT models. The DIFlasso (Schauberger, 2017) and DIFboost (Schauberger, 2016) packages implement penalty approach and boosting techniques in the IRT Rasch model. The GDINA package (Ma and de la Torre, 2019) offers the implementation of generalized deterministic inputs, noisy, and gate model. The mirt package (Chalmers, 2012) provides DIF detection based on various IRT models including those for dichotomous as well as ordinal or nominal items. The lordif package (Choi et al., 2016) brings an iterative method based on a mixture of ordinal logistic regression and an IRT approach for DIF detection. Finally, the psychotree package (Strobl et al., 2015) offers a method build on model-based recursive partitioning to detect latent groups of subjects exhibiting uniform DIF under the Rasch model. Mentioned packages do not allow for detecting differences in guessing and inattention, do not offer possibility of various matching criteria, and/or they assume IRT models which may become computationally demanding.
A potential gap in DIF methodology can be filled by generalizations of logistic regression models implemented within the difNLR package described here. These nonlinear extensions include estimations of pseudo-guessing parameters as proposed by . However, the difNLR package goes further by introducing a wide range of nonlinear regression models which cover both guessing and the possibility of inattention and allow for the testing of group differences in these parameters or their combinations. Moreover, in the case of ordinal data from rating scales or items allowing for partial credit, the difNLR package is able to detect differential use of the scale among groups by adjacent category logit or cumulative logit models. Finally, in the case of nominal data from multiple choice items, the package offers a multinomial model to test for so-called differential distractor functioning (DDF, Green et al., 1989), a situation when respondents with the same ability or knowledge but from different social groups are attracted differently by offered distractors (incorrect option choices). In addition, the difNLR package provides features such as item purification (see e.g., Lord, 1980) or corrections for multiple comparisons (see e.g., Kim and Oshima, 2013). This paper gives a detailed description of the difNLR package with functional examples, from data generation, through the model fitting, to interpretation of the results and their graphical representation, while separate parts are dedicated to features and troubleshooting. To bring the concepts together, we conclude by illustrative real data example connecting the usage of the main functions of the package.
The main advantage of non-IRT DIF detection methods implemented within the difNLR package is their flexibility with respect to the issue of guessing or inattention, while they retain pleasant properties such as low ratio of convergence issues when compared to IRT models . In addition to that, the difNLR provides DIF and DDF detection methods among ordinal and nominal data and thus offers a wide range of techniques to deal with an important topic of detecting potentially unfair items. For users who are new to R, interactive implementation of difNLR functions within the ShinyItemAnalysis package (Martinková and Drabinová, 2018) with toy datasets may be helpful.

Generalized logistic models for DIF and DDF detection
The class of generalized logistic models described here includes nonlinear regression models for DIF detection in dichotomously scored items, cumulative logit and adjacent category logit models for DIF detection in ordinal items, and a multinomial model for DDF detection in the case when items are nominal (e.g., multiple-choice).

Nonlinear regression models for binary items
Nonlinear regression models for DIF detection among binary items are extensions of the logistic regression method (Swaminathan and Rogers, 1990). These extensions account for the possibility that an item can be guessed, in other words correctly answered without possessing the necessary knowledge, i.e., lower asymptote of the item characteristic curve can be larger than zero . Similarly, these models account for a situation when an item is incorrectly answered by a person with high ability, e.g., due to inattention or lack of time, i.e., upper asymptote of the item characteristic curve can be lower than one.
The probability of a correct answer on item i by person p is then given by where X p is a variable describing observed knowledge or ability of person p and G p stands for their membership to social group (G p = 1 for focal group, usually seen as the disadvantaged one, G p = 0 for the reference group). Terms a iG p , b iG p ∈ R, and c iG p , d iG p ∈ [0, 1] represent discrimination, difficulty, probability of guessing, and a parameter related to the probability of inattention in item i, while they can differ for the reference and the focal group (denoted by the index G p ). Further, we use parametrization a iG p = a i + a iDIF G p , where a iDIF is a difference between the focal and the reference group in discrimination (analogously for other parameters). In other words, a i0 = a i for the reference group and a i1 = a i + a iDIF for the focal group. Thus, there are eight parameters for each item in total. For simplicity of the formulas, we stick with the notation of a iG p , b iG p , c iG p , and d iG p .
The class of models determined by equation (1) contains a wide range of methods for DIF detection. For instance, with d iG p = 1, we get a nonlinear model for the detection of DIF allowing for differential guessing in groups as presented by . Assuming moreover c iG p = 0, one can obtain a classical logistic regression model for the detection of uniform and non-uniform DIF as proposed by Swaminathan and Rogers (1990).
In contrast to IRT models, model (1) assumes that knowledge X p is represented by the total test score or its standardized form (Z-score) and thus the described method belongs to the class of non-IRT approaches. As such, model (1) can be seen as a proxy to the 4PL IRT model for DIF detection. While estimation of the asymptote parameters is notoriously challenging in IRT models, it was shown that generalized logistic models require a smaller sample size to be fitted while they keep pleasant properties in terms of power and rejection rates .
The difNLR package offers two approaches to estimate parameters of model (1). The first option is the nonlinear least squares method (Dennis et al., 1981;Ritz and Streibig, 2008), that is, minimization of the residual sums of squares (RSS) for item i with respect to item parameters: where n denotes number of respondents, y ip is the response of respondent p to the item i, x p is their (standardized) test score, and g p is their group membership. The second option is the maximum likelihood method. Likelihood function of item i is maximized with respect to item parameters.

Regression models for ordinal and nominal items
The logistic regression procedure which estimates the probability of the correct answer can be extended to estimate the probability of partial credit scores or option choices. When the responses are ordinal, DIF detection can be performed using the cumulative logit regression model (see e.g., Agresti, 2010, Chapter 3). For K + 1 outcome categories, i.e., Y ∈ {0, 1, . . . , K}, the cumulative probability is for k = 1, . . . , K, where P(Y ip ≥ 0|X p , G p ) = 1 and b ikG p = b ik + b iDIF G p . The parameter b iDIF can be interpreted as the difference in difficulty of item i between the reference and the focal group. The parameter b ik can be seen as category k specific difficulty of item i which is the same for both groups. The category probability is then given by Again, knowledge is represented by observed (standardized) total test score X p and therefore model (2) can be seen as a proxy to a graded response IRT model (Samejima, 1969).
Another approach for DIF detection in ordinal data is the adjacent category logit model (see e.g., Agresti, 2010, Chapter 4), which for K + 1 outcome categories is given by where k = 1, . . . , K and b ikG p = b ik + b iDIF G p . The category probability takes the following form: where k = 1, . . . , K and parameters for k = 0 are set to zero, i.e., a iG p (X p − b i0G p ) = 0 and hence . As such, an adjacent category logit model (3) can be seen as a proxy to the rating scale IRT model (Andrich, 1978). Both models, cumulative logit and adjacent category logit, are estimated by iteratively re-weighted least squares.
When responses are nominal (e.g., multiple choice), DDF detection can be performed with the multinomial model. Considering that k = 0, . . . , K possible option choices are offered, with k = 0 being the correct answer and other ones distractors, the probability of choosing distractor k is given by while for the correct answer k = 0 the parameters are set to zero, i.e., a i0G p (X p − b i0G p ) = 0 and thus . In contrast to ordinal models (2) and (3), discrimination a ikG p and difficulty b ikG p are category-specific and they may vary between groups. The parameters are estimated via neural networks.

Implementation in examples
The difNLR package provides implementation of the methods to detect DIF and DDF based on generalized logistic models. Specifically, a nonlinear regression model (1), cumulative logit model (2), adjacent category model (3), and multinomial model (4) are available within functions difNLR(), difORD(), and ddfMLR() (see Table 1). All three functions were designed to correspond to one of the most widely used R packages for DIF detection -difR (Magis et al., 2010).
Function Description difNLR() Performs DIF detection procedure for dichotomous data based on a nonlinear regression model (generalized logistic regression) and either likelihood-ratio or F test of the submodel. difORD() Performs DIF detection procedure for ordinal data based on either an adjacent category logit model or a cumulative logit model and likelihood ratio test of the submodel. ddfMLR() Performs DDF detection procedure for nominal data based on a multinomial log-linear regression model and likelihood ratio test of the submodel.

DIF detection
In this part, we discuss implementation and usage of the difNLR() function which offers a wide range of methods for DIF detection among dichotomous data based on a generalized logistic regression model (1). The full syntax of the difNLR() function is difNLR( Data, group, focal.name, model, constraints, type = "all", method = "nls", match = "zscore", anchor = NULL, purify = FALSE, nrIter = 10, test = "LR", alpha = 0.05, p.adjust.method = "none", start, initboot = TRUE, nrBo = 20 ) To detect DIF using the difNLR() function, the user always needs to provide four pieces of information: 1. the binary data set, 2. the group membership vector, 3. the indication of the focal group, and 4. the model.
Data. Data is a matrix or a data.frame with rows representing dichotomously scored respondents' answers (1 correct, 0 incorrect) and columns which correspond to the items. In addition, Data may contain the vector of group membership. If so, the group is a column identifier of the Data. Otherwise, the group must be a dichotomous vector of the same length as the number of rows (respondents) in Data. The name of the focal group is specified in focal.name argument.

Data generation.
To run a simulation study or to create an illustrative example, the difNLR package contains a data generator genNLR(), which can be used to generate dichotomous, ordinal, or nominal data. The type of items to be generated can be specified via itemtype argument: itemtype = "dich" for dichotomous items, "ordinal" for ordinal items, and "nominal" for nominal items.
For the generation of dichotomous items, discrimination and difficulty parameters need to be specified within a and b arguments in the form of matrices with two columns. The first column stands for the reference group and the second one for the focal group. Each row of matrices corresponds to one item. Additionally, one can provide guessing and inattention parameters via arguments c and d in the same way as for discriminations and difficulties. By default, values of guessing parameters are set to 0 in both groups, and the values of inattention parameters to 1 in both groups.
Distribution of the underlying latent trait is considered to be Gaussian. The user can specify its mean and standard deviation via arguments mu and sigma respectively. By default, mean is 0 and standard deviation is 1 and they are the same for both groups.
Furthermore, the user needs to provide a sample size (N) and the ratio of respondents in the reference and focal group (ratio). The latent trait for both groups is then generated and together with item parameters is used to generate item data. Output of the genNLR() function is a data.frame with items represented by columns and responses to them represented by rows. The last column is a group indicator, where 0 stands for a focal group and 1 indicates a reference group.
To illustrate generation of dichotomously scored items and to exemplify basic DIF detection with a difNLR() function, we create an example dataset. We choose discrimination a, difficulty b, guessing c, and inattention d parameters for 15 items. Parameters are then set the same for both groups. We generate dichotomous data with 500 observations in the reference group and 500 in the focal group. We assume that an underlying latent trait comes from a standard normal distribution for both groups (default setting). The output is a data.frame where the first 15 columns are dichotomously scored answers of 1,000 respondents and the last column is a group membership variable.

# discrimination
Model. The last necessary input of the difNLR() function is specification of the model to be estimated. This can be made by model argument. There are several predefined models, all of them based on the 4PL model stated in equation (1) (see Table 2).
We are now able to perform the basic DIF detection with a 4PL model for all the items on a generated example dataset DataDIF.

Model annotation Description
"4PL" 4PL model "4PLcdg", "4PLc" 4PL model with an inattention parameter set equal for the two groups "4PLcgd", "4PLd" 4PL model with a guessing parameter set equal for the two groups "4PLcgdg" 4PL model with a guessing and an inattention parameters set equal for the two groups "3PLd" 3PL model with an inattention parameter and c = 0 "3PLc", "3PL" 3PL model with a guessing parameter and d = 1 "3PLdg" 3PL model with an inattention parameter set equal for the two groups "3PLcg" 3PL model with a guessing parameter set equal for the two groups "2PL" Logistic regression model, i.e. c = 0 and d = 1 "1PL" 1PL model with a discrimination parameter set equal for the two groups "Rasch" 1PL model with a discrimination parameter fixed on value of 1 for the two groups (fit1 <-difNLR(DataDIF, groupDIF, focal.name = 1, model = "4PL")) Detection of all types of differential item functioning using generalized logistic regression model The output returns values of the test statistics for DIF detection, corresponding p-values, and set of items which are detected as functioning differently. All items (5, 8, 11, and 15) are correctly identified.
Estimates of parameters can be viewed with coef() method. Method coef() returns a list of parameters, which can be simplified to a matrix by setting simplify = TRUE. Each row then corresponds to one item and columns indicate parameters of the estimated model. round(coef(fit1, simplify = TRUE), 3) a b c d aDif bDif cDif dDif Item1 1.484 1.294 0.049 1.000 0.000 0.000 0.000 0.000 Item2 1.176 0.153 0.000 1.000 0.000 0.000 0.000 0.000 Item3 1.281 1.766 0.001 1.000 0.000 0.000 0.000 0.000 Item4 1.450 0.421 0.000 1.000 0.000 0.000 0.000 0.000 Item5 1.965 -1.147 0.000 0.868 -0.408 0.769 0.023 -0.006 Item6 1.458 -0.527 0.000 0.954 0.000 0.000 0.000 0.000 Item7 0.888 1.392 0.000 1.000 0.000 0.000 0.000 0.000 Item8 1.162 1.407 0.000 0.866 -0.117 0.974 0.007 0.134 Item9 1.482 -1.337 0.000 0.928 0.000 0.000 0.000 0.000 Item10 1.375 -0.570 0.007 0.967 0.000 0.000 0.000 0.000 Item11 1.071 -1.027 0.000 0.969 1.173 -0.499 0.000 0.011 Item12 1.051 1.560 0.080 1.000 0.000 0.000 0.000 0.000 Item13 1.009 1.348 0.084 1.000 0.000 0.000 0.000 0.000 Item14 1.093 1.659 0.141 1.000 0.000 0.000 0.000 0.000 Item15 0.875 -0.565 0.000 0.945 0.205 0.348 0.000 -0.142 The user can also print standard errors of the estimates using an option SE = TRUE. For example, estimated difference in difficulty between the reference and the focal groups in item 5 is 0.769 with standard error of 0.483. The difNLR() function provides a visual representation of the item characteristic curves using the ggplot2 package (Wickham, 2016) and its graphical environment. Curves are always based on results of a DIF detection procedure -when an item displays DIF, two curves are plotted, one for the reference and one for the focal group. Curves are accompanied by points representing empirical probabilities, i.e., proportions of correct answers with respect to the ability level and group membership. Size of the points is determined by the number of respondents at this ability level. Characteristic curves may simply be rendered with method plot() and by specifying items to be plotted. We show here characteristic curves for DIF items only ( Figure 1). plot(fit1, item = fit1$DIFitems) Besides predefined models (see Table 2), all parameters of the model can be further constrained using argument constraints specifying which parameters should be set equally for the two groups. For example, choice "ac" in 4PL model means that discrimination parameter a and pseudo-guessing parameter c are set equally for the two groups while the remaining parameters (b and d) are not. In addition, both arguments model and constraints are item-specific, meaning that a single value for all items can be introduced as well as a vector specifying the setting for each item. While the model specification can be challenging, this offers a wide range of models for DIF detection which goes hand in hand with the complexity of the offered method.
Furthermore, via type argument one can specify which type of DIF to test. Default option type = "all" allows one to test the difference in any parameter which is not constrained to be the same for both groups. Uniform DIF (difference in difficulty b only) can be tested by setting type = "udif", while nonuniform DIF (difference also in discrimination a) by setting type = "nudif". With the argument type = "both", the differences in both parameters (a and b) are tested. Moreover, to identify DIF in more detail, one can determine in which parameters the difference should be tested. The argument type is also item-specific.
The difNLR() function offers two techniques to estimate parameters of a generalized logistic regression model (1). With a default option method = "nls", nonlinear least square estimation is applied using a nls() function from the stats package. With an option method = "likelihood", the maximum likelihood method is used via an optim() function again from the stats package. Moreover, with an argument test, the user can specify what test of a submodel should be used to analyze DIF. The default option is the likelihood-ratio test.
Fitted values and residuals can be standardly calculated with methods fitted() and residuals(), again for all items or for those specified via item argument. This also holds for predicted values and method predict(). Predictions for any new respondents can be obtained by group and match arguments representing group membership and the value of matching criterion (e.g., standardized total score) of the new respondent. For example, with fit1 in item 5, new respondents with average performance (match = 0) have approximately a 22% lower probability of a correct answer if they come from a focal rather than reference group. predict(fit1, item = 5, group = c(0, 1), match = 0) item match group prob 1 Item5 0 0 0.7851739 2 Item5 0 1 0.5624883 This can also be observed when comparing item characteristic curves for the reference and the focal group in item 5 (see upper left Figure 1).

DIF detection among ordinal data
Here we show implementation and usage of the difORD() function which offers two models -cumulative logit model (2) and adjacent category logit model (3) to detect DIF among ordinal data. The full syntax of the difORD() function is difORD( Data, group, focal.name, model = "adjacent", type = "both", match = "zscore", anchor = NULL, purify = FALSE, nrIter = 10, p.adjust.method = "none", parametrization = "irt", alpha = 0.05 ) To detect DIF among ordinal data using the difORD() function, the user needs to provide four pieces of information: 1. the ordinal data set, 2. the group membership vector, 3. the indication of the focal group, and 4. the model to be fitted.

Data.
Data takes a similar format as used for the difNLR() function, however, rows represent ordinally scored respondents' answers instead of dichotomous. Specifications of group and focal.name remain the same.
Data generation. Data generator genNLR() is able to generate ordinal data using an adjacent category logit model (3) by setting itemtype = "ordinal". For polytomous items (ordinal or nominal), a and b have the form of matrices as well as for dichotomous data but each column now represents parameters of partial scores (or distractors). For example, to generate an item with 4 partial scores (i.e., 0-3), the user needs to provide 3 sets of discrimination and difficulty parameters. The parameters for minimal partial scores (i.e., 0; or correct answer in the case of nominal data) do not need to be specified because their probabilities are calculated as a complement to the sums of the partial scores probabilities. Guessing and inattention parameters are disregarded.
To illustrate usage of the difORD() function, we created an example ordinal dataset with 5 items, each scored with a range of 0-4. We first generated discrimination parameters a and difficulties b from a uniform distribution for partial scores k = 1, . . . , 4 for each item. In an adjacent category logit model (3), parameter b ik corresponds to an ability level for which the response categories k and k − 1 intersect in item i. For this reason and to create well-functioning items, parameters b ik are sorted so that b ik < b ik+1 . The parameters are set the same for both the reference and focal group.
set.seed(42) # discrimination a <-matrix(rep(runif(5, 0.25, 1), 8), ncol = 8) # difficulty b <-t(sapply(1:5, function(i) rep(sort(runif(4, -1, 1)), 2))) For the first two items we introduce uniform and non-uniform DIF respectively. Using parameters a and b of an adjacent category logit model, we generated ordinal data with a total sample size of 1,000 (500 observations per group). The first 5 columns of dataset DataORD represent ordinally scored items, while the last column represents a group membership variable. Model. The last input of the difORD() function which needs to be specified is model. It offers two possibilities. With an option model = "cumulative" a cumulative logit model (2) is fitted, while with an option model = "adjacent" (default) DIF detection is performed using an adjacent category logit model (3). The parameters for both models are estimated via vgam() function from the VGAM package (Yee, 2010).
DIF detection with cumulative logit model. In this part we exemplify usage of the difORD() function to fit a cumulative logit model for DIF detection among ordinal data. The group argument is introduced here by specifying the name of the group membership variable in DataORD dataset, i.e., group = "group". Knowledge is represented by observed standardized total score, i.e., standardized sum of all item scores.
(fit5 <-difORD(DataORD, group = "group", focal.name = 1, model = "cumulative")) Detection of both types of Differential Item Functioning for ordinal data using cumulative logit regression model Output provides test statistics for the likelihood ratio test, corresponding p-values, and the set of items which were detected as functioning differently. Items 1 and 2 are correctly identified as DIF items.
Similarly as for the difNLR() function, characteristic curves can be imaged with the plot() method. Besides characteristic curves, the method plot() for the cumulative logit model also offers the plot of cumulative probabilities. This can be achieved using plot.type = "cumulative", while with plot.type = "category" characteristic curves are shown. The plot of cumulative probabilities shows only 4 partial scores and does not show the cumulative probability of P(Y ip ≥ 0) since it is always equal to 1. Note that category probability of the highest score corresponds to its cumulative probability (Figure 3). plot(fit5, item = "Item1", plot.type = "cumulative") plot(fit5, item = "Item1", plot.type = "category") DIF detection with adjacent logit model. We illustrate here the fitting of an adjacent category logit model for DIF detection using the difORD() function. The group argument is now introduced by specifying the identifier of a group membership variable in Data (i.e., group = 6).
(fit6 <-difORD(DataORD, group = 6, focal.name = 1, model = "adjacent")) Detection of both types of Differential Item Functioning for ordinal data using adjacent category logit model Output again provides test statistics for the likelihood ratio test, corresponding p-values, and the set of items which were detected as functioning differently. Items 1 and 2 are correctly identified as DIF items. Characteristic curves can again be rendered using the plot() method (Figure 4).
Function difORD() offers the possibility to specify parametrization of regression coefficients with an argument parametrization. By default, IRT parametrization as stated in (2) and (3) is utilized using parametrization = "irt", but also classical intercept-slope parametrization (parametrization = "classic") with the effect of group membership and its interaction with matching criterion as covariates may be applied, i.e., b 0kG p + b 1G p X p instead of a iG p (X p − b ikG p )). The DIF detection is the same as with IRT parametrization, the only difference can be found in parameter estimates: Note that estimated discrimination for the reference group (parameter a) corresponds to the effect of matching criterion x, and in both cases their value is 1.776 for item 1. The same holds for the difference in discrimination and the effect of interaction between the matching criterion and group membership.

DDF detection among nominal data
Function ddfMLR() offers DDF detection among nominal data with the multinomial model (4). Here we illustrate its implementation and usage on a generated example. The full syntax of the ddfMLR() function is: ddfMLR( Data, group, focal.name, key, type = "both", match = "zscore", anchor = NULL, purify = FALSE, nrIter = 10, p.adjust.method = "none", parametrization = "irt", alpha = 0.05 ) To detect DDF among nominal data using the ddfMLR() function, the user needs to provide four pieces of information: 1. the unscored data set, 2. the key of correct answers, 3. the group membership vector, and 4. the indication of the focal group. The parameters are estimated via multinom() function from the nnet package (Venables and Ripley, 2002).
Data. The format of Data argument is similar to previously described functions. However, rows here represent respondents' unscored answers (e.g., in ABCD format). The group and focal.name is specified as in difNLR() or difORD() functions. Data generation. Data generator genNLR() can be used to generate nominal data using a multinomial model (4) by setting itemtype = "nominal". Specification of arguments a and b is the same as for ordinal items, however, it now represents parameters for distractors (incorrect answers).
To create an illustrative example dataset of nominal data, we must first generate discrimination a and difficulty b parameters from a uniform distribution for distractors of 10 items. The parameters are set the same for the reference and the focal group. For the first 5 items, we consider only two distractors (i.e., three item choices in total). For the last 5 items, we consider three distractors (i.e., four item choices in total). set.seed(42) # discrimination a <-matrix(rep(runif(30, -2, -0.5), 2), ncol = 6) a[1:5, c(3, 6)] <-NA # difficulty b <-matrix(rep(runif(30, -3, 1), 2), ncol = 6) b[1:5, c(3, 6)] <-NA For item 1, we introduce DDF by difference in discrimination and for item 6 by difference in difficulty. Finally, we generate nominal data with 500 observations in each group, i.e. 1,000 in total. The first 10 columns of the generated dataset DataDDF represent the unscored answers of respondents and the last column describes a group membership variable.
The correct answers in the generated dataset are denoted by A for each item; the key is hence a vector of As with a length of 10.
Now we have all the necessary inputs to fit a multinomial model (4) using the ddfMLR() function. The group argument is introduced here by specifying the name of group membership variable in Data (i.e., group = "group"). For the generated data, the total score is calculated as number of correct answers (i.e., number of As on a given row) and the matching criterion is then its standardized value (Z-score).
(fit7 <-ddfMLR(DataDDF, group = "group", focal.name = 1, key = rep("A", 10))) Detection of both types of Differential Distractor Functioning using multinomial log-linear regression model The output again summarizes the statistics of likelihood ratio test of a submodel, corresponding p-values, and the set of items identified as DDF. As expected, items 1 and 6 are detected as DDF.
Their characteristic curves can be displayed with a plot() method while the name of the reference and focal group can be modified via group.names argument ( Figure 5). This option is also available for functions difNLR() and difORD() and their plotting methods. plot(fit7, item = fit7$DDFitems, group.names = c("Group 1", "Group 2")) Similarly as for difNLR() and difORD(), item fit measures are offered via AIC(), BIC(), and logLik() S3 methods. Parameter estimates can be obtained using the method coef().

Further features
The difNLR covers user-friendly features that are common in standard DIF software -various matching criteria, anchor items, item purification, and p-value adjustments. These features are available for all main functions -difNLR(), difORD(), and ddfMLR().
Matching criterion. The models covered in the difNLR package are extensions of the logistic regression model described by Swaminathan and Rogers (1990) who considered the total test score as an observed ability X p . While this approach is well rooted in the psychometric research, note that it may lead to contradictions, e.g. when a nonzero item score is predicted for a respondent with a zero total test score. In the difNLR package, the default observed ability considered in all three main functions is the standardized total score. However, this estimate can be changed using the match argument. Besides default option "zscore" (standardized total score), it can also be the total test score (match = "score") or any numeric vector of the same length as the number of respondents. It is hence possible to use, for instance, latent trait estimates provided by some IRT models, or to use a pre-test score instead of the total score of the current test to be examined and analyze DIF in the context of item validity.

Anchor items and item purification.
Including DIF items into the calculation of matching criterion can lead to a potential bias and misidentification of DIF and non-DIF items. With an argument anchor, one can specify which items are supposed to be used for the calculation of matching criterion.
In the following examples, we go back to our generated dichotomous dataset DataDIF and the difNLR() function. For illustration, we take only items 1-6 and we apply some features with the 4PL model. Matching criterion is now calculated as a total test score based on items 1-6. Similar examples can also be illustrated with functions difORD() and ddfMLR().
We start with not specifying the anchor items. This indicates that any item can be considered as DIF one.
The process of including and omitting DIF and potentially unfair items could be demanding and time consuming. However, this process can be applied iteratively and automatically. This procedure is called item purification (Lord, 1980;Marco, 1977) and it has been shown that it can improve DIF detection. Item purification can be accessed with a purify argument. This can only be done when the matching criterion is either the total score or Z-score. The maximal number of iterations is determined by the nrIter argument, where the default value is 10.
In the initial step, items 5 and 6 were identified as DIF as it was shown with fit8a. The matching criterion was then calculated as the sum of the correct answers in items 1-4 as demonstrated by fit8b.
In the next step, only item 5 was identified as DIF and the matching criterion was based on items 1-4 and 6. The result of the DIF detection procedure was the same in the next step and the item purification process thus ended.

Multiple comparison corrections.
As the DIF detection procedure is done item by item, corrections for multiple comparisons may be considered (see Kim and Oshima, 2013). For example, applying Holm's adjustment (Holm, 1979) results in item 5 being detected as DIF.

Troubleshooting
In this section, we focus on several issues which can be encountered when fitting generalized logistic regression models and using the features offered in the difNLR package.

Convergence issues.
First, there is no guarantee that the estimation process in the difNLR() function will always end successfully. For instance, in the case of a small sample size, convergence issues may appear.
The easiest way to fix such issues is to specify different starting values. Various starting values can be applied via a start argument as a list with named numeric vectors as its elements. Each element needs to include values for the parameters a, b, c, and d of the reference group and the differences between reference and focal groups denoted by aDif, bDif, cDif, and dDif. However, there is no need to determine initial values manually. In the instance of convergence issues, the initial values are by default automatically re-calculated based on bootstrapped samples and applied only to models that failed to converge. This is also performed when starting values were initially introduced via a start argument. This feature can be turned off by setting initboot = FALSE. In such a case, no estimates are obtained for items that failed to converge. To demonstrate described situations, we now use a sample of our original simulated data set.
# sampled data set.seed(42) sam <-sample(1:1000, 420) # using re-calculation of starting values With an option initboot = TRUE in fit12a, starting values were re-calculated and no convergence issue occurred. When setting initboot = FALSE in fit12b we observed convergence failures in items 3 and 14.
The re-calculation process is by default performed up to twenty times, but the number of runs can be increased via the nrBo argument.
Another option is to apply the maximum likelihood method instead of nonlinear least squares to estimate parameters. There is no convergence issue in fit13 using the maximum likelihood estimation in contrast to fit12b and nonlinear least squares option.
Item purification. Issues may also occur when applying an item purification process. Although this is rare in practice, there is no guarantee that the process will end successfully. This can be observed, for instance, when we use DataDIF with the first 12 items only.
The maximum number of item purification iterations can be increased via the nrIter argument. However, in our example this would not necessarily lead to success as the process was not able to decide whether or not to include item 1 in the calculation of matching criterion. In this context, we advise considering such items as DIF to be on the safe side. As a general rule, any suspicious item should be reviewed by content experts. Not every DIF item is necessarily unfair, however even in such a case, understanding the reasons behind DIF may inform educators and help provide the best assessment and learning experience to all individuals involved.

Real data example
To illustrate the work with the difNLR package, we present a real data example from a study exploring the effect of student tracking on gains in learning competence . The LearningToLearn dataset presented here is available in the ShinyItemAnalysis package (Martinková and Hladká, 2020). It contains dichotomously scored responses of 782 students in total; 391 from a basic school (BS) track and 391 from a selective academic school (AS) track, whereas each student answered exactly the same test questions in Grade 6 and Grade 9. The sample was created from a larger dataset of 2,917 + 1,322 students using propensity score matching with the aim of obtaining similar baseline achievement and socio-economic characteristics of the students in both tracks. In addition, the matching algorithm required an exact match of the LtL total score in Grade 6 (i.e., exactly the same total score value in each matched pair of BS and AS students).
We demonstrate the three main functions of the difNLR package using items from the most difficult subscale 6, called Mathematical concepts, which consists of 8 multiple-choice items (6A, 6B, . . . , 6H) with four response options. We first test for DIF in Grade 6 using the difNLR() function to provide a detailed analysis of between-group differences on the baseline. Then, to get a more detailed insight into the differences in student gains after three years of education in two different tracks, we perform an analysis of the differential item functioning in change (DIF-C;  using student item responses in Grade 9 and then also using the nominal and ordinal data of item changes from Grade 6 to Grade 9. To identify items functioning differently on the baseline, we use a reduced dataset LtL6_gr6 which includes only student responses to eight subscale 6 items collected in Grade 6 and a group membership variable track. We further use a standardized total test score achieved in Grade 6 as an estimate of the students' ability. data("LearningToLearn", package = "ShinyItemAnalysis") # dichotomous items for Grade 6 LtL6_gr6 <-LearningToLearn [,c("track",paste0("Item6",LETTERS[1:8], "_6"))] head (LtL6_gr6)  track Item6A_6 Item6B_6 Item6C_6 Item6D_6 Item6E_6 Item6F_6 Item6G_6 Item6H_6  3453  AS  0  0  0  0  0  0  0  1  3456  AS  0  0  0  1  1  1  0  0  3458  AS  0  0  0  0  0  1  0  1  3464  AS  0  0  0  1  0  0  0  1  3474  AS  0  1  0  0  1  1  0  1  3491 AS 0 0 1 0 0 1 0 0 # standardized total score achieved in Grade 6 zscore6 <-LearningToLearn$score_6  hypothesized that, on the baseline, DIF found in this difficult subscale may possibly be caused by differences in guessing. To test this hypothesis, we can use the 3PL generalized logistic model which allows for detecting group differences in item difficulty, item discrimination, as well as the probability of guessing. fitex1 <-difNLR(Data = LtL6_gr6, group = "track", focal.name = "AS", model = "3PL", match = zscore6) fitex1$DIFitems [1] 8 Using the 3PL model, the eighth item (i.e., item 6H) was identified as DIF, while this DIF may have been caused by difference in any of the three item parameters.
In this eighth item (ninth column in the LtL_gr6 dataset), we now test a more specific hypothesis that DIF is caused by differences in guessing . This can be done by adding type = "c". fitex2 <-difNLR (Data = LtL6_gr6[,c(1,9)], group = "track", focal.name = "AS", model = "3PL", type = "c", match = zscore6) fitex2$DIFitems [1] 1 When specifically analyzing differences in the pseudo-guessing parameter c, we again detected a significant DIF in this item (6H). According to the model, for BS the estimated probability of guessing was close to 0.25, which is to be expected in multiple-choice items providing four response options. In contrast, for AS the probability of guessing was close to 0 ( Figure 6). plot(fitex2, item = fitex2$DIFitems) Figure 6: Characteristic curves of item 6H with DIF caused by differences in guessing between tracks in Grade 6.
We further perform DIF-C analysis to get a more detailed insight into the differences between AS and BS students in item gains. In contrast to DIF analysis, we use pre-test achievement in DIF-C analysis as the matching criterion while testing for group differences in the post-test, or in changes from pre-test to post-test.
In what follows, we demonstrate use of the difNLR(), difORD(), and ddfMLR() functions to provide evidence of DIF-C as described above. We first create a dataset LtL6_gr9 which consists of student responses in Grade 9 to subscale 6 items and the track variable.
fitex3 <-difNLR(Data = LtL6_gr9, group = "track", focal.name = "AS", model = "3PL", match = zscore6) fitex3$DIFitems [1] 1 2 The first two easier items of the subscale (i.e., items 6A and 6B) were identified as functioning differently in Grade 9 when accounting for the standardized total score in Grade 6. The differences between tracks can be illustrated by comparing the probabilities of answering item 6A correctly in Grade 9 by students with various baseline abilities: A low-performing student, average-performing student, and student performing above average, specifically, students with standardized total scores in Grade 6 equal to -1, 0, and 1. Method predict() also offers calculation of confidence intervals using delta method approach. Note that applying delta method in this case can result in confidence intervals slightly exceeding probability bounds 0 and 1. . Note that confidence intervals may overlap even in case when differential item functioning is significant among the two groups.
A second method for testing DIF-C proposed in  is based upon a multinomial model (4). To demonstrate the usage of ddfMLR() we use a restricted dataset of variables combining information about responses in Grades 6 and 9 into four response patterns: "00" indicating that the student did not answer correctly in either of the grades (did not improve); "10" meaning that the student answered correctly in Grade 6 but incorrectly in Grade 9 (deteriorated); "01" standing for a situation where the student answered incorrectly in Grade 6 but correctly in Grade 9 (improved); "11" describing the case where the student answered correctly in both grades (did not deteriorate). These variables denoted as "changes" can be extracted directly from the LearningToLearn dataset.
# nominal data for changes between 6th and 9th grade LtL6_change <-LearningToLearn[, c("track", paste0("Item6", LETTERS[1:8], "_changes"))] summary( We then fit a multinomial model (4) with the ddfMLR() function on LtL6_change data, again using the standardized total score achieved in Grade 6 as an estimate of a student's baseline ability, and by using the vector of patterns "11" as the argument key.
fitex4 <-ddfMLR(Data = LtL6_change, group = "track", focal.name = "AS", key = rep("11", 8), match = zscore6) fitex4$DDFitems [1] 2 5 Using the multinomial model, we identified items 2 and 5 (i.e., items 6B and 6E) as DIF-C items. In both items, the AS track was favoured in patterns "01" and "11", meaning that AS students had a higher probability of improving (pattern "01") or not deteriorating (pattern "11") than students with the same baseline standardized total score from the BS track. On the contrary, students with the same baseline standardized total score had a higher probability of not improving (pattern "00") or even deteriorating (pattern "10") in the BS track than in the AS track (Figure 7). plot(fitex4, item = fitex4$DDFitems) Finally, we introduce a third method for testing DIF-C using ordinal regression models. For this purpose, we create an ordinal dataset combining information about student responses in Grade 6 and Grade 9 into three categories: Score "0" means that the student deteriorated, i.e., answered correctly in Grade 6 but incorrectly in Grade 9; score "1" indicates no change in the accuracy of the answer, i.e., either correct or incorrect in both grades; and the best score "2" stands for improvement, i.e., the student did not answer correctly in Grade 6 but answered correctly in Grade 9.
fitex5 <-difORD(Data = LtL6_change_ord, group = "track", focal.name = "AS", model = "adjacent", match = zscore6) fitex5$DIFitems [1] 2 4 5 Using the ordinal model, we identified items 2, 4, and 5 (i.e., items 6B, 6D, and 6E) as functioning differently. In all the items, AS track was favoured in the improvement category "2", while BS students had a higher probability of deteriorating (category "0") and no change (category "1") when compared to AS students with the same baseline ability. The probability of belonging to category "2" was decreasing with an increasing baseline ability in items 6B and 6D, which is reasonable since the students with a higher baseline ability were more likely to have already answered these items correctly while in Grade 6. On the contrary, the probability of falling into the deteriorating category "0" was slightly increasing with baseline ability in these items, which can be explained as a tendency for over-thinking or inattention in very well-performing students. Trends were different in item 6E, where category probabilities seemed to be stable for all baseline ability levels, however, differences between the BS and AS track remained the same as for items 6B and 6D, favouring the AS track in the improving category "2" (Figure 8).
plot(fit5, item = fit5$DIFitems) Figure 8: Characteristic curves of DIF-C items 6B, 6D, and 6E using ordinal data and adjacent category logit model. All three approaches, using functions difNLR(), difORD(), and ddfMLR() jointly confirmed DIF-C in item 6B. Conclusions regarding other items were somewhat ambiguous, nevertheless, this is to be expected given the fact that each approach used different data (binary, nominal, and ordinal) which contained slightly different information.
Similar analysis can be conducted using an interactive environment of the ShinyItemAnalysis application (Martinková and Drabinová, 2018) which offers complex psychometric analysis including some functionalities implemented via difNLR package (see Figure 9). The datasets used in the real data example are included in the application; moreover, the user may upload and analyze interactively their own datasets as well.

Generalized logistic regression
Generalized logistic regression models are extensions of logistic regression method which account for possibility of guessing by allowing for nonzero lower asymptote -pseudoguessing (Drabinova & Martinkova, 2017) or upper asymptote lower than one -inattention . Similarly to logistic regression, its extensions also provide detection of uniform and non-uniform DIF by letting the difficulty parameter (uniform) and the discrimination parameter (non-uniform) differ for groups and by testing for difference in their values. Moreover, these extensions allow for testing differences in pseudo-guessing and inattention parameters and they can be seen as proxies of 3PL and 4PL IRT models for DIF detection.

Method specification
Here you can specify the assumed model. In 3PL and 4PL models, the abbreviations or mean that parameters or are assumed to be the same for both groups, otherwise they are allowed to differ. With type you can specify the type of DIF to be tested by choosing the parameters in which difference between groups should be tested. You can also select correction method for multiple comparison or item purification.
Finally, you may change the DIF matching variable. While matching on standardized total score is typical, upload of other DIF matching variable is possible in section Data. Using a pre-test (standardized) total score allows for testing differential item functioning in change (DIF-C) to provide proofs of instructional sensitivity (Martinkova et al., 2020), also see L e a r n i n g T o L e a r n 9 toy dataset. For selected item you can display plot of its characteristic curves and table of its estimated parameters with standard errors.
Plot with estimated DIF generalized logistic curve Points represent proportion of correct answer (empirical probabilities) with respect to the DIF matching variable. Their size is determined by count of respondents who achieved given level of DIF matching variable with respect to the group membership.

Correction method
Item purification  Figure 9: Generalized logistic model for DIF detection using ShinyItemAnalysis.

Summary
This article introduced the R package difNLR version 1.3.5 for DIF and DDF detection with extensions of a logistic regression model. The release version of the difNLR package is hosted on CRAN and the newest development version on GitHub, which can be accessed by devtools::install_github(" adelahladka/difNLR"). The current paper offered a description of the implementation of its three main functions difNLR(), difORD(), and ddfMLR() using simulated examples as well as a real data example with a longitudinal dataset LearningToLearn from the ShinyItemAnalysis package.
While the difR package already offers a classical logistic regression model for DIF detection via its function difLogistic(), the difNLR offers a wide range of methods based on extensions of the original model. The package is user-friendly since its model-fitting functions were designed to behave like the analogous functions in the difR package. Hence users who are already proficient in DIF detection among dichotomous items using R will find it easy to follow DIF and DDF detection via difNLR.
In addition, the difNLR package offers various S3 methods, including graphical representation of characteristic curves using the ggplot2 environment, various fit measures, extraction of fitted values and residuals as well as a prediction for a generalized logistic model.
Described models are also accessible via shiny application ShinyItemAnalysis which provides not only an interactive environment with a number of toy datasets including the real dataset analyzed here, but also a selected R code which can serve as a springboard for those new to R.
Future plans for the difNLR package include extensions of its functionality. In the difNLR() function, this comprises more options for estimation algorithms, as well as fine-tuning of the starting values to prevent convergence issues and to improve accuracy of the nonlinear models. We also plan to implement methods to extract fitted values and residuals for the difORD() and ddfMLR() functions together with their prediction via predict() method. Further, we would like to offer confidence and prediction intervals and their graphical representation for all three main functions. Finally, considered extensions include generalizations for longitudinal data with more time points.
In summary, the difNLR package provides various methods for DIF and DDF detection and can thus serve as a complex tool for detection of between-group differences in various contexts of educational and psychological tests and their items.