HW11 // Логистическая регрессия. Пуассоновская регрессия

Данные

Сегодня мы работаем с данными о факторам сердечно-сосудистых заболеваний. Ниже представлено описание данных из источника.

Introduction World Health Organization has estimated 12 million deaths occur worldwide, every year due to Heart diseases. Half the deaths in the United States and other developed countries are due to cardio vascular diseases. The early prognosis of cardiovascular diseases can aid in making decisions on lifestyle changes in high risk patients and in turn reduce the complications. This research intends to pinpoint the most relevant/risk factors of heart disease as well as predict the overall risk using logistic regression Data Preparation

Source The dataset is publically available on the Kaggle website, and it is from an ongoing cardiovascular study on residents of the town of Framingham, Massachusetts. The classification goal is to predict whether the patient has 10-year risk of future coronary heart disease (CHD).The dataset provides the patients’ information. It includes over 4,000 records and 15 attributes.

Variables Each attribute is a potential risk factor. There are both demographic, behavioral and medical risk factors.

  • Demographic:
    • male: male or female (Nominal)
    • age: Age of the patient; (Continuous - Although the recorded ages have been truncated to whole numbers, the concept of age is continuous)
  • Behavioral:
    • currentSmoker: whether or not the patient is a current smoker (Nominal)
    • cigsPerDay: the number of cigarettes that the person smoked on average in one day (can be considered continuous as one can have any number of cigarettes, even half a cigarette.)
  • Medical (history):
    • BPMeds: whether or not the patient was on blood pressure medication (Nominal)
    • prevalentStroke: whether or not the patient had previously had a stroke (Nominal)
    • prevalentHyp: whether or not the patient was hypertensive (Nominal)
    • diabetes: whether or not the patient had diabetes (Nominal)
  • Medical (current):
    • totChol: total cholesterol level (Continuous)
    • sysBP: systolic blood pressure (Continuous)
    • diaBP: diastolic blood pressure (Continuous)
    • BMI: Body Mass Index (Continuous)
    • heartRate: heart rate (Continuous — In medical research, variables such as heart rate though in fact discrete, yet are considered continuous because of large number of possible values.)
    • glucose: glucose level (Continuous)
  • Predict variable (desired target):
    • TenYearCHD: 10 year risk of coronary heart disease CHD (binary: “1”, means “Yes”, “0” means “No”)

Необходимо выяснить, какие факторы влияют на риск сердечнососудистых заболеваний в перспективе 10 лет. Переменную TenYearCHD мы будем рассматривать как целевую, все остальные — в качестве предикторов.

Основные задания

#1

  1. Загрузите датасет в переменую heart. Изучите их структуру.
  2. Предобработайте данные, если это необходимо: удалите пропущенные значения, скорректируйте типы переменных. Перезапишите датасет.

В качестве ответа для самопроверки введите количество строк в датасете, полученном после удаления пропущенных значений.

Подсказка

Удалить все пропущенные значения можно с помощью функции drop_na().

#2

Сегодня мы поработаем с изученной моделью в более «машиннообученском». подходе. Начнем в разбиения выборки на обучающую (тренировочную) и тестовую.

Разбейте данные на обучающую и тестовую выборки в соотношении 7/3 — 70% данных должно оказать в обучающей и 30% в тестовой. Сохраните обучающий датасет в переменную heart_train, а тестовый — в heart_test.

Для воспроизводимых результатов используйте set.seed(616).

В качестве ответа для самопроверки введите среднее по переменной TenYearCHD обучающей выборки. В качестве десятичного разделителя используйте точку.

Подсказки

Чтобы разбить датасет на тренировочную и тестовую части, необходимо:

  • сгенерировать вектор случайных индексов из номеров строк общего датасета длиной, равной 70% от числа строк датасета
  • с помощью индексации (или функции из tidyverse) вытащить строки со сгенерированными номерами и сохранить результат в heart_train
  • с помощью «отрицательной» индексации (или функции из tidyverse) вытащить строки со всеми номерами, кроме сгенерированных, и сохранить результат в heart_test

