Получение P-значений нуля в регрессии Кокса: R

Question

Я студент, проводящий анализ выживания по экспрессии генов в R. У меня есть данные по экспрессии для 249 пациентов, и я использую 6000 генов, а также время их выживания без событий и жизненное состояние в качестве переменных ответа. Когда я попытался запустить регрессию Кокса на моем наборе данных, я получил чрезвычайно странные результаты (значения р 0,00 и странная опасность крысы ios). Я проверял свой код несколько раз, но не смог уловить свою ошибку (когда я пытался ранее с одним геном, он работал нормально, но когда я пытаюсь протестировать несколько генов с помощью функции '.', Я не получение результатов порпера). Я был бы очень признателен за любую помощь и приложил мой код и вывод! Дайте мне знать, если потребуется дополнительная информация.

library(survival)
options(expressions = 5e5)
firstSplitData <- read.delim("/Users/menon/OneDrive/Desktop/csrsef files/FirstSplitDataFrame.txt")
firstInitialData <- data.frame(firstSplitData)
firstEventFreeTime <- firstInitialData[ , c("EFST")] 
firstVitalStatus <- firstInitialData[, c("Status")]
#create a temporary object to use in the final object in order to be able to use '.'
temporaryObj <- Surv(as.numeric(firstEventFreeTime), firstVitalStatus == 2)
firstFinalData <- data.frame(SurvObj = temporaryObj)
#bind the two together for the final data 
firstFinalData <- cbind(firstFinalData, firstInitialData[, 2:ncol(firstInitialData)])
#create final cox model
firstCox <- coxph(SurvObj ~ ., data =  firstFinalData)
summary(firstCox)$coefficients

А вот (некоторые из) мой вывод:

> summary(firstCox)$coefficients
                     coef     exp(coef)     se(coef)             z      Pr(>|z|)
EFST         3.644083e-03  1.003651e+00 0.0001340611    27.1822581 1.052851e-162
Status      -2.926090e+00  5.360625e-02 0.3182658189    -9.1938542  3.790122e-20
AADACL3      1.502153e+02  1.728460e+65 0.3665374081   409.8224582  0.000000e+00
AADACL4      5.857192e+01  2.738174e+25 0.3681708023   159.0889828  0.000000e+00
ACADM        2.455978e+02 4.589695e+106 0.2175220391  1129.0710334  0.000000e+00
ACAP3        4.093913e+02 6.256964e+177 0.2756635268  1485.1121632  0.000000e+00
ACOT11       1.940976e+01  2.688751e+08 0.3251033140    59.7033512  0.000000e+00
ACOT7       -2.841794e+02 3.823403e-124 0.3139848504  -905.0736377  0.000000e+00
ACTB        -5.562202e+01  6.976896e-25 0.3173481100  -175.2713234  0.000000e+00
ACTL8       -4.017414e+02 3.356676e-175 0.3435128215 -1169.5093020  0.000000e+00
ACTRT2      -7.613568e+01  8.603881e-34 0.2861088372  -266.1074036  0.000000e+00
ADC         -1.244476e+02  8.976070e-55 0.3201452217  -388.7223972  0.000000e+00
ADPRHL2      4.887427e+01  1.681998e+21 0.2895110526   168.8165913  0.000000e+00
AGMAT        7.266946e+02           Inf 0.4295874196  1691.6104194  0.000000e+00
AGO1         3.352041e+02 3.778188e+145 0.2633158947  1273.0111995  0.000000e+00
...

И вот что dput(firstFinalData[1:10, 1:10]) производит:

