Diskriminanzanalyse

Einleitung

Die (deskriptive) Diskriminanzanalyse geht der entgegengesetzten Fragestellung der MANOVA auf den Grund. Mit ihr können wir (deskriptiv) untersuchen, ob Gruppenzugehörigkeiten durch die AVs der MANOVA vorhergesagt werden können (siehe bspw. Pituch und Stevens, 2016, Kapitel 10 sowie Eid, Gollwitzer & Schmitt, 2017, Kapitel 15.4). Wir wollen uns wieder das fiktive Datenbeispiel (Datensatz Therapy aus dem gleichnamigen .rda File Therapy.rda) ansehen, den wir bereits in der MANOVA-Sitzung untersucht haben. Sie können den Datensatz “Therapy.rda” hier herunterladen.

Daten laden

Wir laden zunächst die Daten: entweder lokal von Ihrem Rechner:

load("C:/Users/Musterfrau/Desktop/Therapy.rda")

oder wir laden sie direkt über die Website:

load(url("https://pandar.netlify.app/daten/Therapy.rda"))

Nun sollte in R-Studio oben rechts in dem Fenster unter der Rubrik “Data” unser Datensatz mit dem Namen “Therapy” erscheinen.

Übersicht über die Daten

Wir wollen uns einen Überblick über die Daten verschaffen:

head(Therapy)

##   Lebenszufriedenheit Arbeitsbeanspruchung Depressivitaet Arbeitszufriedenheit
## 1                   7                    4              7                    5
## 2                   5                    5              8                    3
## 3                   8                    7              6                    6
## 4                   6                    4              5                    5
## 5                   6                    9              8                    5
## 6                   8                    7              8                    6
##     Intervention Geschlecht
## 1 Kontrollgruppe          0
## 2 Kontrollgruppe          1
## 3 Kontrollgruppe          0
## 4 Kontrollgruppe          1
## 5 Kontrollgruppe          1
## 6 Kontrollgruppe          1

levels(Therapy$Intervention)

## [1] "Kontrollgruppe"              "VT Coaching"                
## [3] "VT Coaching + Gruppenuebung"

levels(Therapy$Geschlecht)

## [1] "0" "1"

Die abhängigen Variablen sind Lebenszufriedenheit, Arbeitsbeanspruchung Depressivitaet und Arbeitszufriedenheit. Die Variable Intervention hat drei Stufen: eine Kontrollgruppe, ein verhaltenstherapiebasiertes Coaching, sowie das verhaltenstherapiebasierte Coaching inklusive einer Gruppenübung. Das Geschlecht ist 0-1 kodiert, wobei 0 für männlich und 1 für weiblich steht. Insgesamt sind die Variablennamen der AVs recht lang. Wir wollen diese kürzen:

colnames(Therapy) # Spaltennamen ansehen

## [1] "Lebenszufriedenheit"  "Arbeitsbeanspruchung" "Depressivitaet"      
## [4] "Arbeitszufriedenheit" "Intervention"         "Geschlecht"

colnames(Therapy) <- c("LZ", "AB", "Dep", "AZ", "Intervention", "Geschlecht") # Spaltennamen neu zuordnen
head(Therapy)

##   LZ AB Dep AZ   Intervention Geschlecht
## 1  7  4   7  5 Kontrollgruppe          0
## 2  5  5   8  3 Kontrollgruppe          1
## 3  8  7   6  6 Kontrollgruppe          0
## 4  6  4   5  5 Kontrollgruppe          1
## 5  6  9   8  5 Kontrollgruppe          1
## 6  8  7   8  6 Kontrollgruppe          1

So - schon viel übersichtlicher!

Pakete laden

Nachdem wir neue Pakete installiert haben (install.packages), laden wir diese:

library(MASS) # für lineare Diskrimianzanalys (lda)
library(ggplot2) # Grafiken

Ziel der linearen Diskriminanzanalyse

