GLM u R: Generalizirani linearni model s primjerom

ล to je logistiฤka regresija?

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 mjeri sljedeฤ‡om funkcijom vjerojatnosne veze tzv sigmoidni zbog svog oblika slova S.:

Logistiฤka regresija

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

Logistiฤka regresija

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

Logistiฤka regresija

Kako stvoriti generalizirani model linije (GLM)

Koristimo se odrasla osoba skup podataka za ilustraciju logistiฤke regresije. "Odrasli" je izvrstan skup podataka za zadatak klasifikacije. Cilj je predvidjeti hoฤ‡e li godiลกnji prihod u dolarima pojedinca premaลกiti 50.000. Skup podataka sadrลพi 46,033 XNUMX opaลพanja i deset znaฤajki:

  • 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:

  • 1: Nacrtajte distribuciju sati po tjednu
  • 2: 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:

Provjerite kontinuirane varijable

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 gornja 2 posto percentila

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

Code Objaลกnjenje

  • kvantil(data_adult$hours.per.week, .99): Izraฤunajte vrijednost 99 posto radnog vremena

Izlaz:

## 99% 
##  80

98 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
  1. 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 vjeลกtiji. ลฝelite iscrtati trakasti grafikon za svaki stupac u faktoru okvira podataka. Pogodnije je automatizirati proces, posebno u situaciji 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]]

Provjerite varijable faktora

## ## [[2]]

Provjerite varijable faktora

## ## [[3]]

Provjerite varijable faktora

## ## [[4]]

Provjerite varijable faktora

## ## [[5]]

Provjerite varijable faktora

## ## [[6]]

Provjerite varijable faktora

Napomena: Koristite sljedeฤ‡i gumb za navigaciju do sljedeฤ‡eg grafikona

Provjerite varijable faktora

Korak 3) Inลพenjering znaฤajki

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:

Stara razina Nova razina
Predลกkolski maloljetnik koji se ispisao iz ลกkole
10th maloljetnik koji se ispisao iz ลกkole
11th maloljetnik koji se ispisao iz ลกkole
12th maloljetnik koji se ispisao iz ลกkole
1.-4 maloljetnik koji se ispisao iz ลกkole
5th-6th maloljetnik koji se ispisao iz ลกkole
7th-8th maloljetnik koji se ispisao iz ลกkole
9th maloljetnik koji se ispisao iz ลกkole
HS-Grad HighGrad
Neki fakulteti Zajednica
izv.-acdm Zajednica
izv. voc Zajednica
Celibat Celibat
Masters Masters
Prof-ลกkola Masters
Doktorat Dr.
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:

Stara razina Nova razina
Nikad oลพenjen Neoลพenjen
Oลพenjen-supruลพnik-odsutan Neoลพenjen
Oลพenjen-AF-braฤni drug Oลพenjen
Oลพenjen-graฤ‘anka-supruลพnik
Odvojen Odvojen
Rastavljen
udovice Udovica
# 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:

Sumarna statistika

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:

Sumarna statistika

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:

Sumarna statistika

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:

Sumarna statistika

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:

Nelinearnost

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:

Korelacija

Korak 5) Skup za obuku/testiranje

Svaki pod nadzorom stroj za uฤenje Zadatak zahtijeva dijeljenje podataka izmeฤ‘u skupa vlakova i testnog skupa. Moลพete upotrijebiti "funkciju" koju ste izradili u drugim vodiฤima za nadzirano uฤenje za izradu skupa za treniranje/testiranje.

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 paket glm(). The 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.

Matrica zabune

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 zabune predstavlja stvarni cilj, dok svaki stupac predstavlja predviฤ‘eni cilj. Prvi redak ove matrice uzima u obzir dohodak niลพi od 50 tisuฤ‡a (klasa laลพnih): 6241 je ispravno klasificirano kao pojedinci s dohotkom niลพim od 50 tisuฤ‡a (Istinski negativan), dok je preostali pogreลกno klasificiran kao iznad 50 tisuฤ‡a (Laลพno pozitivno). U drugom redu su primanja iznad 50k, pozitivna klasa je 1229 (Prava pozitiva), dok Istinski negativan bio 1074.

