Литература

Логистичка регресија

У линеарном регресионом моделу смо моделовали зависност средње вриједности зависне промјенљиве \(E(Y)\) (односно \(E(Y|X)\)) од предиктора. Поставља се питање да ли то можемо да урадимо у случају да \(Y\) нема нормалну расподјелу или њој сличну. Одговор је потврдан у великом броју случајева, али уз мале модификације. У овом дијелу ћемо размотрити случај када је \(Y\) индикатор - неко обиљежје које узима само двије вриједности \(0\) и \(1\), и приказати како се може направити модел за \(p_i = P\{Y_i=1\}\). Само неки од примјера су:

Претпоставимо да је \(Y\) нека зависна промјенљива чија се средња вриједност може моделирати линеарним моделом са нормално расподјељеним грешкама, односно да посматрамо модел \[Y_i = \beta_0+\beta_1x_i+\epsilon_i.\]

На примјер, нека је промјенљива \(Y\) ниво ријечног корита. Умјесто узорака \(Y_1\), \(Y_2\),…, \(Y_n\) на располагању имамо само информацију да ли је вриједност \(Y_i\) прешла критични ниво \(c\), након чега се покреће ванредна одбрана од поплава. Дакле, имамо узорак \(Y_1^c\), \(Y_2^c\),…, \(Y_n^c\), гдје је \[Y_i^c = \begin{cases} 1, & Y_i>c; \\ 0, & Y_i \leq c. \end{cases}\]

Желимо да направимо модел којим ћемо оцијенити вјероватноћу да је \(Y_i\) веће од нивоа \(c\), у зависности од неког предиктора \(x_i\), на примјер просјечне дневне температуре у претходних \(5\) дана. Тада је

\[p_i = P\{Y_i^c=1\} = P\{\epsilon_i > c - \beta_0-\beta_1x_i\} = 1 - \Phi \left(\frac{c-\beta_0-\beta_1x_i}{\sigma}\right) = \Phi \left(\frac{-c+\beta_0+\beta_1x_i}{\sigma}\right)\]

гдје је \(\sigma^2 = D(\epsilon_i)\). Одавде је

\[\Phi^{-1}(p_i) = \frac{-c+\beta_0+\beta_1x_i}{\sigma}=A+Bx_i, \]

гдје је \(A=\frac{-c+\beta_0}{\sigma}\) и \(B=\frac{\beta_1}{\sigma}\). Дакле, трансформација коју смо примијенили је \(\Phi^{-1}(\cdot)\).

Тип регресије приказане у претходном разматрању (када примјењујемо трансформацију инверз функције расподјеле нормално расподјељене случајне величине, односно \(\Phi^{-1}(\cdot)\)), назива се пробит регресије. У случају да се ради о стандардној логистичкој расподјели, односно \(F(x) = \frac{1}{1+e^{-x}}\), за \(x\in \mathbb{R}\), чија је инверзна функција \(F^{-1}(x) = \log \frac{x}{1-x}\), за \(x\in (0,1)\), модел

\[\log \frac{p_i}{1-p_i}=\beta_0+\beta_1x_i\]

се назива логистички регресиони модел. Функција

\[\lambda (x) = \log \frac{p(x)}{1-p(x)} \]

се назива логит трансформација. Количник \(\frac{p(x)}{1-p(x)}\) се назива кота (енг. odds ration). На примјер, уколико је квота 4, то значи да је, за фиксирану вриједност предиктора \(X=x\), вјероватноћа да је \(Y=1\), \(4\) пута већа од вјероватноће да је \(Y=0\).

Управо због јасне интерпретације утицаја предиктора на квоту, логистички регресиони модел се најчешће и користи у случају бинарне зависности промјенљиве. Наиме, уколико се предиктор повећа за 1, квота се промијени \(e^{\beta_1}\) пута. Зато оцјењујући коефицијент \(\beta_1\) можемо да видимо какав је утицај промјене предиктора на саму промјену квоте. Примијетимо да ако је \(\beta_1<0\) онда је \(e^{\beta_1}<1\), па долази до смањења квота, док у супротном долази до њеног повећања.

Из логит трансформације добија се да је

\[p(x) = \frac{1}{1+e^{-\lambda (x)}}.\]

Параметре \(\beta_0\) и \(\beta_1\) оцјењујемо методом максималне вјеродостојности. Функција вјеродостојности дата је са

\[L(\beta_0, \beta_1) = \prod_{i=1}^n p_i^{y_i}(1-p_i)^{1-y_i}.\]

