library(MASS)
caith
caith.m <- as.matrix(caith)
names(dimnames(caith)) <- c("Eyes","Hair")

mosaicplot(caith.m,color=T)
chisq.test(caith)

a <- dim(M)[1]
M <- as.matrix(caith)
M
a <- dim(M)[1]
b <- dim(M)[2]
n.. <- sum(M)
col.sums <- colSums(M)
col.proportions <- M%*%diag(1/col.sums)
row.sums <- rowSums(M)

D <- matrix(0,b,b)
for (i in 1:b){
	for (j in 1:b){
	D[i,j] <- n..*sum(((col.proportions[,i]-col.proportions[,j])^2)/row.sums)
	}
}

A <- -0.5*D
J <- matrix(1,b,b)
c.M <- diag(b) - J/b
B <- c.M%*%A%*%c.M
E <- eigen(B)

Z <- E$vectors[,1:3]%*%diag(sqrt(E$values[1:3]))

x <- c(Z[,1],Z[,1],Z[,2])
y <- c(Z[,2],Z[,3],Z[,3])
panl <- rep(c("1 vs 2","1 vs 3","2 vs 3"),each=b)
library(lattice)
xyplot(y~x|panl,col=1:b,pch=16)
legend(x="topright",legend=colnames(M),pch=16,col=1:b)

biplot(corresp(M,nf=2))


####################################################### Donner Party ###################################

Data <- read.table(file="C:/Documents and Settings/steeleb/My Documents/Math 543/DonnerParty.txt",header=T)
head(Data)
 

Age <- rep("Old",times=45)
Age[Data[,2]<31] <- "Young"
Age <- as.factor(Age)
Gender <- as.factor(Data[,3])
Survival <- as.factor(Data[,4])
DP <- data.frame(Age,Gender,Survival)

head(DP)
tabl <- table(DP)
chisq.test(tabl)

mosaicplot(tabl,color=T)

################## set up matrix for MCA ################
X <- cbind(as.integer(Age =="Old"),as.integer(Age =="Young"))
X <- cbind(X,as.integer(Gender == "MALE"),as.integer(Gender == "FEMALE"))
X <- cbind(X,as.integer(Survival=="DIED"),as.integer(Survival=="SURVIVED"))
head(X)
qr(X)$rank

a<-c(1,1,-1,-1,0,0)
X%*%a

b<-c(1,1,0,0,-1,-1)
X%*%b

d <- c(0,0,1,1,-1,-1)
X%*%d

cbind(a,b,d)
cbind(b-a,d)

colnames(X) <-c("Old","Young","MALE","FEMALE","DIED","SURVIVED")
cor(X)
p <- dim(X)[2]
S <- t(X)%*%X
S
p <- dim(S)[1]
e <- eigen(S)
m <- cbind(1:p,e$values,100*e$values/sum(e$values),100*cumsum(e$values)/sum(e$values))
colnames(m) <- c("e-pair","eigenvalue","%acc't by","% cumulative acc't by")
options(digits = 3)
m

vectrs <- e$vectors
rownames(vectrs) <-colnames(X)
colnames(vectrs) <- c("PC 1","PC 2","PC 3","PC 4","PC 5","PC 6")

Z <- S%*%vectrs[,1:3]%*%diag(1/sqrt(e$values[1:3]))
Z <- S%*%vectrs[,1:3]
par(mfrow=c(2,2))
mosaicplot(tabl,color=T)
plot(Z[,1:2],type="n")
text(Z[,1:2],labels=colnames(X))

plot(Z[,c(1,3)],type="n")
text(Z[,c(1,3)],labels=colnames(X))

plot(Z[,2:3],type="n")
text(Z[,2:3],labels=colnames(X))

################################### R function for MCA ##################
library(MASS)
par(mfrow=c(1,1))
plot(mca(DP))