Moลพete izraฤunati model toฤnost zbrajanjem stvarnog pozitivnog + istinskog negativnog ukupnog promatranja

Matrica zabune

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, precjenjuje broj laลพno negativnih rezultata. Ovo se zove paradoks testa toฤnosti. Naveli smo da je toฤnost omjer toฤnih predviฤ‘anja prema ukupnom broju sluฤajeva. Moลพemo imati relativno visoku toฤnost, ali beskoristan model. To se dogaฤ‘a kada postoji dominantna klasa. Ako pogledate unatrag na matricu zabune, moลพete vidjeti da je veฤ‡ina sluฤajeva klasificirana kao stvarno negativna. Zamislite sada, model je klasificirao sve klase kao negativne (tj. niลพe od 50k). Imali biste toฤnost od 75 posto (6718/6718+2257). Vaลก model ima bolju izvedbu, ali muฤi se razlikovati pravo pozitivno od pravog 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

Kada model kaลพe da se radi o osobi iznad 50 tisuฤ‡a, toฤan je u samo 54 posto sluฤajeva, a moลพe potraลพivati โ€‹โ€‹osobe iznad 50 tisuฤ‡a u 72 posto sluฤajeva.

Moลพete stvoriti Preciznost protiv prisjeฤ‡anja rezultat na temelju preciznosti i prisjeฤ‡anja. The Preciznost protiv prisjeฤ‡anja je harmonijska sredina ove dvije metrike, ลกto znaฤi da daje veฤ‡u teลพinu niลพim vrijednostima.

Preciznost protiv prisjeฤ‡anja

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 iscrtavanje ROC krivulje moramo instalirati biblioteku koja se zove RORC. Moลพemo pronaฤ‡i u kondi knjiลพnica. Moลพete upisati kod:

conda install -cr r-rocr โ€“da

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:

ROC krivulja

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.

Morate upotrijebiti test rezultata da biste usporedili 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

Ocjena je neลกto viลกa od prethodne. Moลพete nastaviti raditi na podacima pokuลกavajuฤ‡i nadmaลกiti rezultat.

Rezime

Moลพemo saลพeti funkciju za treniranje logistiฤke regresije u tablici u nastavku:

Paket Cilj funkcija Argument
- Stvorite skup podataka o vlaku/testiranju create_train_set() podaci, veliฤina, vlak
glm Uvjeลพbajte generalizirani linearni model glm() formula, podaci, obitelj*
glm Saลพmite model Saลพetak() ugraฤ‘eni model
baza Napravite predviฤ‘anje predvidjeti() prilagoฤ‘eni model, skup podataka, tip = 'odgovor'
baza Napravite matricu zabune stol() y, predvidi()
baza Stvorite ocjenu toฤnosti zbroj(dijag(tablica())/zbroj(tablica()
ROCR Stvorite ROC: Korak 1 Stvorite predviฤ‘anje predviฤ‘anje() predvidjeti(), y
ROCR Stvorite ROC: 2. korak Stvorite izvedbu izvoฤ‘enje() predviฤ‘anje(), 'tpr', 'fpr'
ROCR Stvorite ROC: korak 3 Iscrtajte grafikon zemljiลกte() izvoฤ‘enje()

Drugi GLM vrste modela su:

โ€“ binom: (link = โ€œlogitโ€)

โ€“ gaussov: (link = โ€œidentitetโ€)

โ€“ Gama: (link = โ€œobrnutoโ€)

โ€“ inverse.gaussian: (link = โ€œ1/mu^2โ€)

โ€“ poisson: (link = โ€œlogโ€)

โ€“ kvazi: (link = โ€œidentitetโ€, varijanca = โ€œkonstantaโ€)

โ€“ kvazibinom: (link = โ€œlogitโ€)

โ€“ quasipoisson: (link = โ€œlogโ€)

Saลพmite ovu objavu uz: