#' Input types and model arguments of the pooling functions

#'

#' @description

#' Shared documentation of \code{\link{km_pooling}} and \code{\link{cox_pooling}}.

#' Both functions fit one model per imputed dataset and then pool the results with Rubin's rules.

#' This page describes which inputs they accept and which arguments are passed on to the

#' per-dataset model.

#'

#' @name pooling_inputs

#' @keywords internal

#' @md

#'

#' @param x The imputed datasets to analyze. One of:

#' - a list of data frames (one per imputed dataset), e.g., the output of

#'   \code{\link{raking_weights}} or of \code{MatchThem::complete(..., action = "all")}

#' - a weighted (\code{wimids}) or matched (\code{mimids}) object created with \pkg{MatchThem}

#' - an imputed (\code{mids}) object created with \pkg{mice}

#'

#' See the section "Input types and default model arguments" for how each input is handled.

#' @param ... Additional arguments passed on to the per-dataset model call

#' (\code{\link[ggsurvfit]{survfit2}} for \code{km_pooling()}, \code{\link[survival]{coxph}} for

#' \code{cox_pooling()}). They are used in addition to, and take precedence over, the default

#' arguments of the input type (see the section "Input types and default model arguments").

#' Only \code{formula} and \code{data} are set internally and cannot be supplied here.

#' Column names are written without quotes and are looked up in each dataset, e.g.

#' \code{weights = w} or \code{id = patientid}.

#'

#' @section Input types and default model arguments:

#'

#' The datasets are first brought into the form of a list of data frames. Depending on the type of

#' \code{x}, default model arguments are set:

#'

#' | Input \code{x} | How the datasets are obtained | \code{weights} | \code{cluster} |

#' |:---|:---|:---|:---|

#' | \code{wimids} | \code{MatchThem::complete(x, action = "all", all = FALSE, include = FALSE)} | \code{weights} | none |

#' | \code{mimids} | same as \code{wimids} | \code{weights} | \code{subclass} |

#' | \code{mids} | \code{mice::complete(x, action = "all", include = FALSE)} | none (unweighted) | none |

#' | list of data frames | used as supplied | none | none |

#'

#' **Lists of data frames are never interpreted automatically.** A list can come from many sources

#' (weighting, matching, raking, your own code) so the function cannot know which columns are meant

#' as weights or clusters. You state this explicitly through \code{...}, for example:

#' - weighted data: \code{weights = weights}

#' - matched data: \code{weights = weights, cluster = subclass}

#' - weights stored in a column with another name: \code{weights = w}

#'

#' If the datasets contain a column named \code{weights} or \code{subclass} that is not used

#' because no \code{weights} or \code{cluster} (or \code{id}) argument was specified, a warning

#' is raised that names the argument to add. Specify \code{weights = NULL} or \code{cluster = NULL}

#' to confirm that the column should deliberately be ignored.

#'

#' \strong{Important}: when you build the list from a \code{mimids} or \code{wimids} object yourself,

#' use \code{MatchThem::complete(x, action = "all", all = FALSE, include = FALSE)}. With

#' \code{all = TRUE}, unmatched patients and patients with zero weight are included in the models.

#'

#' @section Precedence of arguments in \code{...}:

#'

#' Default arguments of the input type are applied first, then the arguments you supply via

#' \code{...}. If both contain the same argument, \strong{your argument wins}. For example, for a

#' \code{mimids} object \code{cluster = NULL} fits the model without clustering on the matched sets

#' and \code{cluster = pair_id} clusters on another column. \code{formula} and \code{data} are the only

#' reserved arguments.

#'

#' \code{robust} is set automatically: unless you specify it, \code{robust = TRUE} is used whenever

#' weights, a cluster or an id are part of the model (the model-based variance is not valid in these

#' cases), otherwise the default of the model function is kept. An explicitly supplied \code{robust}

#' always takes precedence.

#'

#' @section Transparency:

#'

#' Every call prints a message that lists how the datasets were obtained and exactly which

#' arguments were used for the models, including where each argument comes from (default for the

#' input type, derived, or specified by you). Use \code{suppressMessages()} to silence it.

NULL

# prepare_pooling_input() -------------------------------------------------------

#' Turn the input of a pooling function into a list of datasets plus default arguments

#'

#' S3 generic with one small method per supported input type. Every method returns the same

#' three elements so that the pooling functions never need to look at the class of \code{x} again:

#' \itemize{

#'   \item \code{datasets}: plain list of data frames, one per imputed dataset

#'   \item \code{args}: named list of default model arguments as unevaluated expressions

#'     (e.g., \code{list(weights = quote(weights))}); column names stay symbols so that

#'     they are looked up in each dataset when the model is fitted

#'   \item \code{note}: one-line description of the input that is shown to the user

#' }

#' The methods have no side effects (no messages), which keeps them easy to test.

#'

#' @param x mimids, wimids or mids object, or a list of data frames

#'

#' @return list with elements \code{datasets}, \code{args} and \code{note}

#'

#' @noRd

prepare_pooling_input <- function(x) {

}

#' @noRd

#' @exportS3Method

prepare_pooling_input.wimids <- function(x) {

    # MatchThem always names the weight column "weights"

  )

}

