Purpose

This workflow analyzes the time until a doctoral student reaches a program milestone. It demonstrates right censoring, Kaplan-Meier estimation, group comparison, Cox proportional hazards regression, model diagnostics, adjusted survival curves, and results reporting.

The included data are entirely synthetic. Run the complete example before adapting it. Define the event, time origin, censoring rules, and unit of analysis in writing before fitting any survival model.

Survival analysis foundations

  • Time origin: when follow-up begins for every participant.
  • Event: the outcome whose timing is modeled.
  • Observed time: time from origin to the event or censoring.
  • Right censoring: follow-up ended before the event was observed.
  • Survival function: probability of remaining event-free beyond time t.
  • Hazard: instantaneous event rate among participants still at risk.
  • Hazard ratio: relative hazard associated with a predictor, conditional on the model.

A hazard ratio is not a risk ratio, odds ratio, or difference in median survival time.

User settings

input_file <- "Data/synthetic_survival_data.csv"
generate_synthetic_example <- TRUE
id_variable <- "participant_id"
time_variable <- "time_months"
event_variable <- "event"
group_variable <- "support_program"
analysis_horizon <- 36
landmark_times <- c(0, 6, 12, 18, 24, 30, 36)

Simulate and import the example

if (generate_synthetic_example) {
  synthetic_data <- simulate_survival_data()
  readr::write_csv(synthetic_data, input_file, na = "")
}

raw_data <- readr::read_csv(input_file, show_col_types = FALSE) %>%
  mutate(
    support_program = factor(
      support_program,
      levels = c("Standard support", "Enhanced support")
    ),
    part_time = factor(
      part_time,
      levels = c("Full-time", "Part-time")
    )
  )

glimpse(raw_data)
Rows: 450
Columns: 8
$ participant_id    <chr> "D001", "D002", "D003", "D004", "D005", "D006", "D00~
$ time_months       <dbl> 16.73, 25.13, 9.49, 12.26, 4.18, 15.88, 9.81, 5.73, ~
$ event             <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 0, 1, 1, 1~
$ support_program   <fct> Enhanced support, Enhanced support, Enhanced support~
$ part_time         <fct> Full-time, Full-time, Part-time, Full-time, Full-tim~
$ age               <dbl> 30, 35, 37, 23, 34, 26, 48, 32, 28, 44, 29, 29, 51, ~
$ baseline_progress <dbl> 61, 61, 70, 84, 45, 62, 65, 70, 47, 65, 73, 37, 55, ~
$ advisor_support   <dbl> 28, 33, 31, 31, 30, 23, 25, 33, 20, 36, 29, 29, 30, ~

Data audit

Event and follow-up summary

event_summary <- raw_data %>%
  summarise(
    participants = n(),
    events = sum(.data[[event_variable]] == 1, na.rm = TRUE),
    censored = sum(.data[[event_variable]] == 0, na.rm = TRUE),
    event_percent = round(100 * mean(.data[[event_variable]] == 1, na.rm = TRUE), 1),
    median_follow_up = round(median(.data[[time_variable]], na.rm = TRUE), 1),
    maximum_follow_up = max(.data[[time_variable]], na.rm = TRUE)
  )

knitr::kable(event_summary)
participants events censored event_percent median_follow_up maximum_follow_up
450 291 159 64.7 16.5 36

Missingness

missingness <- tibble(
  variable = names(raw_data),
  missing_n = map_int(raw_data, ~ sum(is.na(.x))),
  missing_percent = round(100 * map_dbl(raw_data, ~ mean(is.na(.x))), 1)
) %>%
  arrange(desc(missing_n), variable)

knitr::kable(missingness)
variable missing_n missing_percent
baseline_progress 8 1.8
advisor_support 7 1.6
age 0 0.0
event 0 0.0
part_time 0 0.0
participant_id 0 0.0
support_program 0 0.0
time_months 0 0.0

This example uses complete cases for the adjusted Cox model so the effect of missing predictors is visible. For substantive research, justify complete-case analysis or use a suitable multiple-imputation strategy that includes the event indicator and survival information.

Validity checks

stopifnot(!anyDuplicated(raw_data[[id_variable]]))
stopifnot(all(raw_data[[time_variable]] >= 0, na.rm = TRUE))
stopifnot(all(raw_data[[event_variable]] %in% c(0, 1)))
stopifnot(all(raw_data[[time_variable]] <= analysis_horizon, na.rm = TRUE))

