GLM R:ssä: Yleistetty lineaarinen malli ja logistinen regressio

Mikä on yleistetty lineaarinen malli (GLM) R:ssä?

A Yleistetty lineaarinen malli (GLM) laajentaa tavallista lineaarista regressiota siten, että vastemuuttuja voi noudattaa muuta kuin normaalijakaumaa. R, asennat sellaisen, jossa on sisäänrakennettu glm() funktio tilastopaketista.

Jokainen GLM määritellään kolmen komponentin perusteella:

• Satunnainen komponentti: vastemuuttujan todennäköisyysjakauma, joka on otettu eksponentiaalisesta perheestä (binomi-, Poisson-, Gamma-, Gauss- ja muut).
• Systemaattinen komponentti: lineaarinen ennustaja eli selittävien muuttujien painotettu yhdistelmä.
• Linkkitoiminto: funktio, joka yhdistää vasteen keskiarvon lineaariseen ennustajaan, esimerkiksi binääridatan logit-linkki tai lukumäärien log-linkki.

Tuo rakenne tekee mallista "yleistettyä". Normaalijakauman pakottamisen sijaan oikea jakauma määritetään käyttämällä perhe argumentti ja R arvioi kertoimet suurimman uskottavuuden perusteella. Logistinen regressio on yksinkertaisesti GLM, jossa on binomiparveke ja logit-linkki, joten se on luonnollinen lähtökohta.

Mikä on logistinen regressio R:ssä?

Logistista regressiota käytetään ennustamaan luokkaa eli todennäköisyyttä. Logistinen regressio voi ennustaa binäärituloksen tarkasti.

Kuvittele, että haluat ennustaa, onko laina evätty/hyväksytty monien ominaisuuksien perusteella. Logistinen regressio on muotoa 0/1. y = 0, jos laina hylätään, y = 1, jos laina hyväksytään.

Logistinen regressiomalli eroaa lineaarisesta regressiomallista kahdella tavalla.

• Ensinnäkin logistinen regressio hyväksyy vain dikotomisen (binäärisyötteen) riippuvaiseksi muuttujaksi (eli vektorin 0 ja 1).
• Toiseksi, lopputulos kuvataan todennäköisyyspohjaisen linkkifunktion, jota kutsutaan sigmoidi (logistinen) funktio S-muotonsa vuoksi:

Toiminnon lähtö on aina välillä 0 ja 1. Tarkista alla oleva kuva

Sigmoidifunktio palauttaa arvot välillä 0–1. Luokittelutehtävää varten tarvitsemme erillisen tulosteen 0 tai 1.

Jatkuvan virtauksen muuttamiseksi diskreetiksi arvoksi voimme asettaa päätösrajaksi 0.5. Kaikki tämän kynnyksen ylittävät arvot luokitellaan 1:ksi

Nyt kun linkkifunktio on selvä, vertaa yleistettyä mallia jo tuntemaasi tavalliseen lineaariseen malliin.

GLM vs. lineaarinen regressio: R:n keskeiset erot

Ennen kuin kirjoitat koodia, on hyödyllistä tietää tarkalleen, milloin standardi lineaarinen regressio Funktio lm() ei ole enää sopiva ja sen tilalle tulisi ottaa glm().

Lyhyesti sanottuna, valitse lm(), kun tulos on normaalijakautunut mittaustulos, ja glm(), kun tulos on kyllä/ei-päätös, tehtävä, jonka muuten voisit antaa jollekin toiselle. päätöspuu luokittelija, tapahtumien määrä tai ehdottoman positiivinen suure, jonka hajonta kasvaa keskiarvon mukana.

Yleistetyn lineaarisen mallin (GLM) luominen R:ssä

Teorian selvitettyä loppuosa tästä tutoriaalista soveltaa binomista GLM:ää päästä päähän todelliseen tietojoukkoon.

