General information

Below you can find all the R commands given in Appendix B of Correspondence Analysis in Practice, Second Edition, starting from page 215.

You can find the data sets as Excel files in the section 'data sets'

R code

# Chapter 2: Profiles and the Profile Space

# The data file 'travels.xls' contains the original count data
# so first calculate the row profiles in Excel by dividing each
# row by its row total.  Then copy the profiles with their row
# and column names, remembering to leave the top left hand cell
# empty.
# Alternatively, copy the five lines below 

               Holidays HalfDays FullDays
Norway            0.333    0.056    0.611
Canada            0.067    0.200    0.733
Greece            0.138    0.862    0.000
France/Germany    0.083    0.083    0.833

# To read the data into R the following statement must be typed 
# in the R command window and not copied and pasted (since the 
# data is in the "clipboard")

profiles <- read.table("clipboard")
profiles

# You need to have installed and loaded the rgl package now

rgl.lines(c(0,1.2), c(0,0),   c(0,0))
rgl.lines(c(0,0)  , c(0,1.2), c(0,0))
rgl.lines(c(0,0)  , c(0,0),   c(0,1.2))

rgl.lines(c(0,0), c(0,1), c(1,0), size = 2)
rgl.lines(c(0,1), c(1,0), c(0,0), size = 2)
rgl.lines(c(0,1), c(0,0), c(1,0), size = 2)

rgl.points(profiles[,3], profiles[,1], profiles[,2], size = 4)
rgl.texts(profiles[,3], profiles[,1], profiles[,2], text = row.names(profiles))

# Now expand the 3-d graphics window and spin the points around using the
# left mouse button; the right button or mouse wheel allows zooming

# Chapter 3: Masses and Centroids

# Copy data in file 'readers.xls' or copy following lines

   C1 C2 C3
E1  5  7  2
E2 18 46 20
E3 19 29 39
E4 12 40 49
E5  3  7 16

# Again, TYPE the following statement

table <- read.table("clipboard")

table

table.pro <- table/apply(table,1,sum)
table.x   <- 1-table.pro[,1]-table.pro[,3]/2
table.y   <- table.pro[,3]*sqrt(3)/2
plot.new()
lines(c(0,1,0.5,0), c(0,0,sqrt(3)/2,0), col = "gray")
text(c(0,1,0.5), c(0,0,sqrt(3)/2), labels = colnames(table))
text(table.x, table.y, labels = rownames(table))

# Chapter 4: Chi-square Distance and Inertia

# calculate chi-square statistic

table.rowsum <- apply(table, 1, sum)
table.colsum <- apply(table, 2, sum)
table.sum    <- sum(table)
table.exp    <- table.rowsum%o%table.colsum/table.sum
chi2         <- sum((table-table.exp)^2/table.exp)
chi2

# chi-square distances to centroid

table.colmass <- table.colsum/table.sum
chidist       <- 0
for(j in 1:dim(table)[2]){
  chidist <- chidist+(table.pro[5,j]-table.colmass[j])^2/table.colmass[j]
  }
chidist

apply((t(table.pro)-table.colmass)^2/table.colmass, 2, sum)

tablec.pro <- rbind(table.pro,table.colmass)
rownames(tablec.pro)[6] <- "ave"
dist(sweep(tablec.pro, 2, sqrt(table.colmass), FUN = "/"))

# Chapter 5: Plotting Chi-square Distance

table.wt <- sqrt(table.colmass)

rgl.lines(c(0,1.2/table.wt[2]), c(0,0), c(0,0))
rgl.lines(c(0,0), c(0,1.2/table.wt[3]), c(0,0))
rgl.lines(c(0,0), c(0,0), c(0,1.2/table.wt[1]))

rgl.lines(c(0,0), c(0,1/table.wt[3]), c(1/table.wt[1],0), size = 2)
rgl.lines(c(0,1/table.wt[2]), c(1/table.wt[3],0), c(0,0), size = 2)
rgl.lines(c(0,1/table.wt[2]), c(0,0), c(1/table.wt[1],0), size = 2)