#3

  1. Постройте модель на обучающей выборке. Целевой переменной выступит TenYearCHD, все остальные — в качестве предикторов. Так как в модели достаточно много предикторов, взаимодействия мы включать не будем.
  2. Протестируйте общую статистическую значимость модели с помощью теста отношения правдоподобий. Проинтерпретируйте результаты.

В качестве ответа для самопроверки введите значение девиансы построенной модели (deviance), округленное до целого.

Подсказки

Чтобы проестировать гипотезу об общей статистической значимости модели, необходимо:

  • построить нулевую модель, аналогичную заданной, но только с одним предиктором — интерсептом
  • провести тест отношения правдоподобий с помощью функции anova(), указав в аргументе test необходимый тест
  • значение аргумента test соответствует распределению, которому подчиняется распределение значения отношения правдоподобий

#4

  1. Используя тесты отношения превдоподобий, упростите модель, исключив незначимые предикторы.
  2. Сравните с помощью статистического теста и информационных критериев исходную и упрощенную модели.
  3. Выберите, с какой в дальнейшем вы будете работать:
    • если модели различаются, то мы, исключив предикторы, потеряли что-то важное, значит необходимо работать с полной моделью
    • если модели не различаются, то мы, исключив предикторы, не потеряли ничего важного, значит проще работать с сокращенной моделью

В качестве ответа для самопроверки введите значение девиансы сокращенной модели (deviance), округленное до целого.

#5

  1. Постройте confusion matrix по тренировочным данным для выбранной в предыдущем задании модели.
    • Так как мы работаем с медицинскими данными, воспользуемся достаточно либеральным порогом для получения категориальных предсказаний в 0.6. Идея исходит из подхода гипердиагностики — лучше взять и обследовать пациента с меньшим риском, чтобы выявить потенциально опасные патологии на раннем этапе.
  2. Рассчитайте метрики предказательной силы модели — accuracy, precision, recall, F1-меру — по построенной confusion matrix.
  3. Дайте оценку качеству модели по полученным значениям.

Для самопроверки в поле ниже введите значение F1-меры, округленное до тысячных. В качестве десятичного разделителя используйте точку.

Подсказки

Confusion matrix имеет следующую структуру:

