Я работаю с данными о выживании и использую survSplit(), чтобы справиться с непропорциональностью зависящих от времени коэффициентов. Модель основана на статье Терри Терно и др. и др. (2020) (https://cran.r-project.org/web/packages/survival/vignettes/timedep.pdf)
У меня есть факторная переменная с 6 уровнями для представления различных типов протезов коленного сустава. Когда я применяю survSplit() с любыми контрольными точками, коэффициенты для эталонного уровня этого скорректированного по времени коэффициента отображаются как NA в результатах. Коллинеарности нет, и проблема может быть воспроизведена с другими факторными переменными в данных. Кроме того, если я изменяю контрольное значение путем изменения уровней фактора, контрольное значение в любом случае равно NA.
Проблема воспроизводится ниже с переменной фактора celltype в данных veteran:
str(veteran)
'data.frame': 137 obs. of 8 variables:
$ trt : num 1 1 1 1 1 1 1 1 1 1 ...
$ celltype: Factor w/ 4 levels "squamous","smallcell",..: 1 1 1 1 1 1 1 1 1 1 ...
$ time : num 72 411 228 126 118 10 82 110 314 100 ...
$ status : num 1 1 1 1 1 1 1 1 1 0 ...
$ karno : num 60 70 60 60 70 20 40 80 50 70 ...
$ diagtime: num 7 5 3 9 11 5 10 29 18 6 ...
$ age : num 69 64 38 63 65 49 69 68 43 70 ...
$ prior : num 0 10 0 10 10 0 10 0 0 0 ...
library(tidyverse)
library(survival)
library(survminer)
df <- veteran
cox <- coxph(Surv(time, status) ~ celltype + age + prior, data = df)
cox.zph(cox, terms=F)
cox_tdc <- survSplit(Surv(time, status) ~ .,
data= df,
cut=c(150),
zero=0,
episode= "tgroup",
id="id") %>%
dplyr::select(id, tstart, time, status, tgroup, celltype, age, prior)
coxph(Surv(tstart, time, status) ~
celltype:strata(tgroup) + age + prior,
data=cox_tdc)
Call:
coxph(formula = Surv(tstart, time, status) ~ celltype:strata(tgroup) +
age + prior, data = cox_tdc)
coef exp(coef) se(coef) z p
age 0.005686 1.005702 0.009494 0.599 0.549262
prior 0.008592 1.008629 0.020661 0.416 0.677516
celltypesquamous:strata(tgroup)tgroup=1 0.300732 1.350848 0.360243 0.835 0.403828
celltypesmallcell:strata(tgroup)tgroup=1 1.172992 3.231649 0.325177 3.607 0.000309
celltypeadeno:strata(tgroup)tgroup=1 1.232753 3.430660 0.352423 3.498 0.000469
celltypelarge:strata(tgroup)tgroup=1 NA NA 0.000000 NA NA
celltypesquamous:strata(tgroup)tgroup=2 -1.160625 0.313290 0.450989 -2.574 0.010067
celltypesmallcell:strata(tgroup)tgroup=2 -0.238994 0.787420 0.542002 -0.441 0.659252
celltypeadeno:strata(tgroup)tgroup=2 1.455195 4.285319 0.837621 1.737 0.082335
celltypelarge:strata(tgroup)tgroup=2 NA NA 0.000000 NA NA
Likelihood ratio test=34.54 on 8 df, p=3.238e-05
n= 171, number of events= 128
Likelihood ratio test=69.43 on 9 df, p=1.97e-11
n= 272, number of events= 128
Проблема в том, что я не могу проверить Остатки Шенфельда как cox.zph() приводят к ошибке: «Ошибка в solve.default (imat, u): система вычислительно единственная: число взаимных условий = 5.09342e-19». Из-за АН.
Проблема с NA не возникает, если я не использую зависящие от времени коэффициенты (:strata(tgroup))
Кто-нибудь ранее имел дело с этой проблемой? Почему некоторые коэффициенты являются NA? Я очень ценю вашу помощь с этим!
РЕДАКТИРОВАТЬ: пример был изменен, чтобы включить воспроизводимые данные.
РЕДАКТИРОВАТЬ2: Исправлена точка отсечения времени в примере, который приводил к смещенным коэффициентам
РЕДАКТИРОВАТЬ3 : Я спросил Терри Терно об этой проблеме, и по электронной почте ниже:
Терри Терно: Здесь есть две проблемы.
Подпрограмма model.matrix используется lm, glm, coxph и множество других подпрограмм для создания X-матрицы для регрессии. Он пытается быть разумным, чтобы не создавать избыточные столбцы в X; эти столбцы будут иметь коэффициент NA. Это довольно хорошо, но не идеально. Ваш случай модели со стратами (tgroup): факторная переменная - это та, где она оставляет слишком много. Дополнительный NA в распечатке - это неудобство, но я не могу его исправить.
С другой стороны, подпрограмма cox.zph - моя проблема. Он должен игнорировать эти столбцы NA, и не делает. На самом деле есть код для проверки NA, но ваш пример показывает, что он неполон. Я добавлю кейс NA, как и ваш, чтобы мой набор тестов и исправил проблему. (Случай NA работал один раз, но какое-то обновление сломало его.)
Me: Когда я получаю результаты со строками NA, как в этом случае, коэффициенты все еще правильны, и NA представляет справочное значение?
Терри: Да, все коэффициенты верны. SAS, например, не пытается «предварительно исключить» столбцы, а в моделях с факторами всегда есть некоторые пропуски в векторе коэффициентов. Так как они используют "." вместо «NA» при печати пропуски не спрыгивают со страницы. Численно не существует наказания за то или иное действие.