GLM u R-u: Generalizirani linearni model i logistička regresija

Što je generalizirani linearni model (GLM) u R-u?

A Generalizirani linearni model (GLM) proširuje običnu linearnu regresiju tako da varijabla odziva može slijediti distribuciju koja nije normalna. U R, ugradite jedan s ugrađenim glm() funkcija iz paketa statistike.

Svaki GLM definiran je s tri komponente:

• Slučajna komponenta: distribucija vjerojatnosti odzivne varijable, uzeta iz eksponencijalne obitelji (binomna, Poissonova, Gama, Gaussova i druge).
• Sustavna komponenta: linearni prediktor, odnosno ponderirana kombinacija vaših objašnjavajućih varijabli.
• Funkcija veze: funkcija koja povezuje srednju vrijednost odgovora s linearnim prediktorom, na primjer logit veza za binarne podatke ili logaritamska veza za brojanje.

Ta struktura čini model „generaliziranim“. Umjesto nametanja normalno distribuiranog ishoda, deklarirate ispravnu distribuciju putem obitelj argument, a R procjenjuje koeficijente metodom maksimalne vjerojatnosti. Logistička regresija je jednostavno GLM s binomnom obitelji i logit vezom, pa je prirodna polazna točka.

Što je logistička regresija u R-u?

Logistička regresija se koristi za predviđanje klase, tj. vjerojatnosti. Logistička regresija može točno predvidjeti binarni ishod.

Zamislite da želite predvidjeti hoće li zajam biti odbijen/prihvaćen na temelju mnogih atributa. Logistička regresija je oblika 0/1. y = 0 ako je posudba odbijena, y = 1 ako je prihvaćena.

Logistički regresijski model razlikuje se od linearnog regresijskog modela na dva načina.

• Prije svega, logistička regresija prihvaća samo dihotomni (binarni) input kao zavisnu varijablu (tj. vektor od 0 i 1).
• Drugo, ishod se mapira putem probabilističke funkcije povezivanja koja se naziva sigmoidni (logistička) funkcija zbog svog S-oblika:

Izlaz funkcije uvijek je između 0 i 1. Provjerite sliku ispod

Sigmoidna funkcija vraća vrijednosti od 0 do 1. Za zadatak klasifikacije potreban nam je diskretni izlaz od 0 ili 1.

Za pretvaranje kontinuiranog protoka u diskretnu vrijednost, možemo postaviti ograničenje odluke na 0.5. Sve vrijednosti iznad ovog praga klasificiraju se kao 1

Sada kada je funkcija poveznice jasna, usporedite generalizirani model s običnim linearnim modelom koji već poznajete.

GLM vs. linearna regresija: Ključne razlike u R-u

Prije nego što napišete bilo kakav kod, korisno je znati točno kada je standard Linearna regresija Funkcija lm() više nije prikladna i glm() bi je trebao zamijeniti.

Ukratko, odaberite lm() kada je ishod normalno distribuirano mjerenje i glm() kada je ishod odluka da/ne, zadatak koji biste inače mogli prepustiti stablo odluke klasifikator, broj događaja ili strogo pozitivna veličina čiji raspon raste sa srednjom vrijednošću.

Kako stvoriti generalizirani linearni model (GLM) u R-u

Nakon što je teorija riješena, ostatak ovog tutorijala primjenjuje binomni GLM od početka do kraja na stvarni skup podataka.

Koristimo se odrasla osoba skup podataka za ilustraciju logističke regresije. „Odrasla osoba“ je izvrstan skup podataka za zadatak klasifikacije. Cilj je predvidjeti hoće li godišnji prihod pojedinca u američkim dolarima premašiti 50,000. Skup podataka sadrži 48,842 opažanja i deset varijabli:

• starost: starost pojedinca. Numerički
• obrazovanje: Obrazovna razina pojedinca. Faktor.
• bračni.status: Maritalni status pojedinca. Čimbenik tj. nikad u braku, u braku-građanski supružnik, …
• spol: Spol pojedinca. Faktor tj. Muško ili Žensko
• prihod: Target varijabla. Prihod iznad ili ispod 50K. Faktor tj. >50K, <=50K

između ostalog

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

Izlaz:

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

