Checking the proportional hazards assumption in R

To check the proportional hazards assumption in R, fit the Cox model with coxph(), pass it to cox.zph(), and look at the plot it draws. This post does that on the veteran data from the survival package, where the assumption fails for one predictor. It then applies the two standard fixes, stratification and a time-varying effect, and shows what each one changes in the answer.

The data come from a randomised trial of two treatment regimens for lung cancer and are used here only to illustrate the statistics. Nothing below is a statement about any person’s outlook or about treatment. If you are new to the model, the guide to Kaplan-Meier and Cox regression explains censoring and hazards first, and Cox regression in R walks through reading the main output. This post picks up at the check that those two only touch briefly.

What the assumption says

A Cox model describes each person’s hazard, the rate of the event at a moment given that they have not yet had it, as a shared baseline hazard multiplied by a factor that depends on their predictors. The baseline can have any shape. The multiplier does not change with time. So if one group has twice the hazard of another at day 10, the model says it has twice the hazard at day 400 as well. That is the proportional hazards assumption, and it is what lets one hazard ratio per predictor stand for the whole follow-up (Cox 1972).

The assumption is about the ratio, not the hazards themselves. Both hazards can fall steeply over time and the model is still fine, provided they fall in step. It fails when the ratio drifts: a predictor that matters a great deal early and fades, or one whose effect grows. How to read and report a ratio once you have one is covered in how to interpret a hazard ratio.

When the assumption fails, the single hazard ratio the model reports is still a number, but it is an average of a ratio that changed. It describes some weighted mix of early and late follow-up, and the weights depend on how many events happened when. That is why the check is worth running before you report anything.

The data and the question

The veteran data frame has 137 rows, one per person, so the observations are independent. time is the survival time in days and status is 1 for died and 0 for censored. Censoring is right-censoring, and only 9 of the 137 people were censored, so the 128 deaths are almost the whole sample. trt is the treatment, 1 for the standard regimen and 2 for the test regimen, with 69 and 68 people. karno is the Karnofsky performance score, from 10 to 99 in these data, where a higher score means better functional status (100 is good, according to the dataset’s help page). age is in years.

The question for the worked example is: adjusting for the Karnofsky score and age, how is the test treatment associated with the hazard of death, and do the hazards stay proportional over follow-up? Treatment is entered as the number 1 or 2, so its hazard ratio is the test regimen against the standard one. The Karnofsky score is entered as a number too, so its hazard ratio is per one-point step. The model has 128 events and 3 predictors, about 43 events per predictor, which clears the usual rule-of-thumb floor of 10. That rule is a convention, not a guarantee of a stable model.

The script

Here is the hand-written version. It is plain survival code and nothing else needs installing:

library(survival)

fit <- coxph(Surv(time, status) ~ trt + karno + age, data = veteran)
summary(fit)

zph <- cox.zph(fit)   # Schoenfeld residual test, per term and global
print(zph)
plot(zph[2])          # the karno panel: coefficient against time
abline(h = coef(fit)["karno"], lty = 3)   # the fitted, constant coefficient

According to the survival manual, cox.zph() takes a transform argument whose default is "km", the Kaplan-Meier transformation of time, and it reports one test per model term plus a global test. The numbers below come from two sources, and the post says which each time. Numbers labelled “engine run” come from the Celeus run of this model, recorded in its results file, whose R script you can download below the article. Numbers labelled “extra R call” come from running the hand-written code (and variants) with survival 3.8-6 under R 4.6.1, because the engine’s record does not hold every figure.

Step 1: test each term and the whole model

The test of Grambsch and Therneau (1994) works on Schoenfeld residuals. At each death time, a residual compares the predictor values of the person who died with a weighted average over everyone still at risk. If the hazard ratio is constant, those residuals show no trend over time. cox.zph() scales them and tests, for each term, whether a straight-line trend against the transformed time is zero. Engine run:

The Karnofsky score is the problem. Its test gives a chi-square of 12.00 on 1 degree of freedom (p = 0.00053), and the global test, which asks the same question of the whole model at once, gives 19.10 on 3 degrees of freedom (p = 0.00026). Treatment (p = 0.594) and age (p = 0.147) show no evidence of a drift. The engine’s assumption check records the same finding: “evidence the proportional-hazards assumption is violated”. The three per-term p-values are not corrected for multiplicity, but the global p-value is a single joint test, and it agrees.

