# Example 2.4 (Table 2.7) in Davis, 2002

week <- c(2,4,6,8)
m <- matrix(scan("e:\\vhm\\vhm882\\data\\leppik.dat"),ncol=6,byrow=T)
group <- m[,1]
subject <- m[,2]
y <- m[,3:6]
sumy <- c()
mediany <- c()
for (i in 1:59) {
  sumy[i] <- sum(y[i,])
  mediany[i] <- median(y[i,])
  }
median(sumy[group==0])
median(sumy[group==1])
wilcox.test(sumy ~ group)
median(mediany[group==0])
median(mediany[group==1])
wilcox.test(mediany ~ group)

# code for graphics, using nlme library
library(nlme)
leppik <- balancedGrouped( y ~ week|subject, matrix(m[,3:6],nrow=59,ncol=4,dimnames=list(subject,week)),
                           labels=list(y="Number of seizures during two weeks"))
group <- rep(group,rep(4,59))
p1 <- plot(leppik, outer=~as.factor(group), aspect=1)
p1 # screen plot

ps.options(horizontal=FALSE,onefile=FALSE,print.it=FALSE) # p.75 in VR
library(lattice)
trellis.device("postscript",file="c:\\data.avc\\teaching\\vhm882\\01l\\d_eg2_4.eps")
p1 # file plot
dev.off()