Cox regression in R: coxph() with the exact script

To fit a Cox model in R you call coxph() from the survival package on a Surv(time, status) outcome, read the table of hazard ratios it prints, and then check the proportional hazards assumption with cox.zph(). This post does that on one public dataset, line by line, with the whole script in one place. It also shows how a deterministic Celeus run of the same model lines up with the code you would write by hand.

The data are the North Central Cancer Treatment Group (NCCTG) lung cancer data that ship with survival. They are a standard teaching dataset, used here only to illustrate the statistics. Nothing in this post is a statement about any person’s outlook or about treatment. For the ideas behind the model (censoring, the survival curve, what a hazard is) the guide to Kaplan-Meier and Cox regression is the place to start, and the full survival analysis walkthrough puts the Cox model next to Kaplan-Meier curves and an accelerated failure time model. This post stays on the R code and on reading what it prints. How to read and report the hazard ratios themselves is covered in how to interpret a hazard ratio.

The data and the question

The lung data frame has 228 rows, one per person, so each person appears once. People were enrolled at several institutions, a possible clustering the pitfalls section returns to. time is the survival time in days and status is coded 1 for censored and 2 for dead, as the package documentation states. Censoring is right-censoring: a censored person was still alive at their last contact. Three more columns are used:

• age, in years.
• sex, coded 1 for male and 2 for female. It is stored as a number, so its hazard ratio is the contrast of female against male.
• ph.ecog, the ECOG performance score rated by the physician, from 0 (asymptomatic) upwards. It is entered as a number, so its hazard ratio is per one-point step, which assumes every step multiplies the hazard by the same amount.

The question for the worked example is: adjusting for age and sex, how is the ECOG performance score associated with the hazard of death, and do the hazards stay proportional over follow-up? One row has a missing ph.ecog, so the model uses 227 of the 228 rows, and 164 of those 227 people died during follow-up. Events, not people, set the precision of a Cox model, and 164 events for 3 predictors is about 55 events per predictor, comfortably above the usual rule-of-thumb floor of ten.

The script

Here is the hand-written version. It is plain survival code with nothing else to install:

library(survival)

fit <- coxph(Surv(time, status) ~ age + sex + ph.ecog, data = lung)

print(fit)             # short table: coef, exp(coef), se(coef), z, p
summary(fit)           # adds the confidence interval and the three overall tests
exp(confint(fit))      # hazard ratios with 95% Wald intervals

zph <- cox.zph(fit)    # proportional hazards check
print(zph)
plot(zph)              # one panel per term; flat means proportional

Surv(time, status) pairs each time with its status. According to the survival manual, the status can be coded 0/1, TRUE/FALSE or 1/2 with 2 meaning death, which is why the 1/2 coding in lung needs no recoding. The formula lists the predictors, and coxph() drops any row with a missing value in the outcome or in a predictor before fitting (the default na.action comes from options(), per the manual).

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 above (and a few variants) with survival 3.8-6 under R 4.6.1, because the engine’s record does not hold every line that print() shows.

Reading the output line by line

The first thing print(fit) shows is the call, which repeats your formula and data. The next block is one row per predictor, with five columns.

coef is the estimated change in the log hazard for a one-unit increase in the predictor, with the other predictors held fixed. In an extra R call the coefficients were 0.011067 for age, -0.552612 for sex and 0.463728 for ECOG score. The sign is the direction: positive means a higher hazard, negative a lower one. A log scale is hard to read, so the next column undoes it.

exp(coef) is the hazard ratio, the exponential of the coefficient. This is the number you report. The engine run gave 1.011 for age (per year), 0.575 for sex (female against male) and 1.590 for ECOG score (per point).

se(coef) is the standard error of the coefficient on the log scale, 0.009267, 0.167739 and 0.113577 in the same extra call, in the same order.

z is the coefficient divided by its standard error, a Wald statistic. It was 1.194 for age, -3.294 for sex and 4.083 for ECOG score. The last column is the two-sided p-value from comparing z to a standard normal distribution.

summary(fit) repeats that table and adds a second block with exp(coef), exp(-coef) and the lower and upper 95% limits. The limits are Wald intervals: the coefficient plus or minus 1.96 standard errors, exponentiated. They are the intervals to report. The exp(-coef) column is the same ratio turned round (the hazard ratio for a one-unit decrease), which is handy when you would rather say “male against female”.

Below the intervals come the overall lines. Concordance is the share of comparable pairs of people in which the one with the higher predicted risk had the earlier event. It was 0.637 with a standard error of 0.025 in an extra R call. A value of 0.5 would be a coin toss, so this model ranks pairs only modestly better than chance, and concordance describes ranking, not the size of any hazard ratio. Then come three tests of the whole model against a model with no predictors: the likelihood ratio test (30.50 on 3 degrees of freedom in the engine run), the Wald test (29.93) and the score or log-rank test (30.50). They test the same null hypothesis in three ways and agree closely when the sample is large, as it is here. Finally, n= 227, number of events= 164 is the line to check first. It tells you how many rows the model actually used and how many events they held, and a gap between that n and your dataset’s row count is where missing-data surprises show up.