table.chi <- t(t(table.pro)/table.wt)
rgl.points(table.chi[,2], table.chi[,3], table.chi[,1], size = 4)
rgl.texts(table.chi[,2], table.chi[,3], table.chi[,1], text = row.names(table))

# Chapter 6: Reduction of Dimensionality

# First read data from file 'health.xls' as described before into
# data frame health

health.P    <- health/sum(health)
health.r    <- apply(health.P, 1, sum)
health.c    <- apply(health.P, 2, sum)
health.Dr   <- diag(health.r)
health.Dc   <- diag(health.c)
health.Drmh <- diag(1/sqrt(health.r))
health.Dcmh <- diag(1/sqrt(health.c))

health.P   <- as.matrix(health.P)
health.S   <- health.Drmh%*%(health.P-health.r%o%health.c)%*%health.Dcmh
health.svd <- svd(health.S)

health.rsc <- health.Drmh%*%health.svd$u
health.csc <- health.Dcmh%*%health.svd$v
health.rpc <- health.rsc%*%diag(health.svd$d)
health.cpc <- health.csc%*%diag(health.svd$d)

plot(health.rpc[,1], health.rpc[,2], type = "n", pty = "s")
text(health.rpc[,1], health.rpc[,2], label = rownames(health))

# Notice that the R graphics window is square and does not give the 
# correct aspect ratio, so pull out the window so that a vertical 
# unit is approximately equal to a horizontal unit: in this example
# the window would be as wide as the computer screen but take up less
# than half the height of the screen. 

# Chapter 7: Optimal Scaling

health.range <- max(health.csc[,1])-min(health.csc[,1])
health.scale <- (health.csc[,1]-min(health.csc[,1]))*100/health.range
health.scale

# Chapter 8: Symmetry of Row and Column Analyses

health.cor <- t(health.rsc[,1])%*%health.P%*%health.csc[,1]

health.cor^2

t(health.rsc[,1])%*%health.Dr%*%health.rsc[,1]

# Chapter 9: Two-dimensional Maps

# For this example you can read in the data from 'smoke.xls' as
# before, or after installing and loading the ca package, the 
# data are available by simply entering:

data(smoke)
smoke
ca(smoke)
plot(ca(smoke))

# (After the previous statement you might need to adjust the 
#  aspect ratio of the graphics window)

plot(ca(smoke), map = "rowprincipal")

# Chapter 10: Three More Examples

data(author)
plot3d.ca(ca(author), labels = c(2,1), sf = 0.000001)

# Chapter 11: Contributions to Inertia

# Read in data from 'funding.xls' into data frame fund, then:

fund.P    <- as.matrix(fund/sum(fund))
fund.r    <- apply(fund.P,1,sum)
fund.c    <- apply(fund.P,2,sum)
fund.Drmh <- diag(1/sqrt(fund.r))
fund.Dcmh <- diag(1/sqrt(fund.c))
fund.res  <- fund.Drmh%*%(fund.P-fund.r%o%fund.c)%*%fund.Dcmh
round(apply(fund.res^2, 1, sum), 5)

round(apply(fund.res^2, 2, sum), 5)

round(1000*fund.res^2/sum(fund.res^2), 0)

fund.svd <- svd(fund.res)
fund.svd$d^2

fund.F    <- fund.Drmh%*%fund.svd$u%*%diag(fund.svd$d)
fund.rowi <- diag(fund.r)%*%fund.F^2
fund.rowi[,1:4]

round(1000*(fund.rowi/apply(fund.rowi,1,sum))[,1:4], 0)

round(1000*t(t(fund.rowi)/fund.svd$d^2)[,1:4], 0)

summary(ca(fund))

# Chapter 12: Supplementary Points

fund.m     <- c(4,12,11,19,7)/53
fund.Gamma <- fund.Dcmh%*%fund.svd$v
t(fund.m)%*%fund.Gamma[,1:2]

