## data entery in R
getwd()

#setwd("/Users/macbookpro/Bank-Mac/0-Desktop-2017")
#setwd()
d<-file.choose()
d
data1<-read.table(d,header=T)

# to read data frow txt format
data1<-read.table("data14.txt",header=TRUE)
data1

# to read excell files convert to csv and then:
 data1<-read.csv("data14.csv",header=TRUE)
data1
# mean

xbar<-apply(data1,2,mean)
xbar
# Covariance matrix
apply(data1,2,var)

var(data1)

S<-var(data1)
S

dim(S)

E<-eigen(S)
E$values

E$vectors

#S=sum(lambda*e1%*%t(e1))

L<-diag(E$values)
G<-E$vectors
G
G%*%L%*%t(G)

G%*%L%*%t(G)-S

round(G%*%L%*%t(G)-S,13)

# plot multivarate normal in p=2

# generating random sample from multivariate normal

# package MASS
library(MASS)
#?mvrnorm

Sigma<-matrix(c(1,.8,.8,1),2,2)
Sigma
mu<-c(0,0)
mu
samp1<-mvrnorm(n = 1000, mu, Sigma)
samp1

#plots

par(mfrow=c(1,2))
hist(samp1[,1],col="orange")
hist(samp1[,2],col="gray")


# install package mvtnorm
# install.packages("mvtnorm",dep=TRUE)
library(mvtnorm)
samp2<-rmvnorm(n=500, mu, Sigma)
samp2
par(mfrow=c(1,2))
hist(samp2[,1])
hist(samp2[,2])

x1<-samp2[,1]
x2<-samp2[,2]

#scatter plot
plot(x1,x2)

# Contour plots mvn

library(car)

dataEllipse(x1, x2, levels=c(.1,.5,0.95))
dataEllipse(x1, x2, levels=c(.1,.5,0.95),pch=16,cex=.5)
# Calculate kernel density estimate

bivar1 <- kde2d(x1, x2, n = 50) # from MASS in 50*50 grid
image(bivar1)
contour(bivar1,add=TRUE)

persp(bivar1,col = "lightblue",theta=90,phi=20)
#3d plot
rgl::persp3d(bivar1,col = "lightblue",theta=180)

# package fields

library(fields)

image.plot(bivar1)
image.plot(bivar1,col=rainbow(20))
image.plot(bivar1,col=topo.colors(30))
image.plot(bivar1,col=hcl.colors(30, "purples", rev = TRUE))
contour(bivar1,add=TRUE,col="gray60")

# Normality tests:
qqnorm(x1)
qqline(x1)

qqnorm(x2)
qqline(x2)

shapiro.test(x1)
shapiro.test(x2)

####### univariate test
#install.packages("nortest",dep=T)
library(nortest)
ad.test(x1) # Anderson-Darling test for normality 
lillie.test(x1)  #Lilliefors (Kolmogorov-Smirnov) test for normality 
pearson.test(x1) # Pearson chi-square test for normality 

# Multivarate normal test

#Shapiro-Wilk Multivariate Normality Test 
#install.packages("mvnormtest",dep=T)
library(mvnormtest)
mshapiro.test(t(samp2))

## data1
mshapiro.test(t(data1))


qqnorm(data1$x1)
qqline(data1$x1)

qqnorm(data1$x2)
qqline(data1$x2)

qqnorm(data1$x6)
qqline(data1$x6)

shapiro.test(data1$x6)

x2<-data1$x2
x6<-data1$x6

bivar1 <- kde2d(x2, x6, n = 50) 
image(bivar1)
contour(bivar1,add=TRUE)

persp(bivar1,col = "lightblue",theta=40,phi=20)

mshapiro.test(rbind(x2,x6))

#An R Package for Assessing Multivariate Normality 
#mult: mvnTest = c("mardia", "hz", "royston", "dh", "energy")
#uni: c("SW", "CVM", "Lillie", "SF", "AD") 
#install.packages("MVN",dep=T)
library(MVN)

#Mardia’s MVN test 
result <- mvn(samp2, mvnTest = "mardia")
result

#data x2 , x6
mvn(cbind(x2,x6), mvnTest = "mardia")
mvn(cbind(x2,x6), mvnTest = "hz")
mvn(cbind(x2,x6), mvnTest = "royston")
mvn(cbind(x2,x6), mvnTest = "energy")

