#CANONICAL CORRELATION ANALYSIS
#data for prototype from:
#http://www.ats.ucla.edu/stat/r/dae/canonical.htm

getwd()
setwd("c:/DATA/Multivariate")
getwd()

library(CCA)
library(yacca)

M=read.table("UCLAidreCCAEx.txt")
M
summary(M)

psyc=M[,1:3]
acad=M[,4:8]

#using {yacca}
CCA=cca(psyc,acad)
CCA

#standardized dataset:
Ms=as.matrix(scale(M))
summary(Ms)
write.table(Ms,"Mstand.txt")

#CCA standardized analysis:
CCAs=cca(psyc,acad,xscale=TRUE,yscale=TRUE) 
CCAs

#automated plots of cca() objects:
#press escape in R to stop at any point
plot(CCAs)

#helio plots with more control:
helio.plot(CCAs,cv=1,type="correlation")
helio.plot(CCAs,cv=2,type="variance",x.name="psyc",y.name="acad",main="HP Variance")
?helio.plot

#plot() will give also give output...
#plotting CV scores directly:
#looking at standardized scaled coefficients:
CCAs$xcoef
CCAs$ycoef
#looking at the scores:
CCAs$canvarx
CCAs$canvary

#plot for CV 1:
plot(CCAs$canvarx[,1],CCAs$canvary[,1])
#plot for CV 2:
plot(CCAs$canvarx[,2],CCAs$canvary[,2])


F.test.cca(CCA)


+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#clarification of CCA coefficient scaling:
#http://www.statpower.net/Content/312/Lecture%20Slides/CanonicalCorrelation.pdf

#X <- matrix(c( 1,1,3,2,3,2,1,1,1,1,1,2,2,2,3,3,3,2,1,3,2,4,3,5,5,5,5),9,3,byrow=T)

#Y <- matrix(c( 4,4,-1.07846,3,3,1.214359,
#+ 2,2,0.307180,2,3,-0.385641,
#+ 2,1,-0.078461,1,1,1.61436,
#+ 1,2,0.814359,2,1,-0.0641016,
#+ 1,2,1.535900 ),9,3,byrow=T)

X=as.matrix(psyc)
Y=as.matrix(acad)

S.xy <- cov(X,Y)
S.xx <- var(X)
S.yx <- cov(Y,X)
S.yy <- var(Y)

#"raw" coefficients A & B
A <- eigen(solve(S.xx) %*% S.xy %*% solve(S.yy) %*% S.yx)$vectors
A
B <- eigen(solve(S.yy) %*% S.yx %*% solve(S.xx) %*% S.xy)$vectors
B
R <- sqrt(eigen(solve(S.yy) %*% S.yx %*% solve(S.xx) %*% S.xy)$values)
R

# Singly standardized weights (SAS "raw")
A.single <- A %*% solve(sqrt(diag(diag(var(X %*% A)))))
A.single
B.single <- B %*% solve(sqrt(diag(diag(var(Y %*% B)))))
B.single


#Deviation score X,Y
UnitVector<-function(n){matrix(rep(1,n),n,1)}
sympower <- function(x,pow) {
edecomp <- eigen(x)
roots <- edecomp$val
v <- edecomp$vec
d <- roots^pow
if(length(roots)==1) d <- matrix(d,1,1) else d <- diag(d)
sympow <- v %*% d %*% t(v)
sympow
}

P <-function(x) {
y <- x %*% solve(t(x) %*% x) %*% t(x)
y
}

Q <- function(x) {
q <- diag(dim(x)[1]) - P(x)
q
}
X.dev <- Q(UnitVector(length(X[,1]))) %*% X
Y.dev <- Q(UnitVector(length(X[,1]))) %*% Y

#Z-score X,Y
# Create diagonal matrices with standard deviations
# Then invert using solve
D.x <- solve(sqrt(diag(diag(var(X)))))
D.y <- solve(sqrt(diag(diag(var(Y)))))
# Postmultiply the deviation score matrix
# to create Z-scores
Z.x <- X.dev %*% D.x
Z.y <- Y.dev %*% D.y

R.xy <- cor(X,Y)
R.xx <- cor(X)
R.yx <- cor(Y,X)
R.yy <- cor(Y)
A.s <- eigen(solve(R.xx) %*% R.xy %*% solve(R.yy) %*% R.yx)$vectors
B.s <- eigen(solve(R.yy) %*% R.yx %*% solve(R.xx) %*% R.xy)$vectors
A.fully <- A.s %*% solve(sqrt(diag(diag(var(Z.x %*% A.s)))))
A.fully
B.fully <- B.s %*% solve(sqrt(diag(diag(var(Z.y %*% B.s)))))
B.fully