0
$\begingroup$

I am running a survival cox model with multiple measures, here is the code:

library(survival)
# Data Simulation
N <- 250
dat <- data.frame(ID = factor(1:N), age = rnorm(N, mean = 45, sd = 5), sex = sample(0:1, 
                                                                                    N, TRUE), basemort = rnorm(N, sd = 3))

interval <- matrix(sample(2:14, N * 3, replace = TRUE), N)
windows <- t(apply(cbind(0, interval), 1, cumsum))
windows <- rbind(windows[, 1:2], windows[, 2:3], windows[, 3:4])

colnames(windows) <- c("time1", "time2")
dat <- cbind(do.call(rbind, rep(list(dat), 3)), windows)
dat <- dat[order(dat$ID), ]
dat$assessment <- rep(1:3, N)
rownames(dat) <- NULL
mortality <- with(dat, {
 mu <- basemort + (0.05 * age) - (2.5 * sex) + (0.3 * time2)
 lp <- rnorm(N * 3, mean = mu, sd = 1)
 as.integer(lp > median(lp))
})
mortality <- as.integer(ave(mortality, dat$ID, FUN = cumsum) >= 1)
mortality[dat$assessment == 1] <- 0
dat$mortality <- mortality

# Survival Model
m <- coxph(Surv(time1, time2, mortality) ~ age + sex, data = dat)
## summary of the model
summary(m)

The output is

Call:
coxph(formula = Surv(time1, time2, mortality) ~ age + sex, data = dat)

  n= 750, number of events= 313 

        coef exp(coef) se(coef)      z Pr(>|z|)    
age  0.01027   1.01032  0.01177  0.872    0.383    
sex -0.48694   0.61451  0.11870 -4.102 4.09e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

    exp(coef) exp(-coef) lower .95 upper .95
age    1.0103     0.9898    0.9873    1.0339
sex    0.6145     1.6273    0.4870    0.7755

Concordance= 0.586  (se = 0.017 )
Likelihood ratio test= 18.23  on 2 df,   p=1e-04
Wald test            = 17.61  on 2 df,   p=1e-04
Score (logrank) test = 17.94  on 2 df,   p=1e-04

The number of subjects is 250, but as it is seen in the output section, it says n= 750, I assume that cox with considering start time and end time takes subjects into consideration, but how can I prove this, if it was STATA, it would show that n=250 instead of 750. I reported this my supervisor and he said my model was wrong and the output should show n=250.

$\endgroup$
2
  • $\begingroup$ This is more about specific software rather than a statistical issue, so it's probably off-topic on this site. See Section 3.2 of the main survival vignette. The cgd data set example there shows the same behavior: 128 unique individuals, but n=203 in the summary, the total number of rows in the data frame. You can use survcheck() on your model to get the number of unique identifiers. $\endgroup$
    – EdM
    Nov 12 at 21:59
  • $\begingroup$ oh wow! Thank you so much. Survcheck is exactly what I was looking for. But how can I incorporate survcheck in my survival model. $\endgroup$
    – Katie
    Nov 12 at 22:08
3
$\begingroup$

The output shows $n=750$ simply because there are 750 observations. The dataframe is treating each of your three assessments of the same ID as an independent observation, so it doesn't matter that you labeled them with ID.

1. As @EdM mentioned in the comments, survcheck() is what you should use to prove that you have the correct number of subjects.

2. You should likely modify your formula to include cluster(ID), the model will then also give you a robust standard error that takes the lack of independence within ID into account. Without this, it assumes that each row is independent.

Adding in these two recommendations, the last portion of your code should change to:

# Survival Model
m <- coxph(Surv(time1, time2, mortality) ~ age + sex + 
       cluster(ID), data = dat)

# Check for correct subject number, events, and censoring
survcheck(Surv(time1, time2, mortality) ~ age + sex + 
          cluster(ID),
          data=dat,
          id=ID)

# Summary of model
summary(m)
$\endgroup$
2
  • $\begingroup$ Assuming that mortality can only occur in at most one of the records for any individual, it's not necessary to use cluster(ID), though it's mostly harmless to add it. cluster(ID) is needed for recurrent events in one individual or correlated events between individuals. $\endgroup$ Nov 18 at 2:34
  • $\begingroup$ @ThomasLumley Mortality can only occur once in reality, but the way that this dataset was constructed forces 3 assessments regardless of whether mortality already occurred. As such, mortality frequently occurs multiple times for the same individual. In reality, I doubt that the experimenter would continue assessing corpses, but for training purposes it is worthwhile to point to the cluster(ID) addendum to the formula for use in real situations where the survival outcome is less... permanent. $\endgroup$
    – AJV
    Nov 18 at 2:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.