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: Стандартизировать непрерывные переменные
  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.

Мы вычисляем верхний 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
  1. Стандартизируйте непрерывные переменные

Вы можете стандартизировать каждый столбец, чтобы повысить производительность, поскольку ваши данные не имеют одинакового масштаба. Вы можете использовать функцию 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 — доля положительных экземпляров, правильно обнаруженных классификатором;

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

  1. Точность конструкции
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».

Вывод:

Кривая ROC

Шаг 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»)

— Пуассон: (ссылка = «журнал»)

– квази: (ссылка = «идентичность», дисперсия = «константа»)

– квазибиномиальный: (ссылка = «логит»)

— квазипуассон: (ссылка = «журнал»)