GLM в R: обобщенная линейная модель с примером
Что такое логистическая регрессия?
Логистическая регрессия используется для прогнозирования класса, то есть вероятности. Логистическая регрессия может точно предсказать бинарный результат.
Представьте, что вы хотите предсказать, будет ли кредит отклонен/принят, основываясь на многих атрибутах. Логистическая регрессия имеет форму 0/1. y = 0, если кредит отклонен, y = 1, если принят.
Модель логистической регрессии отличается от модели линейной регрессии двумя способами.
- Прежде всего, логистическая регрессия принимает в качестве зависимой переменной только дихотомические (двоичные) входные данные (т. е. вектор 0 и 1).
- Во-вторых, результат измеряется следующей вероятностной функцией связи, называемой сигмоид благодаря своей S-образной форме.:
Выходные данные функции всегда находятся в диапазоне от 0 до 1. Проверьте изображение ниже.
Сигмовидная функция возвращает значения от 0 до 1. Для задачи классификации нам нужен дискретный выходной сигнал 0 или 1.
Чтобы преобразовать непрерывный поток в дискретное значение, мы можем установить границу решения на уровне 0.5. Все значения выше этого порога классифицируются как 1.
Как создать обобщенную модель лайнера (GLM)
Давайте использовать для взрослых набор данных для иллюстрации логистической регрессии. «Взрослый» — отличный набор данных для задачи классификации. Цель состоит в том, чтобы предсказать, превысит ли годовой доход человека в долларах 50.000 46,033. Набор данных содержит XNUMX XNUMX наблюдения и десять функций:
- возраст: возраст человека. Числовой
- образование: Образовательный уровень человека. Фактор.
- семейное положение: 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)
Код Пояснение
- 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. Постройте график распределения часов в неделю.
- 2: Стандартизировать непрерывные переменные
- Постройте распределение
Давайте подробнее рассмотрим распределение часов в неделю.
# 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.
Мы вычисляем верхний 2-процентный процентиль
top_one_percent <- quantile(data_adult$hours.per.week, .99) top_one_percent
Код Пояснение
- quantile(data_adult$hours.per.week, .99): вычислите значение 99 процентов рабочего времени.
Вывод:
## 99% ## 80
98 процентов населения работает менее 80 часов в неделю.
Вы можете отбросить наблюдения выше этого порога. Вы используете фильтр из дплир библиотека.
data_adult_drop <-data_adult %>% filter(hours.per.week<top_one_percent) dim(data_adult_drop)
Вывод:
## [1] 45537 10
- Стандартизируйте непрерывные переменные
Вы можете стандартизировать каждый столбец, чтобы повысить производительность, поскольку ваши данные не имеют одинакового масштаба. Вы можете использовать функцию 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)
Код Пояснение
- 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)
Код Пояснение
- data.frame(select_if(data_adult, is.factor)): мы сохраняем столбцы фактора в факторе в типе фрейма данных. Библиотеке ggplot2 требуется объект фрейма данных.
Вывод:
## [1] 6
Набор данных содержит 6 категориальных переменных.
Второй шаг более квалифицированный. Вы хотите построить гистограмму для каждого столбца фактора фрейма данных. Удобнее автоматизировать процесс, особенно если столбцов много.
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)))
Код Пояснение
- 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 уровней. Это существенно, и на некоторых уровнях наблюдается относительно небольшое количество наблюдений. Если вы хотите увеличить объем информации, которую вы можете получить от этой переменной, вы можете перевести ее на более высокий уровень. А именно, вы создаете более крупные группы со схожим уровнем образования. Например, низкий уровень образования будет конвертироваться в отсев. Высший уровень образования будет изменен на магистратуру.
Вот деталь:
Старый уровень | Новый уровень |
---|---|
дошкольный | отсев |
10 | отсев |
11 | отсев |
12 | отсев |
1st-4th | отсев |
5th-6th | отсев |
7th-8th | отсев |
9 | отсев |
ВШ-Град | ХайГрад |
Несколько колледжей некоторые колледжи | Сообщество |
Ассоц-АКДМ | Сообщество |
Доцент | Сообщество |
Бакалавры | Бакалавры |
Магистратура | Магистратура |
Профшкола | Магистратура |
Докторская степень | Аспирантура (PhD) |
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")))))))
Код Пояснение
- Мы используем глагол 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-статус
Также возможно создание более низких уровней семейного положения. В следующем коде вы меняете уровень следующим образом:
Старый уровень | Новый уровень |
---|---|
Никогда не был женат | Не женат не замужем |
Женат-супруга-отсутствует | Не женат не замужем |
Женат-AF-супруга | Женат |
Женат-гражданский-супруг | |
Разделенный | Разделенный |
Разведенный | |
Вдовы | вдова |
# 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()
Вывод:
Ящиковая диаграмма подтверждает, что распределение рабочего времени соответствует разным группам. На коробчатой диаграмме наблюдения обоих полов не являются однородными.
Проверить плотность еженедельного рабочего времени можно по видам образования. В дистрибутивах есть много разных вариантов. Вероятно, это можно объяснить типом контракта в США.
# Plot distribution working time by education ggplot(recast_data, aes(x = hours.per.week)) + geom_density(aes(color = education), alpha = 0.5) + theme_classic()
Код Пояснение
- 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()
Код Пояснение
- 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")
Код Пояснение
- 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(). Обобщенная линейная модель представляет собой коллекцию моделей. Основной синтаксис:
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)
Код Пояснение
- формула <- доход ~ .: Создайте подходящую модель
- 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"
Каждое значение можно извлечь с помощью знака $, за которым следует имя показателя. Например, вы сохранили модель как logit. Чтобы извлечь критерии 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
Код Пояснение
- 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 тыс. (Ложный класс): 6241 человек были правильно классифицированы как лица с доходом ниже 50 тыс. (Истинный отрицательный), а оставшийся был ошибочно отнесен к категории выше 50 тыс. (Ложно положительный). Во второй строке учитываются доходы выше 50 тыс., положительный класс — 1229 (Истинный положительный), в то время как Истинный отрицательный был 1074.
Вы можете рассчитать модель точность путем суммирования истинно положительных + истинно отрицательных результатов по общему наблюдению
accuracy_Test <- sum(diag(table_mat)) / sum(table_mat) accuracy_Test
Код Пояснение
- sum(diag(table_mat)): Сумма диагонали
- sum(table_mat): сумма матрицы.
Вывод:
## [1] 0.8277339
Похоже, у модели есть одна проблема: она переоценивает количество ложноотрицательных результатов. Это называется Парадокс теста точности. Мы заявили, что точность — это отношение правильных предсказаний к общему числу случаев. Мы можем иметь относительно высокую точность, но бесполезную модель. Это происходит, когда есть доминирующий класс. Если вы снова посмотрите на матрицу путаницы, вы увидите, что большинство случаев классифицируются как истинно отрицательные. Представьте себе, что модель классифицировала все классы как отрицательные (т.е. ниже 50 тыс.). У вас будет точность 75 процентов (6718/6718+2257). Ваша модель работает лучше, но с трудом различает истинное положительное от истинно отрицательного.
В такой ситуации предпочтительнее иметь более краткую метрику. Мы можем посмотреть:
- Точность=TP/(TP+FP)
- Напомним=TP/(TP+FN)
Точность против отзыва
Точность смотрит на точность положительного прогноза. Recall — доля положительных экземпляров, правильно обнаруженных классификатором;
Вы можете создать две функции для вычисления этих двух показателей.
- Точность конструкции
precision <- function(matrix) { # True positive tp <- matrix[2, 2] # false positive fp <- matrix[1, 2] return (tp / (tp + fp)) }
Код Пояснение
- mat[1,1]: возвращает первую ячейку первого столбца кадра данных, т. е. истинное положительное значение.
- мат[1,2]; Вернуть первую ячейку второго столбца кадра данных, т.е. ложное срабатывание.
recall <- function(matrix) { # true positive tp <- matrix[2, 2]# false positive fn <- matrix[2, 1] return (tp / (tp + fn)) }
Код Пояснение
- mat[1,1]: возвращает первую ячейку первого столбца кадра данных, т. е. истинное положительное значение.
- мат[2,1]; Вернуть вторую ячейку первого столбца фрейма данных, т.е. ложноотрицательный результат.
Вы можете протестировать свои функции
prec <- precision(table_mat) prec rec <- recall(table_mat) rec
Вывод:
## [1] 0.712877 ## [2] 0.5336518
Когда модель утверждает, что речь идет о человеке с доходом выше 50 тысяч, она верна только в 54 процентах случаев и может претендовать на людей с доходом выше 50 тысяч в 72 процентах случаев.
Вы можете создать Оценка на основе точности и отзыва.
является гармоническим средним значением этих двух показателей, то есть придает больший вес более низким значениям.
f1 <- 2 * ((prec * rec) / (prec + rec)) f1
Вывод:
## [1] 0.6103799
Компромисс между точностью и отзывом
Невозможно одновременно иметь высокую точность и высокую полноту.
Если мы повысим точность, правильную личность можно будет предсказать лучше, но мы пропустим многие из них (меньшая запоминаемость). В некоторых ситуациях мы предпочитаем более высокую точность, чем полноту. Существует вогнутая связь между точностью и отзывом.
- Представьте, вам нужно предсказать, есть ли у пациента заболевание. Вы хотите быть как можно более точными.
- Если вам нужно обнаружить потенциальных мошенников на улице с помощью распознавания лиц, было бы лучше поймать много людей, помеченных как мошенники, даже если точность низкая. Полиция сможет отпустить немошеннического человека.
Кривая ROC
В Получатель OperaТинг Характеристика Кривая — еще один распространенный инструмент, используемый при бинарной классификации. Она очень похожа на кривую точности/отзыва, но вместо построения графика зависимости точности от полноты кривая ROC показывает уровень истинного положительного результата (т. е. отзыва) в сравнении с уровнем ложноположительного результата. Уровень ложноположительных результатов — это соотношение отрицательных случаев, которые ошибочно классифицируются как положительные. Он равен единице минус истинно отрицательная ставка. Истинно отрицательную ставку еще называют специфичность. Следовательно, графики кривой ROC чувствительность (напомним) против 1-специфичности
Чтобы построить кривую ROC, нам нужно установить библиотеку RORC. Мы можем найти в Конде библиотека. Вы можете ввести код:
conda install -cr r-rocr –да
Мы можем построить 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))
Код Пояснение
- Predict(predict, data_test$income): Библиотеке ROCR необходимо создать объект прогнозирования для преобразования входных данных.
- Performance(ROCRpred, 'tpr','fpr'): возвращает две комбинации, которые нужно создать на графике. Здесь построены tpr и fpr. Чтобы построить график точности и отзыва вместе, используйте «prec», «rec».
Вывод:
Шаг 8) Улучшить модель
Можно попробовать добавить в модель нелинейности за счет взаимодействия между
- возраст и часы в неделю
- пол и количество часов в неделю.
Вам необходимо использовать тест оценки для сравнения обеих моделей.
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
Оценка немного выше предыдущей. Вы можете продолжать работать с данными, пытаясь побить рекорд.
Итого
Мы можем суммировать функцию для обучения логистической регрессии в таблице ниже:
Упаковка | Цель | Функция | Аргумент |
---|---|---|---|
– | Создать набор данных для обучения/тестирования | create_train_set() | данные, размер, поезд |
GLM | Обучите обобщенную линейную модель | глм() | формула, данные, семейство* |
GLM | Подведите итоги модели | резюме() | подогнанная модель |
Использование темпера с изогнутым основанием | Сделать прогноз | предсказать, () | подобранная модель, набор данных, тип = «ответ» |
Использование темпера с изогнутым основанием | Создайте матрицу путаницы | Таблица() | да, предсказать() |
Использование темпера с изогнутым основанием | Создать оценку точности | сумма(диаг(таблица())/сумма(таблица() | |
РПЦЗ | Создание ROC: Шаг 1. Создание прогноза | прогноз() | предсказать(), у |
РПЦЗ | Создание ROC: Шаг 2. Создание производительности | производительность() | предсказание(), 'tpr', 'fpr' |
РПЦЗ | Создание ROC: Шаг 3 Постройте график | участок() | производительность() |
Другой GLM Тип моделей:
– биномиальный: (ссылка = «логит»)
– гауссово: (ссылка = «идентичность»)
– Гамма: (ссылка = «инверсия»)
– inverse.gaussian: (ссылка = «1/mu^2»)
— Пуассон: (ссылка = «журнал»)
– квази: (ссылка = «идентичность», дисперсия = «константа»)
– квазибиномиальный: (ссылка = «логит»)
— квазипуассон: (ссылка = «журнал»)