Postupit ćemo na sljedeći način:

• Korak 1: Provjerite kontinuirane varijable
• Korak 2: Provjerite faktorske varijable
• Korak 3: Inženjering značajki
• Korak 4: Zbirna statistika
• Korak 5: Skup za obuku/testiranje
• Korak 6: Izgradite model
• Korak 7: Procijenite izvedbu modela
• Korak 8: Poboljšajte model

Vaš zadatak je predvidjeti koja će osoba imati prihod veći od 50K.

U ovom vodiču detaljno će biti opisan svaki korak za izvođenje analize na stvarnom skupu podataka.

Korak 1) Provjerite kontinuirane varijable

U prvom koraku možete vidjeti distribuciju kontinuiranih varijabli.

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

Code Objašnjenje

• kontinuirano <- select_if(data_adult, is.numeric): Koristite funkciju select_if() iz biblioteke dplyr za odabir samo numeričkih stupaca
• sažetak (kontinuirano): Ispis sažete statistike

Izlaz:

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

Iz gornje tablice možete vidjeti da podaci imaju potpuno različite skale i sati.po.tjednima imaju velike odstupanja (tj. pogledajte zadnji kvartil i maksimalnu vrijednost).

S tim se možete nositi u dva koraka:

• Nacrtajte raspodjelu sati tjedno
• Standardizirajte kontinuirane varijable

1. Nacrtajte distribuciju

Pogledajmo pobliže distribuciju sati po tjednu

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

Izlaz:

Varijabla ima puno outliera i nije dobro definiranu distribuciju. Ovaj problem možete djelomično riješiti brisanjem gornjih 0.01 posto sati tjedno.

Osnovna sintaksa kvantila:

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.

Izračunavamo 99-ti percentil tjednog radnog vremena.

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

Code Objašnjenje

• kvantil(data_adult$hours.per.week, .99): Izračunaj 99. percentil tjednog radnog vremena

Izlaz:

## 99% 
##  80

99 posto stanovništva radi manje od 80 sati tjedno.

Možete ispustiti opažanja iznad ovog praga. Koristite filtar iz dplyr knjižnica.

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

Izlaz:

## [1] 45537    10

2. Standardizirajte kontinuirane varijable

Možete standardizirati svaki stupac kako biste poboljšali izvedbu jer vaši podaci nemaju istu ljestvicu. Možete koristiti funkciju mutate_if iz biblioteke dplyr. Osnovna sintaksa je:

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

Brojčane stupce možete standardizirati na sljedeći način:

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

Code Objašnjenje

• mutate_if(is.numeric, funs(scale)): uvjet je samo numerički stupac, a funkcija je scale

Izlaz:

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

Korak 2) Provjerite faktorske varijable

Ovaj korak ima dva cilja:

• Provjerite razinu u svakom stupcu kategorije
• Definirajte nove razine

Ovaj korak ćemo podijeliti u tri dijela:

• Odaberite kategoričke stupce
• Pohranite trakasti grafikon svakog stupca na popis
• Ispišite grafikone

Stupce faktora možemo odabrati pomoću donjeg koda:

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

Code Objašnjenje

• data.frame(select_if(data_adult, is.factor)): Stupce faktora pohranjujemo u faktor u tipu okvira podataka. Biblioteka ggplot2 zahtijeva objekt okvira podataka.

Izlaz:

## [1] 6

Skup podataka sadrži 6 kategoričkih varijabli

Drugi korak je zahtjevniji. Želite nacrtati stupčasti grafikon za svaki stupac u faktoru okvira podataka. Praktičnije je automatizirati proces, posebno u situacijama kada postoji mnogo stupaca.

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 Objašnjenje

• lapply(): Koristite funkciju lapply() za prosljeđivanje funkcije u svim stupcima skupa podataka. Izlaz pohranjujete na popis
• funkcija(x): funkcija će biti obrađena za svaki x. Ovdje je x stupac
• ggplot(faktor, aes(get(x))) + geom_bar()+ theme(axis.text.x = element_text(angle = 90)): Napravite trakasti dijagram za svaki x element. Napomena, da vratite x kao stupac, morate ga uključiti unutar get()

Posljednji korak je relativno lak. Želite ispisati 6 grafikona.

