GLM în R: Model liniar generalizat și regresie logistică

Ce este un model liniar generalizat (GLM) în R?

A Model liniar generalizat (GLM) extinde regresia liniară obișnuită, astfel încât variabila de răspuns poate urma o distribuție diferită de cea normală. În R, montezi unul cu sistemul încorporat glm() funcție din pachetul de statistici.

Fiecare GLM este definit de trei componente:

• Componentă aleatorie: distribuția de probabilitate a variabilei de răspuns, preluată din familia exponențială (binomială, de Poisson, Gamma, gaussiană și altele).
• Componentă sistematică: predictorul liniar, adică combinația ponderată a variabilelor explicative.
• Funcția link: funcția care conectează media răspunsului la predictorul liniar, de exemplu legătura logit pentru date binare sau legătura logaritmică pentru numărătoare.

Această structură este cea care face ca modelul să fie „generalizat”. În loc să forțezi un rezultat distribuit normal, declari distribuția corectă prin . argument și R estimează coeficienții prin probabilitate maximă. Regresia logistică este pur și simplu un GLM cu o familie binomială și o legătură logit, deci este punctul de plecare natural.

Ce este regresia logistică în R?

Regresia logistică este folosită pentru a prezice o clasă, adică o probabilitate. Regresia logistică poate prezice cu acuratețe un rezultat binar.

Imaginați-vă că doriți să preziceți dacă un împrumut este refuzat/acceptat pe baza mai multor atribute. Regresia logistică este de forma 0/1. y = 0 dacă un împrumut este respins, y = 1 dacă este acceptat.

Un model de regresie logistică diferă de modelul de regresie liniară în două moduri.

• În primul rând, regresia logistică acceptă doar intrare dihotomică (binară) ca variabilă dependentă (adică, un vector de 0 și 1).
• În al doilea rând, rezultatul este reprezentat printr-o funcție de legătură probabilistică numită sigmoid funcție (logistică) datorită formei sale de S:

Ieșirea funcției este întotdeauna între 0 și 1. Verificați imaginea de mai jos

Funcția sigmoidă returnează valori de la 0 la 1. Pentru sarcina de clasificare, avem nevoie de o ieșire discretă de 0 sau 1.

Pentru a converti un flux continuu în valoare discretă, putem seta o limită de decizie la 0.5. Toate valorile peste acest prag sunt clasificate ca 1

Acum, că funcția de legătură este clară, comparați modelul generalizat cu modelul liniar obișnuit pe care îl cunoașteți deja.

GLM vs. Regresie liniară: Diferențe cheie în R

Înainte de a scrie orice cod, este util să știi exact când este necesar standardul regresie liniara Funcția lm() nu mai este potrivită și glm() ar trebui să îi ia locul.

Pe scurt, alegeți lm() când rezultatul este o măsurătoare distribuită normal și glm() ori de câte ori rezultatul este o decizie de tip da/nu, o sarcină pe care altfel ați putea-o încredința unui arborele de decizie clasificator, un număr de evenimente sau o cantitate strict pozitivă a cărei dispersie crește odată cu media sa.

Cum se creează un model liniar generalizat (GLM) în R

Odată stabilită teoria, restul acestui tutorial aplică un GLM binomial de la capăt la capăt pe un set de date real.

Să folosim adult Set de date pentru ilustrarea regresiei logistice. „Adultul” este un set de date excelent pentru sarcina de clasificare. Obiectivul este de a prezice dacă venitul anual al unei persoane în dolari americani va depăși 50,000. Setul de date conține 48,842 de observații și zece variabile:

• vârsta: vârsta individului. Numeric
• educație: Nivelul educațional al individului. Factor.
• starea civilă: Maristatutul individual al individului. Factorul, adică niciodată căsătorit, căsătorit-civ-soț,...
• gen: Genul individului. Factorul, adică bărbat sau femeie
• sursa de venit: Target variabilă. Venituri peste sau sub 50K. Factor, adică >50K, <=50K

printre altii

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

ieșire:

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...

Vom proceda astfel:

• Pasul 1: Verificați variabilele continue
• Pasul 2: Verificați variabilele factor
• Pasul 3: Ingineria caracteristicilor
• Pasul 4: Rezumat statistic
• Pasul 5: Setul de tren/test
• Pasul 6: Construiți modelul
• Pasul 7: Evaluați performanța modelului
• Pasul 8: Îmbunătățiți modelul

Sarcina ta este să prezici care individ va avea un venit mai mare de 50.

În acest tutorial, fiecare pas va fi detaliat pentru a efectua o analiză pe un set de date real.

Pasul 1) Verificați variabilele continue

În primul pas, puteți vedea distribuția variabilelor continue.

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

Code Explicație

• continuu <- select_if(data_adult, is.numeric): Utilizați funcția select_if() din biblioteca dplyr pentru a selecta numai coloanele numerice
• summary(continuous): Imprimați statistica rezumat

ieșire:

##        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	

Din tabelul de mai sus, puteți vedea că datele au scări total diferite și ore.pe.săptămâni au valori aberante mari (adică uitați-vă la ultimul quartil și la valoarea maximă).

Poți să o rezolvi în doi pași:

• Reprezentați grafic distribuția orelor pe săptămână
• Standardizați variabilele continue

1. Trasează distribuția

Să ne uităm mai atent la distribuția orelor.pe.săptămână

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

ieșire:

Variabila are o mulțime de valori aberante și o distribuție nu este bine definită. Puteți rezolva parțial această problemă ștergând primele 0.01% din orele pe săptămână.

Sintaxa de bază a cuantilei:

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.

Calculăm percentila 99 a orelor de lucru săptămânale.

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

Code Explicație

• quantile(data_adult$hours.per.week, .99): Calculează percentila 99 a timpului de lucru săptămânal

ieșire:

## 99% 
##  80

99% din populație lucrează sub 80 de ore pe săptămână.

Puteți renunța la observațiile peste acest prag. Folosești filtrul de la dplyr bibliotecă.

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

ieșire:

## [1] 45537    10

2. Standardizați variabilele continue

Puteți standardiza fiecare coloană pentru a îmbunătăți performanța deoarece datele dvs. nu au aceeași scară. Puteți folosi funcția mutate_if din biblioteca dplyr. Sintaxa de bază este:

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

Puteți standardiza coloanele numerice după cum urmează:

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

Code Explicație

• mutate_if(is.numeric, funs(scale)): condiția este doar o coloană numerică, iar funcția este scară

ieșire:

##           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

Pasul 2) Verificați variabilele factor

Acest pas are două obiective:

• Verificați nivelul din fiecare coloană categorială
• Definiți noi niveluri

Vom împărți acest pas în trei părți:

• Selectați coloanele categorice
• Stocați diagrama cu bare a fiecărei coloane într-o listă
• Tipăriți graficele

Putem selecta coloanele de factori cu codul de mai jos:

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

Code Explicație

• data.frame(select_if(data_adult, is.factor)): Stocăm coloanele factor în factor într-un tip de cadru de date. Biblioteca ggplot2 necesită un obiect cadru de date.

ieșire:

## [1] 6

Setul de date conține 6 variabile categoriale

Al doilea pas este mai solicitant. Doriți să reprezentați grafic o diagramă cu bare pentru fiecare coloană din factorul cadrului de date. Este mai convenabil să automatizați procesul, mai ales în situația în care există multe coloane.

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 Explicație

• lapply(): Utilizați funcția lapply() pentru a trece o funcție în toate coloanele setului de date. Salvați rezultatul într-o listă
• function(x): Funcția va fi procesată pentru fiecare x. Aici x sunt coloanele
• ggplot(factor, aes(get(x))) + geom_bar()+ theme(axis.text.x = element_text(angle = 90)): Creați o diagramă cu bare pentru fiecare element x. Rețineți, pentru a returna x ca coloană, trebuie să-l includeți în get()

Ultimul pas este relativ ușor. Doriți să imprimați cele 6 grafice.

# Print the graph
graph

ieșire:

## [[1]]

## ## [[2]]

## ## [[3]]

