#' 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 = (θ_RCT - θ_RWE) / sqrt(Var(θ_RCT) + Var(θ_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)

#'

#' data_miss <- simulate_data(

#'   n_total = 3500,

#'   seed = 41,

#'   include_id = FALSE,

#'   imposeNA = TRUE,

#'   propNA = .33

#'   )

#'

#' head(data_miss)

#'

#'}

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)

    # define column numbers for pattern and weight determination

    # define missingness pattern

  }

}

#' 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)

    )

}

#' Manually fit and pool Cox proportional hazards model results from multiple imputed datasets

#'

#' @description

#' This function manually fits and pools results from Cox proportional hazards models using

#' a list of imputed datasets. It leverages the `mice::as.mira` and `mice::pool` functions to

#' ensure proper pooling of results across multiple imputations using Rubin's rules.

#'

#' @details

#' The function requires a list of imputed data frames with weights and optional cluster (matching) information,

#' as well as a formula for the Cox proportional hazards model. The data frames must include a column named

#' \code{weights}, and optionally a column named \code{subclass} (for matched datasets to indicate cluster membership).

#'

#' This function is particularly useful when working with imputed datasets that are not in the form of

#' \code{mimids} or \code{wimids} objects, such as when intermediate steps like raking weights

#' (via \code{raking_weights}) are applied. It provides a flexible way to fit and pool Cox models

#' while ensuring compatibility with Rubin's rules for multiple imputation.

#'

#' The function follows these steps:

#' 1. Fit a Cox proportional hazards model to each imputed dataset. If a \code{subclass} column is present,

#'    it is used as a cluster variable for matched pairs.

#' 2. Convert the list of fitted models into a \code{mira} object using \code{mice::as.mira}.

#' 3. Pool the results using \code{mice::pool}, which applies Rubin's rules for combining estimates

#'    and variances across imputations.

#' 4. Format the pooled results, including exponentiating the hazard ratios and calculating confidence intervals.

#'

#' @seealso

#' \code{\link[survival]{coxph}}, \code{\link[mice]{pool}}, \code{\link[mice]{as.mira}},

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

#'

#' @param x A list of imputed datasets with weights or raking weights (e.g., from \code{\link[encore.analytics]{raking_weights}})

#'          and optional cluster information (for matched datasets).

#' @param surv_formula A formula for the Cox proportional hazards model (default is \code{Surv(fu_itt_months, death_itt) ~ treat}).

#'

#' @return A data frame containing the pooled results, including hazard ratios, confidence intervals, and p-values.

#'

#' @md

#'

#' @export

#'

#' @examples

#' library(encore.analytics)

#' library(mice)

#' library(MatchThem)

#'

#' # Simulate a cohort with 500 patients and 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 patients within each imputed dataset

#' wimids <- weightthem(

#'   formula = fit,

#'   datasets = mids,

#'   approach = "within",

#'   method = "glm",

#'   estimand = "ATO"

#' )

#'

#' # Create a list of imputed and weighted datasets

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

#'

#' # Define a survival model formula

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

#'

#' # Fit and pool Cox proportional hazards model results

#' cox_pooling(wimids_list, surv_formula = cox_fit)

cox_pooling <- function(x,

                        surv_formula = stats::as.formula(survival::Surv(fu_itt_months, death_itt) ~ treat)

                        ){

  # input checks----

  # check if object is a list

  # check if surv_formula is a formula

  # check if weights_col is in the data

  # Step 1: fit Cox model to each imputed dataset----

    # if cluster_col is in the data (matched datasets)

        )

    }else{

        )

      }

  }

  # loop cox_fit function over each imputed dataset

    )

  # Step 2: extract results into a format suitable for pooling

  # convert each model into a "mira" object (required for pooling)

  # Step 3: pool the results using Rubin's rules

  # Step 4: format the pooled results

}

#' Wrapper around \code{gtsummary::tbl_summary()} to create a beautiful Table 1 quickly

#'

#' @description

#' Builds a table 1 using the \code{gtsummary} package. The function takes a one-row-per-patient dataframe

#' and creates a summary table of covariates stratified by treatment group.

#'

#' @details

#' Wrapper for \code{\link[gtsummary]{tbl_summary}}.

#'

#' The function \code{create_table1} is a wrapper around the \code{gtsummary::tbl_summary()} function.

#' It is designed to create a summary table of covariates stratified by treatment group.

#' The function requires a one-row-per-patient dataframe and creates a summary table of covariates stratified by treatment group.

#'

#' @importFrom smd smd

#'

#' @param x dataframe with individual-level patient data in a one-row-per-patient format with treatment stratification variable and covariates to be displayed in the Table 1

#' @param covariates character vector of columns/covariate names to be displayed in Table 1

#' @param covariates_labels named character vector or list of formulas specifying variables labels of covariate-label pairs to display in table

#' @param treat character specifying column name of treatment variable

