GLM R-ben: Általánosított lineáris modell és logisztikus regresszió

Mi az az általánosított lineáris modell (GLM) R-ben?

A Általánosított lineáris modell (GLM) kiterjeszti a közönséges lineáris regressziót, így a válaszváltozó a normálistól eltérő eloszlást követhet. R, beépített glm() függvény a statisztikacsomagból.

Minden GLM-et három összetevő határoz meg:

• Véletlenszerű komponens: a válaszváltozó valószínűségi eloszlása, amelyet az exponenciális családból (binomiális, Poisson-, Gamma-, Gauss-eloszlás és mások) veszünk fel.
• Szisztematikus komponens: a lineáris prediktor, azaz a magyarázó változók súlyozott kombinációja.
• Link funkció: az a függvény, amely a válasz átlagát a lineáris prediktorhoz köti, például a bináris adatok logit kapcsolata vagy a darabszámok logaritmikus kapcsolata.

Ez a struktúra teszi a modellt „általánosítottá”. Ahelyett, hogy normális eloszlású eredményt erőltetnénk, a helyes eloszlást a következőn keresztül deklaráljuk: család argumentum, és R a maximális valószínűséggel becsüli meg az együtthatókat. A logisztikus regresszió egyszerűen egy GLM binomiális családdal és logit kapcsolattal, így ez a természetes kiindulópont.

Mi a logisztikus regresszió R-ben?

A logisztikus regressziót egy osztály, azaz egy valószínűség előrejelzésére használják. A logisztikus regresszió pontosan megjósolhatja a bináris eredményt.

Képzelje el, hogy meg akarja jósolni, hogy a hitel elutasítása/elfogadása számos tulajdonság alapján történik. A logisztikus regresszió 0/1 formájú. y = 0, ha egy kölcsönt elutasítanak, y = 1, ha elfogadják.

A logisztikus regressziós modell két dologban különbözik a lineáris regressziós modelltől.

• Először is, a logisztikus regresszió csak a dichotóm (bináris) bemenetet fogadja el függő változóként (azaz 0 és 1 vektort).
• Másodszor, az eredményt egy valószínűségi kapcsolódási függvény, az úgynevezett szigma alakú (logisztikai) függvény S alakja miatt:

A függvény kimenete mindig 0 és 1 között van. Ellenőrizze az alábbi képet

A szigmoid függvény 0 és 1 közötti értékeket ad vissza. Az osztályozási feladathoz 0 vagy 1 diszkrét kimenetre van szükségünk.

A folytonos áramlás diszkrét értékké alakításához beállíthatunk egy döntési korlátot 0.5-re. Minden e küszöbérték feletti érték 1-esnek minősül

Most, hogy a linkfüggvény világos, hasonlítsd össze az általánosított modellt a már ismert közönséges lineáris modellel.

GLM vs. lineáris regresszió: Főbb különbségek az R-ben

Mielőtt bármilyen kódot írnál, hasznos pontosan tudni, hogy mikor érvényes a szabvány lineáris regresszió Az lm() függvény már nem megfelelő, a helyét a glm() függvénynek kell átvennie.

Röviden, válaszd az lm()-t, ha az eredmény normális eloszlású mérés, és a glm()-t, ha az eredmény igen/nem döntés – egy olyan feladatot, amelyet egyébként egy másik felhasználónak is át lehetne adni. döntési fa osztályozó, egy eseményszámláló, vagy egy szigorúan pozitív mennyiség, amelynek szórása az átlagával együtt növekszik.

Általánosított lineáris modell (GLM) létrehozása R-ben

Miután az elmélet tisztázódott, a tutoriál további része egy binomiális GLM-et alkalmaz végponttól végpontig egy valós adathalmazon.

Használjuk a felnőtt egy adathalmaz a logisztikus regresszió szemléltetésére. A „felnőtt” adathalmaz nagyszerű az osztályozási feladathoz. A cél annak előrejelzése, hogy egy személy éves jövedelme meghaladja-e az 50 000 dollárt amerikai dollárban. Az adathalmaz 48 842 megfigyelést és tíz változót tartalmaz:

• életkor: az egyén életkora. Numerikus
• oktatás: Az egyén iskolai végzettsége. Tényező.
• családi állapot: Mariaz egyén státusza. Tényező pl. Soha nem házas, Házas-civ-házastárs,…
• gender: Az egyén neme. Tényező, azaz Férfi vagy Nő
• jövedelem: Target változó. 50 ezer feletti vagy alatti bevétel. Tényező pl. >50K, <=50K

többek között

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

output:

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

A következőképpen járunk el:

• 1. lépés: Ellenőrizze a folytonos változókat
• 2. lépés: Ellenőrizze a faktorváltozókat
• 3. lépés: Funkciótervezés
• 4. lépés: Összefoglaló statisztika
• 5. lépés: Képzés/tesztkészlet
• 6. lépés: Készítse el a modellt
• 7. lépés: Mérje fel a modell teljesítményét
• 8. lépés: Javítsa a modellt

Az Ön feladata, hogy megjósolja, melyik személynek lesz nagyobb bevétele 50 XNUMX-nál.

Ebben az oktatóanyagban minden lépést részletesen ismertetünk egy valós adatkészlet elemzéséhez.

1. lépés) Ellenőrizze a folytonos változókat

Első lépésben a folytonos változók eloszlását láthatjuk.

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

Code Magyarázat

• folyamatos <- select_if(data_adult, is.numeric): A select_if() függvény segítségével csak a numerikus oszlopokat jelölje ki a dplyr könyvtárból
• összegzés(folyamatos): Az összesítő statisztika nyomtatása

output:

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

A fenti táblázatból láthatja, hogy az adatok teljesen eltérő skálákkal rendelkeznek, és az órákonkénti órákban nagy kiugró értékek vannak (.azaz nézd meg az utolsó kvartilist és a maximális értéket).

Két lépéssel kezelheti:

• Ábrázolja a heti órák eloszlását
• Szabványosítsa a folytonos változókat

1. Ábrázolja az eloszlást

Nézzük meg közelebbről az órák.hetenkénti eloszlását

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

output:

A változónak sok kiugró értéke van, és nem jól definiált eloszlása. Ezt a problémát részben megoldhatja, ha törli a heti munkaórák felső 0.01 százalékát.

A kvantilis alapvető szintaxisa:

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.

Kiszámítjuk a heti munkaórák 99. percentilisét.

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

Code Magyarázat

• quantile(data_adult$hours.per.week, .99): Kiszámítja a heti munkaidő 99. percentilisét

output:

## 99% 
##  80

A lakosság 99 százaléka heti 80 óra alatt dolgozik.

A megfigyeléseket e küszöb fölé dobhatja. Használja a szűrőt a dplyr könyvtár.

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

output:

## [1] 45537    10

2. Szabványosítsa a folytonos változókat

A teljesítmény javítása érdekében minden oszlopot szabványosíthat, mivel az adatok nem azonos léptékűek. Használhatja a mute_if függvényt a dplyr könyvtárból. Az alap szintaxis a következő:

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

A numerikus oszlopokat az alábbiak szerint szabványosíthatja:

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

Code Magyarázat

• mute_if(is.numeric, funs(scale)): A feltétel csak numerikus oszlop, a függvény pedig skála

output:

##           X         age        workclass    education educational.num
## 1 -1.732680 -1.02325949          Private         11th     -1.22106443
## 2 -1.732605 -0.03969284          Private      HS-grad     -0.43998868
## 3 -1.732530 -0.79628257        Local-gov   Assoc-acdm      0.73162494
## 4 -1.732455  0.41426100          Private Some-college     -0.04945081
## 5 -1.732379 -0.34232873          Private         10th     -1.61160231
## 6 -1.732304  1.85178149 Self-emp-not-inc  Prof-school      1.90323857
##       marital.status  race gender hours.per.week income
## 1      Never-married Black   Male    -0.03995944  <=50K
## 2 Married-civ-spouse White   Male     0.86863037  <=50K
## 3 Married-civ-spouse White   Male    -0.03995944   >50K
## 4 Married-civ-spouse Black   Male    -0.03995944   >50K
## 5      Never-married White   Male    -0.94854924  <=50K
## 6 Married-civ-spouse White   Male    -0.76683128   >50K

