HW12 // Обобщенные аддитивные модели. Регуляризованная регрессия

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

PISA

Сегодня мы возьмем данные, полученные в ходе программы PISA (Programme for International Student Assessment, 2006), дополненные рядом показателей из данных ООН. Переменные, которые нас будут интересовать:

  • Overall Science Score (average score for 15 year olds)
  • Interest in science
  • Identifying scientific Issues
  • Explaining phenomena scientifically
  • Support for scientific inquiry
  • Income Index
  • Health Index
  • Education Index
  • Human Development Index (HDI, composed of the Income index, Health Index, and Education Index)

Также в данных есть переменная Country. Так как у нас усредненные показатели по странам, то эта переменная выступает в данном случае в качестве ID (идентификатора наблюдения), поэтому её мы не будем использовать в моделировании.

В качестве целевой переменной нас будет интересовать Overall, остальные мы будм рассматривать как предикторы.

#1

  1. Уже классически, загрузите данные и проверьте их структуру.
  2. Если в данных есть пропущенные значения, удалите наблюдения, которые их содержат.
  3. Визуализируйте взаимосвязи между переменными датасета (за исключением Country). Что можно сказать о характере связей между переменными?

Для визуализации в данном случае хорошо подойдет функция ggpairs() из пакета GGally.

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

Подсказки
  • Чтобы удалить строки с пропущенными значениями, воспользуйтей функцией drop_na().
  • Функция ggpairs() ожидает на вход датасет, состоящий только из тех переменных, взаимосвязи между которыми нас интересуют.

#2

  1. Постройте обобщенную аддитивную модель gam1, в которой целевой переменной будет Overall, а нелинейными предикторами — Interest, Support, Income и Health.
  2. Протестируйте гипотезу о статистической значимости модели в целом.

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

Подсказки
  • В GAM-модель нелинейный предикторы вводятся внутри сплайнов (функция s()).
  • Каждый нелинейный предиктор вводится внутри своего сплайна — s(pred1) + s(pred2) + s(pred3)
  • Тестирование гипотезы о значимости модели в целом проводится аналогично обобщенным линейным моделям — потребуется нулевая модель, содержащая в качестве предиктора только интерсепт.
  • Для сравнения построенной и нулевой модели используется функция anova() с аргументом test = "Chisq".

#3

  1. Протестируйте гипотезы о том, что связи между предикторами и целевой переменной построенной модели нелинейные.
  2. Если
    • результаты тестирования гипотез показывают, что связи между всеми предикторами модели и целевой переменной нелинейные, переходите к следующему заданию
    • результаты тестирования гипотез показывают, что в модели есть предикторы, линейно связанные с целевой переменной, то модифицируйте модель так, чтобы нелинейные предикторы были включены в модель внутри сплайнов, а линейные — вне сплайнов.
      • сохраните результат модификации в объект gam2
      • сравните с помощью статистического теста исходную и модифицированную модель
      • проверьте, изменились ли результаты тестов на нелинейность связи после модификации модели
Подсказки
  • Гипотезы о (не)линейности связей между предикторами и целевой переменной тестируются автоматически при выводе статистик модели с помощью функции summary().
  • GAM-модели позволяют вводить как нелинейные, так и линейные предикторы в одну и ту же модель. Общий синтаксис выглядит так (DV — целевая переменная, IV1 и IV2 — нелинейные предикторы, IV3 и IV4 — линейные предикторы):
gam(DV ~ s(IV1) + s(IV2) + IV3 + IV4, data = ds)
  • Сравнить модели друг с другом позволяет функция anova(..., test = "Chisq"), как и в случае с GLM.

Если в предыдущем задании не потребовалась модификация модели, то обобщенная аддитивная модель, с которой необходимо работать далее — это gam1. Если модификации была проведена, то обобщенная линейная модель, с которой необходимо работать далее — это gam2.

#4

  1. Постройте модель множественной линейной регрессии , в которую включены те же предикторы, что и в обобщенную аддитивную модель. Сравните результаты тестирования статистических гипотез для обобщенной аддитивной модели и модели множественной линейной регресси.
  2. Рассчитайте метрики MAPE и RMSE для обобщенной аддитивной модели и модели множественной линейной регрессии. На основании их значений сделайте вывод о том, какая из моделей лучше описывает закономерности, имеющиеся в данных.

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

Подсказки
  • В модели множественной линейной регресии, исходя из её математической формулировки, все предикторы считаются линейными.
  • Получить значения метрик качества можно с помощью соответствующих функций из пакета Metrics.
    • Функции ожидают в качестве первого аргумента реальные значения целевой переменной из данных, в качестве второго — модельные значения.
  • Модельные значения храняться в элементе fitted.values объекта, содержащего модель — общую линейную или обощенную аддитивную.

