Logistic Regression and Classification

Wczytanie wymaganych pakietów.

library(ISLR)
data(Default)

1 Binaryzacja wartości wyliczeniowych

Domyślna wartość dummy variable.

print(contrasts(Default$default))

    Yes
No    0
Yes   1

print(contrasts(Default$student))

    Yes
No    0
Yes   1

Możliwe jest bezpośrednie zdefiniowanie wartości dummy variable w przypadku, gdy funkcja contrasts zwraca niepożądane wartości.

data<-cbind(Default, data.frame(student.Yes=ifelse(Default$student=="Yes", 1, 0),
                                default.Yes=ifelse(Default$default=="Yes", 1, 0)))

2 Model regresji logistycznej

Funkcja logistyczna: \(p(X)=\frac{e^{\beta_0+\beta_1X_1+\dotsb+\beta_pX_p}}{1+e^{\beta_0+\beta_1X_1+\dotsb+\beta_pX_p}}\).

2.1 Trenowanie modelu Default ~ Balance

est<-glm(default~balance, data=Default, family=binomial(link="logit"))

Call:
glm(formula = default ~ balance, family = binomial(link = "logit"), 
    data = Default)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.2697  -0.1465  -0.0589  -0.0221   3.7589  

Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept) -1.065e+01  3.612e-01  -29.49   <2e-16 ***
balance      5.499e-03  2.204e-04   24.95   <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: 2920.6  on 9999  degrees of freedom
Residual deviance: 1596.5  on 9998  degrees of freedom
AIC: 1600.5

Number of Fisher Scoring iterations: 8

2.2 Współczynniki modelu regresji logistycznej.

coef(est)  # to samo co print(summary(est)$coefficients[,"Estimate"])

  (Intercept)       balance 
-10.651330614   0.005498917 

2.3 Prognoza

type="response" powoduje zwrócenie wartości w postaci \(P(Y=1 \mid X)\).

sprintf("%.2f%%", predict(est, newdata = data.frame(
  balance=c(median(Default$balance), 1500, 1900, max(Default$balance))), type="response") * 100)

[1] "0.22%"  "8.29%"  "44.93%" "98.10%"

2.4 Prawdopodobieństwo, Logit i szansa

2.4.1 Prawdopodobieństwo

Definicja funkcji prawdopodobieństwo \(p(X)\).

B1=coef(est)[["balance"]]
B0=coef(est)[["(Intercept)"]]
p<-function(x) exp(B1*x+B0)/(1+exp(B1*x+B0))

2.4.2 Szansa (ang. odds)

\(\frac{p(X)}{1-p(X)}=e^{\beta_0+\beta_1X_1+\dotsb+\beta_pX_p}\).

Definicja funkcji szansa.

Odds<-function(x) exp(B0+B1*x)

2.4.3 Logit

\(\log\left(\frac{p(X)}{1-p(X)}\right)=\beta_0+\beta_1X_1+\dotsb+\beta_pX_p\).

Definicja funkcji Logit.

Logit<-function(x) B0+B1*x

Wartości prawdopodobieństwa, szansy oraz funkcji Logit dla nowej obserwacji zadłużenie = $1900.

sprintf("p(1900)=%.2f%%", 100*p(1900))
paste("Odds(1990)=", signif(Odds(1900), 5), sep = "")
paste("Logit(1900)=", signif(Logit(1900), 5), sep="")

p(1900)=44.93%
Odds(1990)=0.81596
Logit(1900)=-0.20339

Te same wartości z użyciem funkcji predict.glm.

predict.glm(est, newdata = data.frame(balance=1900), type="response")
exp(predict(est, newdata = data.frame(balance=1900)))
predict.glm(est, newdata = data.frame(balance=1900))

0.4493274
0.8159612
-0.2033884

3 Wizualizacja

Wykres dystrybuanty rozkładu logistycznego. Na osi rzędnych prawdopodobieństwo problemów ze spłatą karty (default) w zależności od zadłużenia (balance) na osi odciętych. Kolorem czerwonym zaznaczone są obserwacje klientów, którzy faktycznie karty nie spłacają (default); czarnym natomiast klienci regulujący zobowiązania.

4 Klasyfikacja

Trafność (ang. accuracy) klasyfikacji dychotomicznej zbioru trenującego przy zadanym progu prawdopodobieństwa 50%.