Wir wollen mit der Lebenszufriedenheit, der Depression, der Arbeitsbeanspruchung und der Arbeitszufriedenheit die Zugehörigkeit zu den Therapiegruppen vorhersagen (ob dies inhaltlich sinnvoll ist, sei jetzt mal dahingestellt!). Die lineare Diskriminanzanalyse (LDA) funktioniert ähnlich wie die PCA. Die PCA maximiert die Varianz auf jeder Hauptkomponente. Die LDA maximiert die Diskrimination entlang der Achsen, also quasi die Varianz zwischen den Gruppen auf den Achsen. Somit sind beide Verfahren de facto Varianzmaximierungsverfahren.

Die LDA hat im Grunde die gleichen Annahmen wie die MANOVA. Die Kovarianzmatrizen über die Gruppen müssen homogen sein. Die Daten sollten multivariat normalverteilt sein (diese Annahme lässt sich bspw. in Fisher’s linearer Diskriminanzanalye lockern, sodass sie nicht so wichtig ist). Die Beobachtungen stammen aus einer independent and identically distributed ($i.i.d.$, deutsch: unabhängig und identisch verteilt) Population (dies bedeutet, dass alle Beobachtungen unabhängig sind und den gleichen Verteilungs- und Modellannahmen unterliegen). Die letzte Annahme können wir nicht prüfen. Sie kann nur über die sinnvolle Wahl des Designs (Randomisierung etc.) angenommen werden. Wie dies zu prüfen ist, hatten wir in der Sitzung zur MANOVA kennengelernt.

Analysen

Der Befehl in R heißt lda für lineare Diskriminanzanalyse. Ihr übergeben wir wieder eine Gleichung, ähnlich der MANOVA-Gleichung der vergangenen Sitzung. Allerdings ist diesmal die Intervention die AV und die AVs der MANOVA (LZ, Dep, AB und AZ) sind die UVs, die zur Prädiktion herangezogen werden. Mit Hilfe von $scaling können wir die Gewichtungskoeffizienten der Variablen auf den Diskriminanzachsen einsehen. Insgesamt sind bei einer linearen Diskriminanzanalyse immer maximal $\min(p, k-1)$ Diskriminanzfunktionen/achsen möglich, wobei $p$ die Anzahl an UVs ist und $k$ die Anzahl an Gruppen beschreibt. Hier sind es $p=4$ UVs und $k=3$ Gruppen. Folglich können maximal $2$ Diskriminanzfunktionen bestimmt werden:

model_DA <- lda(Intervention ~ LZ + Dep + AB + AZ, Therapy)
model_DA

## Call:
## lda(Intervention ~ LZ + Dep + AB + AZ, data = Therapy)
## 
## Prior probabilities of groups:
##              Kontrollgruppe                 VT Coaching VT Coaching + Gruppenuebung 
##                   0.3333333                   0.3333333                   0.3333333 
## 
## Group means:
##                                   LZ      Dep       AB       AZ
## Kontrollgruppe              5.933333 7.133333 6.033333 5.133333
## VT Coaching                 5.933333 5.066667 5.833333 6.833333
## VT Coaching + Gruppenuebung 7.333333 4.766667 6.300000 7.500000
## 
## Coefficients of linear discriminants:
##            LD1         LD2
## LZ   0.1618370 -0.88229519
## Dep  0.7453510  0.09824489
## AB  -0.4155893 -0.35778884
## AZ  -0.1519953  0.48104478
## 
## Proportion of trace:
##    LD1    LD2 
## 0.7892 0.2108

model_DA$scaling # Koeffizienten

##            LD1         LD2
## LZ   0.1618370 -0.88229519
## Dep  0.7453510  0.09824489
## AB  -0.4155893 -0.35778884
## AZ  -0.1519953  0.48104478

Für die erste Achse sind besonders Dep und AB relevant, während auf der zweiten Achse AZ und LZ stärker ins Gewicht fallen. Wir können predict nutzen, um eine Vorhersage mittels dieses Modells (model_DA) zu erhalten. predict gibt hierbei 3 Listen aus: $posterior gibt die Wahrscheinlichkeit an, in die jeweilige Gruppe zu gelangen, $class gibt die Vorhersage für die Gruppenzugehörigkeit an und $x gibt die vorhergesagte Ausprägung auf der jeweiligen Diskriminanzachse an.