# Chapter 13: Correspondence Analysis Biplots

diag(sqrt(fund.c))%*%fund.Gamma[,1:2]

fund.est <- fund.F[,1:2]%*%t(diag(sqrt(fund.c))%*%fund.Gamma[,1:2])
oner     <- rep(1, dim(fund)[1])
round(fund.est%*%diag(sqrt(fund.c))+oner%o%fund.c,3)

round(fund.P/fund.r, 3)

sum(diag(fund.r)%*%(fund.est%*%diag(sqrt(fund.c))+oner%o%fund.c - 
     fund.P/fund.r)^2%*%diag(1/fund.c))

sum(fund.svd$d[3:4]^2)

# Chapter 14: Transition and Regression Relationships

fund.vec  <- as.vector(as.matrix(fund))
fund.Phi  <- fund.Drmh%*%fund.svd$u
fund.vecr <- rep(fund.Phi[,1], 5)
fund.vecc <- as.vector(oner%*%t(fund.Gamma[,1]))
lm(fund.vecr~fund.vecc, weight = fund.vec)

fund.y <- (fund.P[1,]/fund.r[1])/fund.c
summary(lm(fund.y~fund.Gamma[,1]+fund.Gamma[,2], weight = fund.c))

cov.wt(cbind(fund.y,fund.Gamma[,1:2]), wt = fund.c, cor = TRUE)$cor

# Chapter 15: Clustering Rows and Columns

# The food stores data is in file 'stores.xls', to be read 
# into R data frame food; then:

food.exp   <- apply(food,1,sum)%o%apply(food,2,sum)/sum(food)
food.stres <- (food-food.exp)/sqrt(food.exp)
sum(food.stres^2)

food.rpro <- food/apply(food, 1, sum)
food.r    <- apply(food, 1, sum)/sum(food)

# using Murtagh's hierclust() function (see web resource p.261)
# to row profiles in food.rpro, using weights in food.r

food.rclust <- hierclust(food.rpro, food.r)
plot(as.dendrogram(food.rclust))

# Chapter 16: Multiway Tables

# The raw data file of 6371 cases for the Spanish health survey 
# has not been given but the derived files are: 
# 'health.xls': cross-tabulation of age by perceived health
# 'health2.xls': cross-tabulation of age-sex by perceived health
#                (Exhibit 16.2)

# The working women raw data are in 'womenraw.xls'
# Since the top lefthand cell is not empty (there are no row names)
# the following statement needs to be typed to read the copied data

women <- read.table("clipboard", header = TRUE)

# Check size of this data frame as follows:

dim(women)

colnames(women)
attach(women)
table(C, Q3)

# The above table is also given in file 'women.xls'

CG <- 2*(C-1)+G
table(CG, Q3)

# these commands declare gender 9 to be NA

women[,6][G==9] <- NA
attach(women)
CG <- 2*(C-1)+G
table(CG, Q3)

# The above table is also given in file 'women2.xls'

CGA <- 6*(CG-1)+A

# Chapter 17: Stacked Tables

women.stack <- table(C,Q3)
for(j in 6:9){
  women.stack<-rbind(women.stack,table(women[,j],Q3))
  }
women.stack <- women.stack[-c(38,39,47,48),]

# The above table is also given in file 'women3.xls' (Exhibit 17.1)

chisq.test(women.stack[27:32,])$statistic/sum(women.stack[27:32,])

sum(ca(women.stack[27:32,])$sv^2)

# The file which concatenates all cross-tabulations (Exhibit 17.3)
# is given in 'women4.xls' 

# Chapter 18: Multiple Correspondence Analysis

# Assuming raw women data is in data frame women (see above) and that
# attach(women) has been executed...

germany <- C==2|C==3
womenG  <- women[germany,]
dim(womenG)

# Elimination of missing data

attach(womenG)
missing <- G==9|M==9|E==98|E==99
womenG  <- (womenG[!missing,])
dim(womenG)
womenG  <- as.matrix(womenG) 

