52 Modellieren in R

Letzte Änderung am 29. June 2025 um 19:53:14

“Ich weiß nicht weiter; Ich will mich verändern; Doch wie fang ich’s an?” — Tocotronic, Die Unendlichkeit

Stand des Kapitels: Stand-By (seit 01.2025)

Dieses Kapitel ist auf Stand-By gesetzt, da ich die Thematik des Kapitels aktuell nicht in meiner Lehre oder der statistischen Beratung benötige. Archivierte Kapitel werden nicht von mir weiter gepflegt oder ergänzt. Auftretende Fehler werden aber natürlich beseitigt, wenn die Fehler mir auffallen.

Dieses Startkapitel gibt dir nochmal eine Übersicht über das statistische Modellieren in R. Hier liegt vor allem der Fokus auf R. Es gibt eben eine Reihe von zusätzlichen Pakten, die es dir erlauben noch mehr aus einem statistischen Modell rauszuholen. Am Ende wurde es mir dann aber zu detailliert alle Pakete in jedem Kapitel vorzustellen und anzuwenden. Das ist auch nicht immer sinnig. Häufig willst du erstmal nur das Modell rechnen. Später kannst du dann noch tiefer ins Detail gehen oder aber komplexere Verfahren nutzen. Ich tue mich also etwas schwer dieses Kapitel einzuordnen. Entweder packen wir es ans Ende vom statistischen Modellieren und schauen, dann wie wir alles in R machen. Das steht aber etwas der Intuition entgegen, dass wir in jedem Kapitel zum statistischen Modellieren ja schon was selber machen wollen. In R gibt es dazu dann noch sehr gute Pakete, die das Modellieren sehr viel einfacher machen, dabei dann aber auch für den Anfänger etwas komplexer sind. Ich habe mich daher entschieden, diese aktuelle und komplexere Modellierung einmal hier am Anfang vorzustellen und in den folgenden Kapiteln teilweise darauf zu verweisen, wenn ich es sinnig fand. Du kannst alle Modelle auf althergebrachte Art und Weise rechnen ohne was zu verpassen. Aber manchmal möchte man dann auch effizienter Modellieren. Dafür ist dann dieses Kapitel da: Eine erweiterte Idee von der statistischen Modellierung zu erlangen. Fangen wir also erstmal mit der naheliegenden Frage an.

Was modellieren wir eigentlich?: Wir modellieren die Varianz in den Daten oder anders formuliert, wir versuchen die Varianz in den Daten den jeweiligen Quellen zuzuordnen und darüber dann neue Erkenntnisse zu erlangen.

Was ist ein Modell

Statistical Thinking for the 21st Century — What is a model?

Models are about what changes, and what doesn’t

Als nächstes wollen wir uns die Frage stellen, was ist eigentlich das ziel des Modellierens? Wir wollen ja mit der Modellierung der Varianz irgendwas erreichen. In der Abbildung 52.1 siehst du einmal die drei großen Fragefelder, die wir mit einer Modellierung bearbeiten können.

Abbildung 52.1— Verschiedene Ziele und Möglichkeiten des statistischen Modellierens. Grob können die Möglichkeiten in drei große thematische Zusammenhänge eingeteilt werden. **(A)** *Kausales Modell* – Wie verändert sich $y$, wenn sich $x_1$ ändert? Wie ist der numerische Zusammenhang zwischen $y$ und $x$? **(B)** *Prädikitives Modell* – Wenn $y$ und $x$ gemessen wurden, wie sehen dann die Werte von $y$ für neue $x$-Werte aus? Können wir mit $x$ die Werte in $y$ vorhersagen? **(C)** *Clusteranalyse* – Haben wir in unseren Variablen $x_1$ und $x_2$ eine oder mehrere unbekannte Cluster vorliegen? Wie gehören die Beobachtungen gegeben $x_1$ und $x_2$ zusammen? *[Zum Vergrößern anklicken]*

In diesem Kapitel wollen wir uns auch mit dem Thema marginale Effektmodelle (eng. marginal effect models) beschäftigen. Auch hier ist der deutsche Begriff nicht gebräuchlich, so dass ich mich hier dann immer auf die marginal effect models beziehen werde. Ein zweideutiger Aspekt der marginal effect models besteht darin, dass das Wort “marginal” auf zwei verschiedene und etwas entgegengesetzte Arten verwendet wird:

Bei “marginalen Effekten” (eng. marginal effects) beziehen wir uns auf die Auswirkung einer winzigen (marginalen) Änderung des Regressors auf das Ergebnis. Dies ist eine Steigung oder eine Ableitung.
Der Begriff “marginaler Mittelwert” (eng. marginal means) bezieht sich auf den Prozess der Marginalisierung über die Zeilen eines Vorhersagerasters. Dies ist ein Durchschnitt oder ein Integral.

