R/RcppExports.R

In trtswitch: Treatment Switching

# Generated by using Rcpp::compileAttributes() -> do not edit by hand
# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393

ipcwcpp <- function(data, id = "id", stratum = "", tstart = "tstart", tstop = "tstop", event = "event", treat = "treat", swtrt = "swtrt", swtrt_time = "swtrt_time", base_cov = "", numerator = "", denominator = "", logistic_switching_model = 0L, strata_main_effect_only = 1L, firth = 0L, flic = 0L, ns_df = 3L, stabilized_weights = 1L, trunc = 0, trunc_upper_only = 1L, swtrt_control_only = 1L, alpha = 0.05, ties = "efron", boot = 0L, n_boot = 1000L, seed = NA_integer_) {
    .Call(`_trtswitch_ipcwcpp`, data, id, stratum, tstart, tstop, event, treat, swtrt, swtrt_time, base_cov, numerator, denominator, logistic_switching_model, strata_main_effect_only, firth, flic, ns_df, stabilized_weights, trunc, trunc_upper_only, swtrt_control_only, alpha, ties, boot, n_boot, seed)
}

ipecpp <- function(data, id = "id", stratum = "", time = "time", event = "event", treat = "treat", rx = "rx", censor_time = "censor_time", base_cov = "", aft_dist = "weibull", strata_main_effect_only = 1L, low_psi = -2, hi_psi = 2, treat_modifier = 1, recensor = 1L, admin_recensor_only = 1L, autoswitch = 1L, root_finding = "brent", alpha = 0.05, ties = "efron", tol = 1.0e-6, boot = 0L, n_boot = 1000L, seed = NA_integer_) {
    .Call(`_trtswitch_ipecpp`, data, id, stratum, time, event, treat, rx, censor_time, base_cov, aft_dist, strata_main_effect_only, low_psi, hi_psi, treat_modifier, recensor, admin_recensor_only, autoswitch, root_finding, alpha, ties, tol, boot, n_boot, seed)
}

logisregcpp <- function(data, rep = "", event = "event", covariates = "", freq = "", weight = "", offset = "", id = "", link = "logit", init = NA_real_, robust = 0L, firth = 0L, flic = 0L, plci = 0L, alpha = 0.05, maxiter = 50L, eps = 1.0e-9) {
    .Call(`_trtswitch_logisregcpp`, data, rep, event, covariates, freq, weight, offset, id, link, init, robust, firth, flic, plci, alpha, maxiter, eps)
}

msmcpp <- function(data, id = "id", stratum = "", tstart = "tstart", tstop = "tstop", event = "event", treat = "treat", swtrt = "swtrt", swtrt_time = "swtrt_time", base_cov = "", numerator = "", denominator = "", strata_main_effect_only = 1L, firth = 0L, flic = 0L, ns_df = 3L, stabilized_weights = 1L, trunc = 0, trunc_upper_only = 1L, swtrt_control_only = 1L, treat_alt_interaction = 1L, alpha = 0.05, ties = "efron", boot = 0L, n_boot = 1000L, seed = NA_integer_) {
    .Call(`_trtswitch_msmcpp`, data, id, stratum, tstart, tstop, event, treat, swtrt, swtrt_time, base_cov, numerator, denominator, strata_main_effect_only, firth, flic, ns_df, stabilized_weights, trunc, trunc_upper_only, swtrt_control_only, treat_alt_interaction, alpha, ties, boot, n_boot, seed)
}

#' @title Simulation Study to Evaluate Recensoring Rules in RPSFTM
#' 
#' @description 
#' Simulates datasets to evaluate the performance of various recensoring 
#' strategies under the Rank Preserving Structural Failure Time Model 
#' (RPSFTM) for handling treatment switching in survival analysis.
#'
#' @param nsim Number of simulated datasets.
#' @param n Number of subjects per simulation.
#' @param shape Shape parameter of the Weibull distribution for time to 
#'   death.
#' @param scale Scale parameter of the Weibull distribution for time to 
#'   death in the control group.
#' @param gamma Rate parameter of the exponential distribution for random 
#'   dropouts in the control group.
#' @param tfmin Minimum planned follow-up time (in days).
#' @param tfmax Maximum planned follow-up time (in days).
#' @param psi Log time ratio of death time for control vs experimental 
#'   treatment.
#' @param omega Log time ratio of dropout time for control vs experimental 
#'   treatment.
#' @param pswitch Probability of treatment switching at disease progression.
#' @param a Shape parameter 1 of the Beta distribution for time to disease 
#'   progression as a fraction of time to death.
#' @param b Shape parameter 2 of the Beta distribution for time to disease 
#'   progression.
#' @param low_psi Lower bound for the search interval of the causal 
#'   parameter \eqn{\psi}.
#' @param hi_psi Upper bound for the search interval of the causal 
#'   parameter \eqn{\psi}.
#' @param treat_modifier Sensitivity parameter modifying the constant 
#'   treatment effect assumption.
#' @param recensor_type Type of recensoring to apply:
#'   \itemize{
#'     \item 0: No recensoring
#'     \item 1: Recensor all control-arm subjects
#'     \item 2: Recensor only switchers in the control arm
#'     \item 3: Recensor only control-arm switchers whose counterfactual 
#'              survival exceeds the planned follow-up time
#'   }
#' @param admin_recensor_only Logical. If \code{TRUE}, recensoring is 
#'   applied only to administrative censoring times. 
#'   If \code{FALSE}, it is also applied to dropout times.
#' @param autoswitch Logical. If \code{TRUE}, disables recensoring in arms 
#'   without any treatment switching.
#' @param alpha Significance level for confidence interval calculation 
#'   (default is 0.05).
#' @param ties Method for handling tied event times in the Cox model. 
#'   Options are \code{"efron"} (default) or \code{"breslow"}.
#' @param tol Convergence tolerance for root-finding in estimation of 
#'   \eqn{\psi}.
#' @param boot Logical. If \code{TRUE}, bootstrap is used to estimate 
#'   the confidence interval for the hazard ratio. If \code{FALSE}, 
#'   the confidence interval is matched to the log-rank p-value.
#' @param n_boot Number of bootstrap samples, used only if 
#'   \code{boot = TRUE}.
#' @param seed Optional. Random seed for reproducibility. If not provided, 
#'   the global seed is used.
#'
#' @return A data frame summarizing the simulation results, including:
#' \itemize{
#'   \item \code{recensor_type}, \code{admin_recensor_only}: Settings 
#'         used in the simulation.
#'   \item Event rates: \code{p_event_1}, \code{p_dropout_1}, 
#'         \code{p_admin_censor_1}, \code{p_event_0}, 
#'         \code{p_dropout_0}, \code{p_admin_censor_0}.
#'   \item Progression and switching: \code{p_pd_0}, \code{p_swtrt_0},
#'         \code{p_recensored_0}.
#'   \item Causal parameter (\eqn{\psi}) estimates: \code{psi}, 
#'         \code{psi_est}, \code{psi_bias}, 
#'         \code{psi_se}, \code{psi_mse}.
#'   \item Log hazard ratio estimates: \code{loghr}, \code{loghr_est}, 
#'         \code{loghr_se}, \code{loghr_mse}.
#'   \item Hazard ratio metrics: \code{hr}, \code{hr_est} (geometric mean), 
#'         \code{hr_pctbias} (percent bias).
#'   \item Standard errors of log hazard ratio: \code{loghr_se_cox}, 
#'         \code{loghr_se_lr}, \code{loghr_se_boot}.
#'   \item Coverage probabilities: \code{hr_ci_cover_cox}, 
#'         \code{hr_ci_cover_lr}, \code{hr_ci_cover_boot}.
#' }
#'
#' @author Kaifeng Lu, \email{kaifenglu@@gmail.com}
#'
#' @examples
#' \donttest{
#' result <- recensor_sim_rpsftm(
#'   nsim = 10, n = 400, shape = 1.5, scale = exp(6.3169),
#'   gamma = 0.001, tfmin = 407.5, tfmax = 407.5,
#'   psi = log(0.5) / 1.5, omega = log(1), pswitch = 0.7,
#'   a = 2, b = 4, low_psi = -5, hi_psi = 5,
#'   treat_modifier = 1, recensor_type = 1,
#'   admin_recensor_only = TRUE, autoswitch = TRUE,
#'   alpha = 0.05, tol = 1e-6, boot = TRUE,
#'   n_boot = 10, seed = 314159)
#' }
#'
#' @export
recensor_sim_rpsftm <- function(nsim = NA_integer_, n = NA_integer_, shape = NA_real_, scale = NA_real_, gamma = NA_real_, tfmin = NA_real_, tfmax = NA_real_, psi = NA_real_, omega = NA_real_, pswitch = NA_real_, a = NA_real_, b = NA_real_, low_psi = -1, hi_psi = 1, treat_modifier = 1, recensor_type = 1L, admin_recensor_only = 1L, autoswitch = 1L, alpha = 0.05, ties = "efron", tol = 1.0e-6, boot = 1L, n_boot = 1000L, seed = NA_integer_) {
    .Call(`_trtswitch_recensor_sim_rpsftm`, nsim, n, shape, scale, gamma, tfmin, tfmax, psi, omega, pswitch, a, b, low_psi, hi_psi, treat_modifier, recensor_type, admin_recensor_only, autoswitch, alpha, ties, tol, boot, n_boot, seed)
}