Käytetään aikuinen aineisto logistisen regression havainnollistamiseksi. ”Aikuinen” on erinomainen aineisto luokittelutehtävään. Tavoitteena on ennustaa, ylittävätkö yksilön vuositulot Yhdysvaltain dollareissa 50 000. Aineisto sisältää 48 842 havaintoa ja kymmenen muuttujaa:

• ikä: henkilön ikä. Numeerinen
• koulutus: Yksilön koulutustaso. Tekijä.
• siviilisääty: Marihenkilön asemasta. Factor eli ei koskaan naimisissa, Naimisissa-siviili-puoliso,…
• sukupuoli: yksilön sukupuoli. Tekijä, eli mies tai nainen
• tulo: Target muuttuja. Tulot yli tai alle 50 50. Kerroin eli >50K, <=XNUMXK

muiden joukossa

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

lähtö:

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

Jatketaan seuraavasti:

• Vaihe 1: Tarkista jatkuvat muuttujat
• Vaihe 2: Tarkista tekijämuuttujat
• Vaihe 3: Ominaisuussuunnittelu
• Vaihe 4: Yhteenvetotilasto
• Vaihe 5: Harjoittelu/testisarja
• Vaihe 6: Rakenna malli
• Vaihe 7: Arvioi mallin suorituskyky
• Vaihe 8: Paranna mallia

Sinun tehtäväsi on ennustaa, kenen tulot ovat yli 50 XNUMX.

Tässä opetusohjelmassa jokainen vaihe on kuvattu yksityiskohtaisesti todellisen tietojoukon analysoimiseksi.

Vaihe 1) Tarkista jatkuvat muuttujat

Ensimmäisessä vaiheessa näet jatkuvien muuttujien jakauman.

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

Code Selitys

• jatkuva <- select_if(data_adult, is.numeric): Käytä funktiota select_if() dplyr-kirjastosta valitaksesi vain numeeriset sarakkeet
• yhteenveto (jatkuva): Tulosta yhteenvetotilasto

lähtö:

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

Yllä olevasta taulukosta näet, että tiedoilla on täysin erilaiset asteikot ja tunnit.viikkoa kohden on suuria poikkeavuuksia (.eli katso viimeistä kvartiilia ja maksimiarvoa).

Voit käsitellä sen seuraavasti:

• Piirrä viikkotuntien jakauma
• Standardoi jatkuvat muuttujat

1. Piirrä jakelu

Katsotaanpa tarkemmin tuntien.per.viikko-jakaumaa

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

lähtö:

Muuttujalla on paljon poikkeavuuksia eikä tarkka jakautuminen. Voit ratkaista tämän ongelman osittain poistamalla ylimmät 0.01 prosenttia viikon tunneista.

Kvantiilin perussyntaksi:

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.

Laskemme viikoittaisten työtuntien 99. persentiilin.

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

Code Selitys

• quantile(data_adult$hours.per.week, .99): Laskee viikoittaisen työajan 99. persentiilin

lähtö:

## 99% 
##  80

99 prosenttia väestöstä työskentelee alle 80 tuntia viikossa.

Voit pudottaa havainnot tämän kynnyksen yläpuolelle. Käytät suodatinta dplyr kirjasto.

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

lähtö:

## [1] 45537    10

2. Standardoi jatkuvat muuttujat

Voit standardoida jokaisen sarakkeen tehokkuuden parantamiseksi, koska tiedoillasi ei ole samaa mittakaavaa. Voit käyttää funktiota mutate_if dplyr-kirjastosta. Perussyntaksi on:

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

Voit standardoida numeeriset sarakkeet seuraavasti:

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

Code Selitys

• mutate_if(is.numeric, funs(scale)): Ehto on vain numeerinen sarake ja funktio on skaala

lähtö:

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

Vaihe 2) Tarkista tekijämuuttujat

Tällä vaiheella on kaksi tavoitetta:

• Tarkista taso kussakin kategoriassa sarakkeessa
• Määrittele uusia tasoja

Jaamme tämän vaiheen kolmeen osaan:

• Valitse kategorialliset sarakkeet
• Tallenna kunkin sarakkeen pylväskaavio luetteloon
• Tulosta kaaviot

Voimme valita tekijäsarakkeet alla olevalla koodilla:

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

Code Selitys

• data.frame(select_if(data_adult, is.factor)): Tallennamme tekijäsarakkeet kertoimena tietokehystyyppiin. Kirjasto ggplot2 vaatii tietokehysobjektin.

lähtö:

## [1] 6

Tietojoukko sisältää 6 kategorista muuttujaa

Toinen vaihe on vaativampi. Haluat piirtää pylväsdiagrammin jokaiselle datakehyksen tekijän sarakkeelle. Prosessin automatisointi on helpompaa, varsinkin tilanteissa, joissa sarakkeita on paljon.

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 Selitys

• lapply(): Käytä funktiota lapply() välittääksesi funktion kaikissa tietojoukon sarakkeissa. Tallennat tulosteen luetteloon
• function(x): Funktio käsitellään jokaiselle x:lle. Tässä x on sarakkeet
• ggplot(tekijä, aes(get(x))) + geom_bar()+ teema(akseli.teksti.x = elementin_teksti(kulma = 90)): Luo pylväsmerkkikaavio jokaiselle x-elementille. Huomaa, että jos haluat palauttaa x:n sarakkeena, sinun on sisällytettävä se get()

Viimeinen vaihe on suhteellisen helppo. Haluat tulostaa 6 kuvaajaa.

# Print the graph
graph

lähtö:

## [[1]]

## ## [[2]]

## ## [[3]]

## ## [[4]]

## ## [[5]]

## ## [[6]]

Huomautus: Käytä seuraava-painiketta siirtyäksesi seuraavaan kaavioon

Vaihe 3) Ominaisuussuunnittelu

Kahdella kategorisella muuttujalla on enemmän tasoja kuin malli tarvitsee. Ryhmittelet ne uudelleen laajempiin, paremmin täytettyihin kategorioihin.

Uudelleenlaadittu koulutus

Yllä olevasta kaaviosta näet, että muuttuvassa koulutuksessa on 16 tasoa. Tämä on merkittävää, ja joillakin tasoilla on suhteellisen vähän havaintoja. Jos haluat parantaa tästä muuttujasta saatavan tiedon määrää, voit muotoilla sen uudelleen korkeammalle tasolle. Nimittäin luot suurempia ryhmiä, joilla on samanlainen koulutustaso. Esimerkiksi alhainen koulutus muuttuu keskeyttäneiksi. Korkeammat koulutustasot muutetaan mestariksi.

Tässä on yksityiskohta:

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 Selitys

• Käytämme verbiä mutate dplyr-kirjastosta. Muutamme koulutuksen arvoja lauseella ifelse

Alla olevaan taulukkoon luot yhteenvetotilaston, josta näet, kuinka monta vuotta koulutusta (z-arvo) keskimäärin kestää kandidaatin, maisterin tai tohtorin tutkinnon saavuttamiseen.

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

lähtö:

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

uudelleenlaaditun Marital-status

On myös mahdollista luoda alempia tasoja siviilisäätyä varten. Seuraavassa koodissa muutat tasoa seuraavasti:

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

Voit tarkistaa kunkin ryhmän yksilöiden lukumäärän.

table(recast_data$marital.status)

lähtö:

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

Vaihe 4) Yhteenvetotilasto

On aika tarkistaa tilastotietoja tavoitemuuttujistamme. Alla olevassa kaaviossa lasket niiden henkilöiden prosenttiosuuden, jotka ansaitsevat yli 50 XNUMX sukupuolensa perusteella.

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

lähtö:

Tarkista seuraavaksi, vaikuttaako henkilön alkuperä hänen tuloihinsa.

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

lähtö:

Työtuntien määrä sukupuolen mukaan.

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

lähtö:

Laatikkokaavio vahvistaa, että työajanjako sopii eri ryhmiin. Laatikkokaaviossa molemmilla sukupuolilla ei ole homogeenisiä havaintoja.

Voit tarkistaa viikoittaisen työajan tiheyden koulutustyypin mukaan. Jakaumissa on useita erillisiä valintoja. Se voidaan luultavasti selittää koulutustyypin mukaan.tract Yhdysvalloissa.

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

• ggplot(recast_data, aes(x= tuntia.per.viikko)): Tiheyskaavio vaatii vain yhden muuttujan
• geom_density(aes(väri = koulutus), alfa =0.5): geometrinen objekti tiheyden ohjaamiseen

lähtö:

Vahvistaaksesi ajatuksesi voit suorittaa yksisuuntaisen ANOVA testi:

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

lähtö:

##                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-testi vahvistaa ryhmien välisen keskiarvon eron.

Ei-lineaarisuus

Ennen kuin suoritat mallin, voit nähdä, liittyykö työtuntien määrä ikään.

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 Selitys

• ggplot(recast_data, aes(x = ikä, y = tuntia.per.viikko)): Aseta kaavion estetiikka
• geom_point(aes(väri=tulo), koko =0.5): Muodosta pistekaavio
• stat_smooth(): Lisää trendiviiva seuraavilla argumenteilla:
  • method='lm': Piirrä sovitettu arvo, jos lineaarinen regressio
  • kaava = y~poly(x,2): Sovita polynomiregressio
  • se = TOSI: Lisää vakiovirhe
  • aes(väri=tulot): Jaa malli tulojen mukaan

lähtö:

Lyhyesti sanottuna, voit testata vuorovaikutustermejä mallissa saadaksesi selville viikoittaisen työajan ja muiden ominaisuuksien epälineaarisuuden vaikutuksen. On tärkeää havaita, missä olosuhteissa työaika vaihtelee.

Korrelaatio

Seuraava tarkistus on visualisoida muuttujien välinen korrelaatio. Muunnat tekijätason tyypin numeeriseksi, jotta voit piirtää lämpökartan, joka sisältää Spearman-menetelmällä lasketun korrelaatiokertoimen.

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 Selitys

• data.frame(lapply(recast_data,as.integer)): Muunna tiedot numeerisiksi
• ggcorr() piirtää lämpökartan seuraavilla argumenteilla:
  • method: Menetelmä korrelaation laskemiseksi
  • nbreaks = 6: Taukojen lukumäärä
  • hjust = 0.8: Säätää muuttujan nimen sijaintia kaaviossa
  • etiketti = TOSI: Lisää tunnisteet ikkunoiden keskelle
  • label_size = 3: Kokotarrat
  • väri = "grey50"): Tarran väri

lähtö:

Vaihe 5) Harjoittelu/testisarja

Kaikki valvotut koneoppiminen Tehtävä vaatii datan jakamista junajoukon ja testijoukon kesken. Voit käyttää muissa ohjatun oppimisen tutoriaaleissa luomaasi "funktiota" juna-/testijoukon luomiseen.

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)

lähtö:

## [1] 36429     9

dim(data_test)

lähtö:

## [1] 9108    9

Vaihe 6) Rakenna malli

Algoritmin toiminnan näkemiseksi käytät stats-paketin glm()-funktiota. Yleistetty lineaarinen malli on kokoelma malleja. Perussyntaksi on:

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

Olet valmis arvioimaan logistisen mallin tulotason jakamiseksi eri ominaisuuksien kesken.

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

Code Selitys

• kaava <- tulo ~ .: Luo malli sopivaksi
• logit <- glm(kaava, data = data_train, family = 'binomial'): Sovita logistinen malli (family = 'binomial') data_train-tietoihin.
• summary(logit): Tulosta mallin yhteenveto

lähtö:

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

Mallimme yhteenveto paljastaa mielenkiintoisia tietoja. Logistisen regression suorituskykyä arvioidaan tietyillä keskeisillä mittareilla.

• AIC (Akaike Information Criteria): Tämä vastaa R2 logistisessa regressiossa. Se mittaa sovituksen, kun parametrien määrään sovelletaan rangaistusta. Pienempi AIC arvot osoittavat, että malli on lähempänä totuutta.
• Nollapoikkeama: Sopii malliin vain leikkauspisteen kanssa. Vapausaste on n-1. Voimme tulkita sen Chi-neliöarvona (sovitettu arvo, joka eroaa todellisen arvon hypoteesin testauksesta).
• Jäännöspoikkeama: Malli kaikilla muuttujilla. Se tulkitaan myös khin neliön hypoteesin testaamiseksi.
• Fisher Scoring iteraatioiden määrä: Iteraatioiden määrä ennen konvergoimista.

Glm()-funktion tulos tallennetaan luetteloon. Alla oleva koodi näyttää kaikki saatavilla olevat kohteet logit-muuttujassa, jonka rakensimme arvioimaan logistista regressiota.

# Lista on erittäin pitkä, tulosta vain kolme ensimmäistä elementtiä

lapply(logit, class)[1:3]

lähtö:

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

Jokainen arvo voi olla esim.tracja lisää $-merkki ja sen perään mittareiden nimet. Esimerkiksi tallensit mallin nimellä logit. Esim.tracAIC-kriteereissä käytät:

logit$aic

lähtö:

## [1] 27086.65

Vaihe 7) Arvioi mallin suorituskyky

Sekaannusmatriisi

sekaannusmatriisi on parempi valinta arvioida luokituksen tehokkuutta verrattuna eri mittareihin, jotka olet nähnyt aiemmin. Yleisenä ajatuksena on laskea, kuinka monta kertaa todelliset esiintymät luokitellaan vääriksi.

Sekaannusmatriisin laskemiseksi sinulla on ensin oltava joukko ennusteita, jotta niitä voidaan verrata todellisiin tavoitteisiin.

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

Code Selitys

• ennustaa(logit,data_testi, tyyppi = 'response'): Laske testijoukon ennuste. Aseta tyyppi = 'response' vastauksen todennäköisyyden laskemiseksi.
• table(data_test$income, ennustaa > 0.5): Laske sekaannusmatriisi. ennustaa > 0.5 tarkoittaa, että se palauttaa 1, jos ennustetut todennäköisyydet ovat yli 0.5, muuten 0.

lähtö:

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

Sekaannusmatriisin jokainen rivi edustaa todellista tavoitetta, kun taas jokainen sarake edustaa ennustettua tavoitetta. Matriisin ensimmäinen rivi käsittelee alle 50 000 euron tuloja (negatiivinen luokka): 6 310 havaintoa luokiteltiin oikein alle 50 000 euron tulojen omaaviksi yksilöiksi (Tosi negatiivinen), kun taas 495 luokiteltiin virheellisesti yli 50 000:n arvoisiksi (Väärä positiivinen). Toisella rivillä käsitellään yli 50 000 euron tuloja: 1 229 tunnistettiin oikein (Tosi positiivista), kun taas 1 074 jäi väliin (Väärä negatiivinen).

Voit laskea mallin tarkkuus summaamalla todellinen positiivinen + todellinen negatiivinen kokonaishavainto

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

Code Selitys

• summa(diag(taulukko_matto)): Diagonaalin summa
• sum(table_mat): Matriisin summa.

lähtö:

## [1] 0.8277339

Malli näyttää kärsivän yhdestä ongelmasta: se tuottaa liikaa vääriä negatiivisia tuloksia. Tätä kutsutaan tarkkuustestin paradoksiTotesimme, että tarkkuus on oikeiden ennusteiden suhde tapausten kokonaismäärään. Meillä voi olla suhteellisen korkea tarkkuus, mutta malli on hyödytön. Näin tapahtuu, kun on olemassa dominoiva luokka. Jos katsot taaksepäin sekaannusmatriisiin, voit nähdä, että useimmat tapaukset luokitellaan todella negatiivisiksi. Kuvittele nyt, että malli luokittelee jokaisen havainnon negatiiviseksi (eli alle 50 000). Saat silti noin 75 prosentin tarkkuuden (6 805 / 9 108). Mallisi toimii paremmin, mutta sillä on vaikeuksia erottaa todellinen positiivinen todellisesta negatiivisesta.