# Print the graph
graph

Izlaz:

## [[1]]

## ## [[2]]

## ## [[3]]

## ## [[4]]

## ## [[5]]

## ## [[6]]

Napomena: Koristite sljedeći gumb za navigaciju do sljedećeg grafikona

Korak 3) Inženjering značajki

Dvije kategoričke varijable nose više razina nego što je modelu potrebno. Pregrupirat ćete ih u šire, bolje naseljene kategorije.

Preinačiti obrazovanje

Iz gornjeg grafikona možete vidjeti da varijabla obrazovanje ima 16 razina. To je značajno, a neke razine imaju relativno mali broj opažanja. Ako želite poboljšati količinu informacija koju možete dobiti iz ove varijable, možete je preinačiti na višu razinu. Naime, stvarate veće grupe sa sličnim stupnjem obrazovanja. Na primjer, niska razina obrazovanja pretvorit će se u napuštanje škole. Više razine obrazovanja bit će promijenjene u master.

Evo detalja:

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 Objašnjenje

• Koristimo glagol mutirati iz biblioteke dplyr. Izjavom ifelse mijenjamo vrijednosti obrazovanja

U donjoj tablici kreirate sažetak statistike da vidite, u prosjeku, koliko je godina obrazovanja (z-vrijednost) potrebno da se postigne diploma, magisterij ili doktorat.

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

Izlaz:

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

Preoblikovati Marital-status

Također je moguće stvoriti niže razine za bračni status. U sljedećem kodu mijenjate razinu na sljedeći način:

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

Možete provjeriti broj pojedinaca unutar svake grupe.

table(recast_data$marital.status)

Izlaz:

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

Korak 4) Sažeta statistika

Vrijeme je da provjerimo neke statistike o našim ciljnim varijablama. Na donjem grafikonu brojite postotak pojedinaca koji zarađuju više od 50 tisuća s obzirom na njihov spol.

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

Izlaz:

Zatim provjerite utječe li porijeklo pojedinca na njihovu zaradu.

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

Izlaz:

Broj sati rada prema spolu.

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

Izlaz:

Okvirni dijagram potvrđuje da raspodjela radnog vremena odgovara različitim skupinama. U okvirnom dijagramu, oba spola nemaju homogena opažanja.

Možete provjeriti gustoću tjednog radnog vremena prema vrsti obrazovanja. Raspodjela ima mnogo različitih odabira. To se vjerojatno može objasniti vrstom obrazovanjatract u SAD-u.

# 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 Objašnjenje

• ggplot(recast_data, aes( x= hours.per.week)): Dijagram gustoće zahtijeva samo jednu varijablu
• geom_density(aes(color = education), alpha =0.5): Geometrijski objekt za kontrolu gustoće

Izlaz:

Kako biste potvrdili svoje misli, možete izvesti jednosmjerni ANOVA test:

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

Izlaz:

##                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 test potvrđuje razliku u prosjeku između skupina.

Nelinearnost

Prije nego što pokrenete model, možete vidjeti je li broj radnih sati povezan s dobi.

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 Objašnjenje

• ggplot(recast_data, aes(x = age, y = sati.po.tjednu)): Postavite estetiku grafikona
• geom_point(aes(color= prihod), veličina =0.5): Konstruirajte točkasti dijagram
• stat_smooth(): Dodajte liniju trenda sa sljedećim argumentima:
  • method='lm': iscrtajte prilagođenu vrijednost ako je Linearna regresija
  • formula = y~poly(x,2): Uklopi polinomsku regresiju
  • se = TRUE: Dodajte standardnu ​​pogrešku
  • aes(boja= prihod): Razdvojite model prema prihodu

Izlaz:

Ukratko, možete testirati uvjete interakcije u modelu kako biste uočili učinak nelinearnosti između tjednog radnog vremena i drugih značajki. Važno je otkriti pod kojim se uvjetima radno vrijeme razlikuje.

Korelacija

Sljedeća provjera je vizualizacija korelacije između varijabli. Vrstu razine faktora pretvarate u numeričku tako da možete iscrtati toplinsku kartu koja sadrži koeficijent korelacije izračunat Spearmanovom metodom.

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 Objašnjenje

