R/person.posterior.R

In GPCMlasso: Regularized Explanatory Generalized Partial Credit Models

## From a vector of adjacent-category linear predictors, compute category
## probabilities for the adjacent categories model.
responseFun <- function(eta) {
  q <- length(eta)
  
  eta[eta > 10] <- 10
  eta[eta < -10] <- -10
  
  eta.help <- matrix(
    rep(c(0, eta), each = q + 1),
    ncol = q + 1
  )
  eta.help[upper.tri(eta.help)] <- 0
  
  pi <- cumprod(c(1, exp(eta[-q]))) /
    sum(apply(exp(eta.help), 1, prod))
  
  pi <- (pi - 0.5) * 0.99999 + 0.5
  
  pi
}



## Log likelihood contribution of one person, given one theta value.
##
## lin.pred already contains all fixed effects in the correct parametrization.
## The theta contribution is added as sigma_i * theta.
##
## y_offset controls response coding:
##   y_offset = 0: responses are coded as 1, ..., K
##   y_offset = 1: responses are coded as 0, ..., k
loglik_i_theta <- function(theta,
                           lin.pred,
                           sigma,
                           yp,
                           q,
                           y_offset = 0L) {
  
  eta <- lin.pred + rep(sigma, q) * theta
  eta_split <- split(eta, rep(seq_along(yp), q))
  
  loglik <- 0
  
  for (i in seq_along(eta_split)) {
    mu_i <- responseFun(eta_split[[i]])
    mu_i <- c(mu_i, 1 - sum(mu_i))
    
    category_index <- as.integer(yp[i]) + y_offset
    
    if (
      is.na(category_index) ||
      category_index < 1 ||
      category_index > length(mu_i)
    ) {
      return(-Inf)
    }
    
    prob_i <- mu_i[category_index]
    
    if (is.na(prob_i) || prob_i <= 0) {
      return(-Inf)
    }
    
    loglik <- loglik + log(prob_i)
  }
  
  loglik
}



## Create vector of fixed linear predictors for all observations.
##
## If scale_cols = numeric(0), old behavior is used:
##   lin.pred = sigma_i * (des_total %*% coef_short)
##
## If scale_cols is a 0/1 vector:
##   scale_cols == 1: corresponding design column is multiplied by sigma_i
##   scale_cols == 0: corresponding design column is not multiplied by sigma_i
lin_pred <- function(model, coef_short, sigma) {
  
  with(model$design_list, {
    
    des_total <- matrix(
      rep(t(design), n),
      byrow = TRUE,
      ncol = ncol(design)
    )
    
    if (ncol(designX) > 0) {
      des_total <- cbind(des_total, designX)
    }
    
    scale_cols <- model$scale_cols
    
    if (is.null(scale_cols)) {
      scale_cols <- numeric(0)
    }
    
    if (length(scale_cols) == 0) {
      
      des_total_scaled <- rep(rep(sigma, q), n) * des_total
      lin.pred <- des_total_scaled %*% coef_short
      
    } else {
      
      if (length(scale_cols) != ncol(des_total)) {
        stop(
          paste0(
            "scale_cols must have length equal to ncol(des_total). ",
            "length(scale_cols) = ", length(scale_cols),
            ", ncol(des_total) = ", ncol(des_total), "."
          ),
          call. = FALSE
        )
      }
      
      scale_cols <- as.numeric(scale_cols)
      
      if (any(!scale_cols %in% c(0, 1))) {
        stop("scale_cols must only contain 0 and 1.", call. = FALSE)
      }
      
      beta_long <- rep(rep(sigma, q), n)
      
      idx_scaled <- which(scale_cols == 1)
      idx_unscaled <- which(scale_cols == 0)
      
      lin.pred.scaled <- rep(0, nrow(des_total))
      lin.pred.unscaled <- rep(0, nrow(des_total))
      
      if (length(idx_scaled) > 0) {
        lin.pred.scaled <- as.vector(
          des_total[, idx_scaled, drop = FALSE] %*%
            coef_short[idx_scaled]
        )
      }
      
      if (length(idx_unscaled) > 0) {
        lin.pred.unscaled <- as.vector(
          des_total[, idx_unscaled, drop = FALSE] %*%
            coef_short[idx_unscaled]
        )
      }
      
      lin.pred <- beta_long * lin.pred.scaled + lin.pred.unscaled
    }
    
    as.vector(lin.pred)
  })
}



