Letzte Änderung am 25. March 2024 um 17:41:41
“Denn im Formulieren des Problems ist die Lösung schon enthalten!” — Käptn Peng & Die Tentakel von Delphi, Meister und Idiot
Die Selektion von Variablen in einem Modell. Ein schwieriges Thema. Entweder kenne ich mein Experiment und habe das Experiment so geplant, dass nur die bedeutenden Variablen mit in dem Experiment sind oder ich habe keine Ahnung. Gut, dass ist überspitzt und gemein formuliert. Wir wollen uns in diesem Kapitel den Fall anschauen, dass du sehr viele Variablen \(x\) erhoben hast und nun statistisch bestimmen willst, welche Variablen nun mit in das finale Modell sollen. Achtung, ich spreche hier nicht von einem Blockdesign oder aber einem Feldexperiment. Da hat die Variablenselektion nichts zu suchen. Daher tritt der Fall der Variablenselektion eher in dem Feld Verhaltensbiologie oder aber auch Ökologie auf. Ebenfalls kann die Anwendung in automatisch erfassten Daten einen Sinn machen. Wir nutzen dann die Variablenselektion (eng. feature selection) zu Dimensionsreduktion des Datensatzes. Der Datensatz soll damit einfacher sein… ob der Datensatz das damit auch wird, ist wieder eine andere Frage.
In diesem Kapitel prügeln wir aber einen statistischen Engel. Wir werden hier mal schauen müssen, was alles geht und was nicht. Variablen Selektion ist faktisch nicht ein Kapitel sondern ein Regal(kilo)meter voll mit Büchern.
Zu der Frage welches Verfahren denn nun das richtige Verfahren zur Selektion von Variablen ist, gibt es die Standardantwort in der Statistik: Es kommt auf die Fragestellung an…. Oder aber was ist wie gut implementiert, dass wir das Verfahren einigermaßen gut nutzen können. Wir gehen daher von einfach zu kompliziert und du musst dann schauen, was du nutzen willst und kannst. Wir müssen zum Beispiel unterscheiden, welcher Verteilung das Outcome \(y\) folgt. Wenn wir ein normalverteiltes \(y\) haben, dann haben wir andere Möglichkeiten, als wenn wir uns ein poissonverteiltes oder binominalverteiltes \(y\) anschauen.
Das R Paket {olsrr} erlaubt eine weitreichende Variablen Selektion, wenn ein normalverteiltes Outcome \(y\) vorliegt.
Im Folgenden will ich kurz die fünf Mythen der Variablenselektion von Heinze und Dunkler (2017) zusammenfassen. Wie immer ersetzt meine deutsche Zusammenfassung und Auswahl nicht das eigenständige Lesen der englischen Orgnialquelle, wenn du die Informationen in deiner Abschlussarbeit zitieren willst.
Im Weiteren sei auch noch auf Heinze u. a. (2018) und Talbot und Massamba (2019) verwiesen. Beide Veröffentlichungen liefern nochmal einen fundierten Block auf die Variablenselektion. Wiederum ist das natürlich nur ein winziger Ausschnitt aus der Literatur zur Variablenselektion. Im Zweifel einfach einmal bei Google Scholar nach Variablenselektion suchen und schauen was so in dem Feld der eigenen Forschung alles gemacht wird.
“Variable selection should always be accompanied by sensitivity analyses to avoid wrong conclusions.” (Heinze und Dunkler 2017, p. 9)
Nachdem wir Variablen aus unseren Daten entfernt haben, ist es üblich noch eine Sensitivitätsanalysen durchzuführen. Wir Vergleich dann das selektierte Modell mit anderen Modellen. Oder wir wollen die Frage beantworten, was hat eigentlich meine Variablenselektion am Ergebnis geändert? Habe ich eine wichtige Variable rausgeschmissen? Das machen wir dann gesammelt in dem Kapitel 48 zu den Sensitivitätsanalysen.
Wir wollen folgende R Pakete in diesem Kapitel nutzen.
pacman::p_load(tidyverse, magrittr, dlookr,
MASS, ranger, Boruta, broom,
scales, olsrr, gtsummary, parameters,
conflicted)
conflicts_prefer(dplyr::select)
conflicts_prefer(dplyr::filter)An der Seite des Kapitels findest du den Link Quellcode anzeigen, über den du Zugang zum gesamten R-Code dieses Kapitels erhältst.
Um die Variablenselektion einmal durchzuführen nurtzen wir zwei Datensätze. Zum einen den Datensatz zu den Kichererbsen in Brandenburg mit einem normalverteilten Outcome \(y\) mit dryweight. Wir laden wieder den Datensatz in R und schauen uns einmal die Daten in Tabelle 46.1 als Auszug aus dem Tabellenblatt an.
chickpea_tbl <- read_excel("data/chickpeas.xlsx") Wir sehen, dass wir sehr viele Variablen vorleigen haben. Sind denn jetzt alle Variablen notwendig? Oder können auch ein paar Variablen raus aus dem Modell. So viele Beobachtungen haben wir mit \(n = 95\) ja nicht vorliegen. Daher wollen wir an diesem Datensatz die Variablenselektion unter der Annahme eines normalverteilten \(y\) durchgehen.
| temp | rained | location | no3 | fe | sand | forest | dryweight |
|---|---|---|---|---|---|---|---|
| 25.26 | high | north | 5.56 | 4.43 | 63 | >1000m | 253.42 |
| 21.4 | high | northeast | 9.15 | 2.58 | 51.17 | <1000m | 213.88 |
| 27.84 | high | northeast | 5.57 | 2.19 | 55.57 | >1000m | 230.71 |
| … | … | … | … | … | … | … | … |
| 29.04 | low | north | 5.64 | 2.87 | 53.27 | >1000m | 236.07 |
| 24.11 | high | northeast | 4.31 | 3.66 | 63 | <1000m | 259.82 |
| 28.88 | low | northeast | 7.92 | 2 | 65.75 | >1000m | 274.75 |
Was wir auch noch wissen, ist wie die Effekte in den Daten wirklich sind. Die Daten wurden ja künstlich erstellt, deshalb hier die Ordnung der Effektstärke für jede Variable. Im Prinzip müsste diese Reihenfolge auch bei der Variablenselektion rauskommen. Schauen wir mal, was wir erhalten.
\[ y \sim \beta_0 + 3 * sand + 2 * temp + 1.5 * rained - 1.2 * forest + 1.1 * no3 \]
Viele Beispiele laufen immer unter der Annahme der Normalverteilung. Deshalb als zweites Beispiel nochmal die Daten von den infizierten Ferkeln mit einem binomialverteilten Outcome \(y\) mit infected. Auch hier können wir uns den Auszug der Daten in Tabelle 46.2 anschauen.
pig_tbl <- read_excel("data/infected_pigs.xlsx") Das schöne an diesem Datensatz ist jetzt, dass wir mit \(n = 412\) Beobachtungen sehr viele Daten vorliegen haben. Daher können wir auch alle Methoden gut verwenden und haben nicht das Problem einer zu geringen Fallzahl.
| age | sex | location | activity | crp | frailty | bloodpressure | weight | creatinin | infected |
|---|---|---|---|---|---|---|---|---|---|
| 61 | male | northeast | 15.31 | 22.38 | robust | 62.24 | 19.05 | 4.44 | 1 |
| 53 | male | northwest | 13.01 | 18.64 | robust | 54.21 | 17.68 | 3.87 | 1 |
| 66 | female | northeast | 11.31 | 18.76 | robust | 57.94 | 16.76 | 3.01 | 0 |
| … | … | … | … | … | … | … | … | … | … |
| 61 | male | northwest | 15.26 | 23.1 | robust | 57.18 | 15.55 | 3.08 | 1 |
| 59 | female | north | 13.13 | 20.23 | robust | 56.64 | 18.6 | 3.41 | 0 |
| 63 | female | north | 10.01 | 19.89 | robust | 57.46 | 18.6 | 4.2 | 1 |
Auch in diesem Beispiel wurden die Daten von mir mit folgenden Effekten generiert. Schauen wir mal, was die Variablenselektion bei der hohen Fallzahl mit den Variablen macht bzw. welche Sortierung am Ende rauskommt.
\[ y \sim \beta_0 + 2 * crp + 0.5 * sex + 0.5 * frailty + 0.2 * bloodpressure + 0.05 * creatinin + 0.01 * weight \]
Damit haben wir unsere beiden Beispiel und jetzt gehen wir mal eine Auswahl an Methoden zur Variablenselektion durch. Besonders hier, haltet den statistsichen Engel nah bei euch. Es wird leider etwas ruppig für den statistischen Engel.
In den folgenden Abschnitten wollen wir uns eine Reihe an Methoden anschauen um eine Variablenselektion durchzuführen. Dabei gehen wir von einfach nach komplex. Wobei das komplex eher die Methode und nicht die Anwendung meint. Wir nutzen R Pakete und gehen nicht sehr ins Detail wie der Algorithmus jetzt die Auswahl trifft. Für den Hintergrund sind dann die Verweise auf die anderen Kapitel.
Manchmal ist der Anfang auch das Ende. Wir müssen ja gar keinen Algorithmus auf unsere Daten loslassen um eine Variablenselektion durchzuführen. Deshalb möchte ich gleich den ersten Abschnitt mit einem Zitat von Heinze und Dunkler (2017) beginnen.
“Oft gibt es keinen wissenschaftlichen Grund, eine (algorithmische) Variablenauswahl durchzuführen. Insbesondere erfordern Methoden der (algorithmische) Variablenselektion einen viel größeren Stichprobenumfang als die Schätzung eines multiplen Modells mit einem festen Satz von Prädiktoren auf der Grundlage (klinischer) Erfahrung.” (Übersetzt und ergänzt nach Heinze und Dunkler 2017, p. 9)
Fazit dieses kurzen Abschnitts. Wir können auf alles Folgende einfach verzichten und uns überlegen welche Variablen sinnvollerweise mit ins Modell sollen und das mit unserem Expertenwissen begründen. Gut, und was ist, wenn ich kein Experte bin? Oder wir aber wirklich Neuland betreten? Dann können wir eine Reihe anderer Verfahren nutzen um uns algortimisch einer Wahrheit anzunähern.
Und weiter geht es mit Zitaten aus Heinze und Dunkler (2017) zu der Variablenselektion. Dazu musst du wissen, dass die univariate Vorselektion sehr beliebt war und auch noch ist. Denn die univariate Vorselektion ist einfach durchzuführen und eben auch gut darzustellen.
“Obwohl die univariable Vorfilterung nachvollziehbar und mit Standardsoftware leicht durchführbar ist, sollte man sie besser ganz vergessen, da sie für die Erstellung multivariabler Modelle weder Voraussetzung noch von Nutzen ist.” (Übersetzt nach Heinze und Dunkler 2017, p. 8)
Wir sehen also, eigentlich ist die univariate Variablensleketion nicht so das gelbe vom Ei, aber vielleicht musst die Variablenselektion durchführen, so dass her die Lösung in R einmal dargestellt ist. Wir nutzen einmal die gaussian lineare Regression für den Kichererbsendatensatz. Es ist eine ganze Reihe an Code, das hat aber eher damit zu tun, dass wir die Modellausgabe noch filtern und anpassen wollen. Die eigentliche Idee ist simple. Wir nehmen unseren Datensatz und pipen den Datensatz in select und entfernen unser Outcome drymatter. Nun iterieren wir für jede Variable .x im Datensatz mit der Funktion map() und rechnen in jeder Iteration eine gaussian lineare Regression. Dann entferne wir noch den Intercept und sortieren nach den \(p\)-Werten.
chickpea_tbl |>
select(-dryweight) |>
map(~glm(dryweight ~ .x, data = chickpea_tbl, family = gaussian)) |>
map(tidy) |>
map(filter, term != "(Intercept)") |>
map(select, -term, -std.error, -statistic) |>
bind_rows(.id="term") |>
arrange(p.value) |>
mutate(p.value = pvalue(p.value),
estimate = round(estimate, 2))# A tibble: 8 × 3
term estimate p.value
<chr> <dbl> <chr>
1 sand 2.9 <0.001
2 temp 1.65 0.011
3 no3 2.51 0.037
4 fe 2.61 0.091
5 location -7.27 0.130
6 rained 4.95 0.165
7 forest -4.39 0.197
8 location -2.87 0.434 Würden wir nur nach dem Signifikanzniveau von 5% gehen, dann hätten wir die Variablen sand und location selektiert. Bei der selektion mit dem \(p\)-Wert wird aber eher eine Schwelle von 15.7% vorgeschlagen (Heinze und Dunkler 2017, p. 9). Daher würden wir auch noch no3 und temp mit Selektieren und in unser Modell nehmen.
Es gibt ja immer zwei Wege nach Rom. Deshalb hier auch nochmal die Funktion tbl_uvregression() aus dem R Paket {gtsummary}, die es erlaubt die univariaten Regressionen über alle Variablen laufen zu lassen. Wir kriegen dann auch eine schöne Tabelle 46.3 wieder.
Das R Paket {gtsummary} erlaubt es Ergebnisse der Regression in dem Tutorial: tbl_regression gut darzustellen.
chickpea_tbl |>
tbl_uvregression(
method = glm,
y = dryweight,
method.args = list(family = gaussian),
pvalue_fun = ~style_pvalue(.x, digits = 2)
) |>
add_global_p() |> # add global p-value
add_q() |> # adjusts global p-values for multiple testing
bold_p() |> # bold p-values under a given threshold (default 0.05)
bold_p(t = 0.10, q = TRUE) |> # now bold q-values under the threshold of 0.10
bold_labels()tbl_uvregression().
| Characteristic | N | Beta | 95% CI1 | p-value | q-value2 |
|---|---|---|---|---|---|
| temp | 95 | 1.6 | 0.41, 2.9 | 0.009 | 0.032 |
| rained | 95 | 0.16 | 0.23 | ||
| high | — | — | |||
| low | 4.9 | -2.0, 12 | |||
| location | 95 | 0.31 | 0.31 | ||
| north | — | — | |||
| northeast | -2.9 | -10, 4.3 | |||
| west | -7.3 | -17, 2.1 | |||
| no3 | 95 | 2.5 | 0.18, 4.8 | 0.035 | 0.081 |
| fe | 95 | 2.6 | -0.38, 5.6 | 0.087 | 0.15 |
| sand | 95 | 2.9 | 2.5, 3.3 | <0.001 | <0.001 |
| forest | 95 | 0.19 | 0.23 | ||
| <1000m | — | — | |||
| >1000m | -4.4 | -11, 2.2 | |||
| 1 CI = Confidence Interval | |||||
| 2 False discovery rate correction for multiple testing | |||||
Nun führen wir die univariate Regression erneut auf den Ferkeldaten aus. Hier ändern wir nur die family = binomial, da wir hier jetzt eine logistische lineare Regression rechnen müssen. Unser Outcome infected ist ja \(0/1\) codiert. Sonst ändert sich der Code nicht.
pig_tbl |>
select(-infected) |>
map(~glm(infected ~ .x, data = pig_tbl, family = binomial)) |>
map(tidy) |>
map(filter, term != "(Intercept)") |>
map(select, -term, -std.error, -statistic) |>
bind_rows(.id="term") |>
arrange(p.value) |>
mutate(p.value = pvalue(p.value),
estimate = round(estimate, 2))# A tibble: 12 × 3
term estimate p.value
<chr> <dbl> <chr>
1 crp 0.96 <0.001
2 bloodpressure 0.08 0.013
3 creatinin 0.11 0.133
4 location -0.4 0.141
5 sex -0.29 0.188
6 location -0.24 0.435
7 frailty 0.16 0.645
8 activity 0.02 0.698
9 frailty -0.12 0.699
10 location 0.12 0.712
11 weight -0.02 0.765
12 age 0 0.877 In diesem Fall reicht die Schwelle von 15.7% nur für zwei Variablen (Heinze und Dunkler 2017, p. 9). Wir erhalten die Variablen crp und bloodpressure für das Modell selektiert.
In der Tabelle 46.4 sehen wir dann nochmal die Anwendung der Funktion tbl_uvregression() auf den Ferkeldatensatz. Ich musste hier die Option pvalue_fun = ~style_pvalue(.x, digits = 2) entfernen, da sonst die Variable crp keinen \(p\)-Wert erhält. Leider sehe ich den \(p\)-Wert mit \(<0.001\) in meiner Ausgabe in R aber wie du siehst, wird die Tabelle auf der Webseite nicht korrekt angezeigt. Das Problem von automatischen Tabellen. Ein Fluch und Segen zugleich. Du musst immer wieder überprüfen, ob die Optionen dann auch für sich und deine Analyse passen.
pig_tbl |>
tbl_uvregression(
method = glm,
y = infected,
method.args = list(family = binomial),
exponentiate = TRUE
) |>
add_global_p() |> # add global p-value
add_nevent() |> # add number of events of the outcome
add_q() |> # adjusts global p-values for multiple testing
bold_p() |> # bold p-values under a given threshold (default 0.05)
bold_p(t = 0.10, q = TRUE) |> # now bold q-values under the threshold of 0.10
bold_labels()tbl_uvregression().
| Characteristic | N | Event N | OR1 | 95% CI1 | p-value | q-value2 |
|---|---|---|---|---|---|---|
| age | 412 | 276 | 1.00 | 0.95, 1.04 | 0.9 | 0.9 |
| sex | 412 | 276 | 0.2 | 0.4 | ||
| female | — | — | ||||
| male | 0.75 | 0.49, 1.15 | ||||
| location | 412 | 276 | 0.3 | 0.5 | ||
| north | — | — | ||||
| northeast | 1.12 | 0.61, 2.10 | ||||
| northwest | 0.67 | 0.39, 1.14 | ||||
| west | 0.79 | 0.43, 1.44 | ||||
| activity | 412 | 276 | 1.02 | 0.90, 1.16 | 0.7 | 0.9 |
| crp | 412 | 276 | 2.62 | 2.14, 3.27 | <0.001 | <0.001 |
| frailty | 412 | 276 | 0.5 | 0.7 | ||
| frail | — | — | ||||
| pre-frail | 1.17 | 0.59, 2.27 | ||||
| robust | 0.88 | 0.46, 1.64 | ||||
| bloodpressure | 412 | 276 | 1.08 | 1.02, 1.15 | 0.012 | 0.053 |
| weight | 412 | 276 | 0.98 | 0.86, 1.12 | 0.8 | 0.9 |
| creatinin | 412 | 276 | 1.12 | 0.97, 1.30 | 0.13 | 0.4 |
| 1 OR = Odds Ratio, CI = Confidence Interval | ||||||
| 2 False discovery rate correction for multiple testing | ||||||
Neben der Berechnung von univariaten logistischen Regressionen ist auch die Darstellung der Daten in einer Tabelle 46.5 bei Medizinern sehr beliebt. Deshalb an dieser Stelle auch die Tabelle 1 (eng. table 1) für die Zusammenfasung der Daten getrennt nach dem Infektionsstatus zusammen mit dem \(p\)-Wert. Ich nutze hier die Funktion tbl_summary() aus dem R Paket {gtsummary}.
pig_tbl |> tbl_summary(by = infected) |> add_p()| Characteristic | 0, N = 1361 | 1, N = 2761 | p-value2 |
|---|---|---|---|
| age | 59.5 (57.0, 63.0) | 60.0 (57.0, 63.0) | 0.9 |
| sex | 0.2 | ||
| female | 47 (35%) | 114 (41%) | |
| male | 89 (65%) | 162 (59%) | |
| location | 0.3 | ||
| north | 36 (26%) | 85 (31%) | |
| northeast | 23 (17%) | 61 (22%) | |
| northwest | 48 (35%) | 76 (28%) | |
| west | 29 (21%) | 54 (20%) | |
| activity | 13.40 (12.25, 14.34) | 13.24 (12.28, 14.54) | 0.8 |
| crp | 19.12 (18.13, 19.83) | 20.57 (19.77, 21.46) | <0.001 |
| frailty | 0.5 | ||
| frail | 18 (13%) | 37 (13%) | |
| pre-frail | 42 (31%) | 101 (37%) | |
| robust | 76 (56%) | 138 (50%) | |
| bloodpressure | 56.2 (54.3, 58.5) | 57.2 (55.1, 59.6) | 0.021 |
| weight | 18.61 (17.34, 19.41) | 18.32 (17.19, 19.60) | 0.8 |
| creatinin | 4.85 (3.67, 5.93) | 4.86 (4.06, 5.85) | 0.3 |
| 1 Median (IQR); n (%) | |||
| 2 Wilcoxon rank sum test; Pearson’s Chi-squared test | |||
Tja, auch hier ist dann die Frage, wie sortiere ich Variablen. Da es sich bei dem table 1-Stil um eine Übersichtstabelle handelt, ist die Tabelle nach den Variablen sortiert. Auch hier finden wir dann die Variablen crp und bloodpressure wieder. Das Problem hierbei ist natürlich, dass sich die \(p\)-Werte unterscheiden. Das muss ja auch so sein, denn eine logitische Regression ist nun mal kein Wilcoxon rank sum test oder ein Pearson’s Chi-squared test.
Fassen wir als Fazit dieses Abschnitts zusammen wie unsere Modelle nach der Variablenslektion aussehen würde. In unserem Beispiel für die Kichererbsen im sandigen Brandenburg würden wir dann folgendes Modell nehmen.
\[ y \sim \beta_0 + sand + location + no3 + temp \]
Unsere infizierten Ferkel würden dann folgendes selektiertes Modell erhalten.
\[ y \sim \beta_0 + crp + bloodpressure \]
Schauen wir mal, was die anderen Algorithmen noch so finden.
Für die gaussian lineare Regression gibt es mit dem R Paket {oslrr} eine große Auswahl an Variable Selection Methods. Einfach mal die Möglichkeiten anschauen, die dort angeboten werden. Wir nutzen jetzt nicht alles was din oslrr möglich ist, sondern nur eien Auswahl. Zuerst müssen wir wieder unser Modell fitten. Wir nehmen alle Variablen mit rein und nutzen die Funktion lm() für ein lineares Modell mit einem normalverteilten Outcome \(y\) mit dryweight.
chickenpea_fit <- lm(dryweight ~ temp + rained + location + no3 + fe + sand + forest,
data = chickpea_tbl)Nun gibt es wirklich viele Lösungen in dem R Paket {oslrr}. Ich möchte einmal die Variante mit ols_step_all_possible präsentieren. In dem Fall rechnen wir alle Modelle die gehen. Und damit meine ich wirklich alle Modelle. Deshalb filtern wir noch nach dem \(R^2_{adj}\) um nicht von dem Angebot erschlagen zu werden. Darüber hinaus möchte ich nur Modelle sehen, die maximal vier Variablen mit in dem Modell haben. Das ist zufällig von mir gewählt… ich will ja ein kleineres Modell haben.
ols_step_all_possible(chickenpea_fit)$result |>
as_tibble() |>
arrange(desc(adjr)) |>
filter(n <= 4) |>
select(predictors, adjr, aic) # A tibble: 98 × 3
predictors adjr aic
<chr> <dbl> <dbl>
1 temp no3 sand forest 0.893 587.
2 temp rained no3 sand 0.886 592.
3 temp no3 fe sand 0.886 593.
4 temp no3 sand 0.884 593.
5 temp location no3 sand 0.882 596.
6 temp rained sand forest 0.881 597.
7 temp sand forest 0.877 599.
8 temp fe sand forest 0.876 600.
9 temp location sand forest 0.875 602.
10 temp rained fe sand 0.872 603.
# ℹ 88 more rowsEine Alternative ist die Funktion ols_step_backward_aic(), die es erlaubt die Selektion anhand dem \(AIC\)-Wert zu machen. Der \(AIC\)-Wert beschreibt die Güte eines Modells und je kleiner der \(AIC\)-Wert ist, desto besser ist das Modell im Vergleich zu anderen Modellen gleicher Art. Da der \(AIC\)-Wert von den Daten abhängt in denen der \(AIC\)-Wert geschätzt wurde, können verschiedene \(AIC\)-Werte nicht übergreifend vergleichen werden.
chick_step_aic <- ols_step_backward_aic(chickenpea_fit)In Abbildung 46.1 sehen wir einmal den Verlauf der \(AIC\)-Wert durch die Entfernung der jeweiligen Variable. Wenn du auf die y-Achse schaust, ist der Effekt numerisch nicht sehr groß. Davon darf man sich aber nicht beeindrucken lassen., Wir erhalten ein besseres Modell, wenn wir Variablen entfernen. Darüber hinaus sehen wir auch eine Sättigung.
plot(chick_step_aic) 
ols_step_backward_aic mit der Reduktion des AIC-Wertes.
Gut, und wie sieht nun unser finales Modell aus? Dafür müssen wir usn aus dem Objekt chick_step_aic das model raus ziehen. Dafür nutzen wir die Funktion pluck(). Dann noch die Ausgabe in die Funktion model_parameters() gepipt und schon haben wir das finale Modell nach \(AIC\)-Werten in einer gaussian linearen Regression.
pluck(chick_step_aic, "model") |>
model_parameters()Parameter | Coefficient | SE | 95% CI | t(89) | p
-----------------------------------------------------------------------
(Intercept) | -2.09 | 9.44 | [-20.85, 16.67] | -0.22 | 0.825
temp | 2.37 | 0.22 | [ 1.94, 2.81] | 10.88 | < .001
rained [low] | 1.79 | 1.18 | [ -0.56, 4.15] | 1.51 | 0.133
no3 | 1.42 | 0.40 | [ 0.62, 2.22] | 3.53 | < .001
sand | 3.03 | 0.12 | [ 2.80, 3.26] | 26.12 | < .001
forest [>1000m] | -3.17 | 1.11 | [ -5.38, -0.95] | -2.84 | 0.006 Spannenderweise ist location nicht mehr im finalen Modell plus die Variable location flog auch sehr früh raus. Das passt auch besser zu den Daten. Ich hatte die Daten so gebaut, dass der Ort eigentlich keinen Effekt haben sollte. Wir sehen, dass je nach Verfahren was anderes herauskommt. Aber Achtung, das schrittweise Verfahren ist der Auswahl nach \(p\)-Werten auf jeden Fall vorzuziehen!
stepAICWas das R Paket {oslrr} für die gaussian linear Regression kann, kann das R Paket {MASS} mit der Funktion stepAIC für den Rest der möglichen Verteilungen. Da wir mit dem Fekerldatensatz ein binominales Outcome \(y\) mit infected vorliegen haben nutzen wir un die Funktion stepAIC(). Wir hätten auch den Kichererbsendatensatz mit der Funktion bearbeiten können, aber im Falle der Normalverteilung stehen uns dann eben noch andere Algorithmen zu Verfügung. Wie immer müssen wir zuerst das volle Modell mit der Funktion glm() fitten. Wir geben noch die Verteilungsfamilie mit family = binomial noch mit an und definieren so eine logistische lineare Regression.
fit <- glm(infected ~ age + sex + location + activity + crp + frailty + bloodpressure + weight + creatinin,
data = pig_tbl, family = binomial)Nachdem wir das Modell gefittet haben, können wir das Modell direkt in die Funktion stepAIC stecken. Wir nutzen noch die Option direction = "backward" um eine Rückwärtsselektion durchzuführen.
fit_step <- stepAIC(fit, direction = "backward")Start: AIC=421.22
infected ~ age + sex + location + activity + crp + frailty +
bloodpressure + weight + creatinin
Df Deviance AIC
- location 3 398.49 418.49
- frailty 2 396.61 418.61
- weight 1 395.22 419.22
- age 1 395.24 419.24
- activity 1 395.47 419.47
- sex 1 395.93 419.93
- creatinin 1 396.74 420.74
<none> 395.22 421.22
- bloodpressure 1 399.67 423.67
- crp 1 504.90 528.90
Step: AIC=418.49
infected ~ age + sex + activity + crp + frailty + bloodpressure +
weight + creatinin
Df Deviance AIC
- frailty 2 399.68 415.68
- weight 1 398.50 416.50
- age 1 398.50 416.50
- activity 1 398.83 416.83
- sex 1 399.01 417.01
- creatinin 1 400.10 418.10
<none> 398.49 418.49
- bloodpressure 1 403.24 421.24
- crp 1 509.50 527.50
Step: AIC=415.68
infected ~ age + sex + activity + crp + bloodpressure + weight +
creatinin
Df Deviance AIC
- weight 1 399.69 413.69
- age 1 399.73 413.73
- activity 1 400.06 414.06
- sex 1 400.30 414.30
- creatinin 1 401.21 415.21
<none> 399.68 415.68
- bloodpressure 1 404.17 418.17
- crp 1 511.29 525.29
Step: AIC=413.69
infected ~ age + sex + activity + crp + bloodpressure + creatinin
Df Deviance AIC
- age 1 399.74 411.74
- activity 1 400.09 412.09
- sex 1 400.40 412.40
- creatinin 1 401.22 413.22
<none> 399.69 413.69
- bloodpressure 1 404.18 416.18
- crp 1 512.26 524.26
Step: AIC=411.74
infected ~ sex + activity + crp + bloodpressure + creatinin
Df Deviance AIC
- activity 1 400.11 410.11
- sex 1 400.45 410.45
- creatinin 1 401.40 411.40
<none> 399.74 411.74
- bloodpressure 1 404.20 414.20
- crp 1 512.28 522.28
Step: AIC=410.11
infected ~ sex + crp + bloodpressure + creatinin
Df Deviance AIC
- sex 1 400.47 408.47
- creatinin 1 401.86 409.86
<none> 400.11 410.11
- bloodpressure 1 404.72 412.72
- crp 1 513.77 521.77
Step: AIC=408.47
infected ~ crp + bloodpressure + creatinin
Df Deviance AIC
- creatinin 1 402.21 408.21
<none> 400.47 408.47
- bloodpressure 1 406.37 412.37
- crp 1 514.11 520.11
Step: AIC=408.21
infected ~ crp + bloodpressure
Df Deviance AIC
<none> 402.21 408.21
- bloodpressure 1 408.12 412.12
- crp 1 516.27 520.27Jetzt ist die Selektion durchgelaufen und wir sehen in jeden Schritt welche Variable jeweils entfernt wurde und wie sich dann der \(AIC\)-Wert ändert. Wir starten mit einem \(AIC = 425.65\) und enden bei einem \(AIC=415.01\). Schauen wir uns nochmal das finale Modell an.
fit_step |>
model_parameters()Parameter | Log-Odds | SE | 95% CI | z | p
-------------------------------------------------------------------
(Intercept) | -23.54 | 3.14 | [-29.96, -17.61] | -7.49 | < .001
crp | 0.97 | 0.11 | [ 0.76, 1.20] | 8.83 | < .001
bloodpressure | 0.09 | 0.04 | [ 0.02, 0.16] | 2.39 | 0.017 Hier erscheint jetzt noch die Variable sex mit in der Auswahl. Das hat natürlich auch weitreichende Auswirkungen! Es macht schon einen gewaltigen Unterschied, ob wir annehmen das, dass Geschelcht der Ferkel keinen Einfluss auf die Infektion hat oder eben doch. Wir sehen auch hier, dass wir Aufpassen müssen wenn wir eine Variablenselektion durchführen. Aber Achtung, das schrittweise Verfahren ist der Auswahl nach \(p\)-Werten auf jeden Fall vorzuziehen!
rangerIn diesem Abschnitt wollen wir die Variablenselektion mit einem maschinellen Lernverfahren durchführen. Im Bereich des maschinellen Lernens heist die Variablenselektion dann aber Feature Selektion. Wir versuchen jetzt die Selektion auf den Orginaldaten durchzuführen. Eigentlich wird empfohlen die Daten vorher zu normalisieren und dann mit den maschinellen Lernverfahren zu nutzen.
Eine Herausforderung für maschinelle Lernverfahren sind nicht normalisierte Daten. Das heist, dass wir Variablen haben, die kategorial oder kontinuierlich sein können oder aber sehr unterschiedlich von den Einheiten sind. Deshalb wird empfohlen die Daten vorher zu Standardisieren oder zu Normalisieren. In dem Kapitel 19 kannst du mehr über das Transformieren von Daten nachlesen.
Wir nutzen als erstes einen Random Forest Algorithmus wie er in Kapitel 74 beschrieben ist. Es bietet sich hier die Implementation im R Paket {ranger} an. Bevor wir aber einen Random Forest auf unsere Daten laufen lassen, Standardisieren wir unsere Daten nochmal. Überall wo wir einen numerischen Wert als Variableneintrag haben rechnen wir eine \(z\)-Transformation. Wir erhalten dann die standardisierten Daten zurück.
pig_norm_tbl <- pig_tbl |>
mutate(across(where(is.character), as_factor),
across(where(is.numeric), dlookr::transform, "zscore"))Warning: There was 1 warning in `mutate()`.
ℹ In argument: `across(where(is.numeric), dlookr::transform, "zscore")`.
Caused by warning:
! The `...` argument of `across()` is deprecated as of dplyr 1.1.0.
Supply arguments directly to `.fns` through an anonymous function instead.
# Previously
across(a:b, mean, na.rm = TRUE)
# Now
across(a:b, \(x) mean(x, na.rm = TRUE))pig_norm_tbl# A tibble: 412 × 10
age sex location activity crp frailty bloodpressure
<transfrm> <fct> <fct> <transfrm> <transfrm> <fct> <transfrm>
1 0.2156014 male northeast 1.23906607 1.6154076 robust 1.6152124
2 -1.5521155 male northwest -0.15495130 -0.9942285 robust -0.7895759
3 1.3204244 female northeast -1.18531196 -0.9104969 robust 0.3274677
4 -0.2263278 female north 0.03899895 -0.4848611 robust -0.2085934
5 0.6575306 male northwest 0.87540937 1.0502190 robust -0.4391895
6 -1.1101862 male northwest 1.54211333 0.9664873 robust 1.0312351
7 -2.4359739 male west 0.32386336 -0.7569889 pre-frail -0.6488224
8 -1.5521155 male northwest 0.13597407 -0.7569889 robust -0.4122367
9 -0.4472924 female west -0.87014282 1.2735034 robust 0.6838436
10 -0.6682570 male northwest 0.46326510 0.6245832 robust -0.3343731
# ℹ 402 more rows
# ℹ 3 more variables: weight <transfrm>, creatinin <transfrm>,
# infected <transfrm>Den Random Forest rechnen wir mit der Funktion ranger(). Dafür müssen wir wieder unser vollständiges Modell definieren und können dann die Funktion starten. Damit wir eine Variablenwichtigkeit (eng. variable importance) wiederbekommen, müssen wir noch die Option importance = "permutation" verwenden.
fit_raw <- ranger(infected ~ age + sex + location + activity + crp + frailty + bloodpressure + weight + creatinin,
data = pig_tbl, ntree = 1000, importance = "permutation")
pluck(fit_raw, "variable.importance") |>
sort(decreasing = TRUE) |>
round(3) crp creatinin bloodpressure location frailty
0.097 0.007 0.006 0.002 0.001
sex age activity weight
0.000 -0.001 -0.001 -0.002 Wir sehen, dass wir als wichtigste Variable wiederum crp zurückbekommen. Danach wird es schon etwas schwieriger, da die Werte sehr schnell kleiner werden und auch ein Art Plateau bilden. Daher würde man hier nur annehmen, dass crp bedeutend für das Modell ist. Es kann aber auch sein, dass hier eine kontinuierliche Variable sehr vom Algorithmus bevorzugt wurde. Daher schauen wir uns die Sachlage einmal mit den standardisierten Daten an.
fit_norm <- ranger(infected ~ age + sex + location + activity + crp + frailty + bloodpressure + weight + creatinin,
data = pig_norm_tbl, ntree = 1000, importance = "permutation")
pluck(fit_norm , "variable.importance") |>
sort(decreasing = TRUE) |>
round(3) crp creatinin bloodpressure location age
0.415 0.020 0.015 0.005 -0.001
frailty sex weight activity
-0.002 -0.005 -0.011 -0.013 Auch hier erhalten wir ein ähnliches Bild. Audf jeden Fall ist crp bedeutend für den Infektionsstatus. Danach werden die Werte etwas zufällig. Wir können Werte für die variable importance nicht unter Datensätzen vergleichen. Jeder Datensatz hat seine eigene variable importance, die von den Werten in dem Datensatz abhängt.
Wir ziehen als Fazit, dass wir nur crp als bedeutenden Wert für die Klassifikation des Infektionsstatus ansehen würden. Hier stehen wir wirklich etwas wie das Schwein vor dem Uhrwerk, denn was nun richtiger ist, stepAIC oder ranger lässt sich so einfach nicht bewerten. Zum einen wollen wir ja eigentlich mit Random Forest eine Klassifikation durchführen und mit linearen Regressionsmodellen eher kausale Modelle schätzen. Am Ende musst du selber abschätzen, was in das finale Modell soll. Ich kann ja auch den Threshold für den Variablenausschluss selber wählen. Wähle ich einen Threshold von \(0.008\), dann hätte ich crp, weight, bloodpressure und sex mit in dem finalen Modell.
borutaMit dem Boruta Algorithmus steht uns noch eine andere Implemnetierung des Random Forest Algorithmus zu verfügung um Feature Selektion zu betreiben. Wir nutzen wieder den Boruta Algorithmus in seine einfachen Form und gehen nicht tiefer auf alle Optionen ein. Wir nehmen wieder als Beispiel den Datensatz zu den infizierten Ferkeln und nutzen in diesem Fall auch nur die rohen Daten. Über eine Standardisierung könnte man wiederum nachdenken.
Ich empfehle noch das Tutorium Feature Selection in R with the Boruta R Package. Wir gehen hier nicht tiefer auf die Funktionalität von Boruta ein.
Eine Herausforderung für maschinelle Lernverfahren sind nicht normalisierte Daten. Das heist, dass wir Variablen haben, die kategorial oder kontinuierlich sein können oder aber sehr unterschiedlich von den Einheiten sind. Deshalb wird empfohlen die Daten vorher zu Standardisieren oder zu Normalisieren. In dem Kapitel 19 kannst du mehr über das Transformieren von Daten nachlesen.
Um die Funktion Boruta() zu nutzen brauchen wir wieder das Modell und den Datensatz. Sonst ist erstmal nichts weiter anzugeben. Die Funktion läuft dann durch und gibt auch gleich den gewollten Informationshappen. Wichtig ist hierbei, dass wir natürlich noch andere Optionen mit angeben können. Wir können die Anzahl an Iterationen erhöhen und andere Tuning Parameter eingeben. Hier muss man immer schauen was am besten passt. Da würde ich auch immer rumprobieren und auf der Hilfeseite der Funktion ?Boruta einmal nachlesen.
set.seed(20221031)
boruta_output <- Boruta(infected ~ age + sex + location + activity + crp + frailty + bloodpressure + weight + creatinin,
data = pig_tbl)
boruta_outputBoruta performed 99 iterations in 2.946144 secs.
1 attributes confirmed important: crp;
6 attributes confirmed unimportant: activity, age, frailty, location,
sex and 1 more;
2 tentative attributes left: bloodpressure, creatinin;Manchmal ist es super praktisch, wenn eine Funktion einem die Antwort auf die Frage welche Variable bedeutend ist, gleich liefert. Wir erhalten die Information, dass die Variable crp als bedeutsam angesehen wird. Wir können uns den Zusammenhang auch in der Abbildung 46.2 auch einmal anschauen. Die grünen Variablen sind die bedeutenden Variablen.
plot(boruta_output, cex.axis=.7, las=2, xlab="", main="Variable Importance") 
Am Ende finden wir auch hier die Variable crp als einziges als bedeutend wieder. Wenn wir noch Variablen haben die verdächtig oder vorläufig bedeutend sind, angezeigt durch 2 tentative attributes left: bloodpressure, sex, dann können wir noch die Funktion TentativeRoughFix() nutzen. Die Funktion TentativeRoughFix() rechnet die Variablen nochmal nach und versucht alle Variablen in bedeutend oder nicht bedeutend zu klassifizieren. Wir haben ja zwei tentative Variablen in unseren Fall vorliegen, also nutzen wir noch kurz die Funktion um uns auch hier Klarheit zu schaffen.
TentativeRoughFix(boruta_output)Boruta performed 99 iterations in 2.946144 secs.
Tentatives roughfixed over the last 99 iterations.
2 attributes confirmed important: creatinin, crp;
7 attributes confirmed unimportant: activity, age, bloodpressure,
frailty, location and 2 more;Am Ende ist die klare Aussage einer Funktion auch immer ein zweischneidiges Schwert. Wir verlieren jetzt noch die beiden tentative Variablen. Wo wir bei ranger die Qual der Wahl haben, werden wir bei Boruta eher vor vollendete Tatsachen gestellt. Meistens neigt man nach einer Boruta-Analyse nicht dazu noch eine zusätzliche Variable mit ins Modell zu nehmen. Dafür ist dann die Ausgabe zu bestimmt, obwohl die Entscheidung am Ende auch genau so unsicher ist wie von ranger und den anderen Modellen.