2. lépés) Ellenőrizze a faktorváltozókat

Ennek a lépésnek két célja van:

• Ellenőrizze a szintet minden kategóriás oszlopban
• Határozzon meg új szinteket

Ezt a lépést három részre osztjuk:

• Válassza ki a kategorikus oszlopokat
• Tárolja az egyes oszlopok oszlopdiagramját egy listában
• Nyomtassa ki a grafikonokat

A faktor oszlopokat az alábbi kóddal tudjuk kiválasztani:

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

Code Magyarázat

• data.frame(select_if(data_adult, is.factor)): A faktor oszlopokat faktorban tároljuk egy adatkeret típusban. A ggplot2 könyvtár adatkeret objektumot igényel.

output:

## [1] 6

Az adatkészlet 6 kategorikus változót tartalmaz

A második lépés igényesebb. Minden oszlophoz oszlopdiagramot kell készíteni az adatkeret-faktorban. Kényelmesebb automatizálni a folyamatot, különösen sok oszlop esetén.

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 Magyarázat

• lapply(): Használja a lapply() függvényt egy függvény átadásához az adatkészlet összes oszlopában. A kimenetet egy listában tárolja
• függvény(x): A függvény minden x esetén feldolgozásra kerül. Itt x az oszlopok
• ggplot(faktor, aes(get(x))) + geom_bar()+ theme(axis.text.x = element_text(angle = 90)): Hozzon létre egy oszlopdiagramot minden x elemhez. Megjegyzés: az x oszlopként való visszaadásához fel kell vennie a get() függvénybe.

Az utolsó lépés viszonylag egyszerű. Ki szeretné nyomtatni a 6 grafikont.

# Print the graph
graph

output:

## [[1]]

## ## [[2]]

## ## [[3]]

## ## [[4]]

## ## [[5]]

## ## [[6]]

Megjegyzés: A következő gombbal navigálhat a következő grafikonra

3. lépés) Funkciótervezés

Két kategorikus változó több szintet hordoz, mint amennyire a modellnek szüksége van. Ezeket szélesebb, jobban kitöltött kategóriákba kell csoportosítani.

Átdolgozott oktatás

A fenti grafikonon látható, hogy a változó végzettségnek 16 szintje van. Ez jelentős, és egyes szinteken viszonylag alacsony a megfigyelések száma. Ha javítani szeretné az ebből a változóból nyerhető információ mennyiségét, akkor átdolgozhatja azt magasabb szintre. Ugyanis nagyobb, hasonló iskolai végzettségű csoportokat hoz létre. Például az alacsony iskolai végzettséget lemorzsolják. A felsőfokú oktatás mesterképzésre változik.

Íme a részlet:

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 Magyarázat

• A mute igét használjuk a dplyr könyvtárból. Az ifelse állítással megváltoztatjuk az oktatás értékrendjét

Az alábbi táblázatban összefoglaló statisztikát készít, amely megmutatja, hogy átlagosan hány évnyi oktatás (z-érték) szükséges a Bachelor, Master vagy PhD fokozat megszerzéséhez.

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

output:

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

Átdolgozás Marital-státusz

A családi állapothoz alacsonyabb szinteket is lehet létrehozni. A következő kódban a következőképpen módosíthatja a szintet:

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

Az egyes csoportokon belül ellenőrizheti a személyek számát.

table(recast_data$marital.status)

output:

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

4. lépés) Összefoglaló statisztika

Itt az ideje, hogy ellenőrizzünk néhány statisztikai adatot a célváltozóinkról. Az alábbi grafikonon megszámolja azon személyek százalékos arányát, akik nemük alapján 50 ezernél többet keresnek.

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