rpsftmcpp <- function(data, id = "id", stratum = "", time = "time", event = "event", treat = "treat", rx = "rx", censor_time = "censor_time", base_cov = "", psi_test = "logrank", aft_dist = "weibull", strata_main_effect_only = 1L, low_psi = -2, hi_psi = 2, n_eval_z = 101L, treat_modifier = 1, recensor = 1L, admin_recensor_only = 1L, autoswitch = 1L, gridsearch = 0L, root_finding = "brent", alpha = 0.05, ties = "efron", tol = 1.0e-6, boot = 0L, n_boot = 1000L, seed = NA_integer_) {
    .Call(`_trtswitch_rpsftmcpp`, data, id, stratum, time, event, treat, rx, censor_time, base_cov, psi_test, aft_dist, strata_main_effect_only, low_psi, hi_psi, n_eval_z, treat_modifier, recensor, admin_recensor_only, autoswitch, gridsearch, root_finding, alpha, ties, tol, boot, n_boot, seed)
}

#' @title B-Spline Design Matrix 
#' @description Computes the design matrix for B-splines based on the 
#' specified \code{knots} and evaluated at the values in \code{x}. 
#' 
#' @param knots A numeric vector specifying the positions of the knots, 
#'   including both boundary and internal knots. 
#' @param x A numeric vector of values where the B-spline functions 
#'   or their derivatives will be evaluated. The values of \code{x} 
#'   must lie within the range of the "inner" knots, i.e., between 
#'   \code{knots[ord]} and \code{knots[length(knots) - (ord - 1)]}. 
#' @param ord A positive integer indicating the order of the B-spline. 
#'   This corresponds to the number of coefficients in each piecewise 
#'   polynomial segment, where \code{ord = degree + 1}. 
#' @param derivs An integer vector specifying the order of derivatives 
#'   to be evaluated at the corresponding \code{x} values. Each value 
#'   must be between \code{0} and \code{ord - 1}, and the vector is 
#'   conceptually recycled to match the length of \code{x}. 
#'   The default is \code{0}, meaning the B-spline functions themselves 
#'   are evaluated.
#'   
#' @return A matrix with dimensions 
#' \code{c(length(x), length(knots) - ord)}. Each row corresponds 
#' to a value in \code{x} and contains the coefficients of the 
#' B-splines, or the specified derivatives, as defined by the 
#' \code{knots} and evaluated at that particular value of \code{x}. 
#' The total number of B-splines is \code{length(knots) - ord}, 
#' with each B-spline defined by a set of \code{ord} consecutive knots.
#' 
#' @author Kaifeng Lu, \email{kaifenglu@@gmail.com}
#'
#' @examples
#'
#' splineDesigncpp(knots = 1:10, x = 4:7)
#' splineDesigncpp(knots = 1:10, x = 4:7, derivs = 1)
#'
#' @export
splineDesigncpp <- function(knots = NA_real_, x = NA_real_, ord = 4L, derivs = as.integer( c(0))) {
    .Call(`_trtswitch_splineDesigncpp`, knots, x, ord, derivs)
}

#' @title B-Spline Basis for Polynomial Splines 
#' @description Computes the B-spline basis matrix for a given polynomial 
#' spline. 
#' 
#' @param x A numeric vector representing the predictor variable. 
#' @param df Degrees of freedom, specifying the number of columns in the 
#'   basis matrix. If \code{df} is provided, the function automatically 
#'   selects \code{df - degree - intercept} internal knots based on 
#'   appropriate quantiles of \code{x}, ignoring any missing values. 
#' @param knots A numeric vector specifying the internal breakpoints 
#'   that define the spline. If not provided, \code{df} must be specified. 
#' @param degree An integer specifying the degree of the piecewise 
#'   polynomial. The default value is \code{3}, which corresponds to 
#'   cubic splines. 
#' @param intercept A logical value indicating whether to include an 
#'   intercept in the basis. The default is \code{FALSE}. 
#' @param boundary_knots A numeric vector of length 2 specifying the 
#'   boundary points where the B-spline basis should be anchored. 
#'   If not supplied, the default is the range of non-missing values 
#'   in \code{x}. 
#' @param warn_outside A logical value indicating whether a warning 
#'   should be issued if any values of \code{x} fall outside the 
#'   specified boundary knots. 
#'
#' @return A matrix with dimensions \code{c(length(x), df)}. If 
#' \code{df} is provided, the matrix will have \code{df} columns. 
#' Alternatively, if \code{knots} are supplied, the number of columns 
#' will be \code{length(knots) + degree + intercept}. The matrix 
#' contains attributes that correspond to the arguments 
#' passed to the \code{bscpp} function. 
#' 
#' @author Kaifeng Lu, \email{kaifenglu@@gmail.com}
#'
#' @examples
#'
#' bscpp(women$height, df = 5)
#'
#' @export
bscpp <- function(x = NA_real_, df = NA_integer_, knots = NA_real_, degree = 3L, intercept = 0L, boundary_knots = NA_real_, warn_outside = 1L) {
    .Call(`_trtswitch_bscpp`, x, df, knots, degree, intercept, boundary_knots, warn_outside)
}

#' @title Natural Cubic Spline Basis 
#' @description Computes the B-spline basis matrix for a natural cubic 
#' spline. 
#' 
#' @param x A numeric vector representing the predictor variable. 
#'   Missing values are allowed. 
#' @param df Degrees of freedom, specifying the number of columns in 
#'   the basis matrix. If \code{df} is provided, the function selects 
#'   \code{df - 1 - intercept} internal knots based on appropriate 
#'   quantiles of \code{x}, ignoring any missing values. 
#' @param knots A numeric vector specifying the internal breakpoints 
#'   that define the spline. If provided, the number of degrees of 
#'   freedom will be determined by the length of \code{knots}. 
#' @param intercept A logical value indicating whether to include an 
#'   intercept in the basis. The default is \code{FALSE}. 
#' @param boundary_knots A numeric vector of length 2 specifying the 
#'   boundary points where the natural boundary conditions are applied 
#'   and the B-spline basis is anchored. If not supplied, the default 
#'   is the range of non-missing values in \code{x}. 
#'
#' @return A matrix with dimensions \code{c(length(x), df)}, where 
#' \code{df} is either provided directly or computed as 
#' \code{length(knots) + 1 + intercept} when \code{knots} are supplied. 
#' The matrix contains attributes that correspond to the arguments 
#' passed to the \code{nscpp} function. 
#' 
#' @author Kaifeng Lu, \email{kaifenglu@@gmail.com}
#'
#' @examples
#'
#' nscpp(women$height, df = 5)
#'
#' @export
nscpp <- function(x = NA_real_, df = NA_integer_, knots = NA_real_, intercept = 0L, boundary_knots = NA_real_) {
    .Call(`_trtswitch_nscpp`, x, df, knots, intercept, boundary_knots)
}

