Share |

среда, 23 августа 2017 г.

Прогноз индивидуального дожития и тюнинг параметров регрессии Кокса в пакете mlr



Одна из проблем в анализе дожития заключалась в том, что не было четких алгоритмов тюнинга параметров модели. В разных пакетах типа CoxBoost есть варианты кросс валидации, определения числа итераций и т.п.. Но grid_search не было ни в одном из пакетов.
Теперь такая возможность существует в пакете mlr. Это пакет фреймворк типа caret, который является оберткой для алгоритмов машинного обучения. 
Я покажу, как тренируются параметры регресии Кокса в этом пакете. Регрессия Кокса важная HR в процессах управления текучестью персонала, карьерного роста, других процессах, где мы управляем ростом. 
В конкретном примере я опускаю некоторые детали тип препроцессинга и т.п., что сосредоточиться исключительно на тюнинге параметров и прогнозе кривой дожития.
А вы не забывайте кликать на дирет рекламу в качестве благодарности. Также добавлю, что это один из сюжетов Семинар-практикум "HR-Аналитика в R", Москва,16-17 ноября 2017 - я вас там буду этому учить. 

Датасет

Я не стал возиться со сложными датасетами, взял данные, предоставленные мне Максимом Андреенок. Данные включают в себя:
  1. стаж работы в компании, 
  2. событие - уволился на момент анализа или нет, 
  3. наличие наставника, 
  4. пол работник,
  5. пол руководителя.
Всего 427 строк данных. 

Код

поехали.Загружаем пакеты.
library(mlr)
library(survival)
library(pec)
Пропущу загрузку данных, сразу разбиваю на трейн и тест сет.

train = sample(nrow(data), 0.7 * nrow(data))
test = setdiff(seq_len(nrow(data)), train)
train.task = makeSurvTask(data = data[train, ], target = c("stag", "event"))
train.task
test.task = makeSurvTask(data = data[test, ], target = c("stag", "event"))
test.task
Задаем учителя

lrn = makeLearner("surv.coxph", id = "cph")
Перед тюнингом параметров давайте посмотрим, какие параметры вообще можно тренировать в регрессии Кокса.

getParamSet("surv.coxph")
                Type len      Def              Constr Req Tunable Trafo
ties        discrete   -    efron efron,breslow,exact   -    TRUE     -
singular.ok  logical   -     TRUE                   -   -    TRUE     -
eps          numeric   -    1e-09            0 to Inf   -    TRUE     -
toler.chol   numeric   - 1.82e-12            0 to Inf   -    TRUE     -
iter.max     integer   -       20            1 to Inf   -    TRUE     -
toler.inf    numeric   - 1.35e-06            0 to Inf   -    TRUE     -
outer.max    integer   -       10            1 to Inf   -    TRUE     -
model        logical   -    FALSE                   -   -   FALSE     -
x            logical   -    FALSE                   -   -   FALSE     -
y            logical   -     TRUE                   -   -   FALSE     -
Надеюсь, читатели сами разберутся, какие параметры будут тренировать они. Я беру вот такие параметры

surv_param = makeParamSet(
  makeDiscreteParam("ties",  values = c('efron','breslow','exact')),
  makeIntegerParam("iter.max", lower = 1, upper = 150),
  makeIntegerParam("outer.max", lower = 1, upper = 50)
)
rancontrol = makeTuneControlRandom(maxit = 10L)
set_cv = makeResampleDesc("RepCV", folds = 5L, reps = 5L)
Сразу туда же кросс валидацию на пяти фолдах. И самое любимое

surv_tune = tuneParams(learner = lrn, resampling = set_cv, task = train.task, 
                        par.set = surv_param, control = rancontrol)
Не указываю меру качества модели, по умолчанию cindex - конкорданс.
$ties
[1] "breslow"

$iter.max
[1] 139

$outer.max
[1] 15
Заметьте, по дефолту в регрессии Кокса стоит метод efron, а машинка нам показывает, что breslow получшее будет.
surv_tune$y
cindex.test.mean 
        0.752014 
0.752014 - это про наше качество модели. Рекомендую Вам отдельно прочитать, что обозначает параметр C-index, скажу только, что он аналогичен площади под кривой и находится в пределах от 0, 5 (совсем плохо) и 1 (нереально хорошо). В нашем случае даже 0, 75 кажется нереально хорошо, потому что всего три переменные. Делаем модель
surv.tree = setHyperPars(lrn, par.vals = surv_tune$x)
surva = mlr::train(surv.tree, test.task)
getLearnerModel(surva)
model = predict(surva, test.task)
model
И смотрим, что нам говорят тестовые данные

rcorr.cens(-model$data$response, 
           Surv(data[test, ]$stag, data[test, ]$event))["C Index"]
 C Index 
0.6898021 
Можно еще такой метод проверки сделать - на пакете risksetROC и тестовых данных
w.ROC = risksetROC(Stime = data[test, ]$stag,  
                   status = data[test, ]$event, 
                   marker = model$data$response, 
                   predict.time = median(data[test, ]$stag), 
                   method = "Cox", 
                   main = paste("OOB Survival ROC Curve at t=", 
                                median(model$data$truth.time)), 
                   lwd = 3, 
                   col = "red" )

