R/irtc_mml_2pl_mstep_slope.R

In IRTC: Marginal Maximum Likelihood Estimation for Item Response Models

# IRTC
# Copyright (C) 2026 WEIAN DATA TECH (Beijing) Co., Ltd.
# SPDX-License-Identifier: GPL-2.0-or-later
# See inst/COPYRIGHTS for licensing details.

## File Name: irtc_mml_2pl_mstep_slope.R

irtc_mml_2pl_mstep_slope <- function (B_orig, B, B_obs, B.fixed, max.increment, nitems, A,
            AXsi, xsi, theta, nnodes, maxK, itemwt, Msteps, ndim, convM,
            irtmodel, progress, est.slopegroups, E, basispar, se.B,
            equal.categ, B_acceleration, trim_increment="half", iter, eps=1E-10,
            maxcat=NULL, use_rcpp=FALSE, use_rcpp_calc_prob=TRUE, n_threads=1L )
{
    old_increment <- array(max.increment, dim=c(nitems, maxK, ndim))
    xbar2 <- xxf <- xbar <- array(0, dim=c(nitems, maxK))
    xtemp <- matrix(0, nrow=1, ncol=1)
    converge <- FALSE
    Biter <- 1
    mK <- 1L:maxK
    items.temp <- 1L:nitems
    items.conv <- NULL
    items.conv1 <- c(-1)
    increments_msteps <- rep(NA, Msteps)
    B_old <- B
    eps <- 1E-10
    NI <- nitems
    TP <- nrow(theta)
    if (progress) {
        cat("\nM Step Slopes       |")
        utils::flush.console()
    }

    if (irtmodel %in% c("GPCM", "GPCM.groups")) {
        old_increment.temp <- matrix(.3, nrow=NI, ncol=ndim)
    }
    if (irtmodel == "GPCM.design") {
        Nlambda <- ncol(E)
        old_increment.temp <- matrix(.3, nrow=Nlambda, ncol=ndim)
    }
    if (irtmodel == "2PL.groups") {
        ES <- length(unique(est.slopegroups))
        old_increment.temp <- array(.3, dim=c(ES, maxK, ndim))
        use_rcpp <- FALSE
    }

    while (!converge & Biter <= Msteps) {
        probability_result <- irtc_mml_calc_prob(
            iIndex=items.temp, A=A, AXsi=AXsi, B=B, xsi=xsi,
            theta=theta, nnodes=nnodes, maxK=maxK, maxcat=maxcat,
            use_rcpp=use_rcpp_calc_prob
        )
        rprobs <- probability_result$rprobs
        LIT <- length(items.temp)

        for (dd in 1L:ndim) {
            if (irtmodel %in% c("GPCM", "GPCM.design", "GPCM.groups")) {
                xtemp <- matrix(0, nrow=LIT, ncol=TP)
            }
            if (irtmodel == "2PL.groups") {
                xtemp <- array(0, dim=c(LIT, TP, maxK))
            }

            if (!equal.categ) {
                rprobs[is.na(rprobs)] <- 0
            }

            if (use_rcpp) {
                sufficient_statistics <- irtc_rcpp_mml_2pl_mstep_item_slopes_suffstat(
                    rprobs=as.vector(rprobs), items_temp=items.temp,
                    theta=theta, dd=dd-1, LIT=LIT, TP=TP,
                    nitems=nitems, maxcat=maxcat, maxK=maxK,
                    itemwt=itemwt, xxf_=xxf, xbar_=xbar,
                    xbar2_=xbar2, irtmodel=irtmodel, xtemp_=xtemp,
                    items_conv=items.conv1, n_threads=n_threads
                )
            } else {
                sufficient_statistics <- irtc_mml_2pl_mstep_item_slopes_suffstat_R(
                    rprobs=rprobs, maxK=maxK, LIT=LIT, TP=TP,
                    itemwt=itemwt, theta=theta, dd=dd,
                    items.temp=items.temp, items.conv=items.conv,
                    xbar=xbar, xbar2=xbar2, xxf=xxf,
                    xtemp=xtemp, irtmodel=irtmodel
                )
            }
            xbar <- sufficient_statistics$xbar
            xbar2 <- sufficient_statistics$xbar2
            xxf <- sufficient_statistics$xxf
            xtemp <- sufficient_statistics$xtemp

            if (irtmodel %in% c("GPCM", "GPCM.design", "GPCM.groups")) {
                B_obs.temp <- B_obs
                category_scores <- irtc_matrix2(mK - 1, LIT, maxK)
                B_obs.temp[items.temp, , dd] <- (
                    category_scores * B_obs[items.temp, , dd]
                )
                xbar[is.na(xbar)] <- 0
                xbar.temp <- category_scores * xbar
                diff.temp <- rowSums(B_obs.temp[, , dd]) - rowSums(xbar.temp)
                xbar2.temp <- diag(xtemp^2 %*% itemwt[, items.temp])
                xxf.temp <- rowSums(xxf)
                if (!is.null(items.conv)) {
                    diff.temp[items.conv, ] <- 0
                }
                deriv.temp <- xbar2.temp - xxf.temp
            }

            if (irtmodel %in% c("GPCM", "GPCM.groups")) {
                deriv.temp[is.na(deriv.temp)] <- eps

                if (irtmodel == "GPCM.groups") {
                    diff.temp <- irtc_aggregate_derivative_information(
                        deriv=diff.temp, groups=est.slopegroups
                    )
                    deriv.temp <- irtc_aggregate_derivative_information(
                        deriv=deriv.temp, groups=est.slopegroups
                    )
                }

                increment.temp <- diff.temp * abs(1 / (deriv.temp + eps))
                increment.temp <- irtc_trim_increment(
                    increment=increment.temp,
                    max.increment=abs(old_increment.temp[, dd]),
                    trim_increment=trim_increment
                )
                old_increment.temp[, dd] <- increment.temp

                increment <- irtc_outer( increment.temp, mK - 1)
                if (Biter == 1) {
                    se.B[, , dd] <- irtc_outer(
                        sqrt(1 / abs(deriv.temp)), mK - 1
                    )
                }
            }

            if (irtmodel == "GPCM.design") {
                diff.temp <- (t(E) %*% diff.temp)[, 1]
                deriv.temp <- (t(E^2) %*% deriv.temp)[, 1]
                increment.temp <- diff.temp * abs(1 / (deriv.temp + eps))
                increment.temp <- irtc_trim_increment(
                    increment=increment.temp,
                    max.increment=abs(old_increment.temp[, dd]),
                    trim_increment=trim_increment
                )
                old_increment.temp[, dd] <- increment.temp
                basispar[, dd] <- basispar[, dd] + increment.temp
                increment.temp <- (E %*% increment.temp)[, 1]
                increment <- irtc_outer( increment.temp, mK-1)
                d1 <- irtc_outer(1 / abs(deriv.temp), mK - 1)
                LL <- ncol(d1)
                for (ll in 1L:LL) {
                    m1 <- sqrt(diag(
                        E %*% d1[, ll] %*% t(d1[, ll]) %*% t(E)
                    ))
                    if (Biter == 1) {
                        se.B[, ll, dd] <- m1
                    }
                }
                nB <- dim(B)
                B_ind <- 1 * (B_orig != 0)
                for (dd in 1L:nB[3]) {
                    EB <- E %*% basispar[, dd]
                    for (cc in 1L:(nB[2] - 1)) {
                        B[, cc+1, dd] <- cc * EB * B_ind[, cc+1, dd]
                    }
                }
                B00 <- B
            }

            if (irtmodel == "2PL.groups") {
                a1 <- stats::aggregate(
                    B_obs[, , dd] - xbar, list(est.slopegroups), sum
                )
                a2 <- stats::aggregate(
                    xbar2 - xxf, list(est.slopegroups), sum
                )
                deriv.temp <- as.matrix(a2[, -1])
                diff.temp <- as.matrix(a1[, -1])
                increment.temp <- diff.temp * abs(1 / (deriv.temp + eps))
                increment.temp <- irtc_trim_increment(
                    increment=increment.temp,
                    max.increment=abs(old_increment.temp[, , dd]),
                    trim_increment=trim_increment
                )
                old_increment.temp[, , dd] <- increment.temp
                ind.match <- match(est.slopegroups, a1[, 1])
                increment <- increment.temp[ind.match, ]
                if (Biter == 1) {
                    se.B[, , dd] <- sqrt(1 / abs(deriv.temp[ind.match, ]))
                }
            }

            if (irtmodel == "2PL") {
                diff <- B_obs[, , dd] - xbar
                if (!is.null(items.conv)) {
                    diff[items.conv, ] <- 0
                }
                deriv <- xbar2 - xxf
                increment <- diff * abs(1 / (deriv + eps))
                increment <- irtc_trim_increment(
                    increment=increment,
                    max.increment=abs(old_increment[, , dd]),
                    trim_increment=trim_increment
                )
                if (Biter == 1) {
                    se.B[, , dd] <- sqrt(1 / abs(deriv))
                }
            }

            increment[B_orig[, , dd] == 0] <- 0
            old_increment[, , dd] <- increment
            B[, , dd] <- B[, , dd] + increment
            if (irtmodel == "GPCM.design") {
                B <- B00
            }
        }

        if (irtmodel == "2PL") {
            items.temp <- which(apply(
                old_increment, 1,
                FUN=function(ll) { !(max(abs(ll)) < convM) }
            ))
            items.conv <- setdiff(1L:nitems, items.temp)
            if (length(items.conv) > 0) {
                items.conv1 <- items.conv
            }
        }

        max_change <- max(abs(old_increment))
        increments_msteps[Biter] <- max_change
        if (max_change < convM) {
            converge <- TRUE
        }
        Biter <- Biter + 1
        if (progress) {
            cat("-")
            utils::flush.console()
        }
    }

    se.B[B_orig == 0] <- 0

    if (B_acceleration$acceleration != "none") {
        B_acceleration <- irtc_accelerate_parameters(
            xsi_acceleration=B_acceleration, xsi=as.vector(B),
            iter=iter, itermin=3
        )
        B <- array(B_acceleration$parm, dim(B))
    }
    B_change <- max(abs(B - B_old))

    list(
        "B"=B, "basispar"=basispar, "se.B"=se.B, Biter=Biter,
        B_acceleration=B_acceleration, B_change=B_change,
        increments_msteps=increments_msteps
    )
}

Mstep_slope.v2 <- irtc_mml_2pl_mstep_slope