3 Logistische Regression

Die logistische Regression ist Beispiel für den Modellierung von abhängigen Variablen mit unterschiedlicher Skalierung, wie beispielsweise nominal oder ordinal skalierten Daten. Der “Trick” ist einen linearen Prediktor durch eine Funktion so umzuwandeln, dass der Output dieser Funktion wieder zur Skalierung der abhängigen Variable passt. Bei der logistischen Regression ist die abhängige Variable nominal skaliert, genauer hat sie die Ausprägungen 0 und 1. Geht aus einem linearen Regressionsmodell eine Prognose hervor ist dies in der Regel eine metrische Zahl, wird diese durch eine geeignete (Response-)Funktion transformiert, z.B. so dass die Zahl im Intervall [0,1] liegt, kann dies als Prognnose für die Wahrscheinlichkeit eine 1 zu beobachten verwendet werden.

3.1 Responsefunktionen

Im Modell einer binären Zielvariable können verschiedene Responsefunktionen gewählt werden. Die logistische, die Probit oder die des komplementären Log-Log Modells. Je nach Wahl kann es zu unterschiedlichen Resultaten kommen, jedoch sind die Unterschiede in der Regel relativ gering.

#Definition der drei Funktionen

#Sigmoid
sigmoid <- function(x){
  return(1/(1 + exp(-x)))
}

#Log-Log
log_log <- function(x){
  return(1 - exp(-exp(x)))
}

#Probit
probit <- function(x){
  return(pnorm(x))
}


x <- seq(-5, 5, length.out = 200)


p <- plot_ly() %>%
  add_lines(x = x, y = sigmoid(x), name = 'sigmoid', line = list(color = the_grey)) %>% 
  add_lines(x = x, y = log_log(x), name = 'log_log', line = list(color = the_green)) %>% 
  add_lines(x = x, y = probit(x), name = 'probit', line = list(color = the_blue)) %>% 
  layout(
    title = 'Responsefunktion binäres Modell',
    xaxis = list(title = TeX('$x$')),
    yaxis = list(title = TeX('$h(x)$'))
  )

p

3.2 Maximum Likelihood Schätzung

Bei der logistischen Regression ist die abhängige Variable binomial- bzw. Bernoulli-verteilt, also \(Y_i \sim B(1, \pi_i)\). Die bedeutet die Likelihood Funktion ist:

\[ \mathcal{L}(\pi) = \prod_{i = 1} f(y_i | \pi ) = \prod_{i = 1} \pi^{y_i} (1 - \pi)^{1 - y_i} \]

Bei der Regression wird an Stelle eines Wertes \(\pi\) für jeden Datenpunkt \(y_i\) ein spezieller Wert \(\pi_i\) eingesetzt, der von den Prediktoren \(\boldsymbol{x_i}\) abhängt. Genauer bestimmt sich:

\[ \pi_i = h(\eta_i) \] wobei \(h\) eine Responsefunktion mit \(f: \mathbb{R}^n \to [0,1]\) ist und \(\eta_i\) der lineare Prediktor \(\boldsymbol{x}_i^T \boldsymbol{\beta}\) mit \(\boldsymbol{x}_i, \boldsymbol{\beta} \in \mathbb{R}^p\) ist.

Die Likelihood Funktion wird somit zu:

\[ \mathcal{L}(\boldsymbol{\beta}) = \prod_{i = 1} f(y_i | \boldsymbol{x}_i , \boldsymbol{\beta} ) = \prod_{i = 1} \pi_i^{y_i} (1 - \pi_i)^{1 - y_i} \]

Simulieren wir uns ein paar Datenpunkte und bestimmen die Log-Likelihood:

#Wir verwenden zwei fiktive Prediktoren mit den Paramterwerten
beta_0 <- 0.2
beta_1 <- 1.1
beta_2 <- -0.8

#Anzahl der Datenpunkte
n <- 100


#Simuliere fiktive deterministische Werte der Prediktoren
set.seed(42)
x_1 <- rnorm(n)
x_2 <- rnorm(n)

#linearer Prediktor
lin_pred <- beta_0 + beta_1 * x_1 + beta_2 * x_2

#Transformationen zu Wahrscheinlichkeiten pi_i
prob <- sigmoid(lin_pred)

