R/polyreg-outputs.R

In cifmodeling: Visualization and Polytomous Modeling of Survival and Competing Risks

reportEffects <- function(outcome.type,
                          report.nuisance.parameter,
                          report.optim.convergence,
                          report.sandwich.conf,
                          report.boot.conf,
                          nuisance.model,
                          exposure,
                          estimand,
                          alpha_beta_estimated,
                          cov_estimated,
                          cov_bootstrap,
                          out_bootstrap,
                          out_getResults,
                          iteration,
                          converged.by,
                          objective.function,
                          max.absolute.difference,
                          relative.difference,
                          sol,
                          conf.int,
                          nleqslv.method) {

  out_getCoefTerm <- getCoefTerm(report.nuisance.parameter, report.sandwich.conf, report.boot.conf, nuisance.model, outcome.type, exposure, estimand, alpha_beta_estimated, cov_estimated, cov_bootstrap, out_bootstrap, conf.int, out_getResults)
  coef1 <- out_getCoefTerm$coef1
  coef2 <- out_getCoefTerm$coef2
  terms_text <- out_getCoefTerm$terms_text

  median_followup <- paste(round(median(out_getResults$t), 2))
  range_followup <- sprintf("[%.*f, %.*f]",
                        2, min(out_getResults$t),
                        2, max(out_getResults$t))
  #  range_followup <- paste("[",round(min(out_getResults$t), 2), ",",round(max(out_getResults$t), 2), "]")


  tg_event1_tidy <- tidy_df(coef1, terms_text)
  attr(tg_event1_tidy, "coefficient.index") <- out_getCoefTerm$index1
  if (outcome.type == "competing-risk" | outcome.type == "survival" | outcome.type == "binomial") {
    effect1_text <- paste(estimand$effect.measure1, "at", estimand$time.point)
  } else if (outcome.type == "proportional-survival" | outcome.type == "proportional-competing-risk") {
    effect1_text <- paste(estimand$effect.measure1, "of", exposure, "over time")
  }

  tg_event1_glance <- glance_df(
    report.optim.convergence,
    effect_text                  = effect1_text,
    event_text                   = paste(sum(out_getResults$y_1), "in N =", length(out_getResults$y_1)),
    median_followup              = median_followup,
    range_followup               = range_followup,
    param_len                    = length(alpha_beta_estimated),
    iteration                    = iteration,
    converged.by                 = converged.by,
    objective.function           = objective.function,
    max.absolute.difference      = max.absolute.difference,
    relative.difference          = relative.difference,
    nleqslv.value.at.solution    = max(sol$fvec),
    nleqslv.iteration            = sol$iter,
    nleqslv.message              = sol$message
  )

  if (outcome.type == "competing-risk" | outcome.type == "proportional-competing-risk") {
    if (outcome.type == "competing-risk") {
      effect2_text <- paste(estimand$effect.measure2, "at", estimand$time.point)
    } else if (outcome.type == "proportional-competing-risk") {
      effect2_text <- paste(estimand$effect.measure2, "of", exposure, "over time")
    }
    tg_event2_tidy <- tidy_df(coef2, terms_text)
    attr(tg_event2_tidy, "coefficient.index") <- out_getCoefTerm$index2
    tg_event2_glance <- glance_df(
      report.optim.convergence,
      effect_text                  = effect2_text,
      event_text                   = paste(sum(out_getResults$y_2), "in N =", length(out_getResults$y_2)),
      median_followup              = "-",
      range_followup               = "-",
      param_len                    = "-",
      iteration                    = "-",
      converged.by                 = "-",
      objective.function           = "-",
      max.absolute.difference      = "-",
      relative.difference          = "-",
      nleqslv.value.at.solution    = "-",
      nleqslv.iteration            = "-",
      nleqslv.message              = "-"
    )
    tg <- list(
      event1 = list(tidy = tg_event1_tidy, glance = tg_event1_glance),
      event2 = list(tidy = tg_event2_tidy, glance = tg_event2_glance)
    )
    class(tg$event1) <- class(tg$event2) <- "modelsummary_list"
    return(tg)
  }

  tg1 <- list(
    tidy   = tg_event1_tidy,
    glance = tg_event1_glance
  )
  class(tg1) <- "modelsummary_list"
  tg <- list("event 1 (no competing risk)" = tg1)
  return(tg)
}

tidy_df <- function(coef, text) {
  data.frame(
    term = text,
    estimate = coef$coef,
    std.error = coef$coef_se,
    p.value = coef$p_value,
    conf.low = coef$conf_low,
    conf.high = coef$conf_high
  )
}

