Обобщенная линейная модель в R: обобщенная линейная модель и логистическая регрессия

Что такое обобщенная линейная модель (GLM) в R?

A Обобщенная линейная модель (GLM) Расширяет возможности обычной линейной регрессии, позволяя зависимой переменной следовать распределению, отличному от нормального. Rвы устанавливаете один со встроенным глм() функция из пакета stats.

Каждая обобщенная линейная модель (GLM) определяется тремя компонентами:

• Случайная составляющая: Вероятностное распределение зависимой переменной, взятое из семейства экспоненциальных распределений (биномиальное, Пуассона, гамма-распределение, гауссово распределение и другие).
• Систематический компонент: линейный предиктор, то есть взвешенная комбинация ваших объясняющих переменных.
• Функция ссылки: Функция, связывающая среднее значение отклика с линейным предиктором, например, логит-связь для бинарных данных или логарифмическая связь для количества событий.

Именно эта структура делает модель «обобщенной». Вместо того чтобы принудительно задавать результат, имеющий нормальное распределение, вы определяете правильное распределение посредством семья В качестве аргумента используется R, а коэффициенты оцениваются методом максимального правдоподобия. Логистическая регрессия — это просто обобщенная линейная модель с биномиальным распределением и логит-связью, поэтому она является естественной отправной точкой.

Что такое логистическая регрессия в R?

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

Представьте, что вы хотите предсказать, будет ли кредит отклонен/принят, основываясь на многих атрибутах. Логистическая регрессия имеет форму 0/1. y = 0, если кредит отклонен, y = 1, если принят.

Модель логистической регрессии отличается от модели линейной регрессии двумя способами.

• Прежде всего, логистическая регрессия принимает в качестве зависимой переменной только дихотомические (двоичные) входные данные (т. е. вектор 0 и 1).
• Во-вторых, результат отображается с помощью вероятностной функции связи, называемой сигмоид (логистическая) функция благодаря своей S-образной форме:

Выходные данные функции всегда находятся в диапазоне от 0 до 1. Проверьте изображение ниже.

Сигмовидная функция возвращает значения от 0 до 1. Для задачи классификации нам нужен дискретный выходной сигнал 0 или 1.

Чтобы преобразовать непрерывный поток в дискретное значение, мы можем установить границу решения на уровне 0.5. Все значения выше этого порога классифицируются как 1.

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

Обобщенная линейная модель против линейной регрессии: ключевые различия в R²

Прежде чем писать код, полезно точно знать, когда применяется стандарт. линейная регрессия Функция lm() больше не подходит, и её место должна занять функция glm().

Короче говоря, выбирайте lm(), когда результат представляет собой нормально распределенное измерение, и glm(), когда результат представляет собой решение «да/нет» — задачу, которую в противном случае можно было бы поручить другой функции. Древо решений классификатор, количество событий или строго положительная величина, разброс которой увеличивается пропорционально её среднему значению.

Как создать обобщенную линейную модель (GLM) в R

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

Давайте использовать для взрослых Набор данных для иллюстрации логистической регрессии. Набор данных «взрослый» отлично подходит для задачи классификации. Цель состоит в том, чтобы предсказать, превысит ли годовой доход человека в долларах США 50 000. Набор данных содержит 48 842 наблюдения и десять переменных:

• возраст: возраст человека. Числовой
• образование: Образовательный уровень человека. Фактор.
• семейное положение: MariТальный статус личности. Фактор, т. е. никогда не состоявший в браке, женатый гражданский супруг,…
• пол: Пол человека. Фактор, т.е. мужской или женский
• доход: Target переменная. Доход выше или ниже 50 тыс. Фактор, т.е. >50К, <=50К

среди других

library(dplyr)
data_adult <-read.csv("https://raw.githubusercontent.com/guru99-edu/R-Programming/master/adult.csv")
glimpse(data_adult)

Выход:

Observations: 48,842
Variables: 10
$ x               <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,...
$ age             <int> 25, 38, 28, 44, 18, 34, 29, 63, 24, 55, 65, 36, 26...
$ workclass       <fctr> Private, Private, Local-gov, Private, ?, Private,...
$ education       <fctr> 11th, HS-grad, Assoc-acdm, Some-college, Some-col...
$ educational.num <int> 7, 9, 12, 10, 10, 6, 9, 15, 10, 4, 9, 13, 9, 9, 9,...
$ marital.status  <fctr> Never-married, Married-civ-spouse, Married-civ-sp...
$ race            <fctr> Black, White, White, Black, White, White, Black, ...
$ gender          <fctr> Male, Male, Male, Male, Female, Male, Male, Male,...
$ hours.per.week  <int> 40, 50, 40, 40, 30, 30, 40, 32, 40, 10, 40, 40, 39...
$ income          <fctr> <=50K, <=50K, >50K, >50K, <=50K, <=50K, <=50K, >5...

Мы будем действовать следующим образом:

• Шаг 1. Проверьте непрерывные переменные
• Шаг 2: Проверьте переменные фактора
• Шаг 3. Разработка функций
• Шаг 4. Сводная статистика
• Шаг 5: Обучение/тестовый набор
• Шаг 6: Постройте модель
• Шаг 7: Оцените производительность модели
• Шаг 8. Улучшите модель

Ваша задача — предсказать, какой человек будет иметь доход выше 50 тысяч.

В этом руководстве будет подробно описан каждый шаг для выполнения анализа реального набора данных.

Шаг 1) Проверьте непрерывные переменные

На первом этапе вы можете увидеть распределение непрерывных переменных.

continuous <-select_if(data_adult, is.numeric)
summary(continuous)

Code объяснение

• Continuous <- select_if(data_adult, is.numeric): используйте функцию select_if() из библиотеки dplyr, чтобы выбрать только числовые столбцы.
• summary(continual): печать сводной статистики.

Выход:

##        X              age        educational.num hours.per.week 
##  Min.   :    1   Min.   :17.00   Min.   : 1.00   Min.   : 1.00  
##  1st Qu.:11509   1st Qu.:28.00   1st Qu.: 9.00   1st Qu.:40.00  
##  Median :23017   Median :37.00   Median :10.00   Median :40.00  
##  Mean   :23017   Mean   :38.56   Mean   :10.13   Mean   :40.95  
##  3rd Qu.:34525   3rd Qu.:47.00   3rd Qu.:13.00   3rd Qu.:45.00  
##  Max.   :46033   Max.   :90.00   Max.   :16.00   Max.   :99.00	

Из приведенной выше таблицы видно, что данные имеют совершенно разные масштабы, а часы.в.недели имеют большие выбросы (т.е. посмотрите на последний квартиль и максимальное значение).

С этим можно справиться, выполнив два шага:

• Постройте график распределения часов в неделю.
• Стандартизируйте непрерывные переменные

1. Постройте распределение

Давайте подробнее рассмотрим распределение часов в неделю.

# Histogram with kernel density curve
library(ggplot2)
ggplot(continuous, aes(x = hours.per.week)) +
    geom_density(alpha = .2, fill = "#FF6666")

Выход:

Переменная имеет множество выбросов и нечеткое распределение. Частично решить эту проблему можно, удалив верхние 0.01 процента часов в неделю.

Основной синтаксис квантиля:

quantile(variable, percentile)
arguments:
-variable:  Select the variable in the data frame to compute the percentile
-percentile:  Can be a single value between 0 and 1 or multiple value. If multiple, use this format:  `c(A,B,C, ...)
- `A`,`B`,`C` and `...` are all integer from 0 to 1.

Мы вычисляем 99-й процентиль еженедельного рабочего времени.

top_one_percent <- quantile(data_adult$hours.per.week, .99)
top_one_percent

Code объяснение

• quantile(data_adult$hours.per.week, .99): Вычисляет 99-й процентиль еженедельного рабочего времени.

Выход:

## 99% 
##  80

99 процентов населения работает менее 80 часов в неделю.

Вы можете отбросить наблюдения выше этого порога. Вы используете фильтр из дплир библиотека.

data_adult_drop <-data_adult %>%
filter(hours.per.week<top_one_percent)
dim(data_adult_drop)

Выход:

## [1] 45537    10

2. Стандартизируйте непрерывные переменные

Вы можете стандартизировать каждый столбец, чтобы повысить производительность, поскольку ваши данные не имеют одинакового масштаба. Вы можете использовать функцию mutate_if из библиотеки dplyr. Основной синтаксис:

mutate_if(df, condition, funs(function))
arguments:
-`df`: Data frame used to compute the function
- `condition`: Statement used. Do not use parenthesis
- funs(function):  Return the function to apply. Do not use parenthesis for the function

Вы можете стандартизировать числовые столбцы следующим образом:

data_adult_rescale <- data_adult_drop %>%
	mutate_if(is.numeric, funs(as.numeric(scale(.))))
head(data_adult_rescale)

Code объяснение

• mutate_if(is.numeric, funs(scale)): условие — только числовой столбец, а функция — масштаб.

Выход:

##           X         age        workclass    education educational.num
## 1 -1.732680 -1.02325949          Private         11th     -1.22106443
## 2 -1.732605 -0.03969284          Private      HS-grad     -0.43998868
## 3 -1.732530 -0.79628257        Local-gov   Assoc-acdm      0.73162494
## 4 -1.732455  0.41426100          Private Some-college     -0.04945081
## 5 -1.732379 -0.34232873          Private         10th     -1.61160231
## 6 -1.732304  1.85178149 Self-emp-not-inc  Prof-school      1.90323857
##       marital.status  race gender hours.per.week income
## 1      Never-married Black   Male    -0.03995944  <=50K
## 2 Married-civ-spouse White   Male     0.86863037  <=50K
## 3 Married-civ-spouse White   Male    -0.03995944   >50K
## 4 Married-civ-spouse Black   Male    -0.03995944   >50K
## 5      Never-married White   Male    -0.94854924  <=50K
## 6 Married-civ-spouse White   Male    -0.76683128   >50K

Шаг 2) Проверьте факторные переменные

Этот шаг преследует две цели:

• Проверьте уровень в каждом категориальном столбце
• Определите новые уровни

Мы разделим этот шаг на три части:

• Выберите категориальные столбцы
• Сохраните гистограмму каждого столбца в списке.
• Распечатайте графики

Мы можем выбрать столбцы коэффициентов с помощью кода ниже:

# Select categorical column
factor <- data.frame(select_if(data_adult_rescale, is.factor))
	ncol(factor)

Code объяснение

• data.frame(select_if(data_adult, is.factor)): мы сохраняем столбцы фактора в факторе в типе фрейма данных. Библиотеке ggplot2 требуется объект фрейма данных.

Выход:

## [1] 6

Набор данных содержит 6 категориальных переменных.

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

library(ggplot2)
# Create graph for each column
graph <- lapply(names(factor),
    function(x) 
	ggplot(factor, aes(get(x))) +
		geom_bar() +
		theme(axis.text.x = element_text(angle = 90)))

Code объяснение

• lapply(): используйте функцию lapply() для передачи функции во все столбцы набора данных. Вы сохраняете вывод в списке
• function(x): функция будет обрабатываться для каждого x. Здесь х — столбцы
• ggplot(factor, aes(get(x))) + geom_bar()+ theme(axis.text.x = element_text(angle = 90)): создайте гистограмму для каждого элемента x. Обратите внимание: чтобы вернуть x как столбец, вам нужно включить его в get().

Последний шаг относительно прост. Вы хотите распечатать 6 графиков.

# Print the graph
graph

Выход:

## [[1]]

## ## [[2]]

## ## [[3]]

## ## [[4]]

## ## [[5]]

## ## [[6]]

Примечание. Используйте кнопку «Далее», чтобы перейти к следующему графику.

Шаг 3) Разработка функций

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

Переделать образование

На графике выше видно, что переменная education имеет 16 уровней. Это существенно, и на некоторых уровнях наблюдается относительно небольшое количество наблюдений. Если вы хотите увеличить объем информации, которую вы можете получить от этой переменной, вы можете перевести ее на более высокий уровень. А именно, вы создаете более крупные группы со схожим уровнем образования. Например, низкий уровень образования будет конвертироваться в отсев. Высший уровень образования будет изменен на магистратуру.

Вот деталь:

recast_data <- data_adult_rescale %>%
	select(-X) %>%
	mutate(education = factor(ifelse(education == "Preschool" | education == "10th" | education == "11th" | education == "12th" | education == "1st-4th" | education == "5th-6th" | education == "7th-8th" | education == "9th", "dropout", ifelse(education == "HS-grad", "HighGrad", ifelse(education == "Some-college" | education == "Assoc-acdm" | education == "Assoc-voc", "Community",
    ifelse(education == "Bachelors", "Bachelors",
        ifelse(education == "Masters" | education == "Prof-school", "Master", "PhD")))))))

Code объяснение

• Мы используем глагол mutate из библиотеки dplyr. Мы меняем ценности образования с помощью утверждения ifelse

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

recast_data %>%
	group_by(education) %>%
	summarize(average_educ_year = mean(educational.num),
		count = n()) %>%
	arrange(average_educ_year)

Выход:

## # A tibble: 6 x 3
## education average_educ_year count			
##      <fctr>             <dbl> <int>
## 1   dropout       -1.76147258  5712
## 2  HighGrad       -0.43998868 14803
## 3 Community        0.09561361 13407
## 4 Bachelors        1.12216282  7720
## 5    Master        1.60337381  3338
## 6       PhD        2.29377644   557

Переделать Marital-статус

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

# Change level marry
recast_data <- recast_data %>%
	mutate(marital.status = factor(ifelse(marital.status == "Never-married" | marital.status == "Married-spouse-absent", "Not_married", ifelse(marital.status == "Married-AF-spouse" | marital.status == "Married-civ-spouse", "Married", ifelse(marital.status == "Separated" | marital.status == "Divorced", "Separated", "Widow")))))

Вы можете проверить количество людей в каждой группе.

table(recast_data$marital.status)

Выход:

## ##     Married Not_married   Separated       Widow
##       21165       15359        7727        1286

Шаг 4) Сводная статистика

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

# Plot gender income
ggplot(recast_data, aes(x = gender, fill = income)) +
    geom_bar(position = "fill") +
    theme_classic()

Выход:

Далее проверьте, влияет ли происхождение человека на его заработок.

# Plot origin income
ggplot(recast_data, aes(x = race, fill = income)) +
    geom_bar(position = "fill") +
    theme_classic() +
    theme(axis.text.x = element_text(angle = 90))

Выход:

Количество часов работы по полу.

# box plot gender working time
ggplot(recast_data, aes(x = gender, y = hours.per.week)) +
    geom_boxplot() +
    stat_summary(fun.y = mean,
        geom = "point",
        size = 3,
        color = "steelblue") +
    theme_classic()

Выход:

Ящиковая диаграмма подтверждает, что распределение рабочего времени соответствует разным группам. На коробчатой ​​диаграмме наблюдения обоих полов не являются однородными.

Вы можете проверить плотность распределения еженедельного рабочего времени по типу образования. Распределения имеют множество различных пиков. Вероятно, это можно объяснить типом образования.tracт в США.

# Plot distribution working time by education
ggplot(recast_data, aes(x = hours.per.week)) +
    geom_density(aes(color = education), alpha = 0.5) +
    theme_classic()

Code объяснение

• ggplot(recast_data, aes( x= часы.в.неделю)): для графика плотности требуется только одна переменная.
• geom_density(aes(color = education), Alpha =0.5): геометрический объект для управления плотностью.

Выход:

Для подтверждения своих мыслей можно выполнить одностороннее ANOVA-тест:

anova <- aov(hours.per.week~education, recast_data)
summary(anova)

Выход:

##                Df Sum Sq Mean Sq F value Pr(>F)    
## education       5   1552  310.31   321.2 <2e-16 ***
## Residuals   45531  43984    0.97                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Тест ANOVA подтверждает разницу в средних показателях между группами.

Нелинейность

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

library(ggplot2)
ggplot(recast_data, aes(x = age, y = hours.per.week)) +
    geom_point(aes(color = income),
        size = 0.5) +
    stat_smooth(method = 'lm',
        formula = y~poly(x, 2),
        se = TRUE,
        aes(color = income)) +
    theme_classic()

Code объяснение

• ggplot(recast_data, aes(x = age, y =hours.per.week)): установите внешний вид графика.
• geom_point(aes(color= доход), размер =0.5): построить точечный график.
• stat_smooth(): добавьте линию тренда со следующими аргументами:
  • метод='lm': постройте подобранное значение, если линейная регрессия
  • формула = y~poly(x,2): соответствует полиномиальной регрессии.
  • se = TRUE: добавить стандартную ошибку.
  • aes(color= доход): Разбить модель по доходу.

Выход:

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

Корреляция

Следующая проверка — визуализация корреляции между переменными. Вы преобразуете тип уровня фактора в числовой, чтобы можно было построить тепловую карту, содержащую коэффициент корреляции, вычисленный с помощью метода Спирмена.

library(GGally)
# Convert data to numeric
corr <- data.frame(lapply(recast_data, as.integer))
# Plot the graphggcorr(corr,
    method = c("pairwise", "spearman"),
    nbreaks = 6,
    hjust = 0.8,
    label = TRUE,
    label_size = 3,
    color = "grey50")

Code объяснение

• data.frame(lapply(recast_data,as.integer)): Преобразовать данные в числовые.
• ggcorr() строит тепловую карту со следующими аргументами:
  • метод: метод вычисления корреляции
  • nbreaks = 6: количество перерывов
  • hjust = 0.8: контроль положения имени переменной на графике.
  • label = TRUE: добавить метки в центре окон.
  • label_size = 3: Размер меток
  • color = «grey50»): Цвет метки.

Выход:

Шаг 5) Обучение/тестовый набор

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

set.seed(1234)
create_train_test <- function(data, size = 0.8, train = TRUE) {
    n_row = nrow(data)
    total_row = size * n_row
    train_sample <- 1: total_row
    if (train == TRUE) {
        return (data[train_sample, ])
    } else {
        return (data[-train_sample, ])
    }
}
data_train <- create_train_test(recast_data, 0.8, train = TRUE)
data_test <- create_train_test(recast_data, 0.8, train = FALSE)
dim(data_train)

Выход:

## [1] 36429     9

dim(data_test)

Выход:

## [1] 9108    9

Шаг 6) Постройте модель

Чтобы оценить производительность алгоритма, используйте функцию glm() из пакета stats. Обобщенная линейная модель представляет собой коллекцию моделей. Основной синтаксис:

glm(formula, data=data, family=linkfunction()
Argument:
- formula:  Equation used to fit the model- data: dataset used
- Family:     - binomial: (link = "logit")			
- gaussian: (link = "identity")			
- Gamma:    (link = "inverse")			
- inverse.gaussian: (link = "1/mu^2")			
- poisson:  (link = "log")			
- quasi:    (link = "identity", variance = "constant")			
- quasibinomial:    (link = "logit")			
- quasipoisson: (link = "log")	

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

formula <- income~.
logit <- glm(formula, data = data_train, family = 'binomial')
summary(logit)

Code объяснение

• формула <- доход ~ .: Создайте подходящую модель
• logit <- glm(formula, data = data_train, Family = 'binomial'): Сопоставьте логистическую модель (семейство = 'binomial') с данными data_train.
• summary(logit): Распечатать сводку модели.

Выход:

## 
## Call:
## glm(formula = formula, family = "binomial", data = data_train)
## ## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.6456  -0.5858  -0.2609  -0.0651   3.1982  
## 
## Coefficients:
##                           Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                0.07882    0.21726   0.363  0.71675    
## age                        0.41119    0.01857  22.146  < 2e-16 ***
## workclassLocal-gov        -0.64018    0.09396  -6.813 9.54e-12 ***
## workclassPrivate          -0.53542    0.07886  -6.789 1.13e-11 ***
## workclassSelf-emp-inc     -0.07733    0.10350  -0.747  0.45499    
## workclassSelf-emp-not-inc -1.09052    0.09140 -11.931  < 2e-16 ***
## workclassState-gov        -0.80562    0.10617  -7.588 3.25e-14 ***
## workclassWithout-pay      -1.09765    0.86787  -1.265  0.20596    
## educationCommunity        -0.44436    0.08267  -5.375 7.66e-08 ***
## educationHighGrad         -0.67613    0.11827  -5.717 1.08e-08 ***
## educationMaster            0.35651    0.06780   5.258 1.46e-07 ***
## educationPhD               0.46995    0.15772   2.980  0.00289 ** 
## educationdropout          -1.04974    0.21280  -4.933 8.10e-07 ***
## educational.num            0.56908    0.07063   8.057 7.84e-16 ***
## marital.statusNot_married -2.50346    0.05113 -48.966  < 2e-16 ***
## marital.statusSeparated   -2.16177    0.05425 -39.846  < 2e-16 ***
## marital.statusWidow       -2.22707    0.12522 -17.785  < 2e-16 ***
## raceAsian-Pac-Islander     0.08359    0.20344   0.411  0.68117    
## raceBlack                  0.07188    0.19330   0.372  0.71001    
## raceOther                  0.01370    0.27695   0.049  0.96054    
## raceWhite                  0.34830    0.18441   1.889  0.05894 .  
## genderMale                 0.08596    0.04289   2.004  0.04506 *  
## hours.per.week             0.41942    0.01748  23.998  < 2e-16 ***
## ---## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## ## (Dispersion parameter for binomial family taken to be 1)
## ##     Null deviance: 40601  on 36428  degrees of freedom
## Residual deviance: 27041  on 36406  degrees of freedom
## AIC: 27087
## 
## Number of Fisher Scoring iterations: 6

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

• AIC (информационные критерии Акаике): это эквивалент R2 в логистической регрессии. Он измеряет соответствие, когда к количеству параметров применяется штраф. Меньший АПК значения указывают, что модель ближе к истине.
• Нулевое отклонение: соответствует модели только с пересечением. Степень свободы n-1. Мы можем интерпретировать его как значение хи-квадрат (подогнанное значение, отличное от фактического значения, полученного при проверке гипотезы).
• Остаточное отклонение: Модель со всеми переменными. Это также интерпретируется как проверка гипотезы Хи-квадрат.
• Количество итераций оценки Фишера: количество итераций до сходимости.

Вывод функции glm() сохраняется в списке. В приведенном ниже коде показаны все элементы, доступные в переменной logit, которую мы создали для оценки логистической регрессии.

# Список очень длинный, выведите только первые три элемента

lapply(logit, class)[1:3]

Выход:

## $coefficients
## [1] "numeric"
## 
## $residuals
## [1] "numeric"
## 
## $fitted.values
## [1] "numeric"

Каждое значение может бытьtracОбозначается знаком $, за которым следует название метрики. Например, вы сохранили модель как logit. Например,tracВ соответствии с критериями AIC, вы используете:

logit$aic

Выход:

## [1] 27086.65

Шаг 7) Оцените производительность модели

Матрица путаницы

матрица путаницы — лучший выбор для оценки эффективности классификации по сравнению с различными показателями, которые вы видели ранее. Общая идея состоит в том, чтобы подсчитать, сколько раз истинные экземпляры классифицируются как ложные.

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

predict <- predict(logit, data_test, type = 'response')
# confusion matrix
table_mat <- table(data_test$income, predict > 0.5)
table_mat

Code объяснение

• Predict(logit,data_test, type = 'response'): Вычислить прогноз на тестовом наборе. Установите type = 'response', чтобы вычислить вероятность ответа.
• table(data_test$income, прогноз > 0.5): вычислите матрицу путаницы. Predict > 0.5 означает, что он возвращает 1, если прогнозируемые вероятности превышают 0.5, в противном случае — 0.

Выход:

##        
##         FALSE TRUE
##   <=50K  6310  495
##   >50K   1074 1229	

Каждая строка в матрице ошибок представляет собой фактическую цель, а каждый столбец — прогнозируемую цель. Первая строка этой матрицы рассматривает доход ниже 50 000 (отрицательный класс): 6,310 наблюдений были правильно классифицированы как лица с доходом ниже 50 000 (Истинный отрицательный), при этом 495 были ошибочно классифицированы как имеющие более 50 тысяч (Ложно положительныйВо второй строке рассматривается доход свыше 50 000: 1,229 случаев были определены правильно (Истинный положительный), при этом 1,074 были пропущены (Ложный негатив).

Вы можете рассчитать модель точность путем суммирования истинно положительных + истинно отрицательных результатов по общему наблюдению

accuracy_Test <- sum(diag(table_mat)) / sum(table_mat)
accuracy_Test

Code объяснение

• sum(diag(table_mat)): Сумма диагонали
• sum(table_mat): сумма матрицы.

Выход:

## [1] 0.8277339

Модель, по-видимому, страдает от одной проблемы: она выдает слишком много ложных отрицательных результатов. Это называется Парадокс теста точностиМы утверждали, что точность — это отношение количества правильных прогнозов к общему числу случаев. Мы можем иметь относительно высокую точность, но при этом бесполезную модель. Это происходит, когда существует доминирующий класс. Если вы посмотрите на матрицу ошибок, вы увидите, что большинство случаев классифицируются как истинно отрицательные. Представьте теперь, что модель классифицирует каждое наблюдение как отрицательное (т.е. меньше 50 тыс.). Вы все равно получите точность около 75 процентов (6,805 / 9,108). Ваша модель работает лучше, но с трудом отличает истинно положительные результаты от истинно отрицательных.

В такой ситуации предпочтительнее иметь более краткую метрику. Мы можем посмотреть:

• Точность=TP/(TP+FP)
• Напомним=TP/(TP+FN)

Точность против отзыва

Точность смотрит на точность положительного прогноза. Recall — доля положительных экземпляров, правильно обнаруженных классификатором;

Вы можете создать две функции для вычисления этих двух показателей.

1. Точность конструкции

precision <- function(matrix) {
	# True positive
    tp <- matrix[2, 2]
	# false positive
    fp <- matrix[1, 2]
    return (tp / (tp + fp))
}

Code объяснение

• mat[1,1]: возвращает первую ячейку первого столбца кадра данных, т. е. истинное положительное значение.
• мат[1,2]; Вернуть первую ячейку второго столбца кадра данных, т.е. ложное срабатывание.

recall <- function(matrix) {
# true positive
    tp <- matrix[2, 2]# false positive
    fn <- matrix[2, 1]
    return (tp / (tp + fn))
}

Code объяснение

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

Вы можете протестировать свои функции

prec <- precision(table_mat)
prec
rec <- recall(table_mat)
rec

Выход:

## [1] 0.712877
## [2] 0.5336518

Внимательно прочтите эти два числа. Точность (Precision) составляет 0.71, поэтому, когда модель утверждает, что человек зарабатывает более 50 000 долларов, это верно в 71% случаев. Полнота (Recall) составляет 0.53, поэтому модель обнаруживает только 53% людей, которые действительно зарабатывают более 50 000 долларов.

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

f1 <- 2 * ((prec * rec) / (prec + rec))
f1

Выход:

## [1] 0.6103799

Компромисс между точностью и отзывом

Невозможно одновременно иметь высокую точность и высокую полноту.

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

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

Кривая ROC

Получатель OperaТинг Характеристика Кривая — еще один распространенный инструмент, используемый при бинарной классификации. Она очень похожа на кривую точности/отзыва, но вместо построения графика зависимости точности от полноты кривая ROC показывает уровень истинного положительного результата (т. е. отзыва) в сравнении с уровнем ложноположительного результата. Уровень ложноположительных результатов — это соотношение отрицательных случаев, которые ошибочно классифицируются как положительные. Он равен единице минус истинно отрицательная ставка. Истинно отрицательную ставку еще называют специфичность. Следовательно, графики кривой ROC чувствительность (напомним) против 1-специфичности

Для построения ROC-кривой необходимо установить пакет ROCR. Его можно найти в репозитории conda. библиотеки. Вы можете ввести код:

conda install -c r r-rocr --yes

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

library(ROCR)
ROCRpred <- prediction(predict, data_test$income)
ROCRperf <- performance(ROCRpred, 'tpr', 'fpr')
plot(ROCRperf, colorize = TRUE, text.adj = c(-0.2, 1.7))

Code объяснение

• Predict(predict, data_test$income): Библиотеке ROCR необходимо создать объект прогнозирования для преобразования входных данных.
• Performance(ROCRpred, 'tpr','fpr'): возвращает две комбинации, которые нужно создать на графике. Здесь построены tpr и fpr. Чтобы построить график точности и отзыва вместе, используйте «prec», «rec».

Выход:

Шаг 8) Улучшите модель

Можно попробовать добавить в модель нелинейности за счет взаимодействия между

• возраст и часы в неделю
• пол и количество часов в неделю.

Затем вы сравниваете показатель F1 обеих моделей.

formula_2 <- income~age: hours.per.week + gender: hours.per.week + .
logit_2 <- glm(formula_2, data = data_train, family = 'binomial')
predict_2 <- predict(logit_2, data_test, type = 'response')
table_mat_2 <- table(data_test$income, predict_2 > 0.5)
precision_2 <- precision(table_mat_2)
recall_2 <- recall(table_mat_2)
f1_2 <- 2 * ((precision_2 * recall_2) / (precision_2 + recall_2))
f1_2

Выход:

## [1] 0.6109181

Результат F1 немного выше предыдущего. Вы можете продолжить работу с данными и попытаться превзойти этот результат.

Как интерпретировать коэффициенты GLM и отношения шансов в R

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

Выполните следующие четыре шага.

1. Возведите коэффициенты в степень. Примените функцию exp() ко всем оценкам, чтобы логарифм отношения шансов стал мультипликативным отношением шансов.
2. Добавьте доверительный интервал. Чтобы получить 95-процентный интервал в той же шкале шансов, оберните функцию confint() в функцию exp().
3. Сравните каждое значение с 1. Коэффициент вероятности выше 1 увеличивает вероятность принадлежности к положительному классу, значение ниже 1 уменьшает ее, а значение, близкое к 1, означает, что предиктор мало что добавляет.
4. Проверьте статистическую значимость. Интерпретируйте только те предикторы, значение p которых в сводном отчете ниже выбранного вами порогового значения, обычно 0.05.

# Convert log-odds coefficients into odds ratios
odds_ratio <- exp(coef(logit))
round(odds_ratio, 3)

# Odds ratios with 95% confidence intervals
exp(cbind(OddsRatio = coef(logit), confint(logit)))

Чтение вывода. Коэффициент для количества часов в неделю в нашей модели равен 0.41942. Возведя его в степень, получаем exp(0.41942) = 1.52, что означает, что увеличение количества рабочих часов в неделю на одно стандартное отклонение увеличивает вероятность заработка более 50 000 примерно в 1.5 раза, при условии, что все остальные переменные остаются неизменными.

Отрицательные коэффициенты работают аналогично. marital.statusNot_married равно -2.50346, поэтому exp(-2.50346) = 0.08: вероятность того, что человек, не состоящий в браке, составляет около 8 процентов от вероятности того, что человек состоит в браке. Поскольку непрерывные предикторы были стандартизированы на шаге 1, описывайте изменения в единицах стандартного отклонения, а не в часах.

Примечание для других семей: Экспоненцированные коэффициенты представляют собой отношения шансов только в рамках биномиального распределения с логит-связью. При family = “poisson” и логарифмической связи те же значения exp() считываются как отношения скоростей.

GLM в R: Краткий справочник функций

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

Другой GLM В рамках семейного аргумента доступны следующие варианты семей:

• биномиальное распределение: (ссылка = “логит”)
• гауссова функция: (ссылка = “идентичность”)
• Гамма: (ссылка = “обратная”)
• inverse.gaussian: (link = “1/mu^2”)
• poisson: (link = “log”)
• квази: (ссылка = “тождественность”, дисперсия = “константа”)
• квазибиномиальное распределение: (ссылка = “логит”)
• квазипуассон: (ссылка = “лог”)

Часто задаваемые вопросы (FAQ)