Предсказания: \(0\) Предсказания: \(1\)
Данные: \(0\) \(\text{TN}\) \(\text{FP}\)
Данные: \(1\) \(\text{FN}\) \(\text{TP}\)
  • True Positive (\(\text{TP}\)) — верное предсказанные единицы
  • True Negative (\(\text{TN}\)) — верно предсказанные нули
  • False Positive (\(\text{FP}\)) — ложноположительные предсказания, ошибочно предсказанные единицы
  • False Negative ($ — ложноотрицательные предсказания, ошибочно предсказанные нули

Формулы для расчета метрик таковы:

\[ \begin{split} & \text{accuracy} = \frac{\text{TP} + \text{TN}}{\text{TP} + \text{TN} + \text{FP} + \text{FN}} \\ & \text{precision} = \frac{\text{TP}}{\text{TP} + \text{FP}} \\ & \text{recall} = \frac{\text{TP}}{\text{TP} + \text{FN}} \\ & \text{F1} = 2 \cdot \frac{\text{precision} \cdot \text{recall}}{\text{precision} + \text{recall}} \end{split} \]

#6

  1. Возьмите несколько более строгих значений порогов — 0.7, 0.8 и 0.9 — и посмотрите, справляется ли модель с предсказаниями лучше, если использовать их.
  2. Выберите, при каком пороге модель наиболее хорошо описывает данные. Сравнивая модели, неободимо посмотреть на все метрики, однако для принятия окончательного решения можно воспользоваться F1-мерой, как обобщенным показателем precision и recall.

Возможно, проще будет сначала сделать задание 10. Это сократит код.

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

#7

С помощью имеющейся модели и выбранного порога получите предсказания для тестовых данных.

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

Подсказки
  • Чтобы получить предсказания для тестовых данных, воспользуйтесь функцией predict(), в которую необходимо передать модель и новые данные, а также указать в аргументе type значение "response", чтобы получить предсказания вероятностей.
  • Далее, используя выбранный порог, получите дискретные предсказания модели.

#8

  1. Постройте confusion matrix для предсказаний на тестовой выборке и рассчитайте метрики предсказательной силы.
  2. Сопоставьте результаты с метриками, полученными на тренировочных данных.

Для самопроверки в поле ниже введите значение F1-меры, округленное до тысячных. В качестве десятичного разделителя используйте точку.

Подсказки
  • При сопоставлении метрик необходимо понять, выше или ниже они получились на тестовой выборке по сравнению с обучающей (тренировочной).
    • Если они ниже, то проблем с переобучением модели не возникло — можно пытаться экстраполировать результаты на генеральную совокупность.
    • Если они выше, то результаты моделирования зависят больше от данных, чем от самой модели — есть трудности с обобщением результатов на генеральную сопокупность.

#9

Проверьте основные допущения рассматриваемой модели — линейность связи между предикторами и целевой переменной и отсутствие сверхдисперсии. Проинтерпретируйте полученные результаты.

Для проверки линейности достаточно графической диагностики, для проверси сверхдисперсии используйте функцию ниже (она же использовалась на практике):

# https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#testing-for-overdispersioncomputing-overdispersion-factor
overdisp_fun <- function(model) {
  rdf <- df.residual(model)
  if (inherits(model, 'negbin'))
    rdf <- rdf - 1
  rp <- residuals(model, type = 'pearson')
  Pearson.chisq <- sum(rp ^ 2)
  prat <- Pearson.chisq / rdf
  pval <-
    pchisq(Pearson.chisq, df = rdf, lower.tail = FALSE)
  c(
    chisq = Pearson.chisq,
    ratio = prat,
    rdf = rdf,
    p = pval
  )
}

Для самопроверки в поле ниже введите значение ratio, полученное при проверке допущения об отсутствии сверхдисперсии, округленное до тысячных. В качестве десятичного разделителя используйте точку.

Подсказки
  • Для проверки линейности связи необходимо построить график зависимости пирсоновских остатков от предсказанных моделью значений.
    • Предсказания получаются из модели с помощью функции predict() с указанием type = "response" без использования новых данных, так как нас интересуют значения предсказания для данных, на которых модель была подобрана.
    • Остатки получаются аналогично с помощью функции resid() с указанием type = "pearson", чтобы получить стандартизированные значения остатков.
    • Для отображения тренда используется geom_smooth() с методом "loess", так как линейный тренд здесь не работает из-за использования функции связи внутри модели.
  • При проверке сверхдисперсии нулевая гипотеза теста гласит, что в данных не наблюдается сверхдисперсии, а альтернативная — что сверхдисперсия присутствует. Соответствующим образом интерпретируется и результат тестирования.

#10

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

Напишите функцию prediction_metrics(), которая будет рассчитывать метрики предсказательной силы модели на основе значений целевой переменной из данных и предсказанных вероятностей. Функция должна:

  • принимать на вход:
    • вектор значений бинарной переменной из данных (data)
    • вектор предсказанных моделью вероятностей (response)
    • порог, использующийся для перевода непрерывных предсказаний в дискретные (threshold)
  • возвращать:
    • именованный вектор со значениями accuracy, precision, recall и F1-меры.

Пример работы функции на полной модели представлен ниже.

prediction_metrics(data = heart_train$TenYearCHD,
                   response = predict(model1, type = "response"),
                   threshold = 0.6)
  accuracy  precision     recall         F1 
0.84946646 0.92857143 0.03194103 0.06175772
Подсказки
  • Чтобы написать функцию, необходимо взять код, используемый для расчета метрика предсказательной силы модели и заменить используемые в нем переменные на соответствующие аргументы функции.
  • Чтобы создать именованный вектор, нужно задать соответствие между именем элемента вектора и его значением: c(name1 = value1, name2 = value2, ...).