##-- Statistik-Peter
##-- Wolfgang Peter
##-- 18-05-2010

##-- Literatur  Angewandte Statistik : Methodensammlung mit R;Lothar Sachs; Jürgen Hedderich;Berlin : Springer Berlin, 2009. 13., aktualis. u. erw. Aufl
##-- Seite 543 bis 547
##-- http://gribblelab.org/2009/03/09/repeated-measures-anova-using-r/
#--  library(ez)
#--  http://bm2.genes.nig.ac.jp/RGM2/pkg.php?p=ez
library(reshape)
reset<-   options()
options(contrasts=c("contr.sum","contr.poly"))    #-- Type III (Berechnung der Quadratsummen wie SPSS)
 
 
 Lines <- "ID   Therapie RMT0 RMT1 RMT2 RMT3 RMT4
1      Nein  400 1200 1950  300  150
2      Nein  790  590  800  400  150
3      Nein 1000 1100  800  500  180
4        Ja 1000  650  530  450  600
5        Ja  150  400  100  150  150
6        Ja  580  530  140  300   20
7      Nein  100   50  300  340  200
8        Ja  160  900  370  200  100
9        Ja  250  330  160  150  115
10     Nein  400  200  400  500  300
11       Ja  100  390  310  120  120
12     Nein  200  120  770  650 1200
13       Ja  270  260  250  250   50
14       Ja   50  450  800  100   80
15       Ja   60  350  350  100   40
16     Nein  200  200  220  100  100
17     Nein  240  560 1700  420  420
18       Ja 1612  400  420  300  100
19     Nein   50  230  200  600   80
20     Nein  250  500  370  280   50
21       Ja  700  100  500  550  550
22     Nein  100  200 1100  300  200
23       Ja 1150  240  900  300  100
24     Nein  250  600  400  100  100
25     Nein  140  450  800  590  400
26       Ja   50  100  580  380  300
27       Ja   50  700  250  250  250
28     Nein  500  950  500  500  150
29     Nein  330  500  100  300  200
30       Ja  450  600  600  450  200
31     Nein  480  630  150  800  680
32     Nein  250  310  310  280  100
33       Ja  250  800  500  300  100
34     Nein  700 1300  610  600  670
35       Ja  970  800  650  600  500
36     Nein   50 1000 1500  400  300
37     Nein  900  320  500  400  300
38       Ja   50  420   10  200  300
39       Ja  240    0  500  200  200
40     Nein  500  400  200   20  230
41     Nein  400  600  400  300  190
42     Nein  400  800  400  800  500
43     Nein  450  100 1250  200  180
44     Nein  100   50  650  500  400
45       Ja  350  500  400  200  100
46     Nein  350  600  290  290  290
47     Nein  350 1200  200  210  200"
DF <- read.table(textConnection(Lines), header = TRUE, na.strings="<NA>")

Zeit <- c(0, 1, 2, 3, 4)
Restmengen<-c("RMT0", "RMT1" ,"RMT2", "RMT3", "RMT4")

#--  Area Under Curve  Quelle Sachs  http://www.springer.com/statistics/statistical+theory+and+methods/book/978-3-540-88901-4
AUC   <- function(y, t, n) {
         F <- rep(NA, (n-1))
         for (i in 1:(n-1)) F[i] <- (t[i+1]-t[i])*(y[i]+y[i+1])
         sum(F)/2          }
#-- Regressions-Koeffizient Quelle Sachs http://www.springer.com/statistics/statistical+theory+and+methods/book/978-3-540-88901-4
REGR <- function(y, t, n) {sum((y-mean(y))*(t-mean(t)))/sum((t-mean(t))^2)}

#--
DF$Max   <- apply(DF[,Restmengen[-1]], 1, max, na.rm = T)   # Maximum
DF$AUC   <- apply(DF[,Restmengen], 1, AUC,  t=Zeit, n=5)    # AUC
DF$REGR  <- apply(DF[,Restmengen], 1, REGR, t=Zeit)         # Regression


attach(DF)
 # Maximum
max <- with(t.test(Max ~ Therapie, data = DF), c(estimate,statistic,p.value) )
 # AUC
auc <-  with(t.test(AUC ~ Therapie, data = DF), c(estimate, statistic, p.value) )
 # Aenderung
reg <- with(t.test(REGR ~ Therapie, data = DF), c(estimate, statistic, p.value) )
 # erste Wert
fst <- with(t.test(RMT1 ~ Therapie, data = DF), c(estimate,statistic,p.value) )
# letzter Wert
lst <- with(t.test(RMT4 ~ Therapie, data = DF), c(estimate, statistic, p.value) )

#tab <- as.data.frame(rbind(max,auc,reg,fst,lst), row.names = c("Max","AUC","REGR","Startwert", "Endwert"))

Berechnung <-  as.data.frame(rbind(max,auc,reg,fst,lst), row.names = c("Max","AUC","REGR","Startwert", "Endwert") )
names(Berechnung)  <-  c("Therapie Ja","Therapie Nein","T-Wert","p-Wert")
Berechnung
#write.csv2(Berechnung,file="DF.Berechnung.csv")

#-- Restrukturieren der Datentabelle
DF.r  <- reshape(DF, varying=list(Restmengen),
             v.names=c("Messwert"), timevar="Zeit",
             times=Zeit, idvar="Proband", direction="long")
DF.r<-transform(DF.r
            ,Zeit    = as.factor(Zeit)
            ,Therapie = as.factor(Therapie)
            ,Proband  = as.factor(Proband))
#-- Mittelwerte
cast(DF.r[-length(DF.r)],Zeit~Therapie,c(length,mean,sd))
#write.csv2(cast(DF.r[-length(DF.r)],Zeit~Therapie,c(length,mean,sd)),file="DF.cast.csv")


par(mfrow=c(1,1), lwd=1.7, font.axis=2, bty="l", ps=12)
with(DF.r,
        interaction.plot(Zeit, Therapie, Messwert, lty=c(1,12), lwd=3
        ,ylab="Mittelwert", xlab="Zeit", trace.label="Therapie")        )

#-- ANOVA für Messwiederholungen


# Type I
#rep.aov <-with(DF.r,
# aov(Messwert ~ Therapie * Zeit )   )
#summary(rep.aov)


rep.aov <-with(DF.r,
                    aov(Messwert ~ Therapie * Zeit + Error(Proband))   )
summary(rep.aov)




#--------------------------------------------------
#-- Alternativer wer für mwsswiederholung   + Post Hoc
#-- Entnommen Paul Gribble Repeated Measures ANOVA using R.
mlm1 <- lm(as.matrix(DF[,Restmengen]) ~ Therapie + 1)
mlm1

Zeit <- factor(Restmengen)

 
library(car)

mlm1.aov <- Anova(mlm1, idata = data.frame(Zeit), idesign = ~ Zeit, type="III")
summary(mlm1.aov, multivariate=FALSE)

#-------------------------------------------------------------------------------
# + Post Hoc
require(nlme)
am2 <- lme(Messwert ~ Zeit +Therapie, random = ~1|ID/Zeit, data=DF.r)
 summary(am2)
  anova(am2)
 require(multcomp)
 summary(glht(am2,linfct=mcp(Zeit="Tukey")))
#-- Kontraste Händisch 
contrastD <- DF$RMT2  -DF$RMT4 
t.test(contrastD)