analysis_data <- raw_data %>%
  drop_na(time_months, event, support_program, part_time, age, baseline_progress, advisor_support)

analysis_flow <- tibble(
  stage = c("Imported records", "Complete cases for adjusted model", "Records excluded for missing predictors"),
  n = c(nrow(raw_data), nrow(analysis_data), nrow(raw_data) - nrow(analysis_data))
)

knitr::kable(analysis_flow)
stage n
Imported records 450
Complete cases for adjusted model 435
Records excluded for missing predictors 15

Visualize observed follow-up

follow_up_plot_data <- raw_data %>%
  arrange(desc(time_months), event) %>%
  mutate(display_order = row_number())

ggplot(follow_up_plot_data, aes(x = time_months, y = display_order, color = factor(event))) +
  geom_segment(aes(x = 0, xend = time_months, yend = display_order), color = "#D8C9D0", linewidth = .45) +
  geom_point(size = 1.8) +
  scale_color_manual(
    values = c("0" = dsh_colors[["mauve"]], "1" = dsh_colors[["coral"]]),
    labels = c("Censored", "Event"),
    name = NULL
  ) +
  scale_x_continuous(limits = c(0, analysis_horizon), breaks = landmark_times) +
  labs(
    title = "Observed follow-up for each participant",
    subtitle = "Dots mark the event or last observed follow-up",
    x = "Months since study entry",
    y = "Participants ordered by follow-up time"
  ) +
  theme(axis.text.y = element_blank(), axis.ticks.y = element_blank())

Kaplan-Meier estimates

Overall estimate

survival_object <- with(raw_data, Surv(time_months, event))
km_overall <- survfit(survival_object ~ 1, data = raw_data, conf.type = "log-log")
km_overall
Call: survfit(formula = survival_object ~ 1, data = raw_data, conf.type = "log-log")

       n events median 0.95LCL 0.95UCL
[1,] 450    291   18.9    17.4    20.8

By support program

km_group <- survfit(Surv(time_months, event) ~ support_program, data = raw_data, conf.type = "log-log")
km_group
Call: survfit(formula = Surv(time_months, event) ~ support_program, 
    data = raw_data, conf.type = "log-log")

                                   n events median 0.95LCL 0.95UCL
support_program=Standard support 237    138   21.3    19.4    24.2
support_program=Enhanced support 213    153   16.4    13.6    18.2

Modern survival curves

km_tidy <- broom::tidy(km_group) %>%
  mutate(
    group = sub("support_program=", "", strata),
    group = factor(group, levels = c("Standard support", "Enhanced support"))
  )

ggplot(km_tidy, aes(x = time, y = estimate, color = group)) +
  geom_step(linewidth = 1.1) +
  geom_step(aes(y = conf.low), linewidth = .55, linetype = "dashed", alpha = .75) +
  geom_step(aes(y = conf.high), linewidth = .55, linetype = "dashed", alpha = .75) +
  scale_color_manual(values = c(dsh_colors[["plum"]], dsh_colors[["coral"]]), name = "Support condition") +
  scale_x_continuous(limits = c(0, analysis_horizon), breaks = landmark_times) +
  scale_y_continuous(limits = c(0, 1), labels = scales::label_percent()) +
  labs(
    title = "Kaplan-Meier event-free survival",
    subtitle = "Dashed lines show pointwise 95% confidence intervals",
    x = "Months since study entry",
    y = "Probability of remaining event-free"
  )

Cumulative incidence

ggplot(km_tidy, aes(x = time, y = 1 - estimate, color = group)) +
  geom_step(linewidth = 1.1) +
  scale_color_manual(values = c(dsh_colors[["plum"]], dsh_colors[["coral"]]), name = "Support condition") +
  scale_x_continuous(limits = c(0, analysis_horizon), breaks = landmark_times) +
  scale_y_continuous(limits = c(0, 1), labels = scales::label_percent()) +
  labs(
    title = "Cumulative incidence of the milestone",
    subtitle = "This is 1 minus the Kaplan-Meier estimate in a single-event setting",
    x = "Months since study entry",
    y = "Estimated cumulative incidence"
  )

Number at risk