Default.pred<-rep("No", nrow(Default))
Default.pred[Default.probs > .5]<-"Yes"
table(Default.pred)

Default.pred
  No  Yes 
9858  142 

Macierz pomyłek (ang. confusion matrix / truth table).

print(table(Default.pred, Default$default))

            
Default.pred   No  Yes
         No  9625  233
         Yes   42  100

Trafność.

sprintf("%.2f%%", 100 * mean(Default.pred==Default$default))

[1] "97.25%"

4.1 Podział na zbiór trenujący i walidujący metodą hold-out

Trenowanie modelu.

set.seed(30112)                              # ustawienie ziarna generatora liczb losowych.
library(caret)
partition<-createDataPartition(
  Default$default, p=.9, list=FALSE)
est<-glm(default~balance, data=Default,
         family=binomial, subset=partition)

Alternatywnie.

Default.train<-rep(FALSE, nrow(Default))
Default.train[sample(1:nrow(Default), 0.9*nrow(Default))]<-TRUE
est<-glm(default~balance, data=Default, family=binomial, subset=Default.train)

Call:
glm(formula = default ~ balance, family = binomial, data = Default, 
    subset = partition)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.2980  -0.1432  -0.0569  -0.0211   3.7852  

Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept) -10.813875   0.387893  -27.88   <2e-16 ***
balance       0.005596   0.000236   23.71   <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: 2630.7  on 9000  degrees of freedom
Residual deviance: 1418.4  on 8999  degrees of freedom
AIC: 1422.4

Number of Fisher Scoring iterations: 8

4.2 Macierz pomyłek

Prognoza / Rzeczywista wartość

	Act.False	Act.True
Est.False	True Negative	False Negative / II
Est.True	False Positive / I	True Positive

Wektory prognozy prawdopodobieństwa oraz klasyfikacji przy zadanym progu prawdopodobieństwa 50%.

Default.probs<-predict(est, type="response", newdata = Default[-partition,])
Default.pred<-rep("No", length(Default.probs))
Default.pred[Default.probs > .5]<-"Yes"
valVect<-Default[-partition, "default"]

Alternatywnie.

Default.probs<-predict(est, type="response", newdata = Default[!Default.train,])
Default.pred<-rep("No", length(Default.probs))
Default.pred[Default.probs > .5]<-"Yes"
valVect<-Default[!Default.train,"default"]

Macierz pomyłek.

confusionMatrix<-table(Default.pred, valVect, dnn=c("Predicted","Actual"))
rownames(confusionMatrix)<-c("Predicted.No", "Predicted.Yes")
colnames(confusionMatrix)<-c("Actual.No", "Actual.Yes")
print(confusionMatrix)

               Actual
Predicted       Actual.No Actual.Yes
  Predicted.No        960         25
  Predicted.Yes         6          8

4.3 Współczynniki klasyfikacji

P<-sum(confusionMatrix[,"Actual.Yes"])
N<-sum(confusionMatrix[,"Actual.No"])
TP<-confusionMatrix["Predicted.Yes", "Actual.Yes"]
FP<-confusionMatrix["Predicted.Yes", "Actual.No"]
FN<-confusionMatrix["Predicted.No", "Actual.Yes"]
TN<-confusionMatrix["Predicted.No", "Actual.No"]

Alternatywnie P i N.

P<-TP+FN
N<-FP+TN

4.3.1 Trafność (ang. accuracy)

Rate.ACC<-(TP+TN)/sum(confusionMatrix)

Alternatywnie.

Rate.ACC<-(TP+TN)/(P+N)

[1] 0.968969

4.3.2 Błąd klasyfikacji (ang. error)

Rate.Err<-(FP+FN)/sum(confusionMatrix)

Alternatywnie.

Rate.Err<-(FP+FN)/(P+N)
Rate.Err<-1-ACC

[1] 0.03103103

4.3.3 Czułość / zwrot (ang. sensitivity / recall)

Rate.TPR<-TP/P

[1] 0.2424242

4.3.4 False Negative rate

Rate.FNR<-FN/P  # to samo co Rate.FNR<-1-Rate.TPR

[1] 0.7575758

4.3.5 Specyficzność (ang. specificity)

Rate.TNR<-TN/N

[1] 0.9937888

4.3.6 False Positive rate / Fall-out

Rate.FPR<-FP/N  # to samo co Rate.FPR<-1-Rate.TNR

[1] 0.00621118

4.3.7 Precyzja (ang. precision)