The hazard ratios in a table

Putting the three rows in one place, from the engine run:

Term	Hazard ratio	95% CI	p
Age (per year)	1.011	0.993 to 1.030	0.232
Sex (female vs male)	0.575	0.414 to 0.799	< 0.001
ECOG score (per point)	1.590	1.273 to 1.986	< 0.001

One point about this engine record: its headline line reports the first predictor in the formula, which is age here, so the headline hazard ratio of 1.011 is the age row and not the ECOG row the question is about. Read all three rows from the coefficient table. Reordering the predictors changes which row is the headline, not any estimate.

The hazard ratio sentence

A hazard ratio is a ratio of event rates at each moment among people still at risk. It is not a ratio of probabilities and it is not a ratio of survival times. For the question asked, a defensible sentence is:

Adjusting for age and sex, each one-point higher ECOG performance score was associated with a hazard of death 1.59 times as high (hazard ratio 1.59, 95% CI 1.27 to 1.99, p < 0.001, 164 events among 227 people).

Three habits keep that sentence honest. Name the unit (“per point”, “per year”) because a hazard ratio without one means nothing. Give the interval, which says which ratios the data are compatible with, from about 1.3 times to about 2 times as high here. And say “associated with”: this is an observational dataset, and adjusting for two other variables does not turn a coefficient into a causal effect.

The age row shows the other side. Its interval, 0.993 to 1.030, includes 1, so the data are compatible with no association with age after adjusting for the other two terms and also with a modest one. The p-value of 0.232 does not show that age is unrelated to survival here, only that these data cannot tell the possibilities apart.

Test each term, not each coefficient

The Wald p-values test one coefficient each, which is right for a numeric predictor and wrong for a factor with several levels, where you want one test for the whole factor. The engine reports a likelihood ratio test per term, which drops the term and compares the two fits. In the hand-written script the equivalent is drop1(fit, test = "Chisq"). The engine run gave 1.447 for age (p = 0.229), 11.44 for sex and 16.62 for ECOG score, each on 1 degree of freedom, and the sex and ECOG p-values were below 0.001. Because every predictor here has a single coefficient, these agree with the Wald results.

Ties and other defaults

Cox’s model assumes event times are all different. Real data have ties, and coxph() needs a way to handle them. The ties argument takes "efron", "breslow" or "exact", and according to the survival manual the Efron approximation is the default: “it is more accurate when dealing with tied death times, and is as efficient computationally”. If no death times are tied, all the methods are equivalent.

This dataset does have ties. In an extra R call, the 164 deaths fell on 138 distinct days: 24 days had more than one death, accounting for 50 of the deaths. Refitting with ties = "breslow" moved the age coefficient from 0.011067 to 0.011041 and the ECOG coefficient from 0.463728 to 0.462947, so the choice barely matters at this level of ties. The point is not to tune it but to know which one you ran, and to state it in a methods section.

Two other defaults are worth knowing. summary() reports 95% intervals unless you change conf.int, as the manual’s usage line shows. And cox.zph() has its own defaults, covered next.

Checking proportional hazards with cox.zph

The model assumes the hazard ratio for each predictor is constant over follow-up. Grambsch and Therneau’s test looks for a trend of each term’s scaled Schoenfeld residuals against time, and cox.zph() implements it. The manual’s usage line is cox.zph(fit, transform="km", terms=TRUE, singledf=FALSE, global=TRUE), so by default time is transformed with the Kaplan-Meier transformation before testing, you get one test per model term, and you get a global test as well. The manual notes that versions before survival 3.0 computed an approximation and later versions compute the actual score test.

Engine run, per term and global:

Term	Chi-square	df	p
Age	0.188	1	0.665
Sex	2.305	1	0.129
ECOG score	2.054	1	0.152
Global	4.464	3	0.216

No term and not the global test gives evidence against proportional hazards, and the engine’s assumption check records the same: “no evidence against proportional hazards”. That is weaker than “the assumption holds”. A test with 164 events can miss a modest drift, and the usual advice is to look at the plot as well. plot(zph) draws the estimated coefficient as a function of transformed time, and the manual says that if proportional hazards holds the true curve is a horizontal line. Judge by whether the whole band could sit around a flat line, and be wary of a conclusion that rests on the right-hand end, where few people remain.

Two cautions apply to any such check. With a large sample, cox.zph() can flag a drift too small to matter, and Stensrud and Hernán argue that a test of proportional hazards is a poor basis for choosing between analyses. What matters more is whether a single hazard ratio is a fair summary of the question you asked. And the three per-term p-values are not corrected for multiplicity (the global line is the single joint test), so treat borderline values as a prompt to look at the plot.

