lazer <- read.table(
url("http://stat.unm.edu/~fletcher/LLM/DATA/TAB15-1.DAT"),
# "C:\\E-drive\\Books\\ANREG2\\newdata\\tab5-6.dat",
#"C:\\E-drive\\Books\\LOGLIN3\\DATA\\TAB15-1.dat",
		  sep="",col.names=c("O","Rel","Occ"))
attach(lazer)
lazer
laz <- xtabs(O~Rel+Occ)
laz

# fit indep/homogen model to table
fit <- chisq.test(laz,correct=FALSE)
fit
fit$expected
fit$residual

# Singular Value Decomp of sqrt(n) x residual matrix
caa=svd(fit$residual)
L=diag(caa$d)
L
U=caa$u
U
V=caa$v
V
#This should reproduce fit$residual
PR=U%*%L%*%t(V)
PR

# Pearson Chi-square
X2=PR%*%t(PR)
sum(diag(X2))

#Look at how closely 2-dim SVD reproduces residual matrix
PRca=U[,1:2]%*%L[1:2,1:2]%*%t(V[,1:2])
PRca

#Figure 1
par(mfrow=c(2,1))
#SVD Plot
Ut=U%*%sqrt(L)
Vt=V%*%sqrt(L)
# If done correctly the next line should be
# the residual matrix.
Ut%*%t(Vt)
plot(Ut[,1],Ut[,2],
    xlim=c(-2.5,2.5),ylim=c(-2.5,2.5),
     xlab="Factor 1",ylab="Factor 2",main="SVD Plot")
text(Ut[,1],Ut[,2]-.4,labels=c("Protest","RomCath","Jewish"))
lines(Vt[,1],Vt[,2],pch=15,type="p")
text(Vt[,1]-.125,Vt[,2],labels=c("A","B","C","D"))
text(0,0,labels=c("+"))

#Correspondence Analysis Plot
Pr=rowSums(laz)/sum(laz)
Pc=colSums(laz)/sum(laz)
Dr=diag(1/sqrt(Pr))
Dc=diag(1/sqrt(Pc))
Fr=Dr%*%U%*%L/sqrt(sum(laz))
Fc=Dc%*%V%*%L/sqrt(sum(laz))
plot(Fr[,1],Fr[,2],
    xlim=c(-.7,.3),ylim=c(-.3,.3),
     xlab="Factor 1",ylab="Factor 2",main="CA Plot")
text(Fr[,1],Fr[,2]-.06,labels=c("Protest","RomCath","Jewish"))
lines(Fc[,1],Fc[,2],pch=15,type="p")
text(Fc[,1]-.03,Fc[,2],labels=c("A","B","C","D"))
text(0,0,labels=c("+"))
par(mfrow=c(1,1))


# Figure 2
#SVD Plot
Ut=U%*%sqrt(L)
Vt=V%*%sqrt(L)
# If done correctly the next line should be
# the residual matrix.
Ut%*%t(Vt)
plot(Ut[,1],Ut[,2],
    xlim=c(-2.5,2.5),ylim=c(-2.5,2.5),
     xlab="Factor 1",ylab="Factor 2",main="SVD Plot")
text(Ut[,1],Ut[,2]-.15,labels=c("Protest","RomCath","Jewish"))
lines(Vt[,1],Vt[,2],pch=15,type="p")
text(Vt[,1]-.15,Vt[,2],labels=c("A","B","C","D"))
text(0,0,labels=c("+"))

# Figure 3 is from ANREG2.

oc=c("A","B","C","D")
i=seq(1,4)
P=c(.185,.244,.224,.347)
RC=c(.157,.216,.196,.431)
J=c(.252,.420,.210,.119)
plot(i,P,type="b",xlab="Occupation",ylab="Proportion",at=F,ylim=c(0.10,0.47))
lines(i,RC,type="b",lty=2)
lines(i,J,type="b",lty=3)
axis(side=1,at=i,labels=oc)
axis(side=2)
legend("topleft",c("Protest","RomCath","Jewish"),lty=c(1,2,3))


# Figure 4
plot(Fr[,1],Fr[,2],
    xlim=c(-.575,.23),ylim=c(-.575,.23),
     xlab="Factor 1",ylab="Factor 2",main="CA Rows Plot")
text(Fr[,1],Fr[,2]-.03,labels=c("Protest","RomCath","Jewish"))
text(0,0,labels=c("+"))

#Figure 5
plot(Fc[,1],Fc[,2],
    xlim=c(-.25,.25),ylim=c(-.25,.25),
     xlab="Factor 1",ylab="Factor 2",main="CA Columns Plot")
text(Fc[,1]-.015,Fc[,2],labels=c("A","B","C","D"))
text(0,0,labels=c("+"))

#Figure 6
plot(Fr[,1],Fr[,2],
      xlim=c(-.575,.23),ylim=c(-.575,.23),
     xlab="Factor 1",ylab="Factor 2",main="CA Plot")
text(Fr[,1],Fr[,2]-.03,labels=c("Protest","RomCath","Jewish"))
lines(Fc[,1],Fc[,2],pch=15,type="p")
text(Fc[,1]-.03,Fc[,2],labels=c("A","B","C","D"))

# Or you could just use the ca package
#install("ca")
library(ca)
fitt=ca(laz)
plot(fitt)
par(mfrow=c(1,1))