# ---------------------------------------------------------------
# Beispiele aus NDK-Skript "Repeated Measures". Autor: H.-R. Roth
# ---------------------------------------------------------------

  source("ftp://stat.ethz.ch/NDK/Source-NDK-7/va2-fkt.R")

# ---------------------------------------------------------------
# Beispiel 1 
# ---------------------------------------------------------------
  clotting <- matrix(c(   
     8.4,  9.4,  9.8, 12.2,
    12.8, 15.2, 12.9, 14.4,
     9.6,  9.1, 11.2,  9.8,
     9.8,  8.8,  9.9, 12.0,
     8.4,  8.2,  8.5,  8.5,
     8.6,  9.9,  9.8, 10.9,
     8.9,  9.0,  9.2, 10.4,
     7.9,  8.1,  8.2, 10.0),8,4,byrow=T)

  d1 <- data.frame(
     y      = as.vector(t(clotting)),
     person = factor(rep(1:8,each=4),labels=c("p1","p2","p3","p4","p5","p6","p7","p8")),
     behand = factor(rep(1:4,length=32),labels=c("A","B","C","D")))

#
# Residuenanalyse für die Zielgrösse y ruft nach einer Transformation
#
  fity <- aov(y ~ person + behand, data=d1)
  par(mfrow=c(1,2))
  plot(fity$fitted, fity$residuals, ylab="residuals", xlab="fitted", main="T-A Plot") 
  abline(h=0,lty=2)
  qqnorm(fity$residuals); qqline(fity$residuals)

#
# neue Zielgrösse ist 100/y
# -------------------------

  par(mfrow=c(1,1))
  interaction.plot(d1$behand, d1$person, 100/d1$y, main="Beispiel 1", 
                   xlab="Behandlung", ylab="Gerinnungsrate [%]", trace.label="Person")

  fit <- aov(100/y ~ person + behand, data=d1)
  par(mfrow=c(1,2))
  plot(fit$fitted, fit$residuals, ylab="residuals", xlab="fitted", main="T-A Plot")
  abline(h=0,lty=2)
  qqnorm(fit$residuals); qqline(fit$residuals)
  
  summary(fit)

  model.tables(fit,type="means",se=T)

#
#  paarweise Vergleiche mit der Funktion TukeyHSD
#
  TukeyHSD(aov(100/y ~ person + behand, data=d1))[2]
  
#
#  paarweise Vergleiche mit Funktion aus dem package "multcomp"
#

#  Das Package multcomp gehört nicht zu den standard Packages. Sie finden
#  es auf der Homepage http://stat.ethz.ch/CRAN/ und müssen es selber installieren.
#
  library(multcomp)
  
  options(contrasts=c("contr.treatment","contr.poly"))
  
  fit  <- aov(100/y ~ behand + person, data=d1)      # Faktor behand zuerst
  sulm <- summary.lm(fit,correlation=T)

  param <- sulm$coefficients[,1][1:4]
  cov   <- sulm$sigma^2*sulm$cov.unscaled[1:4,1:4]

  param[1] <- 0    # da fit mit "treatment"-Kontrast berechnet
  cov[1,]  <- 0
  cov[,1]  <- 0

  con <- contrMat(rep(8,4),type="Tukey")    # alle Paarvergleiche
  
  summary(csimint(estpar=param, df=fit$df, covm=cov, cmatrix=con))
  
# ---------------------------------------------------------------
# Beispiel 2 
# ---------------------------------------------------------------
  Y <- matrix(c(
    43,  90,  51,  67,
    87,  36,  12,  14,
    18,  56,  22,  68,
    34,  73,  34,  87,
    81,  55,  29,  54,
    45,  58,  62,  44,
    16,  35,  71,  37,
    43,  47,  87,  27,
    22,  91,  37,  78,
    10,  81,  43,  33,
    58,  84,  35,  43,
    26,  49,  55,  84,
    18,  30,  49,  44,
    13,  14,  25,  45,
    12,   8,  40,  48,
     9,  55,  10,  30,
    31,  45,   9,  66),17,4, byrow=T)
    
