--- format: pdf: fontsize: 11pt number-sections: true include-in-header: ../MV_Header_Seminar.tex lang: de --- \setcounter{section}{11} # Logistische Regression ## Mathematische Grundlagen Im *Modell der Logistischen Regression* (LR) folgt eine Labelvariable $y \in \left\lbrace 0, 1 \right\rbrace$ mit dem assoziierten Zufallsvektor $x \in \mathbb{R}^m$ einer *Bernoulli-Verteilung* \begin{equation} \label{eq:LR} p(y) = \mbox{Bern}\left( y; \, \frac{1}{1 + \exp(-\tilde{x}^T \beta)} \right) \end{equation} mit dem *Parametervektor* $\beta \in \mathbb{R}^{m+1}$ und dem *erweiterten Featurevektor* \begin{equation} \label{eq:x-tilde} \tilde{x} := \begin{pmatrix} 1 \\ x \end{pmatrix} \in \mathbb{R}^{m+1} \; . \end{equation} Wie in der Vorlesung gezeigt wurde, muss für einen gegebenen Datensatz \begin{equation} \label{eq:D} \mathcal{D} = \left\lbrace (x_1, y_1), \ldots, (x_n, y_n) \right\rbrace \end{equation} der *Maximum-Likelihood-Schätzer* des Parametervektors \vspace{-0.5em} \begin{equation} \label{eq:LR-ML} \begin{split} \hat{\beta} &= \operatorname*{arg\,max}_{\beta} \ln \prod_{i=1}^n p(y_i) = \operatorname*{arg\,max}_{\beta} \sum_{i=1}^n \ln p(y_i) \\ &= \operatorname*{arg\,max}_{\beta} \sum_{i=1}^n \left[ y_i \ln (f(\tilde{x}_i^\mathrm{T} \beta)) + (1 - y_i) \ln (1 - f(\tilde{x}_i^\mathrm{T} \beta)) \right] \; , \end{split} \end{equation} wobei $f(x) = 1/\left( 1 + \exp\left( -x \right) \right)$, mithilfe eines Gradientenverfahrens für die *Log-Likelihood-Funktion* des Modells geschätzt werden. Wenn die Parameter eines LR-Modells einmal geschätzt sind, können diese zur *Inferenz*, d.h. zur Berechnung der bedingten Wahrscheinlichkeit $p(y=1|x)$ für einen potentiell neuen Datenpunkt $x \in \mathbb{R}^m$ verwendet werden: \vspace{-0.5em} \begin{equation} \label{eq:LDA-Inf} p(y = 1|x) = \frac{1}{1 + \exp(-\tilde{x}^T\beta)} \; . \end{equation} Als *Klassifikationsregel* für die LR ergibt sich, dass der potentiell neue Datenpunkt $x$ der Klasse 1 zugeschrieben wird, wenn $p(y=1|x) > 0.5$, und der Klasse 0 zugeschrieben wird, wenn $p(y=1|x) \leq 0.5$. In der vorliegenden Übung wollen wir LR-Modellschätzung via *leave-one-out cross-validation* sowie die Prädiktion von Labels auf Grundlage der geschätzten Modellparameter für einen realen Datensatz nachvollziehen. ## Analyse in R Die Datei `FADE_SAME.csv` enthält den in der ersten Seminarsitzung vorgestellten Datensatz. Erklären Sie die Funktion des folgenden R-Codes: \footnotesize ```{r, echo = T, eval = T, results = 'hide'} # Daten einlesen fname = 'FADE_SAME.csv' # Dateiname D = read.csv(fname) # Dataframe # Datenmatrix extrahieren rows = startsWith(D$subject, 'subA') # Personen aus Studie A cols = c('novelty.FADE','novelty.SAME', # Definition Variablen 'memory.FADE', 'memory.SAME') X = t(as.matrix(D[rows,cols])) # Datenmatrix y = t(as.matrix(D$sex[rows])) # gruppendefinierende Variable y = 0*(y == "male") + 1*(y == "female") # Labelvariable # y = t(as.matrix(D$Abitur[rows])) # gruppendefinierende Variable # y = 0*(y == "no") + 1*(y == "yes") # Labelvariable print(dim(X)) # Überprüfung Datenmatrix print(dim(y)) # Überprüfung Labelvariable print(sum(y)) # Anzahl positiver Fälle ``` \normalsize ## Erste Programmieraufgabe Führen Sie eine Klassifikationsanalyse auf Grundlage der Datenmatrix $X$ und der Labelvariable $y$ durch, indem Sie die LR-Modellparameter mit dem Prinzip der *leave-one-out cross-validation* schätzen und zur Vorhersage des jeweils ausgelassenen Datenpunkts verwenden. Gehen Sie dazu wie folgt vor: - Kopieren Sie den R-Code aus Abschnitt 8.4 des Arbeitsblatts (8) *Prädiktive Modellierung*. - Entfernen Sie im Abschnitt "Vorbereitung Datenanalyse" die Definition der Variable \verb|L|. - Berechnen Sie im Abschnitt "Training" innerhalb der Schleife \verb|n_train| als die Anzahl der Spalten von \verb|x_train|. - Trainieren Sie mit \verb|glm()| ein logistisches Regressionsmodell auf den Trainingsdaten (\verb|x_train|, \verb|y_train|) und extrahieren Sie dessen Parameter aus dem Feld \verb|$coefficients| als eine Matrix mit $m+1$ Zeilen und $1$ Spalte (d.h. einen $(m+1)$-dimensionalen Vektor). - Erzeugen Sie im Abschnitt "Test" innerhalb der Schleife mittels \verb|rbind| den erweiteren Featurevektor $\tilde{x}$ anhand der oben angegebenen Formel. - Berechnen Sie die bedingte Wahrscheinlichkeit $p(y=1|x)$ gemäß der oben angegebenen Formel. Speichern das Ergebnis als die Variable \verb|p_y|. - Modifizieren Sie die Klassifikationsregel derart, dass das prädizierte Label des $i$-ten Datenpunkts 0 ist, wenn $p(y=1|x) \leq 0.5$, und 1 anderenfalls. - Geben Sie nach der Schleife zur Überprüfung die Anzahl positiver Klassifikationen (\verb|y_pred[,2]==1|) aus. Sie sollten folgende Ergebnisse erhalten: \footnotesize ```{r, echo = F, eval = T} # Vorbereitung Datenanalyse m = nrow(X) # Anzahl Features n = length(y) # Anzahl Datenpunkte y_pred = matrix(rep(NaN,n*2), nrow = n) # wahre und prädizierte Label # leave-one-out cross-validation for (i in 1:n) { # Iteration über Datenpunkte # Datensatz partitionieren x_train = as.matrix(X[,-i]) # i-ter Featurevektor nicht im Trainingsdatendatz y_train = as.matrix(y[,-i]) # i-tes Label nicht im Trainingsdatendatz x_test = as.matrix(X[, i]) # i-ter Featurevektor als Testdatenpunkt y_test = as.matrix(y[, i]) # i-ites Label als Testdatenpunkt y_pred[i,1] = y_test # wahres Label des i-ten Datenpunkts # "Training": LR aus Trainingsdatensatz schätzen n_train = ncol(x_train) lr = glm(y_train ~ t(x_train), family = 'binomial') beta_hat = as.matrix(lr$coefficients, nrow = m+1) # "Test": Label des Testdatenpunkts prädizieren x_tilde = rbind(1, x_test) # erweiterter Featurevektor p_y = 1/(1 + exp(-t(x_tilde) %*% beta_hat)) if (p_y <= 0.5) {y_pred[i,2] = 0} # prädiziertes Label else {y_pred[i,2] = 1} } print(sum(y_pred[,2]==1)) ``` \normalsize ## Abbildung in R Wir wollen nun die geschätzten Modellparameter $\hat{\beta}$ sowie den Bernoulli-Parameter $\mu$ als Funktion des linearen Prädiktors $\eta = \tilde{x}^\mathrm{T} \beta$ im Vergleich mit den tatsächlich beobachteten Labels $y$ visualisieren. Erklären Sie dazu den folgenden R-Code und die Abbildung, die er erzeugt: \footnotesize ```{r, echo = T, eval = F} # logistische Regressionsparameter lr = glm(t(y) ~ t(X), family = 'binomial') # IWLS-Parameterlernen beta_hat = as.matrix(lr$coefficients, nrow = m+1) # Parameterschätzer X_tilde = rbind(rep(1,n), X) # erweiterte Featurematrix eta_hat = t(beta_hat) %*% X_tilde # linearer Prädiktor (geschätzt) eta = seq(-5, +5, len = 21) # linearer Prädiktor (w.a.u.) mu = 1/(1+exp(-eta)) # Bernoulliparameter (w.a.u.) # Abbildungsparameter library(latex2exp) par( family = "sans", mfcol = c(1,2), pty = "m", bty = "l", lwd = 1, las = 1, mgp = c(2,1,0), xaxs = "i", yaxs = "i", font.main = 1, cex = 1, cex.main = 1) # geschätzter Parametervektor \beta_hat barplot(as.vector(beta_hat), col = "gray60", border = "black", xlim = c(0, length(beta_hat)+1), ylim = c(-1, +1), ylab = TeX("$\\hat{\\beta}$"), cex.names = 0.6, names.arg = c('intercept', cols), main = TeX("$\\hat{\\eta} = \\tilde{x}^T \\hat{\\beta}$")) # Bernoulli-Parameter \mu plot(eta, mu, type = "b", pch = 16, col = "gray60", xlim = c(min(eta)+1, max(eta)-1), ylim = c(-0.1, +1.1), xlab = TeX("$\\hat{\\eta}$"), ylab = TeX("$y$"), main = TeX("$\\hat{\\mu} = 1/(1 + exp(-\\hat{\\eta}))$")) # geschätzter linearer Prädiktor \eta_hat # und tatsächliche Labels y points(eta_hat[y==0], y[y==0], pch = 16, col = "red") points(eta_hat[y==1], y[y==1], pch = 16, col = "blue") grid() # Legende legend("topleft", c(TeX("$y=0$"), TeX("$y=1$"), TeX("$\\mu$")), pch = c(16,16,16), col = c("red","blue","gray60"), bty = "n", cex = 1, x.intersp = 1, y.intersp = 2) # Speichern dev.copy2pdf( file = "Abbildungen/Logistische_Regression_1.pdf", width = 10, height = 5) ``` \normalsize \vspace{-1.5em} ![Parameterschätzer $\hat{\beta}$ eines logistischen Regressionsmodells mit vier Features (links) und bedingte Wahrscheinlichkeit $p(y=1|x)$ als Funktion des linearen Prädiktors $\hat{\eta} = \tilde{x}^T \hat{\beta}$ (rechts).]("Abbildungen/Logistische_Regression_1.pdf"){#logistische-regression-1 fig-align="center" width=80%} \pagebreak ## Zweite Programmieraufgabe Evaluieren Sie die Klassifikationsperformanz, indem Sie die Sensitivität, Spezifizität und (balancierte) Genauigkeit der Klassifikation berechnen. Gehen Sie dazu wie folgt vor: - Orientieren Sie sich für diese Aufgabe an Abschnitt 8.5 des Arbeitsblatts (8) *Prädiktive Modellierung*. - Berechnen Sie die Einträge der $2 \times 2$ Konfusionsmatrix und speichern Sie die Ergebnisse als \verb|TN|, \verb|FP|, \verb|FN| und \verb|TP|. - Berechnen Sie mithilfe der in Vorlesung (8) *Prädiktive Modellierung* angegebenen Formeln die true positive rate (TPR), true negative rate (TNR) sowie accuracy (ACC) und balanced accuracy (BAC). - Geben Sie die Resultate ihrer Analyse aus. Kommentieren Sie die Klassifikationsperformanz im Hinblick auf den Unterschied zwischen TPR und TNR. Sie sollten folgende Ergebnisse erhalten: \footnotesize ```{r, echo = F, eval = T} # Evaluation der Klassifikationsperformanz TN = sum(y_pred[y_pred[,1] == 0, 2] == 0) # true negatives FP = sum(y_pred[y_pred[,1] == 0, 2] == 1) # false positives FN = sum(y_pred[y_pred[,1] == 1, 2] == 0) # false negatives TP = sum(y_pred[y_pred[,1] == 1, 2] == 1) # true positives TPR = (TP)/(TP+FN) # true positive rate (sensitivity) TNR = (TN)/(TN+FP) # true negative rate (specificity) Acc = (TN+TP)/(TN+FP+FN+TP) # accuracy BAcc = (TPR+TNR)/2 # balanced accuracy # Ausgabe der Ergebnisse cat( "true positive rate (sensitivity) : ", round(TPR, digits = 4), "\ntrue negative rate (specificity) : ", round(TNR, digits = 4), "\nAccuracy : ", round(Acc, digits = 4), "\nBalanced Accuracy : ", round(BAcc, digits = 4), "\n") ``` \normalsize - Kommentieren Sie im ersten Code-Segment in Abschnitt 12.1 diejenigen Zeilen, die die Labelvariable "Geschlecht" extrahieren, und entkommentieren Sie diejenigen Zeilen, die die Labelvariable "Abitur" extrahieren. Führen Sie die Analyse erneut durch. Sie sollten folgende Ergebnisse erhalten: \footnotesize ```{r, echo = F, eval = T} # Daten einlesen fname = 'FADE_SAME.csv' # Dateiname D = read.csv(fname) # Dataframe # Datenmatrix extrahieren rows = startsWith(D$subject, 'subA') # Personen aus Studie A cols = c('novelty.FADE','novelty.SAME', # Definition Variablen 'memory.FADE', 'memory.SAME') X = t(as.matrix(D[rows,cols])) # Datenmatrix y = t(as.matrix(D$Abitur[rows])) # gruppendefinierende Variable y = 0*(y == "no") + 1*(y == "yes") # Labelvariable # Vorbereitung Datenanalyse m = nrow(X) # Anzahl Features n = length(y) # Anzahl Datenpunkte y_pred = matrix(rep(NaN, n*2), nrow = n) # wahre und prädizierte Label # leave-one-out cross-validation for (i in 1:n) { # Iteration über Datenpunkte # Datensatz partitionieren x_train = as.matrix(X[,-i]) # i-ter Featurevektor nicht im Trainingsdatendatz y_train = as.matrix(y[,-i]) # i-tes Label nicht im Trainingsdatendatz x_test = as.matrix(X[, i]) # i-ter Featurevektor als Testdatenpunkt y_test = as.matrix(y[, i]) # i-ites Label als Testdatenpunkt y_pred[i,1] = y_test # wahres Label des i-ten Datenpunkts # "Training": logistische Regression schätzen n_train = ncol(x_train) # Anzahl Trainingsdatenpunkte lr = glm(y_train ~ t(x_train), family = 'binomial') # IWLS-Parameterlernen beta_hat = as.matrix(lr$coefficients, nrow = m+1) # Parameterschätzer # "Test": bedingte Wahrscheinlichkeit berechnen x_test_til = rbind(1, x_test) # \tilde{x} p_y = 1/(1+exp(-t(x_test_til) %*% beta_hat)) # p(y = 1|x) if (p_y <= 0.5) {y_pred[i,2] = 0} # prädiziertes Label else {y_pred[i,2] = 1} } # Evaluation der Klassifikationsperformanz TN = sum(y_pred[y_pred[,1] == 0, 2] == 0) # true negatives FP = sum(y_pred[y_pred[,1] == 0, 2] == 1) # false positives FN = sum(y_pred[y_pred[,1] == 1, 2] == 0) # false negatives TP = sum(y_pred[y_pred[,1] == 1, 2] == 1) # true positives TPR = (TP)/(TP+FN) # true positive rate (sensitivity) TNR = (TN)/(TN+FP) # true negative rate (specificity) Acc = (TN+TP)/(TN+FP+FN+TP) # accuracy BAcc = (TPR+TNR)/2 # balanced accuracy # Ausgabe der Ergebnisse cat( "true positive rate (sensitivity) : ", round(TPR, digits = 3), "\ntrue negative rate (specificity) : ", round(TNR, digits = 3), "\naccuracy : ", round(Acc, digits = 3), "\nbalanced accuracy : ", round(BAcc, digits = 3), "\n") ``` \normalsize \pagebreak ## Lückentext Füllen Sie mit den in der Übung gewonnenen Erkenntnissen den folgenden Lückentext aus und präsentieren Sie die Ergebnisse im Seminar: \vspace{1em} **Lückentext:** Für die Klassifikation einer binären Labelvariable mit den Werten 0 und 1 aus einer $m \times n$ Datenmatrix benutzt die Logistische Regression (LR) einen $(m+1)$-dimensionalen Parametervektor $\beta$. Aus dem Parametervektor ergibt sich der skalare \underline{ lineare Prädiktor } $\eta = \tilde{x}^\mathrm{T} \beta$. Aus dem linearen Prädiktor ergibt sich der \underline{ Erwartungswert } $\mu = 1/\left( 1 + \exp\left( -\eta \right) \right)$. Die Labelvariable wird schließlich als \underline{ Bernoulli }-verteilt mit dem Parameter $\mu$ angenommen. Im vorliegenden Datensatz wurden die Labels aus der Variable \underline{ Geschlecht } gewonnen, wobei "negative Fälle" \underline{ männliche VPen waren } und "positive Fälle" \underline{ weibliche VPen waren }. Die Klassifikation der so erzeugten Labelvariable aus \underline{ 4 } Features erfolgt mittels *leave-one-out cross-validation*. Es ergeben sich eine true positive rate (TPR, "Sensitivität") von \underline{ 0.889 } und eine true negative rate (TNR, "Spezifizität") von \underline{ 0.066 }. Die balanced accuracy der Klassifikation ist mit \underline{ 0.477 } kleiner als die accuracy mit \underline{ 0.552 }, weil \underline{ Geschlecht im Datensatz nicht gleichverteilt ist (m: 106; w: 153) und fast alle Datenpunkte } \underline{ als positive Fälle klassifiziert werden }. ## Mögliche Klausurfrage Präsentieren Sie im Seminar folgende Klausurfrage und erklären Sie die richtige Antwort: \vspace{1em} **Frage:** Welche Funktion ist die Mean-Funktion im Modell der Logistischen Regression? \begin{enumerate}[a)] \item $f(\eta) = \eta$ \item $f(\eta) = \ln \left( \frac{\eta}{1-\eta} \right)$ \item $f(\eta) = \frac{1}{1 + \exp\left( -\eta \right)}$ $\checkmark$ \item $f(\eta) = (x-\eta)^2$ \end{enumerate} ## Kinderwitz Wohin geht ein Zyklop mit Sehproblemen? \begin{turn}{180} Antwort: Zum Augearzt. \end{turn}