risk_summary <- summary(km_group, times = landmark_times, extend = TRUE)
risk_table <- tibble(
  group = sub("support_program=", "", risk_summary$strata),
  month = risk_summary$time,
  at_risk = risk_summary$n.risk,
  events = risk_summary$n.event
) %>%
  pivot_wider(names_from = month, values_from = c(at_risk, events), names_glue = "{.value}_{month}m")

knitr::kable(risk_table)
group at_risk_0m at_risk_6m at_risk_12m at_risk_18m at_risk_24m at_risk_30m at_risk_36m events_0m events_6m events_12m events_18m events_24m events_30m events_36m
Standard support 237 214 175 113 62 22 11 0 23 39 29 27 15 5
Enhanced support 213 180 130 78 40 18 7 0 33 50 33 20 16 1

Median survival

median_table <- as.data.frame(summary(km_group)$table) %>%
  rownames_to_column("group") %>%
  transmute(
    group = sub("support_program=", "", group),
    records = records,
    events = events,
    median_months = median,
    median_conf_low = `0.95LCL`,
    median_conf_high = `0.95UCL`
  )

knitr::kable(median_table)
group records events median_months median_conf_low median_conf_high
Standard support 237 138 21.32 19.44 24.24
Enhanced support 213 153 16.41 13.56 18.18

If the survival curve never crosses 0.50, the median is not estimable within follow-up. Report a milestone-specific survival probability or restricted mean survival time instead of inventing a median.

Log-rank comparison

log_rank <- survdiff(Surv(time_months, event) ~ support_program, data = raw_data, rho = 0)
log_rank
Call:
survdiff(formula = Surv(time_months, event) ~ support_program, 
    data = raw_data, rho = 0)

                                   N Observed Expected (O-E)^2/E (O-E)^2/V
support_program=Standard support 237      138      166      4.82      11.3
support_program=Enhanced support 213      153      125      6.42      11.3

 Chisq= 11.3  on 1 degrees of freedom, p= 8e-04 
log_rank_p <- pchisq(log_rank$chisq, df = length(log_rank$n) - 1, lower.tail = FALSE)

The log-rank test compares the full unadjusted survival curves. It does not estimate an adjusted effect and should be interpreted alongside effect estimates and confidence intervals.

Cox proportional hazards models

Unadjusted model

cox_unadjusted <- coxph(
  Surv(time_months, event) ~ support_program,
  data = raw_data,
  ties = "efron",
  x = TRUE
)

broom::tidy(cox_unadjusted, exponentiate = TRUE, conf.int = TRUE) %>%
  select(term, estimate, conf.low, conf.high, p.value) %>%
  knitr::kable(digits = 3)
term estimate conf.low conf.high p.value
support_programEnhanced support 1.481 1.176 1.864 0.001

Adjusted model

cox_adjusted <- coxph(
  Surv(time_months, event) ~ support_program + part_time +
    scale(age) + scale(baseline_progress) + scale(advisor_support),
  data = analysis_data,
  ties = "efron",
  x = TRUE,
  y = TRUE
)

cox_table <- broom::tidy(cox_adjusted, exponentiate = TRUE, conf.int = TRUE) %>%
  transmute(
    predictor = term,
    hazard_ratio = estimate,
    conf_low = conf.low,
    conf_high = conf.high,
    p_value = p.value
  )

knitr::kable(cox_table, digits = 3)
predictor hazard_ratio conf_low conf_high p_value
support_programEnhanced support 1.347 1.054 1.722 0.017
part_timePart-time 0.794 0.615 1.024 0.076
scale(age) 0.899 0.799 1.012 0.077
scale(baseline_progress) 1.125 0.998 1.268 0.054
scale(advisor_support) 1.243 1.096 1.408 0.001 
readr::write_csv(cox_table, "Output/cox_model_results.csv")

Continuous predictors are standardized inside the formula, so their hazard ratios represent a one-standard-deviation increase. Positive coefficients indicate a higher event hazard and typically a shorter time to the event; negative coefficients indicate a lower event hazard and typically a longer time to the event.

Model performance