#' @title Brookmeyer-Crowley Confidence Interval for Quantiles of
#' Right-Censored Time-to-Event Data
#' @description Obtains the Brookmeyer-Crowley confidence
#' interval for quantiles of right-censored time-to-event data.
#'
#' @param time The vector of possibly right-censored survival times.
#' @param event The vector of event indicators.
#' @param cilevel The confidence interval level. Defaults to 0.95.
#' @param transform The transformation of the survival function to use
#'   to construct the confidence interval. Options include
#'   "linear" (alternatively "plain"), "log",
#'   "loglog" (alternatively "log-log" or "cloglog"),
#'   "asinsqrt" (alternatively "asin" or "arcsin"), and "logit".
#'   Defaults to "loglog".
#'
#' @param probs The vector of probabilities to calculate the quantiles.
#'   Defaults to c(0.25, 0.5, 0.75).
#'
#' @return A data frame containing the estimated quantile and
#' confidence interval corresponding to each specified probability.
#' It includes the following variables:
#'
#' * \code{prob}: The probability to calculate the quantile.
#'
#' * \code{quantile}: The estimated quantile.
#'
#' * \code{lower}: The lower limit of the confidence interval.
#'
#' * \code{upper}: The upper limit of the confidence interval.
#'
#' * \code{cilevel}: The confidence interval level.
#'
#' * \code{transform}: The transformation of the survival function to use
#'   to construct the confidence interval.
#'
#' @author Kaifeng Lu, \email{kaifenglu@@gmail.com}
#'
#' @examples
#'
#' survQuantile(
#'   time = c(33.7, 3.9, 10.5, 5.4, 19.5, 23.8, 7.9, 16.9, 16.6,
#'            33.7, 17.1, 7.9, 10.5, 38),
#'   event = c(0, 1, 1, 1, 1, 0, 1, 0, 0, 0, 0, 0, 1, 1),
#'   probs = c(0.25, 0.5, 0.75))
#'
#' @export
survQuantile <- function(time = NA_real_, event = NA_real_, cilevel = 0.95, transform = "loglog", probs = NA_real_) {
    .Call(`_trtswitch_survQuantile`, time, event, cilevel, transform, probs)
}

#' @title Kaplan-Meier Estimates of Survival Curve
#' @description Obtains the Kaplan-Meier estimates of the survival curve.
#'
#' @param data The input data frame that contains the following variables:
#'
#'   * \code{rep}: The replication for by-group processing.
#'
#'   * \code{stratum}: The stratum.
#'
#'   * \code{time}: The possibly right-censored survival time.
#'
#'   * \code{event}: The event indicator.
#'
#' @param rep The name(s) of the replication variable(s) in the input data.
#' @param stratum The name(s) of the stratum variable(s) in the input data.
#' @param time The name of the time variable in the input data.
#' @param event The name of the event variable in the input data.
#' @param conftype The type of the confidence interval. One of "none",
#'   "plain", "log", "log-log" (the default), or "arcsin".
#'   The arcsin option bases the intervals on asin(sqrt(survival)).
#' @param conflev The level of the two-sided confidence interval for
#'   the survival probabilities. Defaults to 0.95.
#' @param keep_censor Whether to retain the censoring time in the output
#'   data frame.
#'
#' @return A data frame with the following variables:
#'
#' * \code{size}: The number of subjects in the stratum.
#'
#' * \code{time}: The event time.
#'
#' * \code{nrisk}: The number of subjects at risk.
#'
#' * \code{nevent}: The number of subjects having the event.
#'
#' * \code{ncensor}: The number of censored subjects.
#'
#' * \code{survival}: The Kaplan-Meier estimate of the survival probability.
#'
#' * \code{stderr}: The standard error of the estimated survival
#'   probability based on the Greendwood formula.
#'
#' * \code{lower}: The lower bound of confidence interval if requested.
#'
#' * \code{upper}: The upper bound of confidence interval if requested.
#'
#' * \code{conflev}: The level of confidence interval if requested.
#'
#' * \code{conftype}: The type of confidence interval if requested.
#'
#' * \code{stratum}: The stratum.
#'
#' * \code{rep}: The replication.
#'
#' @author Kaifeng Lu, \email{kaifenglu@@gmail.com}
#'
#' @examples
#'
#' kmest(data = aml, stratum = "x", time = "time", event = "status")
#'
#' @export
kmest <- function(data, rep = "", stratum = "", time = "time", event = "event", conftype = "log-log", conflev = 0.95, keep_censor = 0L) {
    .Call(`_trtswitch_kmest`, data, rep, stratum, time, event, conftype, conflev, keep_censor)
}

#' @title Estimate of Milestone Survival Difference
#' @description Obtains the estimate of milestone survival difference
#' between two treatment groups.
#'
#' @param data The input data frame that contains the following variables:
#'
#'   * \code{rep}: The replication for by-group processing.
#'
#'   * \code{stratum}: The stratum.
#'
#'   * \code{treat}: The treatment.
#'
#'   * \code{time}: The possibly right-censored survival time.
#'
#'   * \code{event}: The event indicator.
#'
#' @param rep The name of the replication variable in the input data.
#' @param stratum The name of the stratum variable in the input data.
#' @param treat The name of the treatment variable in the input data.
#' @param time The name of the time variable in the input data.
#' @param event The name of the event variable in the input data.
#' @param milestone The milestone time at which to calculate the
#'   survival probability.
#' @param survDiffH0 The difference in milestone survival probabilities
#'   under the null hypothesis. Defaults to 0 for superiority test.
#' @param conflev The level of the two-sided confidence interval for
#'   the difference in milestone survival probabilities. Defaults to 0.95.
#'
#' @return A data frame with the following variables:
#'
#' * \code{rep}: The replication.
#'
#' * \code{milestone}: The milestone time relative to randomization.
#'
#' * \code{survDiffH0}: The difference in milestone survival probabilities
#'   under the null hypothesis.
#'
#' * \code{surv1}: The estimated milestone survival probability for
#'   the treatment group.
#'
#' * \code{surv2}: The estimated milestone survival probability for
#'   the control group.
#'
#' * \code{survDiff}: The estimated difference in milestone survival
#'   probabilities.
#'
#' * \code{vsurv1}: The variance for surv1.
#'
#' * \code{vsurv2}: The variance for surv2.
#'
#' * \code{vsurvDiff}: The variance for survDiff.
#'
#' * \code{survDiffZ}: The Z-statistic value.
#'
#' * \code{survDiffPValue}: The one-sided p-value.
#'
#' * \code{lower}: The lower bound of confidence interval.
#'
#' * \code{upper}: The upper bound of confidence interval.
#'
#' * \code{conflev}: The level of confidence interval.
#'
#' @author Kaifeng Lu, \email{kaifenglu@@gmail.com}
#'
#' @examples
#'
#' df <- kmdiff(data = rawdata, rep = "iterationNumber",
#'              stratum = "stratum", treat = "treatmentGroup",
#'              time = "timeUnderObservation", event = "event",
#'              milestone = 12)
#' head(df)
#'
#' @export
kmdiff <- function(data, rep = "", stratum = "", treat = "treat", time = "time", event = "event", milestone = NA_real_, survDiffH0 = 0, conflev = 0.95) {
    .Call(`_trtswitch_kmdiff`, data, rep, stratum, treat, time, event, milestone, survDiffH0, conflev)
}