output:

Ezután ellenőrizze, hogy az egyén származása befolyásolja-e a keresetét.

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

output:

A munkaórák száma nemek szerint.

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

output:

A dobozdiagram megerősíti, hogy a munkaidő elosztása különböző csoportokhoz illeszkedik. A dobozdiagramban mindkét nemnél nincsenek homogén megfigyelések.

A heti munkaidő sűrűségét ellenőrizheti az oktatás típusa szerint. Az eloszlások sok különböző elemet mutatnak. Ez valószínűleg a képzés típusával magyarázható.tracaz USA-ban.

# 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 Magyarázat

• ggplot(recast_data, aes(x= hours.per.week)): A sűrűségdiagram csak egy változót igényel
• geom_density(aes(szín = oktatás), alfa =0.5): A sűrűséget vezérlő geometriai objektum

output:

Gondolatainak megerősítésére egyirányú eljárást hajthat végre ANOVA teszt:

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

output:

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

Az ANOVA teszt megerősíti a csoportok közötti átlagkülönbséget.

Nemlinearitás

A modell futtatása előtt megnézheti, hogy a ledolgozott órák száma összefüggésben van-e az életkorral.

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 Magyarázat

• ggplot(recast_data, aes(x = életkor, y = óra.hetenként)): A grafikon esztétikájának beállítása
• geom_point(aes(szín=jövedelem), méret =0.5): Készítse el a pontdiagramot
• stat_smooth(): Adja hozzá a trendvonalat a következő argumentumokkal:
  • method='lm': Ábrázolja az illesztett értéket, ha a lineáris regresszió
  • formula = y~poly(x,2): Illesszen egy polinomiális regressziót
  • se = TRUE: Adja hozzá a standard hibát
  • aes(szín=jövedelem): Bontsa fel a modellt jövedelem szerint

output:

Dióhéjban, tesztelheti az interakciós kifejezéseket a modellben, hogy felismerje a heti munkaidő és más jellemzők közötti nemlinearitási hatást. Fontos felismerni, hogy a munkaidő milyen körülmények között tér el.

Összefüggés

A következő ellenőrzés a változók közötti korreláció megjelenítése. A faktorszint-típust numerikussá alakítja, így elkészítheti a Spearman-módszerrel kiszámított korrelációs együtthatót tartalmazó hőtérképet.

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 Magyarázat

• data.frame(lapply(recast_data,as.integer)): Adatok konvertálása numerikussá
• ggcorr() ábrázolja a hőtérképet a következő argumentumokkal:
  • módszer: A korreláció kiszámításának módszere
  • nbreaks = 6: A szünetek száma
  • hjust = 0.8: A változó nevének vezérlési pozíciója a diagramban
  • label = TRUE: Adjon hozzá címkéket az ablakok közepére
  • label_size = 3: Méretcímkék
  • szín = „grey50”): a címke színe

output:

5. lépés) Képzés/tesztkészlet

Bármilyen felügyelt gépi tanulás A feladat megköveteli az adatok felosztását egy vonatkészlet és egy tesztkészlet között. A többi felügyelt tanulási oktatóanyagban létrehozott „függvényt” használhatod egy vonat-/tesztkészlet létrehozásához.

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)

output:

## [1] 36429     9

dim(data_test)

output:

## [1] 9108    9

6. lépés) Építse meg a modellt

Az algoritmus teljesítményének megtekintéséhez használd a stats csomag glm() függvényét. Általánosított lineáris modell modellek gyűjteménye. Az alap szintaxis a következő:

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

Készen áll a logisztikai modell becslésére, hogy a bevételi szintet feloszthassa egy sor jellemző között.

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

Code Magyarázat

• formula <- jövedelem ~ .: Készítse el a modellt, hogy illeszkedjen
• logit <- glm(képlet, adatok = data_train, family = 'binomial'): Illesszen egy logisztikai modellt (family = 'binomial') a data_train adatokkal.
• summary(logit): Kinyomtatja a modell összefoglalóját

