small <- data.frame(
id = 1:5,
time = c(5, 3, 8, 2, 6),
status = c(1, 0, 1, 0, 1)
)
small id time status
1 1 5 1
2 2 3 0
3 3 8 1
4 4 2 0
5 5 6 1Survival Analysis in Base R
01 / A practical introduction to standard survival analysis
Survival analysis concerns the time until an event, such as relapse or death. Follow-up often ends before the event is observed for every participant: some remain event-free at the study end, while others leave earlier. The analysis must account for these incomplete observations while using the information available for each participant.
A standard survival analysis involves defining the outcome, estimating survival curves, comparing groups, fitting regression models, and checking their assumptions. Base R and the survival package provide the tools for these tasks, without requiring tidyverse knowledge.
The standard tools produce usable preliminary results, including console summaries, basic plots, and fitted model objects. Understanding these results is essential for interpreting the estimates and assessing the analysis. Their default presentation often requires additional formatting before inclusion in a report.
Tidy survival tools build on the same statistical methods while simplifying data preparation, extraction of model results, and production of clear figures and tables. A consistent workflow reduces repeated manual formatting and makes reports easier to update when the data or analysis change.
Run the examples in order from the RStudio project folder so that the data paths and objects are available to subsequent commands.
1.1 Survival outcomes and basic concepts
1.1.1 A time origin, an event, and an end to observation
Before choosing a method, specify when the clock starts and what stops it. For overall survival, the origin might be the date of diagnosis and the event is death. For relapse-free survival, the event might be the first of relapse or death. These are different outcomes, even when constructed from the same patient records.
A participant whose event is observed contributes an event time. A participant whose event has not been observed contributes a shorter statement: the event time exceeds the time of last observation. We call the latter observation right censored. Someone followed for 32 months without an event has not demonstrated that their event time is 32 months. They have demonstrated that it is longer than 32 months.
Censoring therefore carries information. We should neither discard a censored participant nor count the censoring time as an event. We also should not treat that person as permanently event-free. Survival methods retain their contribution for the period during which their event status is known.
Let \(T\) be the event time and \(C\) the censoring time. We observe
\[ X = \min(T,C), \qquad \delta = I(T \leq C). \]
The indicator \(I(\cdot)\) equals one when its argument is true and zero otherwise. The pair \((X,\delta)\) records both the length of observation and how it ended. In R, we will usually call these columns time and event, or time and status.
1.1.2 A small example
Consider five people followed from the same kind of starting point. We can create their records with data.frame(), which assembles columns into a rectangular dataset.
The first person had an event at time 5. The second was observed until time 3 without an event, so their event time is greater than 3. The fourth was censored at time 2. The event-free information contributed by these two people is shorter, but still useful.
The function Surv() keeps the two parts of each observation together. Load the package once in the current R session, then pass it the time and event-indicator columns. The $ operator selects a named column from a data frame.
library(survival)
Surv(small$time, small$status)[1] 5 3+ 8 2+ 6 The printed result is 5 3+ 8 2+ 6. The plus sign marks censoring; it does not mean that we should add anything to the time. This compact display represents the same five pairs of numbers in small. Here, 1 means an observed event and 0 means censoring. Other datasets use other conventions, so checking the dictionary comes before constructing a survival object.
A follow-up plot makes the distinction visible. Each horizontal line below represents one subject’s observed follow-up; its endpoint symbol distinguishes an event from censoring. The absence of an event at a line’s endpoint does not tell us what happened afterward.
1.1.3 The assumption behind using censored observations
The usual methods require censoring to be uninformative about future event risk after accounting for the variables used in the analysis. Informally, among comparable people still being followed at a given time, those censored then should not systematically have different subsequent event prospects from those who remain under observation.
A planned study end can make this assumption plausible. Withdrawal caused by worsening health may make it doubtful. The observed data alone generally cannot establish the assumption: the outcomes needed to check it are unobserved after censoring. Knowing how follow-up was collected is part of knowing what the analysis means.
1.1.4 The quantities we want to describe
The survival function answers a probability question:
\[ S(t) = \Pr(T > t). \]
For example, \(S(36)=0.70\) means a 70% probability of remaining free of the defined event beyond 36 time units. If relapse or death is the event, the curve describes remaining alive and relapse-free; it is not an overall-survival curve.
The hazard asks a different question. Among people who have reached time \(t\) without an event, how quickly are events occurring at that time? For a continuous event time,
\[ \lambda(t) = \lim_{\Delta t \to 0} \frac{\Pr(t \leq T < t+\Delta t \mid T \geq t)}{\Delta t}. \]
The fraction restricts attention to people still event-free and considers a short interval ahead. Dividing by the interval length produces a rate. A hazard is not a probability and need not lie between zero and one. When time is measured in months, it is a rate per month among those still at risk.
The cumulative hazard, \(\Lambda(t)=\int_0^t \lambda(u)\,du\), accumulates that rate. For a continuous event-time distribution, \(S(t)=\exp\{-\Lambda(t)\}\). It explains why survival probabilities and hazards both appear in the software: they describe the same distribution from different perspectives.
Keep the interpretation attached to the quantity. A survival probability concerns remaining event-free up to a specified time. A hazard ratio compares event rates among those still event-free. It is not a ratio of survival probabilities.
1.2 Meet the working data
The supplied German breast cancer (GBC) analysis files contain records for 686 patients with node-positive breast cancer. The data dictionary describes their provenance and coding; the original study recruited a larger group.
Place the files in a folder called data inside your RStudio project. With that project open, these paths are relative to its folder. header = TRUE tells read.table() to read column names from the first line.
mortality <- read.table("data/gbc_mort.txt", header = TRUE)
gbc <- read.table("data/gbc.txt", header = TRUE)
head(mortality[, c("id", "time", "status", "hormone", "age")]) id time status hormone age
1 1 74.819672 0 1 38
2 2 65.770492 0 1 52
3 3 47.737705 1 1 47
4 4 4.852459 0 1 40
5 5 61.081967 0 2 64
6 6 63.377049 0 2 49The square brackets select rows and columns. A blank before the comma means all rows; the vector after it names the columns to retain. head() shows just the first six rows, keeping the output manageable.
| File | What a row records | Meaning of status |
|---|---|---|
gbc_mort.txt |
One patient’s death or last mortality follow-up | 0 = censored; 1 = death |
gbc.txt |
A relapse, death, or censoring record; patients can have multiple records | 0 = censored; 1 = relapse; 2 = death |
The baseline variables also need interpretation. hormone is coded 1 for no hormone therapy and 2 for hormone therapy; meno is coded 1 for premenopausal and 2 for postmenopausal. age is in years, size in millimeters, and grade takes values 1, 2, or 3. nodes counts positive lymph nodes. prog and estrg record progesterone and estrogen receptor levels.
A few checks help us avoid counting records as people.
c(
mortality_rows = nrow(mortality),
event_history_rows = nrow(gbc),
patients = length(unique(gbc$id)),
observed_deaths = sum(mortality$status == 1)
) mortality_rows event_history_rows patients observed_deaths
686 985 686 171 There are 686 mortality records but 985 event-history records, representing the same 686 patients. The mortality file contains 171 observed deaths. unique(gbc$id) extracts distinct patient identifiers, and length() counts them. In contrast, nrow(gbc) counts every record, including later records from the same patient.
Inspecting one history makes the distinction concrete.
gbc[gbc$id == 3, c("id", "time", "status")] id time status
5 3 41.93443 1
6 3 47.73770 2Patient 3 has a relapse at about 41.9 months and a death at about 47.7 months. For overall survival, the latter time matters. For relapse-free survival, the event already happened at the earlier time. Counting both rows as independent patients would change the analysis for a purely clerical reason.
1.2.1 Build one record per patient
The outcome for this analysis is time to first relapse or death. Each patient contributes one observation, ending at that first event or at the last known event-free time. We create a separate object, rfs, leaving the source data available for other questions.
ordered_records <- gbc[order(gbc$id, gbc$time, gbc$status == 0), ]
rfs <- ordered_records[!duplicated(ordered_records$id), ]
rfs$event <- as.integer(rfs$status > 0)
head(rfs[, c("id", "time", "status", "event")]) id time status event
1 1 43.836066 1 1
3 2 46.557377 1 1
5 3 41.934426 1 1
7 4 4.852459 0 0
8 5 61.081967 0 0
9 6 63.377049 0 0Read this in three steps:
order(gbc$id, gbc$time, gbc$status == 0)sorts first by patient, then by time. Its third argument places an event before a censoring record when their times tie, becauseFALSEsorts beforeTRUE.!duplicated(ordered_records$id)retains each patient’s first record.duplicated()marks identifiers already encountered;!reverses that result. The blank after the comma keeps all columns.as.integer(rfs$status > 0)creates the event indicator. The comparison asks whether either event occurred, andas.integer()convertsTRUEandFALSEinto 1 and 0.
The tie rule suffices for this composite because either event ends relapse-free survival. It does not establish which event type occurred first when two types share a recorded time. A competing-risk analysis needs its own outcome rules.
c(patients = nrow(rfs), events = sum(rfs$event))patients events
686 299 anyDuplicated(rfs$id)[1] 0The result is 686 patients and 299 first events. Zero from anyDuplicated() confirms that identifiers are not repeated. The event count exceeds the death count because our composite includes relapse. These are observed counts, not estimated probabilities at a common follow-up time.
We also label the categorical variables. A factor tells R to treat a column as categories and records their order.
rfs$hormone <- factor(rfs$hormone, levels = c(1, 2),
labels = c("No", "Yes"))
rfs$meno <- factor(rfs$meno, levels = c(1, 2),
labels = c("Pre", "Post"))
rfs$grade <- factor(rfs$grade, levels = c(1, 2, 3),
labels = c("I", "II", "III"))
table(rfs$hormone, rfs$event)
0 1
No 235 205
Yes 152 94The reference levels are now explicit: no hormone therapy, premenopausal status, and grade I. Treating grade as a factor allows separate comparisons of grades II and III with grade I. Leaving it numerical in a regression would impose an equal change in log hazard for each one-grade increase.
Exercise: change the outcome
Which file would you use to count observed deaths? Why would sum(rfs$event) answer a different question?
The mortality file has one record per patient and codes death as status = 1. The event column in rfs also includes relapse. A matching command is:
sum(mortality$status == 1)[1] 1711.3 Estimate survival and compare groups
1.3.1 Fit a Kaplan–Meier curve
The Kaplan–Meier estimator follows the group through successive event times. At each event time it compares the number of events with the number still under observation and event-free just before that time—the risk set. It combines these conditional survival proportions to estimate survival from the origin.
If five people are at risk and one event occurs, the estimated fraction surviving that step is \(4/5\). If someone is subsequently censored, the curve does not jump at censoring, but later risk sets become smaller. This is how incomplete follow-up contributes without being counted as a failure time.
We ask survfit() for separate curves by hormone-therapy group:
km_fit <- survfit(Surv(time, event) ~ hormone, data = rfs)
km_fitCall: survfit(formula = Surv(time, event) ~ hormone, data = rfs)
n events median 0.95LCL 0.95UCL
hormone=No 440 205 50.1 42.5 59.5
hormone=Yes 246 94 66.2 62.9 NAThe expression around ~ is an R formula. Read the call in three parts:
Surv(time, event)supplies follow-up time and the event indicator on the left of the formula.~ hormonerequests a separate curve for each hormone-therapy group. Use~ 1for one curve combining all patients.data = rfstells R where to find these variables, so we do not need to repeatrfs$before each name.
The assignment stores a fitted object in km_fit. Printing it gives group sizes, observed events, and estimated median survival. Here median survival means the time when the estimated survival curve reaches 0.5. It is not the median of the observed time column, which mixes event and censoring times. An NA for a median or confidence limit can mean that the relevant curve never reaches 0.5 during follow-up.
1.3.2 Ask about specific times
A question about survival at one, three, and five years corresponds to 12, 36, and 60 months in our data.
summary(km_fit, times = c(12, 36, 60))Call: survfit(formula = Surv(time, event) ~ hormone, data = rfs)
hormone=No
time n.risk n.event survival std.err lower 95% CI upper 95% CI
12 379 44 0.897 0.0148 0.868 0.926
36 195 114 0.606 0.0247 0.559 0.656
60 63 41 0.437 0.0297 0.382 0.499
hormone=Yes
time n.risk n.event survival std.err lower 95% CI upper 95% CI
12 223 12 0.950 0.0142 0.922 0.978
36 136 54 0.708 0.0305 0.650 0.770
60 58 20 0.581 0.0362 0.514 0.657Start with time and survival. At 36 months, estimated relapse-free survival is about 0.606 in the no-hormone group and 0.708 in the hormone group: an unadjusted difference of about 10 percentage points in the probability of remaining alive and relapse-free beyond three years.
n.risk counts patients still contributing to the risk set at the stated time. The standard error and confidence limits describe uncertainty in the survival estimate. With successive requested times, n.event reports event counts over the intervals between them, including the first interval from the origin. It does not mean events occurred exactly at 12, 36, or 60 months.
Risk sets shrink with follow-up. A late estimate based on few patients deserves more caution than an early estimate based on most of the sample. Confidence intervals communicate sampling uncertainty, but do not repair informative censoring or a poorly defined outcome.
1.3.3 Make a first plot
Base R already gives us a useful survival plot. We supply labels and line styles so the reader need not translate numerical treatment codes.
plot(km_fit,
xlab = "Months since study entry",
ylab = "Relapse-free survival probability",
ylim = c(0, 1),
col = c("#172d49", "#b0772b"), lty = c(1, 2),
lwd = 2, conf.int = TRUE, mark.time = TRUE)
legend("bottomleft", legend = c("No hormone therapy", "Hormone therapy"),
col = c("#172d49", "#b0772b"), lty = c(1, 2), lwd = 2,
bty = "n")
The curves descend at events and remain level between them. Censoring marks do not cause downward steps. The hormone-therapy curve generally lies above the other curve, indicating better estimated relapse-free survival in that group. This describes a contrast before adjustment for baseline characteristics; it does not explain why the groups differ.
The plotting arguments control different parts of this preliminary display:
xlabandylabname the axes;ylim = c(0, 1)fixes the probability scale.colandltysupply colors and line types in group order. Here, line types 1 and 2 mean solid and dashed.lwd = 2sets line thickness.conf.int = TRUEadds confidence limits, andmark.time = TRUEmarks censoring times.legend()adds the key separately. Repeating the colors and line types keeps it consistent with the curves;bty = "n"removes its border.
This plot is useful for checking the analysis. Preparing a figure for a report may also require a risk table, more control over labels, and a consistent style across analyses. The tidy plotting tools introduced later make these reporting tasks easier to combine.
The confidence limits are pointwise, concerning a given time rather than simultaneous coverage of the entire curve.
1.3.4 Compare groups with a log-rank test
The curves suggest a difference. A log-rank test provides a formal comparison across follow-up. At each event time it compares events observed in each group with those expected if the groups had the same event-time distribution, given the risk sets. It then combines those comparisons.
logrank_fit <- survdiff(Surv(time, event) ~ hormone, data = rfs)
logrank_fitCall:
survdiff(formula = Surv(time, event) ~ hormone, data = rfs)
N Observed Expected (O-E)^2/E (O-E)^2/V
hormone=No 440 205 180 3.37 8.56
hormone=Yes 246 94 119 5.12 8.56
Chisq= 8.6 on 1 degrees of freedom, p= 0.003 The hormone group has fewer events than expected under the null comparison: 94 observed versus about 119 expected. The no-hormone group has the complementary excess. The chi-squared statistic is about 8.6 on one degree of freedom, with a p-value near 0.003. This is evidence against equal survival distributions under the test’s assumptions.
A p-value does not measure the size of the difference; the survival estimates at meaningful times do that more directly. Nor is this an adjustment for age, grade, or other characteristics. The test is often informative when hazards differ consistently over time; opposing differences in different periods can partly cancel. Read the curves along with the test.
If we need the p-value as a separate R object, we can calculate it from the test statistic and number of groups:
pchisq(logrank_fit$chisq,
df = length(logrank_fit$n) - 1,
lower.tail = FALSE)[1] 0.003427282lower.tail = FALSE asks for the probability of a chi-squared value at least as large as the observed statistic. This illustrates a recurring pattern: printed output summarizes a richer R object, and $ retrieves named pieces of that object.
Exercise: compare within menopausal groups
Use a stratified log-rank test to compare hormone groups within menopausal-status strata.
survdiff(Surv(time, event) ~ hormone + strata(meno), data = rfs)Call:
survdiff(formula = Surv(time, event) ~ hormone + strata(meno),
data = rfs)
N Observed Expected (O-E)^2/E (O-E)^2/V
hormone=No 440 205 180 3.52 9.51
hormone=Yes 246 94 119 5.31 9.51
Chisq= 9.5 on 1 degrees of freedom, p= 0.002 strata(meno) forms comparisons separately within menopausal groups before combining them. It does not test a menopausal-status coefficient. It also differs from making four treatment-by-menopause groups and testing all four together.
1.4 Cox regression and prediction
1.4.1 Move from group comparisons to covariates
A Cox model relates the event rate to several explanatory variables at once:
\[ \lambda(t \mid Z) = \lambda_0(t)\exp(\beta^T Z). \]
The baseline hazard \(\lambda_0(t)\) describes the time pattern for the reference covariate values. The remaining term multiplies that hazard according to the covariates. The model leaves the baseline time pattern unspecified but assumes that the hazard ratio for a fixed covariate contrast does not change with time. This is the proportional hazards assumption.
We fit an illustrative model using the covariates in the course example. It is a starting point for learning the workflow, not the result of a complete model-building exercise. For example, the dataset also contains node counts, which are not included in this particular formula.
cox_fit <- coxph(
Surv(time, event) ~ hormone + meno + age + grade + size + prog + estrg,
data = rfs, x = TRUE
)
summary(cox_fit)Call:
coxph(formula = Surv(time, event) ~ hormone + meno + age + grade +
size + prog + estrg, data = rfs, x = TRUE)
n= 686, number of events= 299
coef exp(coef) se(coef) z Pr(>|z|)
hormoneYes -0.3552391 0.7010058 0.1292278 -2.749 0.00598 **
menoPost 0.2683165 1.3077610 0.1840516 1.458 0.14489
age -0.0087937 0.9912449 0.0093783 -0.938 0.34842
gradeII 0.7078888 2.0297017 0.2485544 2.848 0.00440 **
gradeIII 0.8171706 2.2640847 0.2689544 3.038 0.00238 **
size 0.0153213 1.0154392 0.0036538 4.193 2.75e-05 ***
prog -0.0022942 0.9977084 0.0005775 -3.972 7.12e-05 ***
estrg 0.0001788 1.0001789 0.0004686 0.382 0.70274
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
hormoneYes 0.7010 1.4265 0.5442 0.9031
menoPost 1.3078 0.7647 0.9117 1.8758
age 0.9912 1.0088 0.9732 1.0096
gradeII 2.0297 0.4927 1.2470 3.3037
gradeIII 2.2641 0.4417 1.3365 3.8355
size 1.0154 0.9848 1.0082 1.0227
prog 0.9977 1.0023 0.9966 0.9988
estrg 1.0002 0.9998 0.9993 1.0011
Concordance= 0.652 (se = 0.016 )
Likelihood ratio test= 72.76 on 8 df, p=1e-12
Wald test = 62.25 on 8 df, p=2e-10
Score (logrank) test = 62.63 on 8 df, p=1e-10The + signs include separate model terms; they do not request interactions. x = TRUE retains the model matrix for diagnostic calculations. The summary reports 686 observations and 299 events. In other datasets these counts may be reduced by missing covariates, so compare them with the intended analysis sample.
1.4.2 Read the coefficient table
The coef column contains estimated log hazard ratios. exp(coef) exponentiates them to obtain hazard ratios, with one as the reference for no association. The coefficient’s standard error, test statistic, and p-value describe uncertainty and a test of a zero log hazard ratio.
The row hormoneYes compares hormone therapy with the reference level No, at the same modeled values of the other covariates. Its hazard ratio is about 0.70, with a 95% confidence interval of about 0.54 to 0.90. Under this fitted model, the hormone group has an estimated event rate about 30% lower among otherwise comparable patients who remain event-free. This is an adjusted association. It is neither a 30-percentage-point increase in survival nor, by this calculation alone, proof of a causal treatment effect.
The grade rows compare grades II and III separately with grade I. This differs from treating grade as a numerical score and imposing the same step from I to II as from II to III.
Units matter for continuous covariates. The size hazard ratio concerns a one-millimeter increase; the age hazard ratio concerns a one-year increase. For a ten-millimeter size comparison, multiply its log hazard coefficient by ten before exponentiating:
exp(10 * coef(cox_fit)["size"]) size
1.165573 The result is about 1.17: a ten-millimeter larger tumor is associated with roughly a 17% higher hazard, conditional on the other terms and on the linear specification for size. Selecting the coefficient by ["size"] makes the intention explicit without relying on its position in a table.
To extract all hazard ratios and confidence intervals without copying console output, use:
round(cbind(
hazard_ratio = exp(coef(cox_fit)),
exp(confint(cox_fit))
), 3) hazard_ratio 2.5 % 97.5 %
hormoneYes 0.701 0.544 0.903
menoPost 1.308 0.912 1.876
age 0.991 0.973 1.010
gradeII 2.030 1.247 3.304
gradeIII 2.264 1.336 3.836
size 1.015 1.008 1.023
prog 0.998 0.997 0.999
estrg 1.000 0.999 1.001The columns place estimates and limits together. round() changes printed precision, not the fitted model.
1.4.3 Obtain predicted survival probabilities
A hazard ratio is a relative comparison. To estimate three-year relapse-free survival for a particular profile, we also need the estimated baseline survival and specified covariates.
Create a data frame with one profile. Its names must match the model variables, and its factor levels must be consistent with the fitted data.
new_patient <- data.frame(
hormone = factor("No", levels = levels(rfs$hormone)),
meno = factor("Pre", levels = levels(rfs$meno)),
age = 45,
grade = factor("II", levels = levels(rfs$grade)),
size = 20,
prog = 100,
estrg = 100
)
new_patient hormone meno age grade size prog estrg
1 No Pre 45 II 20 100 100This is a hypothetical profile, not an “average patient” automatically supplied by R. Making it explicit lets someone else understand and reproduce the prediction. Different profiles produce different curves.
We use survfit() again, now supplying a fitted Cox model and newdata. The same function name serves different tasks because R selects a method appropriate to the supplied object.
predicted_survival <- survfit(cox_fit, newdata = new_patient)
summary(predicted_survival, times = c(12, 36, 60))Call: survfit(formula = cox_fit, newdata = new_patient)
time n.risk n.event survival std.err lower 95% CI upper 95% CI
12 602 56 0.927 0.0118 0.904 0.951
36 331 168 0.673 0.0327 0.612 0.740
60 121 61 0.515 0.0412 0.440 0.603The times argument in summary() chooses points to report from the predicted curve. Read its survival and confidence-limit columns as model-based estimates for people with this profile, not guarantees about an individual’s future. The risk-set and event counts shown here arise from the fitted data; they are not counts of people with exactly these covariates.
At 36 months, the estimated probability of remaining alive and relapse-free is about 0.673, with a 95% confidence interval from 0.612 to 0.740. We could describe this as an estimated three-year relapse-free survival of 67% for the specified profile. The interval reflects uncertainty in that estimated probability; it is not a range of possible survival times for this person.
plot(predicted_survival,
xlab = "Months since study entry",
ylab = "Predicted relapse-free survival probability",
ylim = c(0, 1), col = "#172d49", lwd = 2,
conf.int = TRUE)
The earlier Kaplan–Meier curves summarized observed treatment groups without covariate adjustment. This curve comes from a regression model for a chosen profile. Its interpretation depends on the model being adequate, which brings us to diagnostics.
1.5 Check what the model assumes
1.5.1 Does a covariate’s association change over time?
The Cox model allows the baseline hazard to change with time. Proportional hazards does not require a constant event rate; it requires a fixed covariate contrast to multiply that changing rate by the same factor throughout follow-up.
cox.zph() checks for time-related changes in coefficients using scaled Schoenfeld residuals. Its printed result gives a test for each term and a joint test labeled GLOBAL.
ph_test <- cox.zph(cox_fit)
ph_test chisq df p
hormone 0.248 1 0.6183
meno 5.171 1 0.0230
age 9.319 1 0.0023
grade 11.894 2 0.0026
size 0.847 1 0.3575
prog 4.934 1 0.0263
estrg 5.674 1 0.0172
GLOBAL 23.663 8 0.0026For this fit the global p-value is about 0.003. Several terms, including age and grade, show evidence of departures from a constant coefficient. The hormone term does not show comparable evidence, but that alone does not validate the model or prove a constant hazard ratio.
Plots reveal patterns that a p-value cannot show. Compare the hormone and age terms:
old_par <- par(mfrow = c(1, 2))
plot(ph_test, var = 1, se = TRUE, col = "#172d49")
plot(ph_test, var = 3, se = TRUE, col = "#172d49")
par(old_par)
The term order follows the diagnostic table. The smooth curves estimate how the coefficients vary over time. Under proportional hazards their underlying patterns would be horizontal, though not necessarily at zero. Uncertainty bands help distinguish a trend from a noisy fluctuation, especially near the ends of follow-up. These are coefficient diagnostics, not survival curves.
1.5.2 Is a straight-line covariate effect reasonable?
A separate assumption concerns the covariate scale. Writing + age gives age a linear association with log hazard. This is different from assuming that age’s effect is constant over time: one assumption concerns variation across ages, the other variation across follow-up times.
An exploratory view plots martingale residuals from a model omitting age against age itself. Each residual contrasts the observed event indicator with the cumulative hazard expected from the fitted model over that person’s observed follow-up. These residuals are asymmetric, so we seek a broad pattern rather than normally distributed errors.
fit_without_age <- update(cox_fit, . ~ . - age)
m <- residuals(fit_without_age, type = "martingale")
plot(rfs$age, m,
xlab = "Age at diagnosis (years)",
ylab = "Martingale residual", pch = 16,
col = grDevices::adjustcolor("#172d49", alpha.f = 0.35))
lines(lowess(rfs$age, m), col = "#925415", lwd = 2)
abline(h = 0, lty = 3, col = "grey50")
update() refits the earlier model after removing age. The dots mean “retain the corresponding parts of the previous formula.” lowess() adds a smooth exploratory trend. Curvature can motivate a more flexible age term, but this plot does not directly estimate age-specific hazard ratios or identify a unique transformation. Other misspecified terms and covariate distributions also influence it.
One possible revision uses a natural spline for age and stratifies the baseline hazard by grade:
revised_fit <- coxph(
Surv(time, event) ~ hormone + meno + splines::ns(age, df = 3) +
size + prog + estrg + strata(grade),
data = rfs, x = TRUE
)
cox.zph(revised_fit) chisq df p
hormone 0.000399 1 0.984
meno 3.282677 1 0.070
splines::ns(age, df = 3) 3.785736 3 0.286
size 0.268522 1 0.604
prog 2.613928 1 0.106
estrg 3.630463 1 0.057
GLOBAL 7.844487 8 0.449splines:: identifies the package providing ns(). This term represents age with a smooth curve using three basis functions, allowing more flexibility than one straight-line coefficient. strata(grade) gives each grade its own baseline hazard, so there is no longer a grade hazard-ratio coefficient to report.
This illustrates available tools rather than an automatic repair. A spline addresses age’s functional form, not a time-varying age effect. The revised diagnostics still need examination, and departures may remain. Arbitrarily splitting age at a cutoff would lose information and would not, by itself, solve either diagnostic problem.
1.6 From preliminary results to tidy reporting
A first-event dataset links the outcome definition to the analysis: Surv() represents the observations, survfit() estimates curves, survdiff() compares groups, and coxph() allows covariate adjustment. Predicted probabilities describe outcomes for specified profiles, while diagnostics assess the assumptions underlying those predictions.
Console summaries and base R plots support preliminary interpretation and model checking. Preparing results for a report also requires selecting coefficients, exponentiating estimates and confidence limits, assigning descriptive labels, and coordinating graphical settings. Repeating these steps manually across several analyses increases the work required to maintain consistent output.
Tidy survival tools organize these tasks into a consistent workflow connecting data, fitted models, and presentation. Their usefulness depends on retaining the outcome definitions and diagnostic checks underlying the analysis. The data operations, graphics, and descriptive tables introduced in Chapter 2 provide the main components of this workflow.
The downloadable R script contains the code in reading order. Start it from the project folder with the data directory available, and read its results alongside the explanations here.
Documentation for further study
See the package documentation for Surv(), Cox-model survival predictions, and cox.zph(). In RStudio, commands such as ?survfit open help for your installed version.