## ## [[4]]

## ## [[5]]

## ## [[6]]

Notă: Folosiți butonul următor pentru a naviga la următorul grafic

Pasul 3) Ingineria caracteristicilor

Două variabile categorice au mai multe niveluri decât este necesar pentru model. Le veți regrupa în categorii mai largi, mai bine populate.

Reformarea educației

Din graficul de mai sus, puteți vedea că variabila educație are 16 niveluri. Acest lucru este substanțial, iar unele niveluri au un număr relativ scăzut de observații. Dacă doriți să îmbunătățiți cantitatea de informații pe care o puteți obține din această variabilă, o puteți transforma la un nivel superior. Și anume, creezi grupuri mai mari cu nivel similar de educație. De exemplu, nivelul scăzut de educație va fi transformat în abandon școlar. Nivelurile superioare de educație vor fi schimbate în master.

Iată detaliul:

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 Explicație

• Folosim verbul mutare din biblioteca dplyr. Schimbăm valorile educației cu afirmația ifelse

În tabelul de mai jos, creați o statistică rezumativă pentru a vedea, în medie, câți ani de educație (valoarea z) sunt necesari pentru a ajunge la licență, master sau doctorat.

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

ieșire:

## # 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

Reformare Marital-status

De asemenea, este posibil să se creeze niveluri inferioare pentru starea civilă. În următorul cod modificați nivelul după cum urmează:

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

Puteți verifica numărul de persoane din fiecare grup.

table(recast_data$marital.status)

ieșire:

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

Pasul 4) Statistică rezumată

Este timpul să verificăm câteva statistici despre variabilele noastre țintă. În graficul de mai jos, numărați procentul de persoane care câștigă mai mult de 50, având în vedere sexul lor.

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

ieșire:

Apoi, verificați dacă originea individului îi afectează câștigurile.

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

ieșire:

Numărul de ore de muncă în funcție de sex.

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

ieșire:

Box plot confirmă faptul că distribuția timpului de lucru se potrivește diferitelor grupuri. În box plot, ambele sexe nu au observații omogene.

Puteți verifica densitatea timpului de lucru săptămânal în funcție de tipul de educație. Distribuțiile au multe opțiuni distincte. Probabil se poate explica prin tipul de educație.tracîn SUA.

# 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 Explicație

• ggplot(recast_data, aes(x= hours.per.week)): Un grafic de densitate necesită doar o variabilă
• geom_density(aes(culoare = educație), alpha =0.5): obiectul geometric pentru a controla densitatea

ieșire:

Pentru a vă confirma gândurile, puteți efectua un singur sens test ANOVA:

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

ieșire:

##                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

Testul ANOVA confirmă diferența de medie între grupuri.

Non-liniaritatea

Înainte de a rula modelul, puteți vedea dacă numărul de ore lucrate este legat de vârstă.

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 Explicație

• ggplot(recast_data, aes(x = age, y = hours.per.week)): setați estetica graficului
• geom_point(aes(culoare=venit), dimensiune =0.5): Construiți graficul cu puncte
• stat_smooth(): Adăugați linia de tendință cu următoarele argumente:
  • method='lm': Trasează valoarea ajustată dacă regresie liniara
  • formula = y~poly(x,2): Potriviți o regresie polinomială
  • se = TRUE: Adăugați eroarea standard
  • aes(culoare= venit): sparge modelul după venit

ieșire:

Pe scurt, puteți testa termenii de interacțiune din model pentru a detecta efectul de neliniaritate dintre timpul de lucru săptămânal și alte caracteristici. Este important să se detecteze în ce condiții diferă timpul de lucru.

Corelație

Următoarea verificare este de a vizualiza corelația dintre variabile. Convertiți tipul de nivel de factor în numeric, astfel încât să puteți reprezenta o hartă termică care conține coeficientul de corelație calculat cu metoda Spearman.

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 Explicație

