Stat4: Lösung 4.2S

Veröffentlichungsdatum

7. November 2023

Musterlösung Übung 4.2S: Multiple logistische Regression (SozWis)

Kommentierter Lösungsweg

df <- read_delim("datasets/stat1-4/Datensatz_novanimal_Uebung_Statistik4.2.csv", delim = ";")

# sieht euch die Verteilung zwischen Mensagänger und Selbstverpfleger an
# sind nicht gleichmässig verteilt, bei der Vorhersage müssen wir das berücksichtigen
table(df$mensa)

  0   1 
282 786 
df |> count(mensa) # alternativ
# A tibble: 2 × 2
  mensa     n
  <dbl> <int>
1     0   282
2     1   786
# definiert das logistische Modell und wendet es auf den Datensatz an

mod0 <- glm(mensa ~ gender + member + age_groups + meat + umwelteinstellung,
  data = df, binomial("logit")
)
summary.lm(mod0) # Umwelteinstellung scheint keinen Einfluss auf die

Call:
glm(formula = mensa ~ gender + member + age_groups + meat + umwelteinstellung, 
    family = binomial("logit"), data = df)

Weighted Residuals:
    Min      1Q  Median      3Q     Max 
-5.6740 -0.8078  0.3712  0.5867  1.2379 

Coefficients:
                             Estimate Std. Error t value Pr(>|t|)    
(Intercept)                  -0.18889    0.40225  -0.470 0.638750    
genderMann                    0.71017    0.16018   4.434 1.02e-05 ***
memberStudent/in             -0.63072    0.29442  -2.142 0.032404 *  
age_groups26- bis 34-jaehrig  1.09429    0.19574   5.591 2.88e-08 ***
age_groups35- bis 49-jaehrig  1.75379    0.45968   3.815 0.000144 ***
age_groups50- bis 64-jaehrig  2.43530    0.78923   3.086 0.002083 ** 
meat                          0.19945    0.05055   3.945 8.49e-05 ***
umwelteinstellung             0.19334    0.18688   1.035 0.301107    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.009 on 1060 degrees of freedom
Multiple R-squared:  0.004042,  Adjusted R-squared:  -0.002536 
F-statistic: 0.6145 on 7 and 1060 DF,  p-value: 0.7443
# Verpflegung zu haben, gegeben die Daten

# neues Modell ohne Umwelteinstellung
mod1 <- update(mod0, ~ . -umwelteinstellung)
summary.lm(mod1)

Call:
glm(formula = mensa ~ gender + member + age_groups + meat, family = binomial("logit"), 
    data = df)

Weighted Residuals:
    Min      1Q  Median      3Q     Max 
-6.0117 -0.8060  0.3584  0.6100  1.2407 

Coefficients:
                             Estimate Std. Error t value Pr(>|t|)    
(Intercept)                   0.03212    0.34053   0.094 0.924860    
genderMann                    0.69697    0.15951   4.369 1.37e-05 ***
memberStudent/in             -0.64418    0.29426  -2.189 0.028806 *  
age_groups26- bis 34-jaehrig  1.11651    0.19458   5.738 1.25e-08 ***
age_groups35- bis 49-jaehrig  1.77409    0.45947   3.861 0.000120 ***
age_groups50- bis 64-jaehrig  2.44683    0.78953   3.099 0.001992 ** 
meat                          0.18070    0.04709   3.837 0.000132 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.01 on 1061 degrees of freedom
Multiple R-squared:  0.003998,  Adjusted R-squared:  -0.001635 
F-statistic: 0.7098 on 6 and 1061 DF,  p-value: 0.6418
# Modeldiagnostik (wenn nicht signifikant, dann OK)
1 - pchisq(mod1$deviance, mod1$df.resid) # Ok
[1] 0.4509591
# Modellgüte (pseudo-R²)
1 - (mod1$dev / mod1$null) # eher kleines pseudo-R2, deckt sich mit dem R-Squared aus dem obigen output summary.lm()
[1] 0.1354244
# Konfusionsmatrix vom  Datensatz
# Model Vorhersage
# hier ein anderes Beispiel:
predicted <- predict(mod1, df, type = "response")

# erzeugt eine Tabelle mit den beobachteten
# Mensagänger/Selbstverpfleger und den Vorhersagen des Modells
km <- table(predicted > 0.5, df$mensa)
# alles was höher/grosser ist als 50% ist
# kommt in die Kategorie Mensagänger

# anpassung der namen
dimnames(km) <- list(
  c("Modell Selbstverpfleger", "Modell Mensagänger"),
  c("Daten Selbstverpfleger", "Daten Mensagänger")
)
km
                        Daten Selbstverpfleger Daten Mensagänger