#
# 4-dim. Zielgrösse
#
  d2m <- data.frame(
    y1     = Y[,1],
    y2     = Y[,2],
    y3     = Y[,3],
    y4     = Y[,4],
    person = factor(1:17),
    group  = factor(rep(c("A","B"),c(9,8))))
    
#
# 1-dim. Zielgrösse
#
  d2u <- data.frame(
    y       = as.vector(t(Y)),
    person  = factor(rep(1:17,each=4)),
    group   = factor(rep(c("A","B"),c(9*4,8*4))),
    problem = factor(rep(1:4,length=68)))
        
#
# graphische Darstellung der Profile
#
  par(mfrow=c(1,1))
  interaction.plot(d2u$problem, d2u$group, d2u$y, fun=median, main="Beispiel 2",
                   xlab="Problem", ylab="Zeit", trace.label="Gruppe")

  par(mfrow=c(2,1))
  for(i in c("A","B")) 
  {
  ss <- subset(d2u, group==i)
  interaction.plot(ss$problem, factor(ss$person), ss$y, ylim=range(d2u$y), 
                 xlab="Problem", ylab="Zeit", trace.label="Person")      
  }  
                 
#
# Mittelwerte tabellarisch
#
  model.tables(aov( y ~ group*problem, data=d2u),type="means",se=T)

  model.tables(aov(log(y)~group*problem, data=d2u),type="means",se=T)

#
# Diverse Lineare Modelle für log(y)
# ----------------------------------
  options(contrasts=c("contr.sum","contr.poly"))
  
#
# Type I SS
#
  fit <- aov(log(y) ~ group*problem + Error(person%in%group), data=d2u)
  summary(fit)
#
# Type III SS
#
  summary(fit,ssType=3)    # nur in S-Plus möglich, R ignoriert einfach!

  fit2 <- aov(log(y) ~ group*problem + person%in%group, data=d2u)
  
  drop1(fit2, .~., test="F")     # une catastrophe !!!


  library(nlme)
  
  options(contrasts=c("contr.sum","contr.poly"))

  lmefit <- lme(log(y) ~ group*problem, data=d2u, random = ~1 | person)
  anova(lmefit)
  
  anova(lmefit,type="marginal")   # nur mit contr.sum richtige Resultate!!!

  intervals(lmefit,which="var-cov")
  random.effects(lmefit)[[1]]

  par(mfrow=c(1,2))
  plot(fitted.lme(lmefit), residuals.lme(lmefit)); abline(h=0,lty=2)
  qqnorm(residuals.lme(lmefit)); qqline(residuals.lme(lmefit))

# --------------------------------------------------------------------------
# Beispiel 3 
# --------------------------------------------------------------------------
  P <- matrix(c(
    1.0,   NA, 10.0, 16.0, 22.0, 20.0, 16.0,   NA, 18.0, 14.0,
    6.5,  5.7,  9.5, 11.6, 17.5, 27.3, 28.5, 22.4, 19.3, 10.0,
    3.0,  4.0,  4.0, 13.0, 15.8, 19.5, 21.2, 17.9, 10.7, 13.4,
    1.0,  2.1,  9.7,   NA, 21.8,   NA, 27.5,   NA, 15.5,  6.2,
    1.0,  1.0,  1.0,  4.2, 22.6, 23.9, 45.5, 42.6, 35.0, 10.6,
    1.0,  1.0,  1.0,  1.0,  3.9, 14.7, 17.6, 16.1,  8.8, 10.8,
    1.0,  1.5,  5.0, 11.0, 16.0, 23.0, 15.0,  9.0,  6.0,  5.0,
    1.0,  1.0,  6.5, 20.0, 22.5, 27.8, 19.0,  9.0,  8.2,  8.0,
    1.0,  1.0,  7.3,  7.5, 18.0, 20.0, 18.9, 12.8,  6.3,  4.8,
    3.0,  2.5,  2.0,  2.7,  3.4,  3.6, 14.0,  7.3,  7.7,  4.7,
    8.3,  7.5,  9.6, 11.0, 11.5, 15.7, 15.2, 15.8, 14.0, 11.5,
    6.2,  5.9,  6.8,  7.7,  9.0,  9.3, 12.1, 12.2, 11.0,  9.0,
    8.4, 10.8,  8.1,  7.8,  8.5, 12.0, 19.8, 22.2, 25.2, 40.5,
    3.5,  3.2,  3.4,  3.3,  8.5,  9.4, 14.5, 12.7, 11.5, 10.2,
    3.5,  4.0,  4.8,  3.5,  3.7, 13.0, 12.5, 15.0, 22.0, 10.5,
    3.7,  3.2,  4.3,  4.5,  5.5,  8.5, 10.3, 11.1,  8.0,  6.0,
    5.0,  5.6,  6.1,  7.2, 13.8, 26.0, 26.1, 25.7, 20.5, 11.0,
    4.5,  5.1, 13.2, 21.0, 26.8, 28.0, 22.0, 17.8, 15.7, 14.0,
    8.4,  6.2,  8.0, 18.5, 33.8, 35.0, 26.2, 23.0, 19.0, 12.6,
    4.2,  3.2,  4.2,  4.8, 10.3, 13.7, 17.1, 18.3, 17.4, 15.8),20,10, byrow=T)

  P <- log(P)