#' @noRd

#' @exportS3Method

prepare_pooling_input.mimids <- function(x) {

    # MatchThem always names the weight column "weights" and the matched sets "subclass"

  )

}

#' @noRd

#' @exportS3Method

prepare_pooling_input.mids <- function(x) {

    # there are no weights or matched sets, so every observation is used

  )

}

#' @noRd

#' @exportS3Method

prepare_pooling_input.list <- function(x) {

  )

    # a list can come from anywhere (weighting, matching, raking, ...), so nothing

    # is assumed here; weights/cluster are specified by the user via `...`

  )

}

#' @noRd

#' @exportS3Method

prepare_pooling_input.mild <- function(x) {

  # the output of complete(action = "all") has class "mild" (not "list"), but is

  # a plain list of data frames, so it is handled like one

}

#' @noRd

#' @exportS3Method

prepare_pooling_input.default <- function(x) {

  # e.g., a single data frame or a vector

  )

}

#' Extract the matched/weighted imputed datasets of a mimids or wimids object

#'

#' @param x mimids or wimids object

#'

#' @return list of data frames, one per imputed dataset

#'

#' @noRd

complete_matched_weighted <- function(x) {

  # MatchThem objects need to contain the columns that the default arguments refer to

  )

    )

  }

  # action = "all"  : produces a mild object (list) of the multiply imputed datasets

  # all = FALSE     : do NOT include observations with a zero estimated weight

  # include = FALSE : do NOT include original data with the missing values

  # Important: the all = parameter needs to be set to FALSE

  # since otherwise unmatched patients or those with zero weight

  # will be included, too!

  )

}

# fit_models() --------------------------------------------------------------------

#' Fit one model per imputed dataset

#'

#' Shared by \code{km_pooling()} and \code{cox_pooling()}. The steps are:

#' \enumerate{

#'   \item merge default arguments of the input type with the arguments of the user

#'   \item derive \code{robust} if the user did not specify it

#'   \item check the arguments against the data and warn about unused weights/subclass columns

#'   \item tell the user how the models are fitted

#'   \item fit the model in each dataset

#' }

#'

#' @param input output of \code{prepare_pooling_input()}

#' @param fit_fun unevaluated name of the model function, e.g., \code{quote(survival::coxph)}

#' @param model_label model name shown to the user, e.g., "Cox proportional hazards model"

#' @param formula model formula

#' @param user_args named list of the unevaluated arguments the user supplied via \code{...},

#'   captured with \code{as.list(substitute(list(...)))[-1]} in the calling function

#'

#' @return list with the fitted models (\code{fits}, one per dataset) and the arguments used (\code{args})

#'

#' @noRd

fit_models <- function(input, fit_fun, model_label, formula, user_args) {

  # 1. merge arguments -----------------------------------------------------------

  # defaults first, then the user's arguments; if both contain the same

  # argument, the user's argument replaces the default

  # 2. derive robust -------------------------------------------------------------

  # a robust (sandwich) variance is required whenever weights, clusters or ids are used.

  # `is.null()` is used because weights = NULL / cluster = NULL mean "not used"

  }

  # 3. checks --------------------------------------------------------------------

  # 4. tell the user ------------------------------------------------------------

  )

  # 5. fit the model in each dataset -------------------------------------------

  # arguments that are NULL (e.g., cluster = NULL) are dropped, they mean "not used"

  # The call is built with as.call()/eval() rather than written literally with a

  # trailing `...`. ggsurvfit::survfit2() reconstructs its call internally

  # (match.call(expand.dots = TRUE)) to dispatch to survival::survfit(); if `...` is

  # spliced in literally, that reconstruction degrades NSE arguments like `id` into

  # positional `..1`-style placeholders that no longer resolve once the model tries to

  # evaluate them as columns of `data = d`. Building the call explicitly keeps

  # `formula` and `data` as plain symbols and the remaining arguments (e.g.,

  # weights = weights) as unevaluated symbols, which the model function then resolves

  # against `data = d`, exactly like a literal call would.

    ))

  }

}

#' Check that column names supplied as weights/cluster/id exist in the datasets

#'

#' @param args merged model arguments

#' @param datasets list of data frames

#'

#' @noRd

check_columns_exist <- function(args, datasets) {

    # only bare column names (symbols) can be checked; NULL or expressions are left alone

        ))

      }

    }

  }

}

#' Warn if a weights or subclass column exists but no matching argument was specified

#'

#' @param args merged model arguments

#' @param datasets list of data frames

#'

#' @noRd

warn_unused_columns <- function(args, datasets) {

  # `%in% names(args)` (rather than is.null) so that weights = NULL / cluster = NULL

  # count as a deliberate choice and silence the warning

  }

  }

    ))

  }

}

#' Tell the user how the models are fitted

#'

#' @inheritParams fit_models

#' @param args merged model arguments (before NULL arguments are dropped)

#' @param derived_robust whether \code{robust} was set automatically

#'

#' @noRd

report_fit_setup <- function(input, model_label, fit_fun, formula, args, user_args, derived_robust) {

  # curly braces have a meaning in cli messages, so escape them in text we insert

  # where does each argument come from?

    } else {

    }

  }

  } else {

      },

    )

  }

  ))

}