## Estimate one person's posterior mean using Gauss-Hermite quadrature.
person_i_gh <- function(person,
                        person.index,
                        all.lin.preds,
                        sigma,
                        Y,
                        q,
                        GHnodes,
                        GHweights,
                        y_offset = 0L) {
  
  yp_all <- Y[person, ]
  obs_items <- !is.na(yp_all)
  
  if (!any(obs_items)) {
    return(NA_real_)
  }
  
  yp <- yp_all[obs_items]
  sigma_obs <- sigma[obs_items]
  q_obs <- q[obs_items]
  
  lin.pred_all <- all.lin.preds[person.index == person]
  
  item_index_long <- rep(seq_along(q), q)
  keep_long <- item_index_long %in% which(obs_items)
  lin.pred <- lin.pred_all[keep_long]
  
  log_lik <- vapply(
    GHnodes,
    loglik_i_theta,
    numeric(1),
    lin.pred = lin.pred,
    sigma = sigma_obs,
    yp = yp,
    q = q_obs,
    y_offset = y_offset
  )
  
  log_w <- log(GHweights) + log_lik
  
  max_log_w <- max(log_w, na.rm = TRUE)
  
  if (!is.finite(max_log_w)) {
    return(NA_real_)
  }
  
  w <- exp(log_w - max_log_w)
  denom <- sum(w)
  
  if (!is.finite(denom) || denom <= 0) {
    return(NA_real_)
  }
  
  sum(GHnodes * w) / denom
}



#' Calculate Posterior Estimates for Trait Parameters
#'
#' Calculates posterior estimates for trait/person parameters for a fitted
#' \code{GPCMlasso} model using the assumed Gaussian distribution of the
#' person parameters.
#'
#' The function computes posterior means of the latent trait parameters by
#' Gauss-Hermite quadrature. If no coefficient vector is supplied, the
#' cross-validation optimal coefficient vector is used when cross-validation
#' was performed; otherwise, the BIC-optimal coefficient vector is used.
#'
#' @param model Object of class \code{GPCMlasso}.
#' @param coefs Optional vector of coefficients. If \code{coefs = NULL}, the
#' coefficients from the BIC-optimal model are used, or, if cross-validation
#' was performed, the coefficients from the cross-validation optimal model are
#' used.
#' @param cores Number of cores used for parallel computation.
#' @param tol Deprecated. Kept for backward compatibility.
#'
#' @return Numeric vector containing posterior estimates of the trait/person
#' parameters.
#'
#' @author Gunther Schauberger\cr \email{gunther.schauberger@@tum.de}
#'
#' @seealso
#' \code{\link{GPCMlasso}},
#' \code{\link{predict.GPCMlasso}}
#'
#' @examples
#' data(tenseness_small)
#'
#' form0 <- as.formula(
#'   paste(
#'     "cbind(",
#'     paste(colnames(tenseness_small)[1:5], collapse = ","),
#'     ") ~ 0"
#'   )
#' )
#'
#' \dontrun{
#' rsm0 <- GPCMlasso(
#'   formula = form0,
#'   data = tenseness_small,
#'   model = "RSM",
#'   control = ctrl_GPCMlasso(cores = 1, trace = FALSE)
#' )
#'
#' theta_hat <- trait.posterior(rsm0, cores = 1)
#' summary(theta_hat)
#' }
#'
#' @export
trait.posterior <- function(model, coefs = NULL, cores = 25, tol = 1e-4) {
  