mjca(womenG[,1:4], lambda = "indicator")

mjca(womenG[,1:4], lambda = "Burt")

sum(mjca(womenG[,1:4], lambda = "indicator")$sv^2)

sum(mjca(womenG[,1:4], lambda = "Burt")$sv^2)

sum(mjca(womenG[,1:4], lambda = "Burt")$subinertia)

16*mjca(womenG[,1:4], lambda = "Burt")$subinertia

summary(mjca(womenG, lambda = "Burt", supcol = 5:9))

# The following computes the Burt matrix, then applies ca() 
# to give the same results as the option lambda="Burt" above

womenG.B <- mjca(womenG[,1:9], lambda = "Burt")$Burt
summary(ca(womenG.B[,1:16], suprow = 17:38))

# Chapter 19: Joint Correspondence Analysis

summary(mjca(womenG[,1:4], lambda = "JCA"))

summary(mjca(womenG[,1:4]))

# The following statements perform the weighted least-squares
# regression of the off-diagonal elements of the Burt matrix
# (in the form of contingency ratios) on the predictors 
# constructed from the standard coordinates. 

# calculate contingency ratios

womenG.mca  <- mjca(womenG[,1:4], lambda = "Burt")
womenG.p    <- as.matrix(womenG.mca$Burt)/sum(womenG.mca$Burt)
womenG.mass <- womenG.mca$colmass
womenG.cr   <- diag(1/womenG.mass)%*%womenG.p%*%diag(1/womenG.mass)-1

# calculate "predictors"

womenG.x1 <- womenG.mca$colcoord[,1]%o%womenG.mca$colcoord[,1]
womenG.x2 <- womenG.mca$colcoord[,2]%o%womenG.mca$colcoord[,2]
womenG.wt <- womenG.mass%o%womenG.mass
for(q in 1:4){
  womenG.cr[(4*(q-1)+1):(4*q),(4*(q-1)+1):(4*q)] <- NA
  womenG.x1[(4*(q-1)+1):(4*q),(4*(q-1)+1):(4*q)] <- NA
  womenG.x2[(4*(q-1)+1):(4*q),(4*(q-1)+1):(4*q)] <- NA
  womenG.wt[(4*(q-1)+1):(4*q),(4*(q-1)+1):(4*q)] <- NA
  }
womenG.cr <- as.vector(womenG.cr[!is.na(womenG.cr)])
womenG.x1 <- as.vector(womenG.x1[!is.na(womenG.x1)])
womenG.x2 <- as.vector(womenG.x2[!is.na(womenG.x2)])
womenG.wt <- as.vector(womenG.wt[!is.na(womenG.wt)])

summary(lm(womenG.cr~womenG.x1+womenG.x2, weight = womenG.wt))

# In the above we see that the regression explains 89.95% of the 
# inertia. As expected this is slightly better than the simplified
# "adjusted" solution and slightly worse than the optimal "JCA" one

# Chapter 20: Scaling Properties of MCA

# Here we use the data set wg93 provided in the ca package

data(wg93)

wg93.mca <- mjca(wg93[,1:4], lambda = "indicator")
plot(wg93.mca, what = c("none","all"))

wg93.F    <- wg93.mca$colcoord%*%diag(sqrt(wg93.mca$sv))
wg93.coli <- diag(wg93.mca$colmass)%*%wg93.F^2
matrix(round(1000*wg93.coli[,1]/wg93.mca$sv[1]^2, 0), nrow = 5)

Ascal <- wg93.mca$colcoord[1:5,1]
Bscal <- wg93.mca$colcoord[6:10,1]
Cscal <- wg93.mca$colcoord[11:15,1]
Dscal <- wg93.mca$colcoord[16:20,1]
As    <- Ascal[wg93[,1]]
Bs    <- Bscal[wg93[,2]]
Cs    <- Cscal[wg93[,3]]
Ds    <- Dscal[wg93[,4]]
AVEs  <- (As+Bs+Cs+Ds)/4