#
# 10-dim. Zielgrösse
#
  d3m <- data.frame(
    p1     = P[,1],
    p2     = P[,2],
    p3     = P[,3],
    p4     = P[,4],
    p5     = P[,5],
    p6     = P[,6],
    p7     = P[,7],
    p8     = P[,8],
    p9     = P[,9],
    p10    = P[,10],
    person = factor(1:20),
    group  = factor(rep(c("A","B","C","D"),c(6,6,4,4))))
    
#
# 1-dim. Zielgrösse
#
  d3u <- data.frame(
    y       = as.vector(t(P)),
    person  = factor(rep(1:20,each=10)),
    group   = factor(rep(c("A","B","C","D"),c(60,60,40,40))),
    zeit    = factor(rep(c(0,1,3,5,10,15,30,45,60,120),20)))
        
#
# graphische Darstellung der Profile
# ----------------------------------

  y.range <- range(d3u$y, na.rm=T)
  par(mfrow=c(2,2))
  for(i in c("A","B","C","D")) 
   {
   ss <- subset(d3u, group==i)
   interaction.plot(ss$zeit, factor(ss$person), ss$y, main="Beispiel 3",
                    ylim=y.range, xlab="Zeit", ylab="Konz.", trace.label="Person")      
   }  
   
  par(mfrow=c(1,1))
  
  f.st <- function(x) { median(x, na.rm=T) }
  
  interaction.plot(d3u$zeit, d3u$group, d3u$y, fun=f.st, main="Beispiel 3", 
                   xlab="Zeit", ylab="Konz.", trace.label="Gruppe")
#
# Person 13 (Ausreisser) wird von allen weiteren Berechnungen ausgeschlossen
#
  ssd3u <- subset(d3u,!is.na(d3u$y) & person!=13)
#
# Mittelwerte tabellarisch
#

  model.tables(aov(y~group*zeit, data=ssd3u),type="means")

#
# AUC für jedes Profil
#
  d3m$auc <- f.auc(P[,1:10], time=c(0,1,3,5,10,15,30,45,60,120))

#
# TTP: Time to Peak für jedes Profil
#
  time    <- c(0,1,3,5,10,15,30,45,60,120)
  pos     <- apply(P[,1:10], 1, which.max)
  d3m$ttp <- time[pos]

#
# Peak für jedes Profil
#
  d3m$peak <- P[cbind(1:length(pos),pos)]

#
# Tabelle und Grafiken der Profil-Eigenschaften
#
  print(d3m[,11:15],digits=4)
  
  par(mfrow=c(2,2))
  plot.default(d3m$group, d3m$auc,  xlab="Gruppe", ylab="AUC")
  plot.default(d3m$group, d3m$ttp,  xlab="Gruppe", ylab="TTP")
  plot.default(d3m$group, d3m$peak, xlab="Gruppe", ylab="PEAK")

