09_Inferenz.qmd

```{r}
#| include: false
library(mosaic)
library(gganimate)
library(tidyr)
library(animation)

set.seed(1896)
theme.fom <- theme_classic(22*1.04)
theme.fom <- theme.fom
theme_set(
  theme.fom  
)

options(OutDec = ',')
```

# Elemente der Inferenzstatistik

## Was Sie lernen &#x1F469;&#x200D;&#x1F3EB; 

- Das Daten zur Schätzung verwendet werden.

- Das Unsicherheiten unvermeidlich sind und was Standardfehler und Konfidenzintervall bedeuten.

- Wie Hypothesen geprüft werden und was der p-Wert ist.

##  Schließen aus Daten :muscle:

Unter 170 Covid-19 Fällen sind [8]{.blue} geimpft und [162]{.violet} nicht. 

Können Sie daraus schließen, dass die Impfung vor der Erkrankung schützt?

- [**Ja**]{.green}

- [**Nein**]{.green}

## Basisratenfehler

Wir müssen wissen, wie hoch der Anteil Geimpfter in der Population ist.

![](img/Base_rate_fallacy_with_vaccines.png){width=70% fig-align="center"}

::: {.aside}
Quelle: <https://de.wikipedia.org/wiki/Prävalenzfehler> auf Basis von <https://twitter.com/MarcRummy/status/1464178903224889345>.
:::

##  Schließen aus Daten -- Forts. :muscle:

In der Studienpopulation gab es genau so viele Geimpfte wie Ungeimpfte. 

Können Sie jetzt daraus schließen, dass die Impfung ursächlich vor der Erkrankung schützt?

- [**Ja**]{.green}

- [**Nein**]{.green}

## RCT

- Randomisiertes Experiment mit $\approx 43660$ Teilnehmenden.

- Multinationale Studie, Teilnehmende über 16 Jahre.

- **Zufällige** (!) Zuordnung zu Impf- oder Placebogruppe, jeweils $\approx 21830$ Personen.

- Für Details der Originalstudie siehe z.B. [https://www.nejm.org/doi/full/10.1056/NEJMoa2034577](https://www.nejm.org/doi/full/10.1056/NEJMoa2034577)

::: {.aside}
*Hinweis*: Aufgrund von Ausschlusskriterien weichen die echten Daten leicht ab.
:::

## Ergebnis

```{r}
#| echo: false
#| out.width: '80%'
#| fig.align: 'center'
ergebnis <- rep(c("I","P"), c(8,162))
plotmatrix <- crossing(x = 1:17, y = 1:10)
plotmatrix <-  plotmatrix |>
  mutate(status = ergebnis) 
impf <- sum(ergebnis=="I")
cols <- c("#0000FF","#DA70D6")
ploterg <- ggplot(plotmatrix, 
                  aes(x = x, y = y, color = status)) +
  geom_tile(color = "white", size = .5, alpha = .1, aes(fill = status)) +
  theme_void() +
  geom_text(size = 8, aes(label = status)) +
  scale_color_manual(values = cols,
                    aesthetics = c("colour", "fill")) +
  labs(title = paste0("Ergebnis: ",impf, " Geimpfte unter den Fällen")) +
  theme(legend.position = "none", plot.title = element_text(hjust = 0.5)) +
  guides(guide = "none")
ploterg
```


## Anzahl und Anteil

- Es gab insgesamt $n=\color{blue}{8}+\color{violet}{162}=170$ Covid-19 Fälle.

- Anzahl aus der [Impfgruppe]{.blue} $\color{blue}{8}$.

- Dann ist der Anteil $\color{orange}p$ der [Impfgruppe]{.blue} unter den Fällen: $$\color{orange}p= \frac{\color{blue}{\text{Fälle in Impfgruppe}}}{\text{Fälle insgesamt}}=\frac{\color{blue}{8}}{\color{blue}{8}+{\color{violet}{162}}}=\frac{\color{blue}{8}}{170}=\color{orange}{`r round(8/170,4)`}$$

## Datengenerierender Prozess

- Wir haben gesehen, dass $\color{blue}{8}$ von $170$ Covid-Fällen aus der [Impfgruppe]{.blue} kamen.

- Aber wie ist es zu dieser Anzahl gekommen? Einige Menschen erkranken (leider), andere nicht.

- Wir stellen uns vor, es gibt eine theoretische Wahrscheinlichkeit $\color{purple}\pi$ (gr.: *pi*), dass in dieser Studie ein Fall aus der [Impfgruppe]{.blue} kommt.


## $\pi$ und $p$

- $\color{purple}{\pi}$ ist der unbekannte, theoretische Anteil im datengenerierenden Prozess. Das ist das was uns *eigentlich* interessiert. Die *Wahrheit*, die wir suchen.

- $\color{orange}{p}$ ist der beobachtete Anteil in unseren Daten.

$$\underbrace{Statistik}_{\color{orange}{p}} = \underbrace{Wert\,Parameter}_{\color{purple}{\pi}} + Verzerrung + Rauschen$$

## Mutmaßlichkeit

Je nachdem, welches $\color{violet}\pi$ vorliegt, desto *mutmaßlicher* (engl.: *likely*) ist eine Anzahl $x=\color{blue}{8}$ bei $n=170$ Fällen:

```{r}
#| out.width: '65%'
#| fig.align: 'center'
pi <- seq(0,1, by = 0.001)
li <- dbinom(8, 170, pi)
gf_line(li ~ pi, linewidth = 2) |>
  gf_labs(x = expression(pi), y = "Likelihood")
```

## Maximum Likelihood

- Je nach dem wie groß das unbekannte, wahre $\color{purple}\pi$ in der Natur ist, desto *mutmaßlicher* ist unsere Beobachtung.

- Als **Punktschätzer** für $\color{purple}\pi$ können wir den Wert verwenden, für den unsere Likelihood maximal ist, hier $\color{orange}p$:

$$\hat{\pi}=p=\frac{8}{170}=`r round(8/170,4)`$$

## $\hat{\pi}$

- $\hat{\pi}$ (sprich: *pi Dach*) ist unser Punktschätzer für $\color{purple}{\pi}$. Wir kennen $\color{purple}{\pi}$ nicht, aber aufgrund unserer Daten würden wir tippen, dass er bei $\color{orange}{p}$ liegt.

```{r}
#| out.width: '65%'
#| fig.align: 'center'
gf_line(li ~ pi, linewidth = 2) |>
  gf_labs(x = expression(pi), y = "Likelihood") |>
  gf_vline(xintercept = 8/170, color = "#FF8811", linewidth = 2)
```


##  Zufall? :muscle:

Wenn die Studie wiederholt werden würde. Muss dann dasselbe Ergebnis, also $\color{orange}{p}=\frac{\color{blue}{8}}{170}$ herauskommen?

- [**Ja**]{.green}

- [**Nein**]{.green}

## Verzerrung und Rauschen

- Die Zuordnung zu Impfung und Placebo erfolgte zufällig. Daher liegt keine (systematische) Verzerrung vor.


$$\underbrace{\color{orange}{`r round(8/170,4)`}}_{Statistik\, \color{orange}{p}} = \underbrace{\color{purple}{\pi}}_{Wert\,Parameter} + \underbrace{0}_{Verzerrung} + Rauschen$$

## Simulierte RCTs

```{r}
#| eval: false
# Einmalig laufen lassen "Run Current Chunk".
# Simulation_Sam.gif dann in den Unterordner img verschieben
simdat <- rbinom(100, size = length(ergebnis), prob = mean(ergebnis=="I"))
SimPlot <- list()
for (i in 1:100) {
  simergebnis <- rep(c("I","P"), c(simdat[i],length(ergebnis)-simdat[i]))
  cols <- c("#0000FF","#DA70D6")
  if (simdat[i]==0) cols <- "#DA70D6"

  plotmatrix <-  plotmatrix |>
    mutate(status = simergebnis) 
  
  SimPlot[[i]] <- ggplot(plotmatrix, 
                         aes(x = x, y = y, color = status)) +
    geom_tile(color = "white", size = .5, alpha = .1, aes(fill = status)) +
    theme_void() +
    geom_text(size = 8, aes(label = status)) +
    scale_color_manual(values = cols,
                       aesthetics = c("colour", "fill")) +
    labs(title = paste0(i, ". Simulation: ", simdat[i], " Geimpfte unter den Fällen")) +
    theme(legend.position = "none", plot.title = element_text(hjust = 0.5)) +
    guides(guide = "none") 
}


saveGIF({for (i in 1:100) plot(SimPlot[[i]])}, 
        movie.name = "Simulation_Sam.gif", interval = 1,
        ani.width = 600, ani.height = 400)
```

![](img/Simulation_Sam.gif){width=95% fig-align="center"}

## Das Problem

:::: {.columns}
::: {.column width="60%"}
- Das vorliegende Studienergebnis ist nur eine der zufälligen Zuordnungen. 

- Das beobachtete $\color{orange}{p}=\color{orange}{`r round(8/170,4)`}$ kann mehr oder weniger vom wahren, aber unbekannten Wert $\color{purple}{\pi}$ abweichen.
:::

::: {.column width="40%"}
![](img/Meme_Unsicherheit.jpg){width=95% fig-align="center"}
:::
::::

## Wichtige Fragen

:::: {.columns}
::: {.column width="50%"}
- Wie präzise ist unsere Schätzung des Parameters $\color{purple}{\pi}$ durch die Statistik $\color{orange}{p}$?

- Wie groß ist unsere Schätzunsicherheit?
:::

::: {.column width="50%"}
![](img/Meme_Sicherheit.jpg){width=80% fig-align="center"}
:::
::::


## Standardfehler

- Um die Präzision der Schätzung, d.h., die Größe des Rauschens abzuschätzen, betrachten wir die Standardabweichung von $p$, den (geschätzten) **Standardfehler** (engl.: *standard error*).

- Er quantifiziert die Streuung des Anteils bei wiederholten zufälligen Daten.
$$se_p=\sqrt{{\frac{p\cdot(1-p)}{n}}}$$

## Standardfehler der Studie

\begin{align*}
\hat{se}_p &= \sqrt{{\frac{\color{orange}{p}\cdot(1-\color{orange}{p})}{n}}} \\
&= \sqrt{{\frac{\color{orange}{\frac{8}{170}}\cdot(1-\color{orange}{\frac{8}{170}})}{170}}}=`r round(sqrt((8/170*(1-8/170))/170),4)`
\end{align*}

## Schätzunsicherheit

- Mit Hilfe des Standarfehlers können wir abschätzen, wie weit $\color{purple}{\pi}$ und $\color{orange}{p}$ üblicherweise voneinander entfernt sind. 

- Häufig wird die **Schätzunsicherheit** wie folgt beschrieben: 
$$p \pm 2 \cdot se_p$$


## Konfidenzintervall

- $p \pm 2 \cdot se_p$ ist eine erste Näherung für das 95%-**Konfidenzintervall**.

- Das Konfidenzintervall ist ein Bereichsschätzer und gibt auf Basis der vorliegenden Daten einen dazu *kompatiblen* Bereich für den (festen, unbekannten) Parameter an.

- Welche Werte für $\color{purple}{\pi}$ sind bei unserem $\color{orange}{p}$ (und $n$) *plausibel*?

::: {.aside}
Für die Impfstudie liegt das exakte Konfidenzintervall bei $[`r confint(binom.test(8,170))$lower`; `r confint(binom.test(8,170))$upper`]$.
:::

##  Konfidenzintervall :muscle:

Bei dieser Studie: Ist ein Anteil von $\color{purple}{\pi}=0,5$ plausibel?

- [**Ja**]{.green}

- [**Nein**]{.green}

## Hypothese

- Die Hypthese *Die Impfung hilft nicht* kann mathematisch als $\color{purple}{\pi}=0,5$ ausgedrückt werden.

- Wenn $\color{purple}{\pi}=0,5$ gelten würde, würden wir um die $170 \cdot \color{purple}{0,5} = 85$ Patient:innen aus der [Impfgruppe]{.blue} erwarten.

- Aber rein zufällig (*Rauschen*) auch mehr oder weniger als $85$. Kleinere Abweichungen sind dabei wahrscheinlicher als große.

## Simulation

```{r}
#| eval: false
# Einmalig laufen lassen "Run Current Chunk".
# Simulation_H0.gif dann in den Unterordner img verschieben

H0Plot <- list()
simdat <- rbinom(100, size = length(ergebnis), prob = 0.5)

for (i in 1:100) {
  simergebnis <- rep(c("I","P"), c(simdat[i],length(ergebnis)-simdat[i]))
  cols <- c("#0000FF","#DA70D6")
  
  plotmatrix <-  plotmatrix |>
    mutate(status = simergebnis) 
  
  H0Plot[[i]] <- ggplot(plotmatrix, 
                         aes(x = x, y = y, color = status)) +
    geom_tile(color = "white", size = .5, alpha = .1, aes(fill = status)) +
    theme_void() +
    geom_text(size = 8, aes(label = status)) +
    scale_color_manual(values = cols,
                       aesthetics = c("colour", "fill")) +
    labs(title = paste0(i, ". Simulation mit \u03C0 = 0,5: ", simdat[i], " Geimpfte unter den Fällen. ")) +
    theme(legend.position = "none", plot.title = element_text(hjust = 0.5)) +
    guides(guide = "none") 
}

saveGIF({for (i in 1:100) plot(H0Plot[[i]])}, 
        movie.name = "Simulation_H0.gif", interval = 1,
        ani.width = 600, ani.height = 400)
```

![](img/Simulation_H0.gif){width=85% fig-align="center"}

## Ergebnis Simulationen

```{r}
#| out.width: '85%'
#| fig.align: 'center'
simdat <- rbinom(1000, size = length(ergebnis), prob = 0.5)

Covid <- data.frame(I = simdat)
yheight <- 75
binwidth <- 0.5
ratio <- 1.5
dotsize <- ratio

dotchart <- ggplot(Covid, aes(x=I), dpi = 800) + 
      geom_dotplot(binwidth=binwidth, method="histodot", dotsize = dotsize, fill="#00998A")
    
    # use coor_fixed(ratio=binwidth*dotsize*max frequency) to setup the right y axis height.
    dotchart <- dotchart + 
      coord_fixed(ratio=binwidth*dotsize*yheight*ratio)
    
        # tweak the theme a little bit
    dotchart <- dotchart + theme(panel.background=element_blank(),
                                 panel.border = element_blank(),
                                 panel.grid.minor = element_blank(),
                                 #plot.margin=unit(c(1,-5,1,-5), "cm"),
                                 axis.line = element_line(colour = "black"),
                                 axis.line.y = element_blank()
    )
    
    # add more tick mark on x axis
    dotchart = dotchart + scale_x_continuous(breaks = seq(0,170,10), limits = c(0,170))
    # add tick mark on y axis to reflect frequencies. Note yheight is max frequency. 
    dotchart = dotchart + scale_y_continuous(limits=c(0, 1/ratio), expand = c(0, 0), 
                                             breaks = seq(0, 1/ratio, 1/yheight*10*1/ratio), 
                                             labels=seq(0,yheight,10))
    dotchart = dotchart + geom_vline(xintercept = 8, color = "#0000FF")
    dotchart = dotchart + annotate("text", x = 30, y = 0.7*1/ratio, 
                                   label = "Tatsächlich beobachtete Fälle")
    dotchart = dotchart + labs(x="Anzahl", y="Häufigkeit", title = "1000 Simulationen mit \u03C0 = 0,5 und n = 170")
    
dotchart
```

##  Interpretation :muscle:

::: {.columns}
:::: {.column width="60%"}
Sind $\color{blue}{8}$ Patient:innen aus der [Impfgruppe]{.blue} bei $n=170$ plausibel, wenn $\color{purple}{\pi}=0,5$ gelten würde?

- [**Ja**]{.green}

- [**Nein**]{.green}
::::

:::: {.column width="40%"}
```{r}
#| out.width: '100%'
#| fig.align: 'center'
dotchart
```
::::
:::

## p-Wert

- Der **p-Wert** beschreibt die Wahrscheinlichkeit, in einer Studie den beobachteten, oder einen noch extremeren Wert zu erhalten, wenn die **Nullhypothese** stimmen würde.

- Die **Nullhypothese** lautet in der Regel, kein Zusammenhang, kein Unterschied, kein Effekt.

$$\underbrace{{\frac{8}{170}}}_{Statistik\, \color{orange}{p}} = \underbrace{0,5}_{Wert\,Parameter\,\color{purple}{\pi_0}} + \underbrace{0}_{Verzerrung} + Rauschen$$

## p-Wert (I/II) :muscle:

::: {.columns}
:::: {.column width="60%"}
Was gilt hier für den p-Wert?

- [**A**]{.green}: Der p-Wert ist klein (<5%)

- [**B**]{.green}: Der p-Wert ist groß (>95%).
::::

:::: {.column width="40%"}
```{r}
#| out.width: '100%'
#| fig.align: 'center'
dotchart
```
::::
:::

## p-Wert (II/II) :muscle:

Wann steigen die Zweifel, dass die Nullhypothese stimmt?

- [**A**]{.green}: Bei einem kleinen p-Wert.

- [**B**]{.green}: Bei einem großen p-Wert.

- [**C**]{.green}: Es gibt keinen Zusammenhang mit dem p-Wert.


## Statistische Signifikanz

- Ist der p-Wert (sehr) klein, z. B. kleiner als $\alpha=5\%$ (gr.: *alpha*), wird das Ergebnis statistisch signifikant genannt.

- Ist der p-Wert nicht (sehr) klein, so wird das Ergebnis als nicht statistisch signifikant bezeichnet.

## Nullhypothese verwerfen

::: {.columns}
:::: {.column width="60%"}
Wenn die beobachten Daten (sehr) unwahrscheinlich sind, wenn die Nullhypothese stimmen würden, dann haben wir mehr als berechtigten Zweifel daran, dass diese stimmt.
::::

:::: {.column width="40%"}
![](img/Meme_H0verwerfen.jpg){fig-align="center" width="85%"}
::::
:::

## Impfstudie

- Der p-Wert ist < 1/1000, daher ist das Ergebnis statistisch signifikant.

- Die Nullhypothese *Die Impfung schützt nicht* wird daher verworfen und auf einen Schutzeffekt der Impfung geschlossen.

## Fallstricke p-Wert

- Der p-Wert misst nicht die Wahrscheinlichkeit, dass die Nullhypothese stimmt.

- Der p-Wert sagt nicht, ob der Effekt groß ist.

- Ein hoher p-Wert bestätigt nicht die Nullhypothese.

## Ablauf Hypothesenprüfung

1. Aufstellen der Nullhypothese.

2. Daten erheben.

3. Daten zusammenfassen.

4. p-Wert bestimmen.

5. Schlussfolgerung.

## Bayes Statistik :nerd_face:

- Wir können berechnen: $Pr(\color{orange}{p} | \color{purple}{\pi})$. Eine *aleatorische* Unsicherheit. Wie wahrscheinlich sind die Daten bei dem Modell?

- Was wir suchen ist (?): $Pr(\color{purple}{\pi} | \color{orange}{p})$. Eine *epistemische* Unsicherheit. Wie wahrscheinlich ist das Modell bei den Daten?

- Was wir dafür brauchen ist der **Satz von Bayes**.

## Zur Erinnerung

::: {.columns}
:::: {.column width="50%"}
<br>
![](img/Meme_pWert.jpg){fig-align="center" width="75%"}
::::

:::: {.column width="50%"}
![](img/Meme_H0akzeptieren.jpg){fig-align="center" width="80%"}
::::
:::

##  Modulziel &#x1F393; :muscle: :pencil:

Grundgedanken von Inferenzstatistiken in medizinischen Studien erfassen.

::: aside
Vergleiche Fragen:

- Was will ich mit Daten machen?

- Was kann ich mit Daten machen?
:::

## Haben Sie gelernt? :eyes:

- Dass Daten zur Schätzung verwendet werden.

- Dass Unsicherheiten unvermeidlich sind und was Standardfehler und Konfidenzintervall bedeuten.

- Wie Hypothesen geprüft werden und was der p-Wert ist.