• data.frame(lapply(recast_data,as.integer)): Convertiți datele în numere
• ggcorr() trasează harta termică cu următoarele argumente:
  • metoda: Metodă de calcul a corelației
  • nbreaks = 6: Număr de pauze
  • hjust = 0.8: Poziția de control a numelui variabilei în grafic
  • label = TRUE: Adăugați etichete în centrul ferestrelor
  • label_size = 3: etichete de dimensiune
  • culoare = „grey50”): culoarea etichetei

ieșire:

Pasul 5) Set de tren/test

Orice supravegheat masina de învățare Sarcina necesită să împărțiți datele între un set de trenuri și un set de teste. Puteți utiliza „funcția” pe care ați creat-o în celelalte tutoriale de învățare supravegheată pentru a crea un set de trenuri/teste.

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)

ieșire:

## [1] 36429     9

dim(data_test)

ieșire:

## [1] 9108    9

Pasul 6) Construiți modelul

Pentru a vedea cum funcționează algoritmul, utilizați funcția glm() din pachetul stats. Model liniar generalizat este o colecție de modele. Sintaxa de bază este:

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

Sunteți gata să estimați modelul logistic pentru a împărți nivelul veniturilor între un set de caracteristici.

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

Code Explicație

• formula <- venit ~ .: Creați modelul care să se potrivească
• logit <- glm(formula, data = data_train, family = 'binom'): Potriviți un model logistic (family = 'binom') cu datele data_train.
• summary(logit): Imprimați rezumatul modelului

ieșire:

## 
## 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

Rezumatul modelului nostru dezvăluie informații interesante. Performanța unei regresii logistice este evaluată cu valori cheie specifice.

• AIC (Akaike Information Criteria): Acesta este echivalentul R2 în regresie logistică. Măsoară potrivirea atunci când se aplică o penalizare numărului de parametri. Mai mic AIC valorile indică faptul că modelul este mai aproape de adevăr.
• Devianța nulă: se potrivește modelului numai cu interceptarea. Gradul de libertate este n-1. O putem interpreta ca o valoare Chi-pătrat (valoare ajustată diferită de testarea ipotezei valorii reale).
• Devianța reziduală: Model cu toate variabilele. De asemenea, este interpretată ca o testare a ipotezei Chi-pătrat.
• Numărul de iterații Fisher Scoring: numărul de iterații înainte de convergență.

Ieșirea funcției glm() este stocată într-o listă. Codul de mai jos arată toate elementele disponibile în variabila logit pe care am construit-o pentru a evalua regresia logistică.

# Lista este foarte lungă, imprimați doar primele trei elemente

lapply(logit, class)[1:3]

ieșire:

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

Fiecare valoare poate fi extracmarcat cu semnul $ urmat de numele metricilor. De exemplu, ați stocat modelul ca logit. Pentru a exemplificatracÎn criteriile AIC, utilizați:

logit$aic

ieșire:

## [1] 27086.65

Pasul 7) Evaluați performanța modelului

Matricea confuziei

matrice de confuzie este o alegere mai bună pentru a evalua performanța clasificării în comparație cu diferitele valori pe care le-ați văzut anterior. Ideea generală este de a număra de câte ori instanțele adevărate sunt clasificate ca fiind false.

Pentru a calcula matricea de confuzie, trebuie mai întâi să aveți un set de predicții, astfel încât acestea să poată fi comparate cu țintele reale.

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

Code Explicație

• predict(logit,data_test, type = 'răspuns'): Calculați predicția pe setul de testare. Setați tip = „răspuns” pentru a calcula probabilitatea de răspuns.
• table(data_test$income, predict > 0.5): Calculați matricea de confuzie. prezice > 0.5 înseamnă că returnează 1 dacă probabilitățile prezise sunt peste 0.5, altfel 0.

ieșire:

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

Fiecare rând dintr-o matrice de confuzie reprezintă o țintă reală, în timp ce fiecare coloană reprezintă o țintă prezisă. Primul rând al acestei matrice ia în considerare venitul mai mic de 50 (clasa negativă): 6,310 observații au fost clasificate corect ca indivizi cu venituri mai mici de 50 (Adevărat negativ), în timp ce 495 au fost clasificate în mod eronat ca fiind peste 50 (Fals pozitiv). Al doilea rând ia în considerare veniturile peste 50: 1,229 au fost identificate corect (Adevărat pozitiv), în timp ce 1,074 au fost ratate (Fals negativ).