head(predict(model_DA)$posterior) # Wahrscheinlichkeit in der jeweiligen Gruppe zu landen

##   Kontrollgruppe VT Coaching VT Coaching + Gruppenuebung
## 1      0.9518348  0.02825584                 0.019909382
## 2      0.9669084  0.02447494                 0.008616686
## 3      0.2868306  0.08677168                 0.626397717
## 4      0.4461846  0.36122866                 0.192586789
## 5      0.5512585  0.12202873                 0.326712765
## 6      0.8868980  0.02379371                 0.089308278

head(predict(model_DA)$class)     # Vorhergesagte Klasse

## [1] Kontrollgruppe              Kontrollgruppe              VT Coaching + Gruppenuebung
## [4] Kontrollgruppe              Kontrollgruppe              Kontrollgruppe             
## Levels: Kontrollgruppe VT Coaching VT Coaching + Gruppenuebung

head(predict(model_DA)$x)         # vorhergesagte Ausprägung auf der jeweiligen Diskriminanzachse

##         LD1        LD2
## 1 2.1797562 -0.3780597
## 2 2.4898346  0.1648972
## 3 0.1974790 -1.9509215
## 4 0.5272172  0.3077457
## 5 0.6853238 -1.1864638
## 6 1.6881810 -1.7544317

Wir wollen uns die Ausprägungen auf den beiden Diskriminanzachsen grafisch ansehen. Dazu speichern wir diese unter DA1 und DA2 jeweils für 1. und 2. Diskriminanzachse ab.

Therapy$DA1 <- predict(model_DA)$x[, 1] # erste DA
Therapy$DA2 <- predict(model_DA)$x[, 2] # zweite DA

Anschließend können wir wieder ggplot verwenden, um die beiden Diskriminanzachsen gegeneinander abzutragen. Außerdem fügen wir noch die Nullpunkte auf den Diskriminanzachsen als gepunktete Linie ein:

ggplot(data = Therapy, aes(x = DA1, y = DA2, color = Intervention)) +
  geom_point()+
  geom_hline(yintercept = 0, lty = 3)+
  geom_vline(xintercept = 0, lty = 3)+
  ggtitle(label = "Diskriminanzachsen", subtitle = "mit Trennlinien")

Was wir nun sehen, ist, dass entlang der x-Achse besonders zwischen blau/grün vs rot unterschieden wird. D.h. die erste Diskriminanzachse hilft uns, zwischen den Interventionsgruppen und der Kontrollgruppe zu unterscheiden. Die zweite Achse trennt eher zwischen den beiden Interventionsgruppen, wobei diese Trennung nicht sehr eindeutig ist.

Wir hätten auch einfach die plot-Funktion auf das model_DA-Objekt anwenden können:

plot(model_DA)

Wenn wir diesem noch Farben zuordnen (hier jeweils 30 mal die selbe Farbe, da die Gruppen so häufig hintereinander im Datensatz standen…), dann sieht diese Grafik der mit ggplot erzeugten Grafik recht ähnlich. Natürlich können wir auch hier die Nullpunkte einfügen:

plot(model_DA, col = c(rep("red", 30), rep("gold3", 30), rep("blue", 30)))
abline(v = 0, lty = 3)
abline(h = 0, lty = 3)

Zu guter Letzt fügen wir noch die Trennlinien, die zwischen den Gruppen unterscheiden, in die erste Grafik ein. Diese liegen immer auf halbem Weg (durchgezogenen schwarzen Linien, hier für Kontrollgruppe vs. VT-Coaching + Gruppenübung) zwischen zwei Mittelwerten pro Gruppe (gestrichelten Linien) auf den jeweiligen Achsen:

# Mittelwerte auf den DAs
Means <- aggregate(cbind(DA1, DA2) ~ Intervention, data = Therapy, FUN = mean)

