1
votes

I have a survival dataset. I would like to run a logrank test for a treatment, broken up into 4 categories. I cannot use the survdiff() command because the asymptotic distribution of those statistics is a chi square, and I need normality (I am doing this in the multiple imputation setting and pooling later). Instead, I want to run a Cox regression and then run the score test, which will be normally distributed.

So, what I would like to do is to take my 4 categories, and then break them down into a few groups to compare them individually. For example

Treatment 2 vs Treatment 3: Is it possible to do this without breaking down the data? Suppose we have the burn dataset from package KMsurv

library(KMsurv)

> summary(coxph(Surv(T1,D1)~factor(Z11),data=burn))
Call:
coxph(formula = Surv(T1, D1) ~ factor(Z11), data = burn)

  n= 154, number of events= 99 

            coef exp(coef) se(coef)      z Pr(>|z|)  
  factor(Z11)2 -0.9820    0.3745   0.4956 -1.982   0.0475 *
  factor(Z11)3 -1.6872    0.1850   0.8029 -2.101   0.0356 *
  factor(Z11)4 -0.4070    0.6656   0.3957 -1.029   0.3037  
 ...
 Likelihood ratio test= 9.17  on 3 df,   p=0.0271
 Wald test            = 7.38  on 3 df,   p=0.06083
 Score (logrank) test = 8  on 3 df,   p=0.04602

This outputs the logrank test for 1 vs 2 vs 3 vs 4, but I only want 2 vs 3. I know I could get it by running before this command

subsetted=subset(burn,Z11==2|Z11==3)
summary(coxph(Surv(T1,D1)~factor(Z11),data=subsetted))

But this will get tedious and hard to debug when we have to do things like compare 1,2 vs 4

So, is there any way to select what groups you want to compare in the coxph command, or is the only way to select groups to presubset them before?

1
The correct method for comparing 1 and 2 together versus 3 would be with contrasts. I suggest using a score statistic when a LR test is available is a faulty approach. The fact that there is "normal theory" strategy for handling it is not a convincing argument for using an inferior statistical test. - IRTFM
Without normality, there is no way for me to pool my results (that I know of, I am currently pooling via Rubin's rules). So I feel that it is kind of a give and take. We lose some optimality in the test, but we gain the ability to pool the results across many imputed datasets. Do you know transformation that would make the LR test normally distributed? I know that it is chi squared, but the square root will be chi distributed, not normal. - RayVelcoro
As I read sites.stat.psu.edu/~jls/mifaq.html#howto the operations are done on the individual estimates and their standard errors, not with a global measure of fit. The parameter estimates from survival analysis are presumed to be normally distributed. - IRTFM

1 Answers

2
votes

Use the subset argument in the coxph function. See ?coxph

But you also need to drop unwanted factor levels from burn$Z11

So you could do

summary(coxph(Surv(T1,D1)~factor(Z11,levels=c('2','3')),data=burn, subset=Z11 %in% c('2','3')))

or bit more convenient like

mylevels <- c('2','3')  #specify factor levels for subset
summary(coxph(Surv(T1,D1)~factor(Z11,levels=mylevels),data=burn, subset=Z11 %in% mylevels))

also see this thread to make a wrapper where subset is defined as argument in wrapper function:

https://stat.ethz.ch/pipermail/r-help/2007-November/145345.html