emission<-read.table('A:/cs63.dat',header=T)
attach(emission)
dim(emission) # 52  6
emission

## I. analysis based on averaged measurements

# averaging over repeats

hcmat<-matrix(hc, byrow=T, ncol=4); hcmat # 13 by 4
hcmean<-apply(hcmat,1, mean); hcmean  # vector, 13 by 1

co2mat<-matrix(co2,byrow=T,ncol=4)
co2mean<-apply(co2mat,1,mean)

comat<-matrix(co,byrow=T,ncol=4)
comean<-apply(comat,1,mean)

# table 6.6
daymean<-1:13
devmean<-c(rep(0,2),rep(1,9),rep(0,2))
table66<-data.frame(daymean,devmean,
  round(hcmean,2),round(co2mean,2),round(comean,2))
table66
 
# graphics

boxplot(hcmean[devmean==0])
boxplot(hcmean[devmean==1])

# two boxplots together
boxplot(split(hcmean,devmean),xlab='Device', ylab='mean hc') 

boxplot(split(co2mean,devmean),xlab='Device', ylab='mean co2')

boxplot(split(comean,devmean),xlab='Device', ylab='mean co')

# t-tests: 

c(mean(hcmean[devmean==0]), mean(hcmean[devmean==1]))
c(sd(hcmean[devmean==0]), sd(hcmean[devmean==1]))

help(t.test)
# paired = FALSE, var.equal = FALSE

t.test(hcmean[devmean==0], hcmean[devmean==1])
# not sig.

c(sd(co2mean[devmean==0]), sd(co2mean[devmean==1]))
t.test(co2mean[devmean==0],co2mean[devmean==1],var.equal=T)
# not sig.

c(sd(comean[devmean==0]), sd(comean[devmean==1]))
t.test(comean[devmean==0],comean[devmean==1],var.equal=T)
t.test(comean[devmean==0],comean[devmean==1])
# not sig.

## II. ANOVA approach

fit1<-lm(hc~dev+factor(day)+factor(rep))
summary(fit1)

fit2<-lm(co2~dev+factor(day)+factor(rep))
summary(fit2)

fit3<-lm(co~dev+factor(day)+factor(rep))
summary(fit3)

## III. MANOVA

y<-cbind(hc,co2,co)
fit4<-manova(y~dev+factor(day)+factor(rep))
summary(fit4,test='Hotelling')

## IV. remarks
# 1. This example, case study 6.3 from Cabrera and McDougall, 
# is for illustrative purposes. In fact, this is not a perfect example
# because of the singularity issue.

# 2. For SAS users, MANOVA method may be implemented by using "proc glm" 
# with a "repeated" statement