#5

Проведите графическую диагностику обобщенной аддитивной модели. Сделайты вывод о её качестве.

Подсказки
  • Для графической диагностики обобщенных аддитивных моделей используется функция gam.check().
  • Необходимо обратить внимание на распределение остатков и на зависимость остатков от предсказанных моделью значений.

#6

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

Интерпретация значений concurvity

На практике у препода произошло помутнение сознания и он нагло наврал о том, как интерпретируются значения concurvity.

Значения concurvity (любой из трёх метрик) варьируются от 0 до 1.

  • Если значение близко к нулю, что проблемы concurvity нет, предикторы независимы друг от друга.
  • Если значение близко к единице, что проблема concurvity есть, один предиктор может быть выражен через другой с помощью какой-либо функции.
Подсказки
  • Для проверки наличия или отсутствия concurvity, используется функция concurvity(), которая в качестве аргумента принимает обобщенную аддитивную модель.
  • Для интерпретации чаще всего используются значения worst и/или observed.

#7

Неким заказчиком нам велено изучить зависимость Overall от Issues, Explain, Evidence и HDI. Мы уже строили графики, по которым можно заподозрить нечто неладное относительно этих переменных, однако давайте для большей наглядности построим корреляционную матрицу между всеми этими переменными — Overall, Issues, Explain, Evidence и HDI. Сделайте вывод о том, можно ли моделировать зависимости в помощью линейных моделей и как это корректнее сделать.

#8

Смоделируйте зависимость Overall от предикторов Issues, Explain, Evidence и HDI с помощью LASSO-регрессии.

  1. Создайте вектор Y, содержащий значения целевой переменной, и матрицу X со значениями предикторов.
  2. Постройте модель LASSO-регрессии и выведите значения коэффициентов при минимальном значении штрафного коэффициента.
  3. Есть ли в модели предикторы, коэффициенты при которых были обнулены в ходе регуляризации?

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

Подсказки
  • Извлечь вектор значений целевой переменной можно через индексацию колонки с помощью оператора $ или с помощью функции pull() из tidyverse.
  • Преобразовать тиббл со значениями предикторов в матрицу можно с помощью функции as.matrix().
  • Модель регуляризованной регрессии строится с помощью функции cv.glmnet().
  • Вид регуляризации задается через аргумент alpha: 0 — ridge-регрессия, 1 — LASSO-регрессия.
  • Вывести значения коэффициентов модели поможет функция coef(), а её аргумент s позволит задать значения штрафного коэффициента — "lambda.min" или "lambda.1se".

#9

Смоделируйте зависимость Overall от Issues, Explain, Evidence и HDI с помощью Ridge-регрессии, оставив только те предикторы, коэффициенты при которых сохранили ненулевое значения в предыдущем задании.

  • Те предикторы, коэффициенты при которых обнулились в LASSO-регрессии не должны быть включены в модель Ridge-регрессии.
  • Если ни один из коэффициентов не был обнулён, то в модель включаются все предикторы.
Подсказки
  • Если какие-либо коэффициент в LASSO-регрессии получили нулевые значения коэффициентов, то их необходимо исключить из матрицы предикторов.
    • Можно создать новую матрицу из исходного датасета аналогично тому, как создавалась матрица X в задании 8a, а можно с помощью индексации убрать колонки из уже созданной матрицы X.
    • Новую матрицу предиктороров нужно передать в функцию cv.glmnet() с указанием аргумента alpha = 0, чтобы была построена именно ridge-регрессия.

#10

  1. Постройте модель множественной линейной регресии (без взаимодействия), соответствующую по структуре модели ridge-регрессии из предыдущего задания.
  2. Рассчитайте метрики качества модели MAPE и RMSE для построенной модели множественной линейной регрессии и регуляризованной регресии (ridge-регрессии) с минимальным значением штрафного коэффициента. Какая из моделей лучше описывает закономерность данных? С чем это может быть связано?

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

Подсказки
  • Для расчета метрик подойдут функции из пакета Metrics.
  • Значения, предсказанные общей линейной моделью можно получить, вытащив элемент fitted.values из объекта, содержащего модель.
  • Для получения предсказаний регуляризованной модели необходимо воспользоваться функцией predict().
    • Так как коэффициенты этой модели зависят от значения штрафного коэффициента — а значит и предсказанные значения также будут зависеть от его значения — в функцию predict() через аргумент s необходимо указать значения штрафного коэффициента, которое будет использовано для вычисления модельных значений. По заданию это будет s = "lambda.min".