model_summary <- summary(cox_adjusted)
model_performance <- tibble(
  analysis_n = cox_adjusted$n,
  events = cox_adjusted$nevent,
  concordance = unname(model_summary$concordance[1]),
  concordance_se = unname(model_summary$concordance[2]),
  likelihood_ratio_chisq = unname(model_summary$logtest[1]),
  likelihood_ratio_df = unname(model_summary$logtest[2]),
  likelihood_ratio_p = unname(model_summary$logtest[3])
)

knitr::kable(model_performance, digits = 3)
analysis_n events concordance concordance_se likelihood_ratio_chisq likelihood_ratio_df likelihood_ratio_p
435 281 0.598 0.018 33.122 5 0 
readr::write_csv(model_performance, "Output/model_performance.csv")

Adjusted survival curves

reference_profiles <- tibble(
  support_program = factor(
    c("Standard support", "Enhanced support"),
    levels = levels(analysis_data$support_program)
  ),
  part_time = factor("Full-time", levels = levels(analysis_data$part_time)),
  age = mean(analysis_data$age),
  baseline_progress = mean(analysis_data$baseline_progress),
  advisor_support = mean(analysis_data$advisor_support)
)

adjusted_fit <- survfit(cox_adjusted, newdata = reference_profiles)
adjusted_matrix <- as.data.frame(adjusted_fit$surv)
names(adjusted_matrix) <- levels(reference_profiles$support_program)

adjusted_tidy <- bind_cols(time = adjusted_fit$time, adjusted_matrix) %>%
  pivot_longer(-time, names_to = "group", values_to = "survival") %>%
  mutate(group = factor(group, levels = c("Standard support", "Enhanced support")))

ggplot(adjusted_tidy, aes(x = time, y = survival, color = group)) +
  geom_step(linewidth = 1.15) +
  scale_color_manual(values = c(dsh_colors[["plum"]], dsh_colors[["coral"]]), name = "Support condition") +
  scale_x_continuous(limits = c(0, analysis_horizon), breaks = landmark_times) +
  scale_y_continuous(limits = c(0, 1), labels = scales::label_percent()) +
  labs(
    title = "Model-adjusted event-free survival",
    subtitle = "Profiles are full-time students at the sample means of continuous covariates",
    x = "Months since study entry",
    y = "Predicted probability of remaining event-free"
  )

Model diagnostics

Proportional hazards

ph_test <- cox.zph(cox_adjusted, transform = "km")
ph_table <- as.data.frame(ph_test$table) %>%
  rownames_to_column("term")

knitr::kable(ph_table, digits = 3)
term chisq df p
support_program 1.305 1 0.253
part_time 0.149 1 0.700
scale(age) 1.981 1 0.159
scale(baseline_progress) 1.740 1 0.187
scale(advisor_support) 0.224 1 0.636
GLOBAL 6.186 5 0.289 
readr::write_csv(ph_table, "Output/proportional_hazards_test.csv")

A small p-value suggests that a coefficient may vary over follow-up. Evaluate the global test, term-specific tests, residual patterns, study design, and statistical power together. If proportional hazards is not plausible, consider a time interaction, stratification, or a different estimand.

Scaled Schoenfeld residuals

schoenfeld_data <- as.data.frame(ph_test$y) %>%
  mutate(event_time = ph_test$time, .before = 1) %>%
  pivot_longer(-event_time, names_to = "term", values_to = "residual")

ggplot(schoenfeld_data, aes(x = event_time, y = residual)) +
  geom_hline(yintercept = 0, color = "#B8A7B6", linewidth = .5) +
  geom_point(color = dsh_colors[["plum"]], alpha = .35, size = 1) +
  geom_smooth(color = dsh_colors[["coral"]], se = TRUE, linewidth = .8) +
  facet_wrap(~ term, scales = "free_y") +
  labs(
    title = "Scaled Schoenfeld residuals",
    subtitle = "A persistent time trend can indicate nonproportional hazards",
    x = "Event time in months",
    y = "Scaled residual"
  )

Functional form

base_for_martingale <- coxph(
  Surv(time_months, event) ~ support_program + part_time,
  data = analysis_data,
  ties = "efron"
)

functional_form_data <- analysis_data %>%
  mutate(martingale = residuals(base_for_martingale, type = "martingale")) %>%
  select(martingale, age, baseline_progress, advisor_support) %>%
  pivot_longer(-martingale, names_to = "predictor", values_to = "value")

