# Demonstration: Zweifaktorielle Varianzanalyse
# (vgl. Welkowitz, Ewen & Cohen, 2002, p. 246 ff.)

# ------------------------------------------------------------------------------
BaseURL="http://www2.jura.uni-hamburg.de/instkrim/kriminologie/Mitarbeiter/Enzmann/Lehre/StatIIKrim/Uebungen/"

# ------------------------------------------------------------------------------
# Gegebenfalls die folgende Anweisung auskommentieren und den Pfadnamen (mit ab-
# schließendem Schrägstrich-Vorwärts) eintragen, z.B.:

# Pfad = 'd:/winword/texte/S2_ss06/StatIIKrim/Uebungen/'

# ------------------------------------------------------------------------------
# Nur bei Internet-Anschluss verwenden:
download.file(paste(BaseURL,'descby.r',sep=""),paste(Pfad,'descby.r',sep=""),mode='wb')

# ------------------------------------------------------------------------------

# Die Datei muss im aktuellen Arbeitsverzeichnis vorliegen!
source(paste(Pfad,'descby.r',sep=""))

# -----------------------------------------------------------------------------

# Beipiel: Messwerte für eine zweifaktorielle Varianzanalyse

y = c(c( 6,15,12,12,13),  #   male, caffeine large
      c(12,10,12,13, 7),  #   male, caffeine moderate
      c(10,13,15,12,10),  #   male, caffeine small
      c( 9,10, 7,12, 7),  #   male, caffeine zero
      c(10,13, 4, 9, 5),  # female, caffeine large
      c( 9, 7,10, 7,13),  # female, caffeine moderate
      c(12,13,15,10,13),  # female, caffeine small
      c( 4, 7, 6, 9, 9))  # female, caffeine zero

# Die folgenden Variablen repräsentieren die Gruppenzugehörigkeit:

male = factor(rep(1:2,each=20))
levels(male)=c("male","female")
caffeine = factor(c(rep(1:4, each=5), rep(1:4, each=5)))
levels(caffeine)=c("large","moderate","small","zero")

# -----------------------------------------------------------------------------
# Deskriptive Statistiken der Zellen:

stats.males=descby(y[male=='male'],caffeine[male=='male'])
stats.males
stats.females=descby(y[male=='female'],caffeine[male=='female'])
stats.females

# Mittelwerte in einer Matrix:

means=matrix(c(stats.males[,1],stats.females[,1]),nrow=2,ncol=4,byrow=T,
             dimnames=list(levels(male),levels(caffeine)))
means

# Mittelwerte für die Faktoren Geschlecht und Koffein:
descby(y,male)
descby(y,caffeine)

# -----------------------------------------------------------------------------
# Plot der y-Werte der Gruppen mit Mittelwertslinien:

windows(record=T)

g = rep(1:8,each=5)

# sortiert  nach Geschlecht:

glabels = c('m-large','m-mod.','m-small','m-zero',
            'f-large','f-mod.','f-small','f-zero')

plot(y~jitter(g,factor=0.5,amount=0),axes=F,xlab='Gruppen',ylab='Testscore',
     main='Werte und Mittelwerte nach Geschlecht')
axis(side=1,at=1:8,lab=glabels)
axis(side=2)

for (i in 1:8)
{
  # Vertikale Linien vom Minimum zum Maximum pro Gruppe:
  lines(c(i,i),range(y[g==i]),lty=2,col="blue")
  # kleine horizontale Mittelwertslinien pro Gruppe:
  lines(c(i-0.1,i+0.1),c(mean(y[g==i]),mean(y[g==i])),lwd=2,col="blue")
}

# Horizontale Linie der Mittelwerte pro Geschlecht:

lines(c(0.5,4.5),c(mean(y[male=='male']),
                   mean(y[male=='male'])),lwd=1,col="blue")
lines(c(4.5,8.5),c(mean(y[male=='female']),
                   mean(y[male=='female'])),lwd=1,col="blue")
legend(x='topright',legend='gender',lty=1,col='blue',cex=0.8)
# Horizontale Linie des Gesamtmittelwerts:

abline(h=mean(y),lty=2,lwd=2,col="red")

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

# sortiert nach Koffein-Bedingung:

glabels = c('large-m','large-f',
            'mod.-m', 'mod.-f',
            'small-m','small-f',
            'zero-m', 'zero-f')

y.cm = y[order(caffeine,male)]

plot(y.cm~jitter(g,factor=0.5,amount=0),axes=F,xlab='Gruppen',ylab='Testscore',
     main='Werte und Mittelwerte nach Koffein-Bedingung')
axis(side=1,at=1:8,lab=glabels)
axis(side=2)

for (i in 1:8)
{
  # Vertikale Linien vom Minimum zum Maximum pro Gruppe:
  lines(c(i,i),range(y.cm[g==i]),lty=2,col="blue")
  # kleine horizontale Mittelwertslinien pro Gruppe:
  lines(c(i-0.1,i+0.1),c(mean(y.cm[g==i]),mean(y.cm[g==i])),lwd=2,col="blue")
}

# Horizontale Linie der Mittelwerte pro Koffein-Bedinung:

lines(c(0.5,2.5),c(mean(y[caffeine=='large']),
                   mean(y[caffeine=='large'])),lwd=1,col="blue")
