## ----echo = FALSE------------------------------------------------------------- library(DiagrammeR) # Define the flowchart in DOT syntax with horizontal layout flowchart <- grViz('digraph wt_rg_S_flowchart { rankdir = LR node [shape = box, style = filled, fillcolor = "#e3e3e3"] A [label = "Start: Call wt.rg.S"] B [label = "Input Validation"] C [label = "Which scheme?"] D [label = "FH: S_left^rho * (1 - S_left)^gamma"] E [label = "Schemper: S_left / Scensor_left"] F [label = "XO: S_left / Ybar"] G [label = "MB: 1 / max(S_left, Shat_tzero)"] H [label = "custom_time: w0.tau or w1.tau by t.tau"] I [label = "fh_exp1: exp(S_left^0.5 * (1 - S_left)^0.5)"] J [label = "fh_exp2: pmin(exp(wt01), wmax)"] K [label = "Return weights"] L [label = "If details=TRUE, return list; else, return weights"] M [label = "End"] A -> B -> C C -> D [label = "fh"] C -> E [label = "schemper"] C -> F [label = "XO"] C -> G [label = "MB"] C -> H [label = "custom_time"] C -> I [label = "fh_exp1"] C -> J [label = "fh_exp2"] D -> K E -> K F -> K G -> K H -> K I -> K J -> K K -> L -> M } ', width = 800, height =400) # Display the flowchart flowchart ## ----message = FALSE, warning = FALSE----------------------------------------- oldpar <- par(no.readonly = TRUE) library(survival) # to install weightedKMplots # library(devtools) # github_install("larry-leon/weightedsurv") library(weightedsurv) library(dplyr) ## ----------------------------------------------------------------------------- library(survival) df_gbsg <- gbsg df_gbsg$tte <- df_gbsg$rfstime / 30.4375 df_gbsg$event <- df_gbsg$status df_gbsg$treat <- df_gbsg$hormon df_gbsg$grade3 <- ifelse(df_gbsg$grade == "3", 1, 0) # Note that these are default names # For alternative names, for example if the # time-to-event outcome (tte) is time_months # then tte.name <- "time_months" tte.name <- "tte" event.name <- "event" treat.name <- "treat" arms <- c("treat", "control") ## ----------------------------------------------------------------------------- # Returns standard log-rank (FH(0,0) via scheme="fh" rho,gamma = 0 in scheme_params) dfcount_gbsg <- df_counting(df=df_gbsg, tte.name = tte.name, event.name = event.name, treat.name = treat.name, arms = arms, by.risk = 12, scheme = "fh", scheme_params = list(rho = 0, gamma =0), lr.digits = 4, cox.digits =3) # test default names #test <- df_counting(df=df_gbsg, by.risk = 12, scheme = "fh", scheme_params = list(rho = 0, gamma =0), lr.digits = 4, cox.digits =3) ## ----eval=FALSE--------------------------------------------------------------- # # MB weighting # #test <- df_counting(df=df_gbsg, tte.name=tte.name, event.name=event.name, treat.name=treat.name, arms=arms, scheme = "MB", scheme_params = list(mb_tstar = 12)) # # 0/1 weighting (0 for time <= t.tau, 1 for time > t.tau) # #test <- df_counting(df=df_gbsg, tte.name=tte.name, event.name=event.name, treat.name=treat.name, arms=arms, scheme = "custom_time", scheme_params = list(t.tau = 6, w0.tau = 0, w1.tau =1)) ## ----fig.width = 10, fig.height=6, warning = FALSE, message = FALSE----------- g <- plot_weight_schemes(dfcount_gbsg) print(g) ## ----fig.width = 10, fig.height=6--------------------------------------------- par(mfrow=c(1,2)) # Mine # Set ymax a little above 1.0 to allow for logrank placement in topleft # topleft is default; xmed.fraction = 0.65 positions the display of the median estimates below the HR estimates ("topright" legend [default]) # For details, please see summary of xmed.fraction in the Appendix below. plot_weighted_km(dfcount=dfcount_gbsg, conf.int=FALSE, show.logrank = TRUE, put.legend.lr = "topleft", ymax = 1.05, xmed.fraction = 0.65) title(main="GBSG (trial) data un-weighted") # Compare with survfit plot_km(df=df_gbsg, tte.name=tte.name, event.name=event.name, treat.name=treat.name) title(main="GBSG (trial) data un-weighted via survfit") ## ----fig.width = 8, fig.height=6---------------------------------------------- plot_weighted_km(dfcount=dfcount_gbsg, conf.int=TRUE, show.logrank = TRUE, ymax = 1.05) title(main="GBSG (trial) data un-weighted") ## ----------------------------------------------------------------------------- gbsg_validate <- within(rotterdam, { rfstime <- ifelse(recur == 1, rtime, dtime) t_months <- rfstime / 30.4375 time_months <- t_months status <- pmax(recur, death) ignore <- (recur == 0 & death == 1 & rtime < dtime) status2 <- ifelse(recur == 1 | ignore, recur, death) rfstime2 <- ifelse(recur == 1 | ignore, rtime, dtime) time_months2 <- rfstime2 / 30.4375 grade3 <- ifelse(grade == "3", 1, 0) treat <- hormon event <- status2 tte <- time_months id <- as.numeric(1:nrow(rotterdam)) SG0 <- ifelse(er <= 0, 0, 1) }) ## ----------------------------------------------------------------------------- df_rotterdam <- subset(gbsg_validate, nodes > 0) ## ----------------------------------------------------------------------------- library(dplyr) rotterdam2 <- df_rotterdam %>% mutate( meno = factor(meno, levels = c(0, 1), labels = c("Pre", "Post")), grade3 = factor(grade3, levels = c(0, 1), labels = c("Not Grade 3", "Grade 3")), chemo = factor(chemo, levels = c(0, 1), labels = c("No", "Yes")), er_negative = factor(ifelse(er == 0, 1, 0), labels = c("ER negative","ER positive")) ) baseline_table <- create_baseline_table( data = rotterdam2, treat_var = "treat", vars_continuous = c("age", "pgr", "er", "nodes"), vars_categorical = c("size","meno","chemo","er_negative"), show_pvalue = TRUE, show_smd = TRUE ) ## ----results = 'asis'--------------------------------------------------------- baseline_table ## ----------------------------------------------------------------------------- fit.ps <- glm(treat ~ age + meno + size + grade3 + nodes + pgr + chemo + er, data=df_rotterdam, family="binomial") pihat <- fit.ps$fitted pihat.null <- glm(treat ~ 1, family="binomial", data=df_rotterdam)$fitted wt.1 <- pihat.null / pihat wt.0 <- (1 - pihat.null) / (1 - pihat) df_rotterdam$sw.weights <- ifelse(df_rotterdam$treat == 1, wt.1, wt.0) # truncated weights df_rotterdam$sw.weights_trunc <- with(df_rotterdam, pmin(pmax(sw.weights, quantile(sw.weights, 0.05)), quantile(sw.weights, 0.95))) ## ----------------------------------------------------------------------------- dfcount_rotterdam_unwtd <- get_dfcounting(df=df_rotterdam, tte.name=tte.name, event.name=event.name, treat.name=treat.name, arms=arms, by.risk=24) dfcount_rotterdam_wtd <- get_dfcounting(df=df_rotterdam, tte.name=tte.name, event.name=event.name, treat.name=treat.name, arms=arms, by.risk=24, weight.name="sw.weights") ## ----fig.width = 10, fig.height=6--------------------------------------------- par(mfrow=c(1,2)) plot_weighted_km(dfcount=dfcount_rotterdam_unwtd, xmed.fraction = 0.65) title(main="Rotterdam un-weighted K-M curves") plot_weighted_km(dfcount=dfcount_rotterdam_wtd, xmed.fraction = 0.65) title(main="Rotterdam weighted K-M curves") ## ----fig.width = 10, fig.height=6--------------------------------------------- par(mfrow=c(1,2)) plot_weighted_km(dfcount=dfcount_gbsg, conf.int = TRUE, xmed.fraction = 0.65, show.logrank = TRUE, ymax = 1.05) title(main="GBSG (trial) data un-weighted K-M curves") plot_weighted_km(dfcount=dfcount_rotterdam_wtd, conf.int = TRUE, xmed.fraction = 0.65, show.logrank = TRUE, ymax = 1.05) title(main="Rotterdam weighted K-M curves") ## ----fig.width = 10, fig.height = 6------------------------------------------- par(mfrow=c(1,2)) temp <- plotKM.band_subgroups( df=df_rotterdam, xlabel="Months", ylabel="Survival difference", yseq_length=5, cex_Yaxis=0.7, risk_cex=0.7, tau_add=42, by.risk = 12, risk_delta=0.075, risk.pad=0.03, ymax.pad = 0.125, tte.name = tte.name, treat.name = treat.name, event.name = event.name, weight.name = "sw.weights", draws.band = 1000, qtau = 0.025, show_resamples = TRUE ) legend("topleft", c("Difference estimate", "95% pointwise CIs"), lwd=c(2,1), col=c("black"), lty=c(1,2), bty="n", cex=0.7) title(main="Rotterdam data propensity score weighted") ## ----fig.width = 10, fig.height = 10------------------------------------------ par(mfrow=c(2,2)) temp <- plotKM.band_subgroups( df=df_rotterdam, xlabel="Months", ylabel="Survival difference", yseq_length=5, cex_Yaxis=0.7, risk_cex=0.7, tau_add=42, by.risk = 12, risk_delta=0.075, risk.pad=0.03, ymax.pad = 0.125, tte.name = tte.name, treat.name = treat.name, event.name = event.name, weight.name = "sw.weights", draws.band = 1000, qtau = 0.025, show_resamples = FALSE ) legend("topleft", c("Difference estimate", "95% pointwise CIs"), lwd=c(2,1), col=c("black"), lty=c(1,2), bty="n", cex=0.7) title(main="Rotterdam data propensity score weighted") get_bands <- cumulative_rmst_bands(df = df_rotterdam, fit = temp$fit_itt, tte.name = tte.name, event.name = event.name, treat.name = treat.name , weight.name = "sw.weights", draws_sb = 1000, xlab="months", rmst_max_cex = 0.7) legend("topleft", c("Cumulative RMST estimate", "95% pointwise CIs"), lwd=c(2,1), col=c("black"), lty=c(1,2), bty="n", cex=0.7) temp <- plotKM.band_subgroups( df=df_gbsg, Maxtau = NULL, xlabel="Months", ylabel="Survival difference", yseq_length=5, cex_Yaxis=0.7, risk_cex=0.7, by.risk = 6, risk_delta=0.075, risk.pad=0.03, tte.name = tte.name, treat.name = treat.name, event.name = event.name, draws.band = 1000, qtau = 0.025, show_resamples = FALSE ) legend("topleft", c("Difference estimate", "95% pointwise CIs"), lwd=c(2,1), col=c("black"), lty=c(1,2), bty="n", cex=0.75) title(main="GBSG data") get_bands <- cumulative_rmst_bands(df = df_gbsg, fit = temp$fit_itt, tte.name = tte.name, event.name = event.name, treat.name = treat.name, draws_sb = 1000, xlab="months", rmst_max_cex = 0.75) legend("topleft", c("Cumulative RMST estimate", "95% pointwise CIs"), lwd=c(2,1), col=c("black"), lty=c(1,2), bty="n", cex=0.75) ## ----------------------------------------------------------------------------- df_gbsg_treat <- subset(df_gbsg, treat == 1) df_rotterdam_control <- subset(df_rotterdam, treat == 0) df_gbsg_treat <- df_gbsg_treat %>% mutate( size = ifelse(size <= 20, "<=20", ifelse(size > 20 & size <=50, "20-50", ">50")), ) # note: GBSG does not have "chemo" variable df_gbsg_treat <- df_gbsg_treat[,c("tte","event","treat","age","meno","size","grade3","nodes","pgr","er")] df_rotterdam_control <- df_rotterdam_control[,c("tte","event","treat","age","meno","size","grade3","nodes","pgr","er")] # cross-trial-comparison (CTC) df_CTC <- rbind(df_gbsg_treat, df_rotterdam_control) fit.ps <- glm(treat ~ age + meno + size + grade3 + nodes + pgr + er, data = df_CTC, family="binomial") pihat <- fit.ps$fitted pihat.null <- glm(treat ~ 1, family="binomial", data = df_CTC)$fitted wt.1 <- pihat.null / pihat wt.0 <- (1 - pihat.null) / (1 - pihat) df_CTC$sw.weights <- ifelse(df_CTC$treat == 1, wt.1, wt.0) # truncated weights df_CTC$sw.weights_trunc <- with(df_CTC, pmin(pmax(sw.weights, quantile(sw.weights, 0.05)), quantile(sw.weights, 0.95))) ## ----------------------------------------------------------------------------- df_CTC$er_negative <- ifelse(df_CTC$er == 0, 1, 0) baseline_table <- create_baseline_table( data = df_CTC, treat_var = "treat", vars_continuous = c("age", "pgr", "er", "nodes"), vars_categorical = c("size","meno","er_negative"), show_pvalue = TRUE, show_smd = TRUE ) ## ----results = 'asis'--------------------------------------------------------- baseline_table ## ----------------------------------------------------------------------------- dfcount_CTC_unwtd <- get_dfcounting(df = df_CTC, tte.name=tte.name, event.name=event.name, treat.name=treat.name, arms=arms, by.risk=24) dfcount_CTC_wtd <- get_dfcounting(df=df_CTC, tte.name=tte.name, event.name=event.name, treat.name=treat.name, arms=arms, by.risk=24, weight.name="sw.weights") ## ----fig.width = 10, fig.height=6--------------------------------------------- par(mfrow=c(2,2)) plot_weighted_km(dfcount=dfcount_gbsg, conf.int = TRUE, xmed.fraction = 0.65, show.logrank = TRUE, ymax = 1.2, show.med = FALSE) title(main="GBSG (trial) data un-weighted K-M curves") plot_weighted_km(dfcount=dfcount_CTC_wtd, conf.int = TRUE, xmed.fraction = 0.65, show.logrank = TRUE, ymax = 1.2, show.med = FALSE) title(main="CTC weighted K-M curves") plot_weighted_km(dfcount=dfcount_CTC_unwtd, conf.int = TRUE, xmed.fraction = 0.65, show.logrank = TRUE, ymax = 1.2, show.med = FALSE) title(main="CTC un-weighted K-M curves") ## ----fig.width = 8, fig.height=6---------------------------------------------- par(mfrow=c(1,1)) ymin <- -1 ymax <- 5 temp <- wlr_cumulative(dfcount_CTC_wtd, scheme_params = list(rho = 0, gamma = 0), scheme = "fh", return_cumulative = TRUE) with(temp, plot(time, z.score, xlab="time", ylab="cumulative log-rank Z", type="s", ylim=c(ymin,ymax)) ) abline(h = qnorm(0.975), lwd=2, col="red") abline(v = 80, lwd=1, col="blue", lty=2) ## ----fig.width = 10, fig.height = 10------------------------------------------ draws <- 1000 par(mfrow=c(2,2)) temp <- plotKM.band_subgroups( df=df_gbsg, Maxtau = NULL, xlabel="Months", ylabel="Survival difference", yseq_length=5, cex_Yaxis=0.7, risk_cex=0.7, by.risk = 6, risk_delta=0.075, risk.pad=0.03, tte.name = tte.name, treat.name = treat.name, event.name = event.name, draws.band = draws, qtau = 0.025, show_resamples = FALSE ) legend("topleft", c("Difference estimate", "95% pointwise CIs"), lwd=c(2,1), col=c("black"), lty=c(1,2), bty="n", cex=0.75) title(main="GBSG data") get_bands <- cumulative_rmst_bands(df = df_gbsg, fit = temp$fit_itt, tte.name = tte.name, event.name = event.name, treat.name = treat.name, draws_sb = draws, xlab="months", rmst_max_cex = 0.75) legend("topleft", c("Cumulative RMST estimate", "95% pointwise CIs"), lwd=c(2,1), col=c("black"), lty=c(1,2), bty="n", cex=0.75) temp <- plotKM.band_subgroups( df=df_CTC, xlabel="Months", ylabel="Survival difference", yseq_length=5, cex_Yaxis=0.7, risk_cex=0.7, tau_add=42, by.risk = 12, risk_delta=0.075, risk.pad=0.03, ymax.pad = 0.125, tte.name = tte.name, treat.name = treat.name, event.name = event.name, weight.name ="sw.weights", draws.band = draws, qtau = 0.025, show_resamples = FALSE ) legend("topleft", c("Difference estimate", "95% pointwise CIs"), lwd=c(2,1), col=c("black"), lty=c(1,2), bty="n", cex=0.7) title(main="CTC data") get_bands <- cumulative_rmst_bands(df = df_CTC, fit = temp$fit_itt, tte.name = tte.name, event.name = event.name, treat.name = treat.name , weight.name = "sw.weights", draws_sb = draws, xlab="months", rmst_max_cex = 0.7) legend("topleft", c("Cumulative RMST estimate", "95% pointwise CIs"), lwd=c(2,1), col=c("black"), lty=c(1,2), bty="n", cex=0.7) ## ----fig.width = 8, fig.height=6---------------------------------------------- par(mfrow=c(1,1)) sg_labs <- c("er == 0","er > 0", "meno == 0", "meno ==1") sg_cols <- c("blue", "brown", "green", "turquoise") # Randomly sample colors from the built-in palette # n_sg <- length(sg_labs) # set.seed(123) # for reproducibility # sg_cols <- sample(colors(), n_sg) temp <- plotKM.band_subgroups( df=df_gbsg, sg_labels = sg_labs, sg_colors = sg_cols, xlabel="Months", ylabel="Survival difference", yseq_length=5, cex_Yaxis=0.7, risk_cex=0.7, tau_add=42, by.risk=6, risk_delta=0.05, risk.pad=0.03, draws.band = 1000, tte.name = tte.name, treat.name = treat.name, event.name = event.name ) legend("topleft", c("ITT", sg_labs), lwd=2, col=c("black", sg_cols), bty="n", cex=0.75) title(main="ITT and subgroups") ## ----fig.width = 12, fig.height=6--------------------------------------------- par(mfrow=c(1,2)) temp <- plotKM.band_subgroups( df = subset(df_gbsg, er > 0), xlabel="Months", ylabel="Survival difference", yseq_length=5, cex_Yaxis=0.7, risk_cex=0.7, by.risk=6, risk_delta=0.075, risk.pad=0.03, tte.name = tte.name, treat.name = treat.name, event.name = event.name, draws.band = 1000, qtau = 0.025, show_resamples = TRUE ) legend("topleft", c("Difference estimate", "95% pointwise CIs"), lwd=c(2,1), col=c("black"), lty=c(1,2), bty="n", cex=0.75) title(main="Estrogen receptor positive (er>0) sub-population") ## ----messages = FALSE, warning = FALSE---------------------------------------- library(data.table) library(gt) dfcount <- dfcount_rotterdam_wtd get_rg <- cox_rhogamma(dfcount = dfcount, scheme = "fh", scheme_params = list(rho = 0, gamma = 0), draws = 1000, verbose = FALSE) # Note: with case-weights (here propensity-score weighting) recommend asymptotic versions of SEs res <- data.table() res$analysis <- "FH(0,0)" res$z <- get_rg$z.score res$hr <- get_rg$hr_ci_asy$hr res$lower <- get_rg$hr_ci_star$lower res$upper <- get_rg$hr_ci_star$upper # res$lower <- get_rg$hr_ci_asy$lower # res$upper <- get_rg$hr_ci_asy$upper # compare with coxph() #cat("Comparison with coxph with weighting, coxph:", dfcount$cox_results$cox_text, "\n") #cat("Mine:", get_rg$cox_text_star, "\n") res_update <- res get_rg <- cox_rhogamma(dfcount = dfcount, scheme = "fh", scheme_params = list(rho = 0, gamma = 1), draws = 1000, verbose = FALSE) res <- data.table() res$analysis <- "FH(0,1)" res$hr <- get_rg$hr_ci_asy$hr res$lower <- get_rg$hr_ci_star$lower res$upper <- get_rg$hr_ci_star$upper res$z <- get_rg$z.score res_update <- rbind(res_update, res) get_rg <- cox_rhogamma(dfcount = dfcount, scheme = "fh_exp2", draws = 1000, verbose = FALSE) res <- data.table() res$analysis <- "FH_exp2" res$hr <- get_rg$hr_ci_asy$hr res$lower <- get_rg$hr_ci_star$lower res$upper <- get_rg$hr_ci_star$upper res$z <- get_rg$z.score res_update <- rbind(res_update, res) get_rg <- cox_rhogamma(dfcount = dfcount, scheme = "MB", scheme_params = list(mb_tstar = 12), draws = 1000, verbose = FALSE) res <- data.table() res$analysis <- "MB(12)" res$hr <- get_rg$hr_ci_asy$hr res$lower <- get_rg$hr_ci_star$lower res$upper <- get_rg$hr_ci_star$upper res$z <- get_rg$z.score res_update <- rbind(res_update, res) get_rg <- cox_rhogamma(dfcount = dfcount, scheme = "custom_time", scheme_params = list(t.tau = 6, w0.tau = 0, w1.tau = 1), draws = 1000, verbose = FALSE) res <- data.table() res$analysis <- "zero_one(6)" res$hr <- get_rg$hr_ci_asy$hr res$lower <- get_rg$hr_ci_star$lower res$upper <- get_rg$hr_ci_star$upper res$z <- get_rg$z.score res_update <- rbind(res_update, res) res_update %>% gt() %>% fmt_number(decimals = 3) ## ----------------------------------------------------------------------------- dfcount_rotterdam_wtd <- get_dfcounting(df=df_rotterdam, tte.name=tte.name, event.name=event.name, treat.name=treat.name, arms=arms, by.risk=24, weight.name="sw.weights", check.seKM = TRUE, draws = 0) ## ----------------------------------------------------------------------------- dfcount_rotterdam_wtd <- get_dfcounting(df=df_rotterdam, tte.name=tte.name, event.name=event.name, treat.name=treat.name, arms=arms, by.risk=24, weight.name="sw.weights", check.seKM = TRUE, draws = 5000) ## ----------------------------------------------------------------------------- # Compare with adjustedCurves library(adjustedCurves) library(pammtools) df_rotterdam$hormon2 <- with(df_rotterdam, factor(hormon, labels=c("No","Yes"))) ps_model <- glm(hormon2 ~ age+meno+size+grade+nodes+pgr+chemo+er, data = df_rotterdam, family="binomial") iptw <- adjustedsurv(data=df_rotterdam, variable="hormon2", ev_time="time_months", event = "status", method = "iptw_km", treatment_model = ps_model, conf_int = TRUE) ## ----fig.width = 8, fig.height = 6-------------------------------------------- plot_weighted_km(dfcount=dfcount_rotterdam_wtd, conf.int = TRUE) title(main="Rotterdam weighted K-M curves") ## ----fig.width = 8, fig.height = 6-------------------------------------------- plot(iptw, conf_int = TRUE, legend.title = "Hormonal therapy") ## ----fig.width = 12, fig.height = 6------------------------------------------- # Check SEs par(mfrow=c(1,2)) df_mine <- dfcount_rotterdam_wtd df_check <- subset(iptw$adj, group == "No") yymax <- max(c(sqrt(df_mine$sig2_surv0[df_mine$idv0.check]), df_check$se)) toget <- df_mine$idv0.check tt <- df_mine$at_points[toget] yy <- sqrt(df_mine$sig2_surv0)[toget] plot(tt,yy, type="s", lty=1, col="lightgrey", lwd=4, ylim=c(0,yymax), xlab="time", ylab="SEs") with(df_check, lines(time, se, type="s", lty=2, lwd=1, col="red")) title(main="Control SEs") df_check <- subset(iptw$adj, group == "Yes") yymax <- max(c(sqrt(df_mine$sig2_surv1[df_mine$idv1.check]), df_check$se)) toget <- df_mine$idv1.check tt <- df_mine$at_points[toget] yy <- sqrt(df_mine$sig2_surv1)[toget] plot(tt,yy, type="s", lty=1, col="lightgrey", lwd=4, ylim=c(0,yymax), xlab="time", ylab="SEs") with(df_check, lines(time, se, type="s", lty=2, lwd=1, col="red")) title(main = "Treat SEs") ## ----------------------------------------------------------------------------- # dfcount_gbsg contains logrank (chisq) stats from specified scheme (default is standard logrank rho=0 and gamma=0) from 3 sources: # (1) wlr_estimates function # (2) survdiff which only provides rho options so gamma !=0 is not available # (3) z_score_calculations function (from a Cox score-test perspective) check_results(dfcount_gbsg) # Check wlr_dhat_estimates which calculates from any scheme based on counting dataset # In addition, calculates dhat's at tzero as well as correlation with log-rank temp <- wlr_dhat_estimates(dfcounting = dfcount_gbsg, scheme_params = list(rho = 0, gamma = 0), scheme = "fh", tzero = 12) temp$z.score^2 ## ----------------------------------------------------------------------------- res <- df_counting(df=df_gbsg, tte.name=tte.name, event.name=event.name, treat.name=treat.name, arms=arms, strata.name="grade") ## ----------------------------------------------------------------------------- with(res,lr_stratified^2/sig2_lr_stratified) # survdiff stratified res$logrank_results$chisq # Cox score stratified with(res,z.score_stratified^2) ## ----------------------------------------------------------------------------- with(res,lr^2/sig2_lr) with(res,z.score^2) # Restore when done par(oldpar)