#' Pooled Kaplan-Meier estimate and survival curve

#'

#' @description

#' Computes pooled median survival Kaplan-Meier estimates using Rubin's rule

#' and outputs corresponding Kaplan-Meier curve across imputed and

#' matched/weighted datasets

#'

#' @details

#' The function takes the imputed datasets (\code{x}) as a list of data frames, or as a

#' mimids, wimids or mids object (see the sections "Input types and default model arguments" and

#' "Precedence of arguments in ..." below). The latter are the output of a workflow that requires

#' imputing multiple (m) datasets using mice or amelia and matching or weighting each imputed dataset

#' via the MatchThem package (see examples). Weights and clusters of lists of data frames are specified

#' explicitly via \code{...}, e.g., \code{weights = weights, cluster = subclass}.

#'

#' The function fits the pre-specified survfit model (\code{surv_formula}, \code{\link[survival]{survfit}} package)

#' to compute survival probabilities at each individual time point according to the Kaplan-Meier method.

#' Survival probabilities can be estimated at all observed event times (default) or at user-specified time points

#' using the \code{times} argument, which is passed directly to \code{\link[ggsurvfit]{tidy_survfit}}.

#' Arguments of the \code{\link[survival]{survfit}} call can be supplied via \code{...},

#' e.g. \code{type = "fleming-harrington"}. For matched and weighted datasets (mimids and wimids objects),

#' weights, cluster membership (matching only) and robust variance estimates are considered

#' by default. \code{...} only affects the underlying Kaplan-Meier fit of each imputed dataset;

#' the complementary log-log transformation, Rubin's rule pooling and confidence interval

#' computation are unaffected.

#'

#' Since survival probabilities typically don't follow normal distributions,

#' these need to be transformed to approximate normality first before pooling

#' across imputed datasets and time points. To that end, the method of Morisot et al. (2015) is

#' implemented, who applied Rubin's rules after a complementary log-log transformation of the

#' survival probabilities, following the recommendation of Marshall et al. (2009) to combine

#' predicted survival probabilities after this transformation. In each imputed dataset

#' \emph{i} and at each time point \emph{t}, the survival probability \eqn{S_i(t)} is transformed to

#' \deqn{Q_i(t) = \log(-\log(1 - S_i(t)))}

#' and the corresponding within-dataset variance is derived from the variance of the survival

#' probability (the squared standard error reported by \code{\link[survival]{survfit}}) with the delta method:

#' \deqn{U_i(t) = \frac{Var[S_i(t)]}{\left[(1 - S_i(t)) \log(1 - S_i(t))\right]^2}}

#'

#' Note that the complementary log-log transformation is applied to the survival probability

#' \eqn{S} itself, as in Morisot et al. (2015). This is not the same transformation as the

#' \eqn{\log(-\log(S))} used for the confidence intervals of \code{\link[survival]{survfit}}

#' with \code{conf.type = "log-log"}. For a single dataset, the point estimates and standard errors

#' therefore agree with \code{\link[survival]{survfit}}, but the confidence intervals differ slightly.

#'

#' To pool the transformed estimates across imputed datasets and time points, the \code{\link[mice]{pool.scalar}}

#' function is used to apply Rubin's rule to combine pooled estimates (qbar) according to

#' formula (3.1.2) Rubin (1987) and to compute the corresponding total variance (t) of the pooled

#' estimate according to formula (3.1.5) Rubin (1987), i.e.,

#' \eqn{T = \bar{U} + (1 + 1/m) B} with the average within-dataset variance \eqn{\bar{U}} and

#' the between-dataset variance \eqn{B} (as in Morisot et al. (2015)). If only a single dataset is supplied (m = 1, e.g., complete

#' data without imputation), the between-imputation variance is 0 by definition, the total variance equals the

#' within-dataset variance and a warning is raised. When \code{times = NULL}, the pooling is performed

#' at each unique combination of event times across all imputed datasets. When specific time points are provided

#' via the \code{times} argument, pooling occurs only at those pre-specified time points, ensuring consistency

#' across datasets and facilitating comparison at clinically relevant follow-up intervals.

#'

#' The pooled survival probabilities are then back-transformed via \code{1-exp(-exp(qbar))}

#' for pooled survival probability estimates and \code{1-exp(-exp(qbar +/- 1.96*sqrt(t)))}

#' for lower and upper 95% confidence intervals. As the formula indicates, the pooled standard

#' error is computed as the square root of the total variance. The vertically stacked table

#' with transformed and backtransformed estimates is returned with the

#' \code{km_survival_table} table.

#'

#' Special scenarios during the transformation, pooling and backtransformation are

#' considered. By convention, Kaplan-Meier survival curves start at time = 0 with survival probability = 1

#' unless a delayed entry survival model is fitted (not supported yet). In the case that there is no

#' observation in the data at at time = 0 with survival probability = 1, the \code{\link[ggsurvfit]{tidy_survfit}}

#' extrapolates this across strata. For the complementary log-log transformation this means, that survival probabilities = 1

#' at time = 0 and corresponding confidence intervals are extrapolated to 1. If the survival proability reaches exactly 0,

