---
title: Lineares Modell, kategoriale Prädiktoren, Vergleich balanced und unbalanced Design
output:
  html_document:
    toc: true
    toc_depth: 4
    toc_float: yes
    theme: paper
header-includes:
 - \usepackage{comment}
params:
 soln: TRUE
 todo: FALSE
csl: ../../apa6.csl
bibliography: ../../mv.bib
---




### [Rmd](lm_cat_unbalanced.Rmd)

# Konzeptionelles

(Aus-)Balancierte varianzanalytische Designs sind Designs mit gleicher Zellbesetzung.
Die Effekte sind dann unabhängig voneinander.
Regressionsgewichte verändern sich nicht mit Aufnahme neuer Prädiktoren. 
Die Aufnahmereihenfolge der Prädiktoren (Faktoren) in das Modell spielt keine Rolle, 
auch nicht bei Quadratsummenbildung vom Typ I.

Unbalancierte varianzanalytische Designs haben ungleiche Zellbesetzungen.
Die Effekte sind abhängig voneinander.
Die Aufnahmereihenfolge der Prädiktoren (Faktoren) spielt eine wichtige Rolle.
Ein hinzugenommener Prädiktor nimmt die maximal ihm zur Verfügung stehende Varianz, 
auch die mit anderen potenziellen Prädiktoren gemeinsame.
Je nach Vorgehen (Aufnahmereihenfolge, verwendeter Quadratsummentyp) werden unterschiedliche Hypothesen geprüft.
SumSq der Residuen im 3-Prädiktor-Modell ist SumSq error-Quadratsumme im Anova-Modell

Regressionsgewichte verändern sich nicht mit Aufnahme neuer Prädiktoren
Faktoren im balancieren Modell sind unabhängig (orthogonal)
Aufnahmereihenfolge der Faktoren (incl. WW) spielt keine Rolle


Die relevanten Ausführungen dazu finden sich in 
[Field (2012, p. 474 - 477)][field] in **Kapitel 11: Analysis of Covariance, ANCOVA**

Etwas über den Hintergrund unter
[lm_cat_unbal_ss_explained.html](lm_cat_unbal_ss_explained.html)

<!-- [ss_explained](ss_explained) -->

