#=====================================================================
# Figure 7.2.R.
# R script to read blood pressure data and draw a profile plot
# of the means under different treatment combinations.
#=====================================================================

#---------------------------------------------------------------------
# 1.  Read the raw numbers into a preliminary data frame having 6 rows
# and 12 columns.  The raw numbers are assumed to be in a space-delimited
# text file named "blood pressure data.txt" in the current working directory.
#---------------------------------------------------------------------
bp.data=read.table("blood pressure data.txt",header=FALSE)

#---------------------------------------------------------------------
# 2.  Stack the columns of the data frame into a vector named bp.
#---------------------------------------------------------------------
bp=c(bp.data[,1],bp.data[,2],bp.data[,3],bp.data[,4],
     bp.data[,5],bp.data[,6],bp.data[,7],bp.data[,8],
     bp.data[,9],bp.data[,10],bp.data[,11],bp.data[,12])

#---------------------------------------------------------------------
# 3.  Use logical statements to build a text vector named biofeed
# containing "y" and "n" on the appropriate lines.  Also, use logical
# statements to build a text vector named drug containing “a”, “b”,
# and “c” on the appropriate lines.
#---------------------------------------------------------------------
n.bp=length(bp)
biofeed=character(n.bp)
drug=biofeed
for (i in ((0:(n.bp/6-1))*6+1)) {
  biofeed[i:(i+2)]="y"
  biofeed[(i+3):(i+5)]="n"
  drug[c(i,i+3)]="a"
  drug[c(i+1,i+4)]="b"
  drug[c(i+2,i+5)]="c"
}

#---------------------------------------------------------------------
# 4.  Use logical statements to build a text vector named diet
# containing “y” and “n” on the appropriate lines.
#---------------------------------------------------------------------
rnum=1:n.bp
diet=character(n.bp)
diet[rnum<=n.bp/2]="n"
diet[rnum>n.bp/2]="y"

#---------------------------------------------------------------------
# 5.  Combine bp, biofeed, diet, drug, into a data frame, and save the
# data frame as a file in the computer system for future use.
#---------------------------------------------------------------------
bp.data.new=data.frame(bp,diet,biofeed,drug)

#---------------------------------------------------------------------
# 6.  Calculate the mean blood pressures within every treatment
# combination of biofeed, drug, and diet.  Put the means in a new data
# frame named bp.means.  Give the variables names.
#---------------------------------------------------------------------
bp.means=aggregate(bp,by=list(diet,biofeed,drug),FUN=mean)

names(bp.means)=c("diet.m","biofeed.m","drug.m","bp.m")
attach(bp.means)

#---------------------------------------------------------------------
# 7.  Graph the means in a profile plot.
#---------------------------------------------------------------------
plot(c(1,2,3),bp.m[(diet.m=="n")&(biofeed.m=="n")],type="o",
  lty=1,pch=1,cex=1.5,ylim=c(160,210),xlab="drug",ylab="blood pressure",
  xaxt="n")
Axis(at=c(1,2,3),side=1,labels=c("a","b","c"))
points(c(1,2,3),bp.m[(diet.m=="y")&(biofeed.m=="n")],type="o",
  lty=1,pch=0)
points(c(1,2,3),bp.m[(diet.m=="n")&(biofeed.m=="y")],type="o",
  lty=2,pch=1,cex=1.5)
points(c(1,2,3),bp.m[(diet.m=="y")&(biofeed.m=="y")],type="o",
  lty=2,pch=0)