cor(cbind(As,Bs,Cs,Ds,AVEs))^2

sum(cor(cbind(As,Bs,Cs,Ds,AVEs))[1:4,5]^2)/4

wg93.mca$sv[1]^2

cov(cbind(As,Bs,Cs,Ds))*870/871

mean(cov(cbind(As,Bs,Cs,Ds))*870/871)

sum(diag(cov(cbind(As,Bs,Cs,Ds))*870/871))

VARs <- ((As-AVEs)^2+(Bs-AVEs)^2+(Cs-AVEs)^2+(Ds-AVEs)^2)/4
mean(VARs)

plot(wg93.mca, map="rowprincipal", labels = c(0,2))

# Note: the above plot may come out inverted compared to Exhibit 20.3

# Chapter 21: Subset Correspondence Analysis

data(author)
vowels     <- c(1,5,9,15,21)
consonants <- c(1:26)[-vowels]
summary(ca(author, subsetcol = consonants))

summary(ca(author, subsetcol = vowels))

plot(ca(author, subsetcol = consonants), map = "rowgreen")

plot(ca(author,subsetcol = vowels), map = "rowgreen")

womenG.B <- mjca(womenG[,1:4])$Burt
subset   <- c(1:16)[-c(4,8,12,16)]
summary(ca(womenG.B[1:16,1:16], subsetrow = subset, subsetcol = subset))

# calculate contingency ratios

womenG.sca  <- ca(womenG.B[1:16,1:16], subsetrow = subset, subsetcol = subset)
womenG.p    <- as.matrix(womenG.mca$Burt)/sum(womenG.mca$Burt)
womenG.mass <- womenG.mca$colmass
womenG.cr   <- diag(1/womenG.mass)%*%womenG.p%*%diag(1/womenG.mass)-1

# calculate "predictors"

womenG.x1 <- womenG.sca$colcoord[,1]%o%womenG.sca$colcoord[,1]
womenG.x2 <- womenG.sca$colcoord[,2]%o%womenG.sca$colcoord[,2]
womenG.wt <- womenG.mass%o%womenG.mass
for(q in 1:4){
  womenG.cr[(4*(q-1)+1):(4*q-1),(4*(q-1)+1):(4*q-1)] <- NA
  womenG.cr[4*q,] <- NA; womenG.cr[,4*q] <- NA 
  womenG.x1[(3*(q-1)+1):(3*q),(3*(q-1)+1):(3*q)] <- NA
  womenG.x2[(3*(q-1)+1):(3*q),(3*(q-1)+1):(3*q)] <- NA
  womenG.wt[(4*(q-1)+1):(4*q-1),(4*(q-1)+1):(4*q-1)] <- NA
  womenG.wt[4*q,] <- NA; womenG.wt[,4*q] <- NA 
  }
womenG.cr <- as.vector(womenG.cr[!is.na(womenG.cr)])
womenG.x1 <- as.vector(womenG.x1[!is.na(womenG.x1)])
womenG.x2 <- as.vector(womenG.x2[!is.na(womenG.x2)])
womenG.wt <- as.vector(womenG.wt[!is.na(womenG.wt)])

summary(lm(womenG.cr~womenG.x1+womenG.x2, weight = womenG.wt))

# The above regression shows that the solution explains 86.7% 
# of the off-diagonal inertia in the subset, compared to the 
# percentage of 62.4% reported in the subset analysis of the 
# Burt matrix

# Chapter 22: Analysis of Square Tables

# Read in data from file 'mobility.xls' into data frame mob

mob  <- as.matrix(mob)
mob2 <- rbind(cbind(mob,t(mob)),cbind(t(mob),mob))
summary(ca(mob2))

# Chapter 23: Data Recoding

# Read in original data from file 'EU.txt' into data frame EU