Mehr dazu dann in dem Abschnitt zu den marginal effect models und dem R Paket {marginaleffects}.

Dafür verweise ich auf Heiss (2022) und das Onlinetutorium Marginalia: A guide to figuring out what the heck marginal effects, marginal slopes, average marginal effects, marginal effects at the mean, and all these other marginal things are. Heiss erklärt dort sehr schön die Zusammenhänge. Dazu dann noch der Verweis auf die Webseite Model to Meaning von Arel-Bundock et al. (2024) um nochmal tiefer in das Modellieren von Daten einzusteigen.

52.1 Genutzte R Pakete

Wir wollen folgende R Pakete in diesem Kapitel nutzen.

R Code [zeigen / verbergen]

pacman::p_load(tidyverse, emmeans, multcomp, conflicted)
conflicts_prefer(dplyr::select)

An der Seite des Kapitels findest du den Link Quellcode anzeigen, über den du Zugang zum gesamten R-Code dieses Kapitels erhältst.

52.2 Theoretisches Modellieren in R

52.2.1 R Pakete

Neben den R Paketen, die wir in den jeweiligen Kapiteln brauchen, kommen noch folgende R Pakete immer wieder dran. Deshalb sind die R Pakete hier schon mal mit den jeweiligen Internetseiten aufgeführt.

Das Buch Tidy Modeling with R gibt nochmal einen tieferen Einblick in das Modellieren in R. Wir immer, es ist ein Vorschlag aber kein Muss.
Das R Paket {parameters} nutzen wir um die Parameter eines Modells aus den Fits der Modelle zu extrahieren. Teilweise sind die Standardausgaben der Funktionen sehr unübersichtich. Hier hilft das R Paket.
Das R Paket {performance} hilft uns zu verstehen, ob die Modelle, die wir gefittet haben, auch funktioniert haben. In einen mathematischen Algorithmus können wir alles reinstecken, fast immer kommt eine Zahl wieder raus.
Das R Paket {tidymodels} nutzen wir als das R Paket um mit Modellen umgehen zu können und eine Vorhersage neuer Daten zu berechnen. Das Paket {tidymodels} ist wie das Paket {tidyverse} eine Sammlung an anderen R Paketen, die wir brauchen werden.

52.2.2 Quantity and tests

Kontrafaktisches Konditional

“Wenn … der Fall wäre, dann wäre …”

Counterfactual comparisons

52.2.3 Generalisierung von `lm()` zu `glm()` und `[g]lmer()`

Die Funktion lm() nutzen wir, wenn das Outcome $y$ einer Normalverteilung folgt.
Die Funktion glm() nutzen wir, wenn das Outcome $y$ einer andere Verteilung folgt.
Die Funktion lmer() nutzen wir, wenn das Outcome $y$ einer Normalverteilung folgt und wir noch einen Block- oder Clusterfaktor vorliegen haben.
Die Funktion glmer() nutzen wir, wenn das Outcome $y$ einer andere Verteilung folgt und wir noch einen Block- oder Clusterfaktor vorliegen haben.

Abbildung 52.2— Übersicht der Namen der Funktionen in R für das lm(), glm() und glmer().

In Abbildung 52.3 sehen wir wie wir den Namen einer Regression bilden. Zuerst entscheiden wir, ob wir nur ein $x$ haben oder mehrere. Mit einem $x$ sprechen wir von einem simplen Modell, wenn wir mehrere $x$ haben wir ein multiples Modell. Im nächsten Schritt benennen wir die Verteilung für das Outcome $y$. Dann müssen wir noch entscheiden, ob wir ein gemischtes Modell vorliegen haben, dann schreiben wir das hin. Sonst lassen wir den Punkt leer. Anschließend kommt noch lineares Modell hinten ran.

Abbildung 52.3— Wie bilde ich den Namen einer Regression? Erst beschreiben wir das $x$, dann das $y$. Am Ende müssen wir noch sagen, ob wir ein gemischtes Modell vorliegen haben oder nicht.

Das R Paket {tidymodels}

Tidy Modeling with R

Veraltet aber manchmal ganz nützlich das R Paket {modelr}

52.2.4 Konfidenzintervall vs. Prädiktionsintervall

Ein Konfidenzintervall gibt den Wertebereich für einen gesuchten Parameter der Grundgesamtheit mit einer bestimmten Wahrscheinlichkeit an. Ein Prognoseintervall gibt den Wertebereich für einen individuellen, zukünftig zu beobachtenden Wert mit einer bestimmten Wahrscheinlichkeit an.

Konfidenzintervall

Text

Konfidenzintervall: Ein Konfidenzintervall gibt den Wertebereich für einen gesuchten Parameter der Grundgesamtheit mit einer bestimmten Wahrscheinlichkeit an.