• data.frame(lapply(recast_data,as.integer)): Pretvori podatke u numeričke
• ggcorr() crta toplinsku kartu sa sljedećim argumentima:
  • metoda: Metoda za izračunavanje korelacije
  • nbreaks = 6: Broj prekida
  • hjust = 0.8: Kontrolni položaj naziva varijable u dijagramu
  • label = TRUE: Dodajte oznake u središte prozora
  • label_size = 3: Oznake veličine
  • boja = “siva50”): Boja naljepnice

Izlaz:

Korak 5) Skup za obuku/testiranje

Svaki pod nadzorom stroj za učenje Zadatak zahtijeva da podijelite podatke između skupa vlakova i skupa testova. Možete koristiti "funkciju" koju ste stvorili u drugim tutorijalima za nadzirano učenje za stvaranje skupa vlakova/testova.

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)

Izlaz:

## [1] 36429     9

dim(data_test)

Izlaz:

## [1] 9108    9

Korak 6) Izgradite model

Da biste vidjeli kako algoritam radi, koristite funkciju glm() iz paketa stats. Generalizirani linearni model je kolekcija modela. Osnovna sintaksa je:

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

Spremni ste procijeniti logistički model kako biste podijelili razinu prihoda između skupa značajki.

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

Code Objašnjenje

• formula <- prihod ~ .: Stvorite model koji odgovara
• logit <- glm(formula, data = data_train, family = 'binomial'): Uklopi logistički model (family = 'binomial') s podacima data_train.
• summary(logit): Ispis sažetka modela

Izlaz:

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

Sažetak našeg modela otkriva zanimljive podatke. Učinkovitost logističke regresije procjenjuje se pomoću specifičnih ključnih metrika.

• AIC (Akaike informacijski kriteriji): Ovo je ekvivalent R2 u logističkoj regresiji. Mjeri prikladnost kada se kazna primijeni na broj parametara. Manji AIC vrijednosti pokazuju da je model bliži istini.
• Nulta devijacija: odgovara modelu samo s presjekom. Stupanj slobode je n-1. Možemo je protumačiti kao vrijednost hi-kvadrat (prilagođena vrijednost koja se razlikuje od testiranja hipoteze stvarne vrijednosti).
• Rezidualna devijacija: Model sa svim varijablama. Također se tumači kao testiranje hi-kvadrat hipoteze.
• Broj ponavljanja Fisherova bodovanja: Broj ponavljanja prije konvergiranja.

Izlaz funkcije glm() pohranjuje se na popis. Kod u nastavku prikazuje sve stavke dostupne u logit varijabli koju smo konstruirali za procjenu logističke regresije.

# Popis je jako dug, ispišite samo prva tri elementa

lapply(logit, class)[1:3]

Izlaz:

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

Svaka vrijednost može biti npr.tracoznačeno znakom $ nakon kojeg slijedi naziv metrike. Na primjer, model ste pohranili kao logit. Na primjertracPrema AIC kriterijima, koristite:

logit$aic

Izlaz:

## [1] 27086.65

Korak 7) Procijenite izvedbu modela

Matrica zabune

The matrica zabune je bolji izbor za procjenu izvedbe klasifikacije u usporedbi s različitim mjernim podacima koje ste vidjeli prije. Opća ideja je brojati koliko su puta True instance klasificirane kao False.

Da biste izračunali matricu zabune, prvo morate imati skup predviđanja kako bi se mogla usporediti sa stvarnim ciljevima.

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

Code Objašnjenje

• predict(logit,data_test, type = 'response'): Izračunajte predviđanje na skupu testova. Postavite type = 'response' za izračunavanje vjerojatnosti odgovora.
• table(data_test$income, predict > 0.5): Izračunajte matricu zabune. predict > 0.5 znači da vraća 1 ako su predviđene vjerojatnosti iznad 0.5, inače 0.

Izlaz:

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