Tällaisessa tilanteessa on parempi käyttää tiiviimpää mittaria. Voimme katsoa:

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

Tarkkuus vs muistaminen

Tarkkuus tarkastelee positiivisen ennusteen tarkkuutta. Palauttaa mieleen on positiivisten tapausten suhde, jotka luokittelija havaitsee oikein;

Voit rakentaa kaksi funktiota näiden kahden metriikan laskemiseksi

1. Rakenna tarkkuus

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

Code Selitys

• mat[1,1]: Palauta tietokehyksen ensimmäisen sarakkeen ensimmäinen solu eli tosipositiivinen
• matto[1,2]; Palauta tietokehyksen toisen sarakkeen ensimmäinen solu eli väärä positiivinen

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

Code Selitys

• mat[1,1]: Palauta tietokehyksen ensimmäisen sarakkeen ensimmäinen solu eli tosipositiivinen
• matto[2,1]; Palauta tietokehyksen ensimmäisen sarakkeen toinen solu eli väärä negatiivinen

Voit testata toimintojasi

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

lähtö:

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

Lue nämä kaksi lukua huolellisesti. Tarkkuus on 0.71, joten kun malli sanoo yksilön ansaitsevan yli 50 000, se pitää paikkansa 71 prosentissa tapauksista. Muistutuskerroin on 0.53, joten malli havaitsee vain 53 prosenttia yksilöistä, jotka todella ansaitsevat yli 50 000.

Voit luoda pisteet tarkkuuden ja muistamisen perusteella. The on näiden kahden mittarin harmoninen keskiarvo, mikä tarkoittaa, että se antaa enemmän painoa alemmille arvoille.

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

lähtö:

## [1] 0.6103799

Precision vs Recall kompromissi

On mahdotonta saavuttaa sekä suurta tarkkuutta että korkeaa muistia.

Jos lisäämme tarkkuutta, oikea yksilö ennustetaan paremmin, mutta monet niistä jäävät paitsi (alempi muistaminen). Joissakin tilanteissa suosittelemme suurempaa tarkkuutta kuin muistamista. Tarkkuuden ja muistamisen välillä on kovera suhde.

• Kuvittele, sinun täytyy ennustaa, onko potilaalla sairaus. Haluat olla mahdollisimman tarkka.
• Jos haluat havaita mahdolliset petolliset ihmiset kadulla kasvojentunnistuksen avulla, olisi parempi saada kiinni monet petollisiksi leimatut ihmiset, vaikka tarkkuus on alhainen. Poliisi voi vapauttaa petollisen henkilön.

ROC-käyrä

Vastaanotin Operating Ominaisuus käyrä on toinen yleinen työkalu, jota käytetään binääriluokituksessa. Se on hyvin samankaltainen kuin tarkkuus/palautuskäyrä, mutta sen sijaan, että piirrettäisiin tarkkuus vs. muistaminen, ROC-käyrä näyttää todellisen positiivisen määrän (eli palauttamisen) väärän positiivisen nopeuden suhteen. Väärien positiivisten tapausten osuus on virheellisesti positiivisiksi luokiteltujen negatiivisten tapausten suhde. Se on yhtä suuri kuin yksi miinus todellinen negatiivinen korko. Todellista negatiivista korkoa kutsutaan myös erityispiirteet. Tästä syystä ROC-käyrä piirretään herkkyys (muistuta) vs. 1-spesifisyys

ROC-käyrän piirtämiseksi meidän on asennettava ROCR-niminen paketti. Löydämme sen condasta kirjasto. Voit kirjoittaa koodin:

conda install -c r r-rocr --yes

Voimme piirtää ROC:n ennuste()- ja performance()-funktioilla.

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 Selitys