The headline of this engine record is the treatment row, so its summary line reports a hazard ratio of 1.209 (95% CI 0.840 to 1.739, p = 0.307). Read that with care. It is a fair summary only if the other terms in the model are well specified, and one of them is not.

Step 2: look at the plot

A p-value says whether a trend is detectable. The plot says what the trend looks like, which decides how you fix it. plot(zph[2]) draws the estimated coefficient for the Karnofsky score as a smooth curve against the transformed time (the axis is labelled in days), with a band of two standard errors around it and the residuals as points. The survival manual says that if proportional hazards holds, the true curve is a horizontal line. The plot does not draw the model’s own coefficient, so the script adds it as a dotted line with abline().

In an extra R call the curve for the Karnofsky score started near -0.06 early in follow-up, rose to about 0 around day 130, and drifted back down at the right-hand end, where the band widens because few people remain. The dotted line at the fitted coefficient, -0.034, is clearly not a fair description of the whole curve.

In words: a lower Karnofsky score was strongly associated with a higher hazard of death in the first few months, and that association faded after that. A more precise statement about the later months is hard to make from 128 events, so be wary of reading the far right of the curve closely. Two checks are worth doing before you fix anything: look for a handful of very late deaths driving the end of the curve, and confirm the predictor is coded the way you think it is.

The result depends on the time transform

The test needs a function of time to trend against. The default is "km", the Kaplan-Meier transformation. You can ask for others, and each one spreads the weight over early and late follow-up differently. In an extra R call, the Karnofsky score’s p-value was 0.031 with transform = "identity" (raw time), 0.0020 with "log", and 0.00038 with "rank". The global p-value was 0.017, 0.00073 and 0.00019 for the same three transforms.

Two lessons follow. The conclusion for the Karnofsky score does not depend on the transform, since all four versions are below 0.05. But a term’s p-value can move a long way with the choice: under the log transform, age had p = 0.062, against 0.147 under the default, so a term nearer the line could cross 0.05 on the transform alone. So choose the transform before you look, state it, and treat a p-value near 0.05 as a prompt to look at the plot rather than as a verdict.

Transform sensitivity is also one reason a formal test is not the whole check. Stensrud and Hernán (2020) argue that testing for proportional hazards and then choosing an analysis by the result is a weak basis for decisions, and that what matters is whether a single hazard ratio answers your question.

Fix 1: stratify on the variable that fails

Stratifying gives each level of a variable its own baseline hazard. The hazard functions of different strata can then differ in any way at all, so the proportionality assumption is no longer needed for that variable. The price is that the model reports no hazard ratio for it. That suits a variable you want to adjust for and not describe. The Karnofsky score is a number, so stratifying needs groups. I chose three: under 60, 60 to 79, and 80 or more, which hold 52, 52 and 33 people. Those cut points are my choice and are arbitrary, so say how you chose yours.

veteran$kgrp <- cut(veteran$karno, c(0, 59, 79, 100),
                    labels = c("<60", "60-79", "80+"))
fit_strat <- coxph(Surv(time, status) ~ trt + age + strata(kgrp),
                   data = veteran)

In an extra R call, the stratified model gave a hazard ratio for treatment of 0.980 (95% CI 0.675 to 1.423) and for age of 0.998 (0.980 to 1.017). Both intervals comfortably include 1. The treatment estimate moved from 1.209 to 0.980 once the Karnofsky score had its own baseline per group, which is the point: leaving a badly fitting term in the model can distort the estimates for the others. Coarse groups also adjust for the score differently from a straight line, so not all of the shift need come from the fix. In this case the data are compatible with no strong association of treatment with the hazard, and with a modest change in either direction.

Stratifying is not a clean bill of health. In the same extra call, cox.zph() on the stratified model gave a global p-value of 0.028, with p = 0.070 for treatment and 0.102 for age. Coarse groups may leave some of the Karnofsky effect inside the strata, and with two terms tested a single global value just under 0.05 is a prompt to look at the plots again, not a reason to stop.

Fix 2: let the effect change with time

If the Karnofsky score is the exposure you care about, you cannot stratify it away. The alternative is to model the change. The cleanest way to see it is to split follow-up into periods and fit a separate coefficient in each. The plot suggested the break is a few months in, and the median observed time in these data (deaths and censored times together) is 80 days (an extra R call), so I split at 90 days.

split <- survSplit(Surv(time, status) ~ ., data = veteran,
                   cut = 90, episode = "period")
split$period <- factor(split$period, labels = c("early", "late"))