Одавде је

\[l(\beta_0, \beta_1) = \sum_{i=1}^n \left(y_i \log \frac{p_i}{1-p_i} +\log (1-p_i)\right) = \sum_{i=1}^n \left( y_i(\beta_0+\beta_1x_i) + \log\frac{1}{1+e^{\beta_0+\beta_1x_i}}\right).\]

Рјешавање система \(\frac{\partial l(\beta_0, \beta_1)}{\partial \beta_0} = \frac{\partial l(\beta_0, \beta_1)}{\partial \beta_1} = 0\) се врши нумерички и нећемо се на томе задржавати.

У \(R\)-у се за прављење модела користи функција glm() (енг. Generalized linear models), која поред аргумента којим се задаје посматрана зависност садржи још један важан елемент којим се спецификује да се ради о бинарној зависној промјенљивој и логит трансформацији: familiy = binomial(link=“logit”).

Након што оцијенимо параметре \(\hat{\beta}_0\) и \(\hat{\beta}_1\), оцијењена логит функција је \[\hat{\lambda}(x) = \hat{\beta}_0 + \hat{\beta}_1x.\]

Њеном трансформацијом добијамо да је оцјена вјероватноће да је \(Y=1\) када је вриједност предиктора \(x\) једнака \[\hat{p}(x) = \frac{1}{1+e^{-(\hat{\beta}_0+\hat{\beta}_1x)}}\] односно \[\hat{p}_i = \frac{1}{1+e^{-(\hat{\beta}_0+\hat{\beta}_1x_i)}}.\]

Примјер 1. У јануару \(1986\). године дошло је до експлозије спејс-шатла \(Challenger\) у коме је живот изгубило \(7\) чланова посаде. као један од разлога катастрофе наводе се лоше перформансе спољних дихтунга на које је утицала изузетно ниска спољашња температура \(31F\). Како би се испитала веза између попуштања дихтунга и спољашње температуре лансирано је \(6\) летјелица у различитим условима. Резултати се могу наћи у бази \(challenger\) (пакет \(mcsm\)). Промјенљива \(oring\) узима вриједност \(0\) уколико није дошло до квара, а \(1\) ако је дошло до квара, а промјенљива \(temp\) означава температуру.

challenger <- read.delim("challenger.txt")
model <- glm(oring~temp, family = "binomial", data = challenger)
summary(model)
## 
## Call:
## glm(formula = oring ~ temp, family = "binomial", data = challenger)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.0611  -0.7613  -0.3783   0.4524   2.2175  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)  15.0429     7.3786   2.039   0.0415 *
## temp         -0.2322     0.1082  -2.145   0.0320 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 28.267  on 22  degrees of freedom
## Residual deviance: 20.315  on 21  degrees of freedom
## AIC: 24.315
## 
## Number of Fisher Scoring iterations: 5
coefficients(model)
## (Intercept)        temp 
##  15.0429016  -0.2321627

Добили смо модел

\[\hat{p}(x) = \frac{1}{1+e^{-(15.04-0.23x)}},\] на основу ког моћемо оцијенити вјероватноћу да дође до отказивања при температури која је била.

predict(model, data.frame(temp=31), type='response')
##         1 
## 0.9996088

Добили смо да је оцијењена вјероватноћа блиска \(1\) што упућује на то да је израда овог истраживања прије самог лансирања спејс-шатла могла спријечити катастрофу.

Вишеструка логистичка регресија

Потпуно аналогно поступамо у случају да имамо више од једног предиктора (што се најчешће дешава у пракси). Тада посматрамо модел:

\[\log \frac{p_i}{1-p_i} = \beta_0+\beta_1 x_{i1}+\beta_2x_{i2}+...+\beta_px_{ip}\]

и параметре оцјењујемо методом максималне вјеродостојности (нумерички). Након тога можемо оцијенити и логит функцију са

\[\hat{\lambda}(x) = \hat{\beta}_0+\hat{\beta}_1x_1+...+\hat{\beta}_px_p,\]

чијом трансформацијом се може добити и оцијењена вјероватноћа.

Примјер 2. Посматрајмо базу податак која садржи резултате студената на пријемним испитима на неком факултету за упис на мастер студије (прве двије колоне), док трећа колона представља индикатор да ли је посматрани студент уписао факултет или не.

