# read in women data from 1994 (concatenated matrix)

women94 <- read.table("clipboard", header=T, row.names=1, check.names=F)
library(ca)
plot(ca(women94))
plot(ca(women94), arrows=c(F,T), map="rowgreen")

# to reverse an axis (e.g., 2nd)
ca.women94 <- ca(women94)
ca.women94$rowcoord[,2] <- -ca.women94$rowcoord[,2]
ca.women94$colcoord[,2] <- -ca.women94$colcoord[,2]
plot(ca.women94, arrows=c(F,T), map="rowgreen")

# read original data for Spain, 2002
women02 <- read.table("clipboard", header=T)

# getting the stacked matrix from the Burt matrix
women02.B <- mjca(women02)$Burt
rownames(women02.B)
colnames(women02.B)

women02.stack <- women02.B[49:69,1:48]
rownames(women02.stack)
colnames(women02.stack)

# checking frequencies of each substantive category
apply(women02.stack,2,sum)/4

# combine H4 and H5 (H5 frequency very low)
women02.stack[,46] <- women02.stack[,46]+women02.stack[,47]
women02.stack <- women02.stack[,-47]

# CA of whole stacked matrix
plot(ca(women02.stack))

# checking frequencies of each demographic category
apply(women02.stack,1,sum)/8

# subset analysis (subsetting out the missing categories)
plot(ca(women02.stack, subsetrow=seq(1,21)[-c(8,15)], subsetcol=seq(1,47)[-c(6,12,18,24,30,36,42,47)]), map="rowprincipal")

# subset analysis, contribution biplot
plot(ca(women02.stack, subsetrow=seq(1,21)[-c(8,15)], subsetcol=seq(1,47)[-c(6,12,18,24,30,36,42,47)]), map="rowgreen", arrows=c(F,T))

summary(ca(women02.stack, subsetrow=seq(1,21)[-c(8,15)], subsetcol=seq(1,47)[-c(6,12,18,24,30,36,42,47)]))


# use wg93 data set on science & environment
data(wg93)

# CA of indicator matrix using mjca function
ca.wg93.Z <- mjca(wg93[,1:4], lambda="indicator")
plot(ca.wg93.Z, labels=c(0,2), map="rowprincipal")  

# inertias and percentages of inertia
ca.wg93.Z$sv^2
100*ca.wg93.Z$sv^2/sum(ca.wg93.Z$sv^2)

# Burt matrix
wg93.B <- mjca(wg93)$Burt

# CA of Burt matrix using ca function
ca.wg93.B <- ca(wg93.B[1:20,1:20])
ca.wg93.B$rowcoord <- - ca.wg93.B$rowcoord
ca.wg93.B$colcoord <- - ca.wg93.B$colcoord
plot(ca.wg93.B, map="rowprincipal")  

# inertias and percentages of inertia
ca.wg93.B$sv^2
100*ca.wg93.B$sv^2/sum(ca.wg93.B$sv^2)

# inertias in the subtables
sum(mjca(wg93[,1:4], lambda="Burt")$subinertia*16)
# total inertia of Burt matrix...
sum(ca.wg93.B$sv^2)
# ...equals average of subinertia tables
sum(mjca(wg93[,1:4], lambda="Burt")$subinertia*16)/16
# average inertia of off-diagonal tables
(sum(mjca(wg93[,1:4], lambda="Burt")$subinertia*16)-16)/12

# adjusted inertias
(16/9)*(0.4574-1/4)^2
(16/9)*(0.4310-1/4)^2
# adjusted parts of inertia
(16/9)*(0.4574-1/4)^2/0.1703
(16/9)*(0.4310-1/4)^2/0.1703

ca.wg93.adj <- ca.wg93.B
ca.wg93.adj$sv[1] <- (4/3)*(0.4574-1/4)
ca.wg93.adj$sv[2] <- (4/3)*(0.4310-1/4)
plot(ca.wg93.adj, map="rowprincipal")

jca.wg93 <- mjca(wg93[,1:4], lambda="JCA")
summary(jca.wg93)
plot(jca.wg9 3, what=c("none","all"))
jca.wg93$rowcoord <- - jca.wg93$rowcoord
jca.wg93$colcoord <- - jca.wg93$colcoord
plot(jca.wg93, what=c("none","all"), map="rowprincipal")

# Cronbach alpha
(4/3)*(1-1/(4*0.4574))

# inertias without variable D
# Cronbach alpha again
mjca(wg93[,1:3], lambda="indicator")$sv^2
(4/3)*(1-1/(4*0.6018))