Simple Survival Analysis of time-to-event data with missing data

my.survival.analysis <- function(dat, study.end.time = NULL)
{
    ## Purpose: Survival Analysis of time-to-event data with missing data. Assume two groups. 
    ## Arguments:
    ##   dat: time-to-event data that look like (variable name is fixed):
    ## 
    ##         time status  arm
    ##           9      1    1
    ##          13      1    0
    ##          13      0    1
    ## 
    ##     The meaning of each variable is as follows: 
    ##       time: the follow-up time (numeric) for right censored data.
    ##       status: the status indicator. 0 = right censored, 1 = event. 
    ##       arm: the group indicator for comparison. 0 = control group, 1 = treatment group. 
    ## 
    ##   study.end.time: truncation time (the end of the time window in which to compare two groups).
    ##                   It needs to be smaller than the default value (minimum of the largest observed time in each of the two groups)
    ## 
    ## Return: Survival Data Analysis Summary and Plots
    ## Author: Feiming Chen
    ## ________________________________________________

    library(survival)

    dat <<- dat                         # so that "cox.zph(f2)" below works. 

    fmla <- Surv(time, status) ~ arm    # formula for survival time analysis

    f <- survfit(formula = fmla, data = dat)

    print(f)

    ## Simple plot of survival curves
    plot(f, lwd = 2, col = c("red", "blue"), mark.time = TRUE, xlab = "Time", ylab = "Survival Probability")
    legend("topright", legend = names(f$strata), bg = "lightyellow", col = c("red", "blue"), title = "Groups", lwd = 2)

    cat("\n------------------------------\nLog-Rank test for comparing two survival curves:\n")
    print(survdiff(fmla, dat))

    cat("\n------------------------------\nSurvival Probability at Study End:\n")
    if (missing(study.end.time)) {
        study.end.time <- max(f$time)
        tau <- NULL
    } else  tau <- study.end.time

    ans <- summary(f, times = study.end.time, extend = TRUE, rmean = study.end.time)
    print(ans)
    ## Post-test disease risk = 1 - survival.probability.at.end

    ## ARR (Absolute Risk Reduction) - Check the direction of the risk reduction manually! 
    ARR <- diff(ans$surv)
    cat("ARR (Absolute Risk Reduction) =", round(ARR, 2), "\n")
    NNT <- ceiling(1 / ARR)
    cat("NNT (Number Needed to Treat) =", NNT, "\n")
    OR <- ans$surv[1] / (1 - ans$surv[1]) / (ans$surv[2] / (1 - ans$surv[2]))
    cat("Odds Ratio:", round(OR, 2), "\n")

    cat("\n------------------------------\nUsing Cox Regression Model:\n")
    f2 <- coxph(formula = fmla, data = dat)
    print(f2)
    cat("Hazard Ratio =", round(exp(f2$coefficients), 2), "\n")

    cat("\n------------------------------\nAssessing Proportional Hazards Assumptions:\n")
    print(cox.zph(f2))

    cat("\n------------------------------\nRMST (Restricted Mean Survival Time) Analysis\n")
    library(survRM2)
    r <- rmst2(time = dat$time, status = dat$status, arm = dat$arm, tau = tau)
    print(r)
    ## plot(r, xlab = "Time", ylab = "Survival Probability")

    f
}

if (F) {                                # Unit Test
    d <- data.frame(time = aml$time, status = aml$status, arm = ifelse(aml$x == "Maintained", 1, 0))
    f <- my.survival.analysis(d, study.end.time = 35)

    d <- survRM2::rmst2.sample.data()
    f <- my.survival.analysis(d, study.end.time = 10)

    ## cannot do it inside a function.
    ## make survival curve with cumulative events/incidences
    f <- survfit(Surv(time, status) ~ x, data = aml)
    library(survminer)
    ggsurvplot(f, data = aml, fun = "event", conf.int = TRUE, pval = TRUE, risk.table = "abs_pct", ggtheme = theme_bw(), palette = c("red", "blue"), risk.table.y.text.col = TRUE, risk.table.y.text = FALSE, ncensor.plot = TRUE)

}

Call: survfit(formula = fmla, data = dat)

       n events median 0.95LCL 0.95UCL
arm=0 12     11     23       8      NA
arm=1 11      7     31      18      NA

------------------------------
Log-Rank test for comparing two survival curves:
Call:
survdiff(formula = fmla, data = dat)

       N Observed Expected (O-E)^2/E (O-E)^2/V
arm=0 12       11     7.31      1.86       3.4
arm=1 11        7    10.69      1.27       3.4

 Chisq= 3.4  on 1 degrees of freedom, p= 0.07 

------------------------------
Survival Probability at Study End:
Call: survfit(formula = fmla, data = dat)

                arm=0 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
     35.0000       2.0000       9.0000       0.1944       0.1219       0.0569       0.6642 

                arm=1 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
      35.000        3.000        6.000        0.368        0.163        0.155        0.875 

ARR (Absolute Risk Reduction) = 0.17 
NNT (Number Needed to Treat) = 6 
Odds Ratio: 0.41 

------------------------------
Using Cox Regression Model:
Call:
coxph(formula = fmla, data = dat)

     coef exp(coef) se(coef)    z    p
arm -0.92      0.40     0.51 -1.8 0.07

Likelihood ratio test=3.4  on 1 df, p=0.066
n= 23, number of events= 18 
Hazard Ratio = 0.4 

------------------------------
Assessing Proportional Hazards Assumptions:
         chisq df    p
arm    0.00788  1 0.93
GLOBAL 0.00788  1 0.93

------------------------------
RMST (Restricted Mean Survival Time) Analysis

The truncation time: tau = 35  was specified. 

Restricted Mean Survival Time (RMST) by arm 
               Est.    se lower .95 upper .95
RMST (arm=1) 27.057 2.906    21.361    32.753
RMST (arm=0) 20.958 3.456    14.185    27.732


Restricted Mean Time Lost (RMTL) by arm 
               Est.    se lower .95 upper .95
RMTL (arm=1)  7.943 2.906     2.247    13.639
RMTL (arm=0) 14.042 3.456     7.268    20.815


Between-group contrast 
                      Est. lower .95 upper .95     p
RMST (arm=1)-(arm=0) 6.098    -2.752    14.949 0.177
RMST (arm=1)/(arm=0) 1.291     0.878     1.899 0.194
RMTL (arm=1)/(arm=0) 0.566     0.238     1.342 0.196