GLM в R: Узагальнена лінійна модель та логістична регресія

Що таке узагальнена лінійна модель (GLM) в R?

A Узагальнена лінійна модель (GLM) розширює звичайну лінійну регресію, тому змінна відгуку може мати розподіл, відмінний від нормального. У R, ви підходите до нього з вбудованим glm() функція з пакета stats.

Кожен GLM визначається трьома компонентами:

• Випадкова складова: розподіл ймовірностей змінної відгуку, взятий з експоненціальної родини (біноміальний, пуассонівський, гамма-, гаусівський та інші).
• Систематичний компонент: лінійний предиктор, тобто зважена комбінація ваших пояснювальних змінних.
• Функція посилання: функція, яка пов'язує середнє значення відповіді з лінійним предиктором, наприклад, логістичний зв'язок для двійкових даних або логарифмічний зв'язок для підрахунків.

Саме ця структура робить модель «узагальненою». Замість того, щоб форсувати нормальний розподіл результату, ви оголошуєте правильний розподіл через сім'я аргумент, а R оцінює коефіцієнти методом максимальної правдоподібності. Логістична регресія — це просто GLM з біноміальною родиною та логістичним зв'язком, тому вона є природною відправною точкою.

Що таке логістична регресія в R?

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

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

Модель логістичної регресії відрізняється від моделі лінійної регресії двома способами.

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

Результат функції завжди між 0 і 1. Перевірте зображення нижче

Сигмоїдна функція повертає значення від 0 до 1. Для завдання класифікації нам потрібен дискретний вихід 0 або 1.

Щоб перетворити безперервний потік у дискретне значення, ми можемо встановити межу рішення на 0.5. Усі значення, що перевищують цей поріг, класифікуються як 1

Тепер, коли функція зв'язку зрозуміла, порівняйте узагальнену модель зі звичайною лінійною моделлю, яку ви вже знаєте.

GLM проти лінійної регресії: ключові відмінності в R

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

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

Як створити узагальнену лінійну модель (GLM) в R

Після того, як теорія визначена, решта цього посібника застосовує біноміальний GLM по черзі до реального набору даних.

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

• вік: вік особи. Числовий
• освіта: Освітній рівень особистості. Фактор.
• сімейний.стан: 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)

Code Пояснення

• безперервний <- 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. Побудуйте графік розподілу

Розглянемо ближче розподіл годин.на.тиждень

# 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 Пояснення

• квантиль(data_adult$hours.per.week, .99): Обчисліть 99-й процентиль тижневого робочого часу

вихід:

## 99% 
##  80

99 відсотків населення працює менше 80 годин на тиждень.

Спостереження вище цього порогу можна відкинути. Ви використовуєте фільтр із dplyr бібліотека

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 категоріальних змінних

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

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(), щоб передати функцію в усі стовпці набору даних. Ви зберігаєте результат у списку
• функція(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 рівнів. Це суттєво, і деякі рівні мають відносно низьку кількість спостережень. Якщо ви хочете покращити кількість інформації, яку ви можете отримати від цієї змінної, ви можете змінити її на вищий рівень. А саме, ви створюєте більші групи з однаковим рівнем освіти. Наприклад, низький рівень освіти перетвориться на відсівання. Вищі рівні освіти замінять на магістерську.

Ось подробиці:

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

Переробити Mariтал-статус

Також можна створити нижчі рівні для сімейного стану. У наступному коді ви змінюєте рівень таким чином:

# 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= 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()

Code Пояснення

• 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")

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'): Підібрати логістичну модель (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"

Кожне значення може бутиtracза допомогою знака $, після якого йде назва метрики. Наприклад, ви зберегли модель як логістичний запис. Щоб навести прикладtracЗгідно з критеріями 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

Code Пояснення

• 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 тисяч (негативний клас): 6,310 спостережень були правильно класифіковані як особи з доходом нижче 50 тисяч (Справжній негатив), тоді як 495 були помилково класифіковані як такі, що перевищують 50 тисяч (Хибно позитивний). У другому рядку враховується дохід понад 50 тисяч: 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)

Точність проти запам'ятовування

Точність дивиться на точність позитивного прогнозу. Згадувати – відношення позитивних випадків, які правильно виявлені класифікатором;

Ви можете створити дві функції для обчислення цих двох показників

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

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

Code Пояснення

• 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))
}

Code Пояснення

• mat[1,1]: повертає першу комірку першого стовпця кадру даних, тобто справжнє позитивне
• mat[2,1]; Повертає другу клітинку першого стовпця кадру даних, тобто хибне негативне значення

Ви можете перевірити свої функції

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

вихід:

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

Уважно прочитайте ці два числа. Точність становить 0.71, тому коли модель каже, що людина заробляє понад 50 тисяч, це правильно у 71 відсотку випадків. Нагадаємо, що модель виявляє лише 53 відсотки людей, які дійсно заробляють понад 50 тисяч.

Ви можете створити оцінка на основі точності та запам'ятовування. The є гармонійним середнім цих двох показників, тобто надає більшої ваги меншим значенням.

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

вихід:

## [1] 0.6103799

Компроміс між точністю та запам’ятовуванням

Неможливо одночасно мати високу точність і високу запам'ятовуваність.

Якщо ми підвищимо точність, правильну людину буде краще передбачити, але ми пропустимо багато з них (нижче запам’ятовування). У деяких ситуаціях ми віддаємо перевагу вищій точності, ніж відкликанню. Існує увігнутий зв’язок між точністю та пригадуванням.

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

Крива ROC

Команда приймач Operaхарактеристика Крива — ще один поширений інструмент, який використовується з бінарною класифікацією. Вона дуже схожа на криву точності/відкликання, але замість побудови графіка точності від відкликання крива ROC показує справжній позитивний коефіцієнт (тобто відкликання) проти хибного позитивного показника. Рівень хибнопозитивних результатів – це співвідношення негативних випадків, які неправильно класифікуються як позитивні. Він дорівнює одиниці мінус справжній негативний коефіцієнт. Справжньою негативною ставкою також називають специфічність. Звідси побудована крива ROC чутливість (пригадування) проти 1-специфічності

Щоб побудувати ROC-криву, нам потрібно встановити пакет під назвою ROCR. Ми можемо знайти його в конді бібліотека. Ви можете ввести код:

conda install -c r 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))

Code Пояснення

• prediction(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. Додайте довірчий інтервал. Оберніть confint() в exp(), щоб отримати 95-відсотковий інтервал на тій самій шкалі шансів.
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 тисяч приблизно на 1.5, зберігаючи всі інші змінні постійними.

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

Примітка для інших сімей: Коефіцієнти, піднесені до степеня, є відношеннями шансів лише для біноміальної родини з логістичним зв'язком. З родиною = "poisson" та логарифмічним зв'язком, ті ж значення exp() зчитуються як відношення швидкостей.

GLM в R: Короткий довідник функцій

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

Інший GLM Сім'ї, доступні через сімейний аргумент:

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

Поширені запитання