# Mittelwerte auf DA1
mDA1_K <- Means[1,2] # Kontrollgruppenmittelwert auf DA1
mDA1_V <- Means[2,2] # Mittelwert VT auf DA1
mDA1_VG <- Means[3,2] # Mittelwert VT + Gruppenuebung auf DA1

# Mittelwerte auf DA2
mDA2_K <- Means[1,3] # Kontrollgruppenmittelwert auf DA2
mDA2_V <- Means[2,3] # Mittelwert VT auf DA2
mDA2_VG <- Means[3,3] # Mittelwert VT + Gruppenuebung auf DA2

ggplot(data = Therapy, aes(x = DA1, y = DA2, color = Intervention)) + geom_point()+
        geom_hline(yintercept = mDA2_K, lty = 2, col = "red")+
        geom_hline(yintercept = mDA2_V, lty = 2, col = "gold3")+
        geom_hline(yintercept = mDA2_VG, lty = 2, col = "blue")+
        geom_hline(yintercept = (mDA2_VG+mDA2_K)/2, lty = 1, col = "black", lwd = 0.2)+
        geom_vline(xintercept = mDA1_K, lty = 2, col = "red")+
        geom_vline(xintercept = mDA1_V, lty = 2, col = "gold3")+
        geom_vline(xintercept = mDA1_VG, lty = 2, col = "blue")+
        geom_vline(xintercept = (mDA1_VG+mDA1_K)/2, lty = 1, col = "black", lwd = 0.2)+
        ggtitle(label = "Diskriminanzachsen", subtitle = "mit Mittelwerten pro Gruppe")

Wir sehen auch an den mit plot erstellten Grafiken, dass sich die Kontrollgruppe von den Interventionsgruppen unterscheiden lässt, während sich die Interventionsgruppen nicht stark voneinander unterscheiden. Dies hatten wir insbesondere auch in der Sitzung zur MANOVA so herausgefunden.

Wie gut ist unsere Gruppenzuordnung?

Um zu prüfen, wie gut wir die Gruppen zugeordnet haben, müssen wir untersuchen, wie häufig richtig zugeordnet wurde. Das können wir bspw. mit table machen und hier die Vorhersagen (predict(model_DA)$class) mit den Beobachtungen (Therapy$Intervention) vergleichen. Wir speichern predict(model_DA)$class zunächst als Variable ab und hängen es an den Datensatz dran:

Therapy$predict_class <- predict(model_DA)$class
table(Therapy$predict_class, Therapy$Intervention)

Die Diagonalelemente dieser Tabelle zeigen die richtig zugeordneten Werte. Wir können diese Werte auch nochmals durch 30 teilen, was die relativen Anteile erzeugt (dies funktioniert hier so einfach, da in jeder Gruppe 30 Probanden waren, normalerweise müssten wir hier die relativen Anteile pro Gruppe betrachten).

table(Therapy$predict_class, Therapy$Intervention)/30

83.33% der Kontrollgruppenprobanden wurden dieser wieder richtig zugeordnet, 50% der VT-Coaching- Gruppe wurden dieser wieder richtig zugeordnet und 66.67% der VT-Coaching + Gruppenübung- Gruppe wurden dieser wieder richtig zugeordnet. Insgesamt wurden

mean(Therapy$predict_class == Therapy$Intervention)*100

## [1] 66.66667

% der Probanden der richtigen Gruppe zugeordnet. Dies zeigt, dass die Prädiktion nicht perfekt war, was höchstwahrscheinlich daran liegt, dass Proband*innen aus den Interventionsgruppen sich kaum in der Lebenszufriedenheit, Arbeitszufriedenheit, Depression und Arbeitsbeanspruchung unterschieden. Ähnliches hatten wir bereits mit der MANOVA herausgefunden.

Trennlinien im ursprünglichen Variablenkoordinatensystem

Nun wollen wir die Trennlinien auch einmal in das ursprüngliche Koordinatensystem einzeichnen. Dazu führen wir eine Diskriminanzanalyse nur mit den beiden Zufriedenheitsmaßen durch.

model_DA2 <- lda(Intervention ~ LZ + AZ, data = Therapy)
model_DA2

