R/difNLR.R

Defines functions residuals.difNLR BIC.difNLR AIC.difNLR logLik.difNLR coef.difNLR predict.difNLR fitted.difNLR plot.difNLR print.difNLR difNLR

Documented in AIC.difNLR BIC.difNLR coef.difNLR difNLR fitted.difNLR logLik.difNLR plot.difNLR predict.difNLR residuals.difNLR

#' DIF detection using non-linear regression method.
#'
#' @aliases difNLR
#'
#' @description
#' Performs DIF detection procedure for dichotomous data based on non-linear
#' regression model (generalized logistic regression) and either
#' likelihood-ratio test, F-test, or Wald's test of a submodel.
#'
#' @param Data data.frame or matrix: dataset which rows represent scored
#'   examinee answers (\code{"1"} correct, \code{"0"} incorrect) and columns
#'   correspond to the items. In addition, \code{Data} can hold the vector of
#'   group membership.
#' @param group numeric or character: a dichotomous vector of the same length as
#'   \code{nrow(Data)} or a column identifier of \code{Data}.
#' @param focal.name numeric or character: indicates the level of \code{group}
#'   which corresponds to focal group.
#' @param model character: generalized logistic regression model to be fitted.
#'   See \strong{Details}.
#' @param constraints character: which parameters should be the same for both
#'   groups. Possible values are any combinations of parameters \code{"a"},
#'   \code{"b"}, \code{"c"}, and \code{"d"}. See \strong{Details}.
#' @param type character: type of DIF to be tested. Possible values are
#'   \code{"all"} for detecting difference in any parameter (default),
#'   \code{"udif"} for uniform DIF only (i.e., difference in difficulty
#'   parameter \code{"b"}), \code{"nudif"} for non-uniform DIF only (i.e.,
#'   difference in discrimination parameter \code{"a"}), \code{"both"} for
#'   uniform and non-uniform DIF (i.e., difference in parameters \code{"a"} and
#'   \code{"b"}), or combination of parameters \code{"a"}, \code{"b"},
#'   \code{"c"}, and \code{"d"}. Can be specified as a single value (for all
#'   items) or as an item-specific vector.
#' @param method character: method used to estimate parameters. The options
#'   are \code{"nls"} for non-linear least squares (default),
#'   \code{"likelihood"} for maximum likelihood method with \code{"L-BFGS-B"}
#'   algorithm, or \code{"irls"} for maximum likelihood method with iteratively
#'   reweighted least squares (available only for \code{model = "2PL"}).
#' @param match numeric or character: matching criterion to be used as an
#'   estimate of trait. Can be either \code{"zscore"} (default, standardized
#'   total score), \code{"score"} (total test score), or vector of the same
#'   length as number of observations in \code{Data}.
#' @param anchor numeric or character: specification of DIF free items. Either
#'   \code{NULL} (default), or a vector of item names (column names of
#'   \code{Data}), or item identifiers (integers specifying the column number)
#'   determining which items are currently considered as anchor (DIF free)
#'   items. Argument is ignored if \code{match} is not \code{"zscore"} or
#'   \code{"score"}.
#' @param purify logical: should the item purification be applied? (default is
#'   \code{FALSE}).
#' @param nrIter numeric: the maximal number of iterations in the item
#'   purification (default is 10).
#' @param test character: test to be performed for DIF detection. Can be either
#'   \code{"LR"} for likelihood ratio test of a submodel (default), \code{"W"}
#'   for Wald test, or \code{"F"} for F-test of a submodel.
#' @param alpha numeric: significance level (default is 0.05).
#' @param p.adjust.method character: method for multiple comparison correction.
#'   Possible values are \code{"holm"}, \code{"hochberg"}, \code{"hommel"},
#'   \code{"bonferroni"}, \code{"BH"}, \code{"BY"}, \code{"fdr"}, and
#'   \code{"none"} (default). For more details see
#'   \code{\link[stats]{p.adjust}}.
#' @param start numeric: initial values for estimation of parameters. If not
#'   specified, starting values are calculated with
#'   \code{\link[difNLR]{startNLR}} function. Otherwise, list with as many
#'   elements as a number of items. Each element is a named numeric vector of
#'   length 8 representing initial values for parameter estimation.
#'   Specifically, parameters \code{"a"}, \code{"b"}, \code{"c"}, and \code{"d"}
#'   are initial values for discrimination, difficulty, guessing, and
#'   inattention for reference group. Parameters \code{"aDif"}, \code{"bDif"},
#'   \code{"cDif"}, and \code{"dDif"} are then differences in these parameters
#'   between reference and focal group.
#' @param initboot logical: in case of convergence issues, should be starting
#'   values re-calculated based on bootstrapped samples? (default is
#'   \code{TRUE}; newly calculated initial values are applied only to
#'   items/models with convergence issues).
#' @param nrBo numeric: the maximal number of iterations for calculation of
#'   starting values using bootstrapped samples (default is 20).
#' @param sandwich logical: should be sandwich estimator used for covariance
#'   matrix of parameters when using \code{method = "nls"}? Default is
#'   \code{FALSE}.
#'
#' @usage
#' 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,
#'   sandwich = FALSE
#' )
#'
#' @details
#' DIF detection procedure based on non-linear regression is the extension of
#' the logistic regression procedure (Swaminathan & Rogers, 1990) accounting for
#' possible guessing and/or inattention when responding (Drabinova & Martinkova,
#' 2017; Hladka & Martinkova, 2020).
#'
#' The unconstrained form of 4PL generalized logistic regression model for
#' probability of correct answer (i.e., \eqn{y = 1}) is \deqn{P(y = 1) = (c +
#' cDif * g) + (d + dDif * g - c - cDif * g) / (1 + exp(-(a + aDif * g) * (x - b
#' - bDif * g))), } where \eqn{x} is by default standardized total score (also
#' called Z-score) and \eqn{g} is a group membership. Parameters \eqn{a},
#' \eqn{b}, \eqn{c}, and \eqn{d} are discrimination, difficulty, guessing, and
#' inattention. Terms \eqn{aDif}, \eqn{bDif}, \eqn{cDif}, and \eqn{dDif} then
#' represent differences between two groups (reference and focal) in relevant
#' parameters.
#'
#' This 4PL model can be further constrained by \code{model} and
#' \code{constraints} arguments. The arguments \code{model} and
#' \code{constraints} can be also combined. Both arguments can be specified as a
#' single value (for all items) or as an item-specific vector (where each
#' element correspond to one item).
#'
#' The \code{model} argument offers several predefined models. The options are as follows:
#' \code{Rasch} for 1PL model with discrimination parameter fixed on value 1 for both groups,
#' \code{1PL} for 1PL model with discrimination parameter fixed for both groups,
#' \code{2PL} for logistic regression model,
#' \code{3PLcg} for 3PL model with fixed guessing for both groups,
#' \code{3PLdg} for 3PL model with fixed inattention for both groups,
#' \code{3PLc} (alternatively also \code{3PL}) for 3PL regression model with guessing parameter,
#' \code{3PLd} for 3PL model with inattention parameter,
#' \code{4PLcgdg} for 4PL model with fixed guessing and inattention parameter for both groups,
#' \code{4PLcgd} (alternatively also \code{4PLd}) for 4PL model with fixed guessing for both groups,
#' \code{4PLcdg} (alternatively also \code{4PLc}) for 4PL model with fixed inattention for both groups,
#' or \code{4PL} for 4PL model.
#'
#' The \code{model} can be specified in more detail with \code{constraints}
#' argument which specifies what parameters should be fixed for both groups. For
#' example, choice \code{"ad"} means that discrimination (parameter \code{"a"})
#' and inattention (parameter \code{"d"}) are fixed for both groups and other
#' parameters (\code{"b"} and \code{"c"}) are not. The \code{NA} value for
#' \code{constraints} means no constraints.
#'
#' Missing values are allowed but discarded for item estimation. They must be
#' coded as \code{NA} for both, \code{Data} and \code{group} arguments.
#'
#' In case that model considers difference in guessing or inattention parameter,
#' the different parameterization is used and parameters with standard errors
#' are re-calculated by delta method. However, covariance matrices stick with
#' alternative parameterization.
#'
#' @return
#' The \code{difNLR()} function returns an object of class \code{"difNLR"}. The
#' output including values of the test statistics, p-values, and items detected
#' as function differently is displayed by the \code{print()} method.
#'
#' Object of class \code{"difNLR"} is a list with the following components:
#' \describe{
#'   \item{\code{Sval}}{the values of the test statistics.}
#'   \item{\code{nlrPAR}}{the estimates of final model.}
#'   \item{\code{nlrSE}}{the standard errors of estimates of final model.}
#'   \item{\code{parM0}}{the estimates of null model.}
#'   \item{\code{seM0}}{the standard errors of estimates of null model.}
#'   \item{\code{covM0}}{the covariance matrices of estimates of null model.}
#'   \item{\code{llM0}}{log-likelihood of null model.}
#'   \item{\code{parM1}}{the estimates of alternative model.}
#'   \item{\code{seM1}}{the standard errors of estimates of alternative model.}
#'   \item{\code{covM1}}{the covariance matrices of estimates of alternative model.}
#'   \item{\code{llM1}}{log-likelihood of alternative model.}
#'   \item{\code{DIFitems}}{either the column identifiers of the items which were detected as DIF, or
#'   \code{"No DIF item detected"} in case no item was detected as function differently.}
#'   \item{\code{model}}{fitted model.}
#'   \item{\code{constraints}}{constraints for the \code{model}.}
#'   \item{\code{type}}{character: type of DIF that was tested. If parameters were specified, the value is \code{"other"}.}
#'   \item{\code{types}}{character: the parameters (specified by user, \code{type} has value \code{"other"}) which were
#'   tested for difference.}
#'   \item{\code{p.adjust.method}}{character: method for multiple comparison correction which was applied.}
#'   \item{\code{pval}}{the p-values of the test.}
#'   \item{\code{adj.pval}}{the adjusted p-values of the test using \code{p.adjust.method}.}
#'   \item{\code{df}}{the degress of freedom of the test.}
#'   \item{\code{test}}{used test.}
#'   \item{\code{purification}}{\code{purify} value.}
#'   \item{\code{nrPur}}{number of iterations in item purification process. Returned only if \code{purify}
#'   is \code{TRUE}.}
#'   \item{\code{difPur}}{a binary matrix with one row per iteration of item purification and one column per item.
#'   \code{"1"} in i-th row and j-th column means that j-th item was identified as DIF in i-th iteration. Returned only
#'   if \code{purify} is \code{TRUE}.}
#'   \item{\code{conv.puri}}{logical: indicating whether item purification process converged before the maximal number
#'   \code{nrIter} of iterations. Returned only if \code{purify} is \code{TRUE}.}
#'   \item{\code{method}}{used estimation method.}
#'   \item{\code{conv.fail}}{numeric: number of convergence issues.}
#'   \item{\code{conv.fail.which}}{the identifiers of the items which did not converge.}
#'   \item{\code{alpha}}{numeric: significance level.}
#'   \item{\code{Data}}{the data matrix.}
#'   \item{\code{group}}{the vector of group membership.}
#'   \item{\code{group.names}}{names of groups.}
#'   \item{\code{match}}{matching criterion.}
#' }
#'
#' For an object of class \code{"difNLR"} several methods are available (e.g.,
#' \code{methods(class = "difNLR")}).
#'
#' @author
#' Adela Hladka (nee Drabinova) \cr
#' Institute of Computer Science of the Czech Academy of Sciences \cr
#' Faculty of Mathematics and Physics, Charles University \cr
#' \email{hladka@@cs.cas.cz} \cr
#'
#' Patricia Martinkova \cr
#' Institute of Computer Science of the Czech Academy of Sciences \cr
#' \email{martinkova@@cs.cas.cz} \cr
#'
#' Karel Zvara \cr
#' Faculty of Mathematics and Physics, Charles University \cr
#'
#' @references
#' Drabinova, A. & Martinkova, P. (2017). Detection of differential item
#' functioning with nonlinear regression: A non-IRT approach accounting for
#' guessing. Journal of Educational Measurement, 54(4), 498--517,
#' \doi{10.1111/jedm.12158}.
#'
#' Hladka, A. (2021). Statistical models for detection of differential item functioning. Dissertation thesis.
#' Faculty of Mathematics and Physics, Charles University.
#'
#' Hladka, A. & Martinkova, P. (2020). difNLR: Generalized logistic regression
#' models for DIF and DDF detection. The R Journal, 12(1), 300--323,
#' \doi{10.32614/RJ-2020-014}.
#'
#' Swaminathan, H. & Rogers, H. J. (1990). Detecting differential item
#' functioning using logistic regression procedures. Journal of Educational
#' Measurement, 27(4), 361--370, \doi{10.1111/j.1745-3984.1990.tb00754.x}
#'
#' @seealso
#' \code{\link[difNLR]{plot.difNLR}} for graphical representation of item characteristic curves and DIF statistics. \cr
#' \code{\link[difNLR]{coef.difNLR}} for extraction of item parameters with their standard errors. \cr
#' \code{\link[difNLR]{predict.difNLR}} for prediction. \cr
#' \code{\link[difNLR]{fitted.difNLR}} and \code{\link[difNLR]{residuals.difNLR}} for extraction of fitted
#' values and residuals. \cr
#' \code{\link[difNLR]{logLik.difNLR}}, \code{\link[difNLR]{AIC.difNLR}}, \code{\link[difNLR]{BIC.difNLR}}
#' for extraction of log-likelihood and information criteria. \cr
#'
#' \code{\link[stats]{p.adjust}} for multiple comparison corrections. \cr
#' \code{\link[stats]{nls}} for nonlinear least squares estimation. \cr
#' \code{\link[difNLR]{startNLR}} for calculation of initial values of fitting algorithms in \code{difNLR()}.
#'
#' @examples
#' # loading data
#' data(GMAT)
#' Data <- GMAT[, 1:20] # items
#' group <- GMAT[, "group"] # group membership variable
#'
#' # testing both DIF effects using likelihood-ratio test and
#' # 3PL model with fixed guessing for groups
#' (x <- difNLR(Data, group, focal.name = 1, model = "3PLcg"))
#' \dontrun{
#' # graphical devices
#' plot(x, item = x$DIFitems)
#' plot(x, item = "Item1")
#' plot(x, item = 1, group.names = c("Group 1", "Group 2"))
#' plot(x, plot.type = "stat")
#'
#' # coefficients
#' coef(x)
#' coef(x, SE = TRUE)
#' coef(x, SE = TRUE, simplify = TRUE)
#'
#' # fitted values
#' fitted(x)
#' fitted(x, item = 1)
#'
#' # residuals
#' residuals(x)
#' residuals(x, item = 1)
#'
#' # predicted values
#' predict(x)
#' predict(x, item = 1)
#'
#' # predicted values for new subjects
#' predict(x, item = 1, match = 0, group = 0)
#' predict(x, item = 1, match = 0, group = 1)
#'
#' # AIC, BIC, log-likelihood
#' AIC(x)
#' BIC(x)
#' logLik(x)
#'
#' # AIC, BIC, log-likelihood for the first item
#' AIC(x, item = 1)
#' BIC(x, item = 1)
#' logLik(x, item = 1)
#'
#' # testing both DIF effects using Wald test and
#' # 3PL model with fixed guessing for groups
#' difNLR(Data, group, focal.name = 1, model = "3PLcg", test = "W")
#'
#' # testing both DIF effects using F test and
#' # 3PL model with fixed guessing for groups
#' difNLR(Data, group, focal.name = 1, model = "3PLcg", test = "F")
#'
#' # testing both DIF effects using
#' # 3PL model with fixed guessing for groups and sandwich estimator
#' # of the covariance matrices
#' difNLR(Data, group, focal.name = 1, model = "3PLcg", sandwich = TRUE)
#'
#' # testing both DIF effects using LR test,
#' # 3PL model with fixed guessing for groups
#' # and Benjamini-Hochberg correction
#' difNLR(Data, group, focal.name = 1, model = "3PLcg", p.adjust.method = "BH")
#'
#' # testing both DIF effects using LR test,
#' # 3PL model with fixed guessing for groups
#' # and item purification
#' difNLR(Data, group, focal.name = 1, model = "3PLcg", purify = TRUE)
#'
#' # testing both DIF effects using 3PL model with fixed guessing for groups
#' # and total score as matching criterion
#' difNLR(Data, group, focal.name = 1, model = "3PLcg", match = "score")
#'
#' # testing uniform DIF effects using 4PL model with the same
#' # guessing and inattention
#' difNLR(Data, group, focal.name = 1, model = "4PLcgdg", type = "udif")
#'
#' # testing non-uniform DIF effects using 2PL model
#' difNLR(Data, group, focal.name = 1, model = "2PL", type = "nudif")
#'
#' # testing difference in parameter b using 4PL model with fixed
#' # a and c parameters
#' difNLR(Data, group, focal.name = 1, model = "4PL", constraints = "ac", type = "b")
#'
#' # testing both DIF effects using LR test,
#' # 3PL model with fixed guessing for groups
#' # using maximum likelihood estimation method with L-BFGS-B algorithm
#' difNLR(Data, group, focal.name = 1, model = "3PLcg", method = "likelihood")
#'
#' # testing both DIF effects using LR test and 2PL model
#' # using maximum likelihood estimation method with iteratively reweighted least squares algorithm
#' difNLR(Data, group, focal.name = 1, model = "2PL", method = "irls")
#' }
#'
#' @keywords DIF
#' @export
difNLR <- function(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, sandwich = FALSE) {
  if (any(type == "nudif" & model == "1PL")) {
    stop("Detection of non-uniform DIF is not possible with 1PL model.", call. = FALSE)
  }
  if (any(type == "nudif" & model == "Rasch")) {
    stop("Detection of non-uniform DIF is not possible with Rasch model.", call. = FALSE)
  }
  # input check
  # model
  if (missing(model)) {
    stop("'model' is missing", call. = FALSE)
  }
  # constraints
  if (missing(constraints)) {
    constraints <- NULL
  }
  # matching criterion
  if (!(match[1] %in% c("score", "zscore"))) {
    if (is.null(dim(Data))) {
      no <- length(Data)
    } else {
      no <- dim(Data)[1]
    }
    if (length(match) != no) {
      stop("Invalid value for 'match'. Possible values are 'zscore', 'score' or vector of the same length as number
           of observations in 'Data'.")
    }
  }
  # purification
  if (purify & !(match[1] %in% c("score", "zscore"))) {
    stop("Purification not allowed when matching variable is not 'zscore' or 'score'.", call. = FALSE)
  }
  # test
  if (!test %in% c("F", "LR", "W") | !is.character(type)) {
    stop("Invalid value for 'test'. Test for DIF detection must be either 'LR', 'W', or 'F'.", call. = FALSE)
  }
  # significance level
  if (alpha > 1 | alpha < 0) {
    stop("Argument 'alpha' must be between 0 and 1.", call. = FALSE)
  }
  # starting values
  if (missing(start)) {
    start <- NULL
  }
  # starting values with bootstrapped samples
  if (initboot) {
    if (nrBo < 1) {
      stop("The maximal number of iterations for calculation of starting values using bootstraped
           samples 'nrBo' need to be greater than 1.", call. = FALSE)
    }
  }
  # estimation method
  if (!(method %in% c("nls", "likelihood", "irls"))) {
    stop("Invalid value for 'method'. Estimation method must be either 'nls', 'likelihood', or 'irls'. ", call. = FALSE)
  }
  # internal NLR function
  internalNLR <- function() {
    if (length(group) == 1) {
      if (is.numeric(group)) {
        GROUP <- Data[, group]
        DATA <- as.data.frame(Data[, (1:dim(Data)[2]) != group])
        colnames(DATA) <- colnames(Data)[(1:dim(Data)[2]) != group]
      } else {
        GROUP <- Data[, colnames(Data) == group]
        DATA <- as.data.frame(Data[, colnames(Data) != group])
        colnames(DATA) <- colnames(Data)[colnames(Data) != group]
      }
    } else {
      GROUP <- group
      DATA <- as.data.frame(Data)
    }
    if (length(levels(as.factor(GROUP))) != 2) {
      stop("'group' must be binary vector", call. = FALSE)
    }
    if (is.matrix(DATA) | is.data.frame(DATA)) {
      if (!all(apply(DATA, 2, function(i) {
        length(levels(as.factor(i))) == 2
      }))) {
        stop("'Data' must be data frame or matrix of binary vectors.", call. = FALSE)
      }
      if (dim(DATA)[1] != length(GROUP)) {
        stop("'Data' must have the same number of rows as is length of vector 'group'.", call. = FALSE)
      }
    } else {
      if (is.vector(DATA)) {
        if (length(DATA) != length(GROUP)) {
          stop("'Data' must have the same number of rows as is length of vector 'group'.", call. = FALSE)
        }
        DATA <- as.data.frame(DATA)
      } else {
        stop("'Data' must be data frame or matrix of binary vectors.", call. = FALSE)
      }
    }
    group.names <- unique(GROUP)[!is.na(unique(GROUP))]
    if (group.names[1] == focal.name) {
      group.names <- rev(group.names)
    }
    GROUP <- as.numeric(as.factor(GROUP) == focal.name)

    if (length(match) == dim(DATA)[1]) {
      df <- data.frame(DATA, GROUP, match, check.names = F)
    } else {
      df <- data.frame(DATA, GROUP, check.names = F)
    }

    df <- df[complete.cases(df), ]

    if (dim(df)[1] == 0) {
      stop("It seems that your 'Data' does not include any subjects that are complete.",
        call. = FALSE
      )
    }

    GROUP <- df[, "GROUP"]
    DATA <- as.data.frame(df[, !(colnames(df) %in% c("GROUP", "match"))])
    colnames(DATA) <- colnames(df)[!(colnames(df) %in% c("GROUP", "match"))]

    if (length(match) > 1) {
      match <- df[, "match"]
    }

    # model
    if (length(model) == 1) {
      model <- rep(model, dim(DATA)[2])
    } else {
      if (length(model) != dim(DATA)[2]) {
        stop("Invalid length of 'model'. Model needs to be specified for each item or
             by single string. ", call. = FALSE)
      }
    }
    if (!all(model %in% c(
      "Rasch", "1PL", "2PL", "3PLcg", "3PLdg", "3PLc", "3PL", "3PLd",
      "4PLcgdg", "4PLcgd", "4PLd", "4PLcdg", "4PLc", "4PL"
    ))) {
      stop("Invalid value for 'model'.", call. = FALSE)
    }

    # type of DIF to be tested
    if (length(type) == 1) {
      type <- rep(type, dim(DATA)[2])
    } else {
      if (length(type) != dim(DATA)[2]) {
        stop("Invalid length of 'type'. Type of DIF needs to be specified for each item or
             by single string. ", call. = FALSE)
      }
    }
    tmptype <- type[!(type %in% c("udif", "nudif", "both", "all"))]
    if (length(tmptype) > 0) {
      tmptypes <- unique(unlist(sapply(
        1:length(tmptype),
        function(i) unique(unlist(strsplit(tmptype[i], split = "")))
      )))
      if (!all(tmptypes %in% letters[1:4])) {
        stop("Type of DIF to be tested not recognized. Only parameters 'a', 'b', 'c' or 'd' can be tested
             or 'type' must be one of predefined options: either 'udif', 'nudif', 'both', or 'all'.", call. = FALSE)
      }
    }

    # constraints
    if (!(is.null(constraints))) {
      if (length(constraints) == 1) {
        constraints <- rep(constraints, dim(DATA)[2])
      } else {
        if (length(constraints) != dim(DATA)[2]) {
          stop("Invalid length of 'constraints'. Constraints need to be specified for each item or
               by single string. ", call. = FALSE)
        }
      }
      constraints <- lapply(1:dim(DATA)[2], function(i) unique(unlist(strsplit(constraints[i], split = ""))))
      if (!all(sapply(1:dim(DATA)[2], function(i) {
        all(constraints[[i]] %in% letters[1:4]) | all(is.na(constraints[[i]]))
      }))) {
        stop("Constraints can be only combinations of parameters 'a', 'b', 'c', and 'd'. ", call. = FALSE)
      }
      if (any(!(type %in% c("udif", "nudif", "both", "all")))) {
        types <- as.list(rep(NA, dim(DATA)[2]))
        wht <- !(type %in% c("udif", "nudif", "both", "all"))
        types[wht] <- unlist(strsplit(type[wht], split = ""))
        if (any(sapply(1:dim(DATA)[2], function(i) (all(!is.na(constraints[[i]])) & (types[[i]] %in% constraints[[i]]))))) {
          stop("The difference in constrained parameters cannot be tested.", call. = FALSE)
        }
      } else {
        types <- NULL
      }
    } else {
      constraints <- as.list(rep(NA, dim(DATA)[2]))
      types <- NULL
    }
    # anchors
    if (!is.null(anchor)) {
      if (is.numeric(anchor)) {
        ANCHOR <- anchor
      } else {
        ANCHOR <- NULL
        for (i in 1:length(anchor)) {
          ANCHOR[i] <- (1:dim(DATA)[2])[colnames(DATA) == anchor[i]]
        }
      }
    } else {
      ANCHOR <- 1:dim(DATA)[2]
    }
    # parameterization
    parameterization <- ifelse(model %in% c(
      "3PLc", "3PL", "3PLd", "4PLcgd", "4PLd",
      "4PLcdg", "4PLc", "4PL"
    ),
    "alternative",
    "classic"
    )
    if (!is.null(start)) {
      if (length(start) != dim(DATA)[2]) {
        stop("Invalid value for 'start'. Initial values must be a list with as many elements
             as number of items in Data.", call. = FALSE)
      }
      if (!all(unique(unlist(lapply(start, names))) %in% c(letters[1:4], paste0(letters[1:4], "Dif")))) {
        stop("Invalid names in 'start'. Each element of 'start' needs to be a numeric vector
             with names 'a', 'b', 'c', 'd', 'aDif', 'bDif', 'cDif' and 'dDif'.", call. = FALSE)
      }
    }
    if (!purify | !(match[1] %in% c("zscore", "score")) | !is.null(anchor)) {
      PROV <- suppressWarnings(NLR(DATA, GROUP,
        model = model, constraints = constraints, type = type, method = method,
        match = match, anchor = ANCHOR, start = start, p.adjust.method = p.adjust.method,
        test = test, alpha = alpha, initboot = initboot, nrBo = nrBo, sandwich = sandwich
      ))
      STATS <- PROV$Sval
      ADJ.PVAL <- PROV$adjusted.pval
      significant <- which(ADJ.PVAL < alpha)

      nlrPAR <- nlrSE <- lapply(
        1:length(PROV$par.m1),
        function(i) {
          if (i %in% PROV$conv.fail.which) {
            structure(rep(NA, length(PROV$par.m1[[i]])),
              names = names(PROV$par.m1[[i]])
            )
          } else {
            structure(rep(0, length(PROV$par.m1[[i]])),
              names = names(PROV$par.m1[[i]])
            )
          }
        }
      )
      for (i in 1:length(PROV$par.m1)) {
        if (!(i %in% PROV$conv.fail.which)) {
          nlrPAR[[i]][names(PROV$par.m0[[i]])] <- PROV$par.m0[[i]]
          nlrSE[[i]][names(PROV$se.m0[[i]])] <- PROV$se.m0[[i]]
        }
      }

      if (length(significant) > 0) {
        DIFitems <- significant
        for (idif in 1:length(DIFitems)) {
          nlrPAR[[DIFitems[idif]]] <- PROV$par.m1[[DIFitems[idif]]]
          nlrSE[[DIFitems[idif]]] <- PROV$se.m1[[DIFitems[idif]]]
        }
      } else {
        DIFitems <- "No DIF item detected"
      }

      if (any(!(type %in% c("udif", "nudif", "both", "all")))) {
        wht <- (!(type %in% c("udif", "nudif", "both", "all")))
        type[wht] <- "other"
      }
      RES <- list(
        Sval = STATS,
        nlrPAR = nlrPAR, nlrSE = nlrSE,
        parM0 = PROV$par.m0, seM0 = PROV$se.m0, covM0 = PROV$cov.m0, llM0 = PROV$ll.m0,
        parM1 = PROV$par.m1, seM1 = PROV$se.m1, covM1 = PROV$cov.m1, llM1 = PROV$ll.m1,
        DIFitems = DIFitems,
        model = model, constraints = constraints, type = type, types = types,
        p.adjust.method = p.adjust.method, pval = PROV$pval, adj.pval = PROV$adjusted.pval,
        df = PROV$df, test = test,
        purification = purify,
        method = method,
        conv.fail = PROV$conv.fail, conv.fail.which = PROV$conv.fail.which,
        alpha = alpha,
        Data = DATA, group = GROUP, group.names = group.names, match = match
      )
    } else {
      nrPur <- 0
      difPur <- NULL
      noLoop <- FALSE
      prov1 <- suppressWarnings(NLR(DATA, GROUP,
        model = model, constraints = constraints, type = type, method = method,
        match = match, start = start, p.adjust.method = p.adjust.method, test = test,
        alpha = alpha, initboot = initboot, nrBo = nrBo, sandwich = sandwich
      ))
      stats1 <- prov1$Sval
      pval1 <- prov1$pval
      significant1 <- which(pval1 < alpha)

      if (length(significant1) == 0) {
        PROV <- prov1
        STATS <- stats1
        DIFitems <- "No DIF item detected"
        nlrPAR <- nlrSE <- structure(data.frame(matrix(0, ncol = dim(PROV$par.m1)[2], nrow = dim(PROV$par.m1)[1])),
          .Names = colnames(PROV$par.m1)
        )
        nlrPAR[, colnames(PROV$par.m0)] <- PROV$par.m0
        nlrSE[, colnames(PROV$par.m0)] <- PROV$se.m0
        noLoop <- TRUE
      } else {
        dif <- significant1
        difPur <- rep(0, length(stats1))
        difPur[dif] <- 1
        repeat {
          if (nrPur >= nrIter) {
            break
          } else {
            nrPur <- nrPur + 1
            nodif <- NULL
            if (is.null(dif)) {
              nodif <- 1:dim(DATA)[2]
            } else {
              for (i in 1:dim(DATA)[2]) {
                if (sum(i == dif) == 0) {
                  nodif <- c(nodif, i)
                }
              }
            }
            prov2 <- suppressWarnings(NLR(DATA, GROUP,
              model = model, constraints = constraints, type = type, method = method,
              match = match, anchor = nodif, start = start, p.adjust.method = p.adjust.method,
              test = test, alpha = alpha, initboot = initboot, nrBo = nrBo, sandwich = sandwich
            ))
            stats2 <- prov2$Sval
            pval2 <- prov2$pval
            significant2 <- which(pval2 < alpha)
            if (length(significant2) == 0) {
              dif2 <- NULL
            } else {
              dif2 <- significant2
            }
            difPur <- rbind(difPur, rep(0, dim(DATA)[2]))
            difPur[nrPur + 1, dif2] <- 1
            if (length(dif) != length(dif2)) {
              dif <- dif2
            } else {
              dif <- sort(dif)
              dif2 <- sort(dif2)
              if (sum(dif == dif2) == length(dif)) {
                noLoop <- TRUE
                break
              } else {
                dif <- dif2
              }
            }
          }
        }
        PROV <- prov2
        STATS <- stats2
        significant1 <- which(PROV$adjusted.pval < alpha)

        nlrPAR <- nlrSE <- lapply(
          1:length(PROV$par.m1),
          function(i) {
            if (i %in% PROV$conv.fail.which) {
              structure(rep(NA, length(PROV$par.m1[[i]])),
                names = names(PROV$par.m1[[i]])
              )
            } else {
              structure(rep(0, length(PROV$par.m1[[i]])),
                names = names(PROV$par.m1[[i]])
              )
            }
          }
        )
        for (i in 1:length(PROV$par.m1)) {
          if (!(i %in% PROV$conv.fail.which)) {
            nlrPAR[[i]][names(PROV$par.m0[[i]])] <- PROV$par.m0[[i]]
            nlrSE[[i]][names(PROV$se.m0[[i]])] <- PROV$se.m0[[i]]
          }
        }

        if (length(significant1) > 0) {
          DIFitems <- significant1
          for (idif in 1:length(DIFitems)) {
            nlrPAR[[DIFitems[idif]]] <- PROV$par.m1[[DIFitems[idif]]]
            nlrSE[[DIFitems[idif]]] <- PROV$se.m1[[DIFitems[idif]]]
          }
        } else {
          DIFitems <- "No DIF item detected"
        }
      }
      if (!is.null(difPur)) {
        rownames(difPur) <- paste0("Step", 0:(dim(difPur)[1] - 1))
        colnames(difPur) <- colnames(DATA)
      }

      if (any(!(type %in% c("udif", "nudif", "both", "all")))) {
        wht <- (!(type %in% c("udif", "nudif", "both", "all")))
        type[wht] <- "other"
      }
      RES <- list(
        Sval = STATS,
        nlrPAR = nlrPAR, nlrSE = nlrSE,
        parM0 = PROV$par.m0, seM0 = PROV$se.m0, covM0 = PROV$cov.m0, llM0 = PROV$ll.m0,
        parM1 = PROV$par.m1, seM1 = PROV$se.m1, covM1 = PROV$cov.m1, llM1 = PROV$ll.m1,
        DIFitems = DIFitems,
        model = model, constraints = constraints, type = type, types = types,
        p.adjust.method = p.adjust.method, pval = PROV$pval, adj.pval = PROV$adjusted.pval,
        df = PROV$df, test = test,
        purification = purify, nrPur = nrPur, difPur = difPur, conv.puri = noLoop,
        method = method,
        conv.fail = PROV$conv.fail, conv.fail.which = PROV$conv.fail.which,
        alpha = alpha,
        Data = DATA, group = GROUP, group.names = group.names, match = match
      )
    }
    if (PROV$conv.fail > 0) {
      warning(paste("Convergence failure in item", PROV$conv.fail.which, "\n"), call. = FALSE)
    }
    if (purify) {
      if (!noLoop) {
        warning(paste0(
          "Item purification process not converged after ",
          nrPur, ifelse(nrPur <= 1, " iteration.", " iterations."), "\n",
          "Results are based on the last iteration of the item purification.\n"
        ),
        call. = FALSE
        )
      }
    }
    class(RES) <- "difNLR"
    return(RES)
  }

  resToReturn <- internalNLR()
  return(resToReturn)
}

#' @export
print.difNLR <- function(x, ...) {
  cat(paste0(
    "Detection of ",
    ifelse(length(unique(x$type)) == 1,
      switch(unique(x$type),
        both = "both types of ",
        udif = "uniform ",
        nudif = "non-uniform ",
        all = "all types of ",
        other = ""
      ),
      ""
    ),
    "differential item functioning\n",
    "using generalized logistic regression model"
  ))
  cat(
    paste0(
      "\n\nGeneralized logistic regression ",
      switch(x$test,
        "F" = "F-test",
        "LR" = "likelihood ratio chi-square",
        "W" = "Wald test"
      ),
      " statistics",
      ifelse(length(unique(x$model)) == 1,
        paste0(
          "\nbased on",
          switch(unique(x$model),
            "Rasch" = "Rasch model",
            "1PL" = "1PL model",
            "2PL" = "2PL model",
            "3PL" = "3PL model",
            "3PLcg" = "3PL model with fixed guessing for groups",
            "3PLdg" = "3PL model with fixed inattention parameter for groups",
            "3PLc" = "3PL model",
            "3PLd" = "3PL model with inattention parameter",
            "4PLcgdg" = "4PL model with fixed guessing and inattention parameter for groups",
            "4PLcgd" = "4PL model with fixed guessing for groups",
            "4PLd" = "4PL model with fixed guessing for groups",
            "4PLcdg" = "4PL model with fixed inattention parameter for groups",
            "4PLc" = "4PL model with fixed inattention parameter for groups",
            "4PL" = "4PL model"
          )
        )
      )
    ),
    ifelse((!all(is.na(x$constraints)) & length(unique(x$constraints)) == 1),
      paste(
        "with constraints on",
        ifelse(length(unique(unlist(x$constraints))) == 1, "parameter", "parameters"),
        paste(unique(unlist(x$constraints)), collapse = ", "), "\n"
      ),
      "\n"
    )
  )
  cat(paste(
    "\nParameters were estimated using",
    switch(x$method,
      "nls" = "non-linear least squares\n",
      "likelihood" = "maximum likelihood method \nwith L-BFGS-B algorithm\n",
      "irls" = "maximum likelihood method \nwith iteratively reweighted least squares algorithm\n"
    )
  ))
  cat(paste0(
    "\nItem purification was",
    ifelse(x$purification, " ", " not "), "applied",
    ifelse(x$purification, paste0(
      " with ", x$nrPur,
      ifelse(x$nrPur <= 1, " iteration.", " iterations.")
    ), ""), "\n"
  ))
  if (x$p.adjust.method == "none") {
    cat("No p-value adjustment for multiple comparisons\n\n")
  } else {
    cat(paste(
      "Multiple comparisons made with",
      switch(x$p.adjust.method,
        holm = "Holm",
        hochberg = "Hochberg",
        hommel = "Hommel",
        bonferroni = "Bonferroni",
        BH = "Benjamini-Hochberg",
        BY = "Benjamini-Yekutieli",
        fdr = "FDR"
      ), "adjustment of p-values\n\n"
    ))
  }
  sign <- ifelse(is.na(x$adj.pval), " ",
    symnum(x$adj.pval,
      c(0, 0.001, 0.01, 0.05, 0.1, 1),
      symbols = c("***", "**", "*", ".", "")
    )
  )
  if (x$p.adjust.method == "none") {
    tab <- format(round(cbind(x$Sval, x$pval), 4))
    tab <- matrix(cbind(tab, sign), ncol = 3)
    colnames(tab) <- switch(x$test,
      "F" = c("F-value", "P-value", ""),
      "LR" = c("Chisq-value", "P-value", ""),
      "W" = c("Chisq-value", "P-value", "")
    )
  } else {
    tab <- format(round(cbind(x$Sval, x$pval, x$adj.pval), 4))
    tab <- matrix(cbind(tab, sign), ncol = 4)
    colnames(tab) <- switch(x$test,
      "F" = c("F-value", "P-value", "Adj. P-value", ""),
      "LR" = c("Chisq-value", "P-value", "Adj. P-value", ""),
      "W" = c("Chisq-value", "P-value", "Adj. P-value", "")
    )
  }

  rownames(tab) <- colnames(x$Data)

  print(tab, quote = FALSE, digits = 4, zero.print = F)
  cat("\nSignif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n")

  if (x$test == "F") {
    if (dim(unique(x$df))[1] == 1) {
      critical <- unique(qf(1 - x$alpha, x$df[, 1], x$df[, 2]))
    } else {
      critical <- NULL
    }
  } else {
    if (length(unique(x$df)) == 1) {
      critical <- unique(qchisq(1 - x$alpha, x$df))
    } else {
      critical <- NULL
    }
  }
  if (!is.null(critical)) {
    cat(paste0("\nDetection thresholds: ", round(critical, 4), " (significance level: ", x$alpha, ")"))
  }

  if (is.character(x$DIFitems)) {
    cat("\nNone of items is detected as DIF \n")
  } else {
    cat("\n\nItems detected as DIF items:")
    cat("\n", paste0(colnames(x$Data)[x$DIFitems], "\n"))
  }
}

#' ICC and test statistics plots for an object of \code{"difNLR"} class.
#'
#' @description Plot method for an object of \code{"difNLR"} class using \pkg{ggplot2}.
#'
#' Two types of plots are available. The first one is obtained by setting \code{plot.type = "cc"}
#' (default). The characteristic curves for an item specified in \code{item} argument are plotted.
#' Plotted curves represent the best model.
#'
#' The second plot is obtained by setting \code{plot.type = "stat"}. The  test statistics
#' (either LR-test, or F-test, depends on argument \code{test}) are displayed on the Y axis,
#' for each converged item. The detection threshold is displayed by a horizontal line and items
#' detected as DIF are printed with the red color. Only parameters \code{size} and \code{title}
#' are used.
#'
#' @param x an object of \code{"difNLR"} class.
#' @param plot.type character: type of plot to be plotted (either \code{"cc"} for characteristic curve
#' (default), or \code{"stat"} for test statistics).
#' @param item numeric or character: either character \code{"all"} to apply for all converged items (default),
#' or a vector of item names (column names of \code{Data}), or item identifiers (integers specifying
#' the column number).
#' @param draw.empirical logical: should empirical probabilities be plotted as points? Default value is \code{TRUE}.
#' @param draw.CI logical: should confidence intervals for predicted values be plotted? Default value is \code{FALSE}.
#' @param group.names character: names of reference and focal group.
#' @param ... other generic parameters for \code{plot()} function.
#'
#' @return For an option \code{plot.type = "stat"}, returns object of class \code{"ggplot"}. In case of
#' \code{plot.type = "cc"}, returns list of objects of class \code{"ggplot"}.
#'
#' Outputs can be edited and modified as standard \code{"ggplot"} object including colours, titles, shapes or linetypes.
#'
#' Note that option \code{draw.CI = TRUE} returns confidence intervals for predicted values as calculated by
#' \code{\link[difNLR]{predict.difNLR}}. Confidence intervals may overlap even in case that item functions differently.
#'
#' @author
#' Adela Hladka (nee Drabinova) \cr
#' Institute of Computer Science of the Czech Academy of Sciences \cr
#' Faculty of Mathematics and Physics, Charles University \cr
#' \email{hladka@@cs.cas.cz} \cr
#'
#' Patricia Martinkova \cr
#' Institute of Computer Science of the Czech Academy of Sciences \cr
#' \email{martinkova@@cs.cas.cz} \cr
#'
#' Karel Zvara \cr
#' Faculty of Mathematics and Physics, Charles University \cr
#'
#' @references
#' Drabinova, A. & Martinkova, P. (2017). Detection of differential item functioning with nonlinear regression:
#' A non-IRT approach accounting for guessing. Journal of Educational Measurement, 54(4), 498--517,
#' \doi{10.1111/jedm.12158}.
#'
#' Hladka, A. & Martinkova, P. (2020). difNLR: Generalized logistic regression models for DIF and DDF detection.
#' The R Journal, 12(1), 300--323, \doi{10.32614/RJ-2020-014}.
#'
#' Swaminathan, H. & Rogers, H. J. (1990). Detecting differential item functioning using logistic regression procedures.
#' Journal of Educational Measurement, 27(4), 361--370, \doi{10.1111/j.1745-3984.1990.tb00754.x}
#'
#' @seealso
#' \code{\link[difNLR]{difNLR}} for DIF detection among binary data using generalized logistic regression model. \cr
#' \code{\link[difNLR]{predict.difNLR}} for prediction.
#' \code{\link[ggplot2]{ggplot}} for general function to plot a \code{"ggplot"} object.
#'
#' @examples
#' \dontrun{
#' # loading data
#' data(GMAT)
#' Data <- GMAT[, 1:20] # items
#' group <- GMAT[, "group"] # group membership variable
#'
#' # testing both DIF effects using likelihood-ratio test and
#' # 3PL model with fixed guessing for groups
#' (x <- difNLR(Data, group, focal.name = 1, model = "3PLcg"))
#'
#' # item characteristic curves
#' plot(x)
#' plot(x, item = x$DIFitems)
#' plot(x, item = 1)
#' plot(x, item = "Item1")
#'
#' # item characteristic curves without empirical probabilities
#' plot(x, item = 1, draw.empirical = FALSE)
#'
#' # item characteristic curves without empirical probabilities but with CI
#' plot(x, item = 1, draw.empirical = FALSE, draw.CI = TRUE)
#'
#' # graphical devices - test statistics
#' plot(x, plot.type = "stat")
#' }
#' @export
plot.difNLR <- function(x, plot.type = "cc", item = "all",
                        group.names, draw.empirical = TRUE, draw.CI = FALSE, ...) {
  plotstat <- function(x) {
    if (x$conv.fail != 0) {
      if (length(x$conv.fail) == sum(!is.na(x$Sval))) {
        switch(x$test,
          "F" = stop(
            "None of items does converge. F-statistic values not plotted",
            call. = FALSE
          ),
          "LR" = stop(
            "None of items does converge. LR-statistic values not plotted",
            call. = FALSE
          )
        )
      } else {
        switch(x$test,
          "F" = warning(paste0("Item ", x$conv.fail.which,
            " does not converge. F-statistic value not plotted",
            collapse = "\n"
          ), call. = FALSE),
          "LR" = warning(paste0("Item ", x$conv.fail.which,
            " does not converge. LR-statistic value not plotted",
            collapse = "\n"
          ), call. = FALSE)
        )
      }
    }

    title <- "Non-linear regression DIF detection with none multiple comparison correction"

    n <- dim(x$Data)[1]
    if (x$test == "F") {
      if (dim(unique(x$df))[1] == 1) {
        Sval_critical <- unique(qf(1 - x$alpha, x$df[, 1], x$df[, 2]))
      } else {
        Sval_critical <- NULL
      }
    } else {
      if (length(unique(x$df)) == 1) {
        Sval_critical <- unique(qchisq(1 - x$alpha, x$df))
      } else {
        Sval_critical <- NULL
      }
    }
    if (is.null(Sval_critical)) {
      stop("Critical values are different for different items. Plot cannot be rendered.")
    }

    items <- setdiff(1:length(x$Sval), x$conv.fail.which)
    g <- as.factor(ifelse(x$Sval > Sval_critical, 1, 0))
    hv <- na.omit(as.data.frame(cbind(1:length(x$Sval), x$Sval, g)))
    hv$g <- as.factor(hv$g)
    hv$V1 <- as.factor(hv$V1)
    plot_stat <- ggplot2::ggplot(data = hv, ggplot2::aes(
      x = .data$V1, y = .data$V2,
      label = .data$V1,
      col = .data$g
    )) +
      # points
      ggplot2::geom_text() +
      ggplot2::scale_color_manual(values = c("black", "red")) +
      # critical value
      ggplot2::geom_hline(yintercept = Sval_critical) +
      # theme
      ggplot2::ggtitle(title) +
      ggplot2::labs(x = "Item", y = switch(x$test,
        "F" = "F-statistic",
        "LR" = "Chisq-statistic"
      )) +
      .plot.theme() +
      ggplot2::theme(
        plot.title = ggplot2::element_text(face = "bold", vjust = 1.5),
        axis.text.x = ggplot2::element_blank(),
        axis.ticks.x = ggplot2::element_blank(),
        legend.position = "none"
      )

    return(plot_stat)
  }

  if (missing(group.names)) {
    group.names <- x$group.names
    if (all(group.names %in% c(0, 1))) group.names <- c("Reference", "Focal")
  } else {
    if (length(group.names) > 2) {
      group.names <- group.names[1:2]
      warning("Only first two values for 'group.names' argument are used.", call. = FALSE)
    } else {
      if (length(group.names) < 2) {
        group.names <- c("Reference", "Focal")
        warning("Argument 'group.names' need to have length of two. Default value is used.", call. = FALSE)
      }
    }
  }

  plotCC <- function(x, item = item) {
    m <- length(x$nlrPAR)
    nams <- colnames(x$Data)

    if (inherits(item, "character")) {
      if (any(item != "all") & !all(item %in% nams)) {
        stop("Invalid value for 'item'. Item must be either character 'all', or
           numeric vector corresponding to column identifiers, or name of the item.",
          call. = FALSE
        )
      }
      if (any(item == "all")) {
        items <- 1:m
      } else {
        items <- which(nams %in% item)
      }
    } else {
      if (!inherits(item, "integer") & !inherits(item, "numeric")) {
        stop("Invalid value for 'item'. Item must be either character 'all', or
           numeric vector corresponding to column identifiers, or name of the item.",
          call. = FALSE
        )
      } else {
        if (!all(item %in% 1:m)) {
          stop("Invalid number for 'item'.",
            call. = FALSE
          )
        } else {
          items <- item
        }
      }
    }

    if (any(x$conv.fail.which %in% items)) {
      if (length(setdiff(items, x$conv.fail.which)) == 0) {
        stop(paste0("Item ", intersect(x$conv.fail.which, items), " does not converge. Characteristic curve not plotted.",
          collapse = "\n"
        ),
        call. = FALSE
        )
      } else {
        warning(paste0("Item ", intersect(x$conv.fail.which, items), " does not converge. Characteristic curve not plotted.",
          collapse = "\n"
        ),
        call. = FALSE
        )
        items <- setdiff(items, x$conv.fail.which)
      }
    }

    if (x$purification) {
      anchor <- c(1:m)[!c(1:m) %in% x$DIFitems]
    } else {
      anchor <- 1:m
    }

    if (length(x$match) > 1) {
      xlab <- "Matching criterion"
      match <- x$match
    } else {
      if (x$match == "score") {
        xlab <- "Total score"
        match <- rowSums(as.data.frame(x$Data[, anchor]))
      } else {
        xlab <- "Standardized total score"
        match <- scale(rowSums(as.data.frame(x$Data[, anchor])))
      }
    }
    xR <- match[x$group == 0]
    xF <- match[x$group == 1]

    max_match <- max(c(xR, xF), na.rm = TRUE)
    min_match <- min(c(xR, xF), na.rm = TRUE)

    plot_CC <- list()

    for (i in items) {
      df_point_0 <- data.frame(cbind(
        as.numeric(levels(as.factor(xR))),
        tapply(
          x$Data[x$group == 0, i],
          as.factor(xR), mean
        )
      ))
      df_point_1 <- data.frame(cbind(
        as.numeric(levels(as.factor(xF))),
        tapply(
          x$Data[x$group == 1, i],
          as.factor(xF), mean
        )
      ))
      df_point <- data.frame(rbind(
        cbind(df_point_0, group = "0"),
        cbind(df_point_1, group = "1")
      ))
      colnames(df_point) <- c("match", "prob", "group")
      rownames(df_point) <- NULL
      df_point$size <- c(table(xR), table(xF))

      TITLE <- colnames(x$Data)[i]

      df_line <- predict.difNLR(x,
        item = i,
        match = rep(
          seq(floor(min_match),
            ceiling(max_match),
            length.out = 100
          ),
          2
        ),
        group = rep(c(0, 1), each = 100),
        interval = "confidence"
      )

      g <- ggplot2::ggplot() +
        ggplot2::geom_line(
          data = df_line, ggplot2::aes(
            x = .data$match, y = .data$prob,
            col = .data$group, linetype = .data$group
          ),
          linewidth = 0.8
        )
      if (draw.empirical) {
        g <- g +
          ggplot2::geom_point(
            data = df_point, ggplot2::aes(
              x = .data$match, y = .data$prob,
              col = .data$group, fill = .data$group,
              size = .data$size
            ),
            shape = 21, alpha = 0.5
          )
      }
      if (draw.CI) {
        g <- g +
          ggplot2::geom_ribbon(
            data = df_line, ggplot2::aes(
              x = .data$match, y = .data$prob,
              ymin = .data$lwr.conf, ymax = .data$upr.conf,
              col = .data$group, fill = .data$group
            ),
            alpha = 0.25
          )
      }

      g <- g +
        # adjusting colours
        ggplot2::scale_fill_manual(values = c("dodgerblue2", "goldenrod2"), labels = group.names) +
        ggplot2::scale_colour_manual(values = c("dodgerblue2", "goldenrod2"), labels = group.names) +
        ggplot2::scale_linetype_manual(values = c("solid", "dashed"), labels = group.names) +
        ggplot2::ggtitle(TITLE) +
        ggplot2::labs(x = xlab, y = "Probability of correct answer") +
        ggplot2::scale_y_continuous(limits = c(0, 1)) +
        # theme
        .plot.theme() +
        # legend
        .plot.theme.legend() +
        ggplot2::guides(
          size = ggplot2::guide_legend(title = "Count", order = 1),
          colour = ggplot2::guide_legend(title = "Group", order = 2),
          fill = ggplot2::guide_legend(title = "Group", order = 2),
          linetype = ggplot2::guide_legend(title = "Group", order = 2)
        )
      plot_CC[[i]] <- g
    }
    plot_CC <- plot_CC[items]
    names(plot_CC) <- nams[items]
    return(plot_CC)
  }
  # checking input
  if (!(plot.type %in% c("cc", "stat"))) {
    stop("Invalid value for 'plot.type'. Possible values of 'plot.type' is 'cc' or 'stat'.",
      call. = FALSE
    )
  } else {
    if (plot.type == "cc") {
      plotCC(x, item = item)
    } else {
      plotstat(x)
    }
  }
}

#' @rdname residuals.difNLR
#' @aliases residuals.difNLR
#' @export
fitted.difNLR <- function(object, item = "all", ...) {
  # checking input
  m <- length(object$nlrPAR)
  nams <- colnames(object$Data)

  if (inherits(item, "character")) {
    if (any(item != "all") & !all(item %in% nams)) {
      stop("Invalid value for 'item'. Item must be either character 'all', or
           numeric vector corresponding to column identifiers, or name of the item.",
        call. = FALSE
      )
    }
    if (any(item == "all")) {
      items <- 1:m
    } else {
      items <- which(nams %in% item)
    }
  } else {
    if (!inherits(item, "integer") & !inherits(item, "numeric")) {
      stop("Invalid value for 'item'. Item must be either character 'all', or
           numeric vector corresponding to column identifiers, or name of the item.",
        call. = FALSE
      )
    } else {
      if (!all(item %in% 1:m)) {
        stop("Invalid number for 'item'.",
          call. = FALSE
        )
      } else {
        items <- item
      }
    }
  }

  if (any(object$conv.fail.which %in% items)) {
    if (length(setdiff(items, object$conv.fail.which)) == 0) {
      stop(paste0("Item ", intersect(object$conv.fail.which, items), " does not converge. ",
        collapse = "\n"
      ),
      call. = FALSE
      )
    } else {
      warning(paste0("Item ", intersect(object$conv.fail.which, items), " does not converge. NA produced.",
        collapse = "\n"
      ),
      call. = FALSE
      )
      ITEMS <- setdiff(items, object$conv.fail.which)
    }
  } else {
    ITEMS <- items
  }

  PAR <- data.frame(
    a = rep(1, m), b = 0, c = 0, d = 1,
    aDif = 0, bDif = 0, cDif = 0, dDif = 0
  )
  for (i in ITEMS) {
    PAR[i, names(object$nlrPAR[[i]])] <- object$nlrPAR[[i]]
  }

  # data
  if (length(object$match) > 1) {
    match <- object$match
  } else {
    if (object$match == "score") {
      match <- c(apply(object$Data, 1, sum))
    } else {
      match <- c(scale(apply(object$Data, 1, sum)))
    }
  }

  # fitted values
  FV <- as.list(rep(NA, m))
  FV[ITEMS] <- lapply(ITEMS, function(i) {
    .gNLR(
      match, object$group,
      PAR[i, "a"], PAR[i, "b"],
      PAR[i, "c"], PAR[i, "d"],
      PAR[i, "aDif"], PAR[i, "bDif"],
      PAR[i, "cDif"], PAR[i, "dDif"]
    )
  })
  FV <- do.call(cbind, FV)
  colnames(FV) <- colnames(object$Data)
  rownames(FV) <- rownames(object$Data)
  FV <- FV[, items]

  return(FV)
}

#' Predicted values for an object of \code{"difNLR"} class.
#'
#' @description S3 method for predictions from the model used in the
#'   object of \code{"difNLR"} class.
#'
#' @param object an object of \code{"difNLR"} class.
#' @param item numeric or character: either character \code{"all"} to
#'   apply for all converged items (default), or a vector of item
#'   names (column names of \code{Data}), or item identifiers
#'   (integers specifying the column number).
#' @param match numeric: matching criterion for new observations.
#' @param group numeric: group membership for new observations.
#' @param interval character: type of interval calculation, either
#'   \code{"none"} (default) or \code{"confidence"} for confidence
#'   interval.
#' @param level numeric: confidence level.
#' @param ... other generic parameters for \code{predict()} function.
#'
#' @author
#' Adela Hladka (nee Drabinova) \cr
#' Institute of Computer Science of the Czech Academy of Sciences \cr
#' Faculty of Mathematics and Physics, Charles University \cr
#' \email{hladka@@cs.cas.cz} \cr
#'
#' Patricia Martinkova \cr
#' Institute of Computer Science of the Czech Academy of Sciences \cr
#' \email{martinkova@@cs.cas.cz} \cr
#'
#' Karel Zvara \cr
#' Faculty of Mathematics and Physics, Charles University \cr
#'
#' @references
#' Drabinova, A. & Martinkova, P. (2017). Detection of differential item functioning with nonlinear regression:
#' A non-IRT approach accounting for guessing. Journal of Educational Measurement, 54(4), 498--517,
#' \doi{10.1111/jedm.12158}.
#'
#' Hladka, A. & Martinkova, P. (2020). difNLR: Generalized logistic regression models for DIF and DDF detection.
#' The R Journal, 12(1), 300--323, \doi{10.32614/RJ-2020-014}.
#'
#' Swaminathan, H. & Rogers, H. J. (1990). Detecting differential item functioning using logistic regression procedures.
#' Journal of Educational Measurement, 27(4), 361--370, \doi{10.1111/j.1745-3984.1990.tb00754.x}
#'
#' @seealso
#' \code{\link[difNLR]{difNLR}} for DIF detection among binary data using generalized logistic regression model. \cr
#' \code{\link[stats]{predict}} for generic function for prediction.
#'
#' @examples
#' \dontrun{
#' # loading data
#' data(GMAT)
#' Data <- GMAT[, 1:20] # items
#' group <- GMAT[, "group"] # group membership variable
#'
#' # testing both DIF effects using likelihood-ratio test and
#' # 3PL model with fixed guessing for groups
#' (x <- difNLR(Data, group, focal.name = 1, model = "3PLcg"))
#'
#' # predicted values
#' summary(predict(x))
#' predict(x, item = 1)
#' predict(x, item = "Item1")
#'
#' # predicted values for new observations - average score
#' predict(x, item = 1, match = 0, group = 0) # reference group
#' predict(x, item = 1, match = 0, group = 1) # focal group
#'
#' # predicted values for new observations - various z-scores and groups
#' new.match <- rep(c(-1, 0, 1), 2)
#' new.group <- rep(c(0, 1), each = 3)
#' predict(x, item = 1, match = new.match, group = new.group)
#'
#' # predicted values for new observations with confidence intervals
#' predict(x, item = 1, match = new.match, group = new.group, interval = "confidence")
#' predict(x, item = c(2, 4), match = new.match, group = new.group, interval = "confidence")
#' }
#'
#' @export
predict.difNLR <- function(object, item = "all", match, group, interval = "none", level = 0.95, ...) {
  # checking input
  m <- length(object$nlrPAR)
  nams <- colnames(object$Data)

  if (inherits(item, "character")) {
    if (any(item != "all") & !all(item %in% nams)) {
      stop("Invalid value for 'item'. Item must be either character 'all', or
           numeric vector corresponding to column identifiers, or name of the item.",
        call. = FALSE
      )
    }
    if (any(item == "all")) {
      items <- 1:m
    } else {
      items <- which(nams %in% item)
    }
  } else {
    if (!inherits(item, "integer") & !inherits(item, "numeric")) {
      stop("Invalid value for 'item'. Item must be either character 'all', or
           numeric vector corresponding to column identifiers, or name of the item.",
        call. = FALSE
      )
    } else {
      if (!all(item %in% 1:m)) {
        stop("Invalid number for 'item'.",
          call. = FALSE
        )
      } else {
        items <- item
      }
    }
  }

  if (missing(match)) {
    if (length(object$match) > 1) {
      match <- object$match
    } else {
      if (object$match == "score") {
        match <- c(apply(object$Data, 1, sum))
      } else {
        match <- c(scale(apply(object$Data, 1, sum)))
      }
    }
  }

  if (missing(group)) {
    group <- object$group
  }
  if (length(match) != length(group)) {
    if (length(match) == 1) {
      match <- rep(match, length(group))
    } else if (length(group) == 1) {
      group <- rep(group, length(match))
    } else {
      stop("Arguments 'match' and 'group' must be of the same length.",
        call. = FALSE
      )
    }
  }

  if (any(object$conv.fail.which %in% items)) {
    if (length(setdiff(items, object$conv.fail.which)) == 0) {
      stop(paste0("Item ", intersect(object$conv.fail.which, items), " does not converge. ",
        collapse = "\n"
      ),
      call. = FALSE
      )
    } else {
      warning(paste0("Item ", intersect(object$conv.fail.which, items), " does not converge. NA produced.",
        collapse = "\n"
      ),
      call. = FALSE
      )
      ITEMS <- setdiff(items, object$conv.fail.which)
    }
  } else {
    ITEMS <- items
  }

  if (!interval %in% c("none", "confidence")) {
    warning(
      "Only confidence interval is supported. ",
      call. = FALSE
    )
    interval <- "none"
  }

  PAR <- data.frame(
    a = rep(1, m), b = 0, c = 0, d = 1,
    aDif = 0, bDif = 0, cDif = 0, dDif = 0
  )
  for (i in ITEMS) {
    PAR[i, names(object$nlrPAR[[i]])] <- object$nlrPAR[[i]]
  }

  # predicted values
  NEW <- lapply(1:m, function(i) {
    if (i %in% ITEMS) {
      .gNLR(
        match, group, PAR[i, "a"], PAR[i, "b"],
        PAR[i, "c"], PAR[i, "d"],
        PAR[i, "aDif"], PAR[i, "bDif"],
        PAR[i, "cDif"], PAR[i, "dDif"]
      )
    } else {
      NA
    }
  })

  if (interval == "confidence") {
    DELTA_new <- as.list(rep(NA, m))
    DELTA_new[ITEMS] <- lapply(ITEMS, function(i) {
      attr(.delta.gNLR(
        match, group, PAR[i, "a"], PAR[i, "b"],
        PAR[i, "c"], PAR[i, "d"],
        PAR[i, "aDif"], PAR[i, "bDif"],
        PAR[i, "cDif"], PAR[i, "dDif"]
      ), "gradient")
    })

    VCOV_par <- matrix(
      0,
      ncol = 8, nrow = 8,
      dimnames = list(
        c("a", "b", "c", "d", "aDif", "bDif", "cDif", "dDif"),
        c("a", "b", "c", "d", "aDif", "bDif", "cDif", "dDif")
      )
    )
    VCOV_par <- lapply(1:m, function(i) VCOV_par)
    VCOV_par[ITEMS] <- lapply(ITEMS, function(i) {
      if (i %in% object$DIFitems) {
        tmp <- colnames(object$covM1[[i]])
        VCOV_par[[i]][tmp, tmp] <- object$covM1[[i]]
      } else {
        tmp <- colnames(object$covM0[[i]])
        VCOV_par[[i]][tmp, tmp] <- object$covM0[[i]]
      }
      return(VCOV_par[[i]])
    })

    SE_new <- const <- lwr <- upp <- lapply(1:m, function(i) c())
    SE_new[ITEMS] <- lapply(ITEMS, function(i) {
      sqrt(diag(DELTA_new[[i]] %*% VCOV_par[[i]] %*% t(DELTA_new[[i]])))
    })

    alpha <- 1 - level
    n <- dim(object$Data)[[1]]
    d <- sapply(1:m, function(i) {
      ifelse(i %in% object$DIFitems, length(object$parM1[[i]]), length(object$parM0[[i]]))
    })
    df <- n - d

    const[ITEMS] <- lapply(ITEMS, function(i) SE_new[[i]] * qt(1 - alpha / 2, df = df)[i])

    lwr[ITEMS] <- lapply(ITEMS, function(i) NEW[[i]] - const[[i]])
    upp[ITEMS] <- lapply(ITEMS, function(i) NEW[[i]] + const[[i]])

    res <- lapply(ITEMS, function(i) {
      data.frame(cbind(
        item = colnames(object$Data)[i],
        match, group,
        prob = NEW[[i]],
        lwr.conf = lwr[[i]],
        upr.conf = upp[[i]]
      ))
    })
    res <- as.data.frame(do.call(rbind, res))
    res$match <- as.numeric(paste(res$match))
    res$prob <- as.numeric(paste(res$prob))
    res$lwr.conf <- as.numeric(paste(res$lwr.conf))
    res$upr.conf <- as.numeric(paste(res$upr.conf))
  } else {
    res <- lapply(ITEMS, function(i) {
      data.frame(cbind(
        item = colnames(object$Data)[i],
        match, group,
        prob = NEW[[i]]
      ))
    })
    res <- as.data.frame(do.call(rbind, res))
    res$match <- as.numeric(paste(res$match))
    res$prob <- as.numeric(paste(res$prob))
  }

  return(res)
}

#' Extract model coefficients from an object of \code{"difNLR"} class.
#'
#' @description S3 method for extracting model coefficients from an
#'   object of \code{"difNLR"} class.
#' @aliases coefficients.difNLR
#'
#' @param object an object of \code{"difNLR"} class.
#' @param SE logical: should the standard errors of estimated
#'   parameters be also returned? (default is \code{FALSE}).
#' @param simplify logical: should the estimated parameters be
#'   simplified to a matrix? (default is \code{FALSE}).
#' @param IRTpars logical: should the estimated parameters be returned
#'   in IRT parameterization? (default is \code{TRUE}).
#' @param CI numeric: level of confidence interval for parameters,
#'   default is \code{0.95} for 95\% confidence interval.
#' @param ... other generic parameters for \code{coef()} function.
#'
#' @author
#' Adela Hladka (nee Drabinova) \cr
#' Institute of Computer Science of the Czech Academy of Sciences \cr
#' Faculty of Mathematics and Physics, Charles University \cr
#' \email{hladka@@cs.cas.cz} \cr
#'
#' Patricia Martinkova \cr
#' Institute of Computer Science of the Czech Academy of Sciences \cr
#' \email{martinkova@@cs.cas.cz} \cr
#'
#' Karel Zvara \cr
#' Faculty of Mathematics and Physics, Charles University \cr
#'
#' @references
#' Drabinova, A. & Martinkova, P. (2017). Detection of differential item functioning with nonlinear regression:
#' A non-IRT approach accounting for guessing. Journal of Educational Measurement, 54(4), 498--517,
#' \doi{10.1111/jedm.12158}.
#'
#' Hladka, A. & Martinkova, P. (2020). difNLR: Generalized logistic regression models for DIF and DDF detection.
#' The R Journal, 12(1), 300--323, \doi{10.32614/RJ-2020-014}.
#'
#' Swaminathan, H. & Rogers, H. J. (1990). Detecting differential item functioning using logistic regression procedures.
#' Journal of Educational Measurement, 27(4), 361--370, \doi{10.1111/j.1745-3984.1990.tb00754.x}
#'
#' @seealso
#' \code{\link[difNLR]{difNLR}} for DIF detection among binary data using generalized logistic regression model. \cr
#' \code{\link[stats]{coef}} for generic function extracting model coefficients.
#'
#' @examples
#' \dontrun{
#' # loading data
#' data(GMAT)
#' Data <- GMAT[, 1:20] # items
#' group <- GMAT[, "group"] # group membership variable
#'
#' # testing both DIF effects using likelihood-ratio test and
#' # 3PL model with fixed guessing for groups
#' (x <- difNLR(Data, group, focal.name = 1, model = "3PLcg"))
#'
#' # estimated parameters
#' coef(x)
#' # includes standard errors
#' coef(x, SE = TRUE)
#' # includes standard errors and simplifies to matrix
#' coef(x, SE = TRUE, simplify = TRUE)
#' # intercept-slope parameterization
#' coef(x, IRTpars = FALSE)
#' # intercept-slope parameterization, simplifies to matrix, turn off confidence intervals
#' coef(x, IRTpars = FALSE, simplify = TRUE, CI = 0)
#' }
#' @export
coef.difNLR <- function(object, SE = FALSE, simplify = FALSE, IRTpars = TRUE, CI = 0.95, ...) {
  if (!inherits(SE, "logical")) {
    stop("Invalid value for 'SE'. 'SE' needs to be logical. ",
      call. = FALSE
    )
  }
  if (!inherits(simplify, "logical")) {
    stop("Invalid value for 'simplify'. 'simplify' needs to be logical. ",
      call. = FALSE
    )
  }
  if (!inherits(IRTpars, "logical")) {
    stop("Invalid value for 'IRTpars'. 'IRTpars' needs to be logical. ",
      call. = FALSE
    )
  }
  if (CI > 1 | CI < 0 | !is.numeric(CI)) {
    stop("Invalid value for 'CI'. 'CI' must be numeric value between 0 and 1. ",
      call. = FALSE
    )
  }

  m <- dim(object$Data)[2]
  nams <- colnames(object$Data)

  if (IRTpars) {
    pars <- object$nlrPAR
    se <- object$nlrSE
  } else {
    nlrPAR <- ifelse(1:m %in% object$DIFitems, object$parM1, object$parM0)
    nlrCOV <- ifelse(1:m %in% object$DIFitems, object$covM1, object$covM0)
    nlrDM <- lapply(1:m, function(i) {
      .deltamethod.NLR.irt2log(
        par = nlrPAR[[i]], cov = nlrCOV[[i]],
        conv = !(i %in% object$conv.fail.which),
        cov_fail = is.null(nlrCOV[[i]])
      )
    })
    pars <- lapply(nlrDM, function(x) x$par)
    se <- lapply(nlrDM, function(x) x$se)
  }
  names(pars) <- nams
  names(se) <- nams

  if (SE) {
    coefs <- lapply(nams, function(i) rbind(pars[[i]], se[[i]]))
    coefs <- lapply(coefs, "rownames<-", c("estimate", "SE"))
    names(coefs) <- nams
  } else {
    coefs <- pars
  }

  if (CI > 0) {
    alpha <- (1 - CI) / 2
    CIlow <- lapply(nams, function(i) pars[[i]] - qnorm(1 - alpha) * se[[i]])
    CIupp <- lapply(nams, function(i) pars[[i]] + qnorm(1 - alpha) * se[[i]])
    names(CIlow) <- names(CIupp) <- nams
    coefs <- lapply(nams, function(i) rbind(coefs[[i]], CIlow[[i]], CIupp[[i]]))
    if (SE) {
      coefs <- lapply(coefs, "rownames<-", c("estimate", "SE", paste0("CI", 100 * alpha), paste0("CI", 100 * (1 - alpha))))
    } else {
      coefs <- lapply(coefs, "rownames<-", c("estimate", paste0("CI", 100 * alpha), paste0("CI", 100 * (1 - alpha))))
    }
    names(coefs) <- nams
  }

  if (simplify) {
    res <- as.data.frame(plyr::ldply(coefs, rbind))
    if (SE) {
      if (CI > 0) {
        resnams <- paste(res[, 1], c("estimate", "SE", paste0("CI", 100 * alpha), paste0("CI", 100 * (1 - alpha))))
      } else {
        resnams <- paste(res[, 1], c("estimate", "SE"))
      }
    } else {
      if (CI > 0) {
        resnams <- paste(res[, 1], c("estimate", paste0("CI", 100 * alpha), paste0("CI", 100 * (1 - alpha))))
      } else {
        resnams <- res[, 1]
      }
    }
    res <- res[, -1]
    rownames(res) <- resnams
    res[is.na(res)] <- 0
  } else {
    res <- coefs
  }

  return(res)
}

#' Log-likelihood and information criteria for an object of \code{"difNLR"} class.
#'
#' @aliases AIC.difNLR BIC.difNLR
#' @rdname logLik.difNLR
#'
#' @description S3 methods for extracting log-likelihood, Akaike's information criterion (AIC) and
#' Schwarz's Bayesian criterion (BIC) for an object of \code{"difNLR"} class.
#'
#' @param object an object of \code{"difNLR"} class.
#' @param item numeric or character: either character \code{"all"} to apply for all converged items (default),
#' or a vector of item names (column names of \code{Data}), or item identifiers (integers specifying
#' the column number).
#' @param ... other generic parameters for S3 methods.
#'
#' @author
#' Adela Hladka (nee Drabinova) \cr
#' Institute of Computer Science of the Czech Academy of Sciences \cr
#' Faculty of Mathematics and Physics, Charles University \cr
#' \email{hladka@@cs.cas.cz} \cr
#'
#' Patricia Martinkova \cr
#' Institute of Computer Science of the Czech Academy of Sciences \cr
#' \email{martinkova@@cs.cas.cz} \cr
#'
#' Karel Zvara \cr
#' Faculty of Mathematics and Physics, Charles University \cr
#'
#' @references
#' Drabinova, A. & Martinkova, P. (2017). Detection of differential item functioning with nonlinear regression:
#' A non-IRT approach accounting for guessing. Journal of Educational Measurement, 54(4), 498--517,
#' \doi{10.1111/jedm.12158}.
#'
#' Hladka, A. & Martinkova, P. (2020). difNLR: Generalized logistic regression models for DIF and DDF detection.
#' The R Journal, 12(1), 300--323, \doi{10.32614/RJ-2020-014}.
#'
#' Swaminathan, H. & Rogers, H. J. (1990). Detecting differential item functioning using logistic regression procedures.
#' Journal of Educational Measurement, 27(4), 361--370, \doi{10.1111/j.1745-3984.1990.tb00754.x}
#'
#' @seealso
#' \code{\link[difNLR]{difNLR}} for DIF detection among binary data using generalized logistic regression model. \cr
#' \code{\link[stats]{logLik}} for generic function extracting log-likelihood. \cr
#' \code{\link[stats]{AIC}} for generic function calculating AIC and BIC.
#'
#' @examples
#' \dontrun{
#' # loading data
#' data(GMAT)
#' Data <- GMAT[, 1:20] # items
#' group <- GMAT[, "group"] # group membership variable
#'
#' # testing both DIF effects using likelihood-ratio test and
#' # 3PL model with fixed guessing for groups
#' (x <- difNLR(Data, group, focal.name = 1, model = "3PLcg"))
#'
#' # AIC, BIC, log-likelihood
#' AIC(x)
#' BIC(x)
#' logLik(x)
#'
#' # AIC, BIC, log-likelihood for the first item
#' AIC(x, item = 1)
#' BIC(x, item = 1)
#' logLik(x, item = 1)
#' }
#' @export
logLik.difNLR <- function(object, item = "all", ...) {
  m <- length(object$nlrPAR)
  nams <- colnames(object$Data)

  if (inherits(item, "character")) {
    if (any(item != "all") & !all(item %in% nams)) {
      stop("Invalid value for 'item'. Item must be either character 'all', or
           numeric vector corresponding to column identifiers, or name of the item.",
        call. = FALSE
      )
    }
    if (any(item == "all")) {
      items <- 1:m
    } else {
      items <- which(nams %in% item)
    }
  } else {
    if (!inherits(item, "integer") & !inherits(item, "numeric")) {
      stop("Invalid value for 'item'. Item must be either character 'all', or
           numeric vector corresponding to column identifiers, or name of the item.",
        call. = FALSE
      )
    } else {
      if (!all(item %in% 1:m)) {
        stop("Invalid number for 'item'.",
          call. = FALSE
        )
      } else {
        items <- item
      }
    }
  }

  if (any(object$conv.fail.which %in% items)) {
    if (length(setdiff(items, object$conv.fail.which)) == 0) {
      stop(paste0("Item ", intersect(object$conv.fail.which, items), " does not converge. ",
        collapse = "\n"
      ),
      call. = FALSE
      )
    } else {
      warning(paste0("Item ", intersect(object$conv.fail.which, items), " does not converge. NA produced.",
        collapse = "\n"
      ),
      call. = FALSE
      )
      ITEMS <- setdiff(items, object$conv.fail.which)
    }
  } else {
    ITEMS <- items
  }

  val <- df <- rep(NA, m)

  val[ITEMS] <- ifelse(ITEMS %in% object$DIFitems,
    object$llM1[ITEMS],
    object$llM0[ITEMS]
  )
  df[ITEMS] <- ifelse(ITEMS %in% object$DIFitems,
    sapply(object$parM1, length)[ITEMS],
    sapply(object$parM0, length)[ITEMS]
  )
  val <- val[items]
  df <- df[items]
  if (length(items) == 1) {
    attr(val, "df") <- df
    class(val) <- "logLik"
  }
  return(val)
}

#' @rdname logLik.difNLR
#' @aliases BIC.difNLR logLik.difNLR
#' @export
AIC.difNLR <- function(object, item = "all", ...) {
  m <- length(object$nlrPAR)
  nams <- colnames(object$Data)

  if (inherits(item, "character")) {
    if (any(item != "all") & !all(item %in% nams)) {
      stop("Invalid value for 'item'. Item must be either character 'all', or
           numeric vector corresponding to column identifiers, or name of the item.",
        call. = FALSE
      )
    }
    if (any(item == "all")) {
      items <- 1:m
    } else {
      items <- which(nams %in% item)
    }
  } else {
    if (!inherits(item, "integer") & !inherits(item, "numeric")) {
      stop("Invalid value for 'item'. Item must be either character 'all', or
           numeric vector corresponding to column identifiers, or name of the item.",
        call. = FALSE
      )
    } else {
      if (!all(item %in% 1:m)) {
        stop("Invalid number for 'item'.",
          call. = FALSE
        )
      } else {
        items <- item
      }
    }
  }

  if (any(object$conv.fail.which %in% items)) {
    if (length(setdiff(items, object$conv.fail.which)) == 0) {
      stop(paste0("Item ", intersect(object$conv.fail.which, items), " does not converge. ",
        collapse = "\n"
      ),
      call. = FALSE
      )
    } else {
      warning(paste0("Item ", intersect(object$conv.fail.which, items), " does not converge. NA produced.",
        collapse = "\n"
      ),
      call. = FALSE
      )
      ITEMS <- setdiff(items, object$conv.fail.which)
    }
  } else {
    ITEMS <- items
  }

  k <- AIC <- rep(NA, m)
  k[ITEMS] <- ifelse(ITEMS %in% object$DIFitems,
    sapply(object$parM1, length)[ITEMS],
    sapply(object$parM0, length)[ITEMS]
  )
  AIC[ITEMS] <- 2 * k[ITEMS] - 2 * logLik(object, item = ITEMS)
  AIC <- AIC[items]

  return(AIC)
}

#' @rdname logLik.difNLR
#' @aliases AIC.difNLR logLik.difNLR
#' @export
BIC.difNLR <- function(object, item = "all", ...) {
  m <- length(object$nlrPAR)
  nams <- colnames(object$Data)

  if (inherits(item, "character")) {
    if (any(item != "all") & !all(item %in% nams)) {
      stop("Invalid value for 'item'. Item must be either character 'all', or
           numeric vector corresponding to column identifiers, or name of the item.",
        call. = FALSE
      )
    }
    if (any(item == "all")) {
      items <- 1:m
    } else {
      items <- which(nams %in% item)
    }
  } else {
    if (!inherits(item, "integer") & !inherits(item, "numeric")) {
      stop("Invalid value for 'item'. Item must be either character 'all', or
           numeric vector corresponding to column identifiers, or name of the item.",
        call. = FALSE
      )
    } else {
      if (!all(item %in% 1:m)) {
        stop("Invalid number for 'item'.",
          call. = FALSE
        )
      } else {
        items <- item
      }
    }
  }

  if (any(object$conv.fail.which %in% items)) {
    if (length(setdiff(items, object$conv.fail.which)) == 0) {
      stop(paste0("Item ", intersect(object$conv.fail.which, items), " does not converge. ",
        collapse = "\n"
      ),
      call. = FALSE
      )
    } else {
      warning(paste0("Item ", intersect(object$conv.fail.which, items), " does not converge. NA produced.",
        collapse = "\n"
      ),
      call. = FALSE
      )
      ITEMS <- setdiff(items, object$conv.fail.which)
    }
  } else {
    ITEMS <- items
  }
  k <- BIC <- rep(NA, m)
  k[ITEMS] <- ifelse(ITEMS %in% object$DIFitems,
    sapply(object$parM1, length)[ITEMS],
    sapply(object$parM0, length)[ITEMS]
  )
  n <- dim(object$Data)[1]
  BIC[ITEMS] <- log(n) * k[ITEMS] - 2 * logLik(object, item = ITEMS)
  BIC <- BIC[items]

  return(BIC)
}

#' Fitted values and residuals for an object of \code{"difNLR"} class.
#'
#' @aliases resid.difNLR fitted.difNLR
#' @rdname residuals.difNLR
#'
#' @description S3 methods for extracting fitted values and residuals for an object of \code{"difNLR"} class.
#'
#' @param object an object of \code{"difNLR"} class.
#' @param item numeric or character: either character \code{"all"} to apply for all converged items (default),
#' or a vector of item names (column names of \code{Data}), or item identifiers (integers specifying
#' the column number).
#' @param ... other generic parameters for S3 methods.
#'
#' @author
#' Adela Hladka (nee Drabinova) \cr
#' Institute of Computer Science of the Czech Academy of Sciences \cr
#' Faculty of Mathematics and Physics, Charles University \cr
#' \email{hladka@@cs.cas.cz} \cr
#'
#' Patricia Martinkova \cr
#' Institute of Computer Science of the Czech Academy of Sciences \cr
#' \email{martinkova@@cs.cas.cz} \cr
#'
#' Karel Zvara \cr
#' Faculty of Mathematics and Physics, Charles University \cr
#'
#' @references
#' Drabinova, A. & Martinkova, P. (2017). Detection of differential item functioning with nonlinear regression:
#' A non-IRT approach accounting for guessing. Journal of Educational Measurement, 54(4), 498--517,
#' \doi{10.1111/jedm.12158}.
#'
#' Hladka, A. & Martinkova, P. (2020). difNLR: Generalized logistic regression models for DIF and DDF detection.
#' The R Journal, 12(1), 300--323, \doi{10.32614/RJ-2020-014}.
#'
#' Swaminathan, H. & Rogers, H. J. (1990). Detecting differential item functioning using logistic regression procedures.
#' Journal of Educational Measurement, 27(4), 361--370, \doi{10.1111/j.1745-3984.1990.tb00754.x}
#'
#' @seealso
#' \code{\link[difNLR]{difNLR}} for DIF detection among binary data using generalized logistic regression model. \cr
#' \code{\link[stats]{fitted}} for generic function extracting fitted values. \cr
#' \code{\link[stats]{residuals}} for generic function extracting residuals.
#'
#' @examples
#' \dontrun{
#' # loading data
#' data(GMAT)
#' Data <- GMAT[, 1:20] # items
#' group <- GMAT[, "group"] # group membership variable
#'
#' # testing both DIF effects using likelihood-ratio test and
#' # 3PL model with fixed guessing for groups
#' (x <- difNLR(Data, group, focal.name = 1, model = "3PLcg"))
#'
#' # fitted values
#' fitted(x)
#' fitted(x, item = 1)
#'
#' # residuals
#' residuals(x)
#' residuals(x, item = 1)
#' }
#' @export
residuals.difNLR <- function(object, item = "all", ...) {
  # checking input
  n <- dim(object$Data)[1]

  m <- length(object$nlrPAR)
  nams <- colnames(object$Data)

  if (inherits(item, "character")) {
    if (any(item != "all") & !all(item %in% nams)) {
      stop("Invalid value for 'item'. Item must be either character 'all', or
           numeric vector corresponding to column identifiers, or name of the item.",
        call. = FALSE
      )
    }
    if (any(item == "all")) {
      items <- 1:m
    } else {
      items <- which(nams %in% item)
    }
  } else {
    if (!inherits(item, "integer") & !inherits(item, "numeric")) {
      stop("Invalid value for 'item'. Item must be either character 'all', or
           numeric vector corresponding to column identifiers, or name of the item.",
        call. = FALSE
      )
    } else {
      if (!all(item %in% 1:m)) {
        stop("Invalid number for 'item'.",
          call. = FALSE
        )
      } else {
        items <- item
      }
    }
  }

  if (any(object$conv.fail.which %in% items)) {
    if (length(setdiff(items, object$conv.fail.which)) == 0) {
      stop(paste("Item", intersect(object$conv.fail.which, items), "does not converge.",
        collapse = "\n"
      ),
      call. = FALSE
      )
    } else {
      warning(paste("Item", intersect(object$conv.fail.which, items), "does not converge. NA produced.",
        collapse = "\n"
      ),
      call. = FALSE
      )
      ITEMS <- setdiff(items, object$conv.fail.which)
    }
  } else {
    ITEMS <- items
  }

  residuals <- matrix(NA, nrow = n, ncol = m)
  residuals[, ITEMS] <- as.matrix(object$Data[, ITEMS] - fitted(object, item = ITEMS))

  residuals <- residuals[, items]

  return(residuals)
}

Try the difNLR package in your browser

Any scripts or data that you put into this service are public.

difNLR documentation built on May 3, 2023, 5:11 p.m.