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

Все знают чем закончилась эта печальная история с трансатлантическим рейсом, который совершал "Титаник" в  апреле 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. Исходный текст