## Question 5.9

## part(a)
Xbar<-c(95.52,164.38,55.69,93.39,17.98,31.13)
Sn<-matrix(c(3266.46,1343.97,731.54,1175.50,162.68, 238.37,
 1343.97,  721.91,  324.25,  537.35,   80.17,  117.73,
  731.54,  324.25,  179.28,  281.17,   39.15,   56.80,
 1175.50,  537.35,  281.17,  474.98,   63.73,   94.85,
  162.68,   80.17,   39.15,   63.73,    9.95,   13.88,
  238.37,  117.73,   56.80,   94.85,   13.88,   21.26), 6, 6)
qchisq6<-qchisq(0.95,6)
n<-61
simultaneCIs<-cbind(Xbar-sqrt(qchisq6)*sqrt(diag(Sn/n)),Xbar+sqrt(qchisq6)*sqrt(diag(Sn/n)))

## part(b)

Sn14<-Sn[c(1,4),c(1,4)]
center<-Xbar[c(1,4)]
x1<-seq(70,120,length=100)
x2<-seq(82,102,length=100)
plot(x1,x2,type="n",asp=1,xlab="weight",ylab="girth")
points(center[1],center[2],col=2,pch=19)
npoints<-1000
n<-61
r<-sqrt(qchisq(0.95,2)/n)
theta<-seq(0, 2*pi, length = npoints)
v<-rbind(r*cos(theta), r*sin(theta))
## transform for points on ellipse
z<-backsolve(chol(solve(Sn14)),v)+center
## plot points
lines(t(z))

## part(c)

adjustalpha<-1-0.05/(2*6)
BonferroniCIs<-cbind(Xbar-qt(adjustalpha,n-1)*sqrt(diag(Sn/n)),Xbar+qt(adjustalpha,n-1)*sqrt(diag(Sn/n)))

## part (d)

Sn14<-Sn[c(1,4),c(1,4)]
center<-Xbar[c(1,4)]
eigenSnvecs<-eigen(Sn14)$vectors
x1<-seq(70,125,length=100)
x2<-seq(80,105,length=100)
plot(x1,x2,type="n",xlab="weight",ylab="girth")
points(center[1],center[2],col=2,pch=19)
npoints<-1000
n<-61
r<-sqrt(qchisq(0.95,2)/n)
theta<-seq(0, 2*pi, length = npoints)
v<-rbind(r*cos(theta), r*sin(theta))
## transform for points on ellipse
z<-backsolve(chol(solve(Sn14)),v)+center
## plot points
lines(t(z))
abline(h=c(BonferroniCIs[4,]),v=c(BonferroniCIs[1,]),lty=2)
abline(h=c(simultaneCIs[4,]),v=c(simultaneCIs[1,]),lty=3,col=4)

## part (e)

adjustalpha<-1-0.05/(2*7)
BonferroniCIs16<-cbind(Xbar-qt(adjustalpha,n-1)*sqrt(diag(Sn/n)),Xbar+qt(adjustalpha,n-1)*sqrt(diag(Sn/n)))
Var65<-t(c(-1,1))%*%Sn[c(5,6),c(5,6)]%*%c(-1,1)
BonferroniCIs7<-c(Xbar[6]-Xbar[5]-qt(adjustalpha,n-1)*sqrt(Var65/n),Xbar[6]-Xbar[5]+qt(adjustalpha,n-1)*sqrt(Var65/n))
BonferroniCIsAll<-rbind(BonferroniCIs16,BonferroniCIs7)

## Question 6.12

## part (b)

n1<-30
n2<-30
xbar1<-c(6.4,6.8,7.3,7.0)
xbar2<-c(4.3,4.9,5.3,5.1)
Spooled<-matrix(c(0.61,0.26,0.07,0.16,0.26,0.64,0.17,0.14,0.07,0.17,0.81,0.03,0.16,0.14,0.03,0.31),4,4)
Cmat<-matrix(c(1,-2,1,0,0,1,-2,1),2,4,byrow=TRUE)
Tn2<-t(xbar1+xbar2)%*%t(Cmat)%*%solve(Cmat%*%Spooled%*%t(Cmat)*(1/n1+1/n2))%*%Cmat%*%(xbar1+xbar2)

## Question 6.20

## part(a)