EUr <- apply(EU,2,rank)-1
EUd <- cbind(EUr,11-EUr)
colnames(EUd) <- c(paste(colnames(EU), "-", sep = ""), 
                   paste(colnames(EU), "+", sep = ""))
EUd

summary(ca(EUd))

plot(ca(EUd))

# Chapter 24: Canonical Correspondence Analysis

# Install and load the vegan package, with function cca()
# Read from file 'benthos.xls' the biological data into data
# frame bio, and the environmental data into data frame env

summary(cca(bio,env))

# Chapter 25: Aspects of Stability and Inference

# Assuming ca package installed load data set author

data(author)
author.ca <- ca(author)
nsim      <- 1000       # number of simulations

# compute row sums

author.rowsum <- apply(author, 1, sum)

# compute nsim simulations of first book

author.sim <- rmultinom(nsim, author.rowsum[1], prob = author[1,])

# compute nsim simulations of other books and column-bind

for (i in 2:12) {
   author.sim <- cbind(author.sim, rmultinom(nsim, author.rowsum[i], 
                                             prob = author[i,]))
   }

# transpose to have same format as original matrix

author.sim  <- t(author.sim)
author.sim2 <- matrix(rep(0,12*nsim*26), nrow = 12*nsim)

# reorganize rows so that matrices are together

for (k in 1:nsim) { 
   for (i in 1:12) {
      author.sim2[(k-1)*12+i,] <- author.sim[k+(i-1)*nsim,]
      }
   }

# get standard coordinates for rows

author.rowsc <- author.ca$rowcoord[,1:2]

# calculate principal coordinates of all replicates using transition formula

author.colsim <- t(t(author.rowsc)%*%author.sim2[1:12,])/
                   apply(author.sim2[1:12,], 2, sum)
for (k in 2:nsim) author.colsim <- rbind(author.colsim,
       t(t(author.rowsc)%*%author.sim2[((k-1)*12+1):(k*12),])/
         apply(author.sim2[((k-1)*12+1):(k*12),], 2, sum))
 

# reorganize rows of coordinates so that letters are together

author.colsim2 <- matrix(rep(0,26*nsim*2), nrow = 26*nsim)
for (j in 1:26) { 
   for (k in 1:nsim) {
      author.colsim2[(j-1)*nsim+k,] <- author.colsim[j+(k-1)*26,]
      }
   }

# plot all points (use first format of coords for labelling...)

plot(author.colsim[,1], -author.colsim[,2], xlab = "dim1", ylab = "dim2", 
     type = "n")
letters <- c("a","b","c","d","e","f","g","h","i","j","k","l","m","n","o",
             "p","q","r","s","t","u","v","w","x","y","z")
text(author.colsim[,1], -author.colsim[,2], letters, cex = 0.8, col = "gray")

# plot convex hulls for each letter
# first calculate principal coordinates of letters for original matrix

author.col <- t(t(author.rowsc)%*%author)/apply(author, 2, sum)
for (j in 1:26) {
   points <- author.colsim2[(nsim*(j-1)+1):(nsim*j),]

  # note we are reversing second coordinate in all these plots

   points[,2] <- -points[,2]
   hpts       <- chull(points)
   hpts       <- c(hpts,hpts[1])
   lines(points[hpts,], lty = 3, lwd = 2)
   text(author.col[j,1], -author.col[j,2], letters[j], font = 2, cex = 1.5)
   }

# peeling a new convex hull (here we peel until just under 5% are removed)

plot(author.colsim2[,1], -author.colsim2[,2], xlab = "dim1", ylab = "dim2", 
     type = "n")
for (j in 1:26) {
   points <- author.colsim2[(nsim*(j-1)+1):(nsim*j),]

   # note we are reversing second coordinate in all these plots

   points[,2] <- -points[,2]
   repeat {
     hpts <- chull(points)
     npts <- nrow(points[-hpts,])   #next number of points in peeled set
     if(npts/nsim < 0.95) break
     points <- points[-hpts,]
     }
   hpts <- c(hpts,hpts[1])
   lines(points[hpts,], lty = 3)
   text(author.col[j,1], -author.col[j,2], letters[j], font = 2)
   }