#Ziehe Realisierungen der abhängigen Variable
y <- rbinom(n, size = 1, prob = prob)

#So würde ein entsprechender Datensatz aussehen
df <- data.frame(y, x_1, x_2)
head(df)

##   y        x_1        x_2
## 1 1  1.3709584  1.2009654
## 2 0 -0.5646982  1.0447511
## 3 1  0.3631284 -1.0032086
## 4 0  0.6328626  1.8484819
## 5 0  0.4042683 -0.6667734
## 6 0 -0.1061245  0.1055138

#Wir kennen in der Realität die Werte beta_0, beta_1, beta_2 nicht, diese wollen wir durch ML schätzen
#Wie bestimmen wir die Likelihood?
#Raten wir Werte für die Parameter
beta_0_guess <- -0.1
beta_1_guess <-  0.5
beta_2_guess <- -0.7

#Bestimmen wir hierzu die zugehörigen geschätzten Werte für pi_i
prob_hat <- sigmoid(beta_0_guess + beta_1_guess * x_1 + beta_2_guess * x_2)

#Berechnen wir die Log-Likelihood
sum(dbinom(y, size = 1, prob = prob_hat, log = TRUE)) #durch log=TRUE werden direkt die logarithmierten Werte ausgegeben

## [1] -56.37434

Da wir die Parameter nicht raten, sondern durch Maximierung der Likelihood schätzen möchten, können wir dies durch eine Funktion, welche die Log-Likelihood in Abhängigkeit der Parameter ausgibt, bestimmen. Diese Funktion muss bei gegebenen Daten maximiert werden.

#Log-Likelihood-Funktion
log_likel <- function(beta, X, y){
  
  #Werte die beider Schätzung eingesetzt werden
  beta_0_hat <- beta[1]
  beta_1_hat <- beta[2]
  beta_2_hat <- beta[3]
  
  lin_pred_hat <- beta_0_hat + beta_1_hat * X[,1] + beta_2_hat * X[,2]
  
  prob_hat <- sigmoid(lin_pred_hat)
  
  ll <- sum(dbinom(y, size = 1, prob = prob_hat, log = TRUE))
  
  return(ll)
  
}

#Test
log_likel(c(beta_0_guess, beta_1_guess, beta_2_guess), cbind(x_1, x_2), y)

## [1] -56.37434

Um die Likelihood Funktion zu maximieren stehen in R eine Reihe numerischer Verfahren zur Verfügung, z.B. die optim-Funktion bei der verschiedene Methoden der numerischen Optimierung verwendet werden können.

opt <- optim(par = c(0,0,0), fn =  log_likel, control = list(fnscale = -1), X = cbind(x_1, x_2), y = y,
      method = 'BFGS')

opt

## $par
## [1]  0.06302419  1.41290438 -0.89443093
## 
## $value
## [1] -50.28128
## 
## $counts
## function gradient 
##       16        8 
## 
## $convergence
## [1] 0
## 
## $message
## NULL

Unter “par” sind hier die geschätzten Werte zu finden. Zusätzlich können wir den Wert der ersten partiellen Ableitungen, also den Gradienten an den geschätzen Werten überprüfen und uns die Hessematrix bzw. die Fisher-Information an den geschätzten Werten ansehen.

#Der Gradient sollte in etwa bei Null liegen
grad(log_likel, x = opt$par, X = cbind(x_1, x_2), y = y)

## [1] -0.0004202148 -0.0004039046  0.0011287774

#Die Fisher-Information sagt etwas über den Informationsgehalt der Variablen aus.
-hessian(log_likel, x = opt$par, X = cbind(x_1, x_2), y = y)

##            [,1]       [,2]       [,3]
## [1,] 16.6164734  0.2213691 -0.9899491
## [2,]  0.2213691 11.3931078  3.4505959
## [3,] -0.9899491  3.4505959 11.7206068

3.3 Logistische Regression in R

Selbstverständlich bestehen in R bereits Funktionen und Pakete zum Schätzen generalisierter linearer Modelle. Für die logistische Regression können wir die glm-Funktion verwenden. Hierzu müssen wir unter dem Attribut family lediglich angeben, wie die abhängige Variable verteilt ist, da neben der logistischen Regresssion weitere Regressionsmodelle mit der Funktion geschätzt werden können.

glm_fit <- glm(y ~ x_1 + x_2, family = binomial(link = 'logit'), data = df)
summary(glm_fit)

## 
## Call:
## glm(formula = y ~ x_1 + x_2, family = binomial(link = "logit"), 
##     data = df)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.3020  -0.7962   0.3596   0.7581   2.3893  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.06301    0.24612   0.256  0.79795    
## x_1          1.41283    0.31065   4.548 5.42e-06 ***
## x_2         -0.89432    0.30702  -2.913  0.00358 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 137.99  on 99  degrees of freedom
## Residual deviance: 100.56  on 97  degrees of freedom
## AIC: 106.56
## 
## Number of Fisher Scoring iterations: 5

Im Output erkennen wir eine ähnliche Struktur, wie wir es vom Output der linearen Regression kennen. Insbesondere interessieren wir uns für die Schätzungen der Parameter, die in der “Estimates” Spalte zu sehen sind. Weiter können wir dem Output Informationen zur statistischen Signifikanz bei der jeweiligen Nullhypothese von keinem Effekt entnehmen.

3.4 Wer überlebt auf der Titanic

Nun wollen wir ein logistisches Regressionsmodell für echte Daten schätzen. Hierzu verwenden wir einen Datensatz zum Untergang der Titanic, mit dem wir analysieren wollen, ob es gewisse Faktoren gab, die zu einer höheren bzw. geringeren Überlebenswahrscheinlichkeit beigetragen haben. Der Datensatz ist bereits in Trainings- und Testdaten aufgeteilt. Mit den Trainingsdaten schätzen wir das Modell, mit den Testdaten können wir überprüfen, inwiefern es dem Modell gelingt, das Überleben zu prognostizieren.

#Einlesen der Daten
train <- read.csv('Data/titanic/train.csv')
test <- read.csv('Data/titanic/test.csv')

#Betrachten wir die Daten
head(train)

##   PassengerId Survived Pclass
## 1           1        0      3
## 2           2        1      1
## 3           3        1      3
## 4           4        1      1
## 5           5        0      3
## 6           6        0      3
##                                                  Name    Sex Age SibSp Parch
## 1                             Braund, Mr. Owen Harris   male  22     1     0
## 2 Cumings, Mrs. John Bradley (Florence Briggs Thayer) female  38     1     0
## 3                              Heikkinen, Miss. Laina female  26     0     0
## 4        Futrelle, Mrs. Jacques Heath (Lily May Peel) female  35     1     0
## 5                            Allen, Mr. William Henry   male  35     0     0
## 6                                    Moran, Mr. James   male  NA     0     0
##             Ticket    Fare Cabin Embarked
## 1        A/5 21171  7.2500              S
## 2         PC 17599 71.2833   C85        C
## 3 STON/O2. 3101282  7.9250              S
## 4           113803 53.1000  C123        S
## 5           373450  8.0500              S
## 6           330877  8.4583              Q

Eine Beschreibung der Daten kann man hier finden. Bevor man mit der Regressionsanalyse beginnt, sollte man sich den Datensatz etwas näher ansehen und vor allem auch auf Besonderheiten und Unregelmäßigkeiten untersuchen. Z.B. können wir damit beginnen, uns anzusehen, welche Datentypen im Datensatz enthalten sind.

str(train)

## 'data.frame':    891 obs. of  12 variables:
##  $ PassengerId: int  1 2 3 4 5 6 7 8 9 10 ...
##  $ Survived   : int  0 1 1 1 0 0 0 0 1 1 ...
##  $ Pclass     : int  3 1 3 1 3 3 1 3 3 2 ...
##  $ Name       : Factor w/ 891 levels "Abbing, Mr. Anthony",..: 109 191 358 277 16 559 520 629 417 581 ...
##  $ Sex        : Factor w/ 2 levels "female","male": 2 1 1 1 2 2 2 2 1 1 ...
##  $ Age        : num  22 38 26 35 35 NA 54 2 27 14 ...
##  $ SibSp      : int  1 1 0 1 0 0 0 3 0 1 ...
##  $ Parch      : int  0 0 0 0 0 0 0 1 2 0 ...
##  $ Ticket     : Factor w/ 681 levels "110152","110413",..: 524 597 670 50 473 276 86 396 345 133 ...
##  $ Fare       : num  7.25 71.28 7.92 53.1 8.05 ...
##  $ Cabin      : Factor w/ 148 levels "","A10","A14",..: 1 83 1 57 1 1 131 1 1 1 ...
##  $ Embarked   : Factor w/ 4 levels "","C","Q","S": 4 2 4 4 4 3 4 4 4 2 ...

Im Datensatz sind mehrere kategoriale Daten wie die Passagierklasse (Pclass) oder das Geschlecht (Sex) enthalten. Um mit diesen Variablen im Modell arbeiten zu können, bringen wir Sie in Dummykodierung, bei der dann eine Variable als Referenzkategorie verwendet wird.

#Funktion zur automatischen Dummykodierung aus dem fastDummies Pakets
train_df <- dummy_cols(train, select_columns = c('Pclass', 'Sex', 'Embarked'))

Und wir können den Datensatz auf NAs untersuchen. NAs sind Datenfelder mit unbekanntem Wert, diese können entweder durch Schätzungen wie, z.B. den Mittelwert der übrigen Beobachtungen, ersetzt oder für die Schätzung ignoriert werden.

apply(train_df, 2, function(x){return(sum(is.na(x)))})

## PassengerId    Survived      Pclass        Name         Sex         Age 
##           0           0           0           0           0         177 
##       SibSp       Parch      Ticket        Fare       Cabin    Embarked 
##           0           0           0           0           0           0 
##    Pclass_1    Pclass_2    Pclass_3  Sex_female    Sex_male   Embarked_ 
##           0           0           0           0           0           0 
##  Embarked_C  Embarked_Q  Embarked_S 
##           0           0           0

#eine gute Übung wäre diese Zeile nachvollziehen zu können

Da wir später auch in der Lage sein, wollen Vorhersagen erzeugen zu können, selbst wenn uns das Alter nicht bekannt ist, ersetzten wir hier die nicht bekannten Werte. Die Frage ist nur mit was? Einfache Varianten wären z.B.

der Mittlerwert des Alters im Datensatz
der Median des Alters im Datensatz

Dies ist aus statistischer Sicht aber eher schlechte Praxis, da hiermit ein Zusammenhang erzeugt wird, wo keiner ist. Wir verwenden daher eine naive Methode. Wir ziehen bei einem NA zufälllig ein Alter aus den bekannten Beobachtungen. Hierfür gibt es in der Statistik selbstverständlich fortgeschrittenere Ansätze.

train_df$Age[1:10]

##  [1] 22 38 26 35 35 NA 54  2 27 14

set.seed(42)
train_df$Age[is.na(train_df$Age)] <- sample(train_df$Age[!is.na(train_df$Age)], size = sum(is.na(train_df$Age)), replace = TRUE)
train_df$Age[1:10]

##  [1] 22 38 26 35 35 39 54  2 27 14

Wir können uns die Daten auch nach Belieben visualisieren, z.B. die abhängige Variable

p_survival <- plot_ly() %>% 
                add_bars(x = c('überlebt', 'nicht überlebt'), y = c(sum(train_df$Survived), dim(train_df)[1] - sum(train_df$Survived)),
                         marker = list(color = the_grey)) %>% 
                layout(title = 'Ueberlebende auf der Titanic')
p_survival

Oder eine Auswahl der Prediktoren

t1 <- list(
  text = "Verteilung Alter",
  xref = "paper",
  yref = "paper",
  yanchor = "bottom",
  xanchor = "center",
  align = "center",
  x = 0.5,
  y = 1,
  showarrow = FALSE
)

t2 <- list(
  text = "Passagierklassen",
  xref = "paper",
  yref = "paper",
  yanchor = "bottom",
  xanchor = "center",
  align = "center",
  x = 0.5,
  y = 1,
  showarrow = FALSE
)

t3 <- list(
  text = "Geschlecht",
  xref = "paper",
  yref = "paper",
  yanchor = "bottom",
  xanchor = "center",
  align = "center",
  x = 0.5,
  y = 1,
  showarrow = FALSE
)

t4 <- list(
  text = "Fahrpreis",
  xref = "paper",
  yref = "paper",
  yanchor = "bottom",
  xanchor = "center",
  align = "center",
  x = 0.5,
  y = 1,
  showarrow = FALSE
)

p_alter <- plot_ly() %>%
  add_histogram(train_df$Age, name = 'Alter', marker = list(color = the_grey)) %>% 
  layout(annotations = t1)

