Kontraste_mehrfaktoriell
Neele Engelmann
22 Juni 2020
Daten einlesen und Pakete laden
#Daten
library(readr)
data <- read_delim("Glastonbury.csv", ";", escape_double = FALSE, trim_ws = TRUE)
#Commander
library(Rcmdr)“music” zum Faktor definieren
# auf richtige Schreibweise der levels achten (groß, klein, Leerzeichen)
data$music <- factor(data$music,levels=c("Crusty","Indie Kid","Metaller","No Musical Affiliation"))
# prüfen, ob levels richtig definiert wurden und in welcher Reihenfolge sie sind
levels(data$music)## [1] "Crusty" "Indie Kid" "Metaller"
## [4] "No Musical Affiliation"Grafik erstellen im commander (Kriterium: day1)
with(data, plotMeans(day1, music, error.bars="conf.int", level=0.95, connect=TRUE))
Helmert-Kontraste anlegen im Commander
contrasts(data$music) <- "contr.helmert"Prüfen, welche Kontrastkoeffizienten dadurch angelegt wurden:
contrasts(data$music)## [,1] [,2] [,3]
## Crusty -1 -1 -1
## Indie Kid 1 -1 -1
## Metaller 0 2 -1
## No Musical Affiliation 0 0 3Wann immer ein Kontrastkoeffizient 0 ist, geht der entsprechende Mittelwert nicht in den Vergleich ein. Koeffizienten mit positivem Vorzeichen werden immer mit Koeffizienten mit negativem Vorzeichen verglichen.
- 3 ist also No Musical Affiliation vs. [Crusty, Indie, Metal]
- 2 ist Metal vs. [Crusty, Indie]
- 1 ist Indie vs. Crusty
Da “No Affiliation” sich laut Graph von allen anderen Gruppen zu unterscheiden scheint und auch inhaltlich als eine Art Baseline-Gruppe gesehen werden kann, können wir die Anordnung so lassen. Falls man sich inhaltlich mehr dafür interessiert, ob sich eine der anderen Gruppen von allen anderen unterscheidet, sollte man die Faktorstufen neu anordnen und diese Gruppe an die letzte Stelle setzen.
Mit dem Code, mit dem wir oben auch den character-Vektor zum Faktor gemacht haben, könnten wir ihn auch erneut überschreiben und die Reihenfolge der levels ändern (hier auskommentiert, weil wir gleich mit den eben angelegten levels weiterrechnen wollen):
#Beispiel: "Crusty" nach hinten stellen
#data$music <- factor(data$music,levels=c("No Musical Affiliation","Indie Kid","Metaller","Crusty"))Rechnen des linearen Modells mit den definierten Kontrasten:
KontrastModell_Helmert <- lm(day1 ~ music, data=data)
summary(KontrastModell_Helmert)##
## Call:
## lm(formula = day1 ~ music, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.82536 -0.46906 0.00884 0.45797 2.07715
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.82763 0.02611 69.986 < 2e-16 ***
## music1 0.05747 0.04160 1.381 0.168
## music2 -0.03360 0.02197 -1.529 0.127
## music3 -0.07160 0.01233 -5.806 9.21e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6803 on 806 degrees of freedom
## Multiple R-squared: 0.04139, Adjusted R-squared: 0.03782
## F-statistic: 11.6 on 3 and 806 DF, p-value: 1.896e-07Interpretation des Outputs:
Intercept: Mittelwert der Mittelwerte
Estimates: Beziehen sich auf die durch die Kontrastkoeffizienten kodierten Vergleiche (wenn größer, dann stärkerer Unterschied zwischen den entsprechenden Gruppen). (Sind hier aber nicht identisch mit den entsprechenden Mittelwertsdifferenzen.)
Rest des Modells (Residuen, R-Squared, F-Test): Zu interpretieren wie immer (Vergleich zum Nullmodell)
Interpretation der Estimates im Detail (siehe VL Kontraste Folie 15)
(Dieser Teil ist Zusatz fürs Verständnis, den Code müsst ihr nicht selbst schreiben können.)
Zur Veranschaulichung legen wir erstmal die Kontrast-Dummy-Variablen sichtbar im Datensatz an (nur damit wir sie nochmal sehen können)
#Zur Erinnerung: wie müssen wir die Koeffizienten vergeben, um die oben im Commander angelegten Kontraste zu reproduzieren?
contrasts(data$music)## [,1] [,2] [,3]
## Crusty -1 -1 -1
## Indie Kid 1 -1 -1
## Metaller 0 2 -1
## No Musical Affiliation 0 0 3#Erster Kontrast: -1 wenn Crusty, 1 wenn Indie Kid, ansonsten 0 (die ifelse-Funktion ist sehr ähnlich aufgebaut wie die "wenn"-Funktion in Excel)
data$music1 <- ifelse(data$music=="Crusty",-1,
ifelse(data$music=="Indie Kid",1,0))
#Zweiter Kontrast: -1 wenn Crusty oder Indie Kid, 2 wenn Metaller, sonst 0
data$music2 <- ifelse(data$music=="Crusty"|data$music=="Indie Kid",-1,
ifelse(data$music=="Metaller",2,0))
#Dritter Kontrast: 3 wenn No Musical Affiliation, sonst -1
data$music3 <- ifelse(data$music=="No Musical Affiliation",3,-1)
#Datensatz ansehen
View(data)Kurzer Check: lineares Modell mit den eben angelegten Kontrasten sollte gleiches Ergebnis liefern wie Helmert-Kontraste über commander
summary(lm(day1~music1+music2+music3,data=data))##
## Call:
## lm(formula = day1 ~ music1 + music2 + music3, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.82536 -0.46906 0.00884 0.45797 2.07715
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.82763 0.02611 69.986 < 2e-16 ***
## music1 0.05747 0.04160 1.381 0.168
## music2 -0.03360 0.02197 -1.529 0.127
## music3 -0.07160 0.01233 -5.806 9.21e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6803 on 806 degrees of freedom
## Multiple R-squared: 0.04139, Adjusted R-squared: 0.03782
## F-statistic: 11.6 on 3 and 806 DF, p-value: 1.896e-07Das statistische Modell ist also:
y =beta0 + beta1 x music1 + beta2 x music2 + beta3 x music3
Die Regressionskoeffizienten sind jetzt also Gewichte von unseren Kontrast-Dummy-Variablen (unseren Prädiktoren).
Da wir wissen, welche Werte die Prädiktoren für jede “music”-Gruppe annehmen (wir haben sie ja schließlich selbst vergeben), können wir uns die Mittelwerte jeder Gruppe mit den ermittelten Regressionskoeffizienten ausgeben lassen.
Crusty:
#kurz checken: bei "Crusty" haben alle drei Kontrastvariablen (music1,music2,music3) den Wert -1
head(subset(data,music=="Crusty"))## # A tibble: 6 x 9
## ticknumb music day1 day2 day3 change music1 music2 music3
## <dbl> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2229 Crusty 0.97 1.41 0.290 -0.68 -1 -1 -1
## 2 2384 Crusty 3.03 NA NA NA -1 -1 -1
## 3 2405 Crusty 0.85 NA NA NA -1 -1 -1
## 4 2490 Crusty 2.29 NA NA NA -1 -1 -1
## 5 2510 Crusty 0.82 0.2 0.47 -0.35 -1 -1 -1
## 6 2521 Crusty 2.79 NA NA NA -1 -1 -1#regressionskoeffizienten als objekte speichern (für übersicht, weil sie jetzt noch mehrmals benutzt werden)
beta0 <- 1.82763
beta1 <- 0.05747
beta2 <- -0.03360
beta3 <- -0.07160
#Einsetzen in die Modellgleichung ergibt den Mittelwert im Kriterium (day1) für die Crusty-Gruppe
beta0 + beta1*(-1) + beta2*(-1) + beta3*(-1) ## [1] 1.87536#kürzer
beta0-beta1-beta2-beta3## [1] 1.87536#Mittelwert der Crusty-Gruppe direkt berechnen zum Vergleich:
mean(subset(data,music=="Crusty")$day1)## [1] 1.875361#mittelwert speichern für später
M_Crusty <- mean(subset(data,music=="Crusty")$day1)Indie Kid:
#kurz checken: Werte der Kontrastvariablen bei "Indie Kid"
head(subset(data,music=="Indie Kid"))## # A tibble: 6 x 9
## ticknumb music day1 day2 day3 change music1 music2 music3
## <dbl> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2467 Indie Kid 1.56 NA NA NA 1 -1 -1
## 2 2478 Indie Kid 3.02 NA NA NA 1 -1 -1
## 3 2565 Indie Kid 2.3 NA NA NA 1 -1 -1
## 4 2602 Indie Kid 0.97 0.38 0.76 -0.21 1 -1 -1
## 5 2677 Indie Kid 1.28 0.38 0.14 -1.14 1 -1 -1
## 6 2688 Indie Kid 2.64 NA NA NA 1 -1 -1#Einsetzen in die Modellgleichung ergibt den Mittelwert im Kriterium (day1) für die Indie Kid-Gruppe
beta0 + beta1*(1) + beta2*(-1) + beta3*(-1) ## [1] 1.9903#kürzer
beta0+beta1-beta2-beta3## [1] 1.9903#Mittelwert direkt berechnen zum Vergleich:
mean(subset(data,music=="Indie Kid")$day1)## [1] 1.990294#speichern
M_Indie <- mean(subset(data,music=="Indie Kid")$day1)Metaller:
#kurz checken: Werte der Kontrastvariablen bei "Metaller"
head(subset(data,music=="Metaller"))## # A tibble: 6 x 9
## ticknumb music day1 day2 day3 change music1 music2 music3
## <dbl> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2111 Metaller 2.65 1.35 1.61 -1.04 0 2 -1
## 2 2514 Metaller 1.41 NA NA NA 0 2 -1
## 3 2533 Metaller 1.91 2.05 NA NA 0 2 -1
## 4 2535 Metaller 2.32 NA NA NA 0 2 -1
## 5 2586 Metaller 1.06 NA NA NA 0 2 -1
## 6 2662 Metaller 0.11 NA NA NA 0 2 -1#Einsetzen in die Modellgleichung ergibt den Mittelwert im Kriterium (day1) für die Indie Kid-Gruppe
beta0 + beta1*(0) + beta2*(2) + beta3*(-1)## [1] 1.83203#kürzer
beta0+2*beta2-beta3## [1] 1.83203#Mittelwert direkt berechnen zum Vergleich:
mean(subset(data,music=="Metaller")$day1)## [1] 1.832034#speichern
M_Metal <- mean(subset(data,music=="Metaller")$day1)No Musical Affiliation:
#kurz checken: Werte der Kontrastvariablen bei "No Musical Affiliation"
head(subset(data,music=="No Musical Affiliation"))## # A tibble: 6 x 9
## ticknumb music day1 day2 day3 change music1 music2 music3
## <dbl> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2338 No Musical Affiliation 0.84 NA NA NA 0 0 3
## 2 2401 No Musical Affiliation 0.88 0.08 NA NA 0 0 3
## 3 2504 No Musical Affiliation 1.11 0.44 0.55 -0.56 0 0 3
## 4 2509 No Musical Affiliation 2.17 NA NA NA 0 0 3
## 5 2515 No Musical Affiliation 1.76 1.64 1.58 -0.18 0 0 3
## 6 2520 No Musical Affiliation 1.38 0.02 NA NA 0 0 3#Einsetzen in die Modellgleichung ergibt den Mittelwert im Kriterium (day1) für die Indie Kid-Gruppe
beta0 + beta1*(0) + beta2*(0) + beta3*(3)## [1] 1.61283#kürzer
beta0+3*beta3## [1] 1.61283#Mittelwert direkt berechnen zum Vergleich:
mean(subset(data,music=="No Musical Affiliation")$day1)## [1] 1.612849#speichern
M_NA <- mean(subset(data,music=="No Musical Affiliation")$day1)Jetzt wissen wir also, wie unsere Gruppenmittelwerte über die Regressionskoeffizienten ausgedrückt werden können. Wir wissen auch, mit welchen Kontrastkoeffizienten wir wiederum die Mittelwerte gewichten müssen, um die Werte unserer drei Kontraste zu erhalten. Zur Erinnerung:
contrasts(data$music)## [,1] [,2] [,3]
## Crusty -1 -1 -1
## Indie Kid 1 -1 -1
## Metaller 0 2 -1
## No Musical Affiliation 0 0 3Wenn wir jetzt die über die Regressionskoeffizienten ausgedrückten Gruppenmittelwerte mit den Kontrastkoeffizienten gewichten und die Gleichungen auflösen, können wir sehen, wie die Regressionskoeffizienten mit den Mittelwertsunterschieden (den Werten der Kontraste) in Verbindung stehen:
Für Kontrast 1 (Indie-Crusty):
#Wert des Kontrasts über Mittelwerte
1*(M_Indie)-1*(M_Crusty)## [1] 0.1149333#Mittelwerte ausgedrückt über betas (siehe oben)
1*(beta0+beta1-beta2-beta3)-1*(beta0-beta1-beta2-beta3)## [1] 0.11494#Gleichung auflösen
2*beta1## [1] 0.11494#music1: 1*Indie Kid - 1*Crusty
1.87536-1.9903## [1] -0.11494#music2: 2*Metaller - Indie Kid - Crusty
2*(1.83203)-1.9903-1.87536## [1] -0.2016#music3: 3*No Musical Affiliation - Crusty - Indie Kid - Metaller
3*1.61283-1.87536-1.9903-1.83203## [1] -0.8592beta1 ist also der halbe Wert des Kontrasts 1.
Für Kontrast 2 (2*Metaller - Indie Kid - Crusty):
#Wert des Kontrasts über Mittelwerte
2*(M_Metal)-1*(M_Indie)-1*(M_Crusty)## [1] -0.2015871#Mittelwerte ausgedrückt über betas (siehe oben)
2*(beta0+2*beta2-beta3)-1*(beta0+beta1-beta2-beta3)-1*(beta0-beta1-beta2-beta3)## [1] -0.2016#Gleichung auflösen
6*beta2## [1] -0.2016beta2 ist also 1/6 des Werts des Kontrasts 2.
Für Kontrast 3 (3*No Musical Affiliation - Crusty - Indie Kid - Metaller):
#Wert des Kontrasts über Mittelwerte
3*(M_NA)-1*(M_Crusty)-1*(M_Indie)-1*(M_Metal)## [1] -0.8591428#Mittelwerte ausgedrückt über betas (siehe oben)
3*(beta0+3*beta3)-1*(beta0-beta1-beta2-beta3)-1*(beta0+beta1-beta2-beta3)-1*(beta0+2*beta2-beta3)## [1] -0.8592#Gleichung auflösen
12*beta3## [1] -0.8592beta3 ist also 1/12 des Werts des Kontrasts 2.
Die Regressionkoeffizienten geben also nicht direkt die durch die Kontraste kodierten Mittelwertsunterschiede an, aber für alle gilt: je größer das beta, desto größer auch der entsprechende Mittelwertsunterschied!