• ennuste(ennustus, data_testi$tulo): ROCR-kirjaston on luotava ennusteobjekti syöttötietojen muuttamiseksi
• suorituskyky (ROCRpred, 'tpr', 'fpr'): Palauta kaksi yhdistelmää, jotka tuotetaan kaaviossa. Tässä konstruoidaan tpr ja fpr. Tot piirrä tarkkuus ja muista yhdessä, käytä "pre", "rec".

lähtö:

Vaihe 8) Paranna mallia

Voit yrittää lisätä malliin epälineaarisuutta vuorovaikutuksella

• ikä ja tunnit viikossa
• sukupuoli ja tunnit.per.viikko.

Sitten vertaat molempien mallien F1-pisteitä.

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

lähtö:

## [1] 0.6109181

F1-pisteet ovat hieman edellistä korkeammat. Voit jatkaa datan parissa työskentelyä ja yrittää parantaa pistemäärää.

GLM-kertoimien ja kerroinsuhteiden tulkitseminen R:ssä

Vaiheessa 6 tulostettu yhteenvetotaulukko raportoi kertoimet log-kertoimet mittakaavassa, jota on vaikea selittää ei-tekniselle yleisölle. Niiden muuntaminen kerroinsuhteiksi tekee mallista paljon helpommin kommunikoitavan.

Noudata näitä neljää vaihetta.

1. Korota kertoimet potenssiin. Sovella exp()-funktiota jokaiseen estimaattiin, jotta log-kertoimista tulee kertovia kerroinsuhteita.
2. Lisää luottamusväli. Sisällytä confint()-funktio exp()-funktioon saadaksesi 95 prosentin välin samalla todennäköisyysasteikolla.
3. Vertaa kutakin arvoa arvoon 1. Yli 1:n kerroinsuhde lisää positiivisen luokan todennäköisyyttä, alle 1:n arvo pienentää sitä ja lähellä 1:tä oleva arvo tarkoittaa, että ennustaja lisää vain vähän.
4. Tarkista tilastollinen merkitsevyys. Tulkitse vain ennustavia tekijöitä, joiden p-arvo yhteenvetotuloksessa on valitsemasi kynnysarvon, tyypillisesti 0.05, alapuolella.

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

Tulosteen lukeminen. Mallissamme kerroin tunteja viikossa on 0.41942. Korottamalla se potenssiin saadaan exp(0.41942) = 1.52, mikä tarkoittaa, että yhden keskihajonnan lisäys viikoittaisissa työtunneissa kertoo yli 50 000 tunnin ansaitsemisen todennäköisyyden noin 1.5:llä, pitäen kaikki muut muuttujat vakioina.

Negatiiviset kertoimet toimivat samalla tavalla. marital.statusNot_married on -2.50346, joten exp(-2.50346) = 0.08: naimattomilla henkilöillä on noin 8 prosenttia naimisissa olevan henkilön todennäköisyydestä. Koska jatkuvat ennustajat standardoitiin vaiheessa 1, kuvaa muutoksia keskihajontayksiköissä, älä raakatunneina.

Huomio muille perheille: Potenssikorotetut kertoimet ovat ristitulosuhteita vain binomimaisessa perheessä, jossa on logit-linkki. Kun perhe = ”poisson” ja linkki on logaritminen, samat exp()-arvot luetaan lukusuhteina.

GLM R:ssä: Funktioiden pikaopas

Pidä tätä taulukkoa vieressäsi koodatessasi. Se listaa kaikki yllä olevissa kahdeksassa vaiheessa käytetyt funktiot sekä funktion toimittavan paketin ja sen odottamat argumentit.

Muut GLM Perheväittelyn kautta saatavilla olevat perheet ovat:

• binomiaalinen: (linkki = “logit”)
• gaussinen: (linkki = “identiteetti”)
• Gamma: (linkki = “käänteinen”)
• käänteinen.gaussian: (linkki = “1/mu^2”)
• poisson: (linkki = “loki”)
• kvasi: (linkki = “identiteetti”, varianssi = “vakio”)
• kvasibinomiaalinen: (linkki = “logit”)
• kvasipoisson: (linkki = “log”)

UKK