p_klasse <- plot_ly() %>% 
  add_bars(x = c('Klasse 1', 'Klasse 2', 'Klasse 3'), y = as.vector(apply(train_df[c('Pclass_1', 'Pclass_2', 'Pclass_3')], 2, sum)),
           marker = list(color = the_grey)) %>% 
  layout(annotations = t2)

p_geschlecht <- plot_ly() %>% 
  add_bars(x = c('Frauen', 'Männer'), y = as.vector(apply(train_df[c('Sex_female', 'Sex_male')], 2, sum)),
           marker = list(color = the_grey)) %>% 
  layout(annotations = t3)

p_ticket <- plot_ly() %>%
  add_histogram(train_df$Fare, name = 'Fahrpreis', marker = list(color = the_grey)) %>% 
  layout(annotations = t4)

subplot(p_alter, p_klasse, p_geschlecht, p_ticket, nrows = 2,
        titleX = TRUE, titleY = TRUE, margin = 0.1) %>% layout(showlegend = FALSE)

3.4.1 Modellfit und Prognose

Nun wollen wir das logistische Regressionsmodell schätzen. Hierzu verwenden wir die glm-Funktion. Die Wahl der Referenzkategorien sind hier beliebig gewählt. Als Referenzkategorie für das Geschlecht dient “Sex_male”, für die Passagierklasse “Pclass_3”.

glm_fit <- glm(Survived ~ Sex_female + Age + SibSp + Parch + Fare + Pclass_1 + Pclass_2 + Embarked_C + Embarked_Q + Embarked_S,
    data = train_df, family = binomial(link = 'logit'))
summary(glm_fit)

## 
## Call:
## glm(formula = Survived ~ Sex_female + Age + SibSp + Parch + Fare + 
##     Pclass_1 + Pclass_2 + Embarked_C + Embarked_Q + Embarked_S, 
##     family = binomial(link = "logit"), data = train_df)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.5006  -0.6371  -0.4030   0.6153   2.5205  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  11.178266 615.839326   0.018  0.98552    
## Sex_female    2.706322   0.200035  13.529  < 2e-16 ***
## Age          -0.029915   0.006822  -4.385 1.16e-05 ***
## SibSp        -0.287712   0.109053  -2.638  0.00833 ** 
## Parch        -0.088104   0.117275  -0.751  0.45249    
## Fare          0.002232   0.002455   0.910  0.36308    
## Pclass_1      2.034785   0.292741   6.951 3.63e-12 ***
## Pclass_2      1.211268   0.235336   5.147 2.65e-07 ***
## Embarked_C  -12.151996 615.839252  -0.020  0.98426    
## Embarked_Q  -12.262893 615.839306  -0.020  0.98411    
## Embarked_S  -12.652307 615.839242  -0.021  0.98361    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1186.66  on 890  degrees of freedom
## Residual deviance:  790.16  on 880  degrees of freedom
## AIC: 812.16
## 
## Number of Fisher Scoring iterations: 13

Auf Basis des P-values können wir erkennen, dass z.B. bei einem Signifikanzniveau die Prediktoren “Sex_female, Age, SibSp, Pclass_1, Pclass_2” einen Effekt mit Blick auf die Überlebenswahrscheinlichkeit haben. Beispielsweise ist die Überlebenswahrscheinlichkeit höher, wenn es sich um eine Frau handelt, genauso für jüngere Personen und für Personen, die in einer besseren Passagierklasse waren. Dies scheint intuitiv auch sinnvoll. Aber wie gut könnten wir mit diesem Modell das Überleben prognostierzen. Um einen Wert zu prognostizieren, bestimmen wir bei gegebenen Daten \(\boldsymbol{x}_i\) und geschätzen Werten \(\hat{\boldsymbol{\beta}}\) den geschätzten linearen Prediktor \(\hat{\eta}_i\) und die zugehörige Wahrscheinlichkeit \(\hat{\pi}_i\). Mit der Wahrscheinlichkeit können wir uns festlegen, ob wir an ein Überleben glauben oder nicht. Naiverweise könnten wir an ein Überleben glauben, wenn \(\hat{\pi}_i > 0.50\). Die Überlebensprognose ist also \(\mathbb{1}_{\hat{\pi}_i > 0.50}\) Die Wahrscheinlichkeit müssen wir nicht manuell bestimmen, sondern wir können die generische predict-Funktion verwenden. Als erste Kennzahl für die Güte der Prognose können wir die Accuracy (\(Acc\)) verwenden. Hierbei vergleichen wir, wie gut im Durchschnitt über alle Datenpunkte das Prognostizierte mit dem Beobachteten übereinstimmt.