ggplot(functional_form_data, aes(x = value, y = martingale)) +
  geom_hline(yintercept = 0, color = "#B8A7B6", linewidth = .5) +
  geom_point(color = dsh_colors[["plum"]], alpha = .3, size = 1) +
  geom_smooth(color = dsh_colors[["coral"]], se = TRUE, linewidth = .9) +
  facet_wrap(~ predictor, scales = "free_x") +
  labs(
    title = "Review continuous-predictor functional form",
    subtitle = "Strong curvature suggests that a linear Cox term may be inadequate",
    x = "Predictor value",
    y = "Martingale residual"
  )

Splines, transformations, or prespecified categories may address nonlinearity. Choose them using subject-matter knowledge and document any data-informed decisions.

Influential observations

dfbeta_matrix <- as.matrix(residuals(cox_adjusted, type = "dfbeta"))
influence_data <- tibble(
  participant_id = analysis_data$participant_id,
  maximum_absolute_dfbeta = apply(abs(dfbeta_matrix), 1, max)
)

influence_threshold <- 2 / sqrt(nrow(analysis_data))

ggplot(influence_data, aes(x = reorder(participant_id, maximum_absolute_dfbeta), y = maximum_absolute_dfbeta)) +
  geom_point(color = dsh_colors[["plum"]], alpha = .65) +
  geom_hline(yintercept = influence_threshold, color = dsh_colors[["coral"]], linetype = "dashed") +
  coord_flip() +
  labs(
    title = "Influence review using DFBETA values",
    subtitle = "The dashed line is a screening heuristic, not an automatic deletion rule",
    x = "Participant",
    y = "Maximum absolute DFBETA"
  ) +
  theme(axis.text.y = element_blank(), axis.ticks.y = element_blank())

influence_review <- influence_data %>%
  filter(maximum_absolute_dfbeta > influence_threshold) %>%
  arrange(desc(maximum_absolute_dfbeta))

readr::write_csv(influence_review, "Output/influence_review.csv")
knitr::kable(head(influence_review, 10), digits = 3)
participant_id maximum_absolute_dfbeta

Investigate influential records for data errors, unusual covariate patterns, and sensitivity. Do not remove a valid observation solely because it is influential.

Event-based planning

For a balanced binary predictor, Schoenfeld’s approximation estimates the number of events needed to detect a target hazard ratio under proportional hazards.

target_hazard_ratio <- 0.75
target_power <- 0.80
alpha <- 0.05
exposed_proportion <- 0.50
anticipated_event_proportion <- 0.65

events_needed <- ceiling(
  (qnorm(1 - alpha / 2) + qnorm(target_power))^2 /
    (exposed_proportion * (1 - exposed_proportion) * log(target_hazard_ratio)^2)
)

approximate_sample_needed <- ceiling(events_needed / anticipated_event_proportion)

planning_summary <- tibble(
  target_hazard_ratio,
  target_power,
  alpha,
  events_needed,
  anticipated_event_proportion,
  approximate_sample_needed
)

knitr::kable(planning_summary)
target_hazard_ratio target_power alpha events_needed anticipated_event_proportion approximate_sample_needed
0.75 0.8 0.05 380 0.65 585

This is a simplified planning calculation. Adjust the design for unequal groups, covariate correlation, attrition, clustering, nonproportional hazards, competing events, and the actual recruitment and follow-up plan.

Dynamic results template

Kaplan-Meier estimates and a Cox proportional hazards model were used to examine time to the doctoral milestone. Among 450 participants, 291 experienced the event and 159 were right censored. The unadjusted log-rank comparison produced chi-square(1) = 11.28, p < .001.

In the adjusted Cox model, enhanced support relative to standard support was associated with a hazard ratio of 1.35, 95% CI [1.05, 1.72], p = .017. The model included enrollment status, age, baseline progress, and advisor support. The concordance statistic was 0.6.

Proportional-hazards assumptions were evaluated using scaled Schoenfeld residuals and the global test. Continuous-predictor functional form and influential observations were also reviewed. Interpret hazard ratios together with adjusted survival curves, absolute event probabilities, diagnostic results, and the study design.

When this template is not enough

  • Competing risks: another event prevents or changes the probability of the event of interest.
  • Time-varying predictors: exposure values change during follow-up.
  • Delayed entry: participants become at risk after the common time origin.
  • Recurrent events: the same participant can experience the event more than once.
  • Clustered data: participants are nested within programs, clinics, schools, or other units.
  • Interval censoring: the event is known only to occur within an interval.

