R/Logistik.R

In difR: Collection of Methods to Detect Dichotomous and Polytomous Differential Item Functioning (DIF)

#' @importFrom stats glm deviance vcov coefficients
#' @export
Logistik <- function(data, member, member.type = "group", match = "score", 
    anchor = 1:ncol(data), type = "both", criterion = "LRT", all.cov = FALSE) 
{
#' @keywords internal
#' @noRd
    R2 <- function(m, n) 1 - (exp(-m$null.deviance/2 + m$deviance/2))^(2/n)
#' @keywords internal
#' @noRd
    R2max <- function(m, n) 1 - (exp(-m$null.deviance/2))^(2/n)
#' @keywords internal
#' @noRd
    R2DIF <- function(m, n) R2(m, n)/R2max(m, n)

    dev <- R2full <- R2simple <- deltaR <- NULL
    mFull <- mSimple <- seFull <- seSimple <- matrix(0, ncol(data), 4)

    cov.matM0 <- cov.matM1 <- if (all.cov) vector("list", ncol(data)) else NULL

    GROUP <- if (member.type == "group") as.factor(member) else member

    for (item in 1:ncol(data)) {
        if (match[1] == "score") {
            data2 <- data[, anchor, drop = FALSE]
            if (!(item %in% anchor)) {
                data2 <- cbind(data2, data[, item])
            }
            SCORES <- rowSums(data2, na.rm = TRUE)
        } else if (match[1] == "restscore") {
            rest_anchor <- setdiff(anchor, item)
            if (length(rest_anchor) == 0) {
                SCORES <- rep(NA, nrow(data))
            } else {
                SCORES <- rowSums(data[, rest_anchor, drop = FALSE], na.rm = TRUE)
            }
        } else {
            SCORES <- match
        }

        ITEM <- data[, item]
        m0 <- switch(type,
                     both = glm(ITEM ~ SCORES * GROUP, family = "binomial"),
                     udif = glm(ITEM ~ SCORES + GROUP, family = "binomial"),
                     nudif = glm(ITEM ~ SCORES * GROUP, family = "binomial"))
        m1 <- switch(type,
                     both = glm(ITEM ~ SCORES, family = "binomial"),
                     udif = glm(ITEM ~ SCORES, family = "binomial"),
                     nudif = glm(ITEM ~ SCORES + GROUP, family = "binomial"))

        if (criterion == "LRT") {
            dev[item] <- deviance(m1) - deviance(m0)
        } else {
            if (criterion != "Wald") stop("'criterion' must be either 'LRT' or Wald'", call. = FALSE)
            coeff <- as.numeric(coefficients(m0))
            covMat <- summary(m0)$cov.scaled
            C <- switch(type,
                        udif = rbind(c(0, 0, 1)),
                        nudif = rbind(c(0, 0, 0, 1)),
                        both = rbind(c(0, 0, 1, 0), c(0, 0, 0, 1)))
            dev[item] <- t(C %*% coeff) %*% solve(C %*% covMat %*% t(C)) %*% C %*% coeff
        }

        R2full[item] <- R2DIF(m0, nrow(data))
        R2simple[item] <- R2DIF(m1, nrow(data))
        deltaR[item] <- R2DIF(m0, nrow(data)) - R2DIF(m1, nrow(data))

        mFull[item, 1:length(m0$coefficients)] <- m0$coefficients
        mSimple[item, 1:length(m1$coefficients)] <- m1$coefficients
        seFull[item, 1:length(m0$coefficients)] <- sqrt(diag(vcov(m0)))
        seSimple[item, 1:length(m1$coefficients)] <- sqrt(diag(vcov(m1)))

        if (all.cov) cov.matM0[[item]] <- vcov(m0)
        if (all.cov) cov.matM1[[item]] <- vcov(m1)
    }

    colnames(mFull) <- colnames(mSimple) <- colnames(seFull) <- colnames(seSimple) <-
        c("(Intercept)", "SCORE", "GROUP", "SCORE:GROUP")

    res <- list(stat = dev, R2M0 = R2full, R2M1 = R2simple, deltaR2 = deltaR,
                parM0 = mFull, parM1 = mSimple, seM0 = seFull, seM1 = seSimple,
                cov.M0 = cov.matM0, cov.M1 = cov.matM1,
                criterion = criterion, member.type = member.type,
                match = ifelse(match[1] == "score", "score", 
                               ifelse(match[1] == "restscore", "restscore", "matching variable")))
    return(res)
}