How to perform a stratified log rank test in R

Question

[ykher92]

Suppose I have two matched sets with n = 50 each. I've arranged them by an ID variable such that each ID variable has 2 subjects. I'd like to compare overall survival with a kaplan meier accounting for their paired nature. I understand the best way to do this is through a stratified log rank test. What is the best way to do this in R? The two methods I've tried:

1) Using the survdiff function along the lines of survdiff(Surv(follow up, event) ~ variable + cluster(ID), data = dataframe)

2) Using the coxPH function along the lines of coxph(Surv(follow up, event) ~ variable + cluster(ID), data = dataframe)

However, I'm getting drastically different results using these two methods. Using a standard non-stratified log-rank test with my data, I get a p value of ~ 0.7. With method 2 (coxPH), I get a similar p value of ~ 0.7 albeit with a different robust standard error. Using method 1 (survdiff), my p value is DRASTICALLY different (p <0.0001). Why is this? Does clustering not work in the survdiff function?

[zy]

I wonder how to perform the stratified log-rank test in R？ Have you got the data or corresponding functions?

[GenoMax]

Code Kevin Blighe posted in his answer is for R.

Answer 1

[Kevin Blighe]

Hey,

No, they should be the same. There are many different survival functions; so, it can become a bit confusing. Here, I'll re-use some code from a previous couples of answers:

• A: survfit(Surv()) P-value interpretation for 3 survival curves?
• A: R survival analysis : surv_pvalue vs fit.coxph for log-rank-test pvalue

1, load AML data

require(survival)
require(survminer)

2, add a third stratum 'SuperMaintained' (log rank test can compare >2 strata / levels)

aml$x <- as.character(aml$x)
aml[10,3] <- 'SuperMaintained'
aml[11,3] <- 'SuperMaintained'
aml[22,3] <- 'SuperMaintained'
aml[23,3] <- 'SuperMaintained'
aml$x <- factor(aml$x, levels = c('Nonmaintained','Maintained','SuperMaintained'))
aml

   time status               x
1     9      1      Maintained
2    13      1      Maintained
3    13      0      Maintained
4    18      1      Maintained
5    23      1      Maintained
6    28      0      Maintained
7    31      1      Maintained
8    34      1      Maintained
9    45      0      Maintained
10   48      1 SuperMaintained
11  161      0 SuperMaintained
12    5      1   Nonmaintained
13    5      1   Nonmaintained
14    8      1   Nonmaintained
15    8      1   Nonmaintained
16   12      1   Nonmaintained
17   16      0   Nonmaintained
18   23      1   Nonmaintained
19   27      1   Nonmaintained
20   30      1   Nonmaintained
21   33      1   Nonmaintained
22   43      1 SuperMaintained
23   45      1 SuperMaintained

3, obtain survdiff log rank p-value

survfit <- survdiff(Surv(time, status) ~ x, data = aml)
1 - pchisq(survfit$chisq, length(survfit$n) - 1)
[1] 0.005309417

4, obtain coxph log rank p-value

coxfit <- coxph(
  Surv(time, status) ~ x,
  data = aml,
  ties = 'exact')

summary(coxfit)

round(summary(coxfit)$sctest[3], digits = 9)

     pvalue 
0.005309417

So, you can see that they are the same.

Kevin