Rate.PPV<-TP/(TP+FP)  # to samo co Rate.PPV<-TP/sum(confusionMatrix["Predicted.Yes",])

[1] 0.5714286

4.3.8 F-miara (ang. F1 score)

Rate.F1<-2*TP/(2*TP+FP+FN)  # to samo co Rate.F1<-2*TP/(sum(Default.pred!=valVect)+2*TP)

[1] 0.3404255

4.4 Zastosowanie pakietu SDMTools

Macierz pomyłek z użyciem funkcji SDMTools::confusion.matrix.

library("SDMTools")
confusion.matrix(ifelse(valVect=="Yes", 1, 0), ifelse(Default.pred=="Yes", 1, 0))

    obs
pred   0  1
   0 960 25
   1   6  8
attr(,"class")
[1] "confusion.matrix"

5 Zbiór walidujący

5.1 Wizualizacja

score_data<-data.frame(balance=Default[-partition, "balance"],
                       p=Default.probs, default=valVect)

Alternatywnie.

score_data<-data.frame(balance=Default[!Default.train, "balance"],
                       p=Default.probs, default=valVect)

Wykres dystrybuanty rozkładu logistycznego dla zbioru walidującego. Na osi rzędnych prawdopodobieństwo problemów ze spłatą karty (default) w zależności od zadłużenia (balance) na osi odciętych. Kolorem czerwonym zaznaczone są obserwacje klientów, którzy faktycznie karty nie spłacają (default); czarnym natomiast klienci regulujący zobowiązania.
Można zaobserwować jak wraz ze zmianą progu prawdopodobieństwa klasyfikacji zmieniać się będą współczynniki precyzji oraz czułości.

5.2 Krzywa ROC (Receiver Operating Characteristic curve)

Jest to wykres, na którym każdy punkt, determinowany przez wartości współczynników (False Positive rate; Czułość) powstał z osobnej macierzy pomyłek, stworzonej z użyciem innej wartości progu prawdopodobieństwa klasyfikacji (odcięcia).
Przekątna odpowiada klasyfikatorowi losowemu.

5.2.1 Stepwise ROC

rocCurve<-data.frame()
for(i in seq(from=1, to=0, by=-0.02)) {
  Default.pred<-rep("No", length(Default.probs))
  Default.pred[Default.probs > i]<-"Yes"
  same<-Default.pred==valVect
  different<-Default.pred!=valVect
  TPR<-sum(Default.pred[same]=="Yes") / sum(valVect=="Yes")
  FPR<-sum(Default.pred[different]=="Yes") / sum(valVect=="No")
  rocCurve<-rbind(rocCurve, data.frame(tpr=TPR, fpr=FPR))
}
spline<-smooth.spline(rocCurve$fpr, rocCurve$tpr, spar = 0.6)
plot(spline, ylim=c(0,1), type="b", main="ROC",
     xlab="False Positive Rate", ylab="True Positive Rate")
plot(function(x){x}, lty=5, add=TRUE)

5.2.2 Samplewise ROC

rocCurve.fnc <- function(labels, scores) {
  labels <- labels[order(scores, decreasing=TRUE)]
  data.frame(tpr=cumsum(labels)/sum(labels), fpr=cumsum(!labels)/sum(!labels), labels)
}
rocCurve<-rocCurve.fnc(valVect=="Yes", Default.probs)
plot(rocCurve$fpr, rocCurve$tpr, ylim=c(0,1), type="o", pch=20, cex=.25, main="ROC",
     xlab="False Positive Rate", ylab="True Positive Rate")
plot(function(x){x}, lty=5, add=TRUE)

5.2.3 Zastosowanie pakietu pROC

library(pROC)
rocCurve<-roc(valVect, Default.probs, direction="<", percent=TRUE)
plot(rocCurve, ylim=c(0,100), xlim=c(100,0), lwd=3)
grid()

5.3 Area Under the Curve

Całka po krzywej ROC. Wartość 1 oznacza klasyfikator idealny, 0.5 losowy; wartości poniżej 0.5 oznaczają klasyfikator gorszy od losowego.

print(rocCurve$auc)

Area under the curve: 92.87%

6 Prognoza

6.1 Default ~ Student

est<-glm(default~student, data=Default, family=binomial)