#' @title Log-Rank Test of Survival Curve Difference
#' @description Obtains the log-rank test using the Fleming-Harrington
#' family of weights.
#'
#' @param data The input data frame that contains the following variables:
#'
#'   * \code{rep}: The replication for by-group processing.
#'
#'   * \code{stratum}: The stratum.
#'
#'   * \code{treat}: The treatment.
#'
#'   * \code{time}: The possibly right-censored survival time.
#'
#'   * \code{event}: The event indicator.
#'
#' @param rep The name of the replication variable in the input data.
#' @param stratum The name of the stratum variable in the input data.
#' @param treat The name of the treatment variable in the input data.
#' @param time The name of the time variable in the input data.
#' @param event The name of the event variable in the input data.
#' @param rho1 The first parameter of the Fleming-Harrington family of
#'   weighted log-rank test. Defaults to 0 for conventional log-rank test.
#' @param rho2 The second parameter of the Fleming-Harrington family of
#'   weighted log-rank test. Defaults to 0 for conventional log-rank test.
#'
#' @return A data frame with the following variables:
#'
#' * \code{uscore}: The numerator of the log-rank test statistic.
#'
#' * \code{vscore}: The variance of the log-rank score test statistic.
#'
#' * \code{logRankZ}: The Z-statistic value.
#'
#' * \code{logRankPValue}: The one-sided p-value.
#'
#' * \code{rho1}: The first parameter of the Fleming-Harrington weights.
#'
#' * \code{rho2}: The second parameter of the Fleming-Harrington weights.
#'
#' * \code{rep}: The replication.
#'
#' @author Kaifeng Lu, \email{kaifenglu@@gmail.com}
#'
#' @examples
#'
#' df <- lrtest(data = rawdata, rep = "iterationNumber",
#'              stratum = "stratum", treat = "treatmentGroup",
#'              time = "timeUnderObservation", event = "event",
#'              rho1 = 0.5, rho2 = 0)
#' head(df)
#'
#' @export
lrtest <- function(data, rep = "", stratum = "", treat = "treat", time = "time", event = "event", rho1 = 0, rho2 = 0) {
    .Call(`_trtswitch_lrtest`, data, rep, stratum, treat, time, event, rho1, rho2)
}

#' @title Estimate of Restricted Mean Survival Time
#' @description Obtains the estimate of restricted means survival time
#' for each stratum.
#'
#' @param data The input data frame that contains the following variables:
#'
#'   * \code{rep}: The replication for by-group processing.
#'
#'   * \code{stratum}: The stratum.
#'
#'   * \code{time}: The possibly right-censored survival time.
#'
#'   * \code{event}: The event indicator.
#'
#' @param rep The name of the replication variable in the input data.
#' @param stratum The name of the stratum variable in the input data.
#' @param time The name of the time variable in the input data.
#' @param event The name of the event variable in the input data.
#' @param milestone The milestone time at which to calculate the
#'   restricted mean survival time.
#' @param conflev The level of the two-sided confidence interval for
#'   the survival probabilities. Defaults to 0.95.
#' @param biascorrection Whether to apply bias correction for the
#'   variance estimate. Defaults to no bias correction.
#'
#' @return A data frame with the following variables:
#'
#' * \code{rep}: The replication.
#'
#' * \code{stratum}: The stratum variable.
#'
#' * \code{size}: The number of subjects in the stratum.
#'
#' * \code{milestone}: The milestone time relative to randomization.
#'
#' * \code{rmst}: The estimate of restricted mean survival time.
#'
#' * \code{stderr}: The standard error of the estimated rmst.
#'
#' * \code{lower}: The lower bound of confidence interval if requested.
#'
#' * \code{upper}: The upper bound of confidence interval if requested.
#'
#' * \code{conflev}: The level of confidence interval if requested.
#'
#' * \code{biascorrection}: Whether to apply bias correction for the
#'   variance estimate.
#'
#' @author Kaifeng Lu, \email{kaifenglu@@gmail.com}
#'
#' @examples
#'
#' rmest(data = aml, stratum = "x",
#'       time = "time", event = "status", milestone = 24)
#'
#' @export
rmest <- function(data, rep = "", stratum = "", time = "time", event = "event", milestone = NA_real_, conflev = 0.95, biascorrection = 0L) {
    .Call(`_trtswitch_rmest`, data, rep, stratum, time, event, milestone, conflev, biascorrection)
}

#' @title Estimate of Restricted Mean Survival Time Difference
#' @description Obtains the estimate of restricted mean survival time
#' difference between two treatment groups.
#'
#' @param data The input data frame that contains the following variables:
#'
#'   * \code{rep}: The replication for by-group processing.
#'
#'   * \code{stratum}: The stratum.
#'
#'   * \code{treat}: The treatment.
#'
#'   * \code{time}: The possibly right-censored survival time.
#'
#'   * \code{event}: The event indicator.
#'
#' @param rep The name of the replication variable in the input data.
#' @param stratum The name of the stratum variable in the input data.
#' @param treat The name of the treatment variable in the input data.
#' @param time The name of the time variable in the input data.
#' @param event The name of the event variable in the input data.
#' @param milestone The milestone time at which to calculate the
#'   restricted mean survival time.
#' @param rmstDiffH0 The difference in restricted mean survival times
#'   under the null hypothesis. Defaults to 0 for superiority test.
#' @param conflev The level of the two-sided confidence interval for
#'   the difference in restricted mean survival times. Defaults to 0.95.
#' @param biascorrection Whether to apply bias correction for the
#'   variance estimate of individual restricted mean survival times.
#'   Defaults to no bias correction.
#'
#' @return A data frame with the following variables:
#'
#' * \code{rep}: The replication number.
#'
#' * \code{milestone}: The milestone time relative to randomization.
#'
#' * \code{rmstDiffH0}: The difference in restricted mean survival times
#'   under the null hypothesis.
#'
#' * \code{rmst1}: The estimated restricted mean survival time for
#'   the treatment group.
#'
#' * \code{rmst2}: The estimated restricted mean survival time for
#'   the control group.
#'
#' * \code{rmstDiff}: The estimated difference in restricted mean
#'   survival times.
#'
#' * \code{vrmst1}: The variance for rmst1.
#'
#' * \code{vrmst2}: The variance for rmst2.
#'
#' * \code{vrmstDiff}: The variance for rmstDiff.
#'
#' * \code{rmstDiffZ}: The Z-statistic value.
#'
#' * \code{rmstDiffPValue}: The one-sided p-value.
#'
#' * \code{lower}: The lower bound of confidence interval.
#'
#' * \code{upper}: The upper bound of confidence interval.
#'
#' * \code{conflev}: The level of confidence interval.
#'
#' * \code{biascorrection}: Whether to apply bias correction for the
#'   variance estimate of individual restricted mean survival times.
#'
#' @author Kaifeng Lu, \email{kaifenglu@@gmail.com}
#'
#' @examples
#'
#' df <- rmdiff(data = rawdata, rep = "iterationNumber",
#'              stratum = "stratum", treat = "treatmentGroup",
#'              time = "timeUnderObservation", event = "event",
#'              milestone = 12)
#' head(df)
#'
#' @export
rmdiff <- function(data, rep = "", stratum = "", treat = "treat", time = "time", event = "event", milestone = NA_real_, rmstDiffH0 = 0, conflev = 0.95, biascorrection = 0L) {
    .Call(`_trtswitch_rmdiff`, data, rep, stratum, treat, time, event, milestone, rmstDiffH0, conflev, biascorrection)
}

liferegcpp <- function(data, rep = "", stratum = "", time = "time", time2 = "", event = "event", covariates = "", weight = "", offset = "", id = "", dist = "weibull", init = NA_real_, robust = 0L, plci = 0L, alpha = 0.05, maxiter = 50L, eps = 1.0e-9) {
    .Call(`_trtswitch_liferegcpp`, data, rep, stratum, time, time2, event, covariates, weight, offset, id, dist, init, robust, plci, alpha, maxiter, eps)
}