#' the survival probability and corresponding confidence intervals are set to 0 as well.

#'

#' Finally, the median survival time is extracted from the \code{km_survival_table} table

#' by determining the time the survival probability drops below .5 for the first time.

#' For this a sub-function of Terry M. Therneau's \code{\link[survival]{print.survfit}} function

#' is used. Therneau also considers some edge cases/nuisances (x = time, y = surv):

#' - Nuisance 1: if one of the y's is exactly .5, we want the mean of the corresponding x and the first x for which y<.5. We need to use the equivalent of all.equal to check for a .5 however: survfit(Surv(1:100)~1) gives a value of .5 + 1.1e-16 due to roundoff error.

#' - Nuisance 2: there may by an NA in the y's

#' - Nuisance 3: if no y's are <=.5, then we should return NA

#' - Nuisance 4: the obs (or many) after the .5 may be censored, giving a stretch of values = .5 +- epsilon

#'

#' \strong{Important Gotcha}: When using the \code{times} argument to specify custom time points, ensure that

#' sufficient time points are requested to capture the range of interest. If too few time points are specified,

#' the \code{km_median_survival} estimates and \code{km_plot} may not represent the true survival trajectory.

#' Specifically, if the median survival time lies between two requested time points or if requested time points

#' do not adequately cover the follow-up period, the median survival estimates will be inaccurate or return \code{NA}.

#' Similarly, the Kaplan-Meier curve will appear sparse or disconnected. For most applications, using the default

#' \code{times = NULL} is recommended to automatically use all observed event times. If using custom time points,

#' ensure they are sufficiently granular to capture the survival dynamics of interest.

#'

#' The function follows the following logic:

#' 1. Fit Kaplan-Meier survival function to each imputed and matched/weighted dataset

#' 2. Transform survival probabilities using complementary log-log transformation

#' 3. Pool transformed survival probabilities and compute total variance using Rubin's rule

#' 4. Back-transform pooled survival probabilities and compute 95% confidence intervals

#' 5. Extract median survival time and corresponding 95% confidence intervals

#' 6. Plot Kaplan-Meier curve with pooled survival probabilities and confidence intervals

#'

#' More references:

#' -  https://stefvanbuuren.name/fimd/sec-pooling.html

#'

#' @references

#' Morisot A, Bessaoud F, Landais P, et al. Prostate cancer: net survival and cause-specific survival

#' rates after multiple imputation. \emph{BMC Med Res Methodol}. 2015;15:54.

#' \doi{10.1186/s12874-015-0048-4}

#'

#' Marshall A, Altman DG, Holder RL, Royston P. Combining estimates of interest in prognostic

#' modelling studies after multiple imputation: current practice and guidelines.

#' \emph{BMC Med Res Methodol}. 2009;9:57. \doi{10.1186/1471-2288-9-57}

#'

#' Rubin DB. \emph{Multiple Imputation for Nonresponse in Surveys}. Wiley; 1987.

#'

#' @seealso

#' \code{\link[survival]{survfit}} \code{\link[mice]{pool.scalar}} \code{\link[MatchThem]{matchthem}} \code{\link[MatchThem]{weightthem}} \code{\link[ggsurvfit]{survfit2}} \code{\link[ggsurvfit]{tidy_survfit}}

#'

#' @inheritParams pooling_inputs

#' @param surv_formula specification of survival model formula to be fitted

#' @param times numeric vector of follow-up time points at which survival probabilities are to be estimated (propagated to \code{\link[ggsurvfit]{tidy_survfit}}).

#' If NULL (default), survival probabilities are estimated and pooled across all observed event times.

#'

#' @inheritSection pooling_inputs Input types and default model arguments

#' @inheritSection pooling_inputs Precedence of arguments in \code{...}

#' @inheritSection pooling_inputs Transparency

#'

#' @return list with pooled median survival estimate and pooled Kaplan-Meier curve

#' \code{km_median_survival}:

#' -  strata = stratum

#' -  t_median = median survival time

#' -  t_lower = lower 95% CI of median survival time

#' -  t_upper = upper 95% CI of median survival time

#'

#' \code{km_survival_table}:

#' -  strata = stratum

#' -  time = observed time point

#' -  m = number of imputed datasets

#' -  qbar = pooled univariate estimate of the complementary log-log transformed survival probabilities \eqn{\log(-\log(1 - S))}, see formula (3.1.2) Rubin (1987)

#' -  t = total variance of the pooled univariate estimate of the complementary log-log transformed survival probabilities, formula (3.1.5) Rubin (1987)

#' -  se = total standard error of the pooled estimate (derived as sqrt(t))

#' -  surv = back-transformed pooled survival probability

#' -  lower = Wald-type lower 95% confidence interval of back-transformed pooled survival probability

#' -  upper = Wald-type upper 95% confidence interval of back-transformed pooled survival probability

#'

#' \code{km_plot}: ggplot2 object with Kaplan-Meier curve

#'

#' @md

#'

#' @export

#'

#' @examples

#'  library(encore.analytics)

#'  library(mice)

#'  library(MatchThem)

#'

#'  # simulate a cohort with 500 patients with 20% missing data

#'  data <- simulate_data(