Prädiktionsintervall

Text

Prädiktionsintervall: Ein Prädiktionsintervall (auch Vorhersageintervall oder Prognoseintervall) gibt den Wertebereich für einen individuellen, zukünftig zu beobachtenden Wert mit einer bestimmten Wahrscheinlichkeit an.

Quantile Regression Forests for Prediction Intervals

The difference between prediction intervals and confidence intervals

P-values for prediction intervals machen keinen Sinn

How NASA didn’t discover the hole in the ozone layer

52.3 Praktisches Modellieren in R

52.3.1 Marginal Effect Models

R Paket {ggeffects}

R Paket {modelbased}

R Paket {marginaleffects}

52.3.2 Visualisierung

R Paket {ggeffects}

R Paket {gtsummary}

Display univariate regression model results in table

52.3.3 …mit `{purrr}`

map() für mehrere Endpunkte wie bei den Spinatdaten in Application gezeigt.

map2() mit Familie für glm zeigen, wenn wir unterschiedliche Outcomes haben.

R Code [zeigen / verbergen]

cutting_tbl <- read_excel("data/multiple_outcomes.xlsx") |> 
  mutate(trt = as_factor(trt),
         block = as_factor(block)) |> 
  mutate_if(is.numeric, round, 2)

Tabelle 52.1— Auszug aus den Daten .

trt	block	shoot	leaf	flower	fruit	ca	drymatter	freshweight	height
control	1	0	136	33	10	0.41	11.59	85.37	25
control	2	0	157	33	9	0.41	10.43	79.78	22
control	1	1	65	10	3	0.4	11.97	78.35	22
control	2	0	88	16	4	0.39	11.24	85.55	27
…	…	…	…	…	…	…	…	…	…
nodium_5th	2	0	137	20	3	0.47	11.11	86.93	28
nodium_5th	1	0	137	23	3	0.43	11.95	89.68	27
nodium_5th	2	1	114	17	9	0.39	11.21	70.82	24
nodium_5th	2	0	192	42	11	0.42	11.29	66.23	22

R Code [zeigen / verbergen]

family_lst <- lst(shoot = binomial(),
                  leaf = quasipoisson(),
                  flower = quasipoisson(),
                  fruit = quasipoisson(),
                  ca = gaussian(), 
                  drymatter = gaussian(), 
                  freshweight = gaussian(), 
                  height = gaussian())

R Code [zeigen / verbergen]

cutting_long_tbl <- cutting_tbl |>
  pivot_longer(cols = shoot:last_col(),
               names_to = "outcome",
               values_to = "rsp") |>
  arrange(outcome, trt, block)

R Code [zeigen / verbergen]

cutting_lst <- cutting_long_tbl |>
  split(~outcome)

R Code [zeigen / verbergen]

glm_lst <- cutting_lst %>% 
  map2(family_lst, 
       ~glm(rsp ~ trt + block + trt:block, 
            data = .x, family = .y))  
  
glm_lst %>% 
  map(pluck, "family")

$ca

Family: binomial 
Link function: logit 


$drymatter

Family: quasipoisson 
Link function: log 


$flower

Family: quasipoisson 
Link function: log 


$freshweight

Family: quasipoisson 
Link function: log 


$fruit

Family: gaussian 
Link function: identity 


$height

Family: gaussian 
Link function: identity 


$leaf

Family: gaussian 
Link function: identity 


$shoot

Family: gaussian 
Link function: identity

R Code [zeigen / verbergen]

glm_lst %>% 
  map(car::Anova)

$ca
Analysis of Deviance Table (Type II tests)

Response: rsp
          LR Chisq Df Pr(>Chisq)
trt       0.222384  3     0.9739
block     0.003979  1     0.9497
trt:block 0.005906  3     0.9999

$drymatter
Analysis of Deviance Table (Type II tests)

Response: rsp
          LR Chisq Df Pr(>Chisq)
trt         5.1410  3     0.1618
block       0.2084  1     0.6480
trt:block   2.0862  3     0.5547

$flower
Analysis of Deviance Table (Type II tests)

Response: rsp
          LR Chisq Df Pr(>Chisq)    
trt        19.3416  3  0.0002323 ***
block       0.5678  1  0.4511477    
trt:block   1.7967  3  0.6156522    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

$freshweight
Analysis of Deviance Table (Type II tests)

Response: rsp
          LR Chisq Df Pr(>Chisq)  
trt         3.6454  3     0.3024  
block       4.0375  1     0.0445 *
trt:block   4.9674  3     0.1742  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

$fruit
Analysis of Deviance Table (Type II tests)

Response: rsp
          LR Chisq Df Pr(>Chisq)  
trt         9.9860  3    0.01869 *
block       0.8118  1    0.36759  
trt:block   1.7163  3    0.63332  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