residuals_liferegcpp <- function(beta, vbeta, data, stratum = "", time = "time", time2 = "", event = "event", covariates = "", weight = "", offset = "", id = "", dist = "weibull", type = "response", collapse = 0L, weighted = 0L) {
    .Call(`_trtswitch_residuals_liferegcpp`, beta, vbeta, data, stratum, time, time2, event, covariates, weight, offset, id, dist, type, collapse, weighted)
}

phregcpp <- function(data, rep = "", stratum = "", time = "time", time2 = "", event = "event", covariates = "", weight = "", offset = "", id = "", ties = "efron", init = NA_real_, robust = 0L, est_basehaz = 1L, est_resid = 1L, firth = 0L, plci = 0L, alpha = 0.05, maxiter = 50L, eps = 1.0e-9) {
    .Call(`_trtswitch_phregcpp`, data, rep, stratum, time, time2, event, covariates, weight, offset, id, ties, init, robust, est_basehaz, est_resid, firth, plci, alpha, maxiter, eps)
}

survfit_phregcpp <- function(p, beta, vbeta, basehaz, newdata, covariates = "", stratum = "", offset = "", id = "", tstart = "", tstop = "", sefit = 1L, conftype = "log-log", conflev = 0.95) {
    .Call(`_trtswitch_survfit_phregcpp`, p, beta, vbeta, basehaz, newdata, covariates, stratum, offset, id, tstart, tstop, sefit, conftype, conflev)
}

residuals_phregcpp <- function(p, beta, vbeta, resmart, data, stratum = "", time = "time", time2 = "", event = "event", covariates = "", weight = "", offset = "", id = "", ties = "efron", type = "schoenfeld", collapse = 0L, weighted = 0L) {
    .Call(`_trtswitch_residuals_phregcpp`, p, beta, vbeta, resmart, data, stratum, time, time2, event, covariates, weight, offset, id, ties, type, collapse, weighted)
}

tsegestcpp <- function(data, id = "id", stratum = "", tstart = "tstart", tstop = "tstop", event = "event", treat = "treat", censor_time = "censor_time", pd = "pd", pd_time = "pd_time", swtrt = "swtrt", swtrt_time = "swtrt_time", base_cov = "", conf_cov = "", low_psi = -2, hi_psi = 2, n_eval_z = 101L, strata_main_effect_only = 1L, firth = 0L, flic = 0L, ns_df = 3L, recensor = 1L, admin_recensor_only = 1L, swtrt_control_only = 1L, gridsearch = 0L, root_finding = "brent", alpha = 0.05, ties = "efron", tol = 1.0e-6, offset = 1, boot = 1L, n_boot = 1000L, seed = NA_integer_) {
    .Call(`_trtswitch_tsegestcpp`, data, id, stratum, tstart, tstop, event, treat, censor_time, pd, pd_time, swtrt, swtrt_time, base_cov, conf_cov, low_psi, hi_psi, n_eval_z, strata_main_effect_only, firth, flic, ns_df, recensor, admin_recensor_only, swtrt_control_only, gridsearch, root_finding, alpha, ties, tol, offset, boot, n_boot, seed)
}

#' @title Simulate Survival Data for Two-Stage Estimation Method Using 
#' g-estimation
#' @description Obtains the simulated data for baseline prognosis, 
#' disease progression, treatment switching, death, and 
#' time-dependent covariates.
#'
#' @param n The total sample size for two treatment arms combined.
#' @param allocation1 The number of subjects in the active treatment group 
#'   in a randomization block.
#' @param allocation2 The number of subjects in the control group in
#'   a randomization block.
#' @param pbprog The probability of having poor prognosis at baseline.
#' @param trtlghr The treatment effect in terms of log hazard ratio.
#' @param bprogsl The poor prognosis effect in terms of log hazard ratio.
#' @param shape1 The shape parameter for the Weibull event distribution 
#'   for the first component.
#' @param scale1 The scale parameter for the Weibull event distribution 
#'   for the first component.
#' @param shape2 The shape parameter for the Weibull event distribution 
#'   for the second component.
#' @param scale2 The scale parameter for the Weibull event distribution 
#'   for the second component.
#' @param pmix The mixing probability of the first component Weibull 
#'   distribution.
#' @param admin The administrative censoring time.
#' @param pcatnotrtbprog The probability of developing metastatic disease
#'   on control treatment with poor baseline prognosis.
#' @param pcattrtbprog The probability of developing metastatic disease
#'   on active treatment with poor baseline prognosis.
#' @param pcatnotrt The probability of developing metastatic disease
#'   on control treatment with good baseline prognosis.
#' @param pcattrt The probability of developing metastatic disease
#'   on active treatment with good baseline prognosis.
#' @param catmult The impact of metastatic disease on shortening remaining 
#'   survival time.
#' @param tdxo Whether treatment crossover depends on time-dependent 
#'   covariates between disease progression and treatment switching.
#' @param ppoor The probability of switching for poor baseline prognosis
#'   with no metastatic disease.
#' @param pgood The probability of switching for good baseline prognosis
#'   with no metastatic disease.
#' @param ppoormet The probability of switching for poor baseline prognosis
#'   after developing metastatic disease.
#' @param pgoodmet The probability of switching for good baseline prognosis
#'   after developing metastatic disease.
#' @param xomult The direct effect of crossover on extending remaining 
#'   survival time.
#' @param milestone The milestone to calculate restricted mean survival 
#'   time.
#' @param outputRawDataset Whether to output the raw data set.
#' @param seed The seed to reproduce the simulation results.
#'   The seed from the environment will be used if left unspecified.
#'
#' @return A list with two data frames.
#' 
#' * \code{sumdata}: A summary data frame with the following variables:
#'
#'     - \code{simtrueconstmean}: The true control group restricted mean 
#'       survival time (RMST).
#'     
#'     - \code{simtrueconstlb}: The lower bound for control group RMST.
#'     
#'     - \code{simtrueconstub}: The upper bound for control group RMST.
#'     
#'     - \code{simtrueconstse}: The standard error for control group RMST.
#'     
#'     - \code{simtrueexpstmean}: The true experimental group restricted 
#'       mean survival time (RMST).
#'     
#'     - \code{simtrueexpstlb}: The lower bound for experimental group RMST.
#'     
#'     - \code{simtrueexpstub}: The upper bound for experimental group RMST.
#'     
#'     - \code{simtrueexpstse}: The standard error for experimental group 
#'       RMST.
#'     
#'     - \code{simtrue_coxwbprog_hr}: The treatment hazard ratio from the 
#'       Cox model adjusting for baseline prognosis.
#'     
#'     - \code{simtrue_cox_hr}: The treatment hazard ratio from the Cox 
#'       model without adjusting for baseline prognosis.
#'
#'     - \code{simtrue_aftwbprog_af}: The average acceleration factor from 
#'       the Weibull AFT model adjusting for baseline prognosis.
#'
#'     - \code{simtrue_aft_af}: The average acceleration factor from 
#'       the Weibull AFT model without adjusting for baseline prognosis.
#'
#' * \code{paneldata}: A counting process style subject-level data frame 
#'   with the following variables:
#'
#'     - \code{id}: The subject ID.
#'
#'     - \code{trtrand}: The randomized treatment arm.
#'
#'     - \code{bprog}: Whether the patient had poor baseline prognosis. 
#'     
#'     - \code{tstart}: The left end of time interval.
#'     
#'     - \code{tstop}: The right end of time interval.
#'     
#'     - \code{event}: Whether the patient died at the end of the interval. 
#'     
#'     - \code{timeOS}: The observed survival time.
#'     
#'     - \code{died}: Whether the patient died during the study. 
#'     
#'     - \code{progressed}: Whether the patient had disease progression. 
#'     
#'     - \code{timePFSobs}: The observed time of disease progression at 
#'       regular scheduled visits.
#'       
#'     - \code{progtdc}: The time-dependent covariate for progression.
#'       
#'     - \code{catevent}: Whether the patient developed metastatic disease.
#'     
#'     - \code{cattime}: When the patient developed metastatic disease.
#'     
#'     - \code{cattdc}: The time-dependent covariate for cat event.
#'     
#'     - \code{xo}: Whether the patient switched treatment. 
#'     
#'     - \code{xotime}: When the patient switched treatment.
#'     
#'     - \code{xotdc}: The time-dependent covariate for treatment 
#'       switching.
#'       
#'     - \code{censor_time}: The administrative censoring time.
#'     
#' @author Kaifeng Lu, \email{kaifenglu@@gmail.com}
#'
#' @references
#' NR Latimer, IR White, K Tilling, and U Siebert.
#' Improved two-stage estimation to adjust for treatment switching in 
#' randomised trials: g-estimation to address time-dependent confounding.
#' Statistical Methods in Medical Research. 2020;29(10):2900-2918.
#'
#' @examples
#'
#' sim1 <- tsegestsim(
#'   n = 500, allocation1 = 2, allocation2 = 1, pbprog = 0.5, 
#'   trtlghr = -0.5, bprogsl = 0.3, shape1 = 1.8, 
#'   scale1 = 360, shape2 = 1.7, scale2 = 688, 
#'   pmix = 0.5, admin = 5000, pcatnotrtbprog = 0.5, 
#'   pcattrtbprog = 0.25, pcatnotrt = 0.2, pcattrt = 0.1, 
#'   catmult = 0.5, tdxo = 1, ppoor = 0.1, pgood = 0.04, 
#'   ppoormet = 0.4, pgoodmet = 0.2, xomult = 1.4188308, 
#'   milestone = 546, outputRawDataset = 1, seed = 2000)
#'
#' @export
tsegestsim <- function(n = 500L, allocation1 = 2L, allocation2 = 1L, pbprog = 0.5, trtlghr = -0.5, bprogsl = 0.3, shape1 = 1.8, scale1 = 360, shape2 = 1.7, scale2 = 688, pmix = 0.5, admin = 5000, pcatnotrtbprog = 0.5, pcattrtbprog = 0.25, pcatnotrt = 0.2, pcattrt = 0.1, catmult = 0.5, tdxo = 1, ppoor = 0.1, pgood = 0.04, ppoormet = 0.4, pgoodmet = 0.2, xomult = 1.4188308, milestone = 546, outputRawDataset = 1L, seed = NA_integer_) {
    .Call(`_trtswitch_tsegestsim`, n, allocation1, allocation2, pbprog, trtlghr, bprogsl, shape1, scale1, shape2, scale2, pmix, admin, pcatnotrtbprog, pcattrtbprog, pcatnotrt, pcattrt, catmult, tdxo, ppoor, pgood, ppoormet, pgoodmet, xomult, milestone, outputRawDataset, seed)
}