\[ Acc = \sum_{i = 1}^n \mathbb{1}_{\hat{\pi}_i > 0.50 ~ \& \text{ Person hat überlebt}} \]

#Um die In-Sample Accuracy zu bestimmen, genügen folgende Zeilen
in_sample_predict <- (predict(glm_fit, train_df, type = 'response') > 0.5) * 1
mean(in_sample_predict == train_df$Survived)

## [1] 0.8013468

Dieser Wert erscheint auf den ersten Blick nicht schlecht, jedoch wurde das Modell auf diese Daten geschätzt und interessant ist es, ob es für noch nicht beobachtete Daten auch so gut funkioniert. Probieren wir es aus.

#Die Testdaten bereiten wir auf die selbe Art auf, wie die Trainingsdaten
test_df <- dummy_cols(test, select_columns = c('Pclass', 'Sex', 'Embarked'))

#Wir schauen auch hier auf fehlende Werte
apply(test_df, 2, function(x){return(sum(is.na(x)))})

## PassengerId      Pclass        Name         Sex         Age       SibSp 
##           0           0           0           0          86           0 
##       Parch      Ticket        Fare       Cabin    Embarked    Pclass_1 
##           0           0           1           0           0           0 
##    Pclass_2    Pclass_3  Sex_female    Sex_male  Embarked_C  Embarked_Q 
##           0           0           0           0           0           0 
##  Embarked_S 
##           0

#Wir stellen fest, dass wieder für das Alter Werte fehlen. Zudem fehlt ein Wert für den Preis. 
#Wir verfahren analog wie bei den Trainingsdaten
set.seed(42)
test_df$Age[is.na(test_df$Age)] <- sample(test_df$Age[!is.na(test_df$Age)], size = sum(is.na(test_df$Age)), replace = TRUE)
test_df$Fare[is.na(test_df$Fare)] <- sample(test_df$Fare[!is.na(test_df$Fare)], size = sum(is.na(test_df$Fare)), replace = TRUE)

#Test
apply(test_df, 2, function(x){return(sum(is.na(x)))})

## PassengerId      Pclass        Name         Sex         Age       SibSp 
##           0           0           0           0           0           0 
##       Parch      Ticket        Fare       Cabin    Embarked    Pclass_1 
##           0           0           0           0           0           0 
##    Pclass_2    Pclass_3  Sex_female    Sex_male  Embarked_C  Embarked_Q 
##           0           0           0           0           0           0 
##  Embarked_S 
##           0

#Nun erstellen wir Prognosen
pred_surv <- (predict(glm_fit, test_df, type = 'response') > 0.5) * 1

#diese Prognose speichern wir zusammen mit der PassengerID
df_pred <- data.frame('PassengerId' = test_df$PassengerId, 'Survived' = pred_surv)

write.csv(df_pred, file = 'titanic_survival.csv', row.names = FALSE)

Wenn wir die erzeugte csv-Datei hier hochladen, wird unsere Accuracy für die Testdaten angezeigt. Diese liegt bei \(0.766\)

3.4.2 Weitere Prognosemaße

Wenn wir das Modell lediglich über die Frequenz an richtigen Zuordnungen im Sinne von richtig-richtig und falsch-falsch machen, so unterscheiden wir nicht zwischen den Fehlern, die bei der Prognose gemacht werden können. Folgende Szenarien sind bei der Prognose möglich:

Table 3.1: Mögliche Szenarien bei der Prognose
	Realität Event	Realität kein Event
Vorhersage lautet Event	true positive	false negative
Vorhersage lautet kein Event	false positive	true negative

In unserem Beispiel:

“true positive” bedeutet, dass wir Survival = 1 vorhersagen und die Person überlebt
“false positive” bedeutet, dass wir Survival = 1 vorhersagen und die Person nicht überlebt
“true negative” bedeutet, dass wir Survival = 0 vorhersagen und die Person nicht überlebt
“false negative” bedeutet, dass wir Survival = 0 vorhersagen und die Person überlebt

