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 тис. Фактор, тобто >50K, <=50K
серед інших
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)
Пояснення коду
- безперервний <- select_if(data_adult, is.numeric): використовуйте функцію select_if() з бібліотеки dplyr, щоб вибрати лише числові стовпці
- summary(continuous): друкує підсумкову статистику
вихід:
## 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 годин на тиждень.
Спостереження вище цього порогу можна відкинути. Ви використовуєте фільтр із dplyr бібліотека
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(), щоб передати функцію в усі стовпці набору даних. Ви зберігаєте результат у списку
- функція(x): функція буде оброблена для кожного 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) Розробка функцій
Переробити освіту
На графіку вище ви можете побачити, що змінна освіта має 16 рівнів. Це суттєво, і деякі рівні мають відносно низьку кількість спостережень. Якщо ви хочете покращити кількість інформації, яку ви можете отримати від цієї змінної, ви можете змінити її на вищий рівень. А саме, ви створюєте більші групи з однаковим рівнем освіти. Наприклад, низький рівень освіти перетвориться на відсівання. Вищі рівні освіти замінять на магістерську.
Ось подробиці:
| Старий рівень | Новий рівень |
|---|---|
| Дошкільна освіта | Викинути |
| 10th | Викинути |
| 11th | Викинути |
| 12th | Викинути |
| 1-4-й | Викинути |
| 5th-6th | Викинути |
| 7th-8th | Викинути |
| 9th | Викинути |
| HS-град | HighGrad |
| Деякий-коледж | Спільнота |
| доц.-акдм | Спільнота |
| доц | Спільнота |
| Бакалавр | Бакалавр |
| майстри | майстри |
| Проф-школа | майстри |
| Докторська ступінь | 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
Переробити Mariтал-статус
Також можна створити нижчі рівні для сімейного стану. У наступному коді ви змінюєте рівень таким чином:
| Старий рівень | Новий рівень |
|---|---|
| Ніколи не одружений | НЕ одружений |
| Одружений-дружина-відсутній | НЕ одружений |
| Одружений-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= hours.per.week)): для графіка щільності потрібна лише одна змінна
- geom_density(aes(колір = освіта), 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(колір= дохід), розмір =0.5): побудуйте точкову діаграму
- stat_smooth(): додайте лінію тренду з такими аргументами:
- method='lm': побудувати підігнане значення, якщо лінійна регресія
- формула = y~poly(x,2): Підгонка поліноміальної регресії
- se = TRUE: додати стандартну помилку
- aes(колір= дохід): Розбити модель за доходом
вихід:
У двох словах, ви можете перевірити умови взаємодії в моделі, щоб виявити ефект нелінійності між тижневим робочим часом та іншими функціями. Важливо визначити, за яких умов робочий час відрізняється.
Кореляція
Наступна перевірка — візуалізація кореляції між змінними. Тип рівня фактора перетворюється на числовий, щоб можна було побудувати теплову карту, яка містить коефіцієнт кореляції, обчислений за методом Спірмена.
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(). The Узагальнена лінійна модель це колекція моделей. Основний синтаксис:
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'): Підібрати логістичну модель (family = '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 (критерії інформації Akaike): це еквівалент 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) Оцініть продуктивність моделі
Матриця плутанини
Команда матриця плутанини є кращим вибором для оцінки ефективності класифікації порівняно з різними показниками, які ви бачили раніше. Загальна ідея полягає в тому, щоб підрахувати кількість разів, коли випадки True класифікуються як False.
Щоб обчислити матрицю плутанини, вам спочатку потрібно мати набір прогнозів, щоб їх можна було порівняти з фактичними цілями.
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, predict > 0.5): обчисліть матрицю плутанини. predict > 0.5 означає, що повертає 1, якщо передбачені ймовірності вище 0.5, інакше 0.
вихід:
## ## FALSE TRUE ## <=50K 6310 495 ## >50K 1074 1229
Кожен рядок у матриці плутанини представляє фактичну ціль, тоді як кожен стовпець представляє прогнозовану ціль. Перший рядок цієї матриці враховує дохід менше 50 тис. (клас False): 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)
Точність проти запам'ятовування
Точність дивиться на точність позитивного прогнозу. Згадувати – відношення позитивних випадків, які правильно виявлені класифікатором;
Ви можете створити дві функції для обчислення цих двох показників
- Точність конструкції
precision <- function(matrix) {
# True positive
tp <- matrix[2, 2]
# false positive
fp <- matrix[1, 2]
return (tp / (tp + fp))
}
Пояснення коду
- mat[1,1]: повертає першу комірку першого стовпця кадру даних, тобто справжнє позитивне
- mat[1,2]; Повертає першу комірку другого стовпця кадру даних, тобто помилковий результат
recall <- function(matrix) {
# true positive
tp <- matrix[2, 2]# false positive
fn <- matrix[2, 1]
return (tp / (tp + fn))
}
Пояснення коду
- mat[1,1]: повертає першу комірку першого стовпця кадру даних, тобто справжнє позитивне
- mat[2,1]; Повертає другу клітинку першого стовпця кадру даних, тобто хибне негативне значення
Ви можете перевірити свої функції
prec <- precision(table_mat) prec rec <- recall(table_mat) rec
вихід:
## [1] 0.712877 ## [2] 0.5336518
Коли модель каже, що це особа вище 50 тис., це правильно лише в 54 відсотках випадків і може претендувати на осіб вище 50 тис. у 72 відсотках випадків.
Ви можете створити 
оцінка на основі точності та запам'ятовування. The є гармонійним середнім цих двох показників, тобто надає більшої ваги меншим значенням.
f1 <- 2 * ((prec * rec) / (prec + rec)) f1
вихід:
## [1] 0.6103799
Компроміс між точністю та запам’ятовуванням
Неможливо одночасно мати високу точність і високу запам'ятовуваність.
Якщо ми підвищимо точність, правильну людину буде краще передбачити, але ми пропустимо багато з них (нижче запам’ятовування). У деяких ситуаціях ми віддаємо перевагу вищій точності, ніж відкликанню. Існує увігнутий зв’язок між точністю та пригадуванням.
- Уявіть, вам потрібно передбачити, чи є у пацієнта захворювання. Ви хочете бути максимально точним.
- Якщо вам потрібно виявити потенційних шахраїв на вулиці за допомогою розпізнавання облич, було б краще зловити багато людей, позначених як шахраї, навіть незважаючи на низьку точність. Поліція зможе відпустити не шахрайську особу.
Крива ROC
Команда приймач Operaхарактеристика Крива — ще один поширений інструмент, який використовується з бінарною класифікацією. Вона дуже схожа на криву точності/відкликання, але замість побудови графіка точності від відкликання крива ROC показує справжній позитивний коефіцієнт (тобто відкликання) проти хибного позитивного показника. Рівень хибнопозитивних результатів – це співвідношення негативних випадків, які неправильно класифікуються як позитивні. Він дорівнює одиниці мінус справжній негативний коефіцієнт. Справжньою негативною ставкою також називають специфічність. Звідси побудована крива ROC чутливість (пригадування) проти 1-специфічності
Щоб побудувати криву ROC, нам потрібно встановити бібліотеку RORC. Ми можемо знайти в conda бібліотека. Ви можете ввести код:
conda install -cr r-rocr –yes
Ми можемо побудувати ROC за допомогою функцій prediction() і performance().
library(ROCR) ROCRpred <- prediction(predict, data_test$income) ROCRperf <- performance(ROCRpred, 'tpr', 'fpr') plot(ROCRperf, colorize = TRUE, text.adj = c(-0.2, 1.7))
Пояснення коду
- prediction(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() | формула, дані, сімейство* |
| glm | Узагальніть модель | резюме() | приталена модель |
| база | Зробіть прогноз | передбачити() | підібрана модель, набір даних, тип = 'відповідь' |
| база | Створіть матрицю плутанини | стіл() | y, predict() |
| база | Створіть оцінку точності | sum(diag(table())/sum(table() | |
| РПЦР | Створення ROC: Крок 1 Створення прогнозу | прогноз() | predict(), y |
| РПЦР | Створення ROC: Крок 2 Створення продуктивності | продуктивність() | prediction(), 'tpr', 'fpr' |
| РПЦР | Створіть ROC: крок 3. Побудуйте графік | сюжет() | продуктивність() |
Інший GLM типи моделей:
– біном: (link = “logit”)
– Гаусс: (link = “ідентичність”)
– Гама: (link = “inverse”)
– inverse.gaussian: (посилання = “1/mu^2”)
– poisson: (link = “log”)
– квазі: (link = «тотожність», variance = «константа»)
– квазібіноміальний: (link = “logit”)
– quasipoisson: (link = “log”)