setwd("/Users/pszhong/Documents/MSU-Teaching/STT-843-Spring-2018")
birdsdatfemale<-read.table(file="birds-dat-table-512.txt")
birdsdatmale<-read.table(file="birds-dat-table-611.txt")
plot(birdsdatmale[,1],birdsdatmale[,2],xlab="Tail Length",ylab="Wing Length")

## part(b)

xbar1<-colMeans(birdsdatmale)
xbar2<-colMeans(birdsdatfemale)
n1<-nrow(birdsdatmale)
n2<-nrow(birdsdatfemale)
Spooled<-(cov(birdsdatfemale)+cov(birdsdatmale))*(44/88)
T2<-(n1*n2/(n1+n2))*t(xbar1-xbar2)%*%solve(Spooled)%*%(xbar1-xbar2)
p<-2
cutoff<-{(n1+n2-2)*p/(n1+n2-p-1)}*qf(0.95,p,n1+n2-p-1)

## removing data point #31 in male birds data

xbar1<-colMeans(birdsdatmale[-31,])
xbar2<-colMeans(birdsdatfemale)
n1<-nrow(birdsdatmale[-31,])
n2<-nrow(birdsdatfemale)
Spooled<-((n1-1)*cov(birdsdatmale[-31,])+(n2-1)*cov(birdsdatfemale))/(n1+n2-2)
T2outlierm<-(n1*n2/(n1+n2))*t(xbar1-xbar2)%*%solve(Spooled)%*%(xbar1-xbar2)
p<-2
cutoff<-{(n1+n2-2)*p/(n1+n2-p-1)}*qf(0.95,p,n1+n2-p-1)
ahat<-solve(Spooled)%*%(xbar1-xbar2)

## replace data point #31 in male birds data

birdsdatmale[31,1]<-184
xbar1<-colMeans(birdsdatmale)
xbar2<-colMeans(birdsdatfemale)
n1<-nrow(birdsdatmale)
n2<-nrow(birdsdatfemale)
Spooled<-(cov(birdsdatfemale)+cov(birdsdatmale))*(44/88)
T2<-(n1*n2/(n1+n2))*t(xbar1-xbar2)%*%solve(Spooled)%*%(xbar1-xbar2)
p<-2
cutoff<-{(n1+n2-2)*p/(n1+n2-p-1)}*qf(0.95,p,n1+n2-p-1)
ahat<-solve(Spooled)%*%(xbar1-xbar2)


## part (c)

xbar1<-colMeans(birdsdatmale[-31,])
xbar2<-colMeans(birdsdatfemale)
n1<-nrow(birdsdatmale[-31,])
n2<-nrow(birdsdatfemale)
Spooled<-((n1-1)*cov(birdsdatmale[-31,])+(n2-1)*cov(birdsdatfemale))/(n1+n2-2)
center<-xbar1-xbar2
x1<-seq(-12,0,length=100)
x2<-seq(-10,10,length=100)
plot(x1,x2,type="n",xlab="tail length difference",ylab="wing length difference")
points(center[1],center[2],col=2,pch=19)
npoints<-1000
r<-sqrt({(n1+n2-2)*p/(n1+n2-p-1)}*qf(0.95,p,n1+n2-p-1)*(n1+n2)/(n1*n2))
theta<-seq(0, 2*pi, length = npoints)
v<-rbind(r*cos(theta), r*sin(theta))
## transform for points on ellipse
z<-backsolve(chol(solve(Spooled)),v)+center
## plot points
lines(t(z))
adjustalpha<-1-0.05/4
BonferroniCIs<-cbind(center-qt(adjustalpha,n1+n2-2)*sqrt(diag((n1+n2)*Spooled/(n1*n2))),center+qt(adjustalpha,n-1)*sqrt(diag((n1+n2)*Spooled/(n1*n2))))
abline(h=c(BonferroniCIs[2,]),v=c(BonferroniCIs[1,]),lty=2)


## Question 6.24

setwd("/Users/pszhong/Documents/MSU-Teaching/STT-843-Spring-2018")
skulldata<-read.table(file="Egyptian-skulls.dat.txt")

## part1
## Obtain the response and design matrix
## Trtment 1: time1; Trtment 2: time2;
## Trtment 3: time3