#' @param explicit_na_categorical logical, should missings in categorical variables be explicitly included as a separate category (default is TRUE)

#'

#' @return object of class "tbl_summary" "gtsummary"

#'

#' @export

#'

#' @examples

#' \dontrun{

#' library(encore.analytics)

#'

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

#' data <- simulate_data(

#'   n = 1000,

#'   imposeNA = TRUE,

#'   propNA = 0.2

#'   )

#'

#' # create a Table 1

#' create_table1(

#'   x = data,

#'   covariates = c("dem_age_index_cont", "dem_sex_cont", "c_smoking_history"),

#'   covariates_labels = list(

#'     "dem_age_index_cont" = "Age",

#'     "dem_sex_cont" = "Sex",

#'     "c_smoking_history" = "Smoking history"

#'     ),

#'   treat = "treat"

#'   )

#'}

create_table1 <- function(x = NULL,

                          covariates = NULL,

                          covariates_labels = NULL,

                          treat = "treat",

                          explicit_na_categorical = TRUE

                          ){

  # check input

  # assign treatment variable

  # and subset columns

    # encode NA explicitly as "Missing" (only applicable to categorical)

      # logicals into categorical

          )

        ) |>

      # for all categorical, assign a "Missing" category

          )

        ) |>

      # drop "Missing" level for covariates which are fully observed/have no NA

  }

      ) |>

      #test_result ~ "**SMD**"

      ) |>

}

#' 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 three types of agreement metrics:

#' - Statistical significance agreement: Same direction and statistical significance

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

#' - SMD agreement: Standardized mean difference below threshold (default 1.96)

#'

#' All estimates should be hazard ratios (HR) and must be positive.

#' Estimates are log-transformed for SMD calculation.

#'

#' @importFrom tidyr pivot_wider pivot_longer

#' @importFrom smd smd

#'

#' @param x A data.frame or tibble containing RCT and RWE results (in wide format, see examples).

#' Must contain these columns:

#'

#' - rct_estimate = Point estimate from the RCT

#' - rct_lower = Lower confidence limit of the RCT result

#' - rct_upper = Upper confidence limit of the RCT result

#' - rwe_estimate = Point estimate from the RWE emulation study

#' - rwe_lower = Lower confidence limit of the RWE result

#' - rwe_upper = Upper confidence limit of the RWE result

#'

#' @param analysis_col Character. Name of column identifying different analyses

#' @param group_col Character. Optional. Name of column for grouping results (e.g., database)

#' @param estimate_label Character. Label for estimates in table header. Default: "HR (95% CI)"

#' @param smd_threshold Numeric. Threshold for SMD agreement. Default: 1.96 (α=0.05)

#'

#' @return A gt table object with formatted agreement metrics

#'

#' @md

#'

#' @examples

#' # Example 1: One analysis per database with grouping

#'

#'x <- tibble::tribble(

#'~Analysis, ~Database, ~rct_estimate, ~rct_lower, ~rct_upper, ~rwe_estimate, ~rwe_lower, ~rwe_upper,

#'"Main analysis", "Database 1", 0.87, 0.78, 0.97, 0.82, 0.76, 0.87,

#'"Main analysis", "Database 2", 0.5, 0.4, 0.6, 2.0, 1.8, 2.2,

#'"Main analysis", "Database 3", 0.8, 0.7, 0.9, 1.5, 1.4, 1.6

#')

#'

#' agreement_metrics(x, analysis_col = "Analysis", group_col = "Database")

#'

#' # Example 2:  Three analyses, one database without grouping

#'

#' x <- tibble::tribble(

#' ~Analysis, ~rct_estimate, ~rct_lower, ~rct_upper, ~rwe_estimate, ~rwe_lower, ~rwe_upper,

#' "Main analysis", 0.87, 0.78, 0.97, 0.82, 0.76, 0.87,

#' "Sensitivity analysis 1", 0.5, 0.4, 0.6, 2.0, 1.8, 2.2,

#' "Sensitivity analysis 2", 0.8, 0.7, 0.9, 1.5, 1.4, 1.6

#' )

#'

#' agreement_metrics(x, analysis_col = "Analysis")

#'

#' # Example 3: Three analyses per database with grouping

#' x <- tibble::tribble(

#' ~Analysis, ~Database, ~rct_estimate, ~rct_lower, ~rct_upper, ~rwe_estimate, ~rwe_lower, ~rwe_upper,

#' "Main analysis", "Database 1", 0.87, 0.78, 0.97, 0.82, 0.76, 0.87,

#' "Main analysis", "Database 2", 0.5, 0.4, 0.6, 2.0, 1.8, 2.2,

#' "Main analysis", "Database 3", 0.8, 0.7, 0.9, 1.5, 1.4, 1.6,

#' "Sensitivity analysis 1", "Database 1", 0.5, 0.4, 0.6, 2.0, 1.8, 2.2,