glance_df <- function(
    report.optim.convergence, effect_text, event_text, median_followup, range_followup,
    param_len, iteration, converged.by, objective.function, max.absolute.difference,relative.difference,
    nleqslv.value.at.solution,nleqslv.iteration,nleqslv.message) {
  if (report.optim.convergence == FALSE) {
    data.frame(
      effect.measure = effect_text,
      n.events = event_text,
      median.follow.up = median_followup,
      range.follow.up = range_followup,
      n.parameters = param_len,
      converged.by = converged.by,
      nleqslv.message = nleqslv.message
    )
  } else {
    data.frame(
      effect.measure = effect_text,
      n.events  = event_text,
      median.follow.up = median_followup,
      range.follow.up = range_followup,
      n.parameters = param_len,
      n.iteration = iteration,
      converged.by = converged.by,
      objective.function = objective.function,
      max.absolute.difference = max.absolute.difference,
      relative.difference = relative.difference,
      nleqslv.value.at.solution = nleqslv.value.at.solution,
      nleqslv.iteration = nleqslv.iteration,
      nleqslv.message = nleqslv.message
    )
  }
}


getCoefTerm <- function(report.nuisance.parameter, report.sandwich.conf, report.boot.conf, nuisance.model, outcome.type, exposure, estimand, alpha_beta_estimated, cov_estimated, cov_bootstrap, out_bootstrap, conf.int, out_getResults) {
  iv <- estimand$index.vector
  code.exposure.ref <- as.character(estimand$code.exposure.ref)

  if (outcome.type == "competing-risk") {
    if (report.nuisance.parameter == TRUE) {
      index1 <- seq_len(iv[3])
      index2 <- seq.int(iv[4], iv[7])
      nuisance_terms <- c("Intercept", attr(terms(nuisance.model), "term.labels"))
    } else {
      index1 <- seq.int(iv[2], iv[3])
      index2 <- seq.int(iv[6], iv[7])
      nuisance_terms <- character(0)
    }
  }
  if (outcome.type == "survival" || outcome.type == "binomial") {
    if (report.nuisance.parameter == TRUE) {
      index1 <- seq_len(iv[3])
      index2 <- NULL
      nuisance_terms <- c("Intercept", attr(terms(nuisance.model), "term.labels"))
    } else {
      index1 <- seq.int(iv[2], iv[3])
      index2 <- NULL
      nuisance_terms <- character(0)
    }
  }
  if (outcome.type %in% c("proportional-survival", "proportional-competing-risk")) {
    layout <- reg_proportional_index(estimand)
    index1 <- layout$exposure
    nuisance_terms <- character(0)
    if (report.nuisance.parameter == TRUE) {
      index1 <- seq_len(layout$block)
      nuisance_terms <- sub("^\\(Intercept\\)$", "Intercept", colnames(out_getResults$x_l))
      if (length(estimand$time.point) > 1L) {
        nuisance_terms <- as.vector(outer(nuisance_terms, estimand$time.point,
          function(term, time) paste(term, "at", time)))
      }
    }
    index2 <- if (outcome.type == "proportional-competing-risk") index1 + layout$block else NULL
  }
  coef1 <- getCoef(index1, alpha_beta_estimated, cov_estimated, cov_bootstrap, outcome.type, report.sandwich.conf, report.boot.conf, out_bootstrap, conf.int)
  coef2 <- getCoef(index2, alpha_beta_estimated, cov_estimated, cov_bootstrap, outcome.type, report.sandwich.conf, report.boot.conf, out_bootstrap, conf.int)
  code.exposure.nonref <- colnames(out_getResults$x_a)
  terms_text_exposure <- if (length(code.exposure.nonref)) {
    paste0(exposure, ", ", sub("^a_", "", code.exposure.nonref), " vs ", code.exposure.ref)
  } else {
    character(0)
  }
  terms_text <- c(nuisance_terms, terms_text_exposure)
  return(list(coef1=coef1, coef2=coef2, terms_text=terms_text, index1=index1, index2=index2))
}