tsesimpcpp <- function(data, id = "id", stratum = "", time = "time", event = "event", treat = "treat", censor_time = "censor_time", pd = "pd", pd_time = "pd_time", swtrt = "swtrt", swtrt_time = "swtrt_time", base_cov = "", base2_cov = "", aft_dist = "weibull", strata_main_effect_only = 1L, recensor = 1L, admin_recensor_only = 1L, swtrt_control_only = 1L, alpha = 0.05, ties = "efron", offset = 1, boot = 1L, n_boot = 1000L, seed = NA_integer_) {
    .Call(`_trtswitch_tsesimpcpp`, data, id, stratum, time, event, treat, censor_time, pd, pd_time, swtrt, swtrt_time, base_cov, base2_cov, aft_dist, strata_main_effect_only, recensor, admin_recensor_only, swtrt_control_only, alpha, ties, offset, boot, n_boot, seed)
}

#' @title Simulate Data for Treatment Switching
#' @description Simulates data for studies involving treatment switching, 
#' incorporating time-dependent confounding. The generated data can be used 
#' to evaluate methods for handling treatment switching in survival 
#' analysis.
#'
#' @param tdxo Logical indicator for timing of treatment switching:
#'   \itemize{
#'     \item 1: Treatment switching can occur at or after disease 
#'              progression.
#'     \item 0: Treatment switching is restricted to the time of disease 
#'              progression.
#'   }
#' @param coxo Logical indicator for arm-specific treatment switching:
#'   \itemize{
#'     \item 1: Treatment switching occurs only in the control arm.
#'     \item 0: Treatment switching is allowed in both arms.
#'   }
#' @param allocation1 Number of subjects in the active treatment group in
#'   a randomization block. Defaults to 1 for equal randomization.
#' @param allocation2 Number of subjects in the control group in
#'   a randomization block. Defaults to 1 for equal randomization.
#' @param p_X_1 Probability of poor baseline prognosis in the experimental 
#'   arm.
#' @param p_X_0 Probability of poor baseline prognosis in the control arm.
#' @param rate_T Baseline hazard rate for time to death.
#' @param beta1 Log hazard ratio for randomized treatment (\code{R}).
#' @param beta2 Log hazard ratio for baseline covariate (\code{X}).
#' @param gamma0 Intercept for the time-dependent covariate model 
#'   (\code{L}).
#' @param gamma1 Coefficient for lagged treatment switching (\code{Alag}) 
#'   in the \code{L} model.
#' @param gamma2 Coefficient for lagged \code{L} (\code{Llag}) in the 
#'   \code{L} model.
#' @param gamma3 Coefficient for baseline covariate (\code{X}) in the 
#'   \code{L} model.
#' @param gamma4 Coefficient for randomized treatment (\code{R}) in the 
#'   \code{L} model.
#' @param zeta0 Intercept for the disease progression model (\code{Z}).
#' @param zeta1 Coefficient for \code{L} in the \code{Z} model.
#' @param zeta2 Coefficient for baseline covariate (\code{X}) in the 
#'   \code{Z} model.
#' @param zeta3 Coefficient for randomized treatment (\code{R}) in the 
#'   \code{Z} model.
#' @param alpha0 Intercept for the treatment switching model (\code{A}).
#' @param alpha1 Coefficient for \code{L} in the \code{A} model.
#' @param alpha2 Coefficient for baseline covariate (\code{X}) in the 
#'   \code{A} model.
#' @param theta1_1 Negative log time ratio for \code{A} for the 
#'   experimental arm.
#' @param theta1_0 Negative log time ratio for \code{A} for the control arm.
#' @param theta2 Negative log time ratio for \code{L}.
#' @param rate_C Hazard rate for random (dropout) censoring.
#' @param accrualTime A vector that specifies the starting time of 
#'   piecewise Poisson enrollment time intervals. Must start with 0, 
#'   e.g., \code{c(0, 3)} breaks the time axis into 2 accrual intervals: 
#'   [0, 3) and [3, Inf).
#' @param accrualIntensity A vector of accrual intensities. One for 
#'   each accrual time interval.
#' @param followupTime	Follow-up time for a fixed follow-up design.
#' @param fixedFollowup Whether a fixed follow-up design is used. 
#'   Defaults to 0 for variable follow-up.
#' @param plannedTime The calendar time for the analysis.
#' @param days Number of days in each treatment cycle.
#' @param n Number of subjects per simulation.
#' @param NSim Number of simulated datasets.
#' @param seed Random seed for reproducibility.
#'
#' @details 
#' The simulation algorithm is adapted from Xu et al. (2022) and simulates
#' disease progression status while incorporating the multiplicative 
#' effects of both baseline and time-dependent covariates on survival time. 
#' The design options \code{tdxo} and \code{coxo} specify
#' the timing of treatment switching and the study arm eligibility 
#' for switching, respectively. Time is measured in days. 
#' 
#' In a fixed follow-up design, all subjects share the same follow-up
#' duration. In contrast, under a variable follow-up design, follow-up
#' time also depends on each subject's enrollment date.  
#' The number of treatment cycles for a subject is determined by
#' dividing the follow-up time by the number of days in each cycle. 
#' \enumerate{
#'   \item At randomization, subjects are assigned to treatment arms 
#'   using block randomization, with \code{allocation1} patients 
#'   assigned to active treatment and \code{allocation2} to control
#'   within each randomized block. A baseline covariate is 
#'   also generated for each subject:
#'   \deqn{X_i \sim \mbox{Bernoulli}(p_{X_1} R_i + p_{X_0} (1-R_i))}
#'         
#'   \item The initial survival time is drawn
#'   from an exponential distribution with hazard:
#'   \deqn{rate_T \exp(\beta_1 R_i + \beta_2 X_i)}
#'   We define the event indicator at cycle \eqn{j} as
#'   \deqn{Y_{i,j} = I(T_i \leq j\times days)}
#'         
#'   \item The initial states are set to
#'   \eqn{L_{i,0} = 0}, \eqn{Z_{i,0} = 0}, \eqn{A_{i,0} = 0},
#'   \eqn{Y_{i,0} = 0}. For each treatment cycle \eqn{j=1,\ldots,J},
#'   we set \eqn{tstart = (j-1) \times days}.
#'         
#'   \item Generate time-dependent covariates:
#'   \deqn{\mbox{logit} P(L_{i,j}=1|\mbox{history}) = 
#'   \gamma_0 + \gamma_1 A_{i,j-1} + \gamma_2 L_{i,j-1} + 
#'   \gamma_3 X_i + \gamma_4 R_i}
#'      
#'   \item If \eqn{T_i \leq j \times days}, set \eqn{tstop = T_i} and 
#'   \eqn{Y_{i,j} = 1}, which completes data generation
#'   for subject \eqn{i}.
#'           
#'   \item If \eqn{T_i > j \times days}, set \eqn{tstop = j\times days},  
#'   \eqn{Y_{i,j} = 0}, and perform the following before proceeding to 
#'   the next cycle for the subject.
#'   
#'   \item Generate disease progression status: 
#'   If \eqn{Z_{i,j-1} = 0},  
#'   \deqn{\mbox{logit} P(Z_{i,j}=1 | \mbox{history}) = \zeta_0 + 
#'   \zeta_1 L_{i,j} + \zeta_2 X_i + \zeta_3 R_i}
#'   Otherwise, set \eqn{Z_{i,j} = 1}.
#'     
#'   \item Generate alternative therapy status:     
#'   If \eqn{A_{i,j-1} = 0}, determine switching eligibility 
#'   based on design options.        
#'   If switching is allowed:
#'   \deqn{\mbox{logit} P(A_{i,j} = 1 | \mbox{history}) = \alpha_0 + 
#'   \alpha_1 L_{i,j} + \alpha_2 X_i}       
#'   If switching is now allowed, set \eqn{A_{i,j} = 0}.      
#'   If \eqn{A_{i,j-1} = 1}, set \eqn{A_{i,j} = 1} (once switched to 
#'   alternative therapy, remain on alternative therapy).
#'          
#'   \item Update survival time based on changes in alternative 
#'   therapy status and time-dependent covariates:
#'   \deqn{T_i = j\times days + (T_i - j\times days) \exp\{
#'   -(\theta_{1,1}R_i + \theta_{1,0}(1-R_i))(A_{i,j} - A_{i,j-1}) 
#'   -\theta_2 (L_{i,j} - L_{i,j-1})\}}   
#' }
#' 
#' Additional random censoring times are generated from an exponential 
#' distribution with hazard rate \eqn{rate_C}.
#' 
#' An extra record is generated when the minimum of the latent survival 
#' time, the random censoring time, and the administrative censoring time 
#' is greater than the number of regular treatment cycles times 
#' days per cycle.
#' 
#' Finally we apply the lag function so that \eqn{Z_{i,j}} and 
#' \eqn{A_{i,j}} represent the PD status and alternative therapy status 
#' at the start of cycle \eqn{j} (and thus remain appplicable for the 
#' entire cycle \eqn{j}) for subject \eqn{i}.
#' 
#' To estimate the true treatment effect in a Cox marginal 
#' structural model, one can set \eqn{\alpha_0 = -\infty}, which 
#' effectively forces \eqn{A_{i,j} = 0} (disabling treatment switching). 
#' The coefficient for the randomized treatment can then be estimated 
#' using a Cox proportional hazards model.
#'
#' @return
#' A list of data frames, each containing simulated longitudinal covariate, 
#' pd status, alternative therapy status, and event history data with the 
#' following variables:
#'
#'  * \code{id}: Subject identifier.
#'  
#'  * \code{arrivalTime}: The enrollment time for the subject.
#'  
#'  * \code{trtrand}: Randomized treatment assignment (0 = control, 
#'    1 = experimental)
#'
#'  * \code{bprog}: Baseline prognosis (0 = good, 1 = poor).
#'
#'  * \code{tpoint}: Treatment cycle index.
#'
#'  * \code{tstart}: Start day of the treatment cycle.
#'
#'  * \code{tstop}: End day of the treatment cycle.
#'
#'  * \code{L}: Time-dependent covariate at \code{tstart} predicting 
#'    survival and switching; affected by treatment switching.
#'
#'  * \code{Llag}: Lagged value of \code{L}.
#'
#'  * \code{Z}: Disease progression status at \code{tstart}.
#'
#'  * \code{A}: Treatment switching status at \code{tstart}.
#'
#'  * \code{Alag}: Lagged value of \code{A}.
#'
#'  * \code{event}: Death indicator at \code{tstop}.
#'
#'  * \code{timeOS}: Observed time to death or censoring.
#'     
#'  * \code{died}: Indicator of death by end of follow-up.
#'     
#'  * \code{progressed}: Indicator of disease progression by end of 
#'    follow-up.
#'       
#'  * \code{timePD}: Observed time to progression or censoring.
#'       
#'  * \code{xo}: Indicator for whether treatment switching occurred.
#'     
#'  * \code{xotime}: Time of treatment switching (if applicable).
#'     
#'  * \code{censor_time}: Administrative censoring time.
#'
#' @author Kaifeng Lu, \email{kaifenglu@@gmail.com}
#'
#' @references
#' 
#' Jessica G. Young, and Eric J. Tchetgen Tchetgen. 
#' Simulation from a known Cox MSM using standard parametric models 
#' for the g-formula.
#' Statistics in Medicine. 2014;33(6):1001-1014. 
#' 
#' NR Latimer, IR White, K Tilling, and U Siebert.
#' Improved two-stage estimation to adjust for treatment switching in 
#' randomised trials: g-estimation to address time-dependent confounding.
#' Statistical Methods in Medical Research. 2020;29(10):2900-2918.
#' 
#' Jing Xu, Guohui Liu, and Bingxia Wang.
#' Bias and type I error control in correcting treatment effect
#' for treatment switching using marginal structural models 
#' in Phse III oncology trials.
#' Journal of Biopharmaceutical Statistics. 2022;32(6):897-914.
#'
#' @examples
#'
#' library(dplyr)
#' 
#' simulated.data <- tssim(
#'   tdxo = 1, coxo = 1, allocation1 = 1, allocation2 = 1,
#'   p_X_1 = 0.3, p_X_0 = 0.3, 
#'   rate_T = 0.002, beta1 = -0.5, beta2 = 0.3, 
#'   gamma0 = 0.3, gamma1 = -0.9, gamma2 = 0.7, gamma3 = 1.1, gamma4 = -0.8,
#'   zeta0 = -3.5, zeta1 = 0.5, zeta2 = 0.2, zeta3 = -0.4, 
#'   alpha0 = 0.5, alpha1 = 0.5, alpha2 = 0.4, 
#'   theta1_1 = -0.4, theta1_0 = -0.4, theta2 = 0.2,
#'   rate_C = 0.0000855, accrualIntensity = 20/30, 
#'   fixedFollowup = FALSE, plannedTime = 1350, days = 30,
#'   n = 500, NSim = 100, seed = 314159)
#'   
#' simulated.data[[1]] %>% filter(id == 1)
#'
#' @export
tssim <- function(tdxo = 0L, coxo = 1L, allocation1 = 1L, allocation2 = 1L, p_X_1 = NA_real_, p_X_0 = NA_real_, rate_T = NA_real_, beta1 = NA_real_, beta2 = NA_real_, gamma0 = NA_real_, gamma1 = NA_real_, gamma2 = NA_real_, gamma3 = NA_real_, gamma4 = NA_real_, zeta0 = NA_real_, zeta1 = NA_real_, zeta2 = NA_real_, zeta3 = NA_real_, alpha0 = NA_real_, alpha1 = NA_real_, alpha2 = NA_real_, theta1_1 = NA_real_, theta1_0 = NA_real_, theta2 = NA_real_, rate_C = NA_real_, accrualTime = 0L, accrualIntensity = NA_real_, followupTime = NA_real_, fixedFollowup = 0L, plannedTime = NA_real_, days = NA_integer_, n = NA_integer_, NSim = 1000L, seed = NA_integer_) {
    .Call(`_trtswitch_tssim`, tdxo, coxo, allocation1, allocation2, p_X_1, p_X_0, rate_T, beta1, beta2, gamma0, gamma1, gamma2, gamma3, gamma4, zeta0, zeta1, zeta2, zeta3, alpha0, alpha1, alpha2, theta1_1, theta1_0, theta2, rate_C, accrualTime, accrualIntensity, followupTime, fixedFollowup, plannedTime, days, n, NSim, seed)
}

