Прогнозировать вероятность события на фактор с помощью модели Кокфа - PullRequest
Новая
       72

Прогнозировать вероятность события на фактор с помощью модели Кокфа

0 голосов
/ 28 января 2020

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

# Generate data
mydata <- data.frame(Site = as.factor(sample(c("SiteA", "SiteB", "SiteC"), 100, replace = TRUE)), 
                     Treatment = as.factor(sample(c("Treat.A", "Treat.B"), 100, replace = TRUE)), 
                     Time = sample(c(1, 2, 3), 100, replace = TRUE), 
                     Surv = sample(c(0, 1), 100, replace = TRUE)) # Alive is 0, death is 1


# Model
mymodel <- coxph(Surv(Time , Surv) ~ Treatment*Site, 
              data = mydata)

Мне нужна вероятность смерти через 3 года для каждого участка и каждого лечения (и соответствующего доверительного интервала). Можно ли извлечь эту информацию?

Исходя из различных форумов, на которых рассматривались подобные вопросы, я бы предположил добавить три столбца в мой набор данных с помощью команды: 

mydata$fit <- survfit(mymodel, newdata=mydata)$surv
mydata$lower <- survfit(mymodel, newdata=mydata)$lower
mydata$upper<- survfit(mymodel, newdata=mydata)$upper

И из этого оставьте только те строки, которые я Я заинтересован в этом. Однако это не работает, и команда генерирует вектор с в 3 раза большим количеством элементов, чем исходный набор данных (в данном примере 300 вместо 100). Есть ли что-то, что я неправильно понял?

Ответы [ 2 ]

0 голосов
/ 29 января 2020

Используйте predict.coxph со значением времени 

testset <-data.frame( Time=3, Surv=1,  # the Surv value is just a placeholder
                      Treatment=factor(rep(c("Treat.A", "Treat.B"),times=3)) , 
                      Site=factor(rep(c("SiteA", "SiteB", "SiteC"), each=2)))

testset$Surv3yr <- exp( -predict(mymodel, newdata=testset, typ="expected") )
testset
  Time Surv Treatment  Site   Surv3yr
1    3    1   Treat.A SiteA 0.1633725
2    3    1   Treat.B SiteA 0.3906895
3    3    1   Treat.A SiteB 0.3432062
4    3    1   Treat.B SiteB 0.2940677
5    3    1   Treat.A SiteC 0.5411742
6    3    1   Treat.B SiteC 0.2047518
0 голосов
/ 28 января 2020

Думаю, у вас возникла эта проблема, потому что элементы surv, lower и upper объекта, возвращаемого survfit, не являются векторами, они являются матрицами. Это дает вам выживание кривых , а не точечные прогнозы. Столбцы в этих матрицах связаны с указанными c комбинациями ковариат, появляющимися в строках фрейма данных, которые вы передали в survfit, тогда как строки этих матриц представляют полный диапазон (последовательных) временных шагов, наблюдаемых в вашем исходные данные. Если вы хотите, чтобы подгонянные значения для определенного c времени, t , вам нужно вытащить t -ую строку этой матрицы, то есть fitted$surv[t,].

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

library(survival)

# Generate data
set.seed(123)
mydata <- data.frame(Site = as.factor(sample(c("SiteA", "SiteB", "SiteC"), 100, replace = TRUE)), 
                     Treatment = as.factor(sample(c("Treat.A", "Treat.B"), 100, replace = TRUE)), 
                     Time = sample(seq(3), 100, replace = TRUE), 
                     Surv = sample(c(0, 1), 100, replace = TRUE)) # Alive is 0, death is 1


# Model
mymodel <- coxph(Surv(Time , Surv) ~ Treatment*Site, data = mydata)

# use expand.grid to get a table with all possible combinations of Site and Treatment
newdata <- with(mydata, expand.grid(Site = unique(Site), Treatment = unique(Treatment)))
# add a vector for your time of interest for clarity's sake; it won't actually factor into survfit
newdata$time = 3

# run survfit on that new table
fitted <- survfit(mymodel, newdata = newdata)

# extract the fitted values for the time slice of interest to you, here 3
newdata$fit <- fitted$surv[3,]
newdata$lower <- fitted$lower[3,]
newdata$upper <- fitted$upper[3,]

# result
print(newdata)
   Site Treatment time       fit      lower     upper
1 SiteA   Treat.B    3 0.3149307 0.15064889 0.6583612
2 SiteC   Treat.B    3 0.1721691 0.04597197 0.6447887
3 SiteB   Treat.B    3 0.3979556 0.18679672 0.8478130
4 SiteA   Treat.A    3 0.6117692 0.37752270 0.9913616
5 SiteC   Treat.A    3 0.3390650 0.15646255 0.7347769
6 SiteB   Treat.A    3 0.3128776 0.13297313 0.7361819
...