baza <- read.csv("http://www.math.rs/p/files/1458635682-70-baza.txt", header=FALSE)
head(baza)
##         V1       V2 V3
## 1 34.62366 78.02469  0
## 2 30.28671 43.89500  0
## 3 35.84741 72.90220  0
## 4 60.18260 86.30855  1
## 5 79.03274 75.34438  1
## 6 45.08328 56.31637  0

Очигледно да је промјенљива \(V3\) категоричка (узима само вриједности \(0\) и \(1\)), па има смисла да примијенимо логистичку регресију.

par(mfrow=c(1,2))
plot(V3~V1+V2, data = baza)

Модел логистичке регресије за вјероватноћу \(p\) да студент упише жељени факултет дат је изразом

\[\log \left(\frac{p}{1-p}\right)=\beta_0+\beta_1V_1+\beta_2V_2.\]

Слично претходним примјерима, оцјењујемо вриједности коефицијената \(\beta_0\), \(\beta_1\) и \(\beta_2\).

model.log <- glm(V3 ~ V1 + V2, family = "binomial", data = baza)
summary(model.log)
## 
## Call:
## glm(formula = V3 ~ V1 + V2, family = "binomial", data = baza)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -2.19287  -0.18009   0.01577   0.19578   1.78527  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -25.16133    5.79836  -4.339 1.43e-05 ***
## V1            0.20623    0.04800   4.297 1.73e-05 ***
## V2            0.20147    0.04862   4.144 3.42e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 134.6  on 99  degrees of freedom
## Residual deviance:  40.7  on 97  degrees of freedom
## AIC: 46.7
## 
## Number of Fisher Scoring iterations: 7

Из претходног добијамо да је тражени модел облика

\[\log \left(\frac{p}{1-p}\right)=−25.16133+0.20623V_1+0.20147V_2.\]

Даље, желимо на основу добијеног модела да направимо предвиђање. Нека је познато да је студент освојио на прим+јемном испиту редом \(V_1=20\) и \(V_2=100\) поена. Занима нас да ли са овим резултатима уочени студент може уписати мастер. На основу датих вриједности имамо да је оцијењена вјероватноћа једнака

\[p = \frac{1}{1+e^{-(-25.16133+0.20623\cdot 20+0.20147\cdot 100)}}=0.29.\]

Претходни резултат добијамо и на сљедећи начин

predict(model.log, newdata =  data.frame(V1 = 20, V2 = 100), type = "response") 
##         1 
## 0.2912049
predict(model.log, newdata =  data.frame(V1 = 20, V2 = 100), type = "response") > 0.5
##     1 
## FALSE

Из претходног добијамо да је предвиђена вриједност једнака \(0\), што значи да упчени студент неће уписати факултет.

Класификација

На основу оцијењеног логистичког регресионог модела можемо вршити класификацију. Једна могућност је да класификујемо податке у двије класе \(0\) и \(1\), ако је \(\hat{p}_i>0.5\) онда је \(\hat{Y}_i = 1\), у супротном је \(\hat{Y}_i=0\). Након тога можемо видјети проценат добро класификованих података. Ово је природан приступ и као такав најпрепоручљивији. Међутим, уколико не добијемо задовољавајућу тачност, а битно нам је да задржимо модел, на примјер због интерпретабилности, можемо да помјеримо праг за класификацију. Односно,

\[\hat{Y}_i = \begin{cases} 0, & \hat{p}_i <C\\ 1, & \hat{p}_i \geq C \end{cases}\]

а праг за класификацију \(C\) бирамо имајући у виду мјере тачности које ћемо описати у наставку.

Резултат предвиђања се може показати такозваном матрицом конфузије.

\(Y\backslash \hat{Y}\) \(0\) \(1\)
\(0\) \(TN - true \, negative\) \(FP - false \, positive\)
\(1\) \(FN - false \, negative\) \(TP - true \, positive\)

Са \(TN\) (енг. true negative) смо означили број оних елемената чија је вриједност \(0\) и који су класификовани као \(0\), \(FP\) (енг. false positive) је број нула које су погрешно класификоване као \(1\), \(FN\) (енг. false negative) је број јединица које су погрешно оцијењене са \(0\), и \(TP\) (енг. true positive) број јединица које су тачно класификоване као \(1\). Тачност класификације (енг. accuracy) је сада \[A=\frac{TN+TP}{TN+FP+FN+TP}.\]