mvn(data1)
mvn(data1, mvnTest = "energy")

mvn(data1[,c(2,4,6)])
mvn(data1[,c(2,4,6)], mvnTest = "energy")

##############################################################
### Example T2 Hotelling in lecture (jozve pivaste1-part2) ###
##############################################################

## Example page 24 (Jozve part2)

#creating data matrix p*n

r1<-c(6,10,8) #first row

r2<-c(9,6,3)  #second row
 
X<-rbind(r1,r2)  #combining rows
X

xbar<-apply(X,1,mean)
xbar


var(r1)
var(r2)
cov(r1,r2)

S<-var(t(X))
S


invS<-solve(S)
invS

# H0: mu=mu0=c(9,5)
# T2 Hotelling T2=n(xbar-mu0)'S^{-1}(xbar-mu0)
n=3 
p=2
mu0=c(9,5)
alpha=0.05
T2<-n*(xbar-mu0)%*%invS%*%as.matrix(xbar-mu0)

T2

c2<-((n-1)*p/(n-p))*qf(1-alpha,p,n-p)

T2>c2

### package ICSNP
#install.packages("ICSNP", "mvtnorm",dep=T)

library("ICSNP","mvtnorm")

?HotellingsT2

HotellingsT2(t(X),mu=c(9,5),test="f")

((n-1)*p/(n-p))^{-1}*T2 > qf(1-alpha,p,n-p)

# computing pvalues?
pval<-1-pf(((n-1)*p/(n-p))^{-1}*T2,p,n-p)
pval
# for large samples chisq approximate
HotellingsT2(t(X),mu=c(9,5),test="chi")


### without package
#?manova
y12 <- cbind(y1 = r1, y2 = r2)
y12
y<-sweep(y12, 2, mu0, "-")

result <- anova(lm(y ~ 1), test = "Hotelling-Lawley")
result

##########################
 ### example page 235 ###
##########################
# confidence region

r1<-c(6,10,8) #first row
r2<-c(9,6,3)  #second row
X<-rbind(r1,r2)  #combining rows
X

xbar<-apply(X,1,mean)
xbar
S<-var(t(X))
S

p=2
n=3

#install.packages("ellipse",dep=T)
library(ellipse)

?ellipse
# confidence 95%
#(x-Xbar)'S^(-1)(x-Xbar)=c^2/n
#c^2=qchisq(0.95, 2) for large data
alpha=0.05
c2<-(n-1)*p/(n-p)*qf(1-alpha,p,(n-p))

t0<-c2/n

plot(ellipse(S, centre = xbar, t=t0), type = "l")
title("Ellipse plot")
points(xbar[1],xbar[2],pch=18,col=2)
text(xbar[1]+20,xbar[2],expression(bar(x)),col=2,cex=1)

# chisq
plot(ellipse(S, centre = xbar, level = 0.95), type = "l")
title("Ellipse plot")
points(r1,r2,pch=16,col=4)
points(xbar[1],xbar[2],pch=18,col=2)
text(xbar[1]+.2,xbar[2],expression(bar(x)),col=2,cex=1)



#########################
## example 3-6 page 288
#########################

xbar1<-c(8.3,4.1)
xbar2<-c(10.2,3.9)
S1<-matrix(c(2,1,1,6),2)
S1
S2<-matrix(c(2,1,1,4),2)
S2
n1=50
n2=50
p=2

Sp<-((n1-1)*S1+(n2-1)*S2)/(n1+n2-2)
Sp

# H0: mu1=mu2 multivariate p=2
diff.x<-xbar1-xbar2
diff.x

Sinv<-solve(Sp)
Sinv

T2<-(n1*n2/(n1+n2))*diff.x%*%Sinv%*%as.matrix(diff.x)
T2
alpha=0.05
c2<-(n1+n2-2)*p/(n1+n2-p-1)*qf(1-alpha,p,(n1+n2-p-1))
c2

T2>c2

# RH0

# confidence region

E<-eigen(Sp)

lambda<-E$values


sqrt(lambda[1])*sqrt((1/n1+1/n2)*c2)

sqrt(lambda[2])*sqrt((1/n1+1/n2)*c2)

t0<-(1/n1+1/n2)*c2
t0
#

#ellipse plot
library(ellipse)

