---
title: "Logistische Regression"
author: "Peter Zezula (pzezula@uni-goettingen.de)"
subtitle: "M.Psy.205, Erklärungen und Beispiele"
output:
  html_document:
    theme: paper
    toc: yes
    toc_depth: 2
    toc_float: yes
  pdf_document: default
  #  latex_engine: xelatex
header-includes:
 - \usepackage{comment}
params:
 soln: TRUE
 todo: FALSE
mainfont: Helvetica
linkcolor: blue
editor_options: 
  chunk_output_type: console
---


[Rmd](glm_logistic.Rmd)

`r if(!params$todo) {"\\begin{comment}"}`

https://ww2.coastal.edu/kingw/statistics/R-tutorials/logistic.html

`r if(!params$todo) {"\\end{comment}"}`


Bei der logistischen Regression ist die Reaktionsvariable binär. Es wird über intervallskalierte oder kategoriale Erklärvariablen versucht, das Eintreffen bzw. Nicht-Eintreffen eines Zielzustandes zu erklären
.

# Generalisiertes Lineares Modell

Beim generalisierten linearen Modell kommt eine Link-Funktion für die Reaktionsvariable mit ins Spiel

\[
f(\pi) = \beta_0 + \beta_1x_1 + \beta_2x_2 + ... + \beta_qx_q 
 \]


glm(formula, family=familytype(link=linkfunction), data=)

Family           | Default Link Function
---------------- | ----------------------------------------
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")

Allgemeine Lineare Modelle (ALM) können ausgedrückt werden als GLM mit der Link-Funktion "identity", 
also einer Funktion, die schlicht den Eingabewert zurückgibt.



# Generalisiertes Lineares Modell: Logistische Regression


Die Reaktionsvariable ist binär (dichotome response variable).
Die erklärendenen Variablen sind stetig oder binär (explanatory variables dichotome or continuous).

Als Reaktionsvariable wird nicht die erhobene binäre Variable direkt geschätzt, 
sondern der Logit der Auftretenswahrscheinlichkeit bzw. die logarithmierten odds der Variablen.
Die Wahrscheinlichkeit, der einen oder der anderen Ausprägung anzgehören, kann dann berechnet werden.
Die Logit-Parameter stehen im Exponenten einer e-Funktion. 
Im Exponenten finden wir dann eine "ganz normale" lineare Funktion.
Die Logit-Funktion (S-Form) wird benutzt, um die zu schätzende Größe (Auftretenswahrscheinlichkeit) in bestimmten Grenzen zu halten (0-1).
Die für diesen Zweck wichtigen Parameter können über eine **lineare** Gleichung beschrieben werden, trotz der nicht-linearen Form.


## Formeln

vgl. [https://md.psych.bio.uni-goettingen.de/mv/unit/glm_logistic/log-reg.pdf]

odds:
\[
odds = \frac{p}{1 - p}
\]


In der logistischen Funktion wird nicht die Wahrscheinlichkeit des Auftretens einer 1 (Zutreffen) einer Beobachtung direkt berechnet/geschätzt, sondern die logit davon.

Es findet also eine Link-Funktion Anwendung.


\[
L = logit(p) = log(odds) = log(\frac{p}{1 - p}) 
 \]
 
Umkehrung [Achtung auf den Vorzeichenwechsel im Exponenten]:
 
\[
p = \frac{e^L}{1 + e^L} = \frac{1}{1 + e^{-L}}
 \]

Umformung erfolgt durch Multiplikation mit $e^{−x}$

$e^{−x} * e^{x} = 1$



bei einem Vorhersagewert:
\[
logit(p) = log(\frac{p}{1 - p}) = \beta_0 + \beta_1x + \epsilon 
 \]

ausgedrückt als Wahrscheinlichkeitsverhältnis (odds) der Wahrscheinlichkeit p zur Gegenwahrscheinlichkeit 1-p

\[
odds = \frac{p}{1 - p} = e^{\beta_0 + \beta_1x + \epsilon} = e^{\beta_0} + e^{\beta_1x} + e^{\epsilon}
 \]


bei q Vorhersagevariablen:

\[
logit(\pi) = log(\frac{\pi}{1 - \pi}) = \beta_0 + \beta_1x_1 + \beta_2x_2 + ... + \beta_qx_q = \beta_0 +\sum_{k=1}^{q} {\beta_kx_k} 
 \]


## Wald Test

Die geschätzten Parameter lassen sich mit einerm Wald-Test darauf testen, ob sie einen signifikanten Einfluss auf das Kriterium ausüben. 
Die Wald-Statistik berechnet sich, indem man eine Parameterschätzung durch ihren Standardfehler teilt:

\[
W = \frac{b}{se(b)} ∼ N(0;1)
\]

Der W-Wert ist standardnormalverteilt, weshalb er mit einem einfachen z-Test getestet werden kann.






## Anmerkungen zur Mathematik 

In der Logistischen Regression werden Logits geschätzt. Logit ist der Logarithmus eines Odd. Odd ist das Verhältnis des Eintreffens eines Ereignisses zum Nicht-Eintreffen des Ereignisses, also eine andere Art, Wahrscheinlichkeiten auszudrücken. Bei einem Odd (Eintreffen/Nicht-Eintreffen) von 1:4 oder von 0.25 ist die Wahrscheinlichkeit des Eintreffens = 1/5 oder 0.2 .

`exp(x)` ist eine andere Schreibweise für $e^x$

`log(exp(x))` ist x, log (Logarithmus naturalis) ist die Umkehrfunktion zu exp.



```{r}
log(0)
log(1)

# prob of event 1 is high: 0.99
pe <- .99
# odds
pe / (1 - pe)
# log odds
log(pe/(1-pe))
# library(psych) knows this as command logit
require(psych)
logit(pe)

# fiffty fiffty
pe <- .5
cat(paste('odds', pe / (1-pe)))
cat(paste('logit or log(odds)', logit(pe)))

# low
pe <- .1
cat(paste('odds', pe / (1-pe)))
cat(paste('logit or log(odds)', logit(pe)))

# very low
pe <- .000000001
cat(paste('odds', pe / (1-pe)))
cat(paste('logit or log(odds)', logit(pe)))
# the known linear model is put to the power of base e in logistic regression
odds <- function(x){
  return( x / (1-x))
  }
```


p |1-p |log(p) |odds(p) = p/(1-p) |logit(p) = log(odds(p)) 
:---|:---|:---|:---|:---
.99 | .01 | -0.01 | 99    |  4.59
.90 | .10 | -0.10 |  9    |  2.20
.70 | .30 | -0.36 |  2.33 |  0.85
.50 | .50 | -0.69 |  1    |  0
.30 | .70 | -1.20 |   .43 | -0.84
.10 | .90 | -2.30 |   .11 | -2.21
.01 | .99 | -4.60 |   .01 | -4.60

- $0 \le p \le 1$ und damit 1-p schwankt zwischen 0 und 1.
- $-\infty \le log(p) \le 0$  
- $0 <= odds(p) \le +\infty$
- $-\infty \le logit(p) \le +\infty$

Die Logit-Werte werden negativ, wenn die Wahrscheinlichkeit eines Ereignisses unter .5 fällt. 

Die logistische Regression schätzt die logitWerte mit einer linearen Funktion. 
Für jede Beobachtung können wir über `predict()` einen persönlichen logitWert berechnen.
Wir können die logitWerte in Wahrscheinlichkeiten umrechnen, der Zielgruppe anzugehören, also der mit 1 als Ausprägung.

Dies geht mit `p = 1/(1 + exp(-logitWert))` bzw. mit `p = exp(logitWert)/(1 + exp(logitWert))`.



# Einführendes Beispiel


Quelle der Idee: [tutorial](https://ww2.coastal.edu/kingw/statistics/R-tutorials/logistic.html)

Kann das Entdecken im Visual-Attention-Test durch die Leistungen im Stroop-Test vorhergesagt werden?

Der "visual attention test"

[Simons Lab](https://www.simonslab.com/videos.html)

oder:
[youtube](https://www.youtube.com/watch?v=vJG698U2Mvo)

Zur besseren Veranschaulichung werden die Daten ein wenig frisiert. 
Die Beobachtungen werden verdoppelt und die Interferenz etwas verstärkt.
(Das Vorgehen findet sich ganz unten unter dem Punkt 'Daten generieren'



## Stroop Effect


[Quelle](https://en.wikipedia.org/wiki/Stroop_effect)

In his experiments, Stroop administered several variations of the same test for which three different kinds of stimuli were created. 
In the first one, names of colors appeared in black ink. 
In the second, names of colors appeared in a different ink than the color named. 
Finally in the third one, there were squares of a given color.

<!--  ![stroop-mat](img/stroop_mat_orig_paper.png) -->
![stroop-mat](https://md.psych.bio.uni-goettingen.de/mv/unit/glm_logistic/img/stroop_mat_orig_paper.png)

Original figure 1 from experiment 2 of the article "Studies of interference in serial verbal reactions" (1935) Stroop J.R. 
Journal of Experimental Psychology 18: 643-622; 
the original article that described the Stroop effect in English. 
The diagram shows that it takes less time to name colors of squares (no conflict) 
than saying the ink color of a when what is written a different color (conflict)

<!-- ![stroop-fig](img/stroop-fig1-exp2.jpg) -->
![stroop-fig](https://md.psych.bio.uni-goettingen.de/mv/unit/glm_logistic/img/stroop-fig1-exp2.jpg)

Es sind viele Varianten im Einsatz. 
<!-- ![stroop-mat-en](img/stroop_mat_en.jpg) -->
![stroop-mat-en](https://md.psych.bio.uni-goettingen.de/mv/unit/glm_logistic/img/stroop_mat_en.jpg)


## Zu den Stroop-Scores


Quelle: [https://www.sciencedirect.com/science/article/pii/S0887617703001215]

In Golden’s (1978) method, 45 s are given to read each page of color words (W) printed in black ink, 
color hues (C) printed as XXXX, 
and color hues printed as competing color words (CW) (e.g., ‘red’ printed in blue ink). 
W is the number of color words read in 45 s. 
C is the number of color hues named in 45 s, 
and CW is the number of color hues named (while ignoring the color words) in 45 s. 
Each of these tasks has its own normative values and table look-up, 
but as one can easily see the CW score will be dependent upon the relative strengths of word reading and of color naming.

Golden’s (1978) ipsative interference score is then based upon an assertion 
that the time to read a CW item is an additive function of the time to read a word plus the time to name a color. 
The addition of the time to read a word (45/W) and the time to name a color (45/C) 
gives the algebraically simplified formula of (W×C)/(W+C) for the number of predicted CW items completed in 45 s ( Golden, 1978). 
In using this function, the investigator thus makes a hidden assumption 
that the brain adds word reading to color naming processes to produce the results on the CW card.

```
Code | Erklärung
---  | ---
C    | Anzahl genannter Farben (keine Worte)
W    | Anzahl gelesener Farbworte in schwarz
CW   | Anzahl genannter Farben von Farbworten mit anderer Bedeutung als Druckfarbe (Interferenz)
```


## Inspektion der Korrelations-Matrix

```{r}
dd <- read.delim("https://md.psych.bio.uni-goettingen.de/mv/data/div/stroop_vis_att_doubled.txt")
head(dd)
### a correlation matrix
cor(dd[2:5])
```

Die Stroop Leistungen sind moderat positiv korreliert, aber keine scheint mit der Erkennens-Leistung korreliert zu sein, zumindest nicht bedeutsam.
Insgesamt ist scheint es hier nicht sonderlich viel zu sehen zu geben.




## Modell mit einer erklärenden Variable: Interferenz-Leistung CW


Start mit der CW (Leistung in der Interferenz-Bedingung) als Prädiktor.
Die normale lineare Regression ist nur bedingt nützlich. 
Sie ergibt u. U. vorhergesagte Werte außerhalb des möglichen Bereichs.


```{r}
dd <- read.delim("https://md.psych.bio.uni-goettingen.de/mv/data/div/stroop_vis_att_doubled.txt")
# predict seen by CW, classical lm
m.lm <- lm(seen ~ CW, data=dd)
summary(m.lm)
m.lm$fitted.values[order(m.lm$fitted.values)]
# extreme values leed to impossible probabilities (predicted values)
m.lm$coefficients
# very low but possible CW: we divide minimum CW by 10 
cat(m.lm$coefficients[1] + 0.1 * min(dd$CW) * m.lm$coefficients[2])
# very high but possible CW: we multiply maximum CW by 2
cat(m.lm$coefficients[1] + 2 * max(dd$CW) * m.lm$coefficients[2])
```

Eine Grafik dazu.

```{r}
dd <- read.delim("https://md.psych.bio.uni-goettingen.de/mv/data/div/stroop_vis_att_doubled.txt")
# model to fit
m.1 = lm(seen ~ CW, data=dd)

# base plot
plot(dd$seen ~ dd$CW, ylim=c(-.5,1.5), xlim=c(-5,110))
abline(m.1$coefficients)

require("ggplot2")
# check the fitted values (predicted seen probabilities)
pplot <- ggplot(data=dd, aes(x=CW, y=seen))
pplot + 
  geom_point() + 
  # stat_smooth(method = "lm", se=TRUE) +
  scale_x_continuous(limits = c(1, 140)) +
  geom_abline(intercept = m.1$coefficients[1], slope = m.1$coefficients[2])

# predictions of values below 0 and above 1 are possible
```

Es ist klar erkennbar, dass die Regressionsgerade Vorhersagen über 1 und unter 0 erlaubt, die nicht möglich sind.



## Besser: Logistische Regression

Die Frage, um wieviel besser ein Modell die Reaktionsvariable besser erklärt als ein Vergleichsmodell, 
kann mit der Modell-Chi-Quadrat-Statistik angegangen werden.
Die Differenz zwischen einem und einem anderen Modell wird ausgedrückt.
Von der Log-Likelihood des kleineren Modells bzw. des Null-Modells wird die Log-Likelihood des größeren Modells abgezogen und ergibt den Modell-Vergleichs-Chi^2 -Wert.
Aus Everitt (2010), S. 119: 
  The lack of fit of a logistic regression model can be measured by a term known as the *deviance*, 
  which is essentially the ratio of the likelihoods of the model of interest to the saturated model that fits the data perfectly.


## Summary


```{r}
dd <- read.delim("https://md.psych.bio.uni-goettingen.de/mv/data/div/stroop_vis_att_doubled.txt")
m.1 = glm(seen ~ CW, family=binomial(), data=dd)
summary(m.1)
```

Im Summary:

- Formel wird wiederholt
- "Deviance Residuals" ist ein Maß für die Güte der Anpassung, also für den Model-Fit.
  Zunächst sehen wir eine Beschreibung der Verteilung der Deviance-Residuen der Beobachtungen, die in dieser Berechnung verwendet wurden.
  Wir werden noch sehen, wie der Model-Fit aus der Deviance ermittelt wird.
- Danach sehen wir die Koeffizienten, ihre Standard-Errors, die z-Statistik (die manchmal "Wald-Z-Statistik" genannt wird) die zugehörigen p-Werte.
  "CW" ist auf dem 5% Niveau statistisch signifikant.
  Die Koeffizienten (Estimate) der logistischen Regression geben uns die Veränderung in logit-Werten (log(odds(p/1-p))) der Outcome-Werte, 
  wenn die Vorhersagevariable um eine Einheit erhöht wird.
  - Eine Erhöhung um einen CW-Punkt lässt den Logit-Wert des Erkennens des Gorilla um 0.04 anwachsen.
- Unter der Tabelle der Koeffizienten finden sich Fit-Indices, wie die Null-Deviance, die Deviance der Residuen und das AIC. 
  Mit Hilfe dieser Werte können wir den Model-Fit einschätzen (s. u.)



## Modellvergleich via `anova()`

Das Berechnen der log-likelihood ratio und das explizite Ermitteln und Prüfen des Chi^2 Wertes kann auch durch den Befehl `anova()` erledigt werden.

Wie beim Vergleich von `lm()`-Ergebnisobjekten können zwei genestete glm-Ergebnisobjekte mit dem Befehl `anova()` verglichen werden.
Wird nur ein Modell übergeben, wird das Null-Modell implizit als Vergleichsmodell angenommen.
Die Chi^2 Prüfgröße wird über `test="Chisq"` angefordert.

```{r}
dd <- read.delim("https://md.psych.bio.uni-goettingen.de/mv/data/div/stroop_vis_att_doubled.txt")
# fit null model and model CW
m.0 = glm(seen ~ 1, family=binomial(), data=dd)
m.1 = glm(seen ~ CW, family=binomial(), data=dd)

# check model explicitly using anova()
anova(m.0, m.1, test="Chisq")

# check model (Null-Model is automatically generated, deviance of null-model can be found in model$null.deviance)
m.1$deviance
m.1$null.deviance
anova(m.1, test="Chisq")

```



## Deviances und Modellvergleiche damit

Deviance eines Modells ist die log(likelihood) des Modells multipliziert mit -2
Die Multiplikation mit -2 ist eine Konvention.
Die Deviances tauchen im Summary als "Null deviance" bzw. "Residual deviance" auf.
Der Befehl `logLik()`, der als Parameter ein glm-Ergebnisobjekt erwartet, gibt uns die log(likelihood(Modell)) des Modells zurück.

Unterschiede in der Deviance können benutzt werden, um genestete Modelle gegeneinander zu testen.
Dies erreichen wir, indem wir die Deviances der zu vergleichenden Modelle voneinander abziehen (Chi^2 Wert) 
sowie die Differenz der jeweiligen Modell-Freiheitsgrade benutzen,
um in der Chi^2 Tabelle den zugehörigen p-Wert nachzuschlagen (Funktion: `pchisq()`)



```{r}
dd <- read.delim("https://md.psych.bio.uni-goettingen.de/mv/data/div/stroop_vis_att_doubled.txt")

# we can see the deviances in the model summary
m.0 = glm(seen ~ 1, family=binomial(), data=dd)
summary(m.0)
# predict seen by CW
m.1 = glm(seen ~ CW, family=binomial(), data=dd)
summary(m.1)

# the likelihood of a model is a very small value
exp(logLik(m.1))
# so we usually work with the log of it
logLik(m.1)
# or more common with -2*log likelihood or deviances
# the deviances are the -2LL
logLik(m.0)
-2*logLik(m.0)
logLik(m.1)
-2*logLik(m.1)

# deviances are properties of the model result object
m.0$deviance
m.1$deviance
m.1$null.deviance

m.1.chi <- m.1$null.deviance - m.1$deviance
m.1.chi
# degrees of freedom are attributes of model-object
m.1.chidf <- m.1$df.null - m.1$df.residual
m.1.chidf
# get the probability of this chisquare
chisq.prob <- 1 - pchisq(m.1.chi, m.1.chidf)
chisq.prob
# anova() does this all for us
anova(m.0, m.1, test="Chisq")
```





## Prüfung der Koeffizienten

z-statistic kann ansatzweise interpretiert werden wie der t-Wert im `lm()`.
Allerdings hat er vor allem bei großen b-Koeffizienten eine Tendenz zum Fehler Typ II.
Potenziell bedeutsame erklärende Variablen werden dann nicht identifiziert.

Die z-Statistik und ihre Prüfung ist zu finden im Summary des Modells.

Wald-Statistik: (z-statistic)^2 

```{r}
dd <- read.delim("https://md.psych.bio.uni-goettingen.de/mv/data/div/stroop_vis_att_doubled.txt")
m.1 = glm(seen ~ CW, family=binomial(), data=dd)
summary(m.1)
```



## R^2 bzw. Pseudo-R^2

Vorschläge zu Analogien von R^2 für die Logistische Regression sind:

* Hosmer and Lemeshow R^2
* Cox and Snell R^2
* Nagelkerke R^2


```{r}
dd <- read.delim("https://md.psych.bio.uni-goettingen.de/mv/data/div/stroop_vis_att_doubled.txt")
# define a function to report R^2 
# cf. Field (2012, p. 334)
pseudoR2 <- function(LogModel){
  dev <- LogModel$deviance
  nullDev <- LogModel$null.deviance
  modelN <- length(LogModel$fitted.values)
  R.l <- 1 - dev / nullDev
  R.cs <- 1 - exp( -(nullDev - dev) / modelN)
  R.n <- R.cs / ( 1 - ( exp (- (nullDev / modelN))))
  cat("Pseudo R^2 for logistic regression\n")
  cat("Hosmer and Lemeshow R^2", round(R.l, 3), "\n")
  cat("Cox and Snell R^2      ", round(R.cs, 3), "\n")
  cat("Nagelkerke R^2         ", round(R.n, 3), "\n")
}
# fit the model
m.1 = glm(seen ~ CW, family=binomial(), data=dd)
# get the parameter tests
pseudoR2(m.1)
# package DescTools offers PseudoR2, take care of capital P
require(DescTools)
DescTools::PseudoR2(m.1)
DescTools::PseudoR2(m.1, which="Nagelkerke")
```




## Vorhersage der Werte {#Vorhersage_der_Werte}

Der Koeffizient für CW `m1.coefficients[1]` kann ähnlich interpretiert werden wie im `lm()`. 
Er drückt aus, um wieviel sich der logit der Reaktionsvariablen ändert, wenn sich die erklärende Variable um eine Einheit ändert.
Logit der Reaktionsvariablen ist der natürliche Logarithmus der odds, dass die Reaktionsvariable 1 wird.

<!--
Eine Erhöhung des CW um 1, also eine im Zeitlimit korrekt genannte Farbe mehr erhöht die logit, also die odds-ratio um `r m.1 = glm(seen ~ CW, family=binomial(), data=dd); exp(m.1$coefficients[2])`.
-->

Einsetzen der errechneten Koeffizienten in die Formel der Auftretenswahrscheinlichkeit.

\[
p(y_i) = \frac{1}{1 + e^{-(b_0 + b_1x_{1i})}}
\]




```{r}
dd <- read.delim("https://md.psych.bio.uni-goettingen.de/mv/data/div/stroop_vis_att_doubled.txt")
# model to fit
m.1 = glm(seen ~ CW, family=binomial(), data=dd)
# coefficients
m.1$coefficients

# predict probability to see of first observation
cat(1 /(1 + exp(-(m.1$coefficients[1] + m.1$coefficients[2] * dd$CW[1]))))
cat(dd$seen[1], dd$CW[1], m.1$fitted.values[1])

# very small value of CW
cat(1 /(1 + exp(-(m.1$coefficients[1] + m.1$coefficients[2] * 0.1   * min(dd$CW)))))
cat(1 /(1 + exp(-(m.1$coefficients[1] + m.1$coefficients[2] * 0.001 * min(dd$CW)))))
# very high value of CW
cat(1 /(1 + exp(-(m.1$coefficients[1] + m.1$coefficients[2] * 2 * max(dd$CW)))))
cat(1 /(1 + exp(-(m.1$coefficients[1] + m.1$coefficients[2] * 5 * max(dd$CW)))))

# check the fitted values (predicted probabilities to see the gorilla)
m.1$fitted.values[order(m.1$fitted.values)]

# change in odds, calculated on base of CW = 50 and 51
p.50 <- 1 /(1 + exp(-(m.1$coefficients[1] + m.1$coefficients[2] * 50)))
p.51 <- 1 /(1 + exp(-(m.1$coefficients[1] + m.1$coefficients[2] * 51)))
# odds.50 is p.50 / (1-p.50)
# odds.51 is p.51 / (1-p.51)
# diff.odds is odds.51/odds.50
cat((p.51 / (1-p.51))/(p.50/(1-p.50)))
cat(m.1$coefficients[2])
cat(exp(m.1$coefficients[2]))
# log() as inverse of exp()
log(exp(m.1$coefficients[2]))

# we get
head(m.1$fitted)
head(m.1$fitted.values)
# we can also use predict() and get the probabilities to see the gorilla, but we have to set type=response
head(predict(m.1, type="response"))
# default is type "link" and returns the log(odds)
head(predict(m.1, type="link"))
# calculate p for the first observations using the log(odds)
head(1/(1 + exp(-predict(m.1, type="link"))))

```



## Eine Grafik hierzu

Nur der eigentliche Range der Ausgangsdaten.
Die S-Form ist erkennbar.

```{r}
dd <- read.delim("https://md.psych.bio.uni-goettingen.de/mv/data/div/stroop_vis_att_doubled.txt")
# model to fit
m.1 = glm(seen ~ CW, family=binomial(), data=dd)

# check the fitted values (predicted seen probabilities)
dd.plot <- data.frame(cw.o = dd$CW[order(dd$CW)],
                       cw.p = m.1$fitted.values[order(m.1$fitted.values)])

ggplot(data=dd.plot, aes(x=cw.o, y=cw.p)) + geom_line() + geom_point()

```

Den angepassten Verlauf über einen größeren Bereich

```{r}
dd <- read.delim("https://md.psych.bio.uni-goettingen.de/mv/data/div/stroop_vis_att_doubled.txt")
# model to fit
m.1 = glm(seen ~ CW, family=binomial(), data=dd)

cat(1 /(1 + exp(-(m.1$coefficients[1] + m.1$coefficients[2] * dd$CW[1]))))

dd.plot <- data.frame(
  a.o = -100:200,
  a.p = 1 /(1 + exp(-(m.1$coefficients[1] + m.1$coefficients[2] * -100:200)))
  )
                    
# check the fitted values (predicted seen probabilities)
dd.plot2 <- data.frame(cw.o = dd$CW[order(dd$CW)],
                       cw.p = m.1$fitted.values[order(m.1$fitted.values)])

ggplot(data=dd.plot, aes(x=a.o, y=a.p)) + 
  geom_line() + 
  geom_point(data=dd.plot2, aes(x=cw.o, y=cw.p))

# including the real answers and a cutoff line at p=0.5
ggplot(data=dd.plot[dd.plot$a.o > 0 & dd.plot$a.o < 100,], aes(x=a.o, y=a.p)) + 
  geom_line() + 
  geom_point(data=dd.plot2, aes(x=cw.o, y=cw.p)) +
  geom_point(data=dd, aes(x=CW, y=seen, color=seen)) +
  geom_vline(xintercept = 55, color="red")

```





## Klassifikation der Personen auf Basis der Ergebnisse der Logistischen Regression

Der Befehl `table()` gibt die Möglichkeit, eine Bedingung mit anzugeben, nach der die Tabelle in zutreffend ('TRUE') und nicht zutreffend ('FALSE') aufgeteilt wird.

```{r}
dd <- read.delim("https://md.psych.bio.uni-goettingen.de/mv/data/div/stroop_vis_att_doubled.txt")
# model to fit
m.1 = glm(seen ~ CW, family=binomial(), data=dd)

cutoff = 0.5
table(m.1$y, fitted(m.1) > cutoff)
```


Bei 56 Beobachtungen, die den Gorilla nicht sehen, ist der vorhergesagte Wert unter dem cutoff, der Vergleich Wert > cutoff also FALSE.
Dies wäre eine korrekte Vorhersage des Nicht-Entdeckens.
Für 4 Personen sagen wir das Entdecken vorher, weil der Vorhersagewert über dem Cutoff liegt, die Vorhersage ist aber falsch.
Bei 8 Beobachtungen, die den Gorilla gesehen haben, ist liegt der Vorhersagewert über dem Cutoff. Für diese Personen würden wir das Entdecken vorhersagen.
Allerdings liegen wir bei 30 Beobachtungen falsch mit unserer Vorhersage.


## Mehr als eine erklärende Variable

Der oben erläuterte Grundansatz, logit-Werte durch ein lineares Modell vorherzusagen, funktioniert mit allen anderen Möglichkeiten, die das lineare Modell ansonsten auch bietet.

Alle Basis Leistungsmaße (CW, W, C) in das Modell aufnehmen.

```{r}
dd <- read.delim("https://md.psych.bio.uni-goettingen.de/mv/data/div/stroop_vis_att_doubled.txt")
# explicit null model
m.0 = glm(seen ~ 1, family=binomial(), data=dd)
summary(m.0)
# predict seen by CW
m.1 = glm(seen ~ CW, family=binomial(), data=dd)
summary(m.1)
require(DescTools)
PseudoR2(m.1)

# predict seen by CW and all derived mesasures of inference
m.2 = glm(seen ~ CW + W + C, family=binomial(), data=dd)
summary(m.2)
PseudoR2(m.2)

# check model (Null-Model is automatically generated, deviance of null-model can be found in model$null.deviance)
anova(m.2, test="Chisq")
# check improvement of expanding model1 to model2
anova(m.1, m.2, test="Chisq")
```




## abgeleitete Stroop-Variablen: Inhibition-Scores


CW und die daraus abgeleiteten Interferenz-Werte in das Modell aufnehmen.
Die Differenz der Farb-Benennleistung und der Interferenzleistung CW - C
Der Quotient Farb-Benennleistung und der Interferenzleistung CW / C


```{r}
dd <- read.delim("https://md.psych.bio.uni-goettingen.de/mv/data/div/stroop_vis_att_doubled.txt")
# explicit null model
m.0 = glm(seen ~ 1, family=binomial(), data=dd)
summary(m.0)
# predict seen by CW
m.1 = glm(seen ~ CW, family=binomial(), data=dd)
summary(m.1)

# predict seen by CW and all derived mesasures of inference
m.2 = glm(seen ~ CW + d.cw.c + q.cw.c, family=binomial(), data=dd)
summary(m.2)

# check model (Null-Model is automatically generated, deviance of null-model can be found in model$null.deviance)
anova(m.2, test="Chisq")
# check improvement of expanding model1 to model2
anova(m.1, m.2, test="Chisq")
require(DescTools)
DescTools::PseudoR2(m.1)
DescTools::PseudoR2(m.2)

```





## Interaktionsterme


Es können auch Interaktionsterme mit eingehen.
Alle Rohmaße und die zweifachen Interaktionen in das Modell aufnehmen.


```{r}
dd <- read.delim("https://md.psych.bio.uni-goettingen.de/mv/data/div/stroop_vis_att_doubled.txt")
# explicit null model
m.0 = glm(seen ~ 1, family=binomial(), data=dd)
summary(m.0)
# predict seen by CW, W, and C
m.1 = glm(seen ~ CW + W + C, family=binomial(), data=dd)
summary(m.1)

# predict seen by CW, W, C and the two way interactions
m.2 = glm(seen ~ CW + W + C + CW:W + CW:C + W:C, family=binomial(), data=dd)
summary(m.2)

# check model (Null-Model is automatically generated, deviance of null-model can be found in model$null.deviance)
anova(m.2, test="Chisq")
# check improvement of expanding model1 to model2
anova(m.1, m.2, test="Chisq")

require(DescTools)
DescTools::PseudoR2(m.1)
DescTools::PseudoR2(m.2)

```





## Stepwise

`stepwise()` (wenn Rcmdr vorhanden ist) oder `stepAIC()` funktioniert auch bei logistischer Regression.

```{r}
dd <- read.delim("https://md.psych.bio.uni-goettingen.de/mv/data/div/stroop_vis_att_doubled.txt")
require("MASS")
m.step <- glm(seen ~ W+C+CW, family = binomial, data = dd)
#stepwise(m.step, trace = FALSE)
#Rcmdr::stepwise(m.step, direction="backward/forward", criterion="AIC")
stepAIC(m.step, direction='both', scope = list(upper = ~ W + C + CW, lower = ~ 1))

#detach("package:MASS", unload=TRUE)
```



Auch Interaktionsterme können mit auftauchen.
Die kompakte model formula mit "*" funktioniert nicht. Die Interaktionsterme müssen explizit mit '+' angehängt werden.


```{r}
dd <- read.delim("https://md.psych.bio.uni-goettingen.de/mv/data/div/stroop_vis_att_doubled.txt")
require("MASS")
# run a model with two way interactions included
m.step <- glm(seen ~ W + C + CW + CW:W + CW:C + C:W, family = binomial, data = dd)
summary(m.step)
#stepwise(m.step, trace = FALSE)
#stepwise(m.step, direction="backward/forward", criterion="AIC")
## Rcmdr::stepwise(m.step, direction="backward/forward", criterion="AIC")
stepAIC(m.step, direction='both', scope = list(upper = ~ W + C + CW + CW:W + CW:C + C:W, lower = ~ 1))
```








## Ausflug: Auswertung mit den Original-Daten


[Quelle:](https://ww2.coastal.edu/kingw/statistics/R-tutorials/logistic.html)

Im Original-Beispiel kamen nur halb so viele Probanden wie in den obigen Beispielen zum Einsatz.

```{r}
# read data
dd <- read.delim("https://md.psych.bio.uni-goettingen.de/mv/data/div/stroop_vis_att.txt")
# create inhibition score
dd['d.cw.c'] <- dd$CW - dd$C
dd['q.cw.c'] <- dd$CW / dd$C

# some descriptive impression
require(psych)
describe(dd)

# full model, only raw scores
m.a = glm(seen ~ W * C * CW, family=binomial(logit), data=dd)
summary(m.a)
# check full model (Null-Model is automatically generated)
anova(m.a, test="Chisq")

# full model, only CW and dependent scores
m.a = glm(seen ~ CW * d.cw.c * q.cw.c , family=binomial(logit), data=dd)
summary(m.a)
# check full model (Null-Model is automatically generated)
anova(m.a, test="Chisq")


```



## Ausflug: Daten "verlängern"


Die Originaldaten wurden etwas "frisiert". 
Zu den Originaldaten wurden 'passende' Beobachtungen hinzugefügt um die Power zu erhöhen.
Für Interessierte hier das Vorgehen:


```{r}
# double the observations with a little noise
dd.n <- read.delim("https://md.psych.bio.uni-goettingen.de/mv/data/div/stroop_vis_att.txt")
dd.n$nr <- dd.n$nr + length(dd.n$nr)
dd.n$W  <- round(dd.n$W  + rnorm(length(dd.n$nr), 0, sd(dd.n$W)/10))
dd.n$C  <- round(dd.n$C  + rnorm(length(dd.n$nr), 0, sd(dd.n$C)/10))
dd.n$CW <- round(dd.n$CW + rnorm(length(dd.n$nr), 0, sd(dd.n$CW)/10))

# 'double' the observations
dd <- rbind(read.delim("https://md.psych.bio.uni-goettingen.de/mv/data/div/stroop_vis_att.txt"), dd.n)

# give CW a little more predictive power
blind <- dd$nr[dd$seen == 0]
s.blind <- sample(blind, length(blind)/2)
# dd[s.blind,]
# dd[s.blind, c("nr", "CW", "seen")]
dd$CW[s.blind] <- round(dd$CW[s.blind] - sd(dd$CW) / 1.5)
# dd[s.blind, c("nr","CW", "seen")]

# create inhibition scores
dd['d.cw.c'] <- dd$CW - dd$C
dd['q.cw.c'] <- dd$CW / dd$C

# a descriptive impression
require(psych)
describe(dd)

# save the data
##write.table(dd, "stroop_vis_att_double.txt", quote=F, sep="\t", row.names=F)

```


## Ausflug: Kodierung und Faktorstufen


Je nach Kodierung und Festlegung der Faktorstufen sind die Interpretationen der Koeffizienten verschieden.
Im Prinzip wird nur die erste Faktorstufe als das Event definiert, dessen Eintreten bzw. dessen log(odds) modelliert werden soll.
Ein Vertauschen der beiden Events invertiert lediglich das Vorzeichen der Koeffizienten.


```{r}
require(tidyverse)
dd <- read_csv("https://md.psych.bio.uni-goettingen.de/mv/data/div/stats_lectures.csv")
dd <- dd %>% 
  rowwise() %>% # wir wollen zeilenweise vorgehen
  mutate(stats_joy = mean(c(stats_joy1, stats_joy2, # mittelwert bilden
                            stats_joy3, stats_joy4, 
                            stats_joy5)))
dd <- dd %>% mutate(aon_a = factor(attend_or_not, levels = c("course attended", "course not attended")))
dd <- dd %>% mutate(aon_n = factor(attend_or_not, levels = c("course not attended", "course attended")))

head(dd$aon_a)
head(dd$aon_n)

m.a <- glm(aon_a ~ stats_joy, data = dd, family = binomial())
m.n <- glm(aon_n ~ stats_joy, data = dd, family = binomial())

head(model.matrix(m.a))
head(model.matrix(m.n))

summary(m.a)
summary(m.n)
```




# Beispiele, Aufgaben


## Everitt GHQ


Vollziehen Sie das in [log-reg.pdf](log-reg.pdf) dargestellte Rechenbeispiel für die GHQ-Daten aus [everitt] nach.
Die Daten und das Beispiel entstammen dem Buch [Everitt (2010)][everitt]
In den [Beispielen](../example) findet sich das in Everitt hierzu vorgestellte GHQ-Beispiel.



## UCLA-Zulassung

UCLA (University of California, Los Angeles)

A researcher is interested in how variables, 
such as GRE (Graduate Record Exam scores), 
GPA (grade point average) and 
prestige of the undergraduate institution, effect admission into graduate school. 
The response variable, admit/don't admit, is a binary variable.

[data](https://www.ats.ucla.edu/stat/dataa/binary.csv)


Variable | Bedeutung
---------|-------------------------------------------------------------------------------------------------------------------------------
admit | Zulassung an UCLA (binär)
gre   | graduate record exmainations (Zulassungsprüfungspunkte)
gpa   | grade point average  (in etwa Notendurchschnitt: hier nach descriptives mit Maximum 4 (Bestnote) (A=4, B=3, C=2, D=1, F=0))
rank  | Bewertung der Qualtiät Institution an der man vorher gelernt hat (vier Stufen, 1..4))

Fragen: 

- Wie gut kann die Zulassung zur UCLA durch die Zulassungs-Prüfungspunkte, den Notendurchschnitt und die Qualität der Herkunftsinstitution vorhergesagt werden?
- Welche Prädiktoren und Kombinationen hieraus sind möglich?
- Welche Prädiktoren und Kombinationen hieraus können wir statistisch gegen Zufall absichern?
- Wie können die Koeffizienten interpretiert werden?
- Ist die 4-stufige Klassifikation der Herkunfts-Institution intervallskaliert? 
  Wie können wir vorgehen, wenn wir das nicht als gegeben ansehen?


```{r}
#df <- read.csv("https://www.ats.ucla.edu/stat/dataa/binary.csv")
df <- read.csv("https://md.psych.bio.uni-goettingen.de/mv/data/div/example_ucla.csv")
```


cf [example_ucla.html](example_ucla.html)



## Beispiel Diabetes 


Beispiel [Diabetes](example_diabetes.html)


## glm() with aggregated data

[source](https://stackoverflow.com/questions/26953585/r-logistic-regression-using-frequency-table-cannot-find-correct-pearson-chi-sq)

If we have frequencies we can also fit a glm(). 
We need a column with proportion and a column with total number of observations for every combination of predictor values.

`glm(..., family=binomial, weights=...)`

```{r}
dd.t <- readr::read_table("y     x1        x2     freq
0     No        0      268
0     No        1       14
0    Yes        0      109
0    Yes        1        1
1     No        0       31
1     No        1        6
1    Yes        0       45
1    Yes        1        6
")
total<-with(dd.t,freq[1:4]+freq[5:8])
logit=glm(freq/total~x1+x2, data=dd.t[1:4,], family=binomial, weights=total)
# pearson chi square value
chi_test<-sum(residuals(logit, type = "pearson")^2);

# so the data are finally
dd <- cbind(dd.t[1:4,], total)
logit2=glm(freq/total~x1+x2, data=dd, family=binomial, weights=total)
summary(logit2)
```

[more generic](https://fukamilab.github.io/BIO202/04-B-binary-data.html)

# Links und Referenzen

Idee für das einführende Beispiel: [https://ww2.coastal.edu/kingw/statistics/R-tutorials/logistic.html]

[example](https://ww2.coastal.edu/kingw/statistics/R-tutorials/logistic.html) including gorilla experiment

[Field (2012, p. 312)][field] **Kapitel 8: Logistic Regression**

<!-- this does not appear in html-page. It's used only to arm links and references. -->
[field]: https://www.uk.sagepub.com/dsur/main.htm "Field et al (2012). Discovering Statistics using R"
[everitt]: https://www.crcpress.com/product/isbn/9781439807699 "Brian S. Everitt (2010). Multivariable Modeling and Multivariate Analysis for the Behavioral Sciences"




# Check

Reaktionsvariable ist dichotom bzw. binär.

## Logistische Regression ist ein Mitglied aus der Familie der Verallgemeinerten Linearen Modelle

## Bedeutung der Koeffizienten

- geschätzt werden die logits der Wahrscheinlichkeit der Zielgruppe anzugehören bzw 
  die logarithmierten Chancen log(odds).

- odds lassen sich in Wahrscheinlichkeiten umrechnen, der Zielgruppe anzugehören
  `1/(1 + exp(-predict(m.1, type="link")))`

- Interpretation `exp(model$coefficients)` ist: 
  Zu welcher Änderung der Chance (odds) führt eine Erhöhung des Prädiktors um 1.
  Die Roh-Koeffizienten in `model$coefficients` beziehen sich also auf die logarithmierten odds.


## Unterschied fitted und predict

`fitted` gibt Wahrscheinlichkeiten aus, zur Zielgruppe zu gehören.
`predict` hat default `predict(type='link')` und gibt die Chancen aus (odds).
Mit  `predict(type='response')` erhalten wir direkt die Wahrscheinlichkeiten, wie mit `fitted`.

Beispiel unter 
[glm_logistic.html#Vorhersage_der_Werte](#Vorhersage_der_Werte)


## Modellvergleiche

Wir vergleichen Modelle mit Maximum Liklihood geschätzten Parametern
Deren Model-Fit wird u. a. ausgedrückt durch auf der Likelihood basierende Kenngrößen -2LL bzw. deviances.
Wir benutzen `anova(..., test = "Chisq")` da die Likelihood-Unterschiede Chi^2 verteilt sind.
Dies müssen wir im Befehl `anova()` mit dem Parameter `test = "Chisq"` angeben.


## p und logit Werte

$p(y_i) = \frac{1}{1 + e^{-(b_0 + b_1x_{1i})}}$

```{r}
library(tidyverse)
dd <- read_csv("https://md.psych.bio.uni-goettingen.de/mv/data/div/stats_lectures.csv")
dd <- dd %>% rowwise() %>% mutate(stats_joy = mean(c(stats_joy1, stats_joy2, stats_joy3, stats_joy4, stats_joy5)))
dd <- dd %>% mutate(attend_or_not.f = factor(attend_or_not, levels = c("course not attended", "course attended")))
dd <- dd %>% mutate(att.bin = ifelse(attend_or_not == "course attended", 1, 0))
m.0 <- glm(att.bin ~ 1, data = dd, family = binomial())
m.1 <- glm(att.bin ~ stats_joy, data = dd, family = binomial())
m.1$coefficients
exp(m.1$coefficients)
# we use our model to see the change in the outcome that is produced by 1 step in the predictor
# the difference is expressed in logit(odds) units.
# we define some data
new.d <- data.frame(stats_joy=c(1, 2, 3, 4))
# and predict the values
predict(m.1, newdata = new.d) -> n.d
n.d
n.d[1] - n.d[2]
n.d[2] - n.d[3]
n.d[3] - n.d[4]
# logit(p) or log(odds) or log(p/1-p) change by 1.0519 when the average joy on statistics increases by 1

# psych package offers commands logit() and logistic() to convert probabilities to logits and vice versa
# logistic function (1/(1+exp(-x)) # logits given, probabilities result
# logit function (log(p/(1-p))     # probabilities given, logits result

# so we can calculate the actual probabilities given the logits
ll <- predict(m.1, type="link")
head(ll)
head(psych::logistic(ll))

pp <- 1/(1 + exp(-predict(m.1, type="link")))
head(pp)
pp <- (predict(m.1, type="response"))
head(pp)
head(psych::logit(pp))


```



#### Modellvergleiche, Deviance und Likelihood, PseudoR^2

[vgl](https://md.psych.bio.uni-goettingen.de/mv/unit/glm_logistic/glm_logistic.html#modellvergleich_via_anova())


```{r}
library(tidyverse)
dd <- read_csv("https://md.psych.bio.uni-goettingen.de/mv/data/div/stats_lectures.csv")
dd <- dd %>% rowwise() %>% mutate(stats_joy = mean(c(stats_joy1, stats_joy2, stats_joy3, stats_joy4, stats_joy5)))
dd <- dd %>% mutate(attend_or_not.f = factor(attend_or_not, levels = c("course not attended", "course attended")))
dd <- dd %>% mutate(att.bin = ifelse(attend_or_not == "course attended", 1, 0))
m.0 <- glm(att.bin ~ 1, data = dd, family = binomial())
m.1 <- glm(att.bin ~ stats_joy, data = dd, family = binomial())
# we want to see the likelihood
exp(logLik(m.1))
# this is very small and not easy to handle. So we use log(likelihood) usually.
logLik(m.1)
# or even more common -2LL
-2*logLik(m.1)
# by the way, this is, what we call deviance
m.1$deviance
# ... and it is
sum(residuals(m.1)^2)
# in model comparison we compare it with a model that uses a constant as predictor
# m.0 <- glm(attend_or_not ~ 1, data = stats_lectures, family = binomial())
m.0$deviance
# this is done automatically and saved to an attribute in the model
m.1$null.deviance
# we compare these deviances to see, whether the expanded model is better
anova(m.0, m.1, test="Chisq")
# AIC is just a modified -2LL, where the number of parameters punish the model fit
AIC(m.1)
m.1$deviance + length(m.1$coefficients) * 2
AIC(m.0)
# BIC accounts for sample size on top of AIC
BIC(m.1)
m.1$deviance + length(m.1$coefficients) * log(nrow(dd))
BIC(m.0)
# the smaller AIC or BIC, the better data are fitted by the model when comparing hierarchical models

# package DescTools offers PseudoR2, take care of capital P
DescTools::PseudoR2(m.0, which="Nagelkerke")
DescTools::PseudoR2(m.1, which="Nagelkerke")

# and the real check of correctly identified attending students
cutoff = 0.5
table(dd$att.bin, fitted(m.1) > cutoff)
chisq.test(table(dd$att.bin, fitted(m.1) > cutoff))
```

#### AIC

AIC ist eine Modifikation von -2LL. Bestraft wird für die Aufnahme von Parametern

$AIC = -2LL + k * 2$ k ist die Anzahl der Prädiktoren

`AIC(m.1)`: `r AIC(m.1)`  = `m.1$deviance + 2*2`: `r m.1$deviance + 2*2` 


#### BIC

Gegenüber AIC bezieht BIC auch noch die Stichprobengröße in die Bewertung des Modells mit ein

`BIC = -2LL + k * log(n)` $BIC = -2LL + k * log(n-1)$ n ist die Stichprobengröße

Bei hierarchischen Modellen ist das Modell mit dem geringeren AIC bzw. BIC das bessere Modell, also das Modell, das weniger "unerklärte Varianz" in der Reaktionsvariable übriglässt.




### Graphics static

```{r}
library(tidyverse)
dd <- read_csv("https://md.psych.bio.uni-goettingen.de/mv/data/div/stats_lectures.csv")
dd <- dd %>% mutate(attend.bin = ifelse(attend_or_not == "course attended", 1, 0))
dd <- dd %>% rowwise() %>% 
  mutate(stats_joy = mean(c(stats_joy1, stats_joy2, stats_joy3, stats_joy4, stats_joy5)))
dd <- dd %>% mutate(attend_or_not = factor(attend_or_not, levels = c("course not attended", "course attended")))
dd <- dd %>% mutate(attend.bin = ifelse(attend_or_not == "course attended", 1, 0))
m.1 <- glm(attend.bin ~ stats_joy, data = dd, family = binomial())

# this can be plotted by a 2 dimensional plot
ggplot(dd, aes(x=stats_joy, y=fitted(m.1))) + geom_point()

# complexer models might be visualized by 3d plots
m.4 <- glm(attend.bin ~ stats_joy + temperature + 
                sun_joy + temperature:sun_joy, 
              data = dd, family = binomial())
dd$fitted <- fitted(m.4)

require(rockchalk)
plotPlane(model = m.4, plotx1 = "temperature", plotx2 = "sun_joy")

plotPlane(model = m.4, plotx1 = "temperature", plotx2 = "sun_joy", theta = 30, phi =30)

```


### Graphics interactive

Graph is not rendered when document is knitted, this makes no sense. 
We therefore set `eval = FALSE` in the R chunk header.

```{r eval=FALSE}
librry(tidyverse)
dd <- read_csv("https://md.psych.bio.uni-goettingen.de/mv/data/div/stats_lectures.csv")
dd <- dd %>% mutate(attend.bin = ifelse(attend_or_not == "course attended", 1, 0))
m.i <- glm(attend.bin ~ temperature + sun_joy + temperature:sun_joy, data = dd, family = binomial())

require(rgl)
open3d()
plot3d(x=dd$temperature, y=dd$sun_joy, z=dd$attend.bin, col="red", size=4)

grd <- expand.grid(temperature = sort(unique(dd$temperature)), sun_joy = sort(unique(dd$sun_joy)) )
# predicted values are logits, we have to convert them to probabilities that we want to display
grd$pred <- 1/(1 + exp(-predict(m.i, type="link", newdata = grd)))

persp3d(x=unique(grd[[1]]), y=unique(grd[[2]]), 
              z=matrix(grd[[3]],length(unique(grd[[1]])),length(unique(grd[[2]]))), add=TRUE)

```




### Screencast 

- Logistic regression - check - concepts and models
  [StudIP](https://studip.uni-goettingen.de/plugins.php/mediacastplugin/media/check/557da707f35cbedf27871677d1e8e982/83b70c9395e666199d8f3586fba3f34d/126?v=log_reg-de-2020-06-11.mp4) - 
[ownCloud](https://owncloud.gwdg.de/index.php/s/eu2e65D13ec6ka6)

- Logistic regression - check - possibilities of visualization
  [StudIP](https://studip.uni-goettingen.de/plugins.php/mediacastplugin/media/check/557da707f35cbedf27871677d1e8e982/4b3e83dd8dd2afab8c4087e5e5418b3a/126?v=log_reg-visualization-de-202.mp4) - 
[ownCloud](https://owncloud.gwdg.de/index.php/s/fEPdRJOLuOiekI8)

`r if(!params$todo) {"\\begin{comment}"}`
`r if(!params$todo) {"\\end{comment}"}`




Version: `r format(Sys.time(), '%d %B, %Y %H:%M')`

<!--
aus Mittelbau Workshop
https://www.the-experimenter.com/wiki?dest=r-space/r_regression


todo:

```
# Logistic Regression
# where F is a binary factor and
# x1-x3 are continuous predictors
fit <- glm(F~x1+x2+x3,data=mydata,family=binomial())
summary(fit) # display results
confint(fit) # 95% CI for the coefficients
exp(coef(fit)) # exponentiated coefficients
exp(confint(fit)) # 95% CI for exponentiated coefficients
predict(fit, type="response") # predicted values
residuals(fit, type="deviance") # residuals 
```


-->