#
# ANOVA's für AUC, TTP und PEAK
# -----------------------------
  par(mfrow=c(3,2))
 
  summary(fit.auc <- aov(d3m$auc ~ group, data=d3m, subset=person!=13))
 
  plot(fit.auc$fitted, fit.auc$residuals, ylab="AUC-residuals"); abline(h=0,lty=2)
  qqnorm(fit.auc$residuals); qqline(fit.auc$residuals)


  summary(fit.ttp <- aov(d3m$ttp ~ group, data=d3m, subset=person!=13))
 
  plot(fit.ttp$fitted, fit.ttp$residuals, ylab="TTP-residuals"); abline(h=0,lty=2)
  qqnorm(fit.ttp$residuals); qqline(fit.ttp$residuals)


  summary(fit.peak <- aov(log(d3m$peak) ~ group, data=d3m, subset=person!=13))
 
  plot(fit.peak$fitted, fit.peak$residuals, ylab="PEAK-residuals"); abline(h=0,lty=2)
  qqnorm(fit.peak$residuals); qqline(fit.peak$residuals)
  par(mfrow=c(1,1))

#
# MANOVA für AUC, TTP und PEAK (ohne Daten von Person 13)
# -------------------------------------------------------
  d3m <- d3m[d3m$person!=13,]
  
  X <- as.matrix(d3m[,13:15])
  summary(manova(X ~ d3m$group))
  

# --------------------------------------------------------------
# Beispiel 4 
# --------------------------------------------------------------
  P <- matrix(c(
    21.0, 20.0, 21.5, 23.0,
    21.0, 21.5, 24.0, 25.5,
    20.5, 24.0, 24.5, 26.0,
    23.5, 24.5, 25.0, 26.5,
    21.5, 23.0, 22.5, 23.5,
    20.0, 21.0, 21.0, 22.5,
    21.5, 22.5, 23.0, 25.0,
    23.0, 23.0, 23.5, 24.0,
    20.0, 21.0, 22.0, 21.5,
    16.5, 19.0, 19.0, 19.5,
    24.5, 25.0, 28.0, 28.0,
    26.0, 25.0, 29.0, 31.0,
    21.5, 22.5, 23.0, 26.5,
    23.0, 22.5, 24.0, 27.5,
    25.5, 27.5, 26.5, 27.0,
    20.0, 23.5, 22.5, 26.0,
    24.5, 25.5, 27.0, 28.5,
    22.0, 22.0, 24.5, 26.5,
    24.0, 21.5, 24.5, 25.5,
    23.0, 20.5, 31.0, 26.0,
    27.5, 28.0, 31.0, 31.5,
    23.0, 23.0, 23.5, 25.0,
    21.5, 23.5, 24.0, 28.0,
    17.0, 24.5, 26.0, 29.5,
    22.5, 25.5, 25.5, 26.0,
    23.0, 24.5, 26.0, 30.0,
    22.0, 21.5, 23.5, 25.0),27,4, byrow=T)

#
# 4-dim. Zielgrösse
#
  d4m <- data.frame(
     p1     = P[,1],
     p2     = P[,2],
     p3     = P[,3],
     p4     = P[,4],
     person = factor(1:27),
     group  = factor(rep(c("G","B"),c(11,16))))
    
#
# 1-dim. Zielgrösse
#
  d4u <- data.frame(
     y      = as.vector(t(P)),
     person = factor(rep(1:27,each=4)),
     group  = factor(rep(c("G","B"),c(11*4,16*4))),
     age    = rep(c(8,10,12,14),length=108))
        
#
# graphische Darstellung der Profile
#
  par(mfrow=c(1,1))
  titel <- "Beispiel 4: mittlere Profile"
  interaction.plot(d4u$age, d4u$group, d4u$y, fun=median, main=titel,
                   xlab="Alter", ylab="Distanz", trace.label="Gruppe")

  par(mfrow=c(1,2))
  titel <- "Beispiel 4: individuelle Profile"
  for(i in c("G","B")) 
  {
  ss <- subset(d4u, group==i)
  interaction.plot(ss$age, factor(ss$person), ss$y, ylim=range(d4u$y),
                 main=titel, xlab="Alter", ylab="Distanz", trace.label="Person")
  }  

#
# Mittelwerte tabellarisch
#
  model.tables(aov( y ~ group*factor(age), data=d4u),type="means",se=F)