#' @title Find Interval Numbers of Indices
#' @description The implementation of \code{findInterval()} in R from
#' Advanced R by Hadley Wickham. Given a vector of non-decreasing
#' breakpoints in v, find the interval containing each element of x; i.e.,
#' if \code{i <- findInterval3(x,v)}, for each index \code{j} in \code{x},
#' \code{v[i[j]] <= x[j] < v[i[j] + 1]}, where \code{v[0] := -Inf},
#' \code{v[N+1] := +Inf}, and \code{N = length(v)}.
#'
#' @param x The numeric vector of interest.
#' @param v The vector of break points.
#' @param rightmost_closed Logical; if true, the rightmost interval, 
#'   `vec[N-1] .. vec[N]` is treated as closed if `left_open` is false, 
#'   and the leftmost interval, `vec[1] .. vec[2]` is treated as 
#'   closed if left_open is true.
#' @param all_inside Logical; if true, the returned indices are coerced
#'   into `1, ..., N-1`, i.e., `0` is mapped to `1` and 
#'   `N` is mapped to `N-1`.
#' @param left_open Logical; if true, all intervals are open at left and 
#'   closedat right. This may be useful, .e.g., in survival analysis.   
#' @return A vector of \code{length(x)} with values in \code{0:N} where
#'   \code{N = length(v)}.
#'
#' @keywords internal
#'
#' @author Kaifeng Lu, \email{kaifenglu@@gmail.com}
#'
#' @examples
#' x <- 2:18
#' v <- c(5, 10, 15) # create two bins [5,10) and [10,15)
#' cbind(x, findInterval3(x, v))
#'
#' @export
findInterval3 <- function(x, v, rightmost_closed = 0L, all_inside = 0L, left_open = 0L) {
    .Call(`_trtswitch_findInterval3`, x, v, rightmost_closed, all_inside, left_open)
}