#' "Sensitivity analysis 1", "Database 2", 0.8, 0.7, 0.9, 1.5, 1.4, 1.6,

#' "Sensitivity analysis 1", "Database 3", 0.87, 0.78, 0.97, 0.82, 0.76, 0.87,

#' "Sensitivity analysis 2", "Database 1", 0.87, 0.78, 0.97, 0.82, 0.76, 0.87,

#' "Sensitivity analysis 2", "Database 2", 0.5, 0.4, 0.6, 2.0, 1.8, 2.2,

#' "Sensitivity analysis 2", "Database 3", 0.8, 0.7, 0.9, 1.5, 1.4, 1.6

#' )

#'

#' agreement_metrics(x, analysis_col = "Analysis", group_col = "Database")

#'

#' @export

agreement_metrics <- function(x,

                              analysis_col,

                              group_col = NULL, # subgroup, e.g, database

                              estimate_label = "HR (95% CI)",

                              smd_threshold = 1.96 # customizable threshold for SMD agreement

                              ){

  # input checks

  # check for non-positive values

  # if no smd_value present in the results table, calculate it

          )

        })

      ) |>

  }

  # calculate agreement metrics

      # significance agreement with additional categories

        # superiority

        # non-inferiority

        # inferiority

        # all other cases

        ),

      # estimate agreement

        ),

      # smd agreement using parameter

        )

      )

  # format table for display

      ) |>

  # conditional grouping

  }

  # create and format gt table

      ) |>

      )

  # apply styles for agreement columns

          } else {

          }

        )

      ) |>

          } else {

          }

        )

      )

  }

  # footnote estimate label

  # Add bold styling for headers

      ) |>

      ) |>

      ) |>

      )

}

#' 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 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 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.

#' For matched and weighted datasets, weights, cluster membership (matching only) and robust

#' variance estimates are considered in the \code{\link[survival]{survfit}} call by default.

#'

#' 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, survival probabilities are first

#' transformed using a complementary log-log transformation (\code{log(-log(1-pr(surv)))})

#' as recommended by multiple sources (Marshall, Billingham, and Bryan (2009)).

#'

#' To pool the transformed estimates across imputed datasets and time points, the 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).

#'

#' 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.

#'

#' 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

#'

#' 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

#' -  https://link.springer.com/article/10.1007/s10198-008-0129-y

#' -  https://bmcmedresmethodol.biomedcentral.com/articles/10.1186/s12874-015-0048-4

#'

#' @seealso

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

#'

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

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

#'

#' @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, 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 1,000 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)

#'

#'  # estimate and pool median survival times and Kaplan-Meier curve

#'  km_out <- km_pooling(

#'    x = wimids,

#'    surv_formula = km_fit

#'    )

#'

#'  # median survival time

#'  km_out$km_median_survival

#'

#'  # KM curve

#'  km_out$km_plot

#'

km_pooling <- function(x = NULL,

                       surv_formula = stats::as.formula(survival::Surv(fu_itt_months, death_itt) ~ treat)

                       ){

  # input checks

  # check if x is a mimids or wimids object

  # check if surv_formula is a formula

  # check weights and subclass

  # check if weights and subclass are in the data

  }

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

  # first, we fit a survival function for ith dataset

    # compute the survival function for dataset i

    ## for weighting

        )

    }

    ## for matching

        )

    }

  }

  # compute KM across all imputed (and weighted/matched) datasets --------------------

  # next, we create a list of datasets to apply the compute_km function 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 zero weight

  # will be included, too!

  # now apply the compute_km function across all imputed and matched/weighted datasets

  # this creates a list with computed survival probabilities

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

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

  # such as survival probabilities, Stef van Buuren

  # recommends to tansform survival probabilities

  # using a complementary log-log transformation

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

  # reference: Marshall, Billingham, and Bryan (2009) https://link.springer.com/article/10.1007/s10198-008-0129-y

  # reference: https://bmcmedresmethodol.biomedcentral.com/articles/10.1186/s12874-015-0048-4

  # that is, the function cloglog performs the

  # transformations on all survival probabilities

  # stored in survival_fit_list

    # convert to data.frame

      ) |>

        # 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)

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

      ) |>

    # we back-transform back and compute CI's

      ) |>

    # arrange by strata 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 {

    }

  }

  # 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 of no confidence intervals are available

    }

    # provide output table

      )

  }

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

      )

    )

}

#' Install package on demand

#'

#' @description

#' This function checks if a package is installed. If not, it

#' asks the user if they want to install it. If the user

#' agrees, it installs the package from CRAN.

#'

#' @param pkg a character string with the name of the package

#' @param quiet logical. If TRUE, suppresses messages during

#' @param ... additional arguments passed to \code{install.packages}

#'

install_on_demand <- function(pkg, quiet = FALSE, ...) {

  # internal function that checks whether package pkg is

  # in the library. If not found, it asks the user permission

  # to install.

  }

  }

}