Svaki redak u matrici konfuzije predstavlja stvarni cilj, dok svaki stupac predstavlja predviđeni cilj. Prvi redak ove matrice uzima u obzir prihod niži od 50 tisuća dolara (negativna klasa): 6,310 opažanja ispravno je klasificirano kao pojedinci s prihodom nižim od 50 tisuća dolara (Istinski negativan), dok je 495 pogrešno klasificirano kao iznad 50 tisuća (Lažno pozitivno). Drugi redak uzima u obzir prihod iznad 50 tisuća: 1,229 je ispravno identificirano (Prava pozitiva), dok je 1,074 propušteno (Lažno negativan).

Možete izračunati model točnost zbrajanjem stvarnog pozitivnog + istinskog negativnog ukupnog promatranja

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

Code Objašnjenje

• sum(diag(table_mat)): Zbroj dijagonale
• sum(table_mat): Zbroj matrice.

Izlaz:

## [1] 0.8277339

Čini se da model pati od jednog problema: proizvodi previše lažno negativnih rezultata. To se naziva paradoks testa točnostiNaveli smo da je točnost omjer točnih predviđanja i ukupnog broja slučajeva. Možemo imati relativno visoku točnost, ali beskoristan model. To se događa kada postoji dominantna klasa. Ako pogledate matricu zbunjenosti, možete vidjeti da je većina slučajeva klasificirana kao istinski negativna. Zamislite sada da je model klasificirao svako opažanje kao negativno (tj. ispod 50 tisuća). I dalje biste dobili točnost od oko 75 posto (6,805 / 9,108). Vaš model radi bolje, ali se muči razlikovati istinski pozitivno od istinski negativnog.

U takvoj situaciji poželjno je imati sažetiju metriku. Možemo pogledati:

• Preciznost=TP/(TP+FP)
• Opoziv=TP/(TP+FN)

Preciznost protiv prisjećanja

Preciznost gleda na točnost pozitivnog predviđanja. Podsjetiti je omjer pozitivnih instanci koje je klasifikator ispravno otkrio;

Možete konstruirati dvije funkcije za izračunavanje ove dvije metrike

1. Konstruirajte preciznost

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

Code Objašnjenje

• mat[1,1]: vraća prvu ćeliju prvog stupca podatkovnog okvira, tj. pravi pozitivan
• prostirka[1,2]; Vrati prvu ćeliju drugog stupca podatkovnog okvira, tj. lažno pozitivno

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

Code Objašnjenje

• mat[1,1]: vraća prvu ćeliju prvog stupca podatkovnog okvira, tj. pravi pozitivan
• prostirka[2,1]; Vrati drugu ćeliju prvog stupca podatkovnog okvira, tj. lažno negativno

Možete testirati svoje funkcije

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

Izlaz:

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

Pažljivo pročitajte ova dva broja. Preciznost je 0.71, pa kada model kaže da pojedinac zarađuje više od 50 tisuća, to je točno u 71 posto slučajeva. Podsjetnik je 0.53, pa model detektira samo 53 posto pojedinaca koji stvarno zarađuju više od 50 tisuća.

Možete stvoriti rezultat na temelju preciznosti i prisjećanja. The je harmonijska sredina ove dvije metrike, što znači da daje veću težinu nižim vrijednostima.

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

Izlaz:

## [1] 0.6103799

Kompromis preciznosti u odnosu na opoziv

Nemoguće je imati i visoku preciznost i visoko pamćenje.

Ako povećamo preciznost, točna će se osoba bolje predvidjeti, ali bismo propustili puno njih (slabije pamćenje). U nekim situacijama više volimo veću preciznost od opoziva. Postoji konkavan odnos između preciznosti i prisjećanja.

• Zamislite, morate predvidjeti ima li pacijent bolest. Želite biti što precizniji.
• Ako trebate otkriti potencijalne prevarante na ulici putem prepoznavanja lica, bilo bi bolje uhvatiti mnogo ljudi označenih kao prevaranti iako je preciznost niska. Policija će moći pustiti osobu koja nije prijevaru.

ROC krivulja

The Prijamnik Operating Karakterističan krivulja je još jedan uobičajeni alat koji se koristi s binarnom klasifikacijom. Vrlo je slična krivulji preciznost/prisjećanje, ali umjesto iscrtavanja preciznosti u odnosu na prisjećanje, ROC krivulja pokazuje pravu pozitivnu stopu (tj. prisjećanje) naspram lažno pozitivne stope. Lažno pozitivna stopa je omjer negativnih slučajeva koji su netočno klasificirani kao pozitivni. Jednaka je jedan minus prava negativna stopa. Prava negativna stopa također se naziva specifičnost. Stoga je prikazana ROC krivulja osjetljivost (prisjećanje) u odnosu na 1-specifičnost

