R/prediction.R

In jpasquier/ad_test: What the Package Does (One Line, Title Case)

#' Compute prediction
#'
#' Compute a prediction with ...
#'
#' @param time individual (and time dependent) times 
#' @param comp individual and time dependent compliances (0/1) 
#' @param id (individual and) time dependent ids 
#' @param w weights of the weighted adherence method 
#' @param cor correlation structure of the GEE model
#' @param k knots for the splines method
#' @param V variance estimate of model coefficients (default="naive")
#'
#' @return a list containing
#' \itemize{
#'   \item predict: dataframe containing prediction and confidence intervals of
#'     the model
#'   \item q: QIC goodness of fit statistics
#' }
#'
#' @examples
#' data(onco)
#' onco$compliance <- as.numeric(onco$observed >= onco$expected)
#' prediction(onco$drel, onco$compliance, onco$id, "ar1",
#'            c(120, 240), rep(1, nrow(onco)))
#'
#' @importFrom splines ns
#' @importFrom geeM geem
#' @importFrom MuMIn QICu
#'
#' @export

prediction <- function(time, comp, id, cor, k, w, V = "naive") {
   
  # To call the QICu function, the variables used in geem have to be in the
  # global environment
  .time <<- time
  .comp <<- comp
  .id <<- id
  .cor <<- cor
  .k <<- k
  .w <<- w
  out <- geeM::geem(.comp ~ splines::ns(.time, knots = .k), id = .id,
                    family = binomial, corstr = .cor, weights = .w)
  q <- MuMIn::QICu(out, typeR = TRUE)
  rm(envir = .GlobalEnv, .time, .comp, .id, .cor, .k, .w)

  t <- unique(time)
  X <- cbind(1, splines::ns(1:length(t), knots = k))
  mu.hat <- X %*% out$beta
  pr <- plogis(mu.hat)

  C <- out$naiv.var
  if (V == "robust") {
    C <- out$var
  }

  if (dim(C)[1] == 1) {
    Var <- as.vector(C)
  } else {
    dia <- matrix(NA, nrow(X), ncol(X))
    for (i in 1:ncol(dia)) {
      #dia[,i] <- dat[,i]^2 * C[i,i]
      dia[, i] <- X[, i]^2 * C[i, i]
    } 
    ndia <- matrix(NA, nrow(X), sum(upper.tri(C)))
    for (i in 1:(ncol(X) - 1)) {
      for (j in ((i + 1):ncol(X))) {
        ndia[, (1:ncol(ndia))[C[upper.tri(C)] == C[i, j]]] <-
          2 * X[, i] * X[, j] * C[i, j]
      }
    }
    Var <- apply(dia, 1, sum) + apply(ndia, 1, sum)
  }

  lower <- plogis(mu.hat - 2 * sqrt(Var))
  upper <- plogis(mu.hat + 2 * sqrt(Var))
  I <- data.frame(time = t, pred = pr, lower = lower, upper = upper)

  res <- list(I, q)
  names(res) <- c("predict", "q")
  res

}