$height
Analysis of Deviance Table (Type II tests)

Response: rsp
          LR Chisq Df Pr(>Chisq)
trt        2.83465  3     0.4178
block      0.12598  1     0.7226
trt:block  0.94488  3     0.8146

$leaf
Analysis of Deviance Table (Type II tests)

Response: rsp
          LR Chisq Df Pr(>Chisq)   
trt        15.5183  3   0.001423 **
block       1.2623  1   0.261207   
trt:block   1.6400  3   0.650353   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

$shoot
Analysis of Deviance Table (Type II tests)

Response: rsp
          LR Chisq Df Pr(>Chisq)    
trt             17  3  0.0007067 ***
block            0  1  1.0000000    
trt:block        2  3  0.5724067    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R Code [zeigen / verbergen]

emm_lst <- glm_lst %>% 
  map(~emmeans(.x, specs = ~ trt, type = "response")) 

emm_lst %>% 
  map(~contrast(.x, method = "pairwise", adjust = "bonferroni")) %>% 
  map(as_tibble) %>% 
  bind_rows(.id = "outcome")

# A tibble: 48 × 11
   outcome   contrast       odds.ratio     SE    df  null z.ratio p.value  ratio
   <chr>     <chr>               <dbl>  <dbl> <dbl> <dbl>   <dbl>   <dbl>  <dbl>
 1 ca        control / str…      0.829 0.978    Inf     1 -0.159    1     NA    
 2 ca        control / nod…      1.30  1.59     Inf     1  0.219    1     NA    
 3 ca        control / nod…      0.780 0.918    Inf     1 -0.211    1     NA    
 4 ca        strong / nodi…      1.57  1.90     Inf     1  0.376    1     NA    
 5 ca        strong / nodi…      0.941 1.10     Inf     1 -0.0523   1     NA    
 6 ca        nodium_3rd / …      0.598 0.719    Inf     1 -0.428    1     NA    
 7 drymatter control / str…     NA     0.0436   Inf     1 -0.425    1      0.981
 8 drymatter control / nod…     NA     0.0479   Inf     1  1.24     1      1.06 
 9 drymatter control / nod…     NA     0.0424   Inf     1 -0.930    1      0.960
10 drymatter strong / nodi…     NA     0.0486   Inf     1  1.67     0.571  1.08 
# ℹ 38 more rows
# ℹ 2 more variables: estimate <dbl>, t.ratio <dbl>

R Code [zeigen / verbergen]

emm_lst %>% 
  map(~cld(.x, Letters = letters, adjust = "bonferroni")) %>% 
  map(as_tibble) %>% 
  bind_rows(.id = "outcome") %>% 
  select(outcome, trt, .group)  %>% 
  print(n = 28)

# A tibble: 32 × 3
   outcome     trt        .group
   <chr>       <fct>      <chr> 
 1 ca          nodium_3rd " a"  
 2 ca          control    " a"  
 3 ca          strong     " a"  
 4 ca          nodium_5th " a"  
 5 drymatter   nodium_3rd " a"  
 6 drymatter   control    " a"  
 7 drymatter   strong     " a"  
 8 drymatter   nodium_5th " a"  
 9 flower      strong     " a " 
10 flower      nodium_3rd " ab" 
11 flower      nodium_5th "  b" 
12 flower      control    "  b" 
13 freshweight strong     " a"  
14 freshweight control    " a"  
15 freshweight nodium_3rd " a"  
16 freshweight nodium_5th " a"  
17 fruit       strong     " a"  
18 fruit       control    " a"  
19 fruit       nodium_3rd " a"  
20 fruit       nodium_5th " a"  
21 height      strong     " a"  
22 height      control    " a"  
23 height      nodium_3rd " a"  
24 height      nodium_5th " a"  
25 leaf        strong     " a " 
26 leaf        nodium_3rd " ab" 
27 leaf        nodium_5th " ab" 
28 leaf        control    "  b" 
# ℹ 4 more rows

52.4 Verschiebebahnhof

R Paket {parsnip}

in diesem Kapitel wollen Fokus auf {parsnip} und dann als Kasten in allen anderen möglichen Modellierungen ergänzen?

Tidymodels, interactions and anova - a short tutorial

XX Regression in {parsnip}

Auch hier können wir die XX Regression in dem R Paket {parsnip} realisieren. Ein Vorteil von {parsnip} ist, dass wir die Funktionen sehr gut mit dem |>-Operator nutzen können. Das ist bei den Funktionen glm() und lm() in der ursprünglichen Implementierung ja leider nur etwas umständlicher möglich. Deshalb hier einmal die bessere Implementierung.

{parsnip}
{base}

R Code [zeigen / verbergen]

linreg_reg_fit <- linear_reg() |> 
  set_engine("lm") |> 
  fit(jump_length ~ grp, data = simple_tbl)

Dann die ANOVA

R Code [zeigen / verbergen]

linreg_reg_fit|>  
  extract_fit_engine() |> 
  anova()

Analysis of Variance Table

Response: jump_length
          Df Sum Sq Mean Sq F value  Pr(>F)  
grp        1 1.5309 1.53089  9.2541 0.01879 *
Residuals  7 1.1580 0.16543                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R Code [zeigen / verbergen]

linreg_reg_fit <- lm(jump_length ~ grp, data = simple_tbl)

R Code [zeigen / verbergen]

linreg_reg_fit|>  
  anova()

Analysis of Variance Table

Response: jump_length
          Df Sum Sq Mean Sq F value  Pr(>F)  
grp        1 1.5309 1.53089  9.2541 0.01879 *
Residuals  7 1.1580 0.16543                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Eine detailliertere Einführung mit mehr Beispielen für die Nutzung vom R Paket {parsnip} findest du im Kapitel Modellieren in R. Hier soll es dann bei der kurzen Gegenüberstellung bleiben.

Referenzen

Arel-Bundock, V., Greifer, N., & Heiss, A. (2024). How to Interpret Statistical Models Using marginaleffects for R and Python. Journal of Statistical Software, 111(9), 1–32. https://doi.org/10.18637/jss.v111.i09

Heiss, A. (2022, Mai 20). Marginalia: A Guide to Figuring Out What the Heck Marginal Effects, Marginal Slopes, Average Marginal Effects, Marginal Effects at the Mean, and All These Other Marginal Things Are. https://doi.org/10.59350/40xaj-4e562