Puteți calcula modelul precizia prin însumarea adevăratului pozitiv + adevăratul negativ peste observația totală

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

Code Explicație

• sum(diag(table_mat)): Suma diagonalei
• sum(table_mat): Suma matricei.

ieșire:

## [1] 0.8277339

Modelul pare să sufere de o singură problemă: produce prea multe rezultate fals negative. Aceasta se numește paradoxul testului de precizieAm afirmat că acuratețea este raportul dintre predicțiile corecte și numărul total de cazuri. Putem avea o acuratețe relativ mare, dar un model inutil. Se întâmplă atunci când există o clasă dominantă. Dacă ne uităm înapoi la matricea de confuzie, putem vedea că majoritatea cazurilor sunt clasificate drept adevărat negativ. Imaginați-vă acum că modelul a clasificat fiecare observație drept negativă (adică mai mică de 50k). Tot ați obține o acuratețe de aproximativ 75% (6,805 / 9,108). Modelul dumneavoastră are performanțe mai bune, dar se chinuie să distingă adevăratul pozitiv de adevăratul negativ.

Într-o astfel de situație, este de preferat să aveți o metrică mai concisă. Ne putem uita la:

• Precizie=TP/(TP+FP)
• Recall=TP/(TP+FN)

Precizie vs rechemare

Precizie se uită la acuratețea predicției pozitive. Rechemare este raportul instanțelor pozitive care sunt detectate corect de către clasificator;

Puteți construi două funcții pentru a calcula aceste două valori

1. Precizie de construcție

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

Code Explicație

• mat[1,1]: Returnează prima celulă a primei coloane a cadrului de date, adică adevăratul pozitiv
• mat[1,2]; Returnează prima celulă din a doua coloană a cadrului de date, adică fals pozitiv

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

Code Explicație

• mat[1,1]: Returnează prima celulă a primei coloane a cadrului de date, adică adevăratul pozitiv
• mat[2,1]; Returnați a doua celulă a primei coloane a cadrului de date, adică negativul fals

Vă puteți testa funcțiile

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

ieșire:

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

Citiți cu atenție aceste două numere. Precizia este de 0.71, deci atunci când modelul spune că o persoană câștigă peste 50 de dolari, este corect în 71% din cazuri. Rata de rechemare este de 0.53, deci modelul detectează doar 53% dintre persoanele care câștigă cu adevărat peste 50 de dolari.

Puteți crea fișierul scor bazat pe precizie și reamintire. The este o medie armonică a acestor două metrici, ceea ce înseamnă că acordă mai multă pondere valorilor mai mici.

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

ieșire:

## [1] 0.6103799

Precizie vs rechemare compromis

Este imposibil să aveți atât o precizie ridicată, cât și o reamintire ridicată.

Dacă creștem precizia, individul corect va fi mai bine prezis, dar am rata foarte multe dintre ele (reamintire mai mică). În unele situații, preferăm o precizie mai mare decât reamintirea. Există o relație concavă între precizie și reamintire.

• Imaginați-vă, trebuie să preziceți dacă un pacient are o boală. Vrei să fii cât mai precis posibil.
• Dacă trebuie să detectați potențialele persoane frauduloase de pe stradă prin recunoașterea facială, ar fi mai bine să prindeți multe persoane etichetate drept frauduloase, chiar dacă precizia este scăzută. Poliția va putea elibera persoana nefrauduloasă.

Curba ROC

Receptor Operating Caracteristică curba este un alt instrument comun utilizat cu clasificarea binară. Este foarte asemănătoare cu curba de precizie/rechemare, dar în loc de a reprezenta un grafic precizie versus reamintire, curba ROC arată rata pozitivă adevărată (adică, reamintirea) față de rata fals pozitivă. Rata fals pozitive este raportul dintre cazurile negative care sunt clasificate incorect drept pozitive. Este egal cu unu minus rata negativă adevărată. Rata adevărată negativă se mai numește specificitate. De aici curba ROC grafică sensibilitate (reamintire) versus 1-specificitate

