R Code [zeigen / verbergen]
::p_load(tidyverse, magrittr, dlookr, broom, modelsummary,
pacman
see, performance, ggpubr, factoextra, FactoMineR, conflicted)
Letzte Änderung am 25. March 2024 um 08:52:02
“Zu tun, worauf man Lust hat, und nicht, was man muss – das ist eines der ältesten und wichtigsten Ziele der Menschheit. Wir arbeiten hart daran, das zu verdrängen.” — Die Not des Müßiggangs, Magazin Brand Eins
Wir brauchen die Sensitivitätsanalyse wenn wir Beobachtungen aus unseren Daten entfernt oder aber hinzugefügt haben. Sensitivitätsanalysen finden eigentlich in dem Kontext von klinischen Studien statt. Der Trend geht aber natürlich auch nicht an den Agrarwissenschaften vorbei und solltest du den Begriff mal hören, weist du wo Sensitivitätsanalyse hingehören. Machen wir also eine Sensitivitätsanalyse. Was haben wir aber davor gemacht? Das heißt du hast entweder eine Variablenselektion wie im Kapitel 46 beschrieben durchgeführt. Oder aber du hast fehlende Werte wie in Kapitel 47 beschrieben imputiert. Es kann auch sein, dass du Ausreißer aus den Daten entfernt oder aber imputiert hast, wie es in Kapitel 45 beschrieben ist. Im Prinzip kannst du auch alles drei gemacht haben, aber meistens beschränkt sich die Veränderung der Daten nur auf eins der drei Möglichkeiten.
Wie immer brauchen wir natürlich auch Fallzahl. Eine Sensitivitätsanalyse kannst du nicht auf zwanzig bis fünfzig Beobachtungen machen. Du brauchst schon eine gute dreistellige Anzahl, damit du hier sauber Modellieren und Darstellen kannst. Wenn du weniger Beobachtungen hast, dann ist ganz natürlich das einzelne Werte einen riesigen Einfluss haben müssen. Im Zweifel frag einfach einmal bei mir nach, dann können wir die Sachlage diskutieren.
Das ist hier natürlich eine Sensitivitätsanalyse für Arme. Wie man es richtig umfangreich macht, findest du in einem sehr gutem und umfangreichen Tutorial zu What Makes a Sensitivity Analysis?
Dieses Kapitel ist relativ übersichtlich. Wir werden die Modelle nach der jeweiligen algorithmischen Veränderung uns nochmal anschauen und dann deskriptive entscheiden, ob wir eine große Veränderung in den Daten sehen. Es gibt zwar auch die Möglichkeit die Modelle untereinander zu vergleichen, aber ist hier die Aussagekraft nicht so stark. Die Idee hinter dem Modellvergleich ist eher die Anzahl an Spalten zu verändern und nicht die Werte in der Datenmatrix. Deshalb machen wir es zwar, genießen die Sache aber mit Vorsicht.
Wir wollen folgende R Pakete in diesem Kapitel nutzen.
::p_load(tidyverse, magrittr, dlookr, broom, modelsummary,
pacman
see, performance, ggpubr, factoextra, FactoMineR, conflicted)
An der Seite des Kapitels findest du den Link Quellcode anzeigen, über den du Zugang zum gesamten R-Code dieses Kapitels erhältst.
In diesem Beispiel betrachten wir wieder die Gummibärchendaten. Auch hier haben wir echte Daten vorliegen, so dass wir Ausreißer entdecken könnten. Da wir hier auch fehlende Werte in den Daten haben, können wir diese fehlenden Werte auch einfach imputieren und uns dann die Effekte anschauen. Das heißt wir haben also einen idealen Datensatz für unsere Sensitivitätsanalysen.
<- read_excel("data/gummibears.xlsx") |>
gummi_tbl select(gender, age, height, semester) |>
mutate(gender = as_factor(gender))
In der Tabelle 45.2 ist der Datensatz gummi_tbl
nochmal für die ersten sieben Zeilen dargestellt. Wir werden später sehen, wie sich die Fallzahl von \(n = 783\) immer wieder ändert, je nachdem wie wir mit den fehlenden Daten und den Variablen umgehen.
gummi_tbl
. Wir betrachten die ersten sieben Zeilen des Datensatzes.
gender | age | height | semester |
---|---|---|---|
m | 35 | 193 | 10 |
w | 21 | 159 | 6 |
w | 21 | 159 | 6 |
w | 36 | 180 | 10 |
m | 22 | 180 | 3 |
NA | NA | NA | NA |
m | 22 | 180 | 3 |
Wir wollen jetzt als erstes das volle Modell schätzen. Das heißt wir packen alle Variablen in das Modell und rechnen dann die lineare Regression. Wir wollen herausfinden in wie weit das Alter, das Geschlecht und das Semester einen Einfluss auf die Körpergröße von Studierenden hat.
\[ height \sim gender + age + semester \]
Wir haben nichts an den Daten geändert und somit dient unser volles Modell als Benchmark für die anderen. Wenn sich einige Werte der Modellgüten im Vergleich zum vollen Modell ändern, dann wissen wir, dass etwas nicht stimmt.
<- lm(height ~ gender + age + semester, data = gummi_tbl) fit_full
Neben dem vollen Modell rechnen wir auch noch das Nullmodel. Das Nullmodell beinhaltet nur den Intercept und sonst keine Einflussvariable. Wir wollen schauen, ob es überhaupt was bringt eine unserer Variablen in das Modell zu nehmen oder ob wir es auch gleich lassen können. Im Prinzip unsere Kontrolle für das Modellieren.
\[ height \sim 1 \]
In R fitten wir das Nullmodell in dem wir keine Variablen mit in das Modell nehmen sondern nur eine 1 schreiben. Wir haben dann nur den Intercept mit in dem Modell und sonst nichts. Was wir schon aus den anderen Kapiteln wissen ist, dass das Nullmodell ein schlechtes Modell sein wird.
<- lm(height ~ 1, data = gummi_tbl) fit_null
Wir schauen uns die Modelle hier nicht weiter an, da es uns nur im Vergleich zu den anderen Modellen interessiert.
Teilweise können wir eine Überprüfung auf Ausreißer nur auf einen Datensatz ohne fehlende Werte durchführen. Hier beißt sich dann die Katze in den Schwanz. Deshalb nutzen wir die Funktion diagnose_outlier()
, die intern die fehlenden Werte entfernt. Das ist natürlich kein richtiges Vorgehen! Aber wir nutzen ja diesen Abschnitt nur als Beispiel. Du findest die Detektion von Ausreißern im Kapitel 45 beschrieben.
diagnose_outlier(gummi_tbl)
# A tibble: 3 × 6
variables outliers_cnt outliers_ratio outliers_mean with_mean without_mean
<chr> <int> <dbl> <dbl> <dbl> <dbl>
1 age 41 5.24 35.7 23.4 22.6
2 height 0 0 NaN 176. 176.
3 semester 67 8.56 8.78 2.55 1.89
Wir sehen, dass wir in der Variable age
und semester
nach der Funktion zu urteilen Ausreißer gefunden haben. Deshalb werden wir jetzt diese Ausreißer durch die Funktion imputate_outlier()
entsprechend ersetzen. Mal schauen, ob wir damit eine substanzielle Änderung in der Modellierung erhalten.
<- gummi_tbl |>
gummi_out_imp_tbl mutate(age = imputate_outlier(gummi_tbl, age, method = "capping"),
semester = imputate_outlier(gummi_tbl, semester, method = "capping"))
Nun modellieren wir noch mit unseren ersetzten und angepassten Daten die Körpergröße und erhalten den Modellfit zurück. Am Ende des Kapitels werden wir dann alle Modelle gegenüberstellen und miteinander vergleichen.
<- lm(height ~ gender + age + semester, data = gummi_out_imp_tbl) fit_outlier
Wir schauen uns das Modell hier nicht weiter an, da es uns nur im Vergleich zu den anderen Modellen interessiert.
Nehmen wir wieder den Gummibärechendatensatz von neuen und imputieren diesmal die fehlenden Werte mit einer univariaten Imputation. Wir machen uns hier nicht die Mühe ein multivariates Verfahren zu nutzen. Das könnte man tun, aber wir wollen hier ja nur den Weg aufzeigen, wie wir den Vergleich der Modelle zur Sensitivitätsanalyse durchführen. Du findest die Imputation von fehlenden Werten im Kapitel 47 beschrieben.
In unserem Fall imputieren wir alle numerischen Variablen mit dem Mittelwert und die kategoriale Variable mit der Methode rpart
. Damit haben wir dann keine fehlenden Werte mehr in den Daten und somit sollte das jetzt auch unserer größter Datensatz für die lineare Regression sein. Nicht vergessen, sobald wir einen fehlenden Wert bei einer Variable in einem Modell haben, fällt die ganze Beobachtung aus dem Modell heraus.
<- gummi_tbl |>
gummi_imp_tbl mutate(age = imputate_na(gummi_tbl, age, method = "mean"),
gender = imputate_na(gummi_tbl, gender, method = "rpart"),
height = imputate_na(gummi_tbl, height, method = "median"),
semester = imputate_na(gummi_tbl, semester, method = "mode"))
Dann rechnen wir noch schnell das Modell für die imputierten Daten. Am Ende des Kapitels werden wir dann alle Modelle gegenüberstellen und miteinander vergleichen.
<- lm(height ~ gender + age + semester, data = gummi_imp_tbl) fit_imp
Wir schauen uns das Modell hier nicht weiter an, da es uns nur im Vergleich zu den anderen Modellen interessiert.
Für die Variablenselektion machen wir es uns sehr einfach. Wir müssen ja nur eine Spalte aus den Daten werfen, mehr ist ja Variablenselektion auch nicht. Wir machen dort nur eine algorithmengetriebene Auswahl. In diesem Fall entscheide ich einfach zufällig welche Variable aus dem Modell muss. Du findest die Variablen Selektion im Kapitel 46 beschrieben. Somit nehmen wir an, wir hätten eine Variablenselektion durchgeführt und die Variable semester
aus dem Modell entfernt.
<- lm(height ~ gender + age, data = gummi_tbl) fit_var_select
Auch dieses Modell schauen wir nicht weiter an, da es uns nur im Vergleich zu den anderen Modellen interessiert.
Kommen wir zu dem eigentlichen Modellvergleich. In Tabelle 48.2 sehen wir den Modellvergleich aller fünf Modelle aus diesem Kapitel. Dazu nutzen wir die Funktion modelsummary()
aus dem R Paket {modelsummary}
. Wir vergleichen die Modelle untereinander aber vor allem mit dem vollen Modell. Das volle Modell basiert ja auf den ursprünglichen nicht veränderten Daten. Den Intercept können wir erstmal ignorieren. Spannend ist, dass sich der Effekt von gender
auf die Körpergröße durch die Imputation um eine Einheit ändert. Der Effekt des Alters verfünffacht sich durch die Outlieranpassung und verdoppelt sich durch die Imputation. Durch die Imputation wird der Effekt des Semesters abgeschwächt.
Wenn wir auf das \(R^2_{adj}\) schauen, dann haben wir eine Verschlechterung durch die Imputation. Sonst bleibt der Wert mehr oder minder konstant. Das ist ein gutes Zeichen, dass wir unser Modell nicht vollkommen an die Wand gefahren haben durch unsere Änderung der Daten. Das \(AIC\) wird folglich für die Imputationsdaten sehr viel schlechter und nähert sich dem Nullmodell an. Das ist wirklcih kein gutes Zeichen für die Imputation. Da haben wir mehr kaputt als heile gemacht. Wir sehen keinen Efdekt bei dem Fehler \(RMSE\), der noch nach dem Fit des Modell übrig bleibt. Aber das kann passieren. Nicht jede Maßzahl muss sich auch ändern. Deshalb haben wir ja mehrere Maßzahlen vorliegen.
modelsummary(lst("Null Modell" = fit_null,
"Volles Modell" = fit_full,
"Outlier" = fit_outlier,
"Imputation" = fit_imp,
"Variablen Selektion" = fit_var_select),
estimate = "{estimate}",
statistic = c("conf.int",
"s.e. = {std.error}",
"t = {statistic}",
"p = {p.value}"))
Null Modell | Volles Modell | Outlier | Imputation | Variablen Selektion | |
---|---|---|---|---|---|
(Intercept) | 176.382 | 184.858 | 187.503 | 183.285 | 184.529 |
[175.596, 177.168] | [182.349, 187.367] | [183.878, 191.128] | [180.821, 185.749] | [182.034, 187.024] | |
s.e. = 0.400 | s.e. = 1.278 | s.e. = 1.846 | s.e. = 1.255 | s.e. = 1.271 | |
t = 440.643 | t = 144.646 | t = 101.557 | t = 145.997 | t = 145.226 | |
p = <0.001 | p = <0.001 | p = <0.001 | p = <0.001 | p = <0.001 | |
genderw | -14.825 | -14.796 | -13.317 | -14.788 | |
[-15.950, -13.700] | [-15.918, -13.674] | [-14.379, -12.255] | [-15.913, -13.663] | ||
s.e. = 0.573 | s.e. = 0.571 | s.e. = 0.541 | s.e. = 0.573 | ||
t = -25.875 | t = -25.892 | t = -24.616 | t = -25.804 | ||
p = <0.001 | p = <0.001 | p = <0.001 | p = <0.001 | ||
age | -0.023 | -0.141 | -0.042 | -0.038 | |
[-0.126, 0.080] | [-0.299, 0.018] | [-0.144, 0.060] | [-0.140, 0.064] | ||
s.e. = 0.052 | s.e. = 0.081 | s.e. = 0.052 | s.e. = 0.052 | ||
t = -0.439 | t = -1.738 | t = -0.802 | t = -0.730 | ||
p = 0.661 | p = 0.083 | p = 0.423 | p = 0.466 | ||
semester | -0.260 | -0.253 | -0.057 | ||
[-0.490, -0.030] | [-0.513, 0.008] | [-0.280, 0.167] | |||
s.e. = 0.117 | s.e. = 0.133 | s.e. = 0.114 | |||
t = -2.223 | t = -1.905 | t = -0.496 | |||
p = 0.027 | p = 0.057 | p = 0.620 | |||
Num.Obs. | 699 | 697 | 697 | 783 | 699 |
R2 | 0.000 | 0.494 | 0.496 | 0.440 | 0.490 |
R2 Adj. | 0.000 | 0.492 | 0.494 | 0.438 | 0.489 |
AIC | 5284.9 | 4801.8 | 4798.7 | 5382.8 | 4817.6 |
BIC | 5294.0 | 4824.6 | 4821.4 | 5406.1 | 4835.8 |
Log.Lik. | -2640.449 | -2395.911 | -2394.337 | -2686.378 | -2404.795 |
F | 204.030 | 334.978 | |||
RMSE | 10.58 | 7.53 | 7.51 | 7.48 | 7.55 |
Das vergleichen von Modellen, die auf unterschiedlichen Daten basieren ist nicht anzuraten. Wir erhalten auch die passende Warnung von der Funktion compare_performance()
aus dem R Paket {performance}
. Dennoch hier einmal der Vergleich. Wir sehen, dass die Modelle mit der Ersetzung der Ausreißer und das volle Modell sich stark ähneln. Das selektierte Modell und das imputierte Modell fallen dagegen ab. Da wir ja hier nicht zeigen wollen, dass sich die Modelle unterscheiden, ist das Ergebnis ähnlich zu der Übersicht. Die Imputation hat so nicht funktioniert.
# Comparison of Model Performance Indices
Name | Model | R2 | R2 (adj.) | RMSE | Sigma | AIC weights | AICc weights | BIC weights | Performance-Score
---------------------------------------------------------------------------------------------------------------------------
fit_outlier | lm | 0.496 | 0.494 | 7.510 | 7.532 | 0.828 | 0.828 | 0.828 | 99.69%
fit_full | lm | 0.494 | 0.492 | 7.527 | 7.549 | 0.172 | 0.172 | 0.171 | 65.42%
fit_var_select | lm | 0.490 | 0.489 | 7.549 | 7.565 | 6.47e-05 | 6.56e-05 | 6.24e-04 | 56.21%
fit_imp | lm | 0.440 | 0.438 | 7.478 | 7.497 | 1.22e-127 | 1.23e-127 | 9.11e-128 | 53.91%
fit_null | lm | 0.000 | 0.000 | 10.575 | 10.583 | 2.17e-106 | 2.24e-106 | 1.98e-103 | 3.42e-102%
Was ist das Fazit aus der Sensitivitätsanalyse für Arme? Nun wir konnten einmal sehen, dass wir auch mit einfachen Werkzeugen Modelle deskriptiv miteinander vergleichen können und dann einen Schluss über die Güte der Detektion von Ausreißern, der Imputation von fehlenden Werten oder aber der Variablenselektion treffen können. Denk immer dran, die Sensitivitätsanalyse findet nach einer sauberen Detektion, Imputation oder Selektion statt und soll nochmal sicherstellen, dass wir nicht künstliche Effekte der Algorithmen modellieren sondern die Effekte in den Daten sehen.