```{r echo = FALSE} #| echo: false #| warning: false #| message: false pacman::p_load(tidyverse, readxl, knitr, kableExtra, see, latex2exp, patchwork, parsnip, conflicted) cbbPalette <- c("#000000", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7") conflicts_prefer(dplyr::filter) ``` # Modellieren in R {#sec-modeling-R} *Letzte Änderung am `r format(fs::file_info("stat-modeling-R.qmd")$modification_time, '%d. %B %Y um %H:%M:%S')`* > *"Ich weiß nicht weiter; Ich will mich verändern; Doch wie fang ich's an?" --- Tocotronic, Die Unendlichkeit* ::: {.callout-caution appearance="simple"} ## Stand des Kapitels: Stand-By (seit 01.2025) Dieses Kapitel ist auf Stand-By gesetzt, da ich die Thematik des Kapitels aktuell nicht in meiner Lehre oder der statistischen Beratung benötige. Archivierte Kapitel werden nicht von mir weiter gepflegt oder ergänzt. Auftretende Fehler werden aber natürlich beseitigt, wenn die Fehler mir auffallen. ::: Dieses Startkapitel gibt dir nochmal eine Übersicht über das statistische Modellieren in R. Hier liegt vor allem der Fokus auf R. Es gibt eben eine Reihe von *zusätzlichen* Pakten, die es dir erlauben noch mehr aus einem statistischen Modell rauszuholen. Am Ende wurde es mir dann aber zu detailliert alle Pakete in jedem Kapitel vorzustellen und anzuwenden. Das ist auch nicht immer sinnig. Häufig willst du erstmal *nur* das Modell rechnen. Später kannst du dann noch tiefer ins Detail gehen oder aber komplexere Verfahren nutzen. Ich tue mich also etwas schwer dieses Kapitel einzuordnen. Entweder packen wir es ans Ende vom statistischen Modellieren und schauen, dann wie wir alles in R machen. Das steht aber etwas der Intuition entgegen, dass wir in jedem Kapitel zum statistischen Modellieren ja schon was selber machen wollen. In R gibt es dazu dann noch sehr gute Pakete, die das Modellieren sehr viel einfacher machen, dabei dann aber auch für den Anfänger etwas komplexer sind. Ich habe mich daher entschieden, diese aktuelle und komplexere Modellierung einmal hier am Anfang vorzustellen und in den folgenden Kapiteln teilweise darauf zu verweisen, wenn ich es sinnig fand. Du kannst alle Modelle auf althergebrachte Art und Weise rechnen ohne was zu verpassen. Aber manchmal möchte man dann auch effizienter Modellieren. Dafür ist dann dieses Kapitel da: Eine erweiterte Idee von der statistischen Modellierung zu erlangen. Fangen wir also erstmal mit der naheliegenden Frage an. Was modellieren wir eigentlich? : Wir modellieren die Varianz in den Daten oder anders formuliert, wir versuchen die Varianz in den Daten den jeweiligen Quellen zuzuordnen und darüber dann neue Erkenntnisse zu erlangen. **Was ist ein Modell** [Statistical Thinking for the 21st Century --- What is a model?](https://statsthinking21.github.io/statsthinking21-core-site/fitting-models.html#what-is-a-model) [Models are about what changes, and what doesn't](https://magesblog.com/post/modelling-change/) Als nächstes wollen wir uns die Frage stellen, was ist eigentlich das ziel des Modellierens? Wir wollen ja mit der Modellierung der Varianz irgendwas erreichen. In der @fig-scatter-modeling-R-01 siehst du einmal die drei großen Fragefelder, die wir mit einer Modellierung bearbeiten können. ```{r} #| echo: false #| message: false #| warning: false #| label: fig-scatter-modeling-R-01 #| fig-align: center #| fig-height: 5 #| fig-width: 15 #| fig-cap: "Verschiedene Ziele und Möglichkeiten des statistischen Modellierens. Grob können die Möglichkeiten in drei große thematische Zusammenhänge eingeteilt werden. **(A)** *Kausales Modell* -- Wie verändert sich $y$, wenn sich $x_1$ ändert? Wie ist der numerische Zusammenhang zwischen $y$ und $x$? **(B)** *Prädikitives Modell* -- Wenn $y$ und $x$ gemessen wurden, wie sehen dann die Werte von $y$ für neue $x$-Werte aus? Können wir mit $x$ die Werte in $y$ vorhersagen? **(C)** *Clusteranalyse* -- Haben wir in unseren Variablen $x_1$ und $x_2$ eine oder mehrere unbekannte Cluster vorliegen? Wie gehören die Beobachtungen gegeben $x_1$ und $x_2$ zusammen? *[Zum Vergrößern anklicken]*" source("images/R/stat-modeling-R-01.R") p11 + p12 + p13 + plot_layout(ncol = 3) + plot_annotation(tag_levels = 'A') ``` In diesem Kapitel wollen wir uns auch mit dem Thema marginale Effektmodelle (eng. *marginal effect models*) beschäftigen. Auch hier ist der deutsche Begriff nicht gebräuchlich, so dass ich mich hier dann immer auf die *marginal effect models* beziehen werde. Ein zweideutiger Aspekt der *marginal effect models* besteht darin, dass das Wort "marginal" auf zwei verschiedene und etwas entgegengesetzte Arten verwendet wird: - Bei "marginalen Effekten" (eng. *marginal effects*) beziehen wir uns auf die Auswirkung einer winzigen (marginalen) Änderung des Regressors auf das Ergebnis. Dies ist eine Steigung oder eine Ableitung. - Der Begriff "marginaler Mittelwert" (eng. *marginal means*) bezieht sich auf den Prozess der Marginalisierung über die Zeilen eines Vorhersagerasters. Dies ist ein Durchschnitt oder ein Integral. Mehr dazu dann in dem Abschnitt zu den *marginal effect models* und dem [R Paket `{marginaleffects}`](https://marginaleffects.com/). Dafür verweise ich auf @heiss2022 und das Onlinetutorium [Marginalia: A guide to figuring out what the heck marginal effects, marginal slopes, average marginal effects, marginal effects at the mean, and all these other marginal things are](https://www.andrewheiss.com/blog/2022/05/20/marginalia/). Heiss erklärt dort sehr schön die Zusammenhänge. Dazu dann noch der Verweis auf die Webseite [Model to Meaning](https://marginaleffects.com/bonus/get_started.html) von @heiss2024 um nochmal tiefer in das Modellieren von Daten einzusteigen. ## Genutzte R Pakete Wir wollen folgende R Pakete in diesem Kapitel nutzen. ```{r echo = TRUE} #| message: false pacman::p_load(tidyverse, emmeans, multcomp, conflicted) conflicts_prefer(dplyr::select) ``` An der Seite des Kapitels findest du den Link *Quellcode anzeigen*, über den du Zugang zum gesamten R-Code dieses Kapitels erhältst. ## Theoretisches Modellieren in R ### R Pakete Neben den R Paketen, die wir in den jeweiligen Kapiteln brauchen, kommen noch folgende R Pakete immer wieder dran. Deshalb sind die R Pakete hier schon mal mit den jeweiligen Internetseiten aufgeführt. - Das Buch [Tidy Modeling with R](https://www.tmwr.org/) gibt nochmal einen tieferen Einblick in das Modellieren in R. Wir immer, es ist ein Vorschlag aber kein Muss. - Das [R Paket `{parameters}`](https://easystats.github.io/parameters/index.html) nutzen wir um die Parameter eines Modells aus den Fits der Modelle zu extrahieren. Teilweise sind die Standardausgaben der Funktionen sehr unübersichtich. Hier hilft das R Paket. - Das [R Paket `{performance}`](https://easystats.github.io/performance/) hilft uns zu verstehen, ob die Modelle, die wir gefittet haben, auch funktioniert haben. In einen mathematischen Algorithmus können wir alles reinstecken, fast immer kommt eine Zahl wieder raus. - Das [R Paket `{tidymodels}`](https://tidymodels.tidymodels.org/) nutzen wir als *das* R Paket um mit Modellen umgehen zu können und eine *Vorhersage neuer Daten* zu berechnen. Das Paket `{tidymodels}` ist wie das Paket `{tidyverse}` eine Sammlung an anderen R Paketen, die wir brauchen werden. ### Quantity and tests [Kontrafaktisches Konditional](https://de.wikipedia.org/wiki/Kontrafaktisches_Konditional?wprov=sfti1) > "Wenn ... der Fall wäre, dann wäre ..." [Counterfactual comparisons](https://marginaleffects.com/chapters/comparisons.html) ### Generalisierung von `lm()` zu `glm()` und `[g]lmer()` - Die Funktion `lm()` nutzen wir, wenn das Outcome $y$ einer Normalverteilung folgt. - Die Funktion `glm()` nutzen wir, wenn das Outcome $y$ einer *andere* Verteilung folgt. - Die Funktion `lmer()` nutzen wir, wenn das Outcome $y$ einer Normalverteilung folgt *und* wir noch einen Block- oder Clusterfaktor vorliegen haben. - Die Funktion `glmer()` nutzen wir, wenn das Outcome $y$ einer *andere* Verteilung folgt *und* wir noch einen Block- oder Clusterfaktor vorliegen haben. ![Übersicht der Namen der Funktionen in R für das lm(), glm() und glmer().](images/glm-01.png){#fig-glm-01 fig-align="center" width="100%"} In @fig-glm-02 sehen wir wie wir den Namen einer Regression bilden. Zuerst entscheiden wir, ob wir nur ein $x$ haben oder mehrere. Mit einem $x$ sprechen wir von einem simplen Modell, wenn wir mehrere $x$ haben wir ein multiples Modell. Im nächsten Schritt benennen wir die Verteilung für das Outcome $y$. Dann müssen wir noch entscheiden, ob wir ein gemischtes Modell vorliegen haben, dann schreiben wir das hin. Sonst lassen wir den Punkt leer. Anschließend kommt noch lineares Modell hinten ran. ![Wie bilde ich den Namen einer Regression? Erst beschreiben wir das $x$, dann das $y$. Am Ende müssen wir noch sagen, ob wir ein gemischtes Modell vorliegen haben oder nicht.](images/glm-02.png){#fig-glm-02 fig-align="center" width="100%"} [Das R Paket `{tidymodels}`](https://www.tidymodels.org/) [Tidy Modeling with R](https://www.tmwr.org/) Veraltet aber manchmal ganz nützlich das [R Paket `{modelr}`](https://modelr.tidyverse.org/) ### Konfidenzintervall vs. Prädiktionsintervall Ein Konfidenzintervall gibt den Wertebereich für einen gesuchten Parameter der Grundgesamtheit mit einer bestimmten Wahrscheinlichkeit an. Ein Prognoseintervall gibt den Wertebereich für einen individuellen, zukünftig zu beobachtenden Wert mit einer bestimmten Wahrscheinlichkeit an. #### Konfidenzintervall {.unnumbered .unlisted} Text Konfidenzintervall : Ein Konfidenzintervall gibt den Wertebereich für einen gesuchten Parameter der Grundgesamtheit mit einer bestimmten Wahrscheinlichkeit an. #### Prädiktionsintervall {.unnumbered .unlisted} Text Prädiktionsintervall : Ein Prädiktionsintervall (auch Vorhersageintervall oder Prognoseintervall) gibt den Wertebereich für einen individuellen, zukünftig zu beobachtenden Wert mit einer bestimmten Wahrscheinlichkeit an. [Quantile Regression Forests for Prediction Intervals](https://www.bryanshalloway.com/2021/04/21/quantile-regression-forests-for-prediction-intervals/) [The difference between prediction intervals and confidence intervals](https://robjhyndman.com/hyndsight/intervals/) [P-values for prediction intervals](https://robjhyndman.com/hyndsight/forecasting-pvalues.html) machen keinen Sinn [How NASA didn’t discover the hole in the ozone layer](https://robjhyndman.com/hyndsight/ozone-hole-anomaly.html) ## Praktisches Modellieren in R ### Marginal Effect Models [R Paket `{ggeffects}`](https://strengejacke.github.io/ggeffects/articles/ggeffects.html) [R Paket `{modelbased}`](https://easystats.github.io/modelbased/) [R Paket `{marginaleffects}`](https://marginaleffects.com/) ### Visualisierung [R Paket `{ggeffects}`](https://strengejacke.github.io/ggeffects/articles/ggeffects.html) [R Paket `{gtsummary}`](https://www.danieldsjoberg.com/gtsummary/index.html) [Display univariate regression model results in table](https://www.danieldsjoberg.com/gtsummary/reference/tbl_uvregression.html) ### ...mit `{purrr}` `map()` für mehrere Endpunkte wie bei den Spinatdaten in Application gezeigt. `map2()` mit Familie für `glm` zeigen, wenn wir unterschiedliche Outcomes haben. ```{r} #| warning: false #| message: false cutting_tbl <- read_excel("data/multiple_outcomes.xlsx") |> mutate(trt = as_factor(trt), block = as_factor(block)) |> mutate_if(is.numeric, round, 2) ``` ```{r} #| echo: false #| message: false #| warning: false #| label: tbl-model-r-multiple-out #| tbl-cap: "Auszug aus den Daten ." cutting_raw_tbl <- cutting_tbl |> mutate_all(as.character) rbind(head(cutting_raw_tbl, 4), rep("...", times = ncol(cutting_raw_tbl)), tail(cutting_raw_tbl, 4)) |> kable(align = "c", "pipe") ``` ```{r} family_lst <- lst(shoot = binomial(), leaf = quasipoisson(), flower = quasipoisson(), fruit = quasipoisson(), ca = gaussian(), drymatter = gaussian(), freshweight = gaussian(), height = gaussian()) ``` ```{r} cutting_long_tbl <- cutting_tbl |> pivot_longer(cols = shoot:last_col(), names_to = "outcome", values_to = "rsp") |> arrange(outcome, trt, block) ``` ```{r} cutting_lst <- cutting_long_tbl |> split(~outcome) ``` ```{r} #| message: false #| warning: false glm_lst <- cutting_lst %>% map2(family_lst, ~glm(rsp ~ trt + block + trt:block, data = .x, family = .y)) glm_lst %>% map(pluck, "family") glm_lst %>% map(car::Anova) emm_lst <- glm_lst %>% map(~emmeans(.x, specs = ~ trt, type = "response")) emm_lst %>% map(~contrast(.x, method = "pairwise", adjust = "bonferroni")) %>% map(as_tibble) %>% bind_rows(.id = "outcome") emm_lst %>% map(~cld(.x, Letters = letters, adjust = "bonferroni")) %>% map(as_tibble) %>% bind_rows(.id = "outcome") %>% select(outcome, trt, .group) %>% print(n = 28) ``` ## Verschiebebahnhof [R Paket `{parsnip}`](https://parsnip.tidymodels.org/) in diesem Kapitel wollen Fokus auf `{parsnip}` und dann als Kasten in allen anderen möglichen Modellierungen ergänzen? [Tidymodels, interactions and anova - a short tutorial](https://www.thomasvanhoey.com/post/2021-10-12-tidymodels-interactions/) :::: callout-note ## XX Regression in `{parsnip}` Auch hier können wir die XX Regression in dem R Paket `{parsnip}` realisieren. Ein Vorteil von `{parsnip}` ist, dass wir die Funktionen sehr gut mit dem `|>`-Operator nutzen können. Das ist bei den Funktionen `glm()` und `lm()` in der ursprünglichen Implementierung ja leider nur etwas umständlicher möglich. Deshalb hier einmal die *bessere* Implementierung. ::: panel-tabset ## `{parsnip}` ```{r} linreg_reg_fit <- linear_reg() |> set_engine("lm") |> fit(jump_length ~ grp, data = simple_tbl) ``` Dann die ANOVA ```{r} linreg_reg_fit|> extract_fit_engine() |> anova() ``` ## `{base}` ```{r} linreg_reg_fit <- lm(jump_length ~ grp, data = simple_tbl) ``` ```{r} linreg_reg_fit|> anova() ``` ::: Eine detailliertere Einführung mit mehr Beispielen für die Nutzung vom [R Paket `{parsnip}`](https://parsnip.tidymodels.org/articles/Examples.html) findest du im Kapitel [Modellieren in R](#sec-modeling-R). Hier soll es dann bei der kurzen Gegenüberstellung bleiben. :::: ## Referenzen {.unnumbered}

52.1 Genutzte R Pakete

52.2 Theoretisches Modellieren in R

52.2.1 R Pakete

52.2.2 Quantity and tests

52.2.3 Generalisierung von lm() zu glm() und [g]lmer()

52.2.4 Konfidenzintervall vs. Prädiktionsintervall

Konfidenzintervall

Prädiktionsintervall

52.3 Praktisches Modellieren in R

52.3.1 Marginal Effect Models

52.3.2 Visualisierung

52.3.3 …mit {purrr}

52.4 Verschiebebahnhof

Referenzen

52.2.3 Generalisierung von `lm()` zu `glm()` und `[g]lmer()`

52.3.3 …mit `{purrr}`