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

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

logisregcpp <- function(data, rep = "", event = "event", covariates = "", freq = "", weight = "", offset = "", id = "", robust = 0L, firth = 0L, flic = 0L, plci = 0L, alpha = 0.05) {
    .Call(`_trtswitch_logisregcpp`, data, rep, event, covariates, freq, weight, offset, id, robust, firth, flic, plci, alpha)
}

rpsftmcpp <- function(data, stratum = "", time = "time", event = "event", treat = "treat", rx = "rx", censor_time = "censor_time", base_cov = "", low_psi = -1, hi_psi = 1, n_eval_z = 100L, treat_modifier = 1, recensor = 1L, admin_recensor_only = 1L, autoswitch = 1L, gridsearch = 0L, alpha = 0.05, ties = "efron", tol = 1.0e-6, boot = 0L, n_boot = 1000L, seed = NA_integer_) {
    .Call(`_trtswitch_rpsftmcpp`, data, stratum, time, event, treat, rx, censor_time, base_cov, low_psi, hi_psi, n_eval_z, treat_modifier, recensor, admin_recensor_only, autoswitch, gridsearch, 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)
}

fsurvci <- function(surv, sesurv, ct, z) {
    .Call(`_trtswitch_fsurvci`, surv, sesurv, ct, z)
}

#' @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.
#'
#' @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{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) {
    .Call(`_trtswitch_kmest`, data, rep, stratum, time, event, conftype, 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", robust = 0L, plci = 0L, alpha = 0.05) {
    .Call(`_trtswitch_liferegcpp`, data, rep, stratum, time, time2, event, covariates, weight, offset, id, dist, robust, plci, alpha)
}

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

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, data, stratum = "", time = "time", time2 = "", event = "event", covariates = "", weight = "", offset = "", id = "", ties = "efron", type = "schoenfeld") {
    .Call(`_trtswitch_residuals_phregcpp`, p, beta, data, stratum, time, time2, event, covariates, weight, offset, id, ties, type)
}

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", swtrt_time_upper = "", base_cov = "", conf_cov = "", low_psi = -3, hi_psi = 3, n_eval_z = 100L, strata_main_effect_only = 1L, firth = 0L, flic = 0L, recensor = 1L, admin_recensor_only = 1L, swtrt_control_only = 1L, alpha = 0.05, ties = "efron", tol = 1.0e-6, 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, swtrt_time_upper, base_cov, conf_cov, low_psi, hi_psi, n_eval_z, strata_main_effect_only, firth, flic, recensor, admin_recensor_only, swtrt_control_only, alpha, ties, tol, 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 swtrt_control_only Whether treatment switching occurred only in
#'   the control group.
#' @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 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{paneldata}: A counting process style 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{died}: Whether the patient died. 
#'     
#'     - \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{catlag}: The lagged value of \code{cattdc}.
#'     
#'     - \code{xo}: Whether the patient switched treatment. 
#'     
#'     - \code{xotime}: When the patient switched treatment.
#'     
#'     - \code{xotdc}: The time-dependent covariate for treatment 
#'       switching.
#'       
#'     - \code{xotime_upper}: The upper bound of treatment switching time.
#'     
#'     - \code{censor_time}: The administrative censoring time.
#'     
#' @author Kaifeng Lu, \email{kaifenglu@@gmail.com}
#'
#' @examples
#'
#' sim1 <- tsegestsim(
#'   n = 500, allocation1 = 2, allocation2 = 1, pbprog = 0.5, 
#'   trtlghr = -0.5, bprogsl = 0.3, shape1 = 1.8, 
#'   scale1 = 0.000025, shape2 = 1.7, scale2 = 0.000015, 
#'   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 = 0.000025, shape2 = 1.7, scale2 = 0.000015, 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, swtrt_control_only = 1L, 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, swtrt_control_only, outputRawDataset, seed)
}

tsesimpcpp <- function(data, 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, 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)
}

set_seed <- function(seed) {
    invisible(.Call(`_trtswitch_set_seed`, 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.
#' @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) {
    .Call(`_trtswitch_findInterval3`, x, v)
}

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

quantilecpp <- function(x, p) {
    .Call(`_trtswitch_quantilecpp`, x, p)
}

c_vectors_i <- function(vec1, vec2) {
    .Call(`_trtswitch_c_vectors_i`, vec1, vec2)
}

c_vectors <- function(vec1, vec2) {
    .Call(`_trtswitch_c_vectors`, vec1, vec2)
}

subset_matrix_by_row <- function(a, q) {
    .Call(`_trtswitch_subset_matrix_by_row`, a, q)
}

c_matrices <- function(a1, a2) {
    .Call(`_trtswitch_c_matrices`, a1, a2)
}

bygroup <- function(data, variables) {
    .Call(`_trtswitch_bygroup`, data, variables)
}

cholesky2 <- function(matrix, n, toler) {
    .Call(`_trtswitch_cholesky2`, matrix, n, toler)
}

chsolve2 <- function(matrix, n, y) {
    invisible(.Call(`_trtswitch_chsolve2`, matrix, n, y))
}

chinv2 <- function(matrix, n) {
    invisible(.Call(`_trtswitch_chinv2`, matrix, n))
}

invsympd <- function(matrix, n, toler) {
    .Call(`_trtswitch_invsympd`, matrix, n, toler)
}

#' @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 survplit.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{stop}: The stopping time of the resulting subrecord.
#'
#' * \code{censor}: Whether the subrecord lies strictly within a record
#'   in the input data.
#'
#' * \code{interval}: The interval number.
#'
#' @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)
}

is_sorted <- function(x) {
    .Call(`_trtswitch_is_sorted`, x)
}

house <- function(x) {
    .Call(`_trtswitch_house`, x)
}

row_house <- function(A, v) {
    invisible(.Call(`_trtswitch_row_house`, A, v))
}

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