Cox model: extremely high Schoenfeld residuals

Question

I have run fixed in time Cox models (crude and covariate adjusted) for the effect of a binary exposure. Sample size is approx. 200,000 and there are approx. 10,000 events and 4000 "yes" values for the exposure.

I am concerned and confused by the very large values of the scaled Schoenfeld residuals, please see plot below from the crude model. The very high (50+ in magnitude) values appear to be from both exposed (5) and non-exposed (50) subjects.

I fitted the model using the coxph() function from the survival package in R. I used the cox.zph() function with identity transform from the same package to produce the plot shown.

UPDATES

library(survival)
library(tidyverse)
fit model
f1 <- coxph(formula = Surv(tstart, time, event) ~ exposed_bin,
    data = test_df,
    robust = TRUE,
    cluster = id1,
    x = TRUE)
check ph
ph_check <- cox.zph(f1)
show some values of the Schoenfeld resids
test_df %>% filter(event == 1) %>% arrange(time) %>% 
 select(id1, id2, exposed_bin, time, event) %>%
 bind_cols(scho_resid = ph_check$y) %>% 
 slice(176:200) %>% 
 print(n = 25)

plot of ph_check with KM transform and no residuals displayed...

Answer 1

I am concerned and confused by the very large values of the scaled Schoenfeld residuals...(emphasis added)

The scaling of Schoenfeld residuals is proportional to the number of events in the dataset. See this page, for example, or this snippet of code from the residuals.coxph() function for the case of scaled Schoenfeld residuals:

if (otype == "scaledsch") {
            ndead <- sum(deaths)
            coef <- ifelse(is.na(object$coefficients), 0, object$coefficients)
            if (nvar == 1) 
                rr <- rr * vv * ndead + coef
            else {
                cname <- colnames(rr)
                rr <- drop(rr %*% vv * ndead + rep(coef, each = nrow(rr)))
                colnames(rr) <- cname
            }
        }

Here, the initial value of rr is the set of unscaled Schoenfeld residuals, and vv is the (non-robust) variance-covariance matrix of the coefficient estimates.

With a large number of deaths, as you have, the absolute values of scaled Schoenfeld residuals can be quite large. What's important is the overall shape of the smoothed plot over time.