hasVariable <- function(df, varName) {
    .Call(`_trtswitch_hasVariable`, df, varName)
}

#' @title Split a survival data set at specified cut points
#' @description For a given survival dataset and specified cut times, 
#' each record is split into multiple subrecords at each cut time. 
#' The resulting dataset is in counting process format, with each 
#' subrecord containing a start time, stop time, and event status.
#' This is adapted from the survsplit.c function from the survival package.
#'
#' @param tstart The starting time of the time interval for 
#'   counting-process data.
#' @param tstop The stopping time of the time interval for 
#'   counting-process data.
#' @param cut The vector of cut points.
#'
#' @return A data frame with the following variables:
#'
#' * \code{row}: The row number of the observation in the input data 
#'   (starting from 0).
#'
#' * \code{start}: The starting time of the resulting subrecord.
#'
#' * \code{end}: The ending time of the resulting subrecord.
#'
#' * \code{censor}: Whether the subrecord lies strictly within a record
#'   in the input data (1 for all but the last interval and 0 for the 
#'   last interval).
#'
#' * \code{interval}: The interval number derived from cut (starting 
#'   from 0 if the interval lies to the left of the first cutpoint).
#'
#' @author Kaifeng Lu, \email{kaifenglu@@gmail.com}
#'
#' @keywords internal
#'
#' @examples
#'
#' survsplit(15, 60, c(10, 30, 40))
#'
#' @export
survsplit <- function(tstart, tstop, cut) {
    .Call(`_trtswitch_survsplit`, tstart, tstop, cut)
}

#' @title QR Decomposition of a Matrix
#' @description Computes the QR decomposition of a matrix.
#'
#' @param X A numeric matrix whose QR decomposition is to be computed.
#' @param tol The tolerance for detecting linear dependencies in the
#'   columns of \code{X}.
#'
#' @details
#' This function performs Householder QR with column pivoting:
#' Given an \eqn{m}-by-\eqn{n} matrix \eqn{A} with \eqn{m \geq n},
#' the following algorithm computes \eqn{r = \textrm{rank}(A)} and
#' the factorization \eqn{Q^T A P} equal to
#' \tabular{ccccc}{
#' | \tab \eqn{R_{11}} \tab \eqn{R_{12}} \tab | \tab \eqn{r} \cr
#' | \tab 0 \tab 0 \tab | \tab \eqn{m-r} \cr
#'   \tab \eqn{r} \tab \eqn{n-r} \tab \tab
#' }
#' with \eqn{Q = H_1 \cdots H_r} and \eqn{P = P_1 \cdots P_r}.
#' The upper triangular part of \eqn{A}
#' is overwritten by the upper triangular part of \eqn{R} and
#' components \eqn{(j+1):m} of
#' the \eqn{j}th Householder vector are stored in \eqn{A((j+1):m, j)}.
#' The permutation \eqn{P} is encoded in an integer vector \code{pivot}.
#'
#' @return A list with the following components:
#'
#' * \code{qr}: A matrix with the same dimensions as \code{X}. The upper
#'   triangle contains the \code{R} of the decomposition and the lower
#'   triangle contains Householder vectors (stored in compact form).
#'
#' * \code{rank}: The rank of \code{X} as computed by the decomposition.
#'
#' * \code{pivot}: The column permutation for the pivoting strategy used
#'   during the decomposition.
#'
#' * \code{Q}: The complete \eqn{m}-by-\eqn{m} orthogonal matrix \eqn{Q}.
#'
#' * \code{R}: The complete \eqn{m}-by-\eqn{n} upper triangular
#'   matrix \eqn{R}.
#'
#' @author Kaifeng Lu, \email{kaifenglu@@gmail.com}
#'
#' @references
#' Gene N. Golub and Charles F. Van Loan.
#' Matrix Computations, second edition. Baltimore, Maryland:
#' The John Hopkins University Press, 1989, p.235.
#'
#' @examples
#'
#' hilbert <- function(n) { i <- 1:n; 1 / outer(i - 1, i, `+`) }
#' h9 <- hilbert(9)
#' qrcpp(h9)
#'
#' @export
qrcpp <- function(X, tol = 1e-12) {
    .Call(`_trtswitch_qrcpp`, X, tol)
}

match3 <- function(id1, v1, id2, v2) {
    .Call(`_trtswitch_match3`, id1, v1, id2, v2)
}

qtpwexpcpp1 <- function(p, piecewiseSurvivalTime, lambda, lowerBound, lowertail, logp) {
    .Call(`_trtswitch_qtpwexpcpp1`, p, piecewiseSurvivalTime, lambda, lowerBound, lowertail, logp)
}