R/model_adequacy.R

In trajeR: Group Based Modeling Trajectory

#' @title Average Posterior Probability of Assignment
#'
#' @description The Average Posterior Probability (AvePP) measures the classification quality. For each group, it is the average of the highest posterior probabilities for all individuals assigned to that group.
#' A value close to 1 indicates a clear classification.
#'
#' @param sol A trajectory object returned by the \code{trajeR} function.
#' @param Y The response variable matrix, as used in the `trajeR` call.
#' @param A The time variable matrix, as used in the `trajeR` call.
#' @param X Optional. A matrix of covariates that modify the group membership probability, if they were used in the model.
#'
#' @return A named numeric vector of the same length as the number of groups, containing the AvePP for each group.
#' @export
#'
#' @examples
#' data <- read.csv(system.file("extdata", "CNORM2gr.csv", package = "trajeR"))
#' data <- as.matrix(data)
#' # Use drop=FALSE for single-column matrices to preserve dimensions
#' sol <- trajeR(Y = data[, 2:6], A = data[, 7:11], Risk = data[, 12, drop = FALSE],
#'               degre = c(2, 2), Model = "CNORM", Method = "EM")
#' AvePP(sol, Y = data[, 2:6], A = data[, 7:11])
AvePP <- function(sol, Y, A, X = NULL) {
  Xt <- X
  prob <- GroupProb(sol, Y, A, X = Xt)
  gr <- as.vector(sapply(1:nrow(prob), function(s) {
    which.max(prob[s, ])
  }))
  res <- c()
  for (i in 1:sol$groups) {
    res <- c(res, mean(apply(prob[gr == i, ], 1, max)))
  }
  names(res) <- paste("Group", 1:sol$groups)
  return(res)
}

#' @title Odds of Correct Classification (OCC)
#'
#' @description The Odds of Correct Classification (OCC) compares the odds of correct assignment using the model's posterior probabilities against the odds of assignment based on the estimated group proportions alone.
#' A high OCC value suggests that the model provides a much better classification than a random assignment.
#'
#' @inheritParams AvePP
#'
#' @return A named numeric vector of the same length as the number of groups, containing the OCC for each group.
#' @export
#'
#' @examples
#' data <- read.csv(system.file("extdata", "CNORM2gr.csv", package = "trajeR"))
#' data <- as.matrix(data)
#' sol <- trajeR(Y = data[, 2:6], A = data[, 7:11], Risk = data[, 12, drop = FALSE],
#'               degre = c(2, 2), Model = "CNORM", Method = "EM")
#' OCC(sol, Y = data[, 2:6], A = data[, 7:11])
OCC <- function(sol, Y, A) {
  tmp <- AvePP(sol, Y, A, X = NULL)
  if (sol$Method == "L" && (length(sol$theta) == sol$groups) > 1) {
    prob <- exp(sol$theta) / sum(exp(sol$theta))
  } else {
    prob <- sol$theta
  }
  return(tmp / (1 - tmp) / (prob / (1 - prob)))
}

#' @title Proportion of Individuals Assigned to Each Group
#' @description Calculate the proportion of individuals in a given group. That is the ratio of
#' the number of individuals in one group and all the individuals.
#' @inheritParams AvePP
#'
#' @return A named numeric vector containing the proportion of individuals assigned to each group.
#' @export
#'
#' @examples
#' data <- read.csv(system.file("extdata", "CNORM2gr.csv", package = "trajeR"))
#' data <- as.matrix(data)
#' sol <- trajeR(Y = data[, 2:6], A = data[, 7:11], Risk = data[, 12, drop = FALSE],
#'               degre = c(2, 2), Model = "CNORM", Method = "EM")
#' propAssign(sol, Y = data[, 2:6], A = data[, 7:11])
propAssign <- function(sol, Y, A) {
  prob <- GroupProb(sol, Y, A, X = NULL)
  gr <- as.vector(sapply(1:nrow(prob), function(s) {
    which.max(prob[s, ])
  }))
  tab <- matrix(
    sapply(1:sol$groups, function(s) {
      sum(gr == s)
    }),
    nrow = 1
  ) /
    sol$Size
  res <- as.vector(tab)
  names(res) <- paste("Group", 1:sol$groups)
  return(res)
}

