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