# confidence ellipses - needs package ellipse (download from CRAN)

plot(author.colsim2[,1], -author.colsim2[,2], xlab = "dim1", ylab = "dim2", 
     type = "n")
for (j in 1:26) {
   points <- author.colsim2[(nsim*(j-1)+1):(nsim*j),]

   # note we are reversing second coordinate in all these plots

   points[,2] <- -points[,2]
   covpoints  <- cov(points)
   meanpoints <- apply(points,2,mean)
   lines(ellipse(covpoints, centre = meanpoints))
   text(author.col[j,1], -author.col[j,2], letters[j], font = 2)
   }

# Delta method, using SPSS output contained in the file
# 'AUTHORDELTA.txt'. This is the file contents:

.007 .008  .000 
.018 .025 -.043 
.019 .024 -.275 
.015 .016  .393 
.009 .014 -.008 
.021 .024 -.041 
.019 .027 -.016 
.012 .021  .066 
.013 .013  .087 
.077 .087  .023 
.027 .033 -.173 
.013 .016 -.004 
.021 .028  .097 
.011 .014  .043 
.011 .020 -.004 
.029 .035 -.080 
.078 .094 -.001 
.015 .017 -.065 
.013 .024  .046 
.006 .012 -.010 
.016 .021  .029 
.032 .040  .027 
.015 .022  .073 
.082 .087  .116 
.039 .027 -.033 
.091 .099 -.075 

# plot all points (use first format of coords for labelling...)

plot(author.colsim[,1], -author.colsim[,2], xlab = "dim1", ylab = "dim2", 
     type = "n")
letters <- c("a","b","c","d","e","f","g","h","i","j","k","l","m","n","o",
             "p","q","r","s","t","u","v","w","x","y","z")

# calculate principal coordinates of column points

author.princol <- t(t(author.rowsc)%*%author)/apply(author, 2, sum)
text(author.princol[,1], -author.princol[,2], letters,font = 2)

# input covariances from SPSS output

covdelta <- read.table("AUTHORDELTA.txt")
for (j in 1:26) {
  covpoints[1,1] <- covdelta[j,1]^2
  covpoints[2,2] <- covdelta[j,2]^2
  covpoints[1,2] <- covdelta[j,3]*covdelta[j,1]*covdelta[j,2]
  covpoints[2,1] <- covpoints[1,2]
  lines(ellipse(covpoints, centre = c(author.princol[j,1], 
                                      -author.princol[j,2])))
  }

# Permutation test of the positions of the books

# calculate the principal coordinates of the books (note: in the
# present version of ca() the standard coordinates only are returned
# and you have to calculate the principal coordinates by multiplying
# the standard coordinates by the respective singular values (square
# roots of eigenvalues, or principal inertias)

pcs <- author.ca$rowcoord%*%diag(author.ca$sv)

# check the order of the books

author.ca$rownames

# check re-ordering in order to bring pairs together

author.ca$rownames[c(1,4,2,11,3,8,5,9,6,7,10,12)]

# reorder the rows of pcs[,] so that pairs of books come together

pcs <- pcs[c(1,4,2,11,3,8,5,9,6,7,10,12),]

# calculate the test-statistic (sum of distances between pairs, 
# calculated in K dimensions)

K    <- 2
stat <- 0
for(i in seq(1, 11, 2)){ 
  dist2 <- 0
  for(k in 1:K) {
    dist2 <- dist2+(pcs[i,k]-pcs[i+1,k])^2
    }
  stat <- stat+sqrt(dist2)
  } 

# generate all permutations (exact null distribution)