Offensichtlich sind die Fehler, die wir machen können, unterschiedlich. Während eine “false negative” Prognose eher Grund zur Freude ist, wirkt eine “false positive” Prognose deutlich schlimmer. In diesem Kontext unterscheidet man zwischen der

Sensitivität: Rate der “true positives” in Relation zu allen “positives”, also “true positives” + “false negatives”
Spezifizität: Rate der “true negatives” in Relation zu allen “negatives”, also “true negatives” + “false positives”

#Wir schreiben eine eigene Funktion zur Bestimmung der möglichen Fehler
cm_fun <- function(pred, real){
  
  pred_eval <- apply(cbind(pred, real), 1, function(x)(paste(x[1], x[2])))
  tp <- sum(pred_eval == "1 1")
  tn <- sum(pred_eval == "0 0")
  fp <- sum(pred_eval == "1 0")
  fn <- sum(pred_eval == "0 1")
  
  cm <- rbind(c(tp, fp), c(fn, tn))
  
  colnames(cm) <- c("Real Event", "Real No Event")
  rownames(cm) <- c("Predict Event", "Predict No Event")
  
  return(cm)
                     
}

#Dies nennt man auch "Confusion Matrix"
cm <- cm_fun(in_sample_predict, train_df$Survived)
cm

##                  Real Event Real No Event
## Predict Event           240            75
## Predict No Event        102           474

#Somit ist die Sensitivität
cm[1,1] / sum(cm[,1])

## [1] 0.7017544

#Somit ist die Spezifizität
cm[2,2] / sum(cm[,2])

## [1] 0.863388

Zuletzt stellt sich noch die Frage, ob die bisherige Zuordnung der Prognose:

\[ \text{Survival} = \begin{cases} 1 \text{, wenn } \hat{\pi_i} > 0.5 \\ 0 \text{, wenn } \hat{\pi_i} \leq 0.5 \end{cases} \]

optimal ist bzw. ist sie auf jeden Fall sehr willkürlich. Würden wir beispielsweise als Grenzwert 0.9 wählen, wäre dies wesentlich strenger. Dies sollte die Anzahl der Falsch-Positive-Rate (1-Spezifizität) reduzierenen, gleichermaßen jedoch auch die Sensitivität, also die Richtig-Positive-Rate, da wir tendenziell wenig Survival = 1 Prognosen machen werden. Reduzieren wir den Grenzwert, sollten sich beide Raten erhöhen. Ist dies nicht der Fall deutet es eher auf ein Modell hin, welches nicht gut differenziert, also eher zufällige Prognosen ausgibt. Dieser Zusammenhang wird oft graphisch in der ROC (Receiving Operating Characteristic) Kurve visualisiert.

#Prognostizierte Wahrscheinlichkeiten
score <- predict(glm_fit, train_df, type = 'response')
#Mögliche Grenzen für die Prognose
threshold <- seq(0.99, 0.01, length.out = 100)

#Ergebnisse für Spezifizität und Sensitivität
spec <- c()
sens <- c()

#Bestimme die Kennzahlen für unterschiedliche Grenzen
for(i in 1:length(threshold)){
  
  pred_temp <- (score > threshold[i]) * 1
  cm_temp <- cm_fun(pred_temp, train_df$Survived)
  
  spec <- c(spec, cm_temp[2,2] / sum(cm_temp[,2]))
  sens <- c(sens, cm_temp[1,1] / sum(cm_temp[,1]))
  
}

plot_ly() %>% 
  add_lines(x = 1 - spec, y = sens, text = as.character(threshold), line = list(color = the_green),  name = 'Modell') %>% 
  add_lines(x = c(0,1), y = c(0,1), line = list(color = the_grey, dash = "dash")) %>% 
  layout(title = "ROC Kurve",
         xaxis = list(title = "Falsch-Positiv-Rate"),
         yaxis = list(title = "Richtig-Positiv-Rate"), 
         showlegend = FALSE)

Die Information der ROC Kurve wird auch manchmal durch die Area Under the ROC curve (AUROC) quantifiziert. Hierbei wird die Fläche unterhalb der grünen Kurve bestimmt. Je näher der Wert an 1, umso besser das Modell.