getCoef <- function(
    index,
    alpha_beta_estimated,
    cov_estimated,
    cov_bootstrap,
    outcome.type,
    report.sandwich.conf,
    report.boot.conf,
    out_bootstrap,
    conf.int) {
  alpha <- 1 - conf.int
  critical_value <- qnorm(1 - alpha / 2)
  if (report.sandwich.conf==TRUE) {
    coef <- alpha_beta_estimated[index]
    coef_se <- sqrt(diag(cov_estimated)[index])
    conf_low <- coef - critical_value * coef_se
    conf_high <- coef + critical_value * coef_se
    p_value <- 2 * (1 - pnorm(abs(coef) / coef_se))
  } else if (report.boot.conf==TRUE && outcome.type %in% c("competing-risk", "survival", "binomial")) {
    coef <- alpha_beta_estimated[index]
    coef_se <- sqrt(diag(cov_bootstrap)[index])
    conf_low <- coef - critical_value * coef_se
    conf_high <- coef + critical_value * coef_se
    p_value <- 2 * (1 - pnorm(abs(coef) / coef_se))
  } else if (report.boot.conf == TRUE) {
    coef <- out_bootstrap$boot.coef[index]
    coef_se <- out_bootstrap$boot.coef_se[index]
    conf_low <- out_bootstrap$boot.conf_low[index]
    conf_high <- out_bootstrap$boot.conf_high[index]
    p_value <- out_bootstrap$boot.p_value[index]
  } else {
    coef     <- alpha_beta_estimated[index]
    coef_se  <- rep(NA_real_, length(index))
    conf_low <- rep(NA_real_, length(index))
    conf_high <- rep(NA_real_, length(index))
    p_value   <- rep(NA_real_, length(index))
  }
  list(coef = coef, coef_se = coef_se, conf_low = conf_low, conf_high = conf_high, p_value = p_value)
}

reg_summarize_bootstrap <- function(fit, index, level, bca) {
  complete <- rowSums(!is.finite(fit$t)) == 0L
  if (sum(complete) < 2L) {
    stop("Fewer than two proportional bootstrap fits succeeded.", call. = FALSE)
  }
  replicates <- fit$t[complete, , drop = FALSE]
  covariance <- stats::cov(replicates)
  center <- colMeans(replicates)
  se <- sqrt(diag(covariance))
  margin <- stats::qnorm((1 + level) / 2) * se
  lower <- center - margin
  upper <- center + margin
  if (isTRUE(bca)) {
    center <- 2 * fit$t0 - center
    lower[] <- upper[] <- NA_real_
    for (j in index[se[index] > 0]) {
      limits <- tryCatch(boot.ci(fit, conf = level, index = j, type = "bca")$bca,
        error = function(e) {
          warning("BCa interval is unavailable for coefficient ", j, ". ",
            conditionMessage(e), call. = FALSE)
          NULL
        })
      if (!is.null(limits)) {
        lower[j] <- limits[1L, 4L]
        upper[j] <- limits[1L, 5L]
      }
    }
  }
  list(
    summary = list(boot.coef = center, boot.coef_se = se,
      boot.p_value = 2 * stats::pnorm(abs(center / se), lower.tail = FALSE),
      boot.conf_low = lower, boot.conf_high = upper),
    covariance = covariance
  )
}

extractOptimizationInfo <- function(sol, method) {
  out <- list(solver = method)
  if (method %in% c("nleqslv","Broyden","Newton")) {
    out$code    <- sol$termcd %||% NA_integer_
    out$message <- sol$message %||% NA_character_
    out$fn.evals <- sol$feval %||% NA_integer_
    out$jac.evals <- sol$jeval %||% NA_integer_
  } else if (method == "multiroot") {
    out$message <- sol$estim.precis %||% NA_character_
    out$iter    <- sol$iter %||% NA_integer_
    out$estim.precis <- sol$estim.precis %||% NA_character_
  } else if (method %in% c("optim","SANN","BFGS")) {
    out$code       <- sol$convergence %||% NA_integer_
    out$message    <- sol$message %||% NA_character_
    if (!is.null(sol$counts)) {
      out$fn.evals <- sol$counts["function"] %||% NA_integer_
      out$gr.evals <- sol$counts["gradient"] %||% NA_integer_
    }
  }
  return(out)
}

append_trace <- function(trace_df, iteration, computation.time.second = NA_real_, nleqslv.method, objective.function, relative.difference,
                         max.absolute.difference, converged.by, nleqslv.info, store_params = FALSE,
                         coefficient = NULL) {
  row <- data.frame(
    iteration = iteration,
    computation.time.second = computation.time.second,
    nleqslv.method = nleqslv.method,
    objective.function = objective.function,
    relative.difference = relative.difference,
    max.absolute.difference = max.absolute.difference,
    converged.by = if (isTRUE(converged.by)) NA_character_ else as.character(converged.by),
    code = nleqslv.info$code %||% NA_integer_,
    msg  = nleqslv.info$message %||% NA_character_,
    fn_evals = nleqslv.info$fn.evals %||% NA_integer_,
    gr_evals = nleqslv.info$gr.evals %||% NA_integer_,
    jac_evals = nleqslv.info$jac.evals %||% NA_integer_,
    stringsAsFactors = FALSE
  )
  if (isTRUE(store_params)) {
    row$coefficient <- list(coefficient)
  }
  if (is.null(trace_df)) return(row)
  rbind(trace_df, row)
}