Как правильно подобрать регрессию Кокса с помощью coxphf? (Метод уменьшения погрешности максимального правдоподобия Фирта в R)

Я пытаюсь подобрать регрессию Кокса (требуется пакет: выживание) и столкнулся с проблемой. Когда я пытаюсь сопоставить «обычную» регрессионную модель Кокса с моими данными, я получаю сообщение об ошибке «Матрица X считается единственной; переменная 9» (и если я удаляю переменную 9, проблема становится переменной 8). Насколько я понимаю проблему, это происходит потому, что слишком много пациентов с этими переменными имеют одно и то же событие (кажется, в другом вопросе это называлось «идеальной классификацией»)

Вот почему я попытался сопоставить модель Кокса с функцией coxphf (одноименный пакет), так как это должно решить проблему, используя «метод уменьшения смещения максимального правдоподобия Фирта со штрафом» для регрессии Кокса. Но это также не работает, пока я не увеличу «максимум» с 50 по умолчанию до 1000 и не удалю из уравнения «неопределенную» переменную. Но если я удалю неопределенную переменную из своего набора данных (это всего лишь 1 человек), модель снова не будет работать.

Итак, мой вопрос: как я могу это решить? Уместно/необходимо ли вообще удалить всю переменную (следовательно, 1 человека) из набора данных? Возможно, я сделал несколько очень очевидных ошибок, но, пожалуйста, потерпите меня, так как у меня нет абсолютно никакого опыта в статистике. Спасибо большое уже заранее.

Я включил следующие образцы данных, а также свои попытки подогнать модель Кокса. Вот как мне удалось заставить модель работать, исключив из модели переменную «Неопределенная»:

# load packages
library("survival")
library ("coxphf")

example<-structure(list(Pat.nr. = c(1L, 2L, 5L, 7L, 8L, 10L, 13L, 14L, 
15L, 17L, 18L, 19L, 20L, 21L, 22L, 23L, 26L, 28L, 29L, 30L, 31L, 
32L, 33L, 34L, 35L, 36L, 37L, 39L, 41L, 42L, 44L, 45L, 46L, 47L, 
48L, 49L, 50L, 52L, 53L, 56L, 57L, 58L, 59L, 60L, 61L, 62L, 63L, 
67L, 68L, 69L), OS = c(1.6, 34.6, 1.5, 35.8, 7.7, 38.6, 37.6, 
8.6, 0.6, 5.7, 0.6, 43.9, 25.8, 7.3, 28.1, 43.8, 12.8, 18.5, 
36.1, 43.1, 15.4, 37.6, 8.6, 2.7, 10.2, 8.1, 37.3, 25.3, 3.7, 
26.1, 41.2, 5.9, 15.5, 56.8, 29.5, 52.1, 5.4, 54.8, 53.5, 16.6, 
49.2, 53.8, 8.5, 56, 7.4, 28, 3.3, 38, 55.7, 0.4), Event = c(1L, 
0L, 0L, 0L, 1L, 0L, 0L, 1L, 1L, 1L, 1L, 0L, 0L, 1L, 1L, 0L, 1L, 
0L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 0L, 1L, 0L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 0L, 0L, 1L, 0L, 0L, 1L, 0L, 1L, 0L, 1L, 0L, 1L, 
1L), Age = c(68.41, 54.9, 44.44, 64.14, 68.86, 62.93, 40.76, 
31.06, 42.97, 69.16, 47.39, 60.14, 27.9, 56.57, 19.63, 47.75, 
45.58, 66.22, 43.73, 45.34, 38.83, 54.46, 48.91, 70.3, 60.51, 
68.55, 63.18, 55.89, 68.27, 57.25, 56.17, 60.83, 74.42, 71.3, 
40.36, 50.85, 59.61, 50.14, 45.77, 19.34, 56.32, 53.38, 70.7, 
55.25, 56.05, 44.06, 51.36, 69.37, 69.71, 75.44), Favorable = c(0L, 
0L, 1L, 0L, 0L, 1L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 
0L, 0L, 1L, 1L, 0L, 1L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 
0L), Intermediate = c(0L, 2L, 0L, 0L, 2L, 0L, 0L, 0L, 0L, 2L, 
0L, 2L, 2L, 0L, 2L, 0L, 2L, 2L, 0L, 0L, 0L, 0L, 2L, 0L, 0L, 0L, 
2L, 0L, 0L, 2L, 2L, 2L, 0L, 2L, 2L, 0L, 0L, 2L, 0L, 2L, 0L, 0L, 
0L, 2L, 2L, 0L, 2L, 0L, 2L, 2L), Adverse = c(1L, 0L, 0L, 1L, 
0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 1L, 0L, 
1L, 0L, 0L, 1L, 1L, 1L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 1L, 0L, 0L, 0L, 0L), Undefined = c(0L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L)), .Names = c("Pat.nr.", "OS", "Event", "Age", "Favorable", 
"Intermediate", "Adverse", "Undefined"), row.names = c(NA, 50L
), class = "data.frame")