Ова мјера није адекватна ако класе нису приближних величина. Уколико је једна класа знатно већа од друге, класификовањем свих елемената тако да припадају тој већој класи постижемо велику тачност, иако је јасно да класификатор није добар. Због тога се дефинишу и друге мјере. Неке од њих су:

У зависности од афинитета, \(C\) се може бирати тако да поменуте мјере имају екстремну вриједност. Наравно, могу се и комбиновати. На примјер, прецизност нам даје податак о броју релевантних јединица (тачно класификованих јединица) јер обично је догађај чију успјешну реализацију означавамо јединицом, тај који посматрамо. С друге стране, одзив нам говори о проценту јединица које су успјешно класификоване. На примјер, ако је \(\hat{Y}=1\) означава да је особа проглашена зараженом неким вирусом, прецизност нам говори колики проценат дијагностификованих су стварно заражени, док одзив пак, колико се добро на основу улазних параметара, добро дијагностификује особа заражена вирусом. Скор представља хармонијску средину ове двије мјере и добро га је користити када су оне равноправне. Уколико то нису може се дефинисати и тежинска хармонијска средина измежу ове двије мјере при чему ћемо оној важнијој задати већу тежину.

За одређивање прага се може одредити и такозвана \(ROC\) крива (енг. receiver operating characteristic curve) која представља зависност између \(TPR\) и \(FPR\), односно на \(y\)-оси је сензитивност, а на \(x\)-оси специфичност. Како нам тачка горњи лијеви угао графика (што се у пракси никада не постиже) говори о идеалном класификатору, једна могућност је да се одабере тачка са криве која је најближа тачки \((0,1)\).

Површина испод \(ROC\) криве \((AUC)\) нам даје информацију о квалитету класификатора. Што је већа површина класификатор боље раздваја категорије. Заправо, \(AUC\) је оцјена вјероватноће да при случајном избору два елемента из различитих класа она из класе означене са \(0\) има мању вриједност (на основу које је подијељена у класе) од оне из класе означене као \(1\).

Примјер 3. Проучавано је како висина утиче на то да ли ће особе склоне кошарци погодити кош са одређене раздаљине. Узет је узорак од \(100\) особа и резултати су приказани у бази \(basketball\). Направићемо одговарајући логистички модел. Како бисмо што непристрасније испитали квалитет модела (и како је крајњи циљ класификација) прво ћемо податке подијелити на тренинг скуп (\(70\%\) података), а затим тестирати квалитет модела на преосталим подацима (\(30\%\)). Уколико се испостави да је модел задовољавајућег квалитета, на комплетном скупу података ћемо направити финални модел.

Најприје, податке дијелимо на тренинг и тест скуп.

basketball <- read.delim("basketball.txt")

set.seed(216)
n <- length(basketball$height)
obim <- floor(n*0.7)
indikator <- sample(1:n, size = obim)

train_basketball <- basketball[indikator,]
test_basketball <- basketball[-indikator,]

Затим на тренинг скупу правимо модел.

modelT <- glm(shot ~ height, family = "binomial", data = train_basketball)
summary(modelT)
## 
## Call:
## glm(formula = shot ~ height, family = "binomial", data = train_basketball)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -2.47890   0.07753   0.19088   0.49124   1.66685  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -44.806     12.638  -3.545 0.000392 ***
## height        24.280      6.697   3.626 0.000288 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 64.14  on 69  degrees of freedom
## Residual deviance: 37.35  on 68  degrees of freedom
## AIC: 41.35
## 
## Number of Fisher Scoring iterations: 6

На основу модела на тренинг скупу можемо оцијенити вјероватноће на тест скупу. Предиктор \(height\) смо задржали у моделу јер је \(p\)-вриједност једнака \(0.000288\). Даље, можемо нацртати \(ROC\) криву и одредити површину испод ње \(AUC\). За то користимо функцију roc() из пакета pROC.

library("pROC")
test_basketball$prob = predict(modelT, newdata = test_basketball, type='response')
basketballROC = roc(shot~height, data=test_basketball)
auc(basketballROC)
## Area under the curve: 0.9062
plot(basketballROC)

Одредићемо и матрицу конфузије.

test_basketball$class = (test_basketball$prob > 0.5)
table(test_basketball$shot, test_basketball$class)
##    
##     FALSE TRUE
##   0     4    2
##   1     0   24

На основу матрице конфузије видимо да је за одабрани праг класификације \((0.5)\) добијена сензитивност \(24/24=1\), док је специфичност \(4/6=0.66\) што говори у прилог квалитету модела.