output:

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

Modellünk összefoglalója érdekes információkat tár fel. A logisztikus regresszió teljesítményét meghatározott kulcsmutatókkal értékelik.

• AIC (Akaike Information Criteria): Ez a megfelelője R2 logisztikus regresszióban. Méri az illeszkedést, amikor büntetést alkalmaznak a paraméterek számára. Kisebb AIC Az értékek azt mutatják, hogy a modell közelebb áll az igazsághoz.
• Null deviance: Csak a metszőponttal illeszkedik a modellhez. A szabadság foka n-1. Értelmezhetjük Khi-négyzet értékként (a tényleges érték hipotézis tesztelésétől eltérő illesztett érték).
• Maradék eltérés: Modell az összes változóval. Khi-négyzet hipotézis tesztelésként is értelmezhető.
• Fisher Scoring iterációk száma: A konvergálás előtti iterációk száma.

A glm() függvény kimenete egy listában tárolódik. Az alábbi kód a logisztikus regresszió értékelésére létrehozott logit változóban elérhető összes elemet mutatja.

# A lista nagyon hosszú, csak az első három elemet nyomtassa ki

lapply(logit, class)[1:3]

output:

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

Minden érték lehet extracA $ jellel jelölt metrikák neve után. Például a modellt logit néven tárolta. PéldáultracAz AIC kritériumok alapján a következőket használja:

logit$aic

output:

## [1] 27086.65

7. lépés) Értékelje a modell teljesítményét

Zavart mátrix

Az zavart mátrix jobb választás az osztályozási teljesítmény értékelésére a korábban látott különböző mérőszámokhoz képest. Az általános ötlet az, hogy megszámoljuk, hogy az Igaz példányok hányszor hamisnak minősülnek.

A zavaros mátrix kiszámításához először rendelkeznie kell egy előrejelzéskészlettel, hogy össze lehessen hasonlítani őket a tényleges célokkal.

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

Code Magyarázat

• előrejelzés(logit,data_teszt, type = 'válasz'): Kiszámítja az előrejelzést a tesztkészleten. Állítsa be a type = 'response' értéket a válasz valószínűségének kiszámításához.
• táblázat(adat_teszt$bevétel, előrejelzés > 0.5): Számítsa ki a zavaros mátrixot. előrejelzés > 0.5 azt jelenti, hogy 1-et ad vissza, ha az előrejelzett valószínűségek 0.5 felett vannak, különben 0.

output:

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

A zavaros mátrix minden sora egy tényleges célt jelöl, míg minden oszlop egy előre jelzett célt. A mátrix első sora az 50 ezer dollárnál alacsonyabb jövedelműeket veszi figyelembe (a negatív osztály): 6,310 megfigyelést helyesen soroltunk be 50 ezer dollárnál alacsonyabb jövedelműként (Igaz negatív), míg 495-öt tévesen 50 ezer feletti kategóriába soroltak (Álpozitív). A második sor az 50 ezer feletti jövedelmet veszi figyelembe: 1,229-et azonosítottak helyesen (Igaz pozitív), míg 1,074-at kihagytak (Hamis negatív).

Kiszámolhatja a modellt pontosság úgy, hogy az igazi pozitív + igaz negatív összeget összeadja a teljes megfigyeléssel

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

Code Magyarázat

• sum(diag(table_mat)): Az átló összege
• sum(table_mat): A mátrix összege.

output:

## [1] 0.8277339

Úgy tűnik, a modellnek egyetlen problémája van: túl sok álnegatív eredményt produkál. Ezt nevezik pontossági vizsgálat paradoxonaAzt állítottuk, hogy a pontosság a helyes előrejelzések és az esetek teljes számának aránya. Viszonylag magas pontosságot érhetünk el, de a modellünk haszontalan. Ez akkor fordul elő, ha van egy domináns osztály. Ha visszanézzük a zavaró mátrixot, láthatjuk, hogy az esetek többségét valódi negatívként osztályozzák. Képzeljük el most, hogy a modell minden megfigyelést negatívként osztályoz (azaz 50 ezer alatt). Még mindig körülbelül 75 százalékos pontosságot kapnánk (6,805 / 9,108). A modellünk jobban teljesít, de nehezen tudja megkülönböztetni a valódi pozitívat a valódi negatívtól.