#row&columns
n_row <- dim(example)[1]
n_col <- dim(example)[2]

#variables:
OS <- c(example[1:n_row,2])
Event <- c(example[1:n_row,3])
age <- c(example[1:n_row,4])
Favorable <- c(example[1:n_row,5])
Intermediate <-  c(example[1:n_row,6])
Adverse <- c(example[1:n_row,7])
Undefined <-  c(example[1:n_row,8])


# dependent and independent variables
y <- Surv(OS, Event)
x <- cbind(age, Favorable, Intermediate, Adverse, Undefined)  
example <- data.frame(cbind(x,y)) 

#coxphf with Firth's Penalized Likelihood --> which doesn't seem to work
cox2<-coxphf(data=example,y~x, firth=TRUE, pl=TRUE, maxit=1000)
summary(cox2)

#coxphf with Firth's Penalized Likelihood (without Variable "Undefined") --> this works
cox2<-coxphf(data=example,y~age+Favorable+Intermediate+Adverse, firth=TRUE, pl=TRUE, maxit=1000)
summary(cox2)

И здесь я изменил набор данных, чтобы он не включал неопределенную переменную (и модель больше не работает):

example1<-structure(list(Pat.nr. = c(1L, 2L, 5L, 7L, 8L, 10L, 13L, 14L, 
15L, 17L, 19L, 20L, 21L, 22L, 23L, 26L, 28L, 29L, 30L, 31L, 32L, 
33L, 34L, 35L, 36L, 37L, 39L, 41L, 42L, 44L, 45L, 46L, 47L, 48L, 
49L, 50L, 52L, 53L, 56L, 57L, 58L, 59L, 60L, 61L, 62L, 63L, 67L, 
68L, 69L, 72L), OS = c(1.6, 34.6, 1.5, 35.8, 7.7, 38.6, 37.6, 
8.6, 0.6, 5.7, 43.9, 25.8, 7.3, 28.1, 43.8, 12.8, 18.5, 36.1, 
43.1, 15.4, 37.6, 8.6, 2.7, 10.2, 8.1, 37.3, 25.3, 3.7, 26.1, 
41.2, 5.9, 15.5, 56.8, 29.5, 52.1, 5.4, 54.8, 53.5, 16.6, 49.2, 
53.8, 8.5, 56, 7.4, 28, 3.3, 38, 55.7, 0.4, 2.8), Event = c(1L, 
0L, 0L, 0L, 1L, 0L, 0L, 1L, 1L, 1L, 0L, 0L, 1L, 1L, 0L, 1L, 0L, 
0L, 0L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 0L, 1L, 0L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 0L, 0L, 1L, 0L, 0L, 1L, 0L, 1L, 0L, 1L, 0L, 1L, 1L, 
1L), Age = c(68.41, 54.9, 44.44, 64.14, 68.86, 62.93, 40.76, 
31.06, 42.97, 69.16, 60.14, 27.9, 56.57, 19.63, 47.75, 45.58, 
66.22, 43.73, 45.34, 38.83, 54.46, 48.91, 70.3, 60.51, 68.55, 
63.18, 55.89, 68.27, 57.25, 56.17, 60.83, 74.42, 71.3, 40.36, 
50.85, 59.61, 50.14, 45.77, 19.34, 56.32, 53.38, 70.7, 55.25, 
56.05, 44.06, 51.36, 69.37, 69.71, 75.44, 71.05), Favorable = c(0L, 
0L, 1L, 0L, 0L, 1L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 1L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 0L, 
0L, 1L, 1L, 0L, 1L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 
1L), Intermediate = c(0L, 2L, 0L, 0L, 2L, 0L, 0L, 0L, 0L, 2L, 
2L, 2L, 0L, 2L, 0L, 2L, 2L, 0L, 0L, 0L, 0L, 2L, 0L, 0L, 0L, 2L, 
0L, 0L, 2L, 2L, 2L, 0L, 2L, 2L, 0L, 0L, 2L, 0L, 2L, 0L, 0L, 0L, 
2L, 2L, 0L, 2L, 0L, 2L, 2L, 0L), Adverse = c(1L, 0L, 0L, 1L, 
0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 1L, 0L, 1L, 
0L, 0L, 1L, 1L, 1L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L)), .Names = c("Pat.nr.", 
"OS", "Event", "Age", "Favorable", "Intermediate", "Adverse"), row.names = c(NA, 
50L), class = "data.frame")