## Call:
## lda(Intervention ~ LZ + AZ, data = Therapy)
## 
## Prior probabilities of groups:
##              Kontrollgruppe                 VT Coaching VT Coaching + Gruppenuebung 
##                   0.3333333                   0.3333333                   0.3333333 
## 
## Group means:
##                                   LZ       AZ
## Kontrollgruppe              5.933333 5.133333
## VT Coaching                 5.933333 6.833333
## VT Coaching + Gruppenuebung 7.333333 7.500000
## 
## Coefficients of linear discriminants:
##           LD1       LD2
## LZ -0.3160376  0.883895
## AZ  0.8736389 -0.337895
## 
## Proportion of trace:
##    LD1    LD2 
## 0.7589 0.2411

model_DA2$scaling # Koeffizienten

##           LD1       LD2
## LZ -0.3160376  0.883895
## AZ  0.8736389 -0.337895

Wir sehen, dass jeweils eine Variable besonders stark auf einer Diskriminanzfunktion diskriminiert (großer Koeffizient). Um jetzt die Trennlinien einzeichnen zu können, brauchen wir ein Koordinatensystem, in dem für alle Kombinationen von LZ und AZ entschieden wird, in welcher Gruppe ein/e Proband*in mit dieser Ausprägung landen würde. Dies übernimmt die expand.grid Funktion für uns.

# Ein Koordinatensystem erstellen von 0 bis 12 auf den beiden Variablen
contour_data <- expand.grid(LZ = seq(0,12, 0.01), AZ = seq(0,12,0.01))
contour_data

$\vdots$

Die predict-Funktion sagt die modellimplizierten Werte vorher. Wir können dieser Funktion auch einen Wert übergeben, der so in den Daten noch nicht vorkam, nämlich dem Argument newdata. Somit bekommen wir die Vorhersage unter unserem Modell für die neuen Datenpunkte. Als letztes kann mit stat_contour die Gruppenzugehörigkeit eingezeichnet werden. Diese übergeben wir der 3. Dimension z in der stat_contour-Funktion. Bei diesem Plot ist es von Vorteil, wenn die Gruppierungsvariable im Originaldatensatz gleich heißt wie die Prädiktion im neuen Datensatz.

# Für das Koordinatensystem für jeden Punkt die Gruppenzugehörigkeit bestimmen
contour_data$Intervention <- as.numeric(predict(object = model_DA2, newdata = contour_data)$class)

head(contour_data$Intervention)

## [1] 1 1 1 1 1 1

# Gruppenzugehörigkeiten in Originalkoordinatensystem einzeichnen
ggplot(data = Therapy, mapping = aes(x = LZ, y = AZ, color = Intervention))+
        geom_point()+
        stat_contour(aes(x = LZ, y = AZ, z = Intervention), data = contour_data)+
        ggtitle("Lebenszufriedenheit vs Arbeitszufriedenheit", subtitle = "inklusive retransformierter Entscheidungslinien\nabgeleitet von den Diskriminanzachsen")

## Warning: The following aesthetics were dropped during statistical transformation: colour
## ℹ This can happen when ggplot fails to infer the correct grouping structure in the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical variable into
##   a factor?

Wir sehen, dass die Trennlinien nicht rechtwinklig zueinander verlaufen. Dies liegt daran, dass die Achsen hier in die Originalskala retransformiert wurden. Die Gruppen werden wie folgt zugeordnet: oben links = VT Coaching, oben rechts = VT Coaching + Gruppenuebung und unten rechts/unten = Kontrollgruppe. Allerdings ist diese Zuordnung nicht sehr genau…

Literatur

Eid, M., Gollwitzer, M., & Schmitt, M. (2017). Statistik und Forschungsmethoden (5. Auflage, 1. Auflage: 2010). Weinheim: Beltz.

Pituch, K. A. & Stevens, J. P. (2016). Applied Multivariate Statistics for the Social Sciences (6th ed.). New York: Taylor & Francis.

• Blau hinterlegte Autor:innenangaben führen Sie direkt zur universitätsinternen Ressource.