structure(list(SurvObj = structure(c(444, 5553, 5296, 922, 205, 
47, 401, 245, 263, 5564, 1, 0, 0, 1, 1, 1, 1, 1, 1, 0), .Dim = c(10L, 
2L), .Dimnames = list(NULL, c("time", "status")), type = "right", class = "Surv"), 
    EFST = c(444L, 5553L, 5296L, 922L, 205L, 47L, 401L, 245L, 
    263L, 5564L), Status = c(2L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 
    2L, 1L), AADACL3 = c(5.52132, 5.64712, 5.45876, 5.71481, 
    5.1269, 5.88764, 5.08912, 4.91729, 5.65387, 5.59824), AADACL4 = c(5.17251, 
    5.41843, 5.10969, 5.23402, 4.60353, 5.70923, 5.02245, 5.1466, 
    4.8355, 4.83986), ACADM = c(7.47834, 7.43494, 7.91155, 7.86337, 
    8.39009, 6.16251, 7.83793, 7.71742, 6.98061, 7.78087), ACAP3 = c(7.80589, 
    8.00354, 7.75014, 7.61566, 7.55267, 7.9449, 7.20561, 7.99776, 
    7.72778, 7.43355), ACOT11 = c(6.75915, 6.30386, 6.38214, 
    6.54392, 6.64743, 6.78981, 6.42641, 6.58761, 6.66693, 6.53731
    ), ACOT7 = c(8.11807, 8.38011, 7.8349, 8.43645, 8.11502, 
    8.0109, 7.6866, 8.55327, 8.17004, 7.44455), ACTB = c(10.8227, 
    11.4556, 11.4216, 11.332, 10.9536, 9.83797, 11.2352, 11.5006, 
    11.1817, 10.895)), row.names = c(NA, 10L), class = "data.frame")

Большое спасибо!

Редактировать:

Я также получаю это предупреждение при запуске firstCox <- coxph(SurvObj ~ ., data = firstFinalData):

In fitter(X, Y, istrat, offset, init, control, weights = weights,  :
  Ran out of iterations and did not converge

Answer 1

Не включайте объект Surv () во фрейм данных.

firstFinalData <- firstFinalData[,-1]
firstCox <- coxph(Surv(EFST, Status) ~ ., data =  firstFinalData)

Он должен работать ( edit : для меньшего числа переменных).

Как говорит Мориц Эверс, использование модели с 6000 предикторов (генов) только для 249 субъектов приведет к проблемам сходимости. Рассмотрите возможность уменьшения количества генов (или получения большего количества пациентов!)

Answer 2

Если вы хотите выполнить несколько моделей регрессии Кокса с одним предиктором, вы можете использовать следующий код, используя ваши опубликованные примеры данных. Сначала я удаляю объект выживания в первом столбце.

myData <- finalData[,-1]

library(survival)
firstCox <- co

coxph(Surv(EFST, Status) ~ ., data =  myData) 

Это возвращает указанное предупреждение (слишком много предикторов)

Warning message:
In fitter(X, Y, istrat, offset, init, control, weights = weights,  :
  Ran out of iterations and did not converge

Чтобы запустить несколько одномерных моделей, сначала создайте список одномерных формул:

formulas <- sapply(names(myData)[3:9], function(x) as.formula(paste('Surv(EFST, Status) ~ ',x)))

Создайте список моделей, используя функцию coxph:

models <- lapply(formulas, function(x) coxph(x, data=myData))

Извлеките крысу опасности ios (exp(coef)) и с вероятностью 95% интервалы:

res <- lapply(models, function(x) return(cbind(HR=exp(coef(x)), exp(confint(x)), Pval=coef(summary(x))[5])))
res

---

$AADACL3
               HR       2.5 %   97.5 %     Pval
AADACL3 0.1858129 0.008579879 4.024119 0.283442

$AADACL4
               HR      2.5 %   97.5 %      Pval
AADACL4 0.8481017 0.02748128 26.17333 0.9249839
...

Answer 3

За исключением первых двух коэффициентов (EFST и Status), коэффициенты для всех других генов либо очень малы, либо чрезвычайно велики, что приводит к очень большой отрицательной / положительной t -статистике, который объясняет p-значения, которые вы видите.

Я не совсем уверен, что понял, что ты делаешь. Разве регрессия на 6000 генов в ваших данных о 249 пациентах не означает, что у вас гораздо больше параметров, чем наблюдений?

В этом случае вы сталкиваетесь с проблемами переоснащения, которые объясняют оценки параметров.