#row&columns
n_row <- dim(example1)[1]
n_col <- dim(example1)[2]

#variables:
OS <- c(example1[1:n_row,2])
Event <- c(example1[1:n_row,3])
age <- c(example1[1:n_row,4])
Favorable <- c(example1[1:n_row,5])
Intermediate <-  c(example1[1:n_row,6])
Adverse <- c(example1[1:n_row,7])




# dependent and independent variables
y <- Surv(OS, Event)
x <- cbind(age, Favorable, Intermediate, Adverse)  
example <- data.frame(cbind(x,y)) 

# dependent and independent variables
y <- Surv(OS, Event)
x <- cbind(age, Favorable, Intermediate, Adverse)  
example1 <- data.frame(cbind(x,y)) 

#coxphf with Firth's Penalized Likelihood

cox2<-coxphf(data=example,y~x, firth=TRUE, pl=TRUE, maxit=1000)
summary(cox2)

r survival-analysis cox-regression multivariate-testing
person Alex    schedule 31.01.2018    source источник

Ответы (1)


arrow_upward
0
arrow_downward

Я пошел немного по-другому с вашими данными, и у меня не было никаких проблем. Прежде всего, я исключил столбец Undefined, так как переменная без каких-либо вариаций ничего не делает для вашей модели (и будет создавать ошибки).

Затем я просто привязал ваш объект Surv() к предоставленным вами данным:

example$survive <- Surv(time = OS, event = Event)

и вызвал coxphf() с настройками по умолчанию (поскольку firth и pl по умолчанию являются TRUE), не изменяя итерацию по умолчанию 50, и использовал правильную спецификацию формулы:

cox2 <- coxphf(survive ~ Age + Favorable + Intermediate + Adverse, data = example)

> summary(cox2)
coxphf(formula = survive ~ Age + Favorable + Intermediate + Adverse, 
    data = example)

Model fitted by Penalized ML
Confidence intervals and p-values by Profile Likelihood 

                     coef   se(coef)  exp(coef)  lower 0.95 upper 0.95     Chisq          p
Age           0.005652315 0.01598045 1.00566832 0.976272961  1.0391381 0.1297188 0.71872373
Favorable    -3.335329337 1.25086821 0.03560286 0.004236556  0.4169422 6.3231629 0.01191709
Intermediate -1.672995547 0.61363779 0.18768401 0.066460220  0.6321023 6.4756709 0.01093610
Adverse      -3.412597661 1.25225958 0.03295548 0.003909284  0.3863706 6.4164893 0.01130655

Likelihood ratio test=6.535659 on 4 df, p=0.1625574, n=50
Wald test = 7.794526 on 4 df, p = 0.09940164

