Pakete
Einführung
Beim allgemeinen linearen Regressionsmodell (ALM) handelt es sich um eine Erweiterung der multiplen linearen Regression, das ebenfalls eine kontinuierliche abhängige Variable vorhersagt, jedoch nicht auf dieses Skalenniveau hinsichtlich der unabhängigen Variablen beschränkt ist. D.h. es können ebenso nominale unabhängige Variablen (dichotom, polychotom) genutzt werden, sowohl exklusiv, als auch in Kombination. Unterteilt nach Skalenniveau und Anzahl der Prädiktoren, handelt es sich um das allgemeine Pendant
- zum t-Test für unabhängige Stichproben (eine dichotome UV)
- zur einfaktoriellen ANOVA (eine polychotome UV)
- zur einfachen linearen Regression (eine metrische UV)
- zur multiplen linearen Regression (mehr als eine metrische UV und/oder codierte nominale UV)
Bortz und Schuster bezeichnen das ALM daher als integrierenden Lösungsansatz (2010, S. 363). Diese Integration wird über die Codierung (s.u.) der UV bewerkstelligt.
Terminologie
Modelle
Während allgemeine lineare Regressionsmodelle (und somit auch deren Spezialfälle, s.o.) von einer normalverteilten abhängigen Variablen ausgehen (bzw. strenggenommen einer Normalverteilung der Störgrößen) und auf dieser Annahme beschränkt sind, lassen sich auch ALMs weiter veralgemeinern. Man spricht dann vom Verallgemeinerten linearen Modell (VLM) oder auch generalisiertes lineares Modell (GLM). Neben der Normalverteilung sind in dieser Klasse weitere Verteilungen wie die Binomial-, Poisson-, oder Gammaverteilung vertreten.
Funktionen
In R kann ein GLM mit der gleichnamigen Funktion glm() gefittet werden (siehe auch ?glm). Da der Paramter family in glm() per Default auf family = gaussian gesetzt ist, kann theoretisch auch stets diese Funktion verwendet werden (ohne die Normalverteilung explizit angeben zu müssen). Wir verwenden nachfolgend jedoch die lm()-Funktion, um die Abgrenzung von ALM und VLM zu verdeutlichen. Die lm()-Funktion ist bereits aus den vorherigen Beiträgen bekannt, anders als in anderen Statistikprogrammen wird in R also nur teilweise zwischen einfachen linaren-, multiplen linearen1-, allgemeinen linearen-2 und verallgemeinerten linearen Regressionsmodellen3 unterschieden.
Unterschiede im Englischen
Spoiler: Dieser Absatz führt zu maximaler Verwirrung, er dient lediglich als Erläuterung, wie welche Modelle und Codierungsvarianten im Englischen anders als bei Schuster benannt werden.
Im Englischen werden allgemeine lineare Modelle auch als Generalized Linear Model (GLM) bezeichnet, was verwirrend sein kann. Entsprechend verwerden wir den Ausdruck ALM.
Zudem werden wir nachfolgend verschiedene Codierungsvarianten betrachten: Dummy-, Effekt- und Kontrastcodierung. Im Englischen werden Codierungen ganz allgmein als Contrasts bezeichnet (was nicht zwangsläufig Kontraste meint, sondern einfach nur Kodierung bzw. Dummycodierung). Kontraste und andere verschiedene Formen der Codierung fallen unter Effects Coding, was nicht zwangsläufig Effektcodierung meint, sondern allgemein den Fall, wenn spezifische Annahmen über die Beziehung zwischen Variablen in die Codierung miteinfließen. Zum Effects Coding zählen Sum Coding (hier als Effektcodierung bezeichnet) sowie User-Defined Contrasts (hier Kontrastcodierung) und weitere Varianten wie z. B. Helmert Coding, die hier aber nicht behandelt werden.
Merke:
- wollen wir die allgmeine Codierung in Erfahrung bringen / bestimmen, nutzten wir die Funktion
contrasts() - wollen wir eine Dummycodierung erzeugen, nutzen wir
contr.treatment(), wobei wir die Referenzgruppe als Basis (base) setzen - wollen wir eine Effektcodierung (engl. Sum Coding) erzeugen, nutzen wir
contr.sum() - und wollen wir eine Kontrastcodierung (engl. User-Defined Contrasts) müssen wir unsere eigene Matrix festlegen (den per Default sind die Kontraste dummycodiert mit der alphabetisch ersten Gruppe als Referenz)
Codierungsvarianten
Es existieren verschiedene Codierungsvarianten, die mit einer unterschiedlichen inhaltlichen Interpretation der resultierenden Regressionskoeffizienten einhergehen, rechnerisch aber zum selben Deteriminationskoeffizienten führen. Es werden der Einfachheit halber ausschließlich balacierte Versuchspläne (gleiche Anzahl von n Merkmalsträgern pro Gruppe) besprochen:
Soll ein p-fach gestufter d.h. nominaler Prädiktor (UV) ein metrisches Kriterium (AV) vorhersagen, so werden \(p-1\) Indikatorveriablen benötigt.
Beispiel: Es soll die Einstellung zur Windkraft (Likert-Skala von 1 = max. Ablehnung bis 10 = max. Zustimmung) durch die Parteipräferenz der Probanden vorhergesagt werden. Sind \(p = 4\) Parteien vorhanden, dann werden \(p-1=3\) Indikatorvariablen benötigt. Es sind folgende Daten von \(N=16\) Personen gegeben, Gesamtmittel \(\bar{G}=4.81\).
Für sämltliche Codierungen muss die Form des Dataframes vom Wide- zum Long Format geändert werden:
df_party %>% pivot_longer(cols = everything(),
names_to = "var",
values_to = "y") %>%
mutate(var = factor(var)) %>%
arrange(var) -> df_long_partyDummycodierung
Je Stufe wird eine eigene Indikatorvariable erstellt, wobei es genau p-1 Indikatorvariablen gibt (für die letzte Stufe bzw. Gruppe ist keine Variable nötig).
- Personen der ersten Stufe erhalten in der ersten Indikatorvariable eine \(1\), für alle weiteren Variablen eine \(0\), alle übrigen Personen werden in der ersten Variable mit \(0\) codiert.
- Personen der zweiten Stufe erhalten für die zweite Indikatorvar. eine \(1\), für alle weiteren eine \(0\), etc.
df_long_party %>% mutate(# x0 = 1,
x1 = if_else(var == "a1", 1, 0),
x2 = if_else(var == "a2", 1, 0),
x3 = if_else(var == "a3", 1, 0)) %>%
select(-var) -> df_dummy_party
paged_table(df_dummy_party)Die Spalte x0 aus Einsen wird zur Berechnung des Y-Achsenabschnitts benötigt. R berechnet diesen jedoch ohnehin, wieso diese Variable hier auskommentiert ist und nachfolgend komplett weggelassen wird.
Interpretation
mdl_dummy_code <- lm(y ~ ., data = df_dummy_party)
# summary(mdl_dummy_code) %$% r.squared
coef(mdl_dummy_code)(Intercept) x1 x2 x3
4.75 2.00 -2.75 1.00 Die Y-Achsenabschnitt (Intercept) entspricht der durchschnittlichen Merkmalsausprägung in der p-ten Gruppe (letzte Gruppe ohne eigenes Gewicht, da durchgängig mit Null codiert). Diese Gruppe kann als Referenzgruppe bezeichnet werden, da die übrigen Gewichte zu ihrem Mittelwert (\(\bar{a}_4=4.75\)) in Relation stehen: eine jede Steigung errechnet sich als Differenz der Mittelwerte \(\bar{a}_i\) zur Referenzgruppe (da \(\hat{y}_i = b_0 + b_i\) für \(i=1, 2,...,p-1\)). Für \(\bar{a}_1=6.75\) ergibt sich entsprechend \(b_1 = 6.75-4.75=2.00\), für \(\bar{a}_2=2.00\) ein \(b_2 = 2.00-4.75=-2.75\) und für \(\bar{a}_3=5.75\) analog \(b_1 = 5.75-4.75=1.00\).
Die vorhergesagten Werte sind die Gruppenmittelwerte:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
6.75 6.75 6.75 6.75 2.00 2.00 2.00 2.00 5.75 5.75 5.75 5.75 4.75 4.75 4.75 4.75 Das Modell klärt mit \(R^2=.73\) ca. \(73\%\) Varianz auf.
Deterimnationskoeffizient und vorhergesagte Werte von der Art der Kodierung unabhängig.
Vorgehen in R
Soll eine Variable Dummycodiert werden, ist kein derartig komplliziertes Verfahren nötig (R nimmt per Default die alphabetisch erste Gruppe als Referenzgruppe, vgl. lm(y ~ var, data = df_long_party)). Die Erstellung von Indikatorvariablen kann ansonsten (insbesondere bei größeren Datensätzen) manuell wie folgt geschehen:
df_long_party %>% rownames_to_column(var = "id") %>%
mutate(id = as.numeric(id),
dummy = 1) %>%
spread(key = var, value = dummy, fill = 0) %>%
column_to_rownames(var = "id") y a1 a2 a3 a4
1 8 1 0 0 0
2 6 1 0 0 0
3 6 1 0 0 0
4 7 1 0 0 0
5 4 0 1 0 0
6 2 0 1 0 0
7 1 0 1 0 0
8 1 0 1 0 0
9 7 0 0 1 0
10 6 0 0 1 0
11 6 0 0 1 0
12 4 0 0 1 0
13 3 0 0 0 1
14 5 0 0 0 1
15 5 0 0 0 1
16 6 0 0 0 1Deutlich einfacher ist aber folgendes Vorgehen: Mit der Funktion contrasts() kann die Designmatrix X, die mit einem Faktor assoziert ist, angezeigt werden:
Stimmt die Refferenzvarialbe nicht mit jener überein, die als Referenz genutzt werden soll, so kann dies mit der Funktion contr.treatment() entsprechend angepasst werden:
contrasts(df_long_party$var) <- contr.treatment(4, base = 4)Das weitere Vorgehen ist in beiden Fällen analog.
lm(y ~ ., data = df_long_party)
Call:
lm(formula = y ~ ., data = df_long_party)
Coefficients:
(Intercept) var1 var2 var3
4.75 2.00 -2.75 1.00 Effektcodierung
Die Effektcodierung gleich der Dummycodierung überwiegend, es sind ebenso \(p-1\) Indikatorvariablen nötig, ihre Kodierung erfolgt analog mit ausnahme der letzten Gruppe (welche oben auf jeder Variable mit null kodiert wurde), diese wird stets mit minus eins codiert:
df_long_party %>% mutate(x1 = case_when(var == "a1" ~ 1,
var == "a4" ~ -1,
TRUE ~ 0),
x2 = case_when(var == "a2" ~ 1,
var == "a4" ~ -1,
TRUE ~ 0),
x3 = case_when(var == "a3" ~ 1,
var == "a4" ~ -1,
TRUE ~ 0)) %>%
select(-var) -> df_effect_party
paged_table(df_effect_party)Interpretation
mdl_effect_code <- lm(y ~ ., data = df_effect_party)
# summary(mdl_effect_code) %$% r.squared
coef(mdl_effect_code)(Intercept) x1 x2 x3
4.8125 1.9375 -2.8125 0.9375 Die Gewichte stehen hier nicht in Referenz zu einer Gruppe sondern zum Gesamtmittelwert \(\bar{G}\) aller Gruppen, hier \(\bar{G}=4.81\). Für \(\bar{a}_1=6.75\) ergibt sich \(b_1 = 6.75-4.81=1.94\), für \(\bar{a}_2=2.00\) ein \(b_2 = 2.00-4.81=-2.81\) und für \(\bar{a}_3=5.75\) analog \(b_1 = 5.75-4.81=0.94\). Der Intetercept ist hier also der Gesamtmittelwert! Vorhergesagte Werte sind dennoch die Mittelwerte.
Vorgehen in R
Das oben erläuterte Vorgehen kann analog fortgesetzt werden:
[,1] [,2] [,3]
a1 1 0 0
a2 0 1 0
a3 0 0 1
a4 -1 -1 -1Mit der gezeigten Designmatrix X ergeben sich dann dieselben Koeffizienten wie mit dem komplizierten manuellen Vorgehen:
lm(y ~ ., data = df_long_party)
Call:
lm(formula = y ~ ., data = df_long_party)
Coefficients:
(Intercept) var1 var2 var3
4.8125 1.9375 -2.8125 0.9375 Kontrastcodierung
Wenn – und nur dann wenn – ein balanciertes Design vorliegt und orthogonale Kontraste gerechnet werden sollen, kann anstelle der Kontraste (Post-Hoc) auch direkt eine entsprechende Codierung gewählt werden (im Fall geplanter Kontraste sind diese ohnehin vor der Berechnung bekannt).
Es gelten dabei die bereits bekannten Regeln für Kontraste (den Gewichten der Treatment-Mittelwerte):
- die Summe der Kontrastgewichte \(c_i\) eines jeden Kontrasts \(D\) muss null betragen: \(\Sigma c_i\stackrel{!}{=}0\)
- das Skalarprodukt jedes möglichen Tupels zweier Gewichte \(D_j\), \(D_{j'}\) muss null sein: \(D_{j}\circ D_{j'}\stackrel{!}{=}0\) (Orthogonalität)
- es sind \(p-1\) orthogonale Kontrase möglich, welche vollständig bestimmt werden müssen
Anstatt wie davor alle Gruppen miteinander zu vergleichen, betrachten wir nun folgende gewichteten Mittelwertsvergleiche (Notation \(D_j(c_1, c_2, c_3, c_4)\):
- \(D_1(\;1\;,\!-1\!\,,\;0\;,\;0\;)\): Vergleich der ersten beiden Gruppen
- \(D_2(\;0\;, \;0\;, \;1\;, \!-1\!\,)\): Vergleich der letzten beiden Gruppen
- \(D_3(\;1\;, \;1\;, \!-1\!\,, \!-1\!\,)\): Vergleich der ersten beiden mit den letzten beiden Gruppen
Für den letzten Kontrast ist es intuitiver \(D_3(0.5,0.5-0.5,-0.5)\) zu wählen
df_long_party %>% mutate(x1 = c(rep(c(1, -1), each=4), rep(0, 8)),
x2 = c(rep(0, 8), rep(c(1, -1), each=4)),
x3 = c(rep(0.5, 8), rep(-0.5, 8))) %>%
select(-var) -> df_contrast_party
paged_table(df_contrast_party) %>% str()Classes 'paged_df' and 'data.frame': 16 obs. of 4 variables:
$ y : num 8 6 6 7 4 2 1 1 7 6 ...
$ x1: num 1 1 1 1 -1 -1 -1 -1 0 0 ...
$ x2: num 0 0 0 0 0 0 0 0 1 1 ...
$ x3: num 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 -0.5 -0.5 ...Interpretation
mdl_contrast_code <- lm(y ~ ., data = df_contrast_party)
# summary(mdl_contrast_code) %$% r.squared
coef(mdl_contrast_code)(Intercept) x1 x2 x3
4.8125 2.3750 0.5000 -0.8750 Der Y-Achsenabschnitt entspricht (wie auch bei der Effektcodierung) dem Gesamtmittelwert: \(b_0=\bar{G}=4.81\)
Kontrast 1: Vergleich der ersten beiden Gruppen, d.h. \(b_1\) ist die halbierte Mittelwertsdifferenz beider Gruppen)
Kontrast 2: Vergleich der letzten beiden Gruppen, d.h. \(b_2\) ist die halbierte Mittelwertsdifferenz beider Gruppen
Kontrast 3: Vergleich der ersten beiden gegen die letzten beiden Gruppen, d.h. \(b_3\) ist die Differenz des Mittelwerts aus der Summe ersten beiden Gruppen gegen jenen Mittelwert der letzten beiden anderen zusammengenommenen Gruppen.
Vorgehen in R
Erstellen der Designmatrix X:
[,1] [,2] [,3]
[1,] 1 0 1
[2,] -1 0 1
[3,] 0 1 -1
[4,] 0 -1 -1Zuweisung als Kontraste und Modellberechnung:
t-Test
Beim Vergleich zweier Gruppen müssen auch diese Codiert werden, z. B. mit \(1\) für Gruppe \(A\) und \(-1\) für Gruppe \(B\) im Faller einer Effektcodierung.
Wir nehmen an, wir würden die Parteien aus dem Beispiel von oben nicht kennen und nur wissen, ob die jeweiligen Personen sich eher links vs. rechts zuordnen würden. Durch den letzten Kontrast kennen wir die Werte strenggenommen eigentlich schon, betrachten diese Daten nun aber für einen t-Test:
df_long_party %<>% mutate(party = case_when(var == "a1" | var == "a2" ~ "links",
var == "a3" | var == "a4" ~ "rechts"))
df_long_party %>% mutate(x = if_else(party == "links", 1, -1)) %>%
select(x, y) -> df_t_party
head(df_t_party)# A tibble: 6 × 2
x y
<dbl> <dbl>
1 1 8
2 1 6
3 1 6
4 1 7
5 1 4
6 1 2
Call:
lm(formula = y ~ x, data = df_t_party)
Residuals:
Min 1Q Median 3Q Max
-3.375 -1.500 0.250 1.625 3.625
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.8125 0.5402 8.908 3.83e-07 ***
x -0.4375 0.5402 -0.810 0.432
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.161 on 14 degrees of freedom
Multiple R-squared: 0.04475, Adjusted R-squared: -0.02348
F-statistic: 0.6558 on 1 and 14 DF, p-value: 0.4316Der Intercept ist auch hier wieder der Mittelwert beider Gruppen und \(b_1\) der die Abweichug des Gruppenmittelwerts aus Gruppe 1, d.h. in Relation dazu (\(b_0-\bar{x}_1=4.81-5.25=-0.44\)). Im Vergleich kommt der klassische t-Test zum selben p-Wert (außerdem: der t-Wert samt Freiheitsgraden des t-Tests für den Steigungskoeffizienten entspricht jenem des zwei-Stichproben-t-Tests)
t.test(y ~ x, data = df_t_party, var.equal = TRUE)
Two Sample t-test
data: y by x
t = 0.80983, df = 14, p-value = 0.4316
alternative hypothesis: true difference in means between group -1 and group 1 is not equal to 0
95 percent confidence interval:
-1.442373 3.192373
sample estimates:
mean in group -1 mean in group 1
5.250 4.375 Für eine Dummy-Codierung ergibt sich analog \(b_0=5.25\) als Mittelwert der Referenzgruppe (0 codiert) und \(b_1=4.37 - 5.25=-0.87\) als Mittelwert der Gruppe 1 (\(\bar{x}_1=4.37\)) in Relation dazu. Diese Kodierung ist deutlich sinnvoller, da aus dem Regressionskoeffizient direkt der Mittelwertsunterschied abgelesen werden kann.
df_long_party %>% mutate(x = if_else(party == "links", 1, 0)) %>%
select(x, y) -> df_t_party
lm(y ~ x, data = df_t_party) %>% summary()
Call:
lm(formula = y ~ x, data = df_t_party)
Residuals:
Min 1Q Median 3Q Max
-3.375 -1.500 0.250 1.625 3.625
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.250 0.764 6.872 7.66e-06 ***
x -0.875 1.081 -0.810 0.432
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.161 on 14 degrees of freedom
Multiple R-squared: 0.04475, Adjusted R-squared: -0.02348
F-statistic: 0.6558 on 1 and 14 DF, p-value: 0.4316Die Übereinstimmung der Werte ist gleich wie oben:
t.test(y ~ x, data = df_t_party, var.equal = TRUE)
Two Sample t-test
data: y by x
t = 0.80983, df = 14, p-value = 0.4316
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
-1.442373 3.192373
sample estimates:
mean in group 0 mean in group 1
5.250 4.375 Eine Kontrastkodierung wird hier nicht besprochen, ist aber möglich wenn inhaltlich belegbar wieso eine Gruppe ein bestimmtes Gewicht erhalten soll.
Der t-Test ist wie oben erläutert nur ein Spezialfall mit nur zwei Gruppen.
ANOVA
Entspricht der oben durchgeführten Regressionsberechnung. Die Quadratsummen der einzelnen Indikatoren ergeben in Summe \(R^2\)
df_effect_party %>% # mutate(x4 = c(rep(0, 12), rep(-1, 4))) %>%
aov(y ~ x1 + x2 + x3, data = .) %>%
summary.aov(intercept = TRUE) Df Sum Sq Mean Sq F value Pr(>F)
(Intercept) 1 370.6 370.6 243.658 2.46e-09 ***
x1 1 8.0 8.0 5.260 0.040671 *
x2 1 37.5 37.5 24.658 0.000328 ***
x3 1 4.7 4.7 3.082 0.104622
Residuals 12 18.3 1.5
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Modellvergleich
\[ F= \frac{R^2-R_b^2}{1-R^2}\times\frac{N-k}{k-k_b} \] Mit Index \(b\) für das beschränktere Modell.
Übungsaufgaben
Aufgabe 1
Teilaufgabe a
df_party %>% pivot_longer(cols = everything(),
names_to = "A",
values_to = "y") %>%
mutate(A = factor(A)) %>%
arrange(A) -> df_long_party
aov(y ~ A, data = df_long_party) %>% summary.aov(intercept = TRUE) Df Sum Sq Mean Sq F value Pr(>F)
(Intercept) 1 370.6 370.6 243.7 2.46e-09 ***
A 3 50.2 16.7 11.0 0.000926 ***
Residuals 12 18.3 1.5
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- Menu: Analysieren -> allgemeines lineares Modell -> univariat
- Zuordnung:
yzu Abhängige Variable undAzu Feste Faktoren
UNIANOVA y BY A
/METHOD=SSTYPE(3)
/INTERCEPT=INCLUDE
/CRITERIA=ALPHA(0.05)
/DESIGN=A.- Munu: Analysieren -> Mittelwerte und Proportionen vergleichen -> einfaktorielle Varianzanalyse findet sich eine Alternative mit gleichem Output
- Zuordnung der AV und UV:
ONEWAY y BY A
/ES=OVERALL
/MISSING ANALYSIS
/CRITERIA=CILEVEL(0.95).Teilaufgabe b
contrasts(df_long_party$A) <- contr.treatment(4, base = 4)
lm(y ~ A, data = df_long_party) %>% summary()
Call:
lm(formula = y ~ A, data = df_long_party)
Residuals:
Min 1Q Median 3Q Max
-1.7500 -0.8125 0.2500 0.5000 2.0000
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.7500 0.6166 7.703 5.52e-06 ***
A1 2.0000 0.8720 2.294 0.04067 *
A2 -2.7500 0.8720 -3.154 0.00832 **
A3 1.0000 0.8720 1.147 0.27383
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.233 on 12 degrees of freedom
Multiple R-squared: 0.7333, Adjusted R-squared: 0.6667
F-statistic: 11 on 3 and 12 DF, p-value: 0.000926- \(b_0=4.75\) ist der Mittelwert \(\bar{a}_4=4.75\) der Referenzgruppe \(a_4\)
- \(b_1=2.00\) ist das Gruppenmittel (\(\bar{a}_1=6.75\)) in Relation zur Referenzgruppe: \(6.75-4.75=2.00\)
- \(b_2=-2.75\) ist das Gruppenmittel (\(\bar{a}_2=2.00\)) in Rel. zur Referenzgr.: \(2.00-4.75=-2.75\)
- \(b_3=1.00\) ist das Gruppenmittel (\(\bar{a}_1=5.75\)) in Rel. zur Referenzgr.: \(5.75-4.75=1.00\)
Teilaufgabe c
Call:
lm(formula = y ~ A, data = df_long_party)
Residuals:
Min 1Q Median 3Q Max
-1.7500 -0.8125 0.2500 0.5000 2.0000
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.8125 0.3083 15.610 2.46e-09 ***
A1 1.9375 0.5340 3.628 0.003462 **
A2 -2.8125 0.5340 -5.267 0.000199 ***
A3 0.9375 0.5340 1.756 0.104622
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.233 on 12 degrees of freedom
Multiple R-squared: 0.7333, Adjusted R-squared: 0.6667
F-statistic: 11 on 3 and 12 DF, p-value: 0.000926Die Koeffizienten sind wie folgt zu interpretieren:
- \(b_1 = \hat{y}_1 - \hat{G} = 6.75-4.81 = 1.94\)
- \(b_2 = \hat{y}_3 - \hat{G} = 2.00-4.81 = -2.81\)
- \(b_3 = \hat{y}_3 - \hat{G} = 5.75-4.81 = 0.94\)
Teilaufgabe d
designmat <- matrix(c( 1, -1, 0, 0,
0, 0, 1, -1,
.5, .5, -.5, -.5),
ncol = 3)
contrasts(df_long_party$A) <- designmat
lm(y ~ A, data = df_long_party) %>% summary()
Call:
lm(formula = y ~ A, data = df_long_party)
Residuals:
Min 1Q Median 3Q Max
-1.7500 -0.8125 0.2500 0.5000 2.0000
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.8125 0.3083 15.610 2.46e-09 ***
A1 2.3750 0.4360 5.447 0.000148 ***
A2 0.5000 0.4360 1.147 0.273831
A3 -0.8750 0.6166 -1.419 0.181334
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.233 on 12 degrees of freedom
Multiple R-squared: 0.7333, Adjusted R-squared: 0.6667
F-statistic: 11 on 3 and 12 DF, p-value: 0.000926Teilaufgabe e
designmat<- matrix(c(.5, -.5, 0, 0,
0, 0, 1, -1,
.5, .5, -.5, -.5),
ncol = 3)
contrasts(df_long_party$A) <- designmat
lm(y ~ A, data = df_long_party) %>% summary()
Call:
lm(formula = y ~ A, data = df_long_party)
Residuals:
Min 1Q Median 3Q Max
-1.7500 -0.8125 0.2500 0.5000 2.0000
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.8125 0.3083 15.610 2.46e-09 ***
A1 4.7500 0.8720 5.447 0.000148 ***
A2 0.5000 0.4360 1.147 0.273831
A3 -0.8750 0.6166 -1.419 0.181334
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.233 on 12 degrees of freedom
Multiple R-squared: 0.7333, Adjusted R-squared: 0.6667
F-statistic: 11 on 3 and 12 DF, p-value: 0.000926Prüfe:
Aufgabe 2
Teilaufgabe a
contr_mat_new <- matrix(c( .75, .5, .25,
-.25, .5, .25,
-.25, -.5, .25,
-.25, -.5, -.75),
byrow = TRUE,
ncol = 3)
df_long_party %>% mutate(x1 = rep(contr_mat_new[,1], each = 4),
x2 = rep(contr_mat_new[,2], each = 4),
x3 = rep(contr_mat_new[,3], each = 4)) %>%
lm(y ~ x1 + x2 + x3, .) %>% summary()
Call:
lm(formula = y ~ x1 + x2 + x3, data = .)
Residuals:
Min 1Q Median 3Q Max
-1.7500 -0.8125 0.2500 0.5000 2.0000
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.8125 0.3083 15.610 2.46e-09 ***
x1 4.7500 0.8720 5.447 0.000148 ***
x2 -3.7500 0.8720 -4.300 0.001031 **
x3 1.0000 0.8720 1.147 0.273831
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.233 on 12 degrees of freedom
Multiple R-squared: 0.7333, Adjusted R-squared: 0.6667
F-statistic: 11 on 3 and 12 DF, p-value: 0.000926Teilaufgabe b
Es lässt sich ein Gleichungssystem aus den Mittelwerten B und der Matrix erstellen.
Wie in der Aufgabenstellung angemerkt muss diese aber zunächst um eine Spalte erweitert werden, die in R in der Modellrechnung nicht verlangt wird und folglich in der Matrx nicht vorhanden ist:
[,1] [,2] [,3] [,4] [,5]
[1,] 1 0 0 0 4.7500
[2,] 0 1 0 0 -3.7500
[3,] 0 0 1 0 1.0000
[4,] 0 0 0 1 4.8125Den Lösungsweg kann man ebenfalls anzeigen lassen und damit nachvollziehen. M.E. nach ist das sinnvoller, als die Gleichungen durch “Hinsehen” zu lösen, denn schafft man dies nicht und spickt in die Lösung, ist die unten aufgeführte Interpretation klar bzw. umgekehrt erweckt dieses Vorgehen den Eindruck, man könne das Gleichungssystem v.A. dadurch “einfach durch Hinsehen” lösen, weil man die Interpretation bereits kennt. Das mathematische Lösen des Gleichungssystems sieht wie folgt aus:
matlib::echelon(A = contr_mat_comp, B = B, verbose = TRUE, fractions = TRUE)
Initial matrix:
[,1] [,2] [,3] [,4] [,5]
[1,] 3/4 1/2 1/4 1 27/4
[2,] -1/4 1/2 1/4 1 2
[3,] -1/4 -1/2 1/4 1 23/4
[4,] -1/4 -1/2 -3/4 1 19/4
row: 1
multiply row 1 by 4/3
[,1] [,2] [,3] [,4] [,5]
[1,] 1 2/3 1/3 4/3 9
[2,] -1/4 1/2 1/4 1 2
[3,] -1/4 -1/2 1/4 1 23/4
[4,] -1/4 -1/2 -3/4 1 19/4
multiply row 1 by 1/4 and add to row 2
[,1] [,2] [,3] [,4] [,5]
[1,] 1 2/3 1/3 4/3 9
[2,] 0 2/3 1/3 4/3 17/4
[3,] -1/4 -1/2 1/4 1 23/4
[4,] -1/4 -1/2 -3/4 1 19/4
multiply row 1 by 1/4 and add to row 3
[,1] [,2] [,3] [,4] [,5]
[1,] 1 2/3 1/3 4/3 9
[2,] 0 2/3 1/3 4/3 17/4
[3,] 0 -1/3 1/3 4/3 8
[4,] -1/4 -1/2 -3/4 1 19/4
multiply row 1 by 1/4 and add to row 4
[,1] [,2] [,3] [,4] [,5]
[1,] 1 2/3 1/3 4/3 9
[2,] 0 2/3 1/3 4/3 17/4
[3,] 0 -1/3 1/3 4/3 8
[4,] 0 -1/3 -2/3 4/3 7
row: 2
multiply row 2 by 3/2
[,1] [,2] [,3] [,4] [,5]
[1,] 1 2/3 1/3 4/3 9
[2,] 0 1 1/2 2 51/8
[3,] 0 -1/3 1/3 4/3 8
[4,] 0 -1/3 -2/3 4/3 7
multiply row 2 by 2/3 and subtract from row 1
[,1] [,2] [,3] [,4] [,5]
[1,] 1 0 0 0 19/4
[2,] 0 1 1/2 2 51/8
[3,] 0 -1/3 1/3 4/3 8
[4,] 0 -1/3 -2/3 4/3 7
multiply row 2 by 1/3 and add to row 3
[,1] [,2] [,3] [,4] [,5]
[1,] 1 0 0 0 19/4
[2,] 0 1 1/2 2 51/8
[3,] 0 0 1/2 2 81/8
[4,] 0 -1/3 -2/3 4/3 7
multiply row 2 by 1/3 and add to row 4
[,1] [,2] [,3] [,4] [,5]
[1,] 1 0 0 0 19/4
[2,] 0 1 1/2 2 51/8
[3,] 0 0 1/2 2 81/8
[4,] 0 0 -1/2 2 73/8
row: 3
multiply row 3 by 2
[,1] [,2] [,3] [,4] [,5]
[1,] 1 0 0 0 19/4
[2,] 0 1 1/2 2 51/8
[3,] 0 0 1 4 81/4
[4,] 0 0 -1/2 2 73/8
multiply row 3 by 1/2 and subtract from row 2
[,1] [,2] [,3] [,4] [,5]
[1,] 1 0 0 0 19/4
[2,] 0 1 0 0 -15/4
[3,] 0 0 1 4 81/4
[4,] 0 0 -1/2 2 73/8
multiply row 3 by 1/2 and add to row 4
[,1] [,2] [,3] [,4] [,5]
[1,] 1 0 0 0 19/4
[2,] 0 1 0 0 -15/4
[3,] 0 0 1 4 81/4
[4,] 0 0 0 4 77/4
row: 4
multiply row 4 by 1/4
[,1] [,2] [,3] [,4] [,5]
[1,] 1 0 0 0 19/4
[2,] 0 1 0 0 -15/4
[3,] 0 0 1 4 81/4
[4,] 0 0 0 1 77/16
multiply row 4 by 4 and subtract from row 3
[,1] [,2] [,3] [,4] [,5]
[1,] 1 0 0 0 19/4
[2,] 0 1 0 0 -15/4
[3,] 0 0 1 0 1
[4,] 0 0 0 1 77/16Die Brüche (oben Dezimalzahlen) in Spalte 5 ergeben genau die Koeffizienten. Die Spalten 1 bis 3 (4 hatten wir hinzugefügt) ergeben genau das Muster der Designmatrix, welche wir als Dummy-Codierung kennengelernt hatten. Inhaltlich ist die Interpretation also, dass der Intercept der Mittelwert von Referenzgruppe vier ist und die anderen dazu in Relation gesetzt werden (muss hier nicht erneut ausgeführt werden). Wichtig ist vielmehr, dass wir die Designmatrix rückwikrend herstellen konnten. Im vorliegenden Fall ist die Interpretation bei Parteien (Mittelwertsdifferenz Gruppe 1 zu 2 bzw. 2 zu 3 bzw. 3 zu 4) wenig sinnvoll. Für Zeitreihenwerte wäre aber genau so eine Cordierung sinnvoll (d.h. Mittelwertsdifferenz von Zeitpunkt 1 zu Zeitpunkt 2 bzw. ZP 2 zu ZP 3 bzw. ZP 3 zu ZP 4). Wir merken uns, dass für repeated measurements4 dies sinnvoll ist.
Aufgabe 3
Call:
lm(formula = y ~ A, data = df_long_party)
Residuals:
Min 1Q Median 3Q Max
-1.7500 -0.8125 0.2500 0.5000 2.0000
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.8125 0.3083 15.610 2.46e-09 ***
A1 4.7500 0.8720 5.447 0.000148 ***
A2 -3.7500 0.8720 -4.300 0.001031 **
A3 1.0000 0.8720 1.147 0.273831
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.233 on 12 degrees of freedom
Multiple R-squared: 0.7333, Adjusted R-squared: 0.6667
F-statistic: 11 on 3 and 12 DF, p-value: 0.000926Button Kontrase, Auswahl in Dropdown: Wiederholt
UNIANOVA y BY A
/CONTRAST(A)=Repeated
/METHOD=SSTYPE(3)
/INTERCEPT=INCLUDE
/CRITERIA=ALPHA(0.05)
/DESIGN=A.Aufgabe 4
Wie sehen die drei Indikatorvariablen (Designmatrix) aus, die verwendet wurden, um die repeated-Kontrastmatrix (in Aufgabenstellung gegeben) zu erstellen? Antwort
MANOVA y BY A(1,4)
/CONTRAST(A)=REPEATED
/PRINT=DESIGN(ONEWAY).