perm     <- rep(0, 12)
nperm    <- 0
perm[1]  <- 1
statnull <- rep(0, 10395)
for(j11 in 2:12){
  perm[2] <- j11
  rem11   <- c(1:12)[-perm[1:2]]
  perm[3] <- rem11[1]
  for(j9 in 2:10){
    perm[4] <- rem11[j9]
    rem9    <- c(1:12)[-perm[1:4]]
    perm[5] <- rem9[1]
    for(j7 in 2:8){
      perm[6] <- rem9[j7]
      rem7    <- c(1:12)[-perm[1:6]]
      perm[7] <- rem7[1]
      for(j5 in 2:6){
        perm[8] <- rem7[j5]
        rem5    <- c(1:12)[-perm[1:8]]
        perm[9] <- rem5[1]
        for(j3 in 2:4){
          perm[10] <- rem5[j3]
          rem3     <- c(1:12)[-perm[1:10]]
          perm[11] <- rem3[1]
          perm[12] <- rem3[2]
          nperm    <- nperm+1
          pcsp     <- pcs[perm,]
          statp    <- 0
          for(i in seq(1,11,2)){ 
            dist2 <- 0
            for(k in 1:K){ 
              dist2 <- dist2+(pcsp[i,k]-pcsp[i+1,k])^2
              }
            statp <- statp+sqrt(dist2)
            } 
          statnull[nperm]<-statp
          }
        }
      }
    }
  }


stat
statsort <- sort(statnull)
statsort[1:100]
hist(statnull, breaks = 50, xlim = c(0.4,1.2))

# Same for consonants

author.ca <- ca(author,subsetcol=consonants)

# calculate the principal coordinates of the books (note: in the
# present version of ca() the standard coordinates only are returned
# and you have to calcalute the principal coordinates by multiplying
# the standard coordinates by the respective singular values (square
# roots of eigenvalues, or principal inertias)

pcs <- author.ca$rowcoord%*%diag(author.ca$sv)

# check the order of the books

author.ca$rownames

# check re-ordering in order to bring pairs together

author.ca$rownames[c(1,4,2,11,3,8,5,9,6,7,10,12)]

# reorder the rows of pcs[,] so that pairs of books come together

pcs <- pcs[c(1,4,2,11,3,8,5,9,6,7,10,12),]

# calculate the test-statistic (sum of distances between pairs, 
# calculated in K dimensions)

K    <- 2
stat <- 0
for(i in seq(1,11,2)){ 
  dist2 <- 0
  for(k in 1:K){
    dist2 <- dist2+(pcs[i,k]-pcs[i+1,k])^2
    }
  stat  <- stat+sqrt(dist2)
  } 

# generate all permutations (exact null distribution)
# (here repeat above nested loops)

statsort <- sort(statnull)
statsort[1:100]
hist(statnull, breaks = 50, xlim = c(0.4,1.2))

# Same for vowels

author.ca <- ca(author, subsetcol = vowels)

# calculate the principal coordinates of the books (note: in the
# present version of ca() the standard coordinates only are returned
# and you have to calcalute the principal coordinates by multiplying
# the standard coordinates by the respective singular values (square
# roots of eigenvalues, or principal inertias)

pcs <- author.ca$rowcoord%*%diag(author.ca$sv)

# check the order of the books

author.ca$rownames

# check re-ordering in order to bring pairs together

author.ca$rownames[c(1,4,2,11,3,8,5,9,6,7,10,12)]

# reorder the rows of pcs[,] so that pairs of books come together

pcs <- pcs[c(1,4,2,11,3,8,5,9,6,7,10,12),]

# calculate the test-statistic (sum of distances between pairs, 
# calculated in K dimensions)

K    <- 2
stat <- 0
for(i in seq(1,11,2)){
  dist2 <- 0
  for(k in 1:K){
    dist2 <- dist2+(pcs[i,k]-pcs[i+1,k])^2
    }
  stat <- stat+sqrt(dist2)
  } 

# generate all permutations (exact null distribution)
# (here repeat above nested loops)

stat
statsort <- sort(statnull)
statsort[1:100]
hist(statnull, breaks = 50, xlim = c(0.13,0.45))