Za crtanje ROC krivulje, potrebno je instalirati paket pod nazivom ROCR. Možemo ga pronaći u conda knjižnica. Možete upisati kod:

conda install -c r r-rocr --yes

Možemo iscrtati ROC s funkcijama 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 Objašnjenje

• predviđanje(predviđanje, data_test$income): ROCR biblioteka treba stvoriti objekt predviđanja za transformaciju ulaznih podataka
• performance(ROCRpred, 'tpr','fpr'): Vrati dvije kombinacije za izradu na grafikonu. Ovdje su konstruirani tpr i fpr. Za precizno iscrtavanje i opoziv zajedno koristite "prec", "rec".

Izlaz:

Korak 8) Poboljšajte model

Možete pokušati dodati nelinearnost modelu s interakcijom između

• dob i sati.po.tjednu
• spol i sati.po.tjednu.

Zatim uspoređujete F1 rezultat oba modela.

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

Izlaz:

## [1] 0.6109181

Rezultat F1 je nešto viši od prethodnog. Možete nastaviti raditi na podacima i pokušati nadmašiti rezultat.

Kako interpretirati GLM koeficijente i omjere šansi u R-u

Sažeta tablica ispisana u koraku 6 prikazuje koeficijente na logaritamske vjerojatnosti skala, što je teško objasniti netehničkoj publici. Pretvaranjem u omjere šansi model je puno lakši za komunikaciju.

Slijedite ova četiri koraka.

1. Potenciraj koeficijente. Primijenite exp() na svaku procjenu tako da logaritamske vjerojatnosti postanu multiplikativni omjeri vjerojatnosti.
2. Dodajte interval pouzdanosti. Omotajte confint() u exp() kako biste dobili interval od 95 posto na istoj skali vjerojatnosti.
3. Usporedite svaku vrijednost s 1. Omjer šansi iznad 1 povećava vjerojatnost pozitivne klase, vrijednost ispod 1 je smanjuje, a vrijednost blizu 1 znači da prediktor malo dodaje.
4. Provjerite statističku značajnost. Interpretirajte samo prediktore čija je p-vrijednost u sažetku ispod odabranog praga, obično 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)))

Čitanje izlaza. Koeficijent za sate tjedno u našem modelu iznosi 0.41942. Njegovo potenciranje daje exp(0.41942) = 1.52, što znači da povećanje tjednog radnog vremena za jednu standardnu ​​devijaciju množi vjerojatnost zarade iznad 50 tisuća dolara s otprilike 1.5, držeći svaku drugu varijablu konstantnom.

Negativni koeficijenti funkcioniraju na isti način. marital.statusNot_married je -2.50346, pa je exp(-2.50346) = 0.08: neoženjene osobe imaju oko 8 posto vjerojatnosti za oženjenu osobu. Budući da su kontinuirani prediktori standardizirani u 1. koraku, opišite promjene u jedinicama standardne devijacije, a ne u sirovim satima.

Napomena za ostale obitelji: Eksponencijalni koeficijenti su omjeri šansi samo pod binomnom obitelji s logit vezom. S obitelji = "poisson" i log vezom, iste vrijednosti exp() čitaju se kao omjeri stopa.

GLM u R-u: Kratki pregled funkcija

Držite ovu tablicu pokraj sebe dok pišete kod. U njoj je navedena svaka funkcija korištena u gornjih osam koraka, zajedno s paketom koji je pruža i argumentima koje očekuje.

Drugi GLM obitelji dostupne putem obiteljskog argumenta su:

• binom: (link = “logit”)
• gaussova: (veza = “identitet”)
• Gama: (link = “inverzno”)
• inverzna.gaussova: (link = “1/mu^2”)
• poisson: (link = “log”)
• kvazi: (veza = „identitet“, varijanca = „konstanta“)
• kvazibinom: (link = “logit”)
• kvazipoisson: (link = “log”)

Pitanja i odgovori