#
# Regressionsparameter für jedes Profil
# -------------------------------------
  for(i in 1:27)
  {
  fit <- lm(y ~ I(age-11), data=d4u, subset=(person==i))
  d4m$alpha[i] <- coef(fit)[[1]]
  d4m$beta[i]  <- coef(fit)[[2]]
  }
 
#
# Tabelle und Grafiken der Regressionskoeffizienten
#
  print(d4m,digits=4)
  
  par(mfrow=c(1,2))
  titel <- "Beispiel 4: individuelle Regressionsparameter"  
  plot(d4m$group,d4m$alpha, xlab="Gruppe", ylab="alpha", main=titel)
  plot(d4m$group,d4m$beta , xlab="Gruppe", ylab="beta", main=titel)

#
# Unterschied zwischen Gruppen mit t-Test
#
  t.test(alpha ~ group, d4m, subset=person!=24, var.equal=T)
  
  t.test(beta ~ group, d4m, subset=person!=24, var.equal=T)

#
# hat Wachstum stattgefunden?
#
  t.test(d4m$beta[1:11], mu=0)                      # Girls
  
  t.test(c(d4m$beta[12:23],d4m$beta[25:27]), mu=0)  # Boys

#
# Repeated Measures ANOVA
# -----------------------

# Kovarianzmatrizen für jede Gruppe separat geschätzt

  M3 <- var(P[1:11,])  # Kovarianzmatrix für Girls
  M4 <- var(P[12:27,]) # Kovarianzmatrix für Boys 
  round(M3,2)
  round(M4,2)
  M5 <- (M3*10 + M4*15) / 25
  S  <- diag(1/sqrt(diag(M5)))
  R5 <- S%*%M5%*%S
  round(M5,2)
  round(R5,3)
  f.epsi(M5,4,2,27)    # compound symmetry check

  d4u$agef <- ordered(d4u$age)

  fit <- aov(y ~ group + agef + group:agef + Error(person%in%group), d4u)
  summary(fit, split=list("group:agef"=list(lin=1,qua=2,cub=3)))

#
# MANOVA
# ------

  M9  <- M5       #  M9  muss nicht unbedingt berechnet werden
  M10 <- var(P)   #  M10 muss nicht unbedingt berechnet werden

#
# Gleichheit der 4-dim Profile
#
  summary(manova(P ~ d4m$group), test="Hotelling-Lawley")

# Parallelität auf 2 Arten getestet
# vgl. C.S. Davis (2002), p60

  c1 <- c(-1,  1,  0, 0)
  c2 <- c( 0, -1,  1, 0)
  c3 <- c( 0,  0, -1, 1)
  C1 <- cbind(c1,c2,c3)
  
  summary(manova(P%*%C1 ~ d4m$group), test="Hotelling-Lawley")
  
  c1 <- c(-3, -1,  1, 3)  # linearer Trend
  c2 <- c( 1, -1, -1, 1)  # quadrat. Trend
  c3 <- c(-1,  3, -3, 1)  # kubisch. Trend
  C2 <- cbind(c1,c2,c3)
  
  mfit <- manova(P%*%C2 ~ d4m$group)
  summary(mfit, intercept=T, test="Hotelling-Lawley")
  
#  
# Unterschied im lin., quadr. und kub. Trend separat getestet
#
  summary.aov(mfit, intercept=T)

#  
# Mixed effects model  bzw. random coefficients model
# ---------------------------------------------------

  library(nlme)
  
  d4u$x <- d4u$age - 11
  
# volles Modell
  fm1 <- lme(y ~ x*group, data=d4u, random = ~ x | person)  
  summary(fm1)
  anova(fm1,type="marginal")
  intervals(fm1)
  
# reduz. Modell (ohne random slope) ist etwas besser
  fm2 <- lme(y ~ x*group, data=d4u, random = ~ 1 | person)    
  summary(fm2)