If the check fails

If a term drifts, you have choices. Stratifying on a variable you do not need an estimate for gives each stratum its own baseline hazard. In an extra R call, coxph(Surv(time, status) ~ age + ph.ecog + strata(sex), data = lung) gave a hazard ratio of 1.588 for ECOG score (95% CI 1.268 to 1.988) and 1.011 for age (0.992 to 1.029). Stratifying removes the sex hazard ratio but leaves the others nearly unchanged here, which fits the passed check. Other options are splitting follow-up into periods with a separate ratio in each, a time-varying coefficient, or a summary that does not need proportional hazards, such as the difference in survival at a stated time. The guide to Kaplan-Meier and Cox regression covers when each fits.

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.

Step	Hand-written R	In the engine run
Outcome	Surv(time, status) reads 1/2	The status column is recoded to 0/1 first (2 becomes the event), a recode the run reports as a note
Missing rows	coxph() drops them by default	A complete-case step drops them and reports it: 227 of 228 rows used, 1 dropped
Model	coxph(Surv(time, status) ~ age + sex + ph.ecog)	The identical formula, fitted by the same survival::coxph()
Proportional hazards	cox.zph(fit)	The same call, with its per-term and global results recorded as an assumption check
Term tests	drop1(fit, test = "Chisq")	The same call, giving the likelihood ratio table above
Printing	print() and summary()	A numbered record plus a results file, so a reader can compare

So the difference is not the statistics but what surrounds them. The downloaded script starts by reading the bundled data file and tidying text cells, recodes the status, drops incomplete rows, and then runs the three calls in the table. Everything after the set-up block (reading, cleaning, recoding, dropping incomplete rows) is the code from the previous section. Run it and the numbers in the results file match the recorded ones, which is the check the run performs before a result is accepted. The engine also adds two notes the hand-written script does not: the events-per-variable check (about 55 here, at or above the rule-of-thumb floor, a convention and not a guarantee) and a descriptive note on values beyond the 1.5 times interquartile range fence in ph.ecog (1 value) and time (10 values). Those notes remove nothing. They only flag values worth a second look.

Since the model is the same, you do not need Celeus to follow this post. The point of a recorded run is the audit trail: the seed, the package versions and the exact script travel together, so a reader can rerun your analysis and reach the same numbers. The guide to reproducible statistical analysis explains why that matters beyond one model.

Pitfalls to avoid

A few mistakes come up often when people run coxph() for the first time.

• Wrong status coding. Surv() reads 0/1 as censored/event and 1/2 as censored/event, so a column where 1 means the event and 2 means censored is read backwards and reverses every hazard ratio. Check table(status) against what the codes mean.
• Treating a coded category as a number. sex here is 1 and 2 and works as a single contrast. A variable with three or more categories coded 1, 2, 3 must be a factor, or the model assumes the hazard rises by equal steps.
• Reading the headline as the question. The first row of a table is not special. Report the row for the exposure you asked about.
• Ignoring clustering. People in this dataset were enrolled at several institutions, and a model that treats them as independent understates the standard errors if people at one site resemble each other. cluster() in the formula gives robust standard errors.
• Too few events. The number of events, 164 here, limits how many predictors you can estimate. Count events, not rows.
• Reporting only a p-value. Always give the hazard ratio and its interval.

Where to go next

For the ideas behind the output, read the guide to Kaplan-Meier and Cox regression, and for a longer walkthrough of the same data see the survival analysis tutorial. 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.

This 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: 227
- statistic: 1.194
- p-value: 0.232
- effect size: hazard ratio = 1.011 (95% CI 0.9929 to 1.03)
- 95% CI: [0.9929, 1.03]
- rows used: 227 of 228 (1 dropped)
- seed: 20261023

## Assumption checks
- proportional hazards (cox.zph): passed=TRUE (no evidence against proportional hazards)
- 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 |
|---|---|---|---|---|---|
| age | LR | 1.447 | 1 | - | 0.229 |
| sex | LR | 11.44 | 1 | - | < 0.001 |
| ph.ecog | LR | 16.62 | 1 | - | < 0.001 |

## 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.31.1  |  catalog c39-d3cd
# Research context: not provided for this run.
# Reproduces the recorded result deterministically.

set.seed(20261023)
.data_path <- "data/cox-regression-in-r.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.
data[["status"]] <- as.integer(factor(data[["status"]], levels = c("1", "2"))) - 1L
.acols <- c("age", "sex", "ph.ecog", "time", "status")
data <- data[stats::complete.cases(data[, .acols, drop = FALSE]), , drop = FALSE]
result <- { .fit <- survival::coxph(survival::Surv(`time`, `status`) ~ `age` + `sex` + `ph.ecog`, data = data); list(summary = summary(.fit), ph = survival::cox.zph(.fit), term_tests = stats::drop1(.fit, test = "Chisq")) }
print(result)