#' @title Bootstrap Confidence Intervals for Group Membership Probabilities
#'
#' @description Calculates confidence intervals for the estimated group membership probabilities using a parametric bootstrap method.
#' It simulates new theta parameters from a normal distribution based on the estimated theta and their standard errors.
#'
#' @inheritParams AvePP
#' @param nb An integer. The number of repetitions in the bootstrap method.
#' @param alpha A number. The significance level for the confidence interval (e.g., 0.05 for a 95\% CI).
#' @return A matrix with two rows (lower and upper bounds) and a column for each group,
#'   containing the confidence intervals for the estimated probabilities.
#' @export
#'
#' @examples
#' data <- read.csv(system.file("extdata", "CNORM2gr.csv", package = "trajeR"))
#' sol <- trajeR(Y = data[, 2:6], A = data[, 7:11], degre = c(2, 2), Model = "CNORM", Method = "L")
#' ConfIntT(sol, Y = data[, 2:6], A = data[, 7:11])
ConfIntT <- function(sol, Y, A, nb = 10000, alpha = 0.98) {
  Xt <- cbind(matrix(rep(1, sol$Size), ncol = 1))
  vtmp <- c(sol$beta, sol$delta, sol$phi, sol$nu, sol$sigma)
  vtmp <- vtmp[!is.na(vtmp)]
  indmin <- length(vtmp) + 1
  indmax <- indmin + length(c(sol$theta)) - 1
  theta <- sol$tab[indmin:indmax, 1]
  sdtheta <- sol$tab[indmin:indmax, 2]
  boottheta <- sapply(1:sol$groups, function(s) {
    stats::rnorm(nb, theta[s], sdtheta[s])
  })
  prob <- boottheta
  if (sol$Method == "L" && (length(sol$theta) == sol$groups) > 1) {
    prob <- exp(boottheta) / rowSums(exp(boottheta))
  }
  sapply(1:sol$groups, function(s) {
    stats::quantile(prob[, s], probs = c((1 - alpha) / 2, 1 - (1 - alpha) / 2))
  })
}

#' @title Model Adequacy Assessment
#'
#' @description Provides a summary table of model adequacy metrics, including estimated group probabilities, confidence intervals, assignment proportions, Average Posterior Probability (AvePP), and Odds of Correct Classification (OCC).
#'
#' @inheritParams AvePP
#' @param nb Integer. The number of repetitions for the bootstrap confidence interval.
#' @param alpha  Real. The degree of confidence of the interval.
#' @return A table of reals. A table with 5 rows: the estimate probabilities, the
#'  two bounds of the confidence interval, the proportion of assignment, the
#'  Average Posterior Probability and the Odds of Correct Classification.
#' @export
#'
#' @examples
#' data <- read.csv(system.file("extdata", "CNORM2gr.csv", package = "trajeR"))
#' data <- as.matrix(data)
#' sol <- trajeR(Y = data[, 2:6], A = data[, 7:11], degre = c(2, 2), Model = "CNORM", Method = "L")
#' adequacy(sol, Y = data[, 2:6], A = data[, 7:11])
adequacy <- function(sol, Y, A, nb = 10000, alpha = 0.98) {
  if (sol$Method == "L" && (length(sol$theta) == sol$groups) > 1) {
    prob_est <- exp(sol$theta) / sum(exp(sol$theta))
  } else {
    prob_est <- sol$theta
  }

  tab <- rbind(
    prob_est,
    ConfIntT(sol, Y, A, nb = 10000, alpha = 0.98),
    propAssign(sol, Y, A),
    AvePP(sol, Y, A),
    OCC(sol, Y, A)
  )
  rownames(tab) <- c(
    "Prob. est.",
    "CI inf.",
    "CI sup.",
    "Prop.",
    "AvePP",
    "OCC"
  )
  return(tab)
}