# ----------------------------------------------------------------
# Beispiel 5 
# ----------------------------------------------------------------
  Y <- matrix(c(
    57, 86, 114, 139, 172,
    60, 93, 123, 146, 177,
    52, 77, 111, 144, 185,
    49, 67, 100, 129, 164,
    56, 81, 104, 121, 151,
    46, 70, 102, 131, 153,
    51, 71,  94, 110, 141,
    63, 91, 112, 130, 154,
    49, 67,  90, 112, 140,
    57, 82, 110, 139, 169,
    59, 85, 121, 156, 191,
    54, 71,  90, 110, 138,
    56, 75, 108, 151, 189,
    59, 85, 116, 148, 177,
    57, 72,  97, 120, 144,
    52, 73,  97, 116, 140,
    52, 70, 105, 138, 171,
    61, 86, 109, 120, 129,
    59, 80, 101, 111, 122,
    53, 79, 100, 106, 133,
    59, 88, 100, 111, 122,
    51, 75, 101, 123, 140,
    51, 75,  92, 100, 119,
    56, 78,  95, 103, 108,
    58, 69,  93, 114, 138,
    46, 61,  78,  90, 107,
    53, 72,  89, 104, 122),27,5, byrow=T)

#
# 5-dim. Zielgrösse
#
  d5m <- data.frame(
    y0     = Y[,1],
    y1     = Y[,2],
    y2     = Y[,3],
    y3     = Y[,4],
    y4     = Y[,5],
    tier   = factor(1:27),
    gruppe = factor(rep(c("G1","G2","G3"),c(10,7,10))))
    
#
# 1-dim. Zielgrösse
#
  d5u <- data.frame(
    y      = as.vector(t(Y)),
    tier   = factor(rep(1:27,each=5)),
    gruppe = factor(rep(c("G1","G2","G3"),c(10*5,7*5,10*5))),
    woche  = rep(-2:2,length=135))

#
# graphische Darstellung der Profile
#
  par(mfrow=c(1,1))
  interaction.plot(d5u$woche, d5u$gruppe, d5u$y, fun=median, main="Beispiel 5",
                   xlab="Woche", ylab="Gewicht", trace.label="Gruppe")

  y.range <- range(d5u$y)
  par(mfrow=c(2,2))
  for(i in c("G1","G2","G3")) 
  {
  ss <- subset(d5u, gruppe==i)
  titel <- paste("Beispiel 5, Gruppe",i)
  interaction.plot(ss$woche, factor(ss$tier), ss$y,  main=titel,
                   ylim=y.range, xlab="Woche", ylab="Gewicht", trace.label="Tier")
  }  

#
# Mittelwerte tabellarisch
#
  model.tables(aov( y ~ gruppe*factor(woche), data=d5u),type="means",se=T)

#
# within groups Kovarianzmatrix für Y, log(Y) und 1. Differenzen
# --------------------------------------------------------------
  M6 <- (var(Y[1:10,])*9 + var(Y[11:17,])*6 + var(Y[18:27,])*9) / 24  

  M7 <- (var(log(Y[1:10,]))*9 + var(log(Y[11:17,]))*6 + var(log(Y[18:27,]))*9) / 24  
  
  D  <- Y[,2:5] - Y[,1:4]
  M8 <- (var(D[1:10,])*9 + var(D[11:17,])*6 + var(D[18:27,])*9) / 24
    
  round(M6,digits=1)
  round(M7,digits=4)
  round(M8,digits=2)
#
# Korrelationsmatrizen
#
  S  <- diag(1/sqrt(diag(M6)))
  R6 <- S%*%M6%*%S
  S  <- diag(1/sqrt(diag(M7)))
  R7 <- S%*%M7%*%S
  S  <- diag(1/sqrt(diag(M8)))
  R8 <- S%*%M8%*%S
  
  round(R6,digits=3)
  round(R7,digits=3)  
  round(R8,digits=3)

#
# compound symmetry checks
#
  f.epsi(M6,5,3,27)    
  f.epsi(M7,5,3,27)    
  f.epsi(M8,4,3,27)  
    
#
# RM-ANOVA mit Split Plot Modell
#
  d5u$wocheF <- ordered(d5u$woche)
  
  fit <- aov(y ~ gruppe*wocheF + Error(tier%in%gruppe), d5u)
  summary(fit, split=list("wocheF"=list(lin=1,qua=2)))
  
#
# Regressionsparameter für jedes Profil
# -------------------------------------

  for(i in 1:27)
  {
  fit <- lm(y ~ I(woche-2) + I((woche-2)^2), data=d5u, subset=(tier==i))
  d5m$beta[i]  <- coef(fit)[[2]]
  d5m$gamma[i] <- coef(fit)[[3]]
  }
 
  par(mfrow=c(1,2))
  plot(d5m$gruppe,d5m$beta, xlab="Gruppe", ylab="beta", main="linearer Trend")
  plot(d5m$gruppe,d5m$gamma , xlab="Gruppe", ylab="gamma", main="quadratischer Trend")

#
# Unterschied zwischen den Gruppen
#
  summary(aov(beta ~ gruppe, d5m))
  
  summary(aov(gamma ~ gruppe, d5m))

#
# MANOVA's
# ---------

  fit.mult <- manova(Y ~ d5m$gruppe)
  summary(fit.mult, test="Wilks")
#
# polynomiale Kontraste
#
  fit.mult <- manova(Y%*%contr.poly(5) ~ d5m$gruppe)
  summary(fit.mult, test="Wilks", intercept=T)
  summary.aov(fit.mult, intercept=T)

# --------------------------------------------------------------
# Beispiel 6 
# --------------------------------------------------------------
  
  y11 <- c(310, 310, 370, 410, 250, 380, 330)
  y12 <- c(270, 260, 300, 390, 210, 350, 365)
  y21 <- c(370, 310, 380, 290, 260,  90) 
  y22 <- c(385, 400, 410, 320, 340, 220)

  p1 <- c(1,4,6,7,10,11,14)
  p2 <- c(2,3,5,9,12,13)

#
# 2-dim. Zielgrösse
#
  d6m <- data.frame(
     y1      = as.vector(c(y11,y21)),
     y2      = as.vector(c(y12,y22)),
     person  = factor(c(p1,p2)),
     sequenz = factor(rep(c("AB","BA"),c(7,6))))
     
  d6m$differenz <- d6m$y1 - d6m$y2
  d6m$summe     <- d6m$y1 + d6m$y2 
         
#
# 1-dim. Zielgrösse
#
  d6u <- data.frame(
     y       = as.vector(c(y11,y12,y21,y22)),
     person  = factor(c(p1,p1,p2,p2)),
     sequenz = factor(rep(c("AB","BA"),c(7*2,6*2))),
     periode = factor(rep(c(1,2,1,2),c(7,7,6,6))),
     behandl = factor(rep(c("A","B","B","A"),c(7,7,6,6))))
        
#
# graphische Darstellung der Profile
#
  par(mfrow=c(1,2))
  titel <- "Beispiel 6: Sequenz"  
  for(i in c("AB","BA")) 
  {
  ss <- subset(d6u, sequenz==i)
  interaction.plot(ss$periode, factor(ss$person), ss$y, ylim=range(d6u$y),
         main=paste(titel,i), xlab="Periode", ylab="PEF", trace.label="Person",cex=0.30)
  }  
                 
#
# Mittelwerte
#
  yq11 <- mean(y11)
  yq12 <- mean(y12)
  yq21 <- mean(y21)
  yq22 <- mean(y22)
  
  print(c(yq11,yq12),digits=4)
  print(c(yq21,yq22),digits=4)
  
#
# geschätzte Effekte
#
  lambda <-  (yq11 + yq12) - (yq21 + yq22)
  pi     <- ((yq11 - yq12) + (yq21 - yq22)) / 2
  alpha  <- ((yq11 - yq12) + (yq22 - yq21)) / 2 
  
  print(c(lambda,pi,alpha),digits=4)
  
#
# Modell 2
#
  summary(aov(y ~ sequenz +Error(person%in%sequenz) + behandl + periode, d6u))
# 
# Modell 3
#
  summary(aov(y ~ behandl:periode +Error(person%in%sequenz) + behandl + periode, d6u)) 

#
# Residuenanalyse  ( aov() ohne Error )
#
  fit <- aov(y ~ sequenz + person%in%sequenz + behandl + periode, d6u)
 
  par(mfrow=c(1,2))
  plot(fit$fitted, fit$residuals); abline(h=0,lty=2)
  qqnorm(fit$residuals); qqline(fit$residuals)