plot(ellipse(Sp, centre = diff.x, t=t0), ylim=c(-1,1),xlim=c(-3,1), type = "l")
title("Ellipse plot")
points(diff.x[1], diff.x[2],pch=18,col=2)
text(diff.x[1], diff.x[2]+.1,expression(bar(x1)-bar(x2)),col=2,cex=.6)
abline(h=0,v=0,lty=2,col=2)

#######################
## example 4-6 page 292
#######################
n1=45
n2=55
p=2

xbar1<-c(204.4,556.6)
xbar1
xbar2<-c(130,355)
xbar2

S1<-matrix(c(13825.3,23823.4,23823.4,73107.4),nrow=2)
S1
S2<-matrix(c(8632,19616.7,19616.7,55964.5),nrow=2)
S2

# equal variance assuming
Sp<-((n1-1)*S1+(n2-1)*S2)/(n1+n2-2)
Sp
alpha=0.05
c2<-(n1+n2-2)*p/(n1+n2-p-1)*qf(1-alpha,p,(n1+n2-p-1))
c2

# simultaneous confidence intervals
#mu1-mu2:?

dif.mu<-xbar1-xbar2


#mu11-mu21 for variable x1
SE.mu1<-sqrt(c2)*sqrt((1/n1+1/n2)*Sp[1,1])
L1<-dif.mu[1]-SE.mu1
L1
U1<-dif.mu[1]+SE.mu1
U1

c(L1,U1)
#mu12-mu22 for variable x2
SE.mu2<-sqrt(c2)*sqrt((1/n1+1/n2)*Sp[2,2])
L2<-dif.mu[2]-SE.mu2
L2
U2<-dif.mu[2]+SE.mu2
U2

c(L2,U2)

###
# confidence region

E<-eigen(Sp)

lambda<-E$values


sqrt(lambda[1])*sqrt((1/n1+1/n2)*c2)

sqrt(lambda[2])*sqrt((1/n1+1/n2)*c2)

t0<-(1/n1+1/n2)*c2
t0
#

#ellipse plot
library(ellipse)

plot(ellipse(Sp, centre = dif.mu, t=t0),ylim=c(0,350),xlim=c(0,200), type = "l")
title("Ellipse plot")
points(dif.mu[1], dif.mu[2],pch=18,col=2)
text(dif.mu[1], dif.mu[2]+10,expression(bar(x1)-bar(x2)),col=2,cex=.6)
abline(h=0,v=0,lty=2,col=2)


#######################
## example 5-6 page 295
#######################
#large samples

S<-S1/n1+S2/n2
S
ch<-qchisq(1-alpha,p)
#mu11-mu21 for variable x1
SE.mu1<-sqrt(ch)*sqrt(S[1,1])
L1<-dif.mu[1]-SE.mu1
L1
U1<-dif.mu[1]+SE.mu1
U1

c(L1,U1)
#mu12-mu22 for variable x2
SE.mu2<-sqrt(ch)*sqrt(S[2,2])
L2<-dif.mu[2]-SE.mu2
L2
U2<-dif.mu[2]+SE.mu2
U2

c(L2,U2)
#######################
## example 1-6 page 276
#######################
# two dependent samples
data16<-read.csv("data-1-6-p1.csv",header=T)

BOD1<-data16$BOD1
SS1<-data16$SS1
BOD2<-data16$BOD2
SS2<-data16$SS2
d1<-BOD1-BOD2
d2<-SS1-SS2
d<-rbind(d1,d2)

dbar<-apply(d,1,mean)

Sd<-var(t(d))
Sd
n=11
T2<-n*dbar%*%solve(Sd)%*%as.matrix(dbar)
T2
p=2
c2<-p*(n-1)/(n-p)*qf(1-alpha,p,n-p)

T2>c2

#######################
## example 9-6 page 308
#######################
# MANOVA
y1<-c(9,6,9,0,2,3,1,2)
y2<-c(3,2,7,4,0,8,9,7)
Factor<-rep(c(1,2,3),c(3,2,3))
Factor
response<-cbind(y1,y2)

Factor<-as.factor(Factor)

fit<-manova(response~Factor)
fit

?summary.manova

summary(fit,test="Wilks")

M<-summary(fit,test="Wilks")

names(M)

M$SS

B<-M$SS$Factor
B

W<-M$SS$Residuals
W

Total=B+W

Total