Ilyen helyzetben célszerű tömörebb mérőszámot használni. Megnézhetjük:

• Pontosság = TP/(TP+FP)
• Visszahívás=TP/(TP+FN)

Precision vs Recall

Pontosság a pozitív előrejelzés pontosságát vizsgálja. visszahívás az osztályozó által helyesen észlelt pozitív esetek aránya;

A két metrika kiszámításához két függvényt hozhat létre

1. Konstrukciós precizitás

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

Code Magyarázat

• mat[1,1]: Visszaadja az adatkeret első oszlopának első celláját, azaz az igazi pozitívat
• mat[1,2]; Adja vissza az adatkeret második oszlopának első celláját, azaz hamis pozitívat

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

Code Magyarázat

• mat[1,1]: Visszaadja az adatkeret első oszlopának első celláját, azaz az igazi pozitívat
• mat[2,1]; Adja vissza az adatkeret első oszlopának második celláját, azaz hamis negatívot

Kipróbálhatja funkcióit

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

output:

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

Olvassa el figyelmesen ezt a két számot. A pontosság 0.71, tehát amikor a modell azt mondja, hogy egy személy 50 ezer dollár felett keres, az az esetek 71 százalékában helyes. A felidézhetőség 0.53, tehát a modell csak az 53 százalékát érzékeli azoknak az egyéneknek, akik valóban 50 ezer dollár felett keresnek.

Létrehozhatod a pontozás a pontosság és a felidézés alapján. A e két metrika harmonikus átlaga, vagyis nagyobb súlyt ad az alacsonyabb értékeknek.

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

output:

## [1] 0.6103799

Precision vs Recall kompromisszum

Lehetetlen egyszerre nagy pontosság és nagy visszahívás.

Ha növeljük a pontosságot, jobban megjósolható lesz a megfelelő egyed, de sok közülük hiányozna (alacsonyabb visszahívás). Bizonyos helyzetekben jobban szeretjük a nagyobb pontosságot, mint a visszahívást. Konkáv kapcsolat van a pontosság és a felidézés között.

• Képzelje el, meg kell jósolnia, ha a betegnek betegsége van. A lehető legpontosabb akar lenni.
• Ha arcfelismeréssel kell felderítenie a potenciális csaló embereket az utcán, jobb lenne sok olyan személyt elkapni, akiket csalónak minősítenek, annak ellenére, hogy a pontosság alacsony. A rendőrség szabadon engedheti a nem csaló személyt.

A ROC görbe

Az Receiver Operating Jellemző A görbe egy másik gyakori eszköz a bináris osztályozáshoz. Nagyon hasonlít a precíziós/visszahívási görbéhez, de ahelyett, hogy a pontosságot a felidézéssel ábrázolná, a ROC görbe a valódi pozitív arányt (azaz a visszahívást) mutatja a hamis pozitív aránnyal szemben. A hamis pozitív arány a hibásan pozitívnak minősített negatív esetek aránya. Ez egyenlő eggyel mínusz a valódi negatív ráta. A valódi negatív rátát is nevezik sajátosság. Ezért a ROC görbe ábrázolja érzékenység (visszahívás) versus 1-specifikusság

A ROC görbe ábrázolásához telepítenünk kell egy ROCR nevű csomagot. A conda fájlban megtaláljuk könyvtár. Beírhatja a kódot:

conda install -c r r-rocr --yes

A ROC-t a predikció() és a performance() függvényekkel ábrázolhatjuk.

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 Magyarázat