Covariance-Matrix:
                       Age   Favorable Intermediate      Adverse
Age           0.0002553747 -0.00220437 -0.001027825 -0.001760528
Favorable    -0.0022043699  1.56467129  0.710597909  1.408715695
Intermediate -0.0010278246  0.71059791  0.376551335  0.705180703
Adverse      -0.0017605280  1.40871570  0.705180703  1.568154050

Могу я спросить, почему столбец Intermediate имеет значения 2 вместо 1? Разве это не просто бинарный индикатор категориальной переменной?

person LAP    schedule 31.01.2018
comment
Благодарю вас! Да, конечно, вы можете спросить, я ошибся при построении выборочного набора данных (поскольку я разделил группы риска, которые ранее были сведены в 1 переменную). Итак, технически можно просто исключить переменную из модели без изменения набора данных? (Это не создает никаких ошибок, так как будет один человек, который не принадлежит ни к одной из трех оставшихся категорий риска?) И что именно делает привязка объекта Surv() к данным? - person Alex; 31.01.2018
comment
Привязка объекта к данным позволяет использовать правильный формат ввода formula, data функции coxphf. Вы можете попробовать исключить один случай без какого-либо наблюдения за одной из трех переменных, чтобы увидеть, есть ли разница в результатах. Как правило, это вполне приемлемо. - person LAP; 31.01.2018
comment
Спасибо, это имеет смысл. Во втором примере (набор данных: пример 1) я попытался исключить этот случай без каких-либо наблюдений за тремя переменными из модели, но это больше не сработало. Что было бы правильнее: исключить случай отсутствия наблюдений или оставить его в модели? И почему модель больше не работает, когда я ее исключаю? - person Alex; 31.01.2018
comment
Да, я понятия не имею, почему возникает эта ошибка. Это может быть проблема самого пакета. Комментарий github на самом деле не проясняет эту ссылку на Github coxphf. - person LAP; 31.01.2018
comment
Спасибо за все что ты сделал для меня. Если вы не возражаете, я спрошу, что бы вы порекомендовали сделать? - person Alex; 31.01.2018
comment
Есть ли у вас качественная информация о деле Undefined? Можно ли перекодировать его в одну из других категорий? Если нет, я бы запустил регрессию для всего набора данных и исключил переменную Undefined. Это не должно иметь значения. - person LAP; 31.01.2018
comment
Спасибо! В этом случае я буду запускать его целиком и исключать переменную Undefined. К сожалению, перекодировать не могу, так как классификация основана на генетических маркерах, а цитогенетического анализа у этого конкретного человека не было. - person Alex; 31.01.2018
comment
Сообщение об ошибке довольно сбивает с толку. Я не могу найти в ваших данных никакой структуры, которая сделала бы этот конкретный случай необходимым для любого вида вычислений. Как только вы что-либо измените, даже установив значение Favorable в 1 для этого случая, а затем оставив его, coxphf больше не будет работать. - person LAP; 31.01.2018
comment
Это действительно очень странно. Возможно, мне следует попытаться найти другой способ соответствовать модели Кокса. На всякий случай в самом пакете ошибка? У вас есть идея для другого метода? - person Alex; 31.01.2018
comment
Есть ли у вас дополнительные объясняющие переменные, которые вы могли бы включить в обычную модель Кокса, чтобы объяснить идеальную классификацию? Прямо сейчас вы прогнозируете только возраст и переменную генетического маркера, которая может быть немного тонкой. - person LAP; 31.01.2018
comment
Вся модель содержит некоторые молекулярные и генетические маркеры, а также пол и количество лейкоцитов. Однако, поскольку все они не вызывают проблемы, я удалил их, чтобы упростить пример в вопросе. К сожалению, когорта состоит только из 103 пациентов, поэтому я думаю, что идеальная классификация объясняется небольшой погрешностью выборки. - person Alex; 31.01.2018
comment
Вы можете задать другой вопрос на странице Cross Validated, чтобы получить больше информации о методе и вашей модели. - person LAP; 31.01.2018