  n <- model$design_list$n
  I <- model$design_list$I
  q <- model$design_list$q
  n_sigma <- model$design_list$n_sigma
  
  if (is.null(coefs)) {
    if (!is.null(model$cv_error)) {
      cat("No coefs are provided, automatically cv-optimal model is chosen", "\n")
      coefs <- model$coefficients[which.min(model$cv_error), ]
    } else {
      cat("No coefs are provided, automatically BIC-optimal model is chosen", "\n")
      coefs <- model$coefficients[which.min(model$BIC), ]
    }
  }
  
  coefs <- as.numeric(coefs)
  
  coef_short <- head(coefs, length(coefs) - n_sigma)
  sigma <- tail(coefs, n_sigma)
  
  if (n_sigma == 1) {
    sigma <- rep(sigma, I)
  }
  
  if (length(sigma) != I) {
    stop(
      paste0(
        "Number of discrimination parameters does not match number of items. ",
        "Expected ", I, " but got ", length(sigma), "."
      ),
      call. = FALSE
    )
  }
  
  person.index <- rep(seq_len(n), each = sum(q))
  
  Y <- matrix(
    as.numeric(as.matrix(model$Y)),
    ncol = ncol(model$Y)
  )
  
  ## Determine response coding.
  ## If at least one observed response is 0, category 0 maps to probability
  ## index 1 in R.
  if (any(Y == 0, na.rm = TRUE)) {
    y_offset <- 1L
  } else {
    y_offset <- 0L
  }
  
  all.lin.preds <- lin_pred(
    model = model,
    coef_short = coef_short,
    sigma = sigma
  )
  
  Q <- NULL
  
  if (!is.null(model$control$Q)) {
    Q <- model$control$Q
  }
  
  if (is.null(Q) && !is.null(model$design_list$Q)) {
    Q <- model$design_list$Q
  }
  
  if (is.null(Q)) {
    Q <- 21
  }
  
  her_poly <- gauss.quad(Q, "hermite")
  GHnodes <- her_poly$nodes
  GHweights <- her_poly$weights * exp(GHnodes^2) * dnorm(GHnodes)
  
  estimates <- person.fit(
    n = n,
    q = q,
    Y = Y,
    sigma = sigma,
    person.index = person.index,
    all.lin.preds = all.lin.preds,
    GHnodes = GHnodes,
    GHweights = GHweights,
    y_offset = y_offset,
    cores = cores
  )
  
  names(estimates) <- rownames(model$data)
  
  estimates
}



person.fit <- function(n,
                       q,
                       Y,
                       sigma,
                       person.index,
                       all.lin.preds,
                       GHnodes,
                       GHweights,
                       y_offset = 0L,
                       cores = 1) {
  
  if (cores > 1) {
    
    cl <- parallel::makeCluster(cores, outfile = "")
    
    on.exit(
      parallel::stopCluster(cl),
      add = TRUE
    )
    
    parallel::clusterExport(
      cl,
      varlist = c(
        "Y",
        "person.index",
        "all.lin.preds",
        "sigma",
        "q",
        "GHnodes",
        "GHweights",
        "y_offset",
        "loglik_i_theta",
        "responseFun",
        "person_i_gh"
      ),
      envir = environment()
    )
    
    estimates <- parallel::parSapply(
      cl,
      seq_len(n),
      person_i_gh,
      person.index = person.index,
      all.lin.preds = all.lin.preds,
      sigma = sigma,
      Y = Y,
      q = q,
      GHnodes = GHnodes,
      GHweights = GHweights,
      y_offset = y_offset
    )
    
  } else {
    
    estimates <- sapply(
      seq_len(n),
      person_i_gh,
      person.index = person.index,
      all.lin.preds = all.lin.preds,
      sigma = sigma,
      Y = Y,
      q = q,
      GHnodes = GHnodes,
      GHweights = GHweights,
      y_offset = y_offset
    )
  }
  
  estimates
}