Xn<-colMeans(skulldata[1:4])
trtment<-skulldata[,5]
responses<-as.matrix(skulldata[,1:4])
trtmentFac<-as.factor(trtment)

##
## Perform an MANOVA test using treatment as factor
##

manova1to4<-manova(responses~trtmentFac)
summary(manova1to4)

overallave<-colMeans(responses)
trtmeans<-matrix(0,3,4)
for (trt in 1:3)
{
  trtIND<-(skulldata[,5]==trt)
  trtmeans[trt,]<-colMeans(responses[trtIND,])
}
centrtmeans<-t(trtmeans)-overallave
B<-30*centrtmeans%*%t(centrtmeans)

residuals<-matrix(0,90,4)
for (i in 1:90)
{
  ithrow<-responses[i,]-trtmeans[skulldata[i,5],]
  residuals[i,]<-ithrow
}
W<-t(residuals)%*%residuals
Lambdastar<-det(W)/det(B+W)
TransLambda<-((90-4-2)/4)*(1-sqrt(Lambdastar))/sqrt(Lambdastar)

## Pairwise simultaneous confidence intervals

SimulCIintervals<-matrix(0,4,6)
alpha<-1-0.05/(4*6)
qtalpha<-qt(alpha,87)
for (k in 1:4)
{
  SimulCIintervals[k,1]<-trtmeans[1,k]-trtmeans[2,k]-qtalpha*sqrt((W[k,k]/87)*(1/15))
  SimulCIintervals[k,2]<-trtmeans[1,k]-trtmeans[2,k]+qtalpha*sqrt((W[k,k]/87)*(1/15))
  SimulCIintervals[k,3]<-trtmeans[1,k]-trtmeans[3,k]-qtalpha*sqrt((W[k,k]/87)*(1/15))
  SimulCIintervals[k,4]<-trtmeans[1,k]-trtmeans[3,k]+qtalpha*sqrt((W[k,k]/87)*(1/15))
  SimulCIintervals[k,5]<-trtmeans[2,k]-trtmeans[3,k]-qtalpha*sqrt((W[k,k]/87)*(1/15))
  SimulCIintervals[k,6]<-trtmeans[2,k]-trtmeans[3,k]+qtalpha*sqrt((W[k,k]/87)*(1/15))
}
SimulCIintervals

## Box’M test for homogeneity of covariances

library(biotools)
CovHomoTest<-boxM(responses,trtmentFac)

## Chi-square plot for normality 

Sn<-cov(residuals)
Sinv<-solve(Sn)
djs<-diag(residuals%*%Sinv%*%t(residuals))
qqplot(qchisq(ppoints(500), df = 4), djs,
 main = expression("Chi-square plot for" ~~ {chi^2}[4]))
qqline(djs,distribution=function(p) qchisq(p, df = 4))


## Question 6.37

setwd("/Users/pszhong/Documents/MSU-Teaching/STT-843-Spring-2018")
turtledata<-read.table(file="turtles-table-609.txt")
library(biotools)
CovHomoTest<-boxM(turtledata[,c(2:4)], turtle[,1])

## Question 6.39

## part (a)

setwd("//sttnas/FacultyReds/pszhong/My Documents/MSU-Teaching/STT-843-2018-Spring/Homework 2")
anacondadata<-read.table(file="anaconda-data-619.txt")
femalesample<-anacondadata[anacondadata[,3]=="F",]
malesample<-anacondadata[anacondadata[,3]=="M",]
xbar1<-colMeans(femalesample[,1:2])
xbar2<-colMeans(malesample[,1:2])
S1<-cov(femalesample[,1:2])
S2<-cov(malesample[,1:2])
Spooled<-(S1+S2)/2
n1<-28
n2<-28
T2<-(n1*n2)*t(xbar1-xbar2)%*%solve(Spooled)%*%(xbar1-xbar2)/(n1+n2)

## part (b)

library(biotools)
CovHomoTest<-boxM(anacondadata[,c(1:2)], as.factor(anacondadata[,3]))

## part (c)

adjustalpha<-1-0.05/(2*2)
Xbar<-xbar1-xbar2
BonferroniCIs<-cbind(Xbar-qt(adjustalpha,n1+n2-2)*sqrt(diag(S1/n1+S2/n2)),Xbar+qt(adjustalpha,n1+n2-2)*sqrt(diag(S1/n1+S2/n2)))