46 Modellieren in R
Letzte Änderung am 09. April 2024 um 09:55:04
“Ich weiß nicht weiter; Ich will mich verändern; Doch wie fang ich’s an?” — Tocotronic, Die Unendlichkeit
Dieses Kapitel wird ist eine Baustelle und wird es vermutlich auch über das Jahr 2024 bleiben. Aktuell weiß ich nämlich noch nicht, wo ich mit dem Kapitel hin will und wo es stehen sollte. Daher können Teile des Codes und des Textes kurzfristig keinen Sinn ergeben oder nicht funktional sein. Es ist geplant eine fertige Version noch 2024 vorliegen zu haben.
Dieses Startkapitel gibt dir nochmal eine Übersicht über das statistische Modellieren in R. Ich tue mich 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 darauf dann zu verweisen. 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.
- 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.
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.
Drei Abbildungen mit explorative, prädiktiv und Cluster
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
{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.
<- linear_reg() |>
linreg_reg_fit set_engine("lm") |>
fit(jump_length ~ grp, data = simple_tbl)
Dann die ANOVA
|>
linreg_reg_fitextract_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
<- lm(jump_length ~ grp, data = simple_tbl) linreg_reg_fit
|>
linreg_reg_fitanova()
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.
46.1 Marginal Effect Models
46.2 Modellieren in R
Veraltet aber manchmal ganz nützlich das R Paket {modelr}
46.3 Visualisierung
46.4 Genutzte R Pakete
Wir wollen folgende R Pakete in diesem Kapitel nutzen.
::p_load(tidyverse, emmeans, multcomp, conflicted)
pacmanconflicts_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.
46.5 Daten
Wir
46.6 Wilde Beispiele
<- read_excel("data/multiple_outcomes.xlsx") |>
cutting_tbl mutate(trt = as_factor(trt),
block = as_factor(block)) |>
mutate_if(is.numeric, round, 2)
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 |
<- lst(shoot = binomial(),
family_lst leaf = quasipoisson(),
flower = quasipoisson(),
fruit = quasipoisson(),
ca = gaussian(),
drymatter = gaussian(),
freshweight = gaussian(),
height = gaussian())
<- cutting_tbl |>
cutting_long_tbl pivot_longer(cols = shoot:last_col(),
names_to = "outcome",
values_to = "rsp") |>
arrange(outcome, trt, block)
<- cutting_long_tbl |>
cutting_lst split(~outcome)
<- cutting_lst %>%
glm_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
%>%
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
<- glm_lst %>%
emm_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 x 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
# i 38 more rows
# i 2 more variables: estimate <dbl>, t.ratio <dbl>
%>%
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 x 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"
# i 4 more rows
46.7 Weitere Ideen
map()
für mehrere Endpunkte wie bei den Spinatdaten in Application gezeigt.
map2()
mit Familie für glm
zeigen, wenn wir unterschiedliche Outcomes haben.
46.8 Genutzte 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.
46.9 95% Konfidenzintervall
46.10 Prädiktionsintervall
Quantile Regression Forests for Prediction Intervals
The difference between prediction intervals and confidence intervals
P-values for prediction intervals machen keinen Sinn
46.11 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.
In Abbildung 46.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.