[https://www.uni-kiel.de/psychologie/dwoll/r/ssTypes.php]
bzw. als [Rmd-File](./lm_cat_unbal_ss_explained.Rmd)

Konzeptioneller Hintergrund und ein empfehlenswertes Tutorial mit Beispiel von Ista Zahn unter
[InteractionsAndTypesOfSS.pdf][zahn] 
[external link](https://sites.psu.edu/stat461psbsp2013/wp-content/uploads/sites/1906/2013/03/InteractionsAndTypesOfSS.pdf)
[lokale Kopie Paper Ista Zahn](https://md.psych.bio.uni-goettingen.de/unit/lm_cat_unbalanced/InteractionsAndTypesOfSS.pdf "Zahn Unbalanced Tutorial")
mit dem zugehörigen Datensatz unter
[Zahn Data](https://md.psych.bio.uni-goettingen.de/data/div/zahn-unbal-paper.txt)





# Versteh-Beispiel Everitt


2x2-Design, je 2 Faktorstufen, Everitt (2010)

Unterschiede zwischen balanced und unbalanced Designs.
Im Beispieldatensatz sind die Faktoren A und B bereits effektkodiert bzw. kontrast-kodiert und können so direkt als Dummy-Variablen benutzt werden, 
um einen regressionsanalytischen Ansatz zu rechnen.

## balanced design

Differenz der Regressions SumSq zwischen 1-Prädiktor-Modell und 2-Prädiktor-Modell ist SumSq für Faktor B in Anova-Modell
Differenz der Regressions SumSq zwischen 2-Prädiktor-Modell und 3-Prädiktor-Modell (Interaktion) ist SumSq für Interaktionsterm in Anova-Modell

Regressionsgewichte verändern sich nicht mit Aufnahme neuer Prädiktoren
Faktoren im balancieren Modell sind unabhängig (orthogonal)
Aufnahmereihenfolge der Faktoren (incl. WW) spielt keine Rolle

die konstante Gesamt-Quadratsumme wird auf immer mehr Effekte aufgeteilt
Residual-Quadratsumme nimmt immer mehr ab.



```{r}
# read datafiles
dd   <- read.delim("https://md.psych.bio.uni-goettingen.de/mv/data/be/e-2x2-bal.txt")
# generate factors
dd$f.A <- factor(dd$A, labels=c('a1', 'a2'))
dd$f.B <- factor(dd$B, labels=c('b1', 'b2'))
#contrasts(dd$f.A) <- contr.sum(2)
#contrasts(dd$f.B) <- contr.sum(2)
#contrasts(dd$f.A) <- contr.treatment(2)
#contrasts(dd$f.B) <- contr.treatment(2)

# frequencies
table(dd$f.A, dd$f.B)

require(Rmisc)
Rmisc::summarySE(data = dd, "y", groupvars=c("f.A", "f.B"))

# explanatory variables are uncorrelated
cor(dd$A, dd$B)

# fit regression with explanatory variable A
m.A <- lm(y ~ f.A, data=dd)
cbind(dd, model.matrix(m.A))
summary(m.A)
# mean(dd$y[dd$f.A == "a1"])
# Rmisc::summarySE(data = dd, "y", groupvars=c("f.A"))
anova(m.A)
anova(m.A)$`Sum Sq`
sum(anova(m.A)$`Sum Sq`)
# add B
m.A.B <- lm(y ~ f.A + f.B, data=dd)
cbind(dd, model.matrix(m.A.B))
m.B.A <- lm(y ~ f.B + f.A, data=dd)
cbind(dd, model.matrix(m.B.A))
summary(m.A.B)
summary(m.B.A)
# in SS type 1 SS are equal, independent of sequence to enter 
anova(m.A.B)
anova(m.B.A)
# mean(dd$y)
# mean(dd$y[dd$f.A == "a1"])
# mean(dd$y[dd$f.A == "a1" & dd$f.B == "b1"])
# Rmisc::summarySE(data = dd, "y", groupvars=c("f.A", "f.B"))
# Rmisc::summarySE(data = dd, "y", groupvars=c("A", "B"))
# coefficients(m.A.B)
# predict(m.A.B)
# dd$y

# model.matrix(m.A.B)
# mean(dd$y[dd$f.B == "b1"])
# mean(c(dd$y[dd$f.A == "a1"], dd$y[dd$f.B == "b1"]))


# difference in SumSq between m.A and m.A.B is SumSq for factor B
anova(m.A, m.A.B)
anova(m.A)$`Sum Sq`
anova(m.A.B)$`Sum Sq`

anova(m.A, m.B.A)
anova(m.B.A)$`Sum Sq`

# add interaction
m.A.B.AB <- lm(y ~ f.A + f.B + f.A:f.B, data=dd)
cbind(dd, model.matrix(m.A.B.AB))
summary(m.A.B.AB)
Rmisc::summarySE(data = dd, "y", groupvars=c("f.A", "f.B"))
anova(m.A.B.AB)
anova(m.A.B, m.A.B.AB)
# difference SumSq of regression between 2 (A, B) and 3 (A, B, A:B) predictor model is SumSq for interaction
# SumSq der Residuen im 3-Prädiktor-Modell ist SumSq error-Quadratsumme im Anova-Modell
# SS resid in 3 predictor model is SS error in Anova (residuals)
# regression weights (coefficients) don't change when new explanatory variables are entered in model
# Regressionsgewichte verändern sich nicht mit Aufnahme neuer Prädiktoren
# factors and their effects are orthogonal
# sequence of entered explanatory variables doesn't matter
summary(m.A.B)
summary(m.B.A)
summary(m.A.B.AB)
anova(m.A.B)
anova(m.B.A)
anova(m.A.B.AB)

# a visualization
require(ggplot2)
ggplot(dd, aes(f.A, y, color=f.B)) +
  # scale_fill_manual(values=c_colors) +
  stat_summary(fun.y=mean, geom="bar", position=position_dodge()) +
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", position=position_dodge(width=0.7), width=0.2) + 
  # put points on top of bars
  stat_summary(fun.data = mean_cl_normal, geom = "point", shape=21, size=5, position=position_dodge(width=0.7)) +
	geom_hline(aes(yintercept=coefficients(m.A.B.AB)[1])) + 
  # add verbose axis labels
  labs(x = "group B", y = "y")



# the constant total SS is divided in more and more effects
# residual SS diminishes with every additional effect entered in the model

```




## unbalanced design (same experimental design)

derselbe Datensatz mit ungleicher Teilstichproben-Größe (teils ein paar Beobachtungen mehr)
Differenz der Regressions SumSq zwischen 1-Prädiktor-Modell und 2-Prädiktor-Modell ist SumSq für Faktor B in Anova-Modell

Wechselwirkung dazu
Reihenfolge der Aufnahme ins Modell verändert die Parameter der Prädiktoren im Modell
sowohl die Regressionsparameter als auch die Quadratsummen
und damit schliesslich auch die statistische Beurteilung.

Die QS der Interaktion kommt hinzu. 
QS der Haupteffekte sind verschieden und abhängig von der Aufnahmereihenfolge.
QS der Interaktion bleibt gleich, unabhängig von der Aufnahmereihenfolge (gilt für die höchste WW).

Berechnung mit `aov()`

```{r}
# unbalanced design
# same experimental design, some observations added, unequal cell size
dd <- read.delim("https://md.psych.bio.uni-goettingen.de/mv/data/be/e-2x2-unbal.txt")

# generate factors
dd$f.A <- factor(dd$A, labels=c('a1', 'a2'))
dd$f.B <- factor(dd$B, labels=c('b1', 'b2'))

# n in cells
table(dd$f.A, dd$f.B)

# descriptives
require(Rmisc)
Rmisc::summarySE(data = dd, "y", groupvars=c("f.A", "f.B"))

# explanatory variables are correlated
cor(dd$A, dd$B)

# fit regression including explanatory variable A
mu.A <- lm(y ~ f.A, data=dd)
unique(model.matrix(mu.A))
summary(mu.A)
anova(mu.A)

# add B to a slightly more complex model
# as we can see, input sequence of explanatory variables matters in unbalanced designs
mu.A.B <- lm(y ~ f.A + f.B, data=dd)
unique(model.matrix(mu.A.B))
mu.B.A <- lm(y ~ f.B + f.A, data=dd)
unique(model.matrix(mu.B.A))
# coefficients keep the same
summary(mu.A.B)
summary(mu.B.A)
# SS differ due to the sequence to enter
anova(mu.A.B)
anova(mu.B.A)
# total SS stays the same
sum(anova(m.A.B)$`Sum Sq`)
sum(anova(m.B.A)$`Sum Sq`)

# difference SS is SS of factor B

# SS of highest interaction stays the same, also in unbalaced designs
mu.A.B.AB <- aov(y ~ f.A + f.B + f.A:f.B, data=dd)
mu.B.A.AB <- aov(y ~ f.B + f.A + f.A:f.B, data=dd)

# sequence to enter alters SS
# A and B add up to the same SS
# interaction stays unchanged
anova(mu.A.B.AB)
anova(mu.B.A.AB)

#c_colors <- c("#66ff66", "#ff6666", "#6666ff")
# sequence to enter alters coefficients and SS
# and therefore the statistical tests

ggplot(dd, aes(f.A, y, color=f.B)) +
  # scale_fill_manual(values=c_colors) +
  stat_summary(fun.y=mean, geom="bar", position=position_dodge()) +
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", position=position_dodge(width=0.7), width=0.2) + 
  # put points on top of bars
  stat_summary(fun.data = mean_cl_normal, geom = "point", shape=21, size=5, position=position_dodge(width=0.7)) +
	geom_hline(aes(yintercept=coefficients(mu.A.B.AB)[1])) + 
  # add verbose axis labels
  labs(x = "group B", y = "y")
# let us control the visualized values
coefficients(mu.A.B.AB)
Rmisc::summarySE(data = dd, "y", groupvars=c("f.A", "f.B"))
aov(mu.A.B.AB)
sum(summary(aov(mu.A.B.AB))[[1]]$`Sum Sq`)

aov(mu.B.A.AB)
sum(summary(aov(mu.B.A.AB))[[1]]$`Sum Sq`)

```





# Beispiel STAI-State Angstgruppe und Geschlecht


In bestimmten Gruppen werden die Beobachtungen 'reduziert' - es wird eine Zufallsauswahl gezogen. 
Hierdurch werden die Zellhäufigkeiten in den Faktorstufen-Kombinationen ungleich.


```{r}
# get data
dd <- read.delim("https://md.psych.bio.uni-goettingen.de/mv/data/virt/v_stait_exam_wide.txt")[c('subj','gender','group', 'neutr')]
length(dd$subj)
# percentages
p.del <- c(.5, .7, .8)
s.del <- c()

# select some women from group 'av' to delete in dataframe
s.a <- dd$subj[dd$gender == 'f' & dd$group == 'av'] 
#s.d <- s.a[!is.na(sample(s.a, length(s.a) * p.del[1]))]
set.seed(42)  # to get reproducable results
s.d <- sample(s.a, length(s.a) * p.del[1])
s.del <- c(s.del, s.d)

# select some trained women to delete in dataframe
s.a <- dd$subj[dd$gender == 'f' & dd$group == 'hi'] 
set.seed(43)  # to get reproducable results
s.d <- sample(s.a, length(s.a) * p.del[2])
s.del <- c(s.del, s.d)

# select some anxious men to delete in dataframe
s.a <- dd$subj[dd$gender == 'm' & dd$group == 'lo'] 
set.seed(44)  # to get reproducable results
s.d <- sample(s.a, length(s.a) * p.del[3])
s.del <- c(s.del, s.d)


# now take all selected observations out of dataframe
dd.u <- dd[-s.del,]
# write.table(dd.u, file="v_stait_exam_wide_unbal.txt", quote=F, sep="\t", row.names=F)

# this way the unbalanced dataframe was created
# we read a reduced dataframe constructed in the above described way
dd.u <- read.delim("https://md.psych.bio.uni-goettingen.de/mv/data/virt/v_stait_exam_wide_unbal.txt")
# control crosstabs
table(dd.u$gender, dd.u$group)

# crosstabulation including test of indepencence
mytable <- xtabs(~gender + group, data=dd.u)
ftable(mytable) # print table
summary(mytable) # chi-square test of indepedence 

# write dummy variables to dataframe
dd.u <- cbind(dd.u, model.matrix(neutr ~ gender + group + gender:group, data=dd.u))
# now we can check correlation of dummy variables
cor(dd.u[c('genderm', 'grouphi', 'grouplo', 'genderm:grouphi', 'genderm:grouplo')])

# descriptives
require(Rmisc)
summarySE(data=dd.u, "neutr", c("group","gender"))
detach("package:Rmisc", unload=TRUE)

# package(car) is required. It defines 'Anova()'. The capital 'A' is significant!
require(car)
# the balanced and unbalanced designs compared
m.b <- (lm(neutr ~ gender*group, data=dd  ))
m.u <- (lm(neutr ~ gender*group, data=dd.u))

# model test
car::Anova(m.b)
car::Anova(m.u)

# a look at sums of squares
car::Anova(m.b)$`Sum Sq`; sum(car::Anova(m.b)$`Sum Sq`)
car::Anova(m.u)$`Sum Sq`; sum(car::Anova(m.u)$`Sum Sq`)

# parameters in detail
# balanced
summary(m.b)
# unbalanced
summary(m.u)
# ... no surprise, results differ due to the loss of data in dd.u (m.u), although it is part of dd (m.b)

# with SS Type I, there should be differences depending on sequence of entrance
summary(aov(lm(neutr ~ gender + group  + gender:group, data=dd.u)))
summary(aov(lm(neutr ~ group  + gender + gender:group, data=dd.u)))

# the highest order interaction is always the same
summary(aov(lm(neutr ~                group + gender + gender:group, data=dd.u)))
summary(aov(lm(neutr ~ gender:group + group + gender,                data=dd.u)))

# the three types of SS result in differences
summary(aov(lm(neutr ~ gender + group + gender:group, data=dd.u)))
Anova(lm(neutr ~ gender + group + gender:group, data=dd.u), type=2)
Anova(lm(neutr ~ gender + group + gender:group, data=dd.u), type=3)

# but coefficients keep the same
coefficients(lm(neutr ~ gender + group + gender:group, data=dd.u))
coefficients(lm(neutr ~ group + gender + gender:group, data=dd.u))
# the highes order interaction is always the same
coefficients(lm(neutr ~ gender:group + group + gender, data=dd.u))

require(Rmisc)
summarySE(data=dd.u, "neutr", c("group","gender"))
detach("package:Rmisc", unload=TRUE)


```



## Grafik dazu


```{r}
dd.u <- read.delim("https://md.psych.bio.uni-goettingen.de/mv/data/virt/v_stait_exam_wide_unbal.txt")
# we want SS type I parameters
m.m <- aov(lm(neutr ~ gender + group + gender:group, data=dd.u))
require(ggplot2)

ggplot(dd.u, aes(x=group, y=neutr, fill=gender)) +
  # define colors
  scale_fill_manual(values=c("#CCCCCC","#FFFFFF")) +
  # create white background etc
  #theme_bw() +
  stat_summary(fun.y=mean, geom="bar", position=position_dodge(width=0.7), colour='black') +
  # add error information
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", position=position_dodge(width=0.7), width=.2) +
  # put big points on top of bars
  stat_summary(fun.data = mean_cl_normal, geom = "point", shape=21, size=5, position=position_dodge(width=0.7)) +
		geom_hline(aes(yintercept=coefficients(m.m)[1])) + 
  # add verbose axis labels
  labs(x = "Anxiety-Group", y = "STAI-Trait")


```





## mit `ezAnova()` der `library(ez)`


Wenn wir wissen, was genau geprüft werden soll und wenn es nicht um das Verständnis der Vorgänge 'hinter den Kulissen' geht, 
ist `ezAnova()` ein komfortabler Weg, varianzanalytische Designs zu rechnen.
Wir können den Quadratsummen-Typ explizit festlegen.


```{r}
dd.u <- read.delim("https://md.psych.bio.uni-goettingen.de/mv/data/virt/v_stait_exam_wide_unbal.txt")
require(ez)
# fit model using ezANOVA
m.f <- ezANOVA(data=dd.u, dv=.(neutr), wid = .(subj), between=.(gender, group), detailed=TRUE, type=2)
m.f
#ezANOVA(data=dd.r, dv=.(anx), wid = .(subj), within = .(sit), detailed=TRUE, type=1)

# compare it to Anova
require(car)
Anova(lm(neutr ~ gender + group + gender:group, data=dd.u), type=2)

# ez offers ezPlot()
ezPlot(data=dd.u, dv=.(neutr), wid = .(subj), between=.(gender, group), x=group)
#ezPlot(data=dd.u, dv=.(neutr), wid = .(subj), between=.(gender, group), x=group, split=.(group))
ezPlot(data=dd.u, dv=.(neutr), wid = .(subj), between=.(gender, group), x=group, split=.(gender))

# descriptives - the ez way
ezStats(data=dd.u, dv=.(neutr), wid = .(subj), between=.(gender, group))

# some visualizations in ggplot
require(ggplot2)
ggplot(dd.u, aes(x=group, y=neutr)) +
  # define colors
  scale_fill_manual(values=c("#FFFFFF", "#CCCCCC", '#AAAAAA')) +
  # create white background etc
  #theme_bw() +
  geom_point() +
  # add verbose axis labels
  labs(x = "Anxiety-Group", y = "STAI-Trait")


ggplot(dd.u, aes(x=group, y=neutr, fill=gender)) +
  # define colors
  scale_fill_manual(values=c("#FFFFFF", "#CCCCCC", '#AAAAAA')) +
  # create white background etc
  #theme_bw() +
  stat_summary(fun.y=mean, geom="bar", position=position_dodge(width=0.7), colour='black') +
  # add error information
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", position=position_dodge(width=0.7), width=0.2) +
  # put big points on top of bars
  stat_summary(fun.data = mean_cl_normal, geom = "point", shape=21, size=5, position=position_dodge(width=0.7)) +
  # add verbose axis labels
  labs(x = "Anxiety-Group", y = "STAI-Trait")

```





# Beispiele Übungen / Exercises

## PACS Data

Führen Sie einfache VAs durch mit den PACS-Daten, beispielsweise, ob sich die verschiedenen Gewichts-Untergruppen hinsichtlich PACS-Werte unterscheiden.



## Zahn, Geschlecht und Bezahlung

Vollziehen Sie das Beispiel von Zahn (2010) nach

[InteractionsAndTypesOfSS.pdf](InteractionsAndTypesOfSS.pdf) 
mit dem zugehörigen 
[Datensatz](https://md.psych.bio.uni-goettingen.de/data/div/zahn-unbal-paper.txt)

[zahn](zahn.html)





## Abwägen von Vor- und Nachteilen

Überlegen Sie sich Situationen, in denen es bei unbalanced Data von Vorteil sein könnte, mit QS Type I zu rechnen d. h. in denen diese Art der Auswertung näher an einer psychologischen Fragestellung wäre, als die Anwendung von QS III


# Links und Referenzen

Konzept und ein Teil der Daten entstammen dem Buch [Everitt (2010)][everitt]

[Field (2012, p. 474 - 477)][field] in **Kapitel 11: Analysis of Covariance, ANCOVA**

[InteractionsAndTypesOfSS.pdf](InteractionsAndTypesOfSS.pdf) 
mit dem zugehörigen 
[Datensatz](https://md.psych.bio.uni-goettingen.de/data/div/zahn-unbal-paper.txt)

<!-- this does not appear in html-page. It's used only to arm links and references. -->
[field]: https://www.uk.sagepub.com/dsur/main.htm "Field et al (2012). Discovering Statistics using R"
[everitt]: https://www.crcpress.com/product/isbn/9781439807699 "Brian S. Everitt (2010). Multivariable Modeling and Multivariate Analysis for the Behavioral Sciences"
[zahn]: https://md.psych.bio.uni-goettingen.de/unit/lm_cat_unbalanced/InteractionsAndTypesOfSS.pdf "Zahn Unbalanced Tutorial"
[ss_explained]: https://md.psych.bio.uni-goettingen.de/unit/lm_cat_unbalanced/lm_cat_unbal_ss_explained.html "Unbalanced SS Tutorial"

[background](https://goanna.cs.rmit.edu.au/~fscholer/anova.php)

ein ganzes [Buch](https://www.ievbras.ru/ecostat/Kiril/R/Biblio/R_eng/Logan%20Biostatistical.pdf)
zu Anova-Designs, auch unbalanced, mit R-Code




# check

- estimated reaction variable is the mean of each subgroup that results of her factor level combination

- ... and is the same for every subgroup member

- unequal subgroup sizes result in different overall probabilities to be member of one group compared to another

- as group membership defines predicted values, some preditions are more probable than others per se, 
  without knowing anything of an observations group membership
  
- various manners to calculate sums of squares try to deal with this problem
  SS type 3 estimate the unique part a factor accounts for
  but does not attribute all explainable variance to effects.
  
- sequence of taking factors into account in our lm formula in combination with type of SS matters for the effect, 
  a factor accounts for  

- type argument defines type of SS in `car::Anova(..., type = ...)` and `ezANOVA(..., type = ...)`
  

```{r check}
```



`r if(!params$todo) {"\\begin{comment}"}`
todo: 
unbalanced einfaktoriell: Problem sind ggf. ungleiche Varianzen in den ungleich großen Gruppen. Weighted means or grand mean.
Ansonsten sind ungleiche Zellbesetzungen bei mehrfaktoriellen Designs ein Problem.

`r if(!params$todo) {"\\end{comment}"}`


Version: `r format(Sys.time(), '%d %B, %Y %H:%M')`