fit_tv <- coxph(Surv(tstart, time, status) ~ trt + age +
                  karno:strata(period), data = split)

survSplit() cuts each person’s follow-up at day 90, so 137 people become 198 rows, with a start time and a stop time. Writing karno:strata(period) gives the Karnofsky score a separate coefficient in each period. Inside an interaction, strata() only labels the periods: it adds no separate baseline hazard and no main effect for period, and none is needed, because the period is a function of time alone. In an extra R call, the hazard ratio per Karnofsky point was 0.951 (95% CI 0.940 to 0.963) in the first 90 days, with 73 deaths, and 1.001 (0.982 to 1.020) afterwards, with 55 deaths.

To read that: in the first 90 days, a person scoring 10 points lower had an estimated hazard 1.65 times as high (1 / 0.951 to the power of 10, a derived figure). After day 90, the data show no association of the score with the hazard, and the interval is compatible with a small effect in either direction. The difference between the two coefficients is itself estimable. Fitting the same model with karno * period gave an interaction hazard ratio of 1.052 (95% CI 1.028 to 1.076, p < 0.001), the factor by which the per-point ratio changes after day 90. In that parameterisation, R prints NA for the period main effect, because a period term is absorbed by the baseline hazard.

How much does the cut point matter? A split at 100 days gave 0.954 early and 0.999 late, and a split at 60 days gave 0.952 and 0.991. The picture is stable: a clear early effect, little later. But the cut is a choice made after seeing the data, so treat the early-period estimate as descriptive and not as a pre-planned finding. There is a smoother alternative, which interacts the term with a function of time. In an extra R call, tt(karno) with x * log(t) gave an interaction coefficient with a hazard ratio of 1.013 (95% CI 1.005 to 1.022, p = 0.0023), which agrees in direction: the effect of a lower score shrinks as follow-up lengthens (Bellera et al. 2010 set out these models and their pitfalls).

Which fix to use

The choice follows the question, not the p-value.

• A nuisance variable that fails: stratify on it. You lose its hazard ratio and keep honest estimates for the rest.
• A variable you want to describe that fails: model the change, with periods or a time interaction, and report the effect by period.
• A treatment effect that fails: say so plainly. A single hazard ratio is then a weighted average, and a difference in survival at a stated time, or the restricted mean survival time, may answer the question better.
• A significant test that you think is trivial: look at the plot and the size of the drift. A large trial can flag a drift too small to matter, and the plot shows that quickly.

Whichever you choose, report what the check showed, the transform you used, and the fix, so a reader can see that the assumption was examined and not assumed.

How the engine run matches your hand-written code

The Celeus run of this model is a recorded, seeded analysis whose reproduction script you can download. It is the same Cox model, and the correspondence is close enough to follow step by step.

The engine run covers the check and the main model. The stratified and period-split fits are extra R calls, not part of the recorded run, so reproduce those from the code in this post. The engine adds one note the hand-written script does not: values beyond 1.5 times the interquartile range in time (12 values), which are descriptive flags and remove nothing. It also raises a warning when the Schoenfeld test is below 0.05, and the warning’s own suggestion is the two fixes above. The guide to reproducible statistical analysis explains why keeping the seed, the versions and the script together matters beyond one model.

Pitfalls to avoid

A few mistakes come up often when checking this assumption.

• Reading a large p-value as proof. No evidence of a trend is not evidence of none. With few events the test has little power, so look at the plot as well.
• Testing and then switching models by the result alone. Decide what you will do before you look, and base the decision on the size and shape of the drift.
• Ignoring the transform. The default is "km". Name the one you used, and do not shop among transforms for the answer you want.
• Stratifying the exposure. You cannot report a hazard ratio for a stratified variable.
• Reading the far right of the plot. Few people remain there, the band is wide, and a few late deaths can bend the curve.
• Treating a coded category as a number. A variable with three or more categories coded 1, 2, 3 must be a factor, or the model assumes equal steps, which is its own form of misspecification.
• Choosing the cut point after the fact and calling it a finding. Say it was chosen after seeing the plot.

Where to go next

For the ideas behind the model, read the guide to Kaplan-Meier and Cox regression, and for a longer walkthrough of the whole workflow see the survival analysis tutorial. To see the same model’s output read line by line, go to Cox regression in R. For the method’s exact options in the Celeus engine, see the reference page for Cox regression. The example script below reproduces every engine number in this post.

