coxph в R, бета зависит от значения фактора?

Question

Сейчас я бегу за рулем. Моя установка: у меня есть справка (без лечения), а затем три разных лечения (A, B и C). У меня также есть взаимодействия A, B и C (например, образцы, обработанные как обработкой A и B, так и A и C и т. Д.). Я создал фиктивные переменные для этих обработок, закодированные как 1 или 2 (1 = получал лечение, 2 = не получал лечение). Я использую as.factor() для загрузки этих переменных.

example:
A<-as.factor(Data$A)

Я могу выполнить это следующим образом и получить результат, показывающий, что получение лечения B (он же B = 1) полезно для продолжительности жизни (коэффициент положительный). Все три важны в некотором роде:

> coxph1<-coxph(Surv(Lifespan,Status)~A+B+C
> summary(coxph1)
Call:
coxph(formula = Surv(Life, Status) ~ A + B + C, data = FlyData, 
    method = "efron")

  n= 162, number of events= 140 

     coef exp(coef) se(coef)      z Pr(>|z|)    
A -0.3486    0.7057   0.1761 -1.980 0.047753 *  
B  0.5911    1.8059   0.1787  3.307 0.000944 ***
C -0.6956    0.4988   0.1815 -3.832 0.000127 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

  exp(coef) exp(-coef) lower .95 upper .95
A    0.7057     1.4170    0.4997    0.9966
B    1.8059     0.5537    1.2722    2.5635
C    0.4988     2.0050    0.3494    0.7119

Concordance= 0.822  (se = 0.095 )
Rsquare= 0.227   (max possible= 1 )
Likelihood ratio test= 41.75  on 3 df,   p=5e-09
Wald test            = 41.35  on 3 df,   p=6e-09
Score (logrank) test = 43.6  on 3 df,   p=2e-09

Но когда я запускаю coxph с терминами взаимодействия, где я хочу знать, имеют ли A: B или A: C и т. Д. Какое-то взаимодействие, отличное от просто A или B, я получаю следующее:

> int.coxph <- coxph(Surv(Life, Status)~A*B*C, data=FlyData, method='efron')

Предупреждающее сообщение: В установщике (X, Y, strats, offset, init, control, weights = weights,: Логлик сходился до переменной 1,2,3,4,5,6,7; бета может быть бесконечной.

> summary(int.coxph)
Call:
coxph(formula = Surv(Life, Status) ~ A * B * C, data = FlyData, 
    method = "efron")

  n= 162, number of events= 140 

            coef  exp(coef)   se(coef)      z Pr(>|z|)
A      3.987e+01  2.066e+17  4.945e+03  0.008    0.994
B      1.856e+01  1.148e+08  2.472e+03  0.008    0.994
C      3.799e+01  3.144e+16  4.945e+03  0.008    0.994
A:B   -1.964e+01  2.967e-09  2.472e+03 -0.008    0.994
A:C   -3.954e+01  6.737e-18  4.945e+03 -0.008    0.994
B:C   -1.874e+01  7.241e-09  2.472e+03 -0.008    0.994
A:B:C  1.962e+01  3.318e+08  2.472e+03  0.008    0.994

      exp(coef) exp(-coef) lower .95 upper .95
A     2.066e+17  4.841e-18         0       Inf
B     1.148e+08  8.714e-09         0       Inf
C     3.144e+16  3.180e-17         0       Inf
A:B   2.967e-09  3.370e+08         0       Inf
A:C   6.737e-18  1.484e+17         0       Inf
B:C   7.241e-09  1.381e+08         0       Inf
A:B:C 3.318e+08  3.014e-09         0       Inf

Concordance= 0.869  (se = 0.095 )
Rsquare= 0.51   (max possible= 1 )
Likelihood ratio test= 115.6  on 7 df,   p=<2e-16
Wald test            = 9.24  on 7 df,   p=0.2
Score (logrank) test = 73.69  on 7 df,   p=3e-13

Итак ... это похоже на некоторые другие вопросы ... но почему бета приближается к бесконечной? Дополнительный поворот, который у меня есть для этого вопроса, состоит в том, что если я перекодирую переменные как 0 или 1 (вместо 1 и 2), то я могу изменить вывод в coxph () взаимодействия. Это перекодирование для coxph:

coxph2<-coxph(Surv(Lifespan, Status)~A2+B2+C2))
summary(coxph2)
Call:
coxph(formula = Surv(Life, Status) ~ A2 + B2 + C2, data = FlyData, 
    method = "efron")

  n= 162, number of events= 140 

      coef exp(coef) se(coef)      z Pr(>|z|)    
A2  0.3486    1.4170   0.1761  1.980 0.047753 *  
B2 -0.5911    0.5537   0.1787 -3.307 0.000944 ***
C2  0.6956    2.0050   0.1815  3.832 0.000127 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

   exp(coef) exp(-coef) lower .95 upper .95
A2    1.4170     0.7057    1.0035     2.001
B2    0.5537     1.8059    0.3901     0.786
C2    2.0050     0.4988    1.4048     2.862

Concordance= 0.822  (se = 0.095 )
Rsquare= 0.227   (max possible= 1 )
Likelihood ratio test= 41.75  on 3 df,   p=5e-09
Wald test            = 41.35  on 3 df,   p=6e-09
Score (logrank) test = 43.6  on 3 df,   p=2e-09

это просто обратное, но таз взаимодействия другой ...

> full.coxph <- coxph(Surv(Life, Status)~A2*B2*C2, data=FlyData, method='efron')
Warning message:
In fitter(X, Y, strats, offset, init, control, weights = weights,  :
  Loglik converged before variable  2,4,6,7 ; beta may be infinite. 
> summary(full.coxph)
Call:
coxph(formula = Surv(Life, Status) ~ A2 * B2 * C2, data = FlyData, 
    method = "efron")

  n= 162, number of events= 140 

               coef  exp(coef)   se(coef)      z Pr(>|z|)
A2       -7.067e-15  1.000e+00  3.204e-01  0.000    1.000
B2       -2.028e+01  1.558e-09  2.472e+03 -0.008    0.993
C2        9.821e-02  1.103e+00  3.204e-01  0.307    0.759
A2:B2     1.960e+01  3.266e+08  2.472e+03  0.008    0.994
A2:C2    -2.991e-01  7.415e-01  4.475e-01 -0.668    0.504
B2:C2     2.050e+01  7.970e+08  2.472e+03  0.008    0.993
A2:B2:C2 -1.962e+01  3.014e-09  2.472e+03 -0.008    0.994

         exp(coef) exp(-coef) lower .95 upper .95
A2       1.000e+00  1.000e+00    0.5337     1.874
B2       1.558e-09  6.417e+08    0.0000       Inf
C2       1.103e+00  9.065e-01    0.5888     2.067
A2:B2    3.266e+08  3.062e-09    0.0000       Inf
A2:C2    7.415e-01  1.349e+00    0.3085     1.782
B2:C2    7.970e+08  1.255e-09    0.0000       Inf
A2:B2:C2 3.014e-09  3.318e+08    0.0000       Inf

Concordance= 0.869  (se = 0.095 )
Rsquare= 0.51   (max possible= 1 )
Likelihood ratio test= 115.6  on 7 df,   p=<2e-16
Wald test            = 9.24  on 7 df,   p=0.2
Score (logrank) test = 73.69  on 7 df,   p=3e-13

Зачем менять числовое значение категориальной переменной? : S Что мне здесь не хватает ... Повторная попытка с нечисловыми переменными («нет» и «да») дает тот же результат, что и при использовании 0 и 1. Например, верхний 0,95 для A равен «1,874», для B - «inf». Точно так же, coxph(Surv()~A+B+C) дает отрицательный коэффициент для B, как и выше.

Answer 1

Я осознал одну из проблем, вызвавших мою проблему: в моих данных о выживании просто не было достаточного разрешения. Я не мог дифференцировать эффекты условий взаимодействия. Если я спроектирую свои данные, чтобы получить ответ, то я смогу получить разумные результаты загрузки модели и условия взаимодействия, которые имеют смысл. В конце дня я планирую использовать комбинированный подход всех трех типов моделей. i.e.:

coxph(Surv(Time, Status)~A+B+C, data=data) #Additive effects
coxph(Surv(Time, Status)~Treatment, data=data) #Base treatment effects
coxph(Surv(Time, Status)~A+B+A:B, data=data) #Test interactions of interest

Базовое понимание аддитивных эффектов дает вам представление о том, как ковариаты глобально способствуют выживанию. Анализ эффектов лечения (т. Е. Фундаментальная переменная интереса) дает вам представление о том, отличаются ли группы, и из этого вы можете вывести закономерности, используя аддитивные эффекты и переменные интереса.

Использование 42-х подходов к исследованию только терминов, представляющих интерес, также весьма полезно при анализе данных. Независимо от того, как я манипулирую данными, даже у данных, которые я разработал для информативности, возникают проблемы, когда вы включаете все термины взаимодействия в трехстороннюю модель. Но использование только взаимодействий, представляющих интерес, может улучшить понимание.

Я полагаю, что такого рода последующий анализ требует независимой проверки из второго эксперимента, сфокусированного на условиях интереса.

Answer 2

У вас, вероятно (почти наверняка на самом деле) есть почти вырожденная "матрица шляпы", которая образуется из матрицы модели с этим взаимодействием. У вас есть все взаимодействия второго порядка, а также взаимодействия третьего порядка. В зависимости от количества уровней в факторах количество терминов, необходимых для полного заполнения матрицы модели, может быть очень большим. Далее я попробую модель с чуть меньшим количеством терминов в модели. Вы можете использовать интерфейс формулы R, чтобы удалить термины третьего порядка и оставить первый и второй члены только одним из двух способов:

int.coxph <- coxph(Surv(Life, Status)~( A+B+C)^2, data=FlyData, method='efron')

Или:

int.coxph <- coxph(Surv(Life, Status)~ A*B*C - A:B:C, data=FlyData, method='efron')

Не уверен, что таким образом вы получите удовлетворение. Вполне возможно, что у вас недостаточно данных, чтобы избежать вырожденности при построении XX ^ t-матрицы, но если ваши результаты не будут взорваны таким очевидным образом, как показано выше, то, возможно, результаты будут значимыми. Другим более безопасным методом было бы сначала посмотреть на сокращенную модель, а затем добавить обратно в конкретные взаимодействия:

 int.coxph.base <- coxph(Surv(Life, Status)~A+B+C,      data=FlyData, method='efron')
int.coxph.intAB <- coxph(Surv(Life, Status)~A+B+C +A:B, data=FlyData, method='efron')

Этот второй вариант имеет дополнительное преимущество, заключающееся в том, что вы можете легко создавать тесты на основе изменения вероятности записи в журнале, а не в зависимости от менее надежных тестов типа Вальда, которые вы видите в распечатках по умолчанию для print.coxph или summary.coxph.