#'    n = 500,

#'    imposeNA = TRUE,

#'    propNA = 0.2

#'    )

#'

#'  # impute the data

#'  set.seed(42)

#'  mids <- mice(data, m = 5, print = FALSE)

#'

#'  # fit a propensity score model

#'  fit <- as.formula(treat ~ dem_age_index_cont + dem_sex_cont + c_smoking_history)

#'

#'  # weight (or alternatively match) patients within each imputed dataset

#'  wimids <- weightthem(

#'    formula = fit,

#'    datasets = mids,

#'    approach = "within",

#'    method = "glm",

#'    estimand = "ATO"

#'    )

#'

#'  # fit a survival model

#'  km_fit <- as.formula(survival::Surv(fu_itt_months, death_itt) ~ treat)

#'

#'  # 1) weighted (wimids) or matched (mimids) object:

#'  # weights (and clusters) are set automatically

#'  km_out <- km_pooling(

#'    x = wimids,

#'    surv_formula = km_fit

#'    )

#'

#'  # median survival time

#'  km_out$km_median_survival

#'

#'  # KM curve

#'  km_out$km_plot

#'

#'  # 2) list of data frames: weights are specified explicitly via `...`

#'  # (all = FALSE excludes observations with zero weight)

#'  datasets <- MatchThem::complete(wimids, action = "all", all = FALSE, include = FALSE)

#'  km_out_list <- km_pooling(

#'    x = datasets,

#'    surv_formula = km_fit,

#'    weights = weights

#'    )

#'

#'  # for matched datasets, additionally specify the matched sets as clusters:

#'  # km_pooling(x = datasets, surv_formula = km_fit, weights = weights, cluster = subclass)

#'

#'  # weights stored in a column with another name, e.g., `w`

#'  datasets_w <- lapply(datasets, function(d) dplyr::rename(d, w = weights))

#'  km_out_w <- km_pooling(

#'    x = datasets_w,

#'    surv_formula = km_fit,

#'    weights = w

#'    )

#'

#'  # 3) imputed (mids) object: unweighted Kaplan-Meier curve

#'  km_out_mids <- km_pooling(

#'    x = mids,

#'    surv_formula = km_fit

#'    )

#'

#'  # 4) pass additional arguments to the underlying survfit2() call:

#'  # column names are resolved in each imputed dataset, e.g. `id = patientid`.

#'  # survfit() uses `id` as the subject identifier for its robust/sandwich

#'  # variance calculation, clustering rows by patient instead of treating

#'  # each row as its own independent cluster (relevant if a patient can

#'  # contribute more than one row, e.g. with time-varying covariates).

#'  km_out_id <- km_pooling(

#'    x = wimids,

#'    surv_formula = km_fit,

#'    id = patientid

#'    )

#'

