###########################
###########################
##
## R commands for drawing the figures in Chapter 2, Plots
## Modeling Longitudinal Data by Robert Weiss
## Published by Springer in 2005.  
##
##
###########################
###########################

## Postscript commands are commented out, along with the dev.off() command.
## Use those to produce postscript files of the figures.  
## If you want to save a copy of a figure, right click on the plot.  



#
# Pick one of the following commands for reading in the data 
# and adapt it as appropriate for your situation.  Read in the data.  
# Big Mice raw data are in long format.  
# 

#mice <- read.table("C:/rob/mystuff/RMbook/rmnotes/comps/mice/mice.dat",header=T)
mice <- read.table("mice.dat",header=T)

x11(width = 4.5, height = 4.5, pointsize = 10)




## Figure 2.1.  Scatterplot of Y_{ij} against t_{ji} for Big Mice data.  

#postscript(file="C:/rob/mystuff/RMbook/rmnotes/plots/graphics/graphics1.ps",
#onefile=FALSE, horizontal=FALSE, height=4.5, width=4.5, pointsize=10)

plot(mice$time, mice$wt, xlab="", ylab="",axes=F)
# next two commands are not needed for this figure, but used in the text.  
#rect(1.5,110,4.5,450)
#rect(14.5,530,17.5,1160)
axis(side=1, at=0:4*5, lab=0:4*5, cex=.8, font=10)
axis(side=2, at=(0:5*200+200), lab=(0:5*200+200), cex=.8, font=10)
mtext("Days", side=1, line=2.15, cex=.8)
mtext("Weight (mg)", side=2, line=2.2, cex=.8)
box()

#dev.off()




## Figure 2.2, Repeated box plots over time for the Big Mice data.  


#postscript(file="C:/rob/mystuff/RMbook/rmnotes/plots/graphics/graphics2.ps", onefile=FALSE, #horizontal=FALSE, height=4.5, width=4.5, pointsize=10)

boxplot(mice$wt ~ mice$time, range=0,xlab="", ylab="",axes=F)
axis(side=1, at=0:4*5, lab=c(0,5,10,15,20), cex=.8,font=10)
axis(side=2, at=(0:5*200+200), lab=(0:5*200+200),cex=.8,font=10)
mtext("Days", side=1, line=2.15, cex=.8)
mtext("Weight (mg)", side=2, line=2.2, cex=.8)
box()

#dev.off()





##  You will find many variations of this function, but almost all are the
##  same.  The variations are minor.  This plots profile plots.  
##  Much of the potential functionality has been removed in places.  
## 

spagplot24 <- function(time, id, y, xlab="Time",ylab="y",
                     lx=min(time),ux=max(time),
                     ly=min(y,na.rm=T),uy=max(y,na.rm=T),...){
  okidx <- !is.na(y)
  y <- y[okidx]
  id <- id[okidx]
  time <- time[okidx]
  plot(c(lx,ux),c(ly,uy),type="n",xlab="",ylab="",axes=F)
#  title(xlab=xlab,ylab=ylab,...)
  n <- length(time)
  subidx <- unique(id)
  for (i in subidx){
    xtmp <- time[id==i]
    ytmp <- y[id==i]
    points(xtmp,ytmp)
    lines(xtmp,ytmp)
  }
}

## Figure 2.3, Profile plot with dots.  


#postscript(file="C:/rob/mystuff/RMbook/rmnotes/plots/graphics/graphics3.ps",
#onefile=FALSE, horizontal=FALSE, height=4.5, width=4.5, pointsize=10)

spagplot24(mice$time, mice$id, mice$wt, xlab="", ylab="")
axis(side=1, at=0:4*5, lab=0:4*5, cex=.8, font=10)
axis(side=2, at=(0:5*200+200), lab=(0:5*200+200), cex=.8, font=10)
mtext("Days", side=1, line=2.15, cex=.8)
mtext("Weight (mg)", side=2, line=2.2, cex=.8)
box()

#dev.off()

########################
##
##
## Figure 2.4, The need for connecting-the-dots.    
## First read in the raw data.  
##
#########################


t1 <- matrix(data=scan(), nrow=5, ncol=6, byrow=T)
-0.410595 -0.203062 -0.0560626 0.308274 0.361065 0.467333 
-0.461461 -0.220357 -0.0795592 0.117201 0.206216 0.405211 
-0.368854 -0.213922 -0.118478 0.229861 0.387449 0.459431 
-0.427479 -0.36239 -0.152453 0.0417976 0.151389 0.36363 
-0.383159 -0.286105 -0.239626 -0.0751478 0.0563148 0.332521

t2 <- matrix(data=scan(), nrow=5, ncol=6, byrow=T)
-0.410595 -0.203062 -0.0560626 0.117201 0.206216 0.405211
-0.368854 -0.213922 -0.118478 0.0417976 0.151389 0.36363
-0.383159 -0.286105 -0.239626 -0.0751478 0.0563148 0.332521
-0.427479 -0.36239 -0.152453 0.229861 0.387449 0.459431
-0.461461 -0.220357 -0.0795592 0.308274 0.361065 0.467333

y1  <- matrix(data=scan(), nrow=5, ncol=6, byrow=T)
1.14483 0.925515 0.97393 -0.00773519 0.0533554 -0.215457
-0.243821 0.0942599 0.103281 0.896782 0.904834 1.22358
0.838395 0.814355 0.708589 0.352087 0.147781 0.174904
0.210492 0.215218 0.321609 0.728796 0.720638 0.769987
0.518491 0.4772 0.417289 0.442167 0.527747 0.451304

y2  <- matrix(data=scan(), nrow=5, ncol=6, byrow=T)
1.14483 0.925515 0.97393 0.896782 0.904834 1.22358
0.838395 0.814355 0.708589 0.728796 0.720638 0.769987 
0.518491 0.4772 0.417289 0.442167 0.527747 0.451304
0.210492 0.215218 0.321609 0.352087 0.147781 0.174904
-0.243821 0.0942599 0.103281 -0.00773519 0.0533554 -0.215457


# The next two commands will close the current plot window, then 
# open up a slightly larger window.  

dev.off()
x11(width = 4, height = 6, pointsize = 10,)


#postscript(file="C:/rob/mystuff/RMbook/rmnotes/plots/graphics/graphics4.ps",
#onefile=FALSE, horizontal=FALSE, height=6, width=4, pointsize=10)

par(mfrow=c(3,1))
par(mar=c(3.1,4.1,2.1,2.1))
m1 <- c(min(t1-.1), max(t1))
m2 <- c(min(y1-.05),max(y1+.05))
plot(m1, m2, type="n", xlab="", ylab="", axes=F)
  for (i in 1:5){
    points(t1,y1)
#    lines(xtmp,ytmp)
  }
mtext("(a)", side=1, line=2.0, cex=.8)
mtext("Time", side=1, line=.85, cex=.8)
mtext("Response", side=2, line=2.2, cex=.8)
box()

plot(m1, m2, type="n", xlab="", ylab="", axes=F)
    points(t2,y2)
  for (i in 1:5){
    lines(t2[i,],y2[i,])
  }
mtext("(b)", side=1, line=2.0, cex=.8)
mtext("Time", side=1, line=.85, cex=.8)
mtext("Response", side=2, line=2.2, cex=.8)
text(x=(t1[,1]-.1),y=y1[,1],labels=c(1,5,2,4,3))
box()


plot(m1, m2, type="n", xlab="", ylab="", axes=F)
    points(t1,y1)
  for (i in 1:5){
    lines(t1[i,],y1[i,])
  }
mtext("(c)", side=1, line=2.0, cex=.8)
mtext("Time", side=1, line=.85, cex=.8)
mtext("Response", side=2, line=2.2, cex=.8)
text(x=(t1[,1]-.1),y=y1[,1],labels=c(1,5,2,4,3))
box()


#dev.off()

##  Figure 2.5 currently.  
##  4 Example profile plots.  
##  Two ways to read the raw data in. 

tempnames = c("observation", "y1", "y2", "y3", "y4", "id", "time")

temp <- scan()
1  -0.71799387  0.84543866 0.1917018  2.5729763  1 1.391201
2  -1.36601557 -0.18348911 0.0400000  3.2434125  2 1.124349
3   1.37196774 -0.71907638 0.9041763  0.9961150  3 1.970888
4   1.44575446  0.84058525 0.1176840  2.4719014  4 1.240704
5  -0.31138223  1.96273185 1.0152547 -0.2474876  5 2.082694
6  -1.61511875  0.09828049 0.9686641 -1.1824746  6 2.483461
7   0.28581649  1.96496427 0.3096998  3.7729474  7 1.343170
8  -1.03798014 -1.62449804 3.3149045  2.7124752  8 1.483506
9  -0.76658209  1.17121745 0.5227769 -0.7188650  1 1.964677
10 -1.57534871 -0.08883809 0.3752078  4.3328321  2 1.377757
11  1.18946950 -0.05874062 2.1886537 -2.4359183  3 3.379799
12  1.46561834  0.88737469 0.2417865 -0.5509173  4 2.243243
13 -0.42457659  2.00860824 1.1586229 -1.0193918  5 2.529357
14 -1.78584864  0.18026145 1.6860050 -1.7162852  6 3.218415
15  0.23684373  2.33576460 1.1690627 -3.1194098  7 3.012688
16 -0.75608475 -1.22096580 3.7680699         NA  8 2.359724
17 -0.87079928  2.27847648 1.1226281 -1.1456607  1 4.257207
18 -1.55595682  0.28566195 0.5887589 -2.1350411  2 2.373881
19  1.27375728  0.60198931 3.1251655 -2.9460690  3 4.547261
20  1.72378839  1.49498012 0.4334495  0.3279331  4 3.564363
21 -0.51124948  2.06684311 1.8412119 -3.2422262  5 2.937632
22 -1.81733270  0.39151351 2.3633574 -1.9842250  6 4.338607
23  0.15796272  2.67480250 1.2354482 -3.5726051  7 3.364642
24 -1.01374920 -0.25462480 3.8609433         NA  8 4.283549
25 -0.44018069  2.37976683 1.4271300  0.1525669  1 4.913104
26 -1.44947801  0.66636568 0.9761347 -1.8593690  2 3.550289
27  1.38276645  0.89361910 4.0715730  1.7583258  3 5.488644
28  1.35405816  1.75688151 0.2458366  0.1334016  4 3.965868
29 -0.66293965  2.92249839 2.6279181 -0.6117954  5 4.233280
30 -2.03553017  0.70981564 2.7571505 -2.2757803  6 4.967419
31  0.04128058  3.03548381 1.7200150 -0.5270822  7 4.219825
32 -0.88441054 -0.10520521 3.3802964         NA  8 5.035887
33 -0.44521788  2.76866305 1.5585535  3.6581928  1 5.493640
34 -1.31832214  1.63630484 1.5313611  4.2995457  2 5.825309
35  1.37877834  1.10426733 4.1169144  4.0721678  3 5.736310
36  1.52497108  2.22968000 0.5971719  3.9570720  4 5.859690
37 -0.34953554  3.14765062 3.1941059 -0.2430971  5 5.226246
38 -1.86254928  1.30876499 2.9818303  2.2261139  6 5.634025
39 -0.05425723  3.19698130 2.3805538 -0.1133163  7 5.056034
40 -0.88113708  0.40796654 3.4519979  2.3698206  8 5.932355

mydata <- matrix(data=temp, nrow=40, ncol=7,byrow=T) 
 colnames(mydata) <- tempnames


# Change the file name to the location of your file.  
# You can load the file directly, rather than reading in the data above.  
load("C:\\rob\\mystuff\\RMbook\\rmnotes\\plots\\graphics\\c2ppexampledata_rsave")

dm<- data.frame(y1=mydata[,2], y2=mydata[,3], y3=mydata[,4], y4=mydata[,5], id=mydata[,6], time1=mydata[,7])


# Yet another version of the profile plot function.  

spagplot <- function(time, id, y, xlab="Time",ylab="y",
                     lx=min(time),ux=max(time),
                     ly=min(y,na.rm=T),uy=max(y,na.rm=T),...){
  okidx <- !is.na(y)
  y <- y[okidx]
  id <- id[okidx]
  time <- time[okidx]
  plot(c(lx,ux),c(ly,uy),type="n",xlab="",ylab="")
  title(xlab=xlab,ylab=ylab,...)
  n <- length(time)
  subidx <- unique(id)
  for (i in subidx){
    xtmp <- time[id==i]
    ytmp <- y[id==i]
    points(xtmp,ytmp)
    lines(xtmp,ytmp)
  }
}

dev.off()
x11(width = 4.5, height = 4.5, pointsize = 10,)


#postscript(file="C:\\rob\\mystuff\\RMbook\\rmnotes\\plots\\graphics\\graphics5.ps",
#onefile=F,print=F,height=4.5,width=4.5,horizontal=FALSE,pointsize=10)



par(mfrow=c(2,2))
par(mar=c(4.1,3.1,1.1,1.1))
spagplot(dm$time1,dm$id,dm$y1, ylab="", xlab="(a)")
mtext("Time", side=1, line=2.0, cex=.8)
mtext("Response", side=2, line=2.2, cex=.8)

spagplot(dm$time1, dm$id, dm$y2,ylab="",xlab="(b)")
mtext("Time", side=1, line=2.0, cex=.8)
mtext("Response", side=2, line=2.2, cex=.8)

spagplot(dm$time1, dm$id, dm$y3, ylab="", xlab="(c)")
mtext("Response", side=2, line=2.2, cex=.8)
mtext("Time", side=1, line=2.0, cex=.8)

spagplot(dm$time1, dm$id, dm$y4,ylab="",xlab="(d)")
mtext("Time", side=1, line=2.0, cex=.8)
mtext("Response", side=2, line=2.2, cex=.8)


#dev.off()


## Here's how the data was constructed.  

## First version of same code for this figure.  
## Note: this code has not been cleaned up.  No guarantees at this time, sorry.  

#fig 2.5(a) -- random intercept, times usable for all plots.  

id    <- rep(1:8, 5)
x1    <- matrix(data=(rep(4 * runif(29)[0:7 * 4 + 1] - 2 ,5) + .15 * rnorm(40) ), nrow=8, ncol=5)
times <- matrix(data=NA, nrow=8,ncol=5)
for (i in 1:8) { times[i,] <- 1 + 5 * sort(runif(9))[c(1,3,5,7,9)] }
y1    <- as.vector(x1)
time1 <- as.vector(times)


# fig 2.5(b) -- random intercept, fixed slope

x2    <- matrix(data=(rep(4 * runif(29)[0:7 * 4 + 1] - 2 ,5) + .15 * rnorm(40) ), nrow=8, ncol=5) 
x2    <- x2 + times * .4 - .2 
y2    <- as.vector(x2) 


# fig 2.5(c) -- random slope + fixed slope, fixed intercept, high outlier

par(mar=c(2.1,2.1,1.1,1.1))

temp2 <- (runif(29)[ 0:7 * 4 + 1 ] -.05 ) * 3
temp3 <- rep( temp2 ,5 )
temp4 <- ((time1-min(time1)) * temp3*5 + .5*temp3 + 2 * rnorm(40)) / 16

x3    <- matrix(data=(temp4 - min(temp4)+ .02), nrow=8, ncol=5) 
x3[8,]<- .25 * rnorm(5) + 3.5  
y3    <- as.vector(x3) + min(x3)


# fig 2.5(d) -- fixed quadratic 

x4    <- matrix(data=( (time1 - median(time1))^2 + rnorm(40) ), nrow=8, ncol=5) 
x4[8,]<- c(  (times[8,1] - median(time1))^2 + rnorm(1), NA, NA, NA,(times[8,5] - median(time1))^2 + rnorm(1))
y4    <- as.vector(x4) + min(x4,na.rm=T)




########################
##
## Figure 2.6 Big Mice.  
## Marginal Means and Marginal Standard Deviations against Day.  
## 
#########################

dev.off()
x11(width = 4.5, height = 2.5, pointsize = 10)

#postscript(file="C:/rob/mystuff/RMbook/rmnotes/plots/graphics/graphics6.ps",
#    onefile=FALSE, horizontal=FALSE, height=2.5, width=4.5, pointsize=10)

micemeans <- by(mice$wt, mice$time,mean)
par(mfrow=c(1,2))
par(mar=c(4.1,4.1,2.6,1.1))
plot(0:20, micemeans, xlab="", ylab="", axes=F)
axis(side=1, at=0:4*5, lab=0:4*5, cex=.8, font=10)
axis(side=2, at=(0:5*200+200), lab=(0:5*200+200), cex=.8, font=10)
mtext("Days", side=1, line=1.7, cex=.8)
mtext("(a)", side=1, line=2.5, cex=.8)
mtext("Mean (mg)", side=2, line=2.2, cex=.8)
box()

micesds <- sqrt(by(mice$wt, mice$time,var))
plot(0:20, micesds, xlab="", ylab="", axes=F)
axis(side=1, at=0:4*5, lab=0:4*5, cex=.8, font=10)
axis(side=2, at=(0:6*20+20), lab=(0:6*20+20), cex=.8, font=10)
mtext("Days", side=1, line=1.7, cex=.8)
mtext("(b)", side=1, line=2.5, cex=.8)
mtext("Sd (mg)", side=2, line=2.2, cex=.8)
box()

#dev.off()





############################################
####
#### Profile plots of Pain data original scale and logged scale  
#### Figure 2.7.  
####
############################################



# read in data from wide table.  

#pain1 <- read.table("C:/rob/RMbook/rmnotes/comps/pain/paindatawithheading.txt",header=T)
pain1 <- read.table("paindatawithheading.txt",header=T)


#create long format data set.  

pain2 <- matrix(0,64*4,5)
pain2[4*pain1$id-3,1] <- pain1$id
pain2[4*pain1$id-2,1] <- pain1$id
pain2[4*pain1$id-1,1] <- pain1$id
pain2[4*pain1$id-0,1] <- pain1$id
pain2[4*pain1$id-3,2] <- pain1$tmt
pain2[4*pain1$id-2,2] <- pain1$tmt
pain2[4*pain1$id-1,2] <- pain1$tmt
pain2[4*pain1$id-0,2] <- pain1$tmt
pain2[4*pain1$id-3,3] <- pain1$cs
pain2[4*pain1$id-2,3] <- pain1$cs
pain2[4*pain1$id-1,3] <- pain1$cs
pain2[4*pain1$id-0,3] <- pain1$cs
pain2[4*pain1$id-3,4] <- pain1$t1
pain2[4*pain1$id-2,4] <- pain1$t2
pain2[4*pain1$id-1,4] <- pain1$t3
pain2[4*pain1$id-0,4] <- pain1$t4
pain2[,5]       <- 0
pain2[4*pain1$id,5]   <- pain1$tmt

#Create data frame version of data set.  

pain <-data.frame(pain2)
names(pain) <- c("id", "tmt", "cs", "paintol", "tmt4")
pain$cs <- as.factor(pain$cs)
pain$tmt <- as.factor(pain$tmt)
pain$time <- rep(c(1,2,3,4),64)
pain$logpaintol <- log2(pain$paintol)
pain$cslabels <-factor(pain$cs, labels=c("A","D"))
pain$tmtlabels <-factor(pain$tmt, labels=c("A","D","N"))
pain$tmt4labels <-factor(pain$tmt4, labels=c("0","A","D","N"))


## There may well be extraneous data management going on here;
## code was set up to serve many purposes.  


#YASPFOT (Yet Another Spaghetti Plot Function Or Two)

spagplot26 <- function(time, id, y, xlab="Time",ylab="y", highoffset=0,
                     lx=min(time),ux=max(time),
                     ly=min(y,na.rm=T),uy=max(y,na.rm=T),...){
  okidx <- !is.na(y)
  y <- y[okidx]
  id <- id[okidx]
  time <- time[okidx]
  plot(c(lx,ux),c(ly,uy+highoffset),type="n",xlab="",ylab="",axes=F)
#  title(xlab=xlab,ylab=ylab,...)
  n <- length(time)
  subidx <- unique(id)
  for (i in subidx){
    xtmp <- time[id==i]
    ytmp <- y[id==i]
#    points(xtmp,ytmp)
    lines(xtmp,ytmp)
  }
}
spagplot27 <- function(time, id, y, xlab="Time",ylab="y", highoffset=0,
                     lx=min(time),ux=max(time),
                     ly=min(y,na.rm=T),uy=max(y,na.rm=T),...){
  okidx <- !is.na(y)
  y <- y[okidx]
  id <- id[okidx]
  time <- time[okidx]
  plot(c(lx,ux),c(ly,uy+highoffset),type="n",xlab="",ylab="",axes=F,log="y")
#  title(xlab=xlab,ylab=ylab,...)
  n <- length(time)
  subidx <- unique(id)
  for (i in subidx){
    xtmp <- time[id==i]
    ytmp <- y[id==i]
#    points(xtmp,ytmp,log="y")
    lines(xtmp,ytmp)
  }
}


dev.off()
x11(width = 4.5, height = 4.0, pointsize = 10)


#postscript(file="C:\\rob\\mystuff\\RMbook\\rmnotes\\plots\\graphics\\graphics8.ps",
#    onefile=FALSE, horizontal=FALSE, height=4.0, width=4.5, pointsize=10)

par(mfrow=c(1,2))
par(mar=c(4.1,4.1,2.6,1.1))



spagplot26(pain$time,pain$id,pain$paintol,xlab="Trial Number", ylab="Tolerance (sec)",highoffset=10)
axis(side=1,labels=1:4,at=1:4)
axis(side=2,labels=0:5*50,at=0:5*50)
mtext("Trial", side=1, line=1.7, cex=.8)
mtext("(a)", side=1, line=2.6, cex=.8)
mtext("Tolerance (sec)", side=2, line=2.2, cex=.8)
box()

spagplot27(pain$time,pain$id,pain$paintol,xlab="Trial Number", ylab="Tolerance (sec)")
axis(side=1,labels=1:4,at=1:4)
axis(side=2,labels=c(7.5,15,30,60,120,240),at=c(7.5,15,30,60,120,240))
mtext("Trial", side=1, line=1.7, cex=.8)
mtext("(b)", side=1, line=2.6, cex=.8)
mtext("Tolerance (sec)", side=2, line=2.2, cex=.8)
box()

#dev.off()




############################################
####
#### Plot illustrating different within-subject variability  
#### Figure 2.8
####
############################################


##
## demonstrate different within subject variability.  
##

load("C:\\rob\\mystuff\\RMbook\\rmnotes\\plots\\graphics\\c2withinsubjectvar_rsave")

dev.off()
x11(width = 4.0, height = 3.0, pointsize = 10)


#postscript(file="C:\\rob\\mystuff\\RMbook\\rmnotes\\plots\\graphics\\graphics9.ps",
#    onefile=FALSE, horizontal=FALSE, height=3.0, width=4.0, pointsize=10)

spagplot26(time1, id, y1, ylab="")
axis(side=1, labels=0:6, at=0:6)
axis(side=2, labels=c(-1:4)*5, at=c(-1:4)*5)
mtext("Time", side=1, line=1.7, cex=.8)
mtext("Response", side=2, line=2.2, cex=.8)
box()


#dev.off()


# Create a version of the data in Figure 2.8.  
# I don't think this is the actual original parameter specification but a variant.  

id    <- rep(1:9, 15)
x1    <- matrix(data=(rep(c(2,2.5,3,5,6,7,10,11.5,13) ,15)), nrow=9, ncol=15)
rnoise <- matrix(data=c(1,1,1,2.5,2.5,2.5,6,6,6),nrow=9,ncol=15) * matrix(data= rnorm(90)*.15,nrow=9,ncol=15)
times <- matrix(data=rep(1:15,9), nrow=9,ncol=15,byrow=T)/3
y1    <- as.vector(x1 + rnoise)
time1 <- as.vector(times)

# Uncomment this to save the data.  
#save(y1,time1,id,
#file="C:\\rob\\mystuff\\RMbook\\rmnotes\\plots\\graphics\\c2withinsubjectvar_rsave")




############################################
####
#### plot mean versus variance for the pain data original and log scale
#### Figure 2.9
####
############################################

## already have pain1

painmatraw <- pain1[,4:7]


fullcases <-1:64
fullcases <- fullcases[!is.na(painmatraw[,1]) & !is.na(painmatraw[,2]) & !is.na(painmatraw[,3]) 
  & !is.na(painmatraw[,4]) & !(painmatraw[,1] > 239) & !(painmatraw[,2] > 239)
 & !(painmatraw[,3] > 239) & !(painmatraw[,4] > 239)]

painmeansraw <- by(pain$paintol, pain$id,mean,na.rm=T)
painsdsraw <- sqrt(by(pain$paintol, pain$id,var,na.rm=T))
painmeanslog <- by(log2(pain$paintol), pain$id,mean,na.rm=T)
painsdslog <- sqrt(by(log2(pain$paintol), pain$id,var,na.rm=T))


dev.off()
x11(width = 4.5, height = 2.5, pointsize = 10)


#postscript(file="C:\\rob\\mystuff\\RMbook\\rmnotes\\plots\\graphics\\graphics10.ps",
#    onefile=FALSE, horizontal=FALSE, height=2.5, width=4.5, pointsize=10)

par(mfrow=c(1,2))
par(mar=c(4.1,4.1,1.6,1.1))

plot(painmeansraw[fullcases],painsdsraw[fullcases],xlab="", ylab="", axes=F)
abline(lm(painsdsraw[fullcases]~painmeansraw[fullcases]))
axis(side=1, labels=1:8*10, at=1:8*10, font=10)
axis(side=2, labels=c(0:4)*20, at=c(0:4)*20, font=10)
mtext("Mean", side=1, line=1.7, cex=.8)
mtext("(a)", side=1, line=2.6, cex=.8)
mtext("Sd", side=2, line=2.2, cex=.8)
box()

plot(painmeanslog[fullcases],painsdslog[fullcases],xlab="", ylab="", axes=F)
abline(lm(painsdslog[fullcases]~painmeanslog[fullcases]))
axis(side=1, labels=3:6, at=3:6, font=10)
axis(side=2, labels=0:4/2, at=0:4/2, font=10)
mtext("Mean", side=1, line=1.7, cex=.8)
mtext("(b)", side=1, line=2.6, cex=.8)
mtext("Sd", side=2, line=2.2, cex=.8)
box()

#dev.off()


############################################
####
#### plot profiles for the pain data separately by attender distracter
#### Figure 2.10
####
############################################

painatt <- pain[pain$cs==1,]
paindis <- pain[pain$cs==2,]


dev.off()
x11(width = 4.5, height = 3.0, pointsize = 10)

#postscript(file="C:\\rob\\mystuff\\RMbook\\rmnotes\\plots\\graphics\\graphics11.ps",
#    onefile=FALSE, horizontal=FALSE, height=3.0, width=4.5, pointsize=10)

par(mfrow=c(1,2))
par(mar=c(4.1,4.1,1.6,1.1))


spagplot27(painatt$time, painatt$id, painatt$paintol, xlab="Trial Number", 
        ylab="Tolerance (sec)", ly=min(pain$paintol,na.rm=T))
axis(side=1, labels=1:4, at=1:4, font=10 )
axis(side=2, labels=c(7.5,15,30,60,120,240), at=c(7.5,15,30,60,120,240), font=10)
mtext("Trial", side=1, line=1.7, cex=.8)
mtext("(a)", side=1, line=2.6, cex=.8)
mtext("Tolerance (sec)", side=2, line=2.2, cex=.8)
box()

spagplot27(paindis$time,paindis$id,paindis$paintol, xlab="Trial Number", 
        ylab="Tolerance (sec)", ly=min(pain$paintol,na.rm=T))
axis(side=1, labels=1:4, at=1:4, font=10 )
axis(side=2, labels=c(7.5,15,30,60,120,240), at=c(7.5,15,30,60,120,240), font=10)
mtext("Trial", side=1, line=1.7, cex=.8)
mtext("(b)", side=1, line=2.6, cex=.8)
mtext("Tolerance (sec)", side=2, line=2.2, cex=.8)
box()


#dev.off()


############################################
############################################
####
#### Ozone data
#### Figure 2.11
####
############################################
############################################


ozonesites <- read.table("C:\\rob\\mystuff\\RMbook\\rmnotes\\comps\\ozone\\ozonestations2.txt",header=T)

sites <- c("RESE",
"BURK",
"SIMI",
"PASA",
"WSLA",
"HAWT",
"LGBH",
"ANAH",
"CELA",
"AZUS",
"POMA",
"FONT",
"SNBO",
"RIVR",
"UPLA",
"CLAR",
"LYNN",
"LAHB",
"WHIT",
"PICO")

dev.off()
x11(width = 3.5, height = 4.0, pointsize = 10)

#postscript(file="C:\\rob\\mystuff\\RMbook\\rmnotes\\plots\\graphics\\graphics12.ps",
#    onefile=FALSE, horizontal=FALSE, height=4.0, width=3.5, pointsize=10)


mm <- c(min(ozonesites$long)-.1, max(ozonesites$long)+.1)
nn <- c(min(ozonesites$lat), max(ozonesites$lat))
plot(-mm, nn, type="n", axes=F, xlab="", ylab="")
text(-ozonesites$long, ozonesites$lat, sites, cex=.7)
axis(1, at = (-118.6 + 0:3 * .4), labels =-(-118.6 + 0:3 * .4) )
axis(2, at = (  33.9 + 0:3 * .1), labels = (  33.9 + 0:3 * .1) )
mtext("Longitude", side=1, line=1.7, cex=.8)
mtext("Latitude", side=2, line=2.2, cex=.8)
box()

#dev.off()


############################################
############################################
####
#### Ozone data
#### Figures 2.12, 2.13, 2.14.
####
############################################
############################################

# Need to load "lattice" package first.  In R, go to menu item "packages" then
# to "load package."  Then select lattice from the list.  

ozone <- read.table("C:\\rob\\mystuff\\RMbook\\rmnotes\\comps\\ozone\\basicozonelong.txt", header=T, col.names=c("hour","o3","site","day"))
ozone$id <- rep(1:60,rep(12,60))
ozone$fday <- factor(ozone$day)
ozone$id <- factor(ozone$id)

#
# Orders the sites by the maximum value for better display.  
#

ozone.max.site <- tapply(ozone$o3,ozone$site,max)
ozone$site2 <- ordered(ozone$site,levels = names(sort(ozone.max.site)))
trellis.par.set("strip.background",rep("white",7))
trellis.par.set("background","white")
trellis.par.get("strip.background")

panel.special <- function(x,y)
    {
     temp.a <- length(x)/12
     for( temp.i in 0:(temp.a-1) ) {temp.begin <- temp.i * 12 + 1
          temp.end <- temp.i * 12 + 12
          llines(x[temp.begin:temp.end], y[temp.begin:temp.end], col="black")}
}


foo <- xyplot(formula=ozone$o3 ~ ozone$hour|ozone$site2,type="l",layout=c(5,4,1)
,col="black",panel=panel.special)
foo <- update(foo, main="",cex=1.5)
foo <- update(foo, xlab="Hour", ylab="Ozone",cex=1.3)
#trellis.device(postscript,onefile=FALSE,print.it=FALSE,
#file="C:\\rob\\mystuff\\RMbook\\rmnotes\\plots\\graphics\\graphics14.ps", 
#horizontal=FALSE, height=4.0, width=4.5, pointsize=7)
dev.off()
trellis.device(x11, height=5.5, width=4, pointsize=7)
trellis.par.set("strip.background",rep("white",7))
trellis.par.set("background","white")
trellis.par.set("strip.background",rep("white",7))
foo 


foo1 <- xyplot(formula=ozone$o3 ~ ozone$hour|factor(ozone$day), type="l", 
     layout=c(3,1,1), col="black", panel=panel.special)
foo1 <- update(foo1, main="", cex=1.5)
foo1 <- update(foo1, xlab="Hour", ylab="Ozone",cex=1.3)
#trellis.device(postscript, onefile=FALSE, print.it=FALSE,
#file="C:\\rob\\mystuff\\RMbook\\rmnotes\\plots\\graphics\\graphics13.ps", 
#horizontal=FALSE, height=3.0, width=4.5, pointsize=7)
trellis.par.set("strip.background",rep("white",7))
trellis.par.set("background","white")
trellis.par.set("strip.background",rep("white",7))
foo1 

#dev.off()


#ozone$site2 * ozone$fday 

foo2 <- xyplot(formula=ozone$o3 ~ ozone$hour|ozone$fday*ozone$site2, type="l", 
      layout=c(3,5,4), col="black", panel=panel.special)
foo2 <- update(foo2, main="", cex=1.5)
foo2 <- update(foo2, xlab="Hour", ylab="Ozone",cex=1.3)
#trellis.device(postscript, onefile=FALSE, print.it=FALSE,
#file="C:\\rob\\mystuff\\RMbook\\rmnotes\\plots\\graphics\\graphics13.ps", 
#horizontal=FALSE, height=3.0, width=4.5, pointsize=7)
trellis.par.set("strip.background",rep("white",7))
trellis.par.set("background","white")
trellis.par.set("strip.background",rep("white",7))
foo2

#dev.off()


############################################
############################################
####
#### Weightloss data
#### Figure 2.15
####
############################################
############################################

weightloss<-
read.table("C:\\rob\\mystuff\\RMbook\\rmnotes\\comps\\weightloss\\weightlosslong.txt", 
header=T)


dev.off()
x11(width = 4.5, height = 4.5, pointsize = 10)


#postscript(file="C:\\rob\\mystuff\\RMbook\\rmnotes\\plots\\graphics\\graphics16.ps",
#    onefile=FALSE, horizontal=FALSE, height=4.5, width=4.5, pointsize=10)

spagplot26(weightloss$week, weightloss$id, weightloss$weight, xlab="", 
        ylab="")
axis(side=1, labels=1:8, at=1:8, font=10)
axis(side=2, labels=0:6*20+140, at=0:6*20+140, font=10)
mtext("Week", side=1, line=2.0, cex=.8)
#mtext("(a)", side=1, line=2.6, cex=.8)
mtext("Weight (lb)", side=2, line=2.2, cex=.8)
box()

#dev.off()


# The old way to make figure 2.16 is to manipulate the figure in the 
# r graphics window to be long and skinny.  New way, is this.  
# New Figure 2.16, weightloss data split among 2 different long thin plots.  

spagplot26a <- function(time, id, y, xlab="Time",ylab="y", highoffset=0,
                     lx=min(time),ux=max(time),
                     ly=min(y,na.rm=T),uy=max(y,na.rm=T),...){
  okidx <- !is.na(y)
  y <- y[okidx]
  id <- id[okidx]
  time <- time[okidx]
  plot(c(lx,ux),c(141,209),type="n",xlab="",ylab="",axes=F)
#  title(xlab=xlab,ylab=ylab,...)
  n <- length(time)
  subidx <- unique(id)
  for (i in subidx){
    xtmp <- time[id==i]
    ytmp <- y[id==i]
#    points(xtmp,ytmp)
    lines(xtmp,ytmp)
  }
}
spagplot26b <- function(time, id, y, xlab="Time",ylab="y", highoffset=0,
                     lx=min(time),ux=max(time),
                     ly=min(y,na.rm=T),uy=max(y,na.rm=T),...){
  okidx <- !is.na(y)
  y <- y[okidx]
  id <- id[okidx]
  time <- time[okidx]
  plot(c(lx,ux),c(201,269),type="n",xlab="",ylab="",axes=F)
#  title(xlab=xlab,ylab=ylab,...)
  n <- length(time)
  subidx <- unique(id)
  for (i in subidx){
    xtmp <- time[id==i]
    ytmp <- y[id==i]
#    points(xtmp,ytmp)
    lines(xtmp,ytmp)
  }
}

dev.off()
x11(width = 4.5, height = 6.5, pointsize = 10)


#postscript(file="C:\\rob\\mystuff\\RMbook\\rmnotes\\plots\\graphics\\graphics17.ps",
#    onefile=FALSE, horizontal=FALSE, height=6.5, width=4.5, pointsize=10)

par(mfrow=c(1,2))

spagplot26a(weightloss$week, weightloss$id, weightloss$weight, xlab="", 
        ylab="")
axis(side=1, labels=1:8, at=1:8, font=10)
axis(side=2, labels=0:6*20+140, at=0:6*20+140, font=10)
mtext("Week", side=1, line=2.0, cex=.8)
#mtext("(a)", side=1, line=2.6, cex=.8)
mtext("Weight (lb)", side=2, line=2.2, cex=.8)
box()

spagplot26b(weightloss$week, weightloss$id, weightloss$weight, xlab="", 
        ylab="")
axis(side=1, labels=1:8, at=1:8, font=10)
axis(side=2, labels=0:6*20+140, at=0:6*20+140, font=10)
mtext("Week", side=1, line=2.0, cex=.8)
#mtext("(a)", side=1, line=2.6, cex=.8)
mtext("Weight (lb)", side=2, line=2.2, cex=.8)
box()

#dev.off()


##################
##
## Third weight loss figure, empirical within-subject residual plot
## Figure 2.17
##
##################

wl.means <- tapply(weightloss$weight,weightloss$id,mean)
weightloss$resids <- weightloss$weight - wl.means[weightloss$id]


dev.off()
x11(width = 4.5, height = 4.5, pointsize = 10)


#postscript(file="C:\\rob\\mystuff\\RMbook\\rmnotes\\plots\\graphics\\graphics18.ps",
#    onefile=FALSE, horizontal=FALSE, height=4.5, width=4.5, pointsize=10)


spagplot26(weightloss$week, weightloss$id, weightloss$resids, 
       xlab="Week", ylab="Weight Residuals (lbs)")
axis(side=1, labels=1:8, at=1:8, font=10)
axis(side=2, labels=-2:2*5, at=-2:2*5, font=10)
mtext("Week", side=1, line=2.0, cex=.8)
#mtext("(a)", side=1, line=2.6, cex=.8)
mtext("Weight (lb)", side=2, line=2.2, cex=.8)
box()


#dev.off()



############################################
############################################
####
#### Mice data empirical population residuals 
#### Figure 2.18
####
############################################
############################################


mice$time2 <- factor(mice$time)
mice.means <- tapply(mice$wt,mice$time2,mean)
mice$upopresids <- mice$wt - mice.means[mice$time+1]

dev.off()
x11(width = 4.0, height = 3.5, pointsize = 10)


#postscript(file="C:/rob/mystuff/RMbook/rmnotes/plots/graphics/graphics19.ps",
#onefile=FALSE, horizontal=F, height=4.0, width=3.5, pointsize=10)

spagplot26(mice$time,mice$id,mice$upopresids,xlab="", ylab="")
#axis(side=1, labels=0:20, at=0:20, font=10)
axis(side=1, at=0:4*5, lab=0:4*5, cex=.8, font=10)
axis(side=2, labels=-3:3*100, at=-3:3*100, font=10)
mtext("Day", side=1, line=2.0, cex=.8)
#mtext("(a)", side=1, line=2.6, cex=.8)
mtext("Weight (mg)", side=2, line=2.2, cex=.8)
box()

#dev.off()

############################################
############################################
####
#### Correlation matrix computations
####
############################################
############################################
## 
## Create Table 2.1.
##
## Pain data correlations and standard errors of the correlations in these
## next two computations.  
temp2 <- cor(painmat,use="pairwise.complete.obs")
temp3 <- matrix(data=c(1,62,58,60,62,1,58,61,58,58,1,58,60,61,58,1),nrow=4,ncol=4)
(1-temp2^2)/sqrt(temp3-3)
# These are temp3 results, the standard errors of the correlation in temp2.  
#           [,1]       [,2]       [,3]       [,4]
#[1,]        NaN 0.06156724 0.04052349 0.08408203
#[2,] 0.06156724        NaN 0.06539382 0.07378227
#[3,] 0.04052349 0.06539382        NaN 0.05608517
#[4,] 0.08408203 0.07378227 0.05608517        NaN


############################################
#### 
####  Calculate Table 2.3, Small Mice correlations.  
####  temp2 are the correlations.  The following calculation are 
####  the estimated, rounded standard errors of the correlations.  
#### 
############################################


micetemp  <- mice[ mice$time==2 | mice$time==5 | mice$time==8 | mice$time==11 | mice$time==14 | mice$time==17 | mice$time==20,]$wt
smallmice <- matrix(data=micetemp,ncol=7,byrow=T)
temp2 <- cor(smallmice, use="pairwise.complete.obs")
round((1-temp2^2)/(sqrt(11)),2)

############################################
############################################
####
#### Pain Tolerance scatterplot matrix
#### For Figure 2.19 
####
############################################
############################################

dev.off()
x11(width = 4.5, height = 4.5, pointsize = 10)

painmat <- log2 (pain1[,4:7])

#postscript(file="C:/rob/mystuff/RMbook/rmnotes/plots/graphics/graphics20.ps",
#onefile=FALSE, horizontal=F, height=4.5, width=4.5, pointsize=10)

pairs(painmat, labels=c("T1", "T2", "T3", "T4"), gap=.4,
row1attop=T,cex.labels=1.3)
mtext("Pain tolerance (log_2 seconds)", side=1, line=3.8, cex=.8)
mtext("Pain tolerance (log_2 seconds)", side=2, line=3, cex=.8)

#dev.off()

############################################
############################################
####
#### Small Mice data scatterplot matrix
#### For Figure 2.20 
####
############################################
############################################


days<- c("Day 0", "Day 1", "Day 2", "Day 3", "Day 4", "Day 5", "Day 6", 
"Day 7", "Day 8", "Day 9", "Day 10", "Day 11", "Day 12", "Day 13", "Day 14", 
"Day 15", "Day 16", "Day 17", "Day 18", "Day 19", "Day 20")


dev.off()
x11(width = 4.5, height = 4.5, pointsize = 10)

micemat <- matrix(data=scan(),byrow=T,ncol=7)
190 388 621 823 1078 1132 1191
218 393 568 729 839 852 1004
141 260 472 662 760 885 878
211 394 549 700 783 870 925
209 419 645 850 1001 1026 1069
193 362 520 530 641 640 751
201 361 502 530 657 762 888
202 370 498 650 795 858 910
190 350 510 666 819 879 929
219 399 578 699 709 822 953
225 400 545 690 796 825 836
224 381 577 756 869 929 999
187 329 441 525 589 621 796
278 471 606 770 888 1001 1105



#postscript(file="C:/rob/mystuff/RMbook/rmnotes/plots/graphics/graphics21.ps",
#onefile=FALSE, horizontal=F, height=4.5, width=4.5, pointsize=10)

pairs(micemat, gap=.4, row1attop=T, labels=days[0:6*3+3], cex.labels=1.3 )
mtext("Weight (mg)", side=1, line=3.8, cex=.8)
mtext("Weight (mg)", side=2, line=3, cex=.8)

#dev.off()


############################################
############################################
####
#### Constructed profile plot illustrating various correlations.  
#### Figure 2.21 
####
############################################
############################################

#save(x1, x2, y1, y2, time1, time2, id1, id2, texty1, texty2, textx1, textx2, corr1, corr2, 
#file="C:\\rob\\mystuff\\RMbook\\rmnotes\\plots\\graphics\\c2ppcorrdata_rsave")

load(file="C:\\rob\\mystuff\\RMbook\\rmnotes\\plots\\graphics\\c2ppcorrdata_rsave")



#R code remaking the original data.   

x1 <- matrix(data=NA, nrow=20, ncol=8)
x2 <- matrix(data=NA, nrow=20, ncol=8)
x1[,1] <- rnorm(20)
x2[,1] <- rnorm(20)
corr1 <- c( .99,  .95,  .8,  .6,  .4,  .2, 0)
corr2 <- c(-.99, -.95, -.8, -.6, -.4, -.2, 0)
corr1text <- c( ".99",  ".95",  ".8",  ".6",  ".4",  ".2", "0")
corr2text <- c("-.99", "-.95", "-.8", "-.6", "-.4", "-.2", "0")

sd1 <- sqrt(1 - corr1^2)
sd2 <- sqrt(1 - corr2^2)

x1[,2] <- x1[,1] * corr1[1] + rnorm(20)*sd1[1]
x1[,3] <- x1[,2] * corr1[2] + rnorm(20)*sd1[2]
x1[,4] <- x1[,3] * corr1[3] + rnorm(20)*sd1[3]
x1[,5] <- x1[,4] * corr1[4] + rnorm(20)*sd1[4]
x1[,6] <- x1[,5] * corr1[5] + rnorm(20)*sd1[5]
x1[,7] <- x1[,6] * corr1[6] + rnorm(20)*sd1[6]
x1[,8] <- x1[,7] * corr1[7] + rnorm(20)*sd1[7]

x2[,2] <- x2[,1] * corr2[1] + rnorm(20)*sd2[1]
x2[,3] <- x2[,2] * corr2[2] + rnorm(20)*sd2[2]
x2[,4] <- x2[,3] * corr2[3] + rnorm(20)*sd2[3]
x2[,5] <- x2[,4] * corr2[4] + rnorm(20)*sd2[4]
x2[,6] <- x2[,5] * corr2[5] + rnorm(20)*sd2[5]
x2[,7] <- x2[,6] * corr2[6] + rnorm(20)*sd2[6]
x2[,8] <- x2[,7] * corr2[7] + rnorm(20)*sd2[7]

#
# Any time you see repetitious code like that, you
# know there is a more efficient way to do it.  
#

time1 <- c(1:8)
time2 <- c(1:8)


y1 <- as.vector(x1)
y2 <- as.vector(x2)
id1 <- rep(1:20,8)
id2 <- rep(1:20,8)
time1 <- rep(1:8,rep(20,8))
time2 <- rep(1:8,rep(20,8))


#spaghetti plot modified with a lower extent for y for this plot.  
spagplot28 <- function(time, id, y, xlab="Time",ylab="y",lowext=.75,
                     lx=min(time),ux=max(time),
                     ly=min(y-lowext,na.rm=T),uy=max(y,na.rm=T),...){
  okidx <- !is.na(y)
  y <- y[okidx]
  id <- id[okidx]
  time <- time[okidx]
  plot(c(lx,ux), c(ly,uy), type="n", xlab="", ylab="", axes=F)
#  title(xlab=xlab,ylab=ylab,...)
  n <- length(time)
  subidx <- unique(id)
  for (i in subidx){
    xtmp <- time[id==i]
    ytmp <- y[id==i]
    points(xtmp,ytmp)
    lines(xtmp,ytmp)
  }
}

dev.off()
x11(width = 4.5, height = 6.0, pointsize = 10)

#postscript(file="C:\\rob\\mystuff\\RMbook\\rmnotes\\plots\\graphics\\graphics22.ps", 
#onefile=F, print=F, horizontal=FALSE, height=6.0, width=4.5, pointsize=10)

par(mfrow=c(2,1))
par(mar=c(4.1,4.1,1.1,1.1))

spagplot28(time1, id1, y1, ylab="", lowext=.5)
axis(side=1, labels=1:8, at=1:8, cex=.7)
axis(side=2, labels=c("",-2:2), at=-3:2, cex=.7)

texty1 <- c(-2.5,-2.75,-2.5,-2.75,-2.5,-2.75,-2.5) +.25
textx1 <- c(1.5,2.5,3.5,4.5,5.5,6.5,7.5)
text(textx1, texty1, corr1text, cex=1.0)

mtext("Time", side=1, line=1.7, cex=.85)
mtext("(a)", side=1, line=2.6, cex=.85)
mtext("Response", side=2, line=2.2, cex=.85)
box()

# second figure part, negative correlations

spagplot28(time2, id2, y2, ylab="", lowext=.6)
axis(side=1, labels=1:8, at=1:8, cex=.7)
axis(side=2, labels=c("",-2:2), at=-3:2, cex=.7)

texty2 <- c(-2.5,-2.75,-2.5,-2.75,-2.5,-2.75,-2.5) +.25
textx2 <- c(1.5,2.5,3.5,4.5,5.5,6.5,7.5)
text(textx2, texty2, corr2text, cex=1.0)

mtext("Time", side=1, line=1.7, cex=.85)
mtext("(b)", side=1, line=2.6, cex=.85)
mtext("Response", side=2, line=2.2, cex=.85)
box()

#dev.off()

#save(x1, x2, y1, y2, time1, time2, id1, id2, texty1, texty2, textx1, textx2, corr1, corr2, 
#file="C:\\rob\\mystuff\\RMbook\\rmnotes\\plots\\chapter2\\c2ppcorrdata_rsave")


############################################
############################################
####
#### Construct correlogram for ozone data.  
#### Figure 2.22 
####
############################################
############################################

# convert to wide format.  
ozonemat <- matrix(data=ozone$o3,byrow=T,ncol=12)

co <- cor(ozonemat)
tdiff <- matrix(data=NA,nrow=12,ncol=12)
for (i in 1:12) {
for (j in 1:12) { tdiff[i,j] <- (i-j)}}
ocols <- matrix(data=rep(1:12,12), nrow=12,ncol=12,byrow=T)
orows <- matrix(data=rep(1:12,12), nrow=12,ncol=12)

ocorrs <- data.frame(cbind(corrs=co[1:144],diff=tdiff[1:144],col=ocols[1:144], 
           row=orows[1:144]))

ocorr <- ocorrs[ocorrs$diff>0,]

# One version of the correlogram... but not the main plot.  
#plot(ocorr$diff, ocorr$corrs, xlab="time difference", 
#         ylab="correlation", main="Ozone Correlogram", 
#         pch=as.character(ocorr$col))

#postscript(file="C:/rob/mystuff/RMbook/rmnotes/comps/ozone/ctwoozonecorrelogram.eps",
#onefile=FALSE, horizontal=FALSE, height=4.5, width=4.5, pointsize=10)

dev.off()
x11(width = 4.5, height = 4.5, pointsize = 10)

#postscript(file="C:\\rob\\mystuff\\RMbook\\rmnotes\\plots\\graphics\\graphics23.ps", 
#onefile=F, print=F, horizontal=FALSE, height=4.5, width=4.5, pointsize=10)

plot(ocorr$diff, ocorr$corrs, xlab="", axes=F,
         ylab="", main="", type="n")
for(j in 1:12) {points(ocorr$diff[ocorr$col==j], ocorr$corrs[ocorr$col==j]) }
for(j in 1:12) {lines(ocorr$diff[ocorr$col==j], ocorr$corrs[ocorr$col==j]) }
axis(side=1, labels=1:5*2, at=1:5*2, cex=.7)
axis(side=2, labels=c("-.2","0",".2",".4",".6",".8","1"), at=-1:5*.2, cex=.7)
mtext("Lag (hours)", side=1, line=1.7, cex=.85)
mtext("Correlation", side=2, line=2.2, cex=.85)
box()

#dev.off()

############################################
############################################
####
#### Empirical fit and predictions for the Mice data
#### Figure 2.23 
####
############################################
############################################


mm  <- tapply(mice$wt, mice$time2, mean, na.rm=TRUE)
msd <- sqrt(tapply(mice$wt, mice$time2, var, na.rm=TRUE))
mn  <- tapply(rep(1,length(mice$time2)), mice$time2, sum)
upperpred <- mm + 1.96*msd
lowerpred <- mm - 1.96*msd
upperfit <- mm + 1.96*msd/sqrt(mn)
lowerfit <- mm - 1.96*msd/sqrt(mn)


dev.off()
x11(width = 4.5, height = 2.5, pointsize = 10)

#postscript(file="C:\\rob\\mystuff\\RMbook\\rmnotes\\plots\\graphics\\graphics24.ps", 
#onefile=F, print=F, horizontal=FALSE, height=2.5, width=4.5, pointsize=10)

par(mfrow=c(1,2))
par(mar=c(4.1,3.1,2.1,1.1))
plot(rep(0:20,2),c(upperpred,lowerpred),type="n",xlab="",ylab="",
          main="")
lines(0:20,mm)
arrows(0:20,lowerfit,0:20,upperfit,angle=90,length=.05, code=3)
axis(side=1, labels=0:4*5, at=0:4*5)
axis(side=2, labels=1:5*200, at=1:5*200)
mtext("Day", side=1, line=1.7, cex=.8)
mtext("(a)", side=1, line=2.6, cex=.8)
mtext("Weight (mg)", side=2, line=2.2, cex=.8)
box()


plot(rep(0:20,2),c(upperpred,lowerpred),type="n",xlab="",ylab="",
          main="")
lines(0:20,mm)
arrows(0:20,lowerpred,0:20,upperpred,angle=90,length=.05, code=3)
axis(side=1, labels=0:4*5, at=0:4*5)
axis(side=2, labels=1:5*200, at=1:5*200)
mtext("Day", side=1, line=1.7, cex=.8)
mtext("(b)", side=1, line=2.6, cex=.8)
mtext("Weight (mg)", side=2, line=2.2, cex=.8)
box()

#dev.off()


############################################
############################################
####
#### Empirical fit and predictions for the pain data separately for attenders
#### and distracters.  
#### Figure 2.24
####
############################################
############################################


painms <- tapply(pain$logpaintol,list(pain$cslabels,pain$time),mean,na.rm = TRUE)
painsds <- sqrt(tapply(pain$logpaintol,list(pain$cslabels,pain$time),var,na.rm = TRUE))
painns <- table(pain$cs[!is.na(pain$paintol)],pain$time[!is.na(pain$paintol)])

painlower <- painms - 1.96 * painsds/sqrt(painns) 
painupper <- painms + 1.96 * painsds/sqrt(painns) 

epainlower <- 2^painlower
epainupper <- 2^painupper
epainms <- 2^painms

dev.off()
X11(width=4.5, height=2.5, pointsize=10)

#postscript(file="C:/rob/mystuff/RMbook/rmnotes/plots/graphics/graphics25.ps",
#onefile=FALSE, horizontal=FALSE, height=2.5, width=4.5, pointsize=10)

par(mfrow=c(1,2))
par(mar=c(4.1,3.1,1.1,1.1))

plot(rep(c(1:4-.15,1:4 + .15),2), c(painlower,painupper), type="n", 
    xlab="", ylab="", axes=F)
lines(1:4-.1,painms[1,])
lines(1:4+.1,painms[2,])

arrows(1:4-.1,painlower[1,],1:4-.1,painupper[1,],angle=90,length=.05, code=3)
arrows(1:4+.1,painlower[2,],1:4+.1,painupper[2,],angle=90,length=.05, code=3)

axis(side=1, at=c(1:4), label=c(1:4))
axis(side=2, at=c(4.2,4.6,5.0,5.4), label=c(4.2,4.6,5.0,5.4))
mtext("Trial", side=1, line=1.7, cex=.8)
mtext("(a)", side=1, line=2.6, cex=.8)
mtext("Log_2 seconds", side=2, line=2.2, cex=.8)
box()

plot(rep(c(1:4-.15,1:4 + .15),2), c(epainlower,epainupper), type="n",
      xlab="", ylab="", axes=F)
lines(1:4-.1,epainms[1,])
lines(1:4+.1,epainms[2,])

arrows(1:4-.1,epainlower[1,],1:4-.1,epainupper[1,],angle=90,length=.05, code=3)
arrows(1:4+.1,epainlower[2,],1:4+.1,epainupper[2,],angle=90,length=.05, code=3)

axis(side=1, at=c(1:4),label=c(1:4))
axis(side=2, at=c(20,30,40),label=c(20,30,40))
mtext("Trial", side=1, line=1.7, cex=.8)
mtext("(b)", side=1, line=2.6, cex=.8)
mtext("Seconds", side=2, line=2.2, cex=.8)
box()

#dev.off()



############################################
############################################
####
#### First, edit vagal tone, then produce plots.  
#### Figure 2.27
####
############################################
############################################



timesall <- scan()
1 2  4 5 NA
4 NA NA NA NA
1 2  4 5 NA
2   4 5 NA NA
1 2 3 4 5
1 2 3 4 5
4 5 NA NA NA
1   3 4 5 NA
3  5 NA NA NA
2 3 4 5 NA
3 4 NA NA NA
1 2 3 NA NA
1 2 4 5 NA
2 3 4 5 NA
1 2 3 4 5
1 2 3 NA NA
2 NA NA NA NA

vtall  <- scan()
1.96 3.04  3.18 3.27 NA
0.97 NA NA NA NA
3.93 3.86  3.59 3.27 NA
2.24  2.55 2.27 NA NA
2.57 2.52 3.92 3.22 1.28
4.44 3.07 1.7 3.43 3.43
1.15 1.03 NA NA NA
1.74 1.23 1.49 0.92 NA
2.14 4.92 NA NA NA
3.06 1.94 2.6 1.18 NA
2.97 5.23 NA NA NA
4.4 3.92 0 NA NA
1.96 1.95 2.51 1.25 NA
3.51 1.88 2.14 2.42 NA
3.63 3.11 1.87 4.46 3.7
2.91 2.91 0.61 NA NA
5.03 NA NA NA NA

vtsep  <- scan()
1.96 3.04 NA NA NA
3.18 3.27 NA NA NA
0.97 NA NA NA NA
3.93 3.86 NA NA NA
3.59 3.27 NA NA NA
2.24 NA NA NA NA
2.55 2.27 NA NA NA
2.57 2.52 3.92 3.22 1.28
4.44 3.07 1.7 3.43 3.43
1.15 1.03 NA NA NA
1.74 NA NA NA NA
1.23 1.49 0.92 NA NA
2.14 NA NA NA NA
4.92 NA NA NA NA
3.06 1.94 2.6 1.18 NA
2.97 5.23 NA NA NA
4.4 3.92 0 NA NA
1.96 1.95 NA NA NA
2.51 1.25 NA NA NA
3.51 1.88 2.14 2.42 NA
3.63 3.11 1.87 4.46 3.7
2.91 2.91 0.61 NA NA
5.03 NA NA NA NA

timessep<- scan()
1 2  NA NA NA
4 5 NA NA NA
4 NA NA NA NA
1 2  NA NA NA
4 5 NA NA NA
2  NA NA NA NA
4 5 NA NA NA
1 2 3 4 5
1 2 3 4 5
4 5 NA NA NA
1 NA NA NA NA
3 4 5 NA NA
3 NA NA NA NA
5 NA NA NA NA
2 3 4 5 NA
3 4 NA NA NA
1 2 3 NA NA
1 2 NA NA NA
4 5 NA NA NA
2 3 4 5 NA
1 2 3 4 5
1 2 3 NA NA
2 NA NA NA NA



vagaltone2 <- data.frame(vt=vtall,times=timesall,id=rep(1:17,rep(5,17)))
vagaltone3 <- data.frame(vt=vtsep,times=timessep,id=rep(1:23,rep(5,23)))

spagplot28 <- function(time, id, y, xlab="",ylab="", highoffset=0,
                     lx=min(time),ux=max(time),
                     ly=min(y,na.rm=T),uy=max(y,na.rm=T),...){
  okidx <- !is.na(y)
  y <- y[okidx]
  id <- id[okidx]
  time <- time[okidx]
  plot(c(lx,ux),c(ly,uy+highoffset),type="n",xlab="",ylab="",axes=F)
#  title(xlab=xlab,ylab=ylab,...)
  n <- length(time)
  subidx <- unique(id)
  for (i in subidx){
    xtmp <- time[id==i]
    ytmp <- y[id==i]
    points(xtmp,ytmp)
    lines(xtmp,ytmp)
  }
}

dev.off()
X11(width=4.5, height=5, pointsize=10)

#postscript(file="C:/rob/mystuff/RMbook/rmnotes/plots/graphics/graphics26.ps",
#onefile=FALSE, horizontal=FALSE, height=5, width=4.5, pointsize=10)

par(mfrow=c(2,1))
par(mar=c(4.1,3.1,1.1,1.1))


spagplot28(vagaltone2$times, vagaltone2$id, vagaltone2$vt, ylab="",xlab="")
axis(side=1, labels=1:5, at=1:5)
axis(side=2, labels=0:5, at=0:5)
mtext("Time", side=1, line=1.7, cex=.8)
mtext("(a)", side=1, line=2.6, cex=.8)
mtext("Vagal tone", side=2, line=2.2, cex=.8)
box()

spagplot28(vagaltone3$times, vagaltone3$id, vagaltone3$vt, ylab="",xlab="")
axis(side=1, labels=1:5, at=1:5)
axis(side=2, labels=0:5, at=0:5)
mtext("Time", side=1, line=1.7, cex=.8)
mtext("(b)", side=1, line=2.6, cex=.8)
mtext("Vagal tone", side=2, line=2.2, cex=.8)
box()

#dev.off()



############################################
####
#### Plot constructed data with and without drop out.  
#### Figure 12.1
####
############################################


spagplot25 <- function(time, id, y, xlab="",ylab="",
                     lx=min(time),ux=max(time),
                     ly=min(y,na.rm=T),uy=max(y,na.rm=T),...){
  okidx <- !is.na(y)
  y <- y[okidx]
  id <- id[okidx]
  time <- time[okidx]
  plot(c(lx,ux), c(ly,uy), type="n", xlab="", ylab="", axes=F)
#  title(xlab=xlab,ylab=ylab,...)
  n <- length(time)
  subidx <- unique(id)
  for (i in subidx){
    xtmp <- time[id==i]
    ytmp <- y[id==i]
    points(xtmp,ytmp)
    lines(xtmp,ytmp)
  }
}


# Figure a random intercept with selective drop out.  
# 

#save(y1,y3,time1,id,usethese,
#file="C:\\rob\\mystuff\\RMbook\\rmnotes\\plots\\chapter2\\c2dropoutdata_rsave")
load("C:\\rob\\mystuff\\RMbook\\rmnotes\\plots\\graphics\\c2dropoutdata_rsave")

# first two pictures calculations: random intercept with and without selective 
# drop out.  
id    <- rep(1:10, 5)
x1    <- matrix(data=(rep(4 * sort(runif(41))[0:9 * 4 + 3] - 2 ,5) + .15 * rnorm(50) ),
 nrow=10, ncol=5)
times <- matrix(data=NA, nrow=10,ncol=5)
for (i in 1:10) { times[i,] <- 1 + 5 * sort(runif(21))[c(3,7,11,15,19)] }
y1    <- as.vector(x1) + 2
time1 <- as.vector(times)

usethese <- c(1:10, 13:20,25:30,37:40,49:50)

# second two pictures calculations
# random slope with and without selective drop out.  
temp2 <- (sort(runif(40))[ 0:9 * 4 + 1 ] -.5 ) * 3
temp3 <- rep( temp2 ,5 )
temp4 <- ((time1-min(time1)) * temp3*5 + .5*temp3 + 2 * rnorm(50)) / 16

x3    <- matrix(data=(temp4 - min(temp4)+ .02), nrow=10, ncol=5) 
y3    <- as.vector(x3) + min(x3)



dev.off()
x11(width = 4.5, height = 4.5, pointsize = 10)



#postscript(file="C:\\rob\\mystuff\\RMbook\\rmnotes\\plots\\graphics\\graphics7.ps",
#    onefile=FALSE, horizontal=FALSE, height=4.5, width=4.5, pointsize=10)

par(mfrow=c(2,2))
par(mar=c(4.1,4.1,2.6,1.1))

# random intercept data
spagplot25(time1, id, y1, ylab="")
axis(side=1, labels=1:6, at=1:6)
axis(side=2, labels=c(0:5), at=0:5)
mtext("Time", side=1, line=1.7, cex=.8)
mtext("(a)", side=1, line=2.6, cex=.8)
mtext("Response", side=2, line=2.2, cex=.8)
box()

# same plot missing the data after dropout.  
spagplot25(time1[usethese], id[usethese], y1[usethese], ylab="")
axis(side=1, labels=1:6, at=1:6)
axis(side=2, labels=c(0:5), at=0:5)
mtext("Time", side=1, line=1.7, cex=.8)
mtext("(b)", side=1, line=2.6, cex=.8)
mtext("Response", side=2, line=2.2, cex=.8)
box()

# figure c.  Random slope, complete data.  
spagplot25(time1,id,y3,ylab="")
axis(side=1,labels=1:6,at=1:6)
axis(side=2,labels=c(0:5),at=0:5)
mtext("Time", side=1, line=1.7, cex=.8)
mtext("(c)", side=1, line=2.6, cex=.8)
mtext("Response", side=2, line=2.2, cex=.8)
box()

# figure d with dropout.  

spagplot25(time1[usethese],id[usethese],y3[usethese],ylab="",lx=min(time1),ux=max(time1),
         ly=min(y3,na.rm=T),uy=max(y3,na.rm=T))
axis(side=1,labels=1:6,at=1:6)
axis(side=2,labels=c(0:5),at=0:5)
mtext("Time", side=1, line=1.7, cex=.8)
mtext("(d)", side=1, line=2.6, cex=.8)
mtext("Response", side=2, line=2.2, cex=.8)
box()

#dev.off()