Saturday, September 28, 2013

Выжить при крушении. Улучшаем прогноз.

В прошлый раз я использовал линейную регрессию для прогнозирования того, выжил ли пассажир при крушении или нет. Точность прогноза составила 0.77, что неплохо для начала. Самыми значимыми переменными (по значению P-Value) были Sex, Fare, Age. Их я и брал в для прогнозирования.

Смотрим чего не хватает

Что делать, если для некоторых значимых переменных отсуствуют значения? Например, для поля возраст в тренировочном наборе отсутствует значение в 177 случаях:
Один из участников соревнования предложил заполнить эти значения спрогнозированными данными. Детальнее можно ознакомиться на форуме Kaggle, я же опишу как я применил его идею.

Добавляем новую колонку

От каких переменных может зависеть возраст? Возможно, что от количества детей (Parch) или братьев-сестер (Sibsp), а может стоимости билета (Fare)? Может можно из данных вытащить дополнительную информацию, которая нам помогла бы спрогнозировать возраст? Например, имя пассажира уже содержит его титул: Mr, Mrs, Miss, Master.

Написать скрипт, который из полного имени вытащит титул из имени не составит особого труда. Я написал его на питоне, точнее на его расширенной версии Anaconda, которая содержит массу нужных пакетов, в том числе модуль для чтения/записи csv файлов. Скрипт находится здесь. После обработки в файлах с тренировочными и тестовыми данными появилась колонка Title.
Считываю полученный csv файл в R:

Определяем переменные с помощью corrgram

Дальше определим, от каких переменных зависит возраст. Для этого воспользуюсь таким графиком, который называется коррелограмма (correlogram). Его я тоже нарыл на том же форуме, ссылка внизу. Реализована эта красота в пакете corrgram.
Поскольку corrgram показывает зависимости между числовыми переменными, добавлю колонку TitleNum, в которой будет числовое значение от Title:
Строю график.
Вот что получилось:
Чем насыщеннее цвет в нижней части, тем более тесная связь между переменными. Синий цвет - связь положительная, красный - отрицательная. Интересные графики в верхней части показывают эти изменения графически (чтобы отобразились графики в верхней части нужно задать параметер upper.panel=panel.ellipse ).
Например, между Pclass и Age есть сильная отрицательная связь, т.к. соответствующий цвет квадратика в нижней части темно-красный. На графике так же видно, что с уменьшением Pclass возраст пассажира увеличивается (сейчас, как и 100 лет назад первым классом в основном путешествовали старперы, студики и прочий нищеброд ютились в эконом классах).

Заполняем недостающие данные

Переменные, наиболее связанные с возрастом определили - это Pclass, SibSp и Title. Осталось заполнить отстутствующие значения и запихнуть в общую модель. Недостающие значения в тренировочных и тестовых данных я проставляю в процедуре, которую вызываю перед тем, как построить модель и делать прогноз. Вот код функции:
```{r}
predictMissing <- function(dataSet){
  #Find missing Age
  age.model = lm(Age~as.factor(Title)  + SibSp + Pclass, data=dataSet)
  fare.model = lm(Fare~Parch+Pclass, data=dataSet)
  for (i in 1:nrow(dataSet)){
    if (is.na(dataSet[i, "Age"])){
      dataSet[i,"Age"] = predict(age.model, newdata=dataSet[i,])
    }
  }
  return(dataSet)
}
```{r}

Финальный результат

После отправки новых прогнозов на Kaggle я увеличил счет на один процент(с 0.77512 до 0.78469). Немного, но главное - это полученный опыт! :)

Ресурсы

1. Пост от Wallace Campbell. https://www.kaggle.com/c/titanic-gettingStarted/forums/t/5232/r-code-to-score-0-79426
2. Anaconda, python visualization and explaration tool. https://store.continuum.io/cshop/anaconda/
3. Скрипт добавления колонки title. https://dl.dropboxusercontent.com/u/22607711/AddTitle.py
4. Пост на Kaggle форуме про корелограммы. https://www.kaggle.com/c/titanic-gettingStarted/prospector#489

Monday, September 02, 2013

Выжить при крушении. Анализ данных катастрофы "Титаник".

Все знают чем закончилась эта печальная история с трансатлантическим рейсом, который совершал "Титаник" в  апреле 1912 года. В результате крушения погибло 1495 человек, спаслось 712. Эта катастрофа стала предметом многих исследований, документальных и художественных фильмов. Создатели онлайн площадки для научного моделирования Kaggle тоже решили использовать эту историю вместе с данными о всех пассажирах Титаника для тренировочного задания.

Коротко о задании

Есть данные о 891м пассажире (в формате csv), которые содержат класс, которым путешествовал, имя, пол, возраст, цену билета и т.д. Также дана информация о том выжил ли пассажир. Это т.н. тренировочный набор данных. Вместе с тренировочными данными прилагается тестовый набор, который содержит такую же информацию о еще 418 пассажирах кроме индикатора "выжил/нет".  Необходимо предсказать кто из 418 выжил.
Более подробно о деталях можно посмотреть здесь https://www.kaggle.com/c/titanic-gettingStarted

Инструменты