Modell Selbstverpfleger                     87                59
Modell Mensagänger                         195               727
#############
### reminder: https://towardsdatascience.com/understanding-confusion-matrix-a9ad42dcfd62
#############

# TP = true positive: you predicted positive and it’s true; hier Vorhersage
# Mensagänger stimmt also (727)

# TN = true negative: you predicted negative and it’s true, hier Vorhersage der
# Selbstverpfleger stimmt (87)

# FP = false positive (fehler 1. art, auch Spezifizität genannt) you predicted
# and it’s false. hier Modell sagt Mensagänger vorher
# (obwohl in Realität Selbstverpfleger) (195)

# FN = false negative (fehler 2. art, auch Sensitivität genannt),
# you predicted negative and it’s false. hier Modell sagt Selbtverpfleger vorher
# (obwohl in Realität Mensagänger) (59)


# es scheint, dass das Modell häufig einen alpha Fehler macht, d.h.
# das Modell weist keine hohe Spezifizität auf: konkret werden viele Mensagänger als
# Selbstverpfleger vorhergesagt resp. klassifiziert. Dafür gibt es mehere Gründe:

# 1) die Kriteriumsvariable ist sehr ungleich verteilt, d.h. es gibt weniger
# Selbstverpfleger als Mensgänger im Datensatz

# 2) nicht adäquates Modell z.B. link mit probit zeigt besserer fit

# 3) Overfitting: wurde hier nicht berücksichtigt, in einem Paper/Arbeit
# müsste noch einen Validierungstest gemacht werden z.B. test-train
# Cross-Validation oder k fold Cross-Validation

# kalkuliert die Missklassifizierungsrate
mf <- 1 - sum(diag(km) / sum(km)) # ist mit knapp 23 %  eher hoch
mf
[1] 0.2378277
# kleiner exkurs: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2636062/
# col wise proportion, da diese die "Realität" darstellt
km_prop <- prop.table(km, 2)

# specificity = a / (a+c) => ability of a test to correctly
spec = km_prop[1] / (km_prop[1] + km_prop[2])
spec
[1] 0.3085106
# sensitivity = d / (b+d) => Sensitivity is the ability of a
sens = km_prop[4] / (km_prop[3] + km_prop[4])
sens
[1] 0.9249364

Methode

In dieser Aufgabe möchten wir untersuchen, ob wir vorhersagen können, ob jemand die Mensa besuchen wird. Wir interessieren uns dabei besonders für sozioökonomische Faktoren, die wahrgenommene Fleischkonsum und die Einstellung zur Umwelt. Die Zielgröße “Mensa-Besuch” ist entweder ja oder nein, daher führen wir eine “multiple logistische Regression” durch. Wir verwenden die Merkmale “Alter”, “Geschlecht”, “Hochschulzugehörigkeit”, “Fleischkonsum” und “Umwelteinstellung” als Vorhersagevariablen.

Wenn du mehr über logistische Regressionen erfahren möchtest, findest du ausführlichere Informationen in den Büchern von Crawley (2015) und Gareth (2016).

Ergebnisse

Daten Selbstverpfleger Daten Mensagänger
Modell Selbstverpfleger 87 59
Modell Mensagänger 195 727

Konfusionsmatrix

Die Ergebnisse des logistischen Modells mit der “logit” Linkfunktion deuten darauf hin, dass das Modell nicht gut zu unseren Daten passt. Das bedeutet, dass es schwierig ist, basierend auf diesem Modell vorherzusagen, ob jemand in Zukunft in der Mensa essen wird oder sein Mittagessen mitbringt. Diese Schlussfolgerung wird durch den kleinen pseudo-R2-Wert (14%) und die hohe Rate der Fehlklassifikation (24%) gestützt. Wenn wir genauer hinsehen, stellen wir fest, dass das Modell häufig einen Fehler vom Typ “Alpha” (1. Fehler) macht, das bedeutet, es sagt oft voraus, dass jemand in die Mensa geht, obwohl sie tatsächlich ihr Essen selbst mitbringen.

Es gibt mehrere Gründe für die schlechte Passung des Modells:

Die Kriteriumsvariable ist in den Daten ungleich verteilt, es gibt also weniger Personen, die ihr Essen selbst mitbringen (26%) im Vergleich zu denen, die in der Mensa essen (74%).

Die Prädiktorvariablen sind alle entweder kategorisch oder ordinal, was dazu führen kann, dass das Modell keine guten Ergebnisse liefert.

Zusammengefasst, es wäre ratsam, nach einem besseren geeigneten Modell zu suchen, insbesondere einem Modell, das mit ordinalen Prädiktorvariablen besser umgehen kann. Beispiele wären:

  • eine andere Link-Funktion für das GLM verwenden z.B. probit
  • polynomiale Kontraste
  • Smooth Splines hier
  • multinomiale Regression z.M. nnet::mulitom() hier