w.ROC$AUC
[1] 0.6648762
Хуже, но тоже неплохо. w.ROC вам нарисует диаграмму под кривой, я не буду ее сюда помещать.

Согласитесь, неплохо для тестовых данных - 0, 69. А почему не 0, 75, как на трейне?

getLearnerModel(surva)
Call:
survival::coxph(formula = f, data = data, ties = "breslow", iter.max = 139L, 
    outer.max = 15L)

             coef exp(coef) se(coef)     z      p
X         0.00234   1.00235  0.00133  1.76 0.0778
coach    -0.95264   0.38572  0.33482 -2.85 0.0044
matchж.м -0.42199   0.65574  0.77277 -0.55 0.5850
matchм.ж -0.03979   0.96099  0.40376 -0.10 0.9215
matchм.м  0.33148   1.39303  0.41926  0.79 0.4292
А потому что у нас значим только один фактор - наличие наставника, остальные факторы - шум. Надо их отбросить, чтобы получить нормальную модель.
И у нас осталась, еще одна задача - самая главная наша задача - прогноз индивидуального риска или дожития. Иначе на практике наши упражнения будут теорией. И вот тут фишка заключается в том, что пакет mlr на сегодня не дает возможности предсказывать индивидуальные риски и дожития (информация в github авторов пакета). На этот счет есть разные подходы, но все они за пределами mlr. Вот как например дают прогноз через пакет pec
# train.coxph:
mod = coxph(Surv(time, status) ~ ., data = data[train, ])

# predict.coxph:
probs = predictSurvProb(mod, newdata = data[test, ], times = data[test, "time"])
timepoints = seq(0, max(data$time), length.out = 100L)
probs = predictSurvProb(mod, newdata = data[test, ], times = timepoints)

pec(probs, Surv(time, status) ~ 1, data = data[test,], exact = F,
    exactness = 99L, maxtime = max(timepoints))
probs
     [,1]      [,2]      [,3]      [,4]      [,5]      [,6]      [,7]      [,8]       [,9]      [,10]
3      1 0.9989732 0.9979358 0.9979358 0.9968463 0.9968463 0.9956966 0.9933632 0.98563169 0.98563169
9      1 0.9989948 0.9979792 0.9979792 0.9969126 0.9969126 0.9957871 0.9935025 0.98593217 0.98593217

Я даю в обобщенном виде, понятно, что в mod вы должны включать полученные параметры (покажу ниже), timepoints я выбрал произвольно, вам рекомендую выбирать по количеству единиц измерения вашего времени. В нашем случае это месяцы, их и можно выбрать.
Я же делаю так

w = coxph(Surv(stag, event) ~ . , data = data[train, ], ties=c("breslow"), iter.max = 139, outer.max = 15)
Я подставляю полученные параметры в регрессию Кокса (ну как будто мы убрали не значимые параметры и заново натренировали).
И тут есть одна важная тонкость: в этой модели w, коэффициенты другие, чем в тренированной.
w$coefficients
          X       coach    matchж.м    matchм.ж    matchм.м 
 0.00382257 -1.74394635  0.58288003 -0.31481756  0.04525168 

Наличие наставника в этой модели идет с коэффициентом -1, 7439, а в той, что тренировали - -0.95264. Не очень приятная разница. Думаю, согласитесь, что второй коэффициент, прошедший сетку и кросс валидацию, является более точным, поэтому нам нужно использовать его, вместо -1, 7439. Мы это исправляем так
w$coefficients[2] = -0.95264
Проверяем
X       coach    matchж.м    matchм.ж    matchм.м 
 0.00382257 -0.95264000  0.58288003 -0.31481756  0.04525168 

Все зашибись. Далее из тест сета я выбираю произвольно троих чуваков и смотрю, как они будут себя вести

newdude = data[test, ][c(4,5,12), ]
e = survfit(w, newdata=newdude)
quantile(e)
Прогноз индивидуального дожития и тюнинг параметров регрессии Кокса в пакете mlr

По оси X у нас стаж работы в компании, по оси Y - вероятность дожития. Обратите внимание особо, что крайне нижнее значение по оси Y равно 0, 4. И совершенно очевидно, что у Петрова наставника не было (ну или в нашем случае мы прогнозируем, что будет с Петровым, если ему не назначить наставника), а у Иванова с Сидоровым они есть. Но даже для Петрова медиана дожития, или если пользоваться HR термином - средний срок жизни в компании - в районе пяти лет. Ну что поделать, такие суровые белорусские реалии.
Это все, что я хотел показать. Понятно, что владея этим инструментом, вы сможете строить более продвинутые модели. 

Понравился пост? 

и Вы захотите выразить мне благодарность, просто покликайте на директ рекламу ниже на странице - у вас это отнимет несколько секунд, а мне принесет немного денег.
На этом все, читайте нас в фейсбуке 

Комментариев нет:

Отправить комментарий

Популярные сообщения