Для построения модели я буду использовать R и RStudio. Ссылки для скачивания внизу. Мне удобнее работать в RStudio, чем в консоли R т.к. там наряду с возможностью выполнять команды можно создать исполняемый документ (Rmd - R Markup Document или "чистый" R). RMD - это Markup document со скриптовыми вставками на языке R, в результате выполнения которого получается HTML.

Загружаем данные

Загрузил тренировочный и тестовый наборы данных. Полный тренировочный набор разделил на данные для кросс-валидации (30% записей) и тренировочный (trainFaith). Это делается для того,чтобы затем можно было оценить эффективность модели. На тренировочных данных (trainFaith) буду "обучать" модель, а данные для кросс-валидации (cvFaith) использую для оценки.
Делал примерно так:

```{r}
set.seed(333)
l=length(train[,1]);
cvSamples <- sample(1:l,size=(l/3),replace=F)
cvFaith <- train[cvSamples,]
trainFaith <- train[-cvSamples,]
```

Определяем переменные

От чего зависит выжил пассажир или нет? Возможно, что сначала спасали женщин и детей, а может пассажиров первого класса спасали первыми? По одноименному фильму вполне можно сделать такие предположения. Тот, кто смотрел фильм, уже может считаться "экспертом" в предметной области. А что делать человеку, который ничего не слышал о Титанике? Можно использовать научный подход к определению первоначальных переменных для дальнейшего построения модели.
Для "научного" подхода я буду использовать P-Values. Очень хорошее объяснение что такое P-Values можно посмотреть на Khan Academy, набор лекций Hypothesis Testing and P-values.
Если коротко, то P-Value - это величина, которая показывает насколько принятая гипотеза не влияет на данные. Например, гипотеза: у женщин больше шансов быть спасенными при крушении. Если эта гипотеза неверна, то P-value для поля Sex по отношению к Survived будет близко к 1. К сожалению, не всегда можно сказать, что если P-value близко к 0, то гипотеза верна, т.к. на данные могут влиять другие факторы.
Самый простой способ получить эти значения в R - это "впихнуть" (по-англ. fit) данные в линейную модель и вызвать на ней функцию summary. Самая простая модель линейной регрессии, ее и возьмем:

```{r}
glm1=glm(Survived ~ Sex + Fare + Age + SibSp + Parch, data=train, family="binomial")
summary(glm1)
```
Здесь family="binomial" указана для построения логистической регрессии.
Вот что получилось. Смотрим сразу на раздел coefficients, колонка Pr:
Для построения модели возьму переменные Sex, Fare, Age и SibSp.

Обучение модели

Тут все просто. Впихиваю тренировочные данные в модель. В качестве набора - trainFaith

```{r}
predictFunction = Survived ~ Sex + Fare + Age + SibSp
glm1 = glm(predictFunction, data=trainFaith, family="binomial")
```

Прогнозирование
Используем обученную модель для предсказаний кто выжил. Встроенная функция predict делает практически все что нам надо. Параметер type="response" конвертирует значения из логарифмической шкалы в обычную. В качестве данных используем cvFaith (помним, что это 30% случайно отобранных данных из полного тренировочн набора).
Далее я создал новый набор, который состоит из 2х колонок: предсказанное значение и PassengerId.

```{r}
rawpredict=predict(glm1, newdata=cvFaith, type="response")
cvPredicted=data.frame(Survived = rawpredict, PassengerId = cvFaith$PassengerId)
head(cvPredicted)
```


Вероятность выживания в колонке Survived:

Оценка

Как известно, что выжить частично очень сложно, поэтому предсказанные значения, которые находятся в диапазоне от 0 до 1 нужно конвертнуть в 0 (не выжил) или 1 (выжил). Тут я решил выпендриться и использовать стороннюю библиотеку. Я хочу посмотреть изменится ли точность моих предугадываний, если в качестве порога я возьму число, отличное от 0.5.

```{r}
library(ROCR)
pred = prediction(cvPredicted$Survived, cvFaith$Survived)
perf <- performance(pred,"acc")
plot(perf)
```
Получился такой график:
По оси Y точность прогнозирования (accuracy), по оси X - порог. Т.е. если в качестве порога я возьму 0.8, то точность модели будет около 0.65. Выбрать 0.5 в качестве порога вполне резонно, т.к. это обеспечит максимальную точность.

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

```{r}
threshold = function (val) as.numeric(val > 0.5)
glmProd = glm(predictFunction, data=train, family="binomial")
rawpredict=predict(glmProd, newdata=test, type="response")
predicted=data.frame(Survived = threshold(rawpredict), PassengerId = test$PassengerId)
```

Запись в файл производится одной командой:
```{r}
write.table(predicted, "D:/Dropbox/kaggle/Titanic/output-R.csv", row.names=FALSE, quote=FALSE, sep=",", append=FALSE)
```

После отправки на kaggle эта модель дала точность 0.77033. Это немногим больше, если прогнозировать только по полю Sex (это дало бы 0.7655). Но для начала пойдет. 


Ресурсы

1. Язык для статистической обработки данных R
2. RStudio - IDE для R
3. P-Values https://www.khanacademy.org/math/probability/statistics-inferential/hypothesis-testing/v/hypothesis-testing-and-p-values
4. Здесь много алгоритмов моделирования : http://www.machinelearning.ru
5. Пакет ROCR http://rocr.bioinf.mpi-sb.mpg.de/
6. Исходный текст