• predikció(predict, data_teszt$bevétel): Az ROCR könyvtárnak létre kell hoznia egy predikciós objektumot a bemeneti adatok átalakításához
• performance (ROCRpred, 'tpr', 'fpr'): Visszaadja a két kombinációt a grafikonban előállítani. Itt tpr és fpr épül fel. A ábrázolás pontosságának és felidézésének együttes használata, használja a „prec”, „rec” kifejezéseket.

output:

8. lépés) Javítsa a modellt

Megpróbálhat nemlinearitást hozzáadni a modellhez a közötti interakcióval

• életkor és óra.hetente
• nem és óra.per.hét.

Ezután összehasonlítod a két modell F1 pontszámát.

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

output:

## [1] 0.6109181

Az F1 pontszám valamivel magasabb, mint az előző. Tovább dolgozhatsz az adatokon, és megpróbálhatod megdönteni a pontszámot.

Hogyan értelmezzük a GLM együtthatókat és az esélyhányadosokat az R-ben?

A 6. lépésben kinyomtatott összefoglaló táblázat a következő együtthatókat mutatja: log-odds skála, amit nehéz elmagyarázni egy nem műszaki közönség számára. Az esélyhányadosokká való átalakításuk sokkal könnyebbé teszi a modell kommunikálását.

Kövesd ezt a négy lépést.

1. Hatványozd az együtthatókat. Alkalmazzuk az exp() függvényt minden becslésre, hogy a log-oddsok multiplikatív esélyhányadosokká váljanak.
2. Adjon hozzá egy konfidenciaintervallumot. A confint() függvényt az exp() függvénybe csomagolva kapjuk meg a 95 százalékos intervallumot ugyanazon az esélyskálán.
3. Hasonlítsa össze az egyes értékeket az 1-gyel. Az 1 feletti esélyhányados növeli a pozitív osztály valószínűségét, az 1 alatti érték csökkenti, az 1-hez közeli érték pedig azt jelenti, hogy a prediktor keveset ad hozzá.
4. Ellenőrizze a statisztikai szignifikanciát. Csak olyan prediktorokat értelmezzen, amelyeknek az összegző kimenetben szereplő p-értéke a kiválasztott küszöbérték, jellemzően 0.05 alatt van.

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

A kimenet olvasása. A modellünkben a heti órák számának együtthatója 0.41942. Ha ezt hatványozzuk, akkor exp(0.41942) = 1.52-t kapunk, ami azt jelenti, hogy a heti munkaórák egy szórásnyi növekedése nagyjából 1.5-szeresére növeli az 50 ezer feletti kereset esélyét, minden más változót állandónak tekintve.

A negatív együtthatók ugyanígy működnek. A martial.statusNot_married értéke -2.50346, tehát exp(-2.50346) = 0.08: a nem házas egyének esélyei körülbelül 8 százalékkal magasabbak, mint a házas egyéneké. Mivel a folytonos prediktorokat az 1. lépésben standardizáltuk, a változásokat szórásegységekben, ne nyers órákban írjuk le.

Megjegyzés más családoknak: A hatványozott együtthatók csak a binomiális család alatti esélyhányadosok logit kapcsolattal. Ha a család = „poisson” és logaritmikus kapcsolat van, akkor ugyanazokat az exp() értékeket olvassuk be arányhányadosként.

GLM R-ben: Gyors függvényreferencia

Tartsd magad mellett ezt a táblázatot kódolás közben. Felsorolja a fenti nyolc lépésben használt összes függvényt, a függvényt szolgáltató csomaggal és a várt argumentumokkal együtt.

A másik GLM A családi érvelésen keresztül elérhető családok a következők:

• binomiális: (link = „logit”)
• gaussian: (link = „identitás”)
• Gamma: (link = „inverz”)
• inverse.gaussian: (link = „1/mu^2”)
• poisson: (link = „napló”)
• kvázi: (kapcsolat = „azonosság”, variancia = „konstans”)
• kvázibinomiális: (link = „logit”)
• kvázipoisson: (link = „napló”)

GYIK