Everything above is research statistics, not clinical decision support: a hazard ratio in a public teaching dataset describes an association in that sample and is not an estimate for any individual.

# Results

- Test: Cox proportional-hazards regression
- n: 137
- statistic: 1.022
- p-value: 0.307
- effect size: hazard ratio = 1.209 (95% CI 0.8402 to 1.739)
- 95% CI (trt): [0.8402, 1.739]
- seed: 20261028

## Assumption checks
- proportional hazards (cox.zph): passed=FALSE (evidence the proportional-hazards assumption is violated)
- events per variable: passed=TRUE (at or above the rule-of-thumb floor (a convention, not a guarantee of a stable model))

## Term tests
| Term | Test | Statistic | df | Residual df | p-value |
|---|---|---|---|---|---|
| trt | LR | 1.042 | 1 | - | 0.307 |
| karno | LR | 42.52 | 1 | - | < 0.001 |
| age | LR | 0.1754 | 1 | - | 0.675 |

## Power
- Power is not calculated for this test.

# Celeus - reproduce.R  (run from this package folder; bundled data in data/)
# Test: Cox proportional-hazards regression  |  engine 0.32.1  |  catalog c39-d3cd
# Research context: not provided for this run.
# Reproduces the recorded result deterministically.

set.seed(20261028)
.data_path <- "data/proportional-hazards-assumption.csv"
if (!file.exists(.data_path)) stop("Run reproduce.R from the package folder (the bundled data is in data/).")
data <- readr::read_csv(.data_path, show_col_types = FALSE)  # needs the readr package
clean_cells <- function (df) 
{
    ctrl_collapse <- "[\r\n\t\v\f]+"
    zero_width_chars <- intToUtf8(c(8203L, 65279L), multiple = TRUE)
    uni_space_chars <- intToUtf8(c(160L, 8239L, 8199L), multiple = TRUE)
    clean_text <- function(x) {
        out <- x
        for (ch in zero_width_chars) out <- gsub(ch, "", out, 
            fixed = TRUE)
        out <- gsub(ctrl_collapse, " ", out)
        for (ch in uni_space_chars) out <- gsub(ch, " ", out, 
            fixed = TRUE)
        trimws(gsub(" +", " ", out))
    }
    changed_cols <- character(0)
    cells_cleaned <- 0L
    retyped <- character(0)
    for (j in seq_along(df)) {
        col <- df[[j]]
        if (!is.character(col)) 
            next
        cleaned <- clean_text(col)
        n_changed <- sum(!is.na(col) & !is.na(cleaned) & col != 
            cleaned, na.rm = TRUE)
        if (n_changed > 0L) {
            changed_cols <- c(changed_cols, names(df)[j])
            cells_cleaned <- cells_cleaned + n_changed
        }
        df[[j]] <- cleaned
        non_blank <- cleaned[!is.na(cleaned) & nzchar(cleaned)]
        if (length(non_blank) && all(grepl("^[+-]?(0|[1-9][0-9]*)([.][0-9]+)?$", 
            non_blank))) {
            parsed <- suppressWarnings(as.numeric(cleaned))
            was_blank <- is.na(cleaned) | !nzchar(cleaned)
            lost <- (is.na(parsed) | is.infinite(parsed)) & !was_blank
            if (!any(lost)) {
                finite <- parsed[!is.na(parsed)]
                all_whole <- all(finite == floor(finite))
                in_range <- all(abs(finite) <= .Machine$integer.max)
                df[[j]] <- if (all_whole && in_range) 
                  as.integer(parsed)
                else parsed
                retyped <- c(retyped, names(df)[j])
            }
        }
    }
    names(df) <- clean_text(names(df))
    if (length(changed_cols)) {
        attr(df, "ingest_clean_columns") <- changed_cols
        attr(df, "ingest_clean_cells_cleaned") <- cells_cleaned
    }
    if (length(retyped)) 
        attr(df, "ingest_clean_retyped") <- retyped
    df
}
data <- clean_cells(data)
# Complete-case analysis: drop rows missing any analysis value.
.acols <- c("trt", "karno", "age", "time", "status")
data <- data[stats::complete.cases(data[, .acols, drop = FALSE]), , drop = FALSE]
result <- { .fit <- survival::coxph(survival::Surv(`time`, `status`) ~ `trt` + `karno` + `age`, data = data); list(summary = summary(.fit), ph = survival::cox.zph(.fit), term_tests = stats::drop1(.fit, test = "Chisq")) }
print(result)