Each situation requires data structures and estimators beyond the basic one-row-per-participant Cox model used here.

Reporting checklist

  • Define the time origin, event, unit of time, and censoring mechanisms.
  • Report the analysis sample, number of events, censoring proportion, and follow-up.
  • Show Kaplan-Meier curves with confidence intervals and numbers at risk.
  • Report hazard ratios with confidence intervals and clarify reference groups and scaling.
  • State the ties method and missing-data approach.
  • Describe proportional-hazards, functional-form, and influence diagnostics.
  • Report absolute survival or cumulative-incidence estimates when they aid interpretation.
  • Document departures from the prespecified model and any sensitivity analyses.

References

  • Cox, D. R. (1972). Regression models and life-tables. Journal of the Royal Statistical Society: Series B, 34(2), 187-220.
  • Kaplan, E. L., & Meier, P. (1958). Nonparametric estimation from incomplete observations. Journal of the American Statistical Association, 53(282), 457-481.
  • Kleinbaum, D. G., & Klein, M. (2012). Survival Analysis: A Self-Learning Text (3rd ed.). Springer.
  • Therneau, T. M., & Grambsch, P. M. (2000). Modeling Survival Data: Extending the Cox Model. Springer.

Exported files

readr::write_csv(raw_data, "Output/synthetic_survival_data.csv", na = "")
readr::write_csv(event_summary, "Output/event_summary.csv")
readr::write_csv(missingness, "Output/missingness_review.csv")

tibble(
  file = c(
    "synthetic_survival_data.csv",
    "event_summary.csv",
    "missingness_review.csv",
    "cox_model_results.csv",
    "model_performance.csv",
    "proportional_hazards_test.csv",
    "influence_review.csv"
  )
) %>%
  knitr::kable()
file
synthetic_survival_data.csv
event_summary.csv
missingness_review.csv
cox_model_results.csv
model_performance.csv
proportional_hazards_test.csv
influence_review.csv

Session information

sessionInfo()
R version 4.6.0 (2026-04-24 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 10 x64 (build 19045)

Matrix products: default
  LAPACK version 3.12.1

locale:
[1] C
system code page: 65001

time zone: America/New_York
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] glue_1.8.1      broom_1.0.12    survival_3.8-6  lubridate_1.9.5
 [5] forcats_1.0.1   stringr_1.6.0   dplyr_1.2.1     purrr_1.2.2    
 [9] readr_2.2.0     tidyr_1.3.2     tibble_3.3.1    ggplot2_4.0.3  
[13] tidyverse_2.0.0

loaded via a namespace (and not attached):
 [1] sass_0.4.10        generics_0.1.4     stringi_1.8.7      lattice_0.22-9    
 [5] hms_1.1.4          digest_0.6.39      magrittr_2.0.5     evaluate_1.0.5    
 [9] grid_4.6.0         timechange_0.4.0   RColorBrewer_1.1-3 fastmap_1.2.0     
[13] jsonlite_2.0.0     Matrix_1.7-5       backports_1.5.1    mgcv_1.9-4        
[17] scales_1.4.0       jquerylib_0.1.4    cli_3.6.6          rlang_1.3.0       
[21] crayon_1.5.3       bit64_4.8.0        splines_4.6.0      withr_3.0.3       
[25] cachem_1.1.0       yaml_2.3.12        otel_0.2.0         parallel_4.6.0    
[29] tools_4.6.0        tzdb_0.5.0         vctrs_0.7.3        R6_2.6.1          
[33] lifecycle_1.0.5    bit_4.6.0          vroom_1.7.1        pkgconfig_2.0.3   
[37] pillar_1.11.1      bslib_0.10.0       gtable_0.3.6       xfun_0.57         
[41] tidyselect_1.2.1   knitr_1.51         farver_2.1.2       nlme_3.1-169      
[45] htmltools_0.5.9    labeling_0.4.3     rmarkdown_2.31     compiler_4.6.0    
[49] S7_0.2.2
Dissertation Stats Helper · Reproducible statistical guidance for doctoral research · Survival models require explicit definitions of time, event, censoring, and the population at risk.