Логранговый тест для определенных групп с помощью регрессии Кокса

У меня есть набор данных выживания. Я хотел бы провести логарифмический тест для лечения, разбитого на 4 категории. Я не могу использовать команду survdiff(), потому что асимптотическое распределение этих статистических данных представляет собой хи-квадрат, и мне нужна нормальность (я делаю это в настройке множественного вменения и объединения позже). Вместо этого я хочу запустить регрессию Кокса, а затем запустить тест оценки, который будет нормально распределен.

Итак, что я хотел бы сделать, так это взять мои 4 категории, а затем разбить их на несколько групп, чтобы сравнить их по отдельности. Например

Обработка 2 против обработки 3: возможно ли сделать это без разбивки данных? Предположим, у нас есть набор данных записи из пакета 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

Это выводит тест logrank для 1 против 2 против 3 против 4, но я хочу только 2 против 3. Я знаю, что могу получить его, выполнив перед этой командой

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

Но это станет утомительным и трудным для отладки, когда нам придется делать такие вещи, как сравнение 1,2 и 4.

Итак, есть ли какой-либо способ выбрать, какие группы вы хотите сравнить в команде coxph, или это единственный способ выбрать группы, чтобы предварительно задать их?


person RayVelcoro    schedule 29.09.2015    source источник
comment
Правильным методом сравнения 1 и 2 вместе с 3 будет использование контрастов. Я предлагаю использовать статистику баллов, когда доступен тест LR, что является ошибочным подходом. Тот факт, что для этого существует нормальная теоретическая стратегия, не является убедительным аргументом в пользу использования некачественного статистического теста.   -  person IRTFM    schedule 30.09.2015
comment
Без нормальности я не могу объединить свои результаты (насколько я знаю, я в настоящее время объединяю по правилам Рубина). Так что я чувствую, что это своего рода давать и брать. Мы теряем некоторую оптимальность в тесте, но получаем возможность объединять результаты по многим вмененным наборам данных. Знаете ли вы преобразование, которое сделало бы тест LR нормально распределенным? Я знаю, что это хи в квадрате, но квадратный корень будет хи распределенным, а не нормальным.   -  person RayVelcoro    schedule 30.09.2015
comment
Когда я прочитал sites.stat.psu.edu/~jls/mifaq.html #howto операции выполняются с отдельными оценками и их стандартными ошибками, а не с глобальной мерой соответствия. Предполагается, что оценки параметров из анализа выживаемости имеют нормальное распределение.   -  person IRTFM    schedule 30.09.2015


Ответы (1)


Используйте аргумент subset в функции coxph. См. ?coxph

Но вам также необходимо исключить нежелательные уровни факторов из burn$Z11

Так что вы могли бы сделать

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

или немного удобнее, как

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

также см. этот поток, чтобы создать оболочку, где подмножество определяется как аргумент в функции-оболочке:

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

person Chris Holbrook    schedule 29.09.2015
comment
Я думаю, что это почти отвечает на него. Что бы я сделал, если бы хотел сравнить группу 1 с группой 2 и 3? - person RayVelcoro; 29.09.2015
comment
mylevels <- c('1','2','3')? или два отдельных сравнения, например. mylevels <- c('1','2') и mylevels <- c('1','3')? Я не уверен, что понимаю, что именно вы спрашиваете. - person Chris Holbrook; 30.09.2015