Добијене резултате можемо упоредити са резултатима на тренинг скупу. Модел, уколико је добар (и ако су скупови за тренинг и тест добро одабрани) требало би да показује приближан квалитет на оба скупа.

train_basketball$prob = fitted(modelT)
basketballROCT = roc(shot~height, data = train_basketball)
auc(basketballROCT)
## Area under the curve: 0.9116
plot(basketballROCT)

Након што смо испитали понашање модела на оба скупа и утврдили да разлике нису велике, иако постоје, направићемо модел на свим подацима и на основу њега приказати оцијењену зависност вјероватноће (укључујући и интервале повјерења) поготка од висине особе.

model = glm(shot~height, family = "binomial", data = basketball)
xheight = seq(from = 1.6, to = 2.2, by = 0.01)
probEstimate = predict(model, newdata = data.frame(height = xheight), type = 'response')

probEstimateL = predict(model, newdata = data.frame(height = xheight), type = 'link', se.fit = TRUE)

Одређујемо горњу и доњу границу интервала повјерења за логит функцију а затим их трансформишемо у границе за вјероватноћу успјеха.

C=qnorm(0.975)
intUpper = plogis(probEstimateL$fit + C*probEstimateL$se.fit)
intLower = plogis(probEstimateL$fit - C*probEstimateL$se.fit)

plot(xheight, probEstimate, type='l', lwd=2, xlab='x', ylab='p')

lines(xheight, intUpper, col="blue", lwd=2)
lines(xheight, intLower, col="blue", lwd=2)

Приликом одређивања интервала користили смо да оцијењена логит функција \(\hat{\lambda}(x)\) има асимптотски нормалну расподјелу.

Из методолошких разлога илустроваћемо и како се \(C\) може одредити коришћењем \(ROC\) криве. Тренинг скуп ћемо подијелити на два скупа, водећи рачуна о односу класа зависне промјенљиве. На првом (већем дијелу) ћемо направити модел, а онда на основу предвићених вјероватноћа на том моделу, нацртати \(ROC\) криву на другом дијелу скупа који се још назива и скуп за валидацију. На тест скупу ћемо одредити матрицу конфузије на основу модела направљеног на комплетном тренинг скупу.

set.seed(15)
library("caTools")

indexTV = sample.split(train_basketball$shot, SplitRatio = 7/10)
basketballTT = train_basketball[indexTV, ]
basketballTestV = train_basketball[!indexTV, ]

modelTT = glm(shot~height, family = "binomial", data = basketballTT)

basketballTestV$prob = predict(modelTT, newdata = basketballTestV, type = 'response')

basketballROCTV = roc(shot~prob, data = basketballTestV)

auc(basketballROCTV)
## Area under the curve: 0.875
plot(basketballROCTV)

Одређујемо тачку најближу горњем лијевом углу.

coords(basketballROCTV, x='best', best.method='closest.topleft', input='threshold')
##   threshold specificity sensitivity
## 1 0.8543604           1   0.7647059

Упоређујемо \(AUC\) на скупу на коме смо развијали модел и на тесту за валидацију.

basketballTT$prob = predict(modelTT, newdata = basketballTT, type = 'response')

basketballROCTT = roc(shot~prob, data = basketballTT)
auc(basketballROCTT)
## Area under the curve: 0.9192

За крај одређујемо матрицу конфузије на тест скупу.

test_basketball$class1 = (test_basketball$prob > 0.8543604)
table(test_basketball$shot, test_basketball$class1)
##    
##     FALSE TRUE
##   0     5    1
##   1     5   19

Задаци за самосталан рад

Задатак 1. У пакету \(mlbench\) налази се база података \(PimaIndiansDiabetes2\) о присуству дијабетеса и неких других клиничких показатеља који потенцијално утичу на присуство шећера у крви. Направити одговарајући модел логистичке регресије који се може користити за предвиђање вјероватноће болести, односно за одговарајућу класификацију.

Задатак 2. У бази \(malaria\) (пакет \(ISwR\)) налазе се подаци о \(100\) случајно одабраних дјечијих испитаника старости \(3\) до \(15\) година из једног села у Гани.

а) Подијелити податке на тренинг и тест скуп.

б) На тренинг скупу направити модел зависности вјероватноће испољавања симптома маларије \((mal)\) од старости \((age)\) и количине антитијела у крви \((ab)\).

в) Одредити \(AUC\) претходно направљеног модела на тест скупу. Шта можете да закључите о класификационој моћи модела?