lines(c(2.5,4.5),c(mean(y[caffeine=='moderate']),
                   mean(y[caffeine=='moderate'])),lwd=1,col="blue")
lines(c(4.5,6.5),c(mean(y[caffeine=='small']),
                   mean(y[caffeine=='small'])),lwd=1,col="blue")
lines(c(6.5,8.5),c(mean(y[caffeine=='zero']),
                   mean(y[caffeine=='zero'])),lwd=1,col="blue")
legend(x='topright',legend='caffeine',lty=1,col='blue',cex=0.8)

# Horizontale Linie des Gesamtmittelwerts:

abline(h=mean(y),lty=2,lwd=2,col="red")

# ------------------------------------------------------------------------------
# Zweifaktorielle Varianzanalye:

summary(aov(y~caffeine*male))

# Hinweis: Vollständige Modelle mit Haupteffekten und Interaktionstermen werden
# in R in der Form f1*f2 (wobei f1 und f2 die jeweiligen Haupteffekte dar-
# stellen) spezifiziert. Man könnte alle Effekte einzeln spezifizieren, wobei
# mit f1:f2 der Interaktioneffekt gekennzeichnet wird):

summary(aov(y~caffeine+male+caffeine:male))

# oder mittels linearer Regression (da die Gruppenvariablen als "Faktoren" de-
# klariert wurden, s.o, werden sie werden automatisch dummy-codiert) und an-
# schließender Prüfung der F-Werte:

mod=lm(y~caffeine*male)
summary(mod)
anova(mod)

# ------------------------------------------------------------------------------
# Post-hoc Tests der Mittelwertsunterschiede (hier nur der Faktor caffeine):

oldpar=par(las=1,mar=c(5,8,4,2))

PostHoc = TukeyHSD(aov(y~caffeine))
PostHoc
plot(PostHoc)

par(oldpar)

# ------------------------------------------------------------------------------
# Plot zur Darstellung des Interaktionseffekts:

# Eine einfache Methode, Interaktionseffekte grafisch darzustellen, bietet die
# Funktion interaction.plot():

interaction.plot(caffeine,male,y,lty=c(1,2),col="blue",main='Interaktionsplot')

# Man kann die Grafik auch mit Hilfe der Funktion plot() erstellen:

plot(1:4,means[1,],ylim=range(means),type='b',col='blue',axes=F,
     ylab="Mean test score",xlab="Caffeine condition",
     main="Interaktionsplot")
par(new=T)
plot(1:4,means[2,],ylim=range(means),type='b',col='blue',lty=2,axes=F,
     ylab="Mean test score",xlab="Caffeine condition",
     main="Interaktionsplot")
axis(side=1,at=1:4,lab=levels(caffeine))
axis(side=2)
legend(x='topleft',legend=levels(male),lty=c(1,2),col=c('blue','blue'),cex=0.8)

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Interaktionsplot ohne Interaktionseffekt (fiktiv):

means.m = c(11.6,10.8,12.0,9.0)
means.f = c(10.0, 9.2,10.4,7.4)

plot(1:4,means.m,ylim=range(means),type='b',col='blue',axes=F,
     ylab="Mean test score",xlab="Caffeine condition",
     main="Interaktionsplot (ohne Interaktionseffekt)")
par(new=T)
plot(1:4,means.f,ylim=range(means),type='b',col='blue',lty=2,axes=F,
     ylab="Mean test score",xlab="Caffeine condition",
     main="Interaktionsplot (ohne Interaktionseffekt)")
axis(side=1,at=1:4,lab=levels(caffeine))
axis(side=2)
legend(x='topleft',legend=levels(male),lty=c(1,2),col=c('blue','blue'),cex=0.8)

for (i in 1:4)
{
  lines(c(i,i),c(range(means.m[i],means.f[i])[1]+.1,
                 range(means.m[i],means.f[i])[2]-.1))
  text(x=c(i,i),y=c(abs(mean(range(means.m[i],means.f[i]))),
                    abs(mean(range(means.m[i],means.f[i])))),
       lab=paste(round(means.m[i]-means.f[i],1)),
       adj=1,cex=0.8)
}

# Interaktionsplot mit Interaktionseffekt (fiktiv):

means.m = c(11.6,10.8,12.0, 9.0)
means.f = c( 8.2, 9.2, 7.0,12.6)

plot(1:4,means.m,ylim=range(means),type='b',col='blue',axes=F,
     ylab="Mean test score",xlab="Caffeine condition",
     main="Interaktionsplot (mit Interaktionseffekt)")
par(new=T)
plot(1:4,means.f,ylim=range(means),type='b',col='blue',lty=2,axes=F,
     ylab="Mean test score",xlab="Caffeine condition",
     main="Interaktionsplot (mit Interaktionseffekt)")
axis(side=1,at=1:4,lab=levels(caffeine))
axis(side=2)
legend(x='topleft',legend=levels(male),lty=c(1,2),col=c('blue','blue'),cex=0.8)

for (i in 1:4)
{
  lines(c(i,i),c(range(means.m[i],means.f[i])[1]+.1,
                 range(means.m[i],means.f[i])[2]-.1))
  text(x=c(i,i),y=c(abs(mean(range(means.m[i],means.f[i]))),
                    abs(mean(range(means.m[i],means.f[i])))),
       lab=paste(round(means.m[i]-means.f[i],1)),
       adj=1,cex=0.8)
}