Pentru a trasa curba ROC, trebuie să instalăm un pachet numit ROCR. Putem găsi în conda bibliotecă. Puteți introduce codul:

conda install -c r r-rocr --yes

Putem reprezenta ROC cu funcțiile prediction() și 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 Explicație

• prediction(predict, data_test$income): biblioteca ROCR trebuie să creeze un obiect de predicție pentru a transforma datele de intrare
• performanță(ROCRpred, 'tpr','fpr'): returnează cele două combinații de produs în grafic. Aici, tpr și fpr sunt construite. Pentru a reprezenta precizia și a reaminti împreună, folosiți „prec”, „rec”.

ieșire:

Pasul 8) Îmbunătățiți modelul

Puteți încerca să adăugați non-liniaritate modelului cu interacțiunea dintre

• vârsta și ore.pe.săptămână
• sex și ore.pe.săptămână.

Apoi comparați scorul F1 al ambelor modele.

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

ieșire:

## [1] 0.6109181

Scorul F1 este puțin mai mare decât cel precedent. Poți continua să lucrezi la date și să încerci să depășești scorul.

Cum se interpretează coeficienții GLM și rapoartele de probabilitate în R

Tabelul sumar tipărit la Pasul 6 prezintă coeficienții pentru log-cote scară, care este dificil de explicat unui public netehnic. Conversia lor în rapoarte de șanse face ca modelul să fie mult mai ușor de comunicat.

Urmați acești patru pași.

1. Exponențiați coeficienții. Aplicați exp() fiecărei estimări astfel încât logaritmul probabilităților să devină rapoarte de probabilitate multiplicative.
2. Adăugați un interval de încredere. Încadrați confint() în exp() pentru a obține intervalul de 95% pe aceeași scală de probabilitate.
3. Comparați fiecare valoare cu 1. Un raport de șanse mai mare de 1 crește probabilitatea clasei pozitive, o valoare sub 1 o scade, iar o valoare apropiată de 1 înseamnă că predictorul adaugă puțin.
4. Verificați semnificația statistică. Interpretați doar predictorii a căror valoare p în rezultatul rezumat este sub pragul ales, de obicei 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)))

Citirea rezultatului. Coeficientul pentru ore pe săptămână în modelul nostru este 0.41942. Prin exponențiere, rezultă exp(0.41942) = 1.52, ceea ce înseamnă că o creștere cu o deviație standard a numărului de ore lucrate săptămânal multiplică probabilitatea de a câștiga peste 50 de dolari cu aproximativ 1.5, menținând constantă fiecare altă variabilă.

Coeficienții negativi funcționează în același mod. marital.statusNot_married este -2.50346, deci exp(-2.50346) = 0.08: persoanele necăsătorite dețin aproximativ 8% din probabilitatea unei persoane căsătorite. Deoarece predictorii continui au fost standardizați în Pasul 1, descrieți modificările în unități de deviație standard, nu în ore brute.

Notă pentru alte familii: Coeficienții exponențiați sunt rapoarte de probabilitate doar în familia binomială cu o legătură logit. Cu family = „poisson” și o legătură logaritmică, aceleași valori exp() sunt citite în schimb ca rapoarte de rată.

GLM în R: Referință rapidă a funcțiilor

Păstrează acest tabel lângă tine în timp ce programezi. Acesta listează fiecare funcție utilizată în cei opt pași de mai sus, împreună cu pachetul care o furnizează și argumentele pe care le așteaptă.

Celălalt GLM Familiile disponibile prin argumentul familial sunt:

• binom: (legătură = „logit”)
• gaussian: (legătură = „identitate”)
• Gamma: (legătură = „inversă”)
• invers.gaussian: (link = „1/mu^2”)
• poisson: (legătură = „buștean”)
• cvasi: (legătură = „identitate”, varianță = „constantă”)
• cvasibinomial: (legătură = „logit”)
• cvasipoisson: (legătură = „log”)

Întrebări frecvente