Call:
glm(formula = default ~ student, family = binomial, data = Default)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.2970  -0.2970  -0.2434  -0.2434   2.6585  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -3.50413    0.07071  -49.55  < 2e-16 ***
studentYes   0.40489    0.11502    3.52 0.000431 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 2920.6  on 9999  degrees of freedom
Residual deviance: 2908.7  on 9998  degrees of freedom
AIC: 2912.7

Number of Fisher Scoring iterations: 6

Prawdopodobieństwo a posteriori wystąpienia zadłużenia przy wiedzy (dowodzie, ang. evidence), że klient jest studentem.

predict(est, newdata = data.frame(student='Yes'), type='response')

[1] "Pr(default=Yes|student=Yes)=4.31%"

Prawdopodobieństwo a posteriori wystąpienia zadłużenia przy wiedzy (dowodzie, ang. evidence), że klient nie jest studentem.

predict(est, newdata = data.frame(student='No'), type='response')

[1] "Pr(default=Yes|student=No)=2.92%"

6.2 Default ~ Student, Balance

est<-glm(default~balance+student, data=Default, family=binomial)

Call:
glm(formula = default ~ balance + student, family = binomial, 
    data = Default)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.4578  -0.1422  -0.0559  -0.0203   3.7435  

Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept) -1.075e+01  3.692e-01 -29.116  < 2e-16 ***
balance      5.738e-03  2.318e-04  24.750  < 2e-16 ***
studentYes  -7.149e-01  1.475e-01  -4.846 1.26e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 2920.6  on 9999  degrees of freedom
Residual deviance: 1571.7  on 9997  degrees of freedom
AIC: 1577.7

Number of Fisher Scoring iterations: 8

Prawdopodobieństwo a posteriori wystąpienia zadłużenia przy wiedzy (dowodzie, ang. evidence), że klient jest studentem, a obciążenia na karcie to $2000.

predict(est, newdata = data.frame(student='Yes', balance=2000), type='response')

[1] "Pr(default=Yes|student=Yes,balance=2000)=50.30%"

Prawdopodobieństwo a posteriori wystąpienia zadłużenia przy wiedzy (dowodzie, ang. evidence), że klient nie jest studentem, a obciążenia na karcie to $2000.

predict(est, newdata = data.frame(student='No', balance=2000), type='response')

[1] "Pr(default=Yes|student=No,balance=2000)=67.41%"

7 Klasyfikacja wieloklasowa

Wartości nominalne zmiennej endogenicznej:

• DH - Disk Hernia
• SL - Spondylolisthesis
• NO - Normal

Wczytanie danych.

data<-read.table('/usr/miswdm/column_3C.dat')
colnames(data)<-c(
  'pelvic incidence',
  'pelvic tilt',
  'lumbar lordosis angle',
  'sacral slope',
  'pelvic radius',
  'grade of spondylolisthesis',
  'class')

Podział na zbiór trenujący i walidujący metodą hold-out.

partition<-createDataPartition(data$class, p=.7, list=FALSE)
data.train<-data[partition,]
data.test<-data[-partition,]

Funkcje pomocnicze.

cls<-function(name) suppressWarnings(glm(class==name~., data=data.train, family=binomial))
pred<-function(row) names(which.max(vapply(classifiers,
                                           function(c) predict(c, row, type='response'), .0)))

Metoda One-vs-All (One-vs-Rest)

classifiers=list(NO=cls('NO'), DH=cls('DH'), SL=cls('SL'))
for (i in 1:nrow(data.test)) data.test[i, 'pred']<-pred(data.test[i,])
cm<-as.matrix(table(Predicted=data.test$pred, Actual=data.test$class))

Macierz pomyłek.

         Actual
Predicted DH NO SL
       DH 10  2  1
       NO  8 28  2
       SL  0  0 42

7.1 Współczynniki klasyfikacji

Wartości pomocnicze.

diag = diag(cm)
rowsums = apply(cm, 1, sum)
colsums = apply(cm, 2, sum)

7.1.1 Trafność (ang. Accuracy)

sum(diag) / sum(cm)

[1] "Accuracy: 86.02%"

7.1.2 Średnia precyzja (ang. Mean precision)

precision = diag / rowsums
mean(precision)

[1] 0.8353576

7.1.3 Średnia czułość (ang. Mean recall)

recall = diag / colsums
mean(recall)

[1] 0.8074074

7.1.4 Średnia F-miara (ang. Mean F1 score)

f1 = 2 * precision * recall / (precision + recall)
mean(f1)

[1] 0.8114026