km_pooling <- function(

  x = NULL,

  surv_formula = stats::as.formula(

    survival::Surv(fu_itt_months, death_itt) ~ treat

  ),

  times = NULL,

  ...

) {

  # input checks ------------------------------------------------------------

  # check if surv_formula is a formula

  )

  # bring the input into the form of a list of datasets plus the default

  # model arguments of the input type (weights, cluster, ...)

    )

  }

  # check that ... does not try to override arguments set internally

  # (...names() reads argument names without forcing/evaluating them, since

  # some, like `id`, are intended to be resolved as columns of the imputed

  # data further down rather than as objects that exist here)

  )

    )

  )

  # capture the *unevaluated* expressions behind `...` (e.g. list(id = patientid)

  # becomes list(id = quote(patientid))). This must not be `list(...)`: some

  # arguments (like `id`) are column names meant to be resolved against

  # `data = d` inside survfit(), not objects that exist in this environment,

  # and `list(...)` would force and fail on them immediately.

  # Fit Kaplan Meier --------------------------------------------------------

  # fit a survival function in each imputed (and weighted/matched) dataset;

  # merging of default and user arguments, checks, the console message and the

  # construction of the survfit2() call are handled in fit_models()

  # with a single dataset (e.g., complete data without imputation) there is no

  # between-imputation variance; the total variance is then the within-dataset

  # variance (see pooling below), which we communicate to the user

    ))

  }

  # pool survival probabilities ---------------------------------------------

  # now to combine/pool quantities with non-normal distributions such as

  # survival probabilities, we apply Rubin's rules after a complementary log-log

  # transformation of the survival probability S, following

  #   Morisot A, Bessaoud F, Landais P, et al. BMC Med Res Methodol 2015;15:54

  #   https://doi.org/10.1186/s12874-015-0048-4

  # who build on the recommendation of

  #   Marshall A, Altman DG, Holder RL, Royston P. BMC Med Res Methodol 2009;9:57

  #   https://doi.org/10.1186/1471-2288-9-57

  # see also https://stefvanbuuren.name/fimd/sec-pooling.html

  #

  # Per dataset i and time t (Morisot et al., 2015):

  #   Q_i = log(-log(1 - S_i))                                 -> q below

  #   U_i = Var[S_i] / [(1 - S_i) * log(1 - S_i)]^2 (delta method) -> u below

  # Note: the transformation is applied to S itself, which is not the

  # log(-log(S)) used by survival::survfit(conf.type = "log-log").

  # that is, the function cloglog performs the

  # transformations on all survival probabilities

  # stored in survival_fit_list

    # tidy up

    ) |>

        # compute transformed qbar

        # compute transformed U

      )

  }

  # we apply this function to all fitted survival curves

    # the resulting list with q and u per stratum and dataset

    # is then vertically stacked to a long data frame (using bind_rows).

    # group by the exposure levels (stratum) and apply the pool.scalar function

    # which pools univariate estimates of m datasets. This returns:

    # m = The number of imputations

    # qbar = The pooled univariate estimate, formula (3.1.2) Rubin (1987)

    # ubar = The average within-imputation variance, formula (3.1.3) Rubin (1987)

    # t = The total variance of the pooled estimated, formula (3.1.5) Rubin (1987).

    ) |>

    # with a single dataset, the between-imputation variance is not defined

    # (pool.scalar returns NA for t); it is 0 by definition, so the total variance

    # is the within-dataset variance (ubar). This is deliberately based on the

    # number of datasets, not on m per time point: with several datasets, a time

    # point observed in only one of them stays NA

    # we back-transform back and compute CI's

    ) |>

    # confidence intervals are both 1 for time = 0 and surv = 1

    ) |>

    # confidence intervals are both 0 for surv = 0

    ) |>

    # arrange by strata, time and descending survival probability

  # Plot a pooled Kaplan-Meier curve ----------------------------------------

  # plot pooled Kaplan-Meier curve based on pooled survival probabilities

    )) +

    ) +

    ) +

  # Extract median survival time --------------------------------------------

  # This section is taken from Terry M. Therneau's

  # print.survfit function in the survival package!

  # to extract the median survival times, we

  # apply the same sub-functions used in

  # print.survfit function of the original survival package

  # which defines the median survival time as the time

  # the survival probability drops below .5 for the first time

  # see https://github.com/cran/survival/blob/master/R/print.survfit.R#L109

  # Therneau defined some edge cases (x = time, y = surv):

  # Nuisance 1: if one of the y's is exactly .5, we want the mean of

  #    the corresponding x and the first x for which y<.5.  We need to

  #    use the equivalent of all.equal to check for a .5 however:

  #    survfit(Surv(1:100)~1) gives a value of .5 + 1.1e-16 due to

  #    roundoff error.

  # Nuisance 2: there may by an NA in the y's

  # Nuisance 3: if no y's are <=.5, then we should return NA

  # Nuisance 4: the obs (or many) after the .5 may be censored, giving

  #   a stretch of values = .5 +- epsilon

    } else {

      } else {

      }

    }

  }

  # apply this across strata listed in the km_survival table

    # median survival time

    # 95% CI of median survival time

    # = time when lower and upper

    # 95% confidence intervals of the survival

    # probability drop below .5

      )

      )

    } else {

      # Therneau assigns a 0 if no confidence intervals are available

    }

    # provide output table

    )

  }

  )

  # return results ----------------------------------------------------------

    )

  )

}

#' Calculate Standardized Mean Difference (SMD) Between RCT and RWE Results

#'

#' @description

#' This function calculates the standardized mean difference (SMD) between results from

#' a randomized controlled trial (RCT) and its real-world evidence (RWE) emulation.

#' The SMD provides a measure of agreement between the two results, taking into account

#' both the point estimates and their corresponding uncertainty.

#'

#' @details

#' The SMD is calculated as the difference between RCT and RWE estimates, divided by

#' the square root of the sum of their variances. Variances are derived from the

#' confidence intervals assuming normal distribution (using 1.96 for 95% CI).

#'

#' The formula used is:

#' SMD = (theta_RCT - theta_RWE) / sqrt(Var(theta_RCT) + Var(theta_RWE))

#'

#' @param rct_estimate Numeric. Point estimate from the RCT (typically log hazard ratio or log odds ratio)

#' @param rct_lower Numeric. Lower bound of the confidence interval from RCT

#' @param rct_upper Numeric. Upper bound of the confidence interval from RCT

#' @param rwe_estimate Numeric. Point estimate from the RWE study

#' @param rwe_lower Numeric. Lower bound of the confidence interval from RWE

#' @param rwe_upper Numeric. Upper bound of the confidence interval from RWE

#'

#' @return

#' Returns a numeric value representing the standardized mean difference.

#' Larger absolute values indicate greater disagreement between RCT and RWE results.

#'

#' @examples

#' # Example with hazard ratios (log scale)

#' smd <- smd_agreement(

#'   rct_estimate = log(0.87),

#'   rct_lower = log(0.78),

#'   rct_upper = log(0.97),

#'   rwe_estimate = log(0.82),

#'   rwe_lower = log(0.76),

#'   rwe_upper = log(0.87)

#'   )

#'

#' @references

#' For methodology on comparing RCT and RWE results using standardized differences.

#'

#' @export

smd_agreement <- function(rct_estimate,

                          rct_lower,

                          rct_upper,

                          rwe_estimate,

                          rwe_lower,

                          rwe_upper

                          ){

  # Input validation

  # Check for missing values

  }

  # Check numeric inputs

  }

  # Check bounds

  }

  }

  # Check if estimate is within bounds

  }

  }

  # numerator

  # variances

  # denominator

  # Final check for valid calculation

  }

  # smd

}

#' Simulates an artifical EHR-derived analysis-ready oncology dataset

#'

#' @description

#' Parameterized function to quickly create an

#' EHR-derived analytic cohort for analytic code development.

#'

#' @details

#' This function simulates a cohort of patients with

#' oncology data. The cohort is simulated using a Weibull distribution for time to event

#' and a logistic distribution for treatment assignment.

#' The function also allows for missingness to be imposed on the data.

#' The function is parameterized to allow for the number of patients

#' to be simulated, the seed for reproducibility, and whether to

#' include a patient id variable.

#'

#' @param n_total integer, number of total patients

#' @param seed integer, seed for reproducibility

#' @param include_id logical, include a generated patientid variable

#' @param imposeNA logical, set covariates to missing

#' @param propNA numeric, proportion of missingness, needs to be between 0 and 1

#'

#' @return data frame with simulated analytic cohort

#'

#' @export

#'

#' @examples

#' \dontrun{

#' library(encore.analytics)

#'

#' # Original uniform missingness

#' data_uniform <- simulate_data(

#'   n_total = 3500,

#'   seed = 41,

#'   include_id = FALSE,

#'   imposeNA = TRUE,

#'   propNA = .33

#'   )

#'

#' # Now creates variable missingness across columns with average of 0.33

#' data_variable <- simulate_data(

#'   n_total = 3500,

#'   seed = 41,

#'   include_id = FALSE,

#'   imposeNA = TRUE,

#'   propNA = .33

#'   )

#'

#' head(data_variable)

#'

#'}

simulate_data <- function(n_total = 3500,

                          seed = 42,

                          include_id = TRUE,

                          imposeNA = TRUE,

                          propNA = NULL

                          ){

  # check input

  # set the seed

  # Build the total cohort

    # patient id

    # age

    # sex/male

    # smoking

    # num met sites

    # c_ecog_cont

    # c_stage_initial_dx_cont

    # dem_race

    # dem_region

    # dem_ses

    # c_hemoglobin_g_dl_cont

    # c_urea_nitrogen_mg_dl_cont

    # c_platelets_10_9_l_cont

    # c_calcium_mg_dl_cont

    # c_glucose_mg_dl_cont

    # c_lymphocyte_leukocyte_ratio_cont

    # c_alp_u_l_cont

    # c_protein_g_l_cont

    # c_alt_u_l_cont

    # c_albumin_g_l_cont

    # c_bilirubin_mg_dl_cont

    # c_chloride_mmol_l_cont

    # c_monocytes_10_9_l_cont

    # c_eosinophils_leukocytes_ratio_cont

    # c_ldh_u_l_cont

    # c_hr_cont

    # c_sbp_cont

    # c_oxygen_cont

    # c_neutrophil_lymphocyte_ratio_cont

    # c_bmi_cont

    # c_ast_alt_ratio_cont

    # c_time_dx_to_index

    # c_de_novo_mets_dx

    # c_height_cont

    # c_weight_cont

    # c_dbp_cont

    # c_year_index

    )

  # assign treatment probabilities

      # treat (to simplify categorical variables don't get betas)

          )

        )

      )

  #summary(cohort$treat)

  # simulate HTE by race

  ## Asian patients

  # assign betas for hazard model

    )

        )

      ) |>

      )

  ## Non-Asian patients

    )

        )

      ) |>

      )

  # combine Asian and non-Asian

  # summary(cohort_outcome$death_itt)

  # summary(cohort_outcome$fu_itt_months)

  # hist(cohort_outcome$fu_itt_months)

    # variables that should always stay complete

    # Create variable missingness proportions per column with average = propNA

    # Safety check - ensure we have variables to make missing

      # No variables to make missing, skip missingness

      }

    # Generate proportions that average to propNA using a simpler approach

    # Start with uniform random values and scale them

    # Scale to have the desired mean

    # Apply adaptive bounds based on target propNA

    # Final adjustment to ensure mean equals propNA

      }

    # Apply missingness column by column

      # Skip if proportion is NA or invalid

      # Randomly select rows to set as missing

        }

      }

    }

}

#' Estimate raking weights for a mimids or wimids object

#'

#' @description

#' Function estimates raking weights for multiple imputed and matched (\code{mimids})

#' or weighted (\code{wimids}) datasets. That is to match the distributions of

#' certain variables in the imputed and matched/weighted datasets to a distribution of

#' a target population (e.g., clinical trial population).

#'

#' @details

#' The function requires an object of class mimids or wimids (\code{x}), which is the output

#' of a workflow that requires imputing multiple (m) datasets using mice or amelia

#' and matching or weighting each imputed dataset via the MatchThem package

#' (see examples).

#'

#' The function additionally requires a list of target distributions (\code{targets}) for each variable

#' that is considered for the raking procedure. The list should contain named vectors with the target

#' distributions for each variable and the names of the vectors should match the variable names in the imputed datasets.

#'

#' In brief, the raking procedure iteratively adjusts the weights to make the weighted sample percentages match the target

#' population percentages for the selected variables.It does this by multiplying the current weight for each case by a

#' factor based on the ratio of the target population proportion to the weighted sample proportion for a given category.

#' This adjustment is performed sequentially for each category of each selected variable. Because adjusting for one

#' variable can disrupt the match for previous variables, the process is repeated through all selected variables in cycles.

#' This iterative process minimizes the Kullback-Leibler (KL) divergence and continues until the weighted sample proportions

#' match the target population proportions for all categories ("full convergence"), or until no further change occurs.

#'

#' The function follows the following logic:

#' 1. Extract the ith imputed dataset from the mimids or wimids object

#' 2. Create a temporary case/patient ID

#' 3. Apply the \code{\link[anesrake]{anesrake}} function to the ith imputed dataset

#' 4. Create a temporary dataframe with the case ID and replace the initial weights with the updated raking weights

#' 5. Merge the temporary dataframe with the ith imputed dataset

#' 6. Drop the temporary case ID

#' 7. Return the ith imputed dataset with the raking weights

#'

#' The function returns a list of data frames with the updated raking weights. These

#' updated raking weights overwrite in each data frame the existing \code{weights} column.

#' This column can then be used in a downstream analysis (e.g., Kaplan-Meier, Cox proportional hazards regression).

#'

#' @seealso

#' \code{\link[anesrake]{anesrake}} \code{\link[MatchThem]{matchthem}} \code{\link[MatchThem]{weightthem}}

#'

#' @param x imputed and matched (mimids) or weighted (wimids) object

#' @param targets list of all target values for the raking procedure

#'

#' @return list of data frames with updated raking weights (\code{raking_weights})

#'

#' @md

#'

#' @export

#'

#' @examples

#'

#'  library(encore.analytics)

#'  library(mice)

#'  library(dplyr)

#'  library(MatchThem)

#'  library(survival)

#'

#'  # simulate a cohort with 1,000 patients with 20% missing data

#'  data <- simulate_data(

#'    n = 500,

#'    imposeNA = TRUE,

#'    propNA = 0.2

#'    ) |>

#'    # anesrake works best with factor variables

#'    mutate(c_smoking_history = factor(ifelse(c_smoking_history == TRUE, "Current/former", "Never")))

#'

#'  # impute the data (create mids object)

#'  set.seed(42)

#'  mids <- mice(data, m = 5, print = FALSE)

#'

#'  # define covariates for propensity score model

#'  covariates <- data |>

#'   select(starts_with("c_"), starts_with("dem_")) |>

#'    colnames()

#'

#'  # define propensity score model

#'  fit <- as.formula(paste("treat ~", paste(covariates, collapse = " + ")))

#'

#'  # match patients within each imputed dataset

#'  mimids <- matchthem(

#'    formula = fit,

#'    datasets = mids,

#'    approach = 'within',

#'    method = 'nearest'

#'    )

#'

#'  smoker_target <- c(.35, .65)

#'  names(smoker_target) <- c("Current/former", "Never")

#'

#'  # summarize target distributions in a named list vector

#'  targets <- list(smoker_target)

#'  names(targets) <- c("c_smoking_history")

#'

#'  # estimate raking weights

#'  mirwds <- raking_weights(

#'    x = mimids,

#'    targets = targets

#'    )

#'

raking_weights <- function(x,

                           targets = NULL

                           ){

  # input checks----

  # check if object is a mimids or wimids object

  # check if targets is a list

  # check if targets is a named list

  # check if targets is a list of named vectors

  # check if targets is a list of named vectors

  # extract imputed and matched/weighted datasets and weights ----

  # create a function to indicate which cases should be re-weighted

  # since we are only interested in matched/weighted cases

    # next, we create a list of datasets to apply the raking procedure one-by-one

    # action = "all"  : produces a mild object (list) of the multiply imputed datasets

    # all = FALSE     : do NOT include observations with a zero estimated weight

    # include = FALSE : do NOT include original data with the missing values

    # Important: the all = parameter needs to be set to FALSE

    # since otherwise unmatched patients or those zero weight

    # will be included, too!

    # create a temporary case ID

    # apply anesrake function to re-weight -----------------------------------

      # target distributions for raking procedure

      # data to be weighted

      # patient identifier

      # other weights that should be accounted for before re-weighting is conducted

      # for matching this is 0 (unmatched) and 1 (matched) and for weighting this is the actual weights

      # no verbosity

      )

    # create a temporary dataframe with id and sampling/re-weighting weights

      # remove "old" weights with new weights

      # now we can drop the temporary id

  }

  # return list of imputed, matched/weighted and re-weighted datasets (mirwds)

    )

}

#' Calculate Agreement Metrics Between RCT and RWE Results

#'

#' @description Creates a formatted table comparing results from randomized controlled trials (RCTs)

#' with real-world evidence (RWE) emulation studies by calculating agreement metrics.

#'

#' @details

#' Calculates up to three types of agreement metrics, controlled via the `metrics`

#' argument:

#' - Statistical significance agreement: RCT and RWE confidence intervals fall on the

#'   same side of the null value of 1 (both entirely above, both entirely below, or

#'   both include the null, i.e. both non-significant)

#' - Estimate agreement: RWE estimate within RCT confidence interval