Statistik und Datenanalyse: Aufbau

10. Sitzung – Logistische Regression – Machine Learning

Benjamin Fretwurst
PDF-Version der Folien

Orga

Photo courtesy of Lisa Vogel

Lernziele

Theoretisch

Übung 4 und Text Zerback Wirz! lesen! (liegt beides auch in OLAT unter Materialien bzw. Texte)

1 Logistische Regression

1.1 Grund und Folgen

RAQ-ML-«S-Kurve» sieht gerade aus

1.2 Multinominale Regression

Multinomial bedeutet, dass es eine kategoriale (multinomiale) AV gibt. Im Prinzip werden mehrere binäre logistische Regressionen durchgeführt und zusätzlich ausgegeben, wie gut das Gesamtmodel ist. Es kommen die Fit-Kennungen AIC und BIC dazu.

Wir schauen uns nächste Woche Outputs an.

2 Machine Learning – ML

Unsupervised Learning

Beim Unsupervised Learning versucht die «Maschiene» ohne Vorgaben, eigenständig Strukturen in den Daten zu finden. Faktorenanalysen (PCA) wird beim Preprocessing der Daten eingesetzt und Clusteranalysen, um Gruppen zu identifizieren.

Supervised Learning

Beim Supervised Learning werden der Maschiene Vorgaben gemacht, indem ihr Trainingsdaten gegeben werden, die Erkennungsmerkmale («Feature» aka UVs) und Ergebnis enthalten. Die Maschiene sucht nach Modellen, die versucht aus den «Features» das Ergebnis abzuleiten. Die Challenge besteht darin, gute Feature zu finden und wird «Feature Engineering» genannt.

Supervised

Trainieren ➭ Testen ➭ Anwenden

Die verfügbaren Daten werden in einen grösseren und einen kleineren Teil getrennt:

Trainingsdaten, an denen das Modell trainiert wird (Modelling)
Testdaten, an denen das Modell evaluiert wird.

Machine Learning

Übung 4

In der Übung beschäftigen wir uns mit Daten des Untergangs der Titanic. Schauen Sie sich jeweils die Befehle an und finden Sie heraus, was die Befehle machen. Machen Sie sich Notizen, wenn Sie etwas nicht verstehen.

Daten einlesen

Legen Sie für diese Übung einen Ordner an, erstellen Sie dort eine qmd und kopieren Sie die Folgenden Daten in einen dortigen Unterordner «data»:

Suchen und nehmen Sie den Datensatz “train.csv”: Downlaod der Daten: https://www.kaggle.com/competitions/titanic/data

# DATEN_titanic <- read_csv("data/titanic/train.csv") # lese beim ersten Mal die Daten ein

DATEN_titanic <- readRDS("data/titanic/train.RDS") # Lese nach dem ersten Mal so die Daten ein.
saveRDS(DATEN_titanic, "data/titanic/train.RDS") # speichere die Daten

DATEN_titanic |> # mache eine Zusammenfassung der Daten
  summary()

PassengerId Survived Pclass Name
Min. : 1.0 Min. :0.0000 Min. :1.000 Length:891
1st Qu.:223.5 1st Qu.:0.0000 1st Qu.:2.000 Class :character
Median :446.0 Median :0.0000 Median :3.000 Mode :character
Mean :446.0 Mean :0.3838 Mean :2.309
3rd Qu.:668.5 3rd Qu.:1.0000 3rd Qu.:3.000
Max. :891.0 Max. :1.0000 Max. :3.000

 Sex                 Age            SibSp           Parch

Length:891 Min. : 0.42 Min. :0.000 Min. :0.0000
Class :character 1st Qu.:20.12 1st Qu.:0.000 1st Qu.:0.0000
Mode :character Median :28.00 Median :0.000 Median :0.0000
Mean :29.70 Mean :0.523 Mean :0.3816
3rd Qu.:38.00 3rd Qu.:1.000 3rd Qu.:0.0000
Max. :80.00 Max. :8.000 Max. :6.0000
NA’s :177
Ticket Fare Cabin Embarked
Length:891 Min. : 0.00 Length:891 S :644
Class :character 1st Qu.: 7.91 Class :character C :168
Mode :character Median : 14.45 Mode :character Q : 77
Mean : 32.20 NA’s: 2
3rd Qu.: 31.00
Max. :512.33

Kinder Age_z Pclass_f Cabin_D
Mode :logical Min. :-29.279 3:491 Min. :0.000
FALSE:643 1st Qu.: -9.574 2:184 1st Qu.:0.000
TRUE :71 Median : -1.699 1:216 Median :0.000
NA’s :177 Mean : 0.000 Mean :0.229
3rd Qu.: 8.301 3rd Qu.:0.000
Max. : 50.301 Max. :1.000
NA’s :177

Die PassengerId ist einefach eine Identifikationsnummer.

Es gibt dann eine Variable, die “Survived” heisst, die ein Minimum von 0 hat und ein Maximum von 1. Das deutet sehr auf eine Dummy hin. Da der Durchschnitt (“Mean”) = 0.38 ist, wissen wir jetzt schon, dass 38 Prozent der Passagiere überlebt haben (der Mittelwert einer Dummy ist immer der Prozentsatz der 1er-Gruppe). Dann kommt noch der Name als Zeichenvariable, das Alter, das von 0.42 bis 80 geht. Von 177 Personen fehlen die Altersangaben.

Informieren Sie sich über die übrigen Variablen auf (https://www.kaggle.com/competitions/titanic/data)[«Kaggle»].

Daten in Trainings- und Testdaten aufteilen

Was passiert hier?

R-Code

# Setze eine Zufallszahl, damit die Ergebnisse replizierbar sind, also nicht jedes Mal eine neue Zufallszahl gesetzt wird und die Ergebnisse (bisschen) abweichen
set.seed(12345)

# Ziehe eine Zufallsstichprobe aus dem Filmdatensatz und bezeichne ihn als "train", also Trainingsdatensatz, mit dem die "Maschine" trainiert wird
train <- DATEN_titanic |> 
  sample_frac(.75)

# Bilde aus dem Rest der nicht für "train" gezogenen Fälle einen Test-Datensatz, indem nach 'id' die Fälle aus "Filme" das Gegenteil (anti) von zusammengetan (join) werden.
test  <- anti_join(DATEN_titanic, train, by = 'PassengerId')

Sehen kann man nach diesem r-Chunk übrigens nichts, weil nur Datensätze im Hintergrund aufgeteilt wurden. Also suchen wir mal nach guten Datenvisualisierungen.

Übung: Experimentieren Sie mit verschiedenen Zahlen in set.seed() und grösseren und kleineren Werten in sample_frac().

Datenvisualisierung

Eine einfache Darstellungsmöglichkeit ist ein sogenannter “Scatterplot” der die Lage von Fällen in ein Koordinatensystem einteilt, das durch zwei Variablen gebildet wird. Im Beispiel (Abbildung @ref(titanic-Datenvisualisierungen)) ist es das Alter der Passagiere “Age” und der Fahrpreis “Fare”. Als Zweites haben wir ein Balkendiagramm für die Passagierklassen. Im letzten sehr aufwendigen Plotgrafik werden die Zweierbeziehungen aller Variablen dargestellt, also wie sie miteinander korrelieren (obere Nebendiagonalen), wie ihre Verteilung ist (auf der Diagonale mit Namen) und wie ihre gemeinsame Streuung ist, also ein Scatterplot in der unteren Nebendiagonalen. Mehr zu diesen SPLOM finden Sie hier: https://cran.r-project.org/web/packages/psych/vignettes/intro.pdf.

R-Code

# Erstelle einen Scatterplot (Punktewolke)

DATEN_titanic |> 
  ggplot(aes(x = Age, y = Fare)) + 
  geom_point() 

# Erstelle ein Balkendiagramm
DATEN_titanic |> 
  ggplot(aes(x=Pclass)) +
  geom_bar()

# Alle auf einmal
psych::pairs.panels(DATEN_titanic)

Modellbildung für den Fahrpreis

Übung: Experimentieren Sie mit dem Syntax! Kopieren Sie sich die Zeile für das Modell, löschen von “Age_z” bis “Kinder” alles heraus und schätzen Sie mal. Schauen Sie sich das Ergebnis gut an und achten Sie darauf, was passiert, wenn Sie die Summanden für I(Age_z^2) usw. wieder in das Modell tun. Am Ende können Sie versuchen das Modell durch weitere Variablen ergänzen und verbessern oder andere Teile wieder herausnehmen. Manche Variablen wurden erst noch erstellt (zB “Kinder” oder “Age_z”). Die entsprechende Datenaufbereitung.qmd können Sie hier abrufen .

Übung: Interpretieren Sie den Output!

R-Code

fit_titanic <- lm(log(Fare + 1) ~ Sex + Age_z + I(Age_z^2) + Survived + Pclass_f + Kinder, data = train)

broom::tidy(fit_titanic) |>
  mutate(across(where(is.numeric), ~round(.,3)))|>
  gt::gt()

term	estimate	std.error	statistic	p.value
(Intercept)	2.657	0.076	35.047	0.000
Sexmale	-0.316	0.070	-4.520	0.000
Age_z	-0.003	0.003	-1.009	0.314
I(Age_z^2)	0.000	0.000	-0.406	0.685
Survived	-0.077	0.073	-1.050	0.294
Pclass_f2	0.492	0.072	6.859	0.000
Pclass_f1	1.778	0.077	23.104	0.000
KinderTRUE	0.607	0.173	3.515	0.000

R-Code


broom::glance(fit_titanic)|> 
  round(2)|>
  gt::gt()

r.squared	adj.r.squared	sigma	statistic	p.value	df	logLik	AIC	BIC	deviance	df.residual	nobs
0.58	0.57	0.64	103.29	0	7	-517.85	1053.7	1092.26	216.7	528	536

Regressionsoutput mit `sjPlot::tab_model`

R-Code

sjPlot::tab_model(fit_titanic,
                  show.std = TRUE, # zeige die standardisierten Koeffizienten
                  show.est = TRUE, # zeige die unstandardisierten estimates
                  show.r2 = TRUE, # zeige R^2
                  show.fstat = FALSE, # zeige die F-Statistik
                  string.est = "b", # beschrifte die Estimator mit b
                  string.std = "std. b" # beschrifte die standardisierten b mit std. b
                  )

	log(Fare+1)
Predictors	b	std. b	CI	standardized CI	p	std. p
(Intercept)	2.66	0.83	2.51 – 2.81	0.80 – 0.87	<0.001	<0.001
Sex [male]	-0.32	-0.06	-0.45 – -0.18	-0.10 – -0.02	<0.001	0.003
Age z	-0.00	-0.02	-0.01 – 0.00	-0.05 – 0.00	0.314	0.063
Age z^2	-0.00	-0.00	-0.00 – 0.00	-0.02 – 0.01	0.685	0.639
Survived	-0.08	-0.00	-0.22 – 0.07	-0.02 – 0.02	0.294	0.914
Pclass f [2]	0.49	0.06	0.35 – 0.63	0.02 – 0.10	<0.001	0.004
Pclass f [1]	1.78	0.44	1.63 – 1.93	0.39 – 0.48	<0.001	<0.001
KinderTRUE	0.61	0.05	0.27 – 0.95	-0.04 – 0.15	<0.001	0.272
Observations	536
R² / R² adjusted	0.578 / 0.572

Als auch nicht ganz schlechte Alternative kann gtsummary::tbl_regression benutzt werden, was durch weitere gepipte Befehle angepasst und ergänzt werden kann.

Regressionsoutput mit `gtsummary::tbl_regression`

Übung: Was für ein Problem sehen Sie im folgenden Output? Bei welcher Befehlszeile müssten Sie das # wegnehmen, um das Problem zu lösen?

R-Code

## Alternative zum sjPlot-Ouptut. Sie können entscheiden, was Sie besser finden. 
fit_titanic |> 
  gtsummary::tbl_regression() |> 
  gtsummary::add_vif() |> 
#  gtsummary::add_global_p() |> # Gebe die Signifikanzen je Variable raus (Varianzaufklärung)
#  gtsummary::modify_header(estimate = "**b**")   |> # Sorgt dafür, dass die b's nicht Beta heissen, sondern b
  gtsummary::add_glance_source_note()

Characteristic	Beta	95% CI¹	p-value	GVIF¹	Adjusted GVIF^2,1
Sex				1.5	1.2
female	—	—
male	-0.32	-0.45, -0.18	<0.001
Age_z	0.00	-0.01, 0.00	0.3	3.1	1.8
I(Age_z^2)	0.00	0.00, 0.00	0.7	2.5	1.6
Survived	-0.08	-0.22, 0.07	0.3	1.7	1.3
Pclass_f				1.4	1.1
3	—	—
2	0.49	0.35, 0.63	<0.001
1	1.8	1.6, 1.9	<0.001
Kinder				3.5	1.9
FALSE	—	—
TRUE	0.61	0.27, 0.95	<0.001
R² = 0.578; Adjusted R² = 0.572; Sigma = 0.641; Statistic = 103; p-value = <0.001; df = 7; Log-likelihood = -518; AIC = 1,054; BIC = 1,092; Deviance = 217; Residual df = 528; No. Obs. = 536
¹ CI = Confidence Interval, GVIF = Generalized Variance Inflation Factor
² GVIF^[1/(2*df)]

Voraussetzungschecks

Übung: Suchen Sie die Befehle und führen Sie sie aus:

Die möglichen Verletzungen der Voraussetzungen sind:

Multikollinearität
Deutliche Abweichung von der Normalverteilung in den Fehlern
Heteroskedastizität
Outlier

Algorithmus für den Fahrpreis

Prognosewerte erstellen:

R-Code

preds <- predict(fit_titanic, test) 

modelEval <- cbind(log(test$Fare), preds) |> 
  as_tibble() 

# Beschrifte die Spalten
colnames(modelEval) <- c('actual', 'predicted')

plot <- ggplot(modelEval, aes(x = predicted, y = actual)) + 
  # Create a diagonal line:
  geom_abline(lty = 2) + 
  geom_point(alpha = 0.5) + 
  labs(y = "Tatsächliche Werte", x = "Vorhergesagte Werte") 

plot

Überlebensprognose

Zweite Analyse: Überlebensprognose (Dummy) mit linearer Regression. Was denken Sie?

DATEN_titanic <- readRDS("data/titanic/train.RDS")

model1 <- lm(Survived ~ Pclass_f + Sex + Age_z + I(Age_z^2) + Kinder, data=train)

broom::tidy(model1)|> gt::gt()

term	estimate	std.error	statistic	p.value
(Intercept)	0.57899125672	0.03753055415	15.4271971	0.00000000000000000000000000000000000000000001312827
Pclass_f2	0.17998008082	0.04200295904	4.2849381	0.00002171949705443183080400514073371454060179530643
Pclass_f1	0.37609658211	0.04287551640	8.7718263	0.00000000000000002416214395714002410899718603234220
Sexmale	-0.50783615340	0.03528670506	-14.3917136	0.00000000000000000000000000000000000000072310102002
Age_z	-0.00292748404	0.00200018508	-1.4636066	0.14389540567839331774990796475321985781192779541016
I(Age_z^2)	-0.00001618071	0.00008789737	-0.1840864	0.85401623963406969597400575366918928921222686767578
KinderTRUE	0.12521089068	0.10287434720	1.2171245	0.22409953431140666735288391464564483612775802612305

broom::glance(model1)|> gt::gt()

r.squared	adj.r.squared	sigma	statistic	p.value	df	logLik	AIC	BIC	deviance	df.residual	nobs
0.4032852	0.3965172	0.3819574	59.58678	0.00000000000000000000000000000000000000000000000000000002851803	6	-241.1569	498.3138	532.5869	77.1766	529	536

Logistische Regression

term	estimate	std.error	statistic	p.value
(Intercept)	0.2353041	0.2209270	1.065076	0.2868414000117959772317988154100021
Pclass_f2	1.1302165	0.2863985	3.946307	0.0000793659122583148654225815299412
Pclass_f1	2.1581173	0.2870682	7.517787	0.0000000000000557111855943479295124
Sexmale	-2.7026219	0.2437692	-11.086808	0.0000000000000000000000000001453835
KinderTRUE	1.1726858	0.3732106	3.142155	0.0016770912214931204143458787214627

null.deviance	df.null	logLik	AIC	BIC	deviance	df.residual	nobs
724.2873	535	-241.9354	493.8708	515.2915	483.8708	531	536

Kennwerte

Übung: Interpretieren Sie den Output

R-Code

model2 <- glm(Survived ~ Pclass_f + Sex +  Kinder,  family=binomial, data=train)

sjPlot::tab_model(model2,title = "Überlebensanalyse zum Titanicunglück mit sjPlot",
                  show.est = TRUE, # zeige die estimates
                  show.r2 = TRUE, # zeige R^2
                  show.fstat = TRUE # zeige die F-Statistik
                  )

Überlebensanalyse zum Titanicunglück mit sjPlot
	Survived
Predictors	Odds Ratios	CI	p
(Intercept)	1.27	0.82 – 1.96	0.287
Pclass f [2]	3.10	1.77 – 5.47	<0.001
Pclass f [1]	8.65	4.99 – 15.41	<0.001
Sex [male]	0.07	0.04 – 0.11	<0.001
KinderTRUE	3.23	1.56 – 6.78	0.002
Observations	536
R² Tjur	0.404

Wieder auch mit `gtsummary::`

R-Code

model2 |> 
  gtsummary::tbl_regression(exponentiate = TRUE) |> 
  gtsummary::add_vif() |> 
#  gtsummary::add_global_p() |> # damit könnte man die Signifikanz ganzer Faktoren testen, statt jede Ausprägung gegen die Referenz
  gtsummary::add_glance_source_note() |> 
  gtsummary::modify_caption("**Überlebensanalyse zum Titanicunglück mit gtsummary**")

**Überlebensanalyse zum Titanicunglück mit gtsummary**
Characteristic	OR¹	95% CI¹	p-value	GVIF¹	Adjusted GVIF^2,1
Pclass_f				1.2	1.0
3	—	—
2	3.10	1.77, 5.47	<0.001
1	8.65	4.99, 15.4	<0.001
Sex				1.1	1.0
female	—	—
male	0.07	0.04, 0.11	<0.001
Kinder				1.1	1.0
FALSE	—	—
TRUE	3.23	1.56, 6.78	0.002
Null deviance = 724; Null df = 535; Log-likelihood = -242; AIC = 494; BIC = 515; Deviance = 484; Residual df = 531; No. Obs. = 536
¹ OR = Odds Ratio, CI = Confidence Interval, GVIF = Generalized Variance Inflation Factor
² GVIF^[1/(2*df)]

Voraussetzungschecks

Übung: Führen Sie die Checks für die Voraussetzungen aus!

Vorhersagetest

Was sagt dieser Test?

preds <- predict(model2, test) 

modelEval <- cbind(test$Survived, preds) |> 
  as.tibble() |> 
  mutate(preds = ifelse(preds > 0.5,1,0))

# Beschrifte die Spalten
colnames(modelEval) <- c('actual', 'predicted')

modelEval |> 
sjmisc::flat_table()

   predicted   0   1

actual
0 100 6 1 31 41


modelEval |> 
  filter(!is.na(predicted)) |> 
  mutate(test = actual == predicted, na.rm = TRUE) |> 
  summarise(Accuracy = mean(test))

3 A tibble: 1 × 1

Accuracy 1 0.792

Take Home – Ausblick – Vokabeln

Take Home

Logistische Regression

Sie können eine logistische Regression in R rechnen und interpretieren
Sie können Machine Learning in R anwenden

Ausblick

Sie lernen Clusteranalysen kennen. Wie man also Gruppen anhand von Eigenschaften identifiziert

Statistik und Datenanalyse: Aufbau

Orga

Lernziele

1 Logistische Regression

1.1 Grund und Folgen

RAQ-ML-«S-Kurve» sieht gerade aus

1.2 Multinominale Regression

2 Machine Learning – ML

Supervised

Übung 4

Daten einlesen

Daten in Trainings- und Testdaten aufteilen

Datenvisualisierung

Modellbildung für den Fahrpreis

Regressionsoutput mit sjPlot::tab_model

Regressionsoutput mit gtsummary::tbl_regression

Voraussetzungschecks

Algorithmus für den Fahrpreis

Überlebensprognose

Logistische Regression

Kennwerte

Wieder auch mit gtsummary::

Voraussetzungschecks

Vorhersagetest

3 A tibble: 1 × 1

Take Home – Ausblick – Vokabeln

Take Home

Ausblick

Regressionsoutput mit `sjPlot::tab_model`

Regressionsoutput mit `gtsummary::tbl_regression`

Wieder auch mit `gtsummary::`