Sample R-code
[2010/09/30]

General information

This page contains the R commands given in "Biplots in Practice".

You can find the data sets in the section 'data sets'

In this page commands that are given in the Computational Appendix of the book are given in brown, just as in the book.
Results of commands are given in green, just like in the book.
Comments about these commands are given in black.
EXTRA means commands given here (but not in the book) to obtain additional results contained in the corresponding chapter.
Both the comments and commands of EXTRA material are in blue, with results given in green.

R code

# Chapter 1: Biplots - the basic idea

# Read in European indicators data set (Exhibit 1.1) into 
# data frame EU2008 (see p.168)
# For example, on Windows PC's, copy table, and then type:  
# EU2008 <- read.table("clipboard")
windows(width = 11, height = 6)
par(mfrow = c(1,2))
plot(EU2008[,2:1], type = "n", cex.axis = 0.7, xlab = "GDP/capita", 
     ylab = "Purchasing power/capita")
text(EU2008[,2:1], labels = rownames(EU2008), col = "forestgreen", 
     font = 2)
plot(EU2008[,2:3], type = "n", cex.axis = 0.7, xlab = "GDP/capita", 
     ylab = "Inflation rate")
text(EU2008[,2:3], labels = rownames(EU2008), col = "forestgreen", 
     font = 2)

# Note: we use "forestgreen" for the texts above above, more like 
# the green in the book.
# Install the R package rgl for 3-d graphics, and then load 
# as follows:
library(rgl)
plot3d(EU2008[,c(2,1,3)], xlab = "GDP", ylab = "Purchasing power", 
       zlab = "Inflation", font = 2, col = "brown", type = "n")
text3d(EU2008[,c(2,1,3)], text = rownames(EU2008), font = 2, 
       col = "forestgreen") 

# Chapter 2: Regression biplots

# 'bioenv' data set with 8 columns (first 8 columns of Exhibit 2.1) 
# in data frame bioenv. d is species d, y is pollution, x is depth
d <- bioenv[,4]
y <- bioenv[,6]
x <- bioenv[,7]

summary(lm(d ~ y + x))
Call:
lm(formula = d ~ y + x)

Residuals:
     Min       1Q   Median       3Q      Max 
-11.7001  -2.4684   0.1749   3.0563   9.1803 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)  6.13518    6.25721   0.980  0.33554   
y           -1.38766    0.48745  -2.847  0.00834 **
x            0.14822    0.06684   2.217  0.03520 * 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Residual standard error: 5.162 on 27 degrees of freedom
Multiple R-squared: 0.4416,     Adjusted R-squared: 0.4003 
F-statistic: 10.68 on 2 and 27 DF,  p-value: 0.0003831 
# two ways of calculating regression coefficients
# 1. by doing the regression on standardized variables
ds <- (d-mean(d)) / sd(d)
ys <- (y-mean(y)) / sd(y)
xs <- (x-mean(x)) / sd(x)
summary(lm(ds ~ ys + xs))
Call:
lm(formula = ds ~ ys + xs)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.75515 -0.37028  0.02623  0.45847  1.37715 

Coefficients:
              Estimate Std. Error  t value Pr(>|t|)   
(Intercept)  2.487e-17  1.414e-01 1.76e-16  1.00000   
ys          -4.457e-01  1.566e-01   -2.847  0.00834 **
xs           3.472e-01  1.566e-01    2.217  0.03520 * 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Residual standard error: 0.7744 on 27 degrees of freedom
Multiple R-squared: 0.4416,     Adjusted R-squared: 0.4003 
F-statistic: 10.68 on 2 and 27 DF,  p-value: 0.0003831 
# 2. by direct calculation
lm(d ~ y + x)$coefficients[2] * sd(y) / sd(d)
         y 
-0.4457286 
lm(d ~ y + x)$coefficients[3] * sd(x) / sd(d)
        x 
0.3471993 
# finding and storing the standardized regression coefficients for 
# all five species - matrix B in (2.2)
B <- lm(bioenv[,1] ~ y+x)$coefficients[2:3] * c(sd(y),sd(x)) / 
        sd(bioenv[,1])
for(j in 2:5) {
  B <- rbind(B,lm(bioenv[,j] ~ y + x)$coefficients[2:3] * 
         c(sd(y),sd(x)) / sd(bioenv[j]))
  }
B
           y           x
B -0.7171713  0.02465266
  -0.4986038  0.22885450
   0.4910580  0.07424574
  -0.4457286  0.34719935
  -0.4750841 -0.39952072
# EXTRA: storing the variance explained for each regression 
# and computing overall variance explained
var.expl <- var(lm(bioenv[,1] ~ y + x)$fitted.values)
print(var(lm(bioenv[,j] ~ y + x)$fitted.values) / var(bioenv[,j]))
[1] 0.2351773
for(j in 2:5) {
  var.expl <- var.expl + var(lm(bioenv[,j] ~ y + x)$fitted.values)
  print(var(lm(bioenv[,j] ~ y + x)$fitted.values) / var(bioenv[,j]))
  }
[1] 0.3912441
[1] 0.2178098
[1] 0.4416404
[1] 0.2351773
# overall R2
var.expl / sum(apply(bioenv[,1:5], 2, var))
[1] 0.4145226
# EXTRA: doing the same as previous 'EXTRA', but for three predictor 
# variables, adding z = temperature
z <- bioenv[,8]
var.expl3 <- var(lm(bioenv[,1] ~ y + x + z)$fitted.values)
print(var(lm(bioenv[,1] ~ y + x + z)$fitted.values) / var(bioenv[,1]))
[1] 0.5569366
for(j in 2:5) {
  var.expl3 <- var.expl3 + 
                 var(lm(bioenv[,j] ~ y + x + z)$fitted.values)
  print(var(lm(bioenv[,j] ~ y + x + z)$fitted.values) / 
          var(bioenv[,j]))
  }
[1] 0.3924837
[1] 0.2181056
[1] 0.4416436
[1] 0.2463643
# overall R2
var.expl3 / sum(apply(bioenv[,1:5], 2, var))
[1] 0.4271042
# Drawing a regression biplot (Exhibit 2.5) using standardized values 
# of x and y.
# species labels put at the point, arrows 5% shorter
# Note: either close the previous graphics window, or open a new  
# one -on Windows machines open with command win.graph()
win.graph()
plot(xs, ys, xlab = "x*(depth)", ylab = "y*(pollution)", type = "n", 
     asp = 1, cex.axis = 0.7)
text(xs, ys, labels = rownames(bioenv), col = "forestgreen")
text(B[,2:1], labels = colnames(bioenv[,1:5]), col = "chocolate", 
     font = 4)
arrows(0, 0, 0.95*B[,2], 0.95*B[,1], col = "chocolate", angle = 15, 
       length = 0.1)

# (Notice that we have varied the color definitions to be more like 
# those in the book)

# Chapter 3: Generalized linear model biplots

# d0 is the fourth-root transformed d; do regression of d0 on 
# ys and xs
d0 <- d^0.25
summary(lm(d0 ~ ys + xs))
Call:
lm(formula = d0 ~ ys + xs)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.50933 -0.04458  0.13950  0.25483  0.85149 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.63908    0.09686  16.923 6.71e-16 ***
ys          -0.28810    0.10726  -2.686   0.0122 *  
xs           0.05959    0.10726   0.556   0.5831    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Residual standard error: 0.5305 on 27 degrees of freedom
Multiple R-squared: 0.2765,     Adjusted R-squared: 0.2229 
F-statistic: 5.159 on 2 and 27 DF,  p-value: 0.01266 
# EXTRA: saving all the coefficients for the double square-root 
# regressions (Exhibit 3.1)
# Also accumulating variance explained in each and computing overall
dsqrt_abcde       <- bioenv[,1:5]^0.25
dsqrt_abcde_coefs <- lm(dsqrt_abcde[,1] ~ ys + xs)$coefficients
dsqrt.var.expl    <- var(lm(dsqrt_abcde[,1] ~ ys + xs)$fitted.values)
print(var(lm(dsqrt_abcde[,j] ~ ys + xs)$fitted.values) / 
        var(dsqrt_abcde[,j]))
[1] 0.2284871
for(j in 2:5){
  dsqrt_abcde_coefs <- rbind(dsqrt_abcde_coefs,
                           lm(dsqrt_abcde[,j] ~ ys + xs)$coefficients)
  dsqrt.var.expl <- dsqrt.var.expl + 
                      var(lm(dsqrt_abcde[,j] ~ ys + xs)$fitted.values)
  print(var(lm(dsqrt_abcde[,j] ~ ys + xs)$fitted.values) / 
          var(dsqrt_abcde[,j]))
  }
[1] 0.3621872
[1] 0.1591791
[1] 0.2764936
[1] 0.2284871
# overall R2
dsqrt.var.expl / sum(apply(dsqrt_abcde[,1:5], 2, var))
[1] 0.3393836
# EXTRA: drawing the biplot of Exhibit 3.2
plot(xs, ys, xlab = "x*(depth)", ylab = "y*(pollution)", type = "n", 
     asp = 1, cex.axis = 0.7)
text(xs, ys, labels = rownames(bioenv), col = "forestgreen")
text(dsqrt_abcde_coefs[,3:2], labels = colnames(bioenv[,1:5]), 
     col = "chocolate", font = 4)
for(j in 1:5){ 
  arrows(0, 0, 0.95*dsqrt_abcde_coefs[j,3], 
         0.95*dsqrt_abcde_coefs[j,2], col = "chocolate", angle = 15, 
         length = 0.1)
  }

# performing Poisson regression for species d:
summary(glm(d ~ ys + xs, family = poisson))
Call:
glm(formula = d ~ ys + xs, family = poisson)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-4.4123  -0.7585   0.1637   0.8262   2.7565  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  2.29617    0.06068  37.838  < 2e-16 ***
ys          -0.33682    0.07357  -4.578 4.69e-06 ***
xs           0.19963    0.06278   3.180  0.00147 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 144.450  on 29  degrees of freedom
Residual deviance:  88.671  on 27  degrees of freedom
AIC: 208.55

Number of Fisher Scoring iterations: 5
# to get error deviance for Poisson regression of species d:
poisson.glm <- glm(d~ys+xs, family=poisson)
poisson.glm$deviance/poisson.glm$null.deviance
[1] 0.6138564
# EXTRA: saving all the coefficients for the Poisson regressions 
# (Exhibit 3.4)
pois_abcde_coefs <- glm(abcde[,1] ~ ys + xs, 
                        family = poisson)$coefficients
for(j in 2:5){ 
  pois_abcde_coefs <- rbind(pois_abcde_coefs, 
              glm(abcde[,j] ~ ys + xs, family=poisson)$coefficients)
  }
pois_abcde_coefs
                 (Intercept)         ys          xs
pois_abcde_coefs   2.1786692 -1.1254080 -0.06714289
                   1.8525985 -0.8122742  0.18286062
                   2.0411089  0.4174381  0.05321664
                   2.2961652 -0.3368244  0.19963421
                   0.8277491 -0.8234332 -0.56805285
# performing logistic regression for species d, after converting to 
# presence/absence:
d01 <- d > 0
summary(glm(d01 ~ ys + xs, family = binomial))
Call:
glm(formula = d01 ~ ys + xs, family = binomial)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.2749   0.2139   0.2913   0.3884   1.2874  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept)   2.7124     0.8533   3.179  0.00148 **
ys           -1.1773     0.6522  -1.805  0.07105 . 
xs           -0.1369     0.7097  -0.193  0.84708   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 19.505  on 29  degrees of freedom
Residual deviance: 15.563  on 27  degrees of freedom
AIC: 21.563

Number of Fisher Scoring iterations: 6
# to get error deviance for logistic regression of species d:
logistic.glm <- glm(d01 ~ ys + xs, family = binomial)
logistic.glm$deviance / logistic.glm$null.deviance
[1] 0.7979165
# EXTRA: saving all the coefficients for the logistic regressions 
# (Exhibit 3.5)
abcde <- bioenv[,1:5]
abcde01 <- abcde > 0
logit_abcde_coefs <- glm(abcde01[,1] ~ ys + xs, 
                         family = binomial)$coefficients
for(j in 2:5){ 
  logit_abcde_coefs <- rbind(logit_abcde_coefs, 
        glm(abcde01[,j] ~ ys + xs, family = binomial)$coefficients)
  }
logit_abcde_coefs
                  (Intercept)         ys         xs
logit_abcde_coefs   2.3842273 -2.8885410  0.8631339
                    1.2731575 -1.4183387 -0.1431720
                    0.8309807  0.9726882  0.3148828
                    2.7124406 -1.1772989 -0.1368686
                    0.2532229 -1.2796484 -0.7858995
# EXTRA: another example of a plot (Example 3.6)
# This time the arrows are drawn to the exact point and labels 
# are offset
plot(xs, ys, type = "n", xlab = "x* (depth)", ylab = "y* (pollution)", 
     ylim = c(-3,2.6), cex.axis = 0.7, asp = 1)
text(xs, ys, labels = rownames(bioenv), col = "forestgreen", cex=0.8)
arrows(0, 0, logit_abcde_coefs[,3], logit_abcde_coefs[,2], 
       col = "chocolate", length = 0.1, angle = 15, lwd = 1)
text(1.05*logit_abcde_coefs[,3:2], labels = colnames(bioenv[,1:5]), 
     col = "chocolate", font = 4, cex = 1.1)

# Chapter 4: Multidimensional scaling biplots

# 'countries' data set -- dissimilarities between 13 countries 
# (Exhibit 4.1) -- 
# assumed to have been read into data frame MT_matrix 
# R function cmdscale performs classical multidimensional scaling
# Exhibit 4.2 obtained as follows
plot(cmdscale(MT_matrix), type = "n", xlab = "dim 1", ylab = "dim 2", 
     asp = 1)
text(cmdscale(MT_matrix), labels = colnames(MT_matrix), 
     col = "forestgreen")

# EXTRA (*only* for MDS experts -- there is a bug in the cmdscale 
# function... read on...)
# The eigenvalues of the MDS (notice that warnings are given that 
# some are negative)
cmdscale(MT_matrix, eig = TRUE, k = 1)$eig

# Function cmdscale does not allow obtaining the last (13th) 
# eigenvalue, and there is thus a bug in the computation of the 
# goodness-of-fit (GOF), which doesn't take last dimension into 
# account
cmdscale(MT_matrix, eig = TRUE)$GOF
[1] 0.6177913 0.6500204         
# In the above the first reported GOF value is wrong. It is based on 
# the first 12 eigenvalues:
sum(abs(cmdscale(MT_matrix, eig = TRUE, k = 12)$eig))
[1] 264.5464 
# (it is not possible to set k=13)
# Here is the correct total, which includes the last (negative) 
# eigenvalue
MT_sum <- sum(abs(eigen(-0.5*cmdscale(MT_matrix, eig = TRUE, k = 12, 
           x.ret = TRUE)$x)$values))
MT_sum
[1] 288.0903
# The following verifiess the second GOF value, which is correct:
cmdscale(MT_matrix, eig = TRUE, k = 12)$eig[1:2] / 
  sum(cmdscale(MT_matrix, eig = TRUE, k = 12)$eig[1:8])
[1] 0.3821252 0.2678952
# (0.6500204 = 0.3821252 + 0.2678952)
# And the following revised calculation gives what the first GOF 
# should have been:
cmdscale(MT_matrix, eig = TRUE, k = 12)$eig[1:2] / MT_sum
[1] 0.3334984 0.2338045
# that is, the sum for the two-dimensional solution is not 0.6177913, 
# but rather
sum(cmdscale(MT_matrix, eig = TRUE, k = 12)$eig[1:2] / MT_sum)
[1] 0.5673029
# 'atributes' data set -- ratings of countries on 6 attributes 
# (Exhibit 4.3) -- 
# assumed to have been read into data frame MT_ratings

# keep coordinates of MDS plot in MT_dims
MT_dims <- cmdscale(MT_matrix, eig = TRUE, k = 2)$points
colnames(MT_dims) <- c("dim1","dim2")

# calculate the regression coefficients of the 6 attributes on the 
# coordinates (Exhibit 4.4)
MT_coefs <- lm(MT_ratings[,1] ~ MT_dims[,1] + 
                                MT_dims[,2])$coefficients
for(j in 2:6){ 
  MT_coefs <- rbind(MT_coefs, MT_coefs <- 
          lm(MT_ratings[,j] ~ MT_dims[,1] + MT_dims[,2])$coefficients)
  }

# plot the regression coefficients on the MDS plot (Exhibit 4.5)
plot(cmdscale(MT_matrix), type = "n", xlab = "dim 1", ylab = "dim 2", 
     asp = 1)
text(cmdscale(MT_matrix), labels = colnames(MT_matrix), 
     col = "forestgreen")
arrows(0, 0, MT_coefs[,2], MT_coefs[,3], col = "chocolate", 
       length = 0.1, angle = 10)
text(1.2*MT_coefs[,2:3], labels = colnames(MT_ratings), 
     col = "chocolate", font = 4, cex = 0.8)

# Function to calculate chi-square distances between row or column
# profiles of a matrix
# Usage: chidist(N,1) calculates the chi-square distances between row 
# profiles
#        (for row profiles, chidist(N) is sufficient)
#        chidist(N,2) calculates the chi-square distances between 
#        column profiles
chidist <- function(mat, rowcol = 1) {
  if(rowcol == 1) {
    prof <- mat / apply(mat, 1, sum)
    rootaveprof <- sqrt(apply(mat, 2, sum) / sum(mat))
    }
  if(rowcol == 2) {
    prof <- t(mat) / apply(mat, 2, sum)
    rootaveprof <- sqrt(apply(mat, 1, sum) / sum(mat))
    }
  dist(scale(prof, FALSE, rootaveprof))  
  }

# Note: functions defined in this way are added as objects to your 
# workspace, and are saved if the workspace is saved, otherwise they
# are lost

# Compute chi-square distances between row profiles, perform 
# classical MDS and plot
abcde <- bioenv[,1:5]
abcde_mds <- cmdscale(chidist(abcde), eig = TRUE, k = 4)
100 * abcde_mds$eig / sum(abcde_mds$eig)
[1] 52.364447 22.014980 16.187784  9.432789
plot(abcde_mds$points[,1:2], type = "n", xlab = "dim 1", 
     ylab = "dim 2", xlim = c(-1.2,1.6), ylim = c(-1.1,1.8), asp = 1) 
text(abcde_mds$points[,1:2], labels = rownames(abcde), 
     col = "forestgreen") 

# Now add columns (species) as biplot vectors
# First convert to row profiles and standardize in the chi-squared way
abcde_prof <- abcde / apply(abcde, 1, sum)
abcde_prof_stand <- t(t(abcde_prof) / sqrt(apply(abcde, 2, sum) / 
                                           sum(abcde)))

# perform regressions and save coefficients for plotting
mds_coefs <- lm(abcde_prof_stand[,1] ~ abcde_mds$points[,1] + 
                abcde_mds$points[,2])$coefficients
for(j in 2:5){
  mds_coefs <- rbind(mds_coefs, MT_coefs <- 
                 lm(abcde_prof_stand[,j] ~ abcde_mds$points[,1] + 
                    abcde_mds$points[,2])$coefficients)
  }
arrows(0, 0, mds_coefs[,2], mds_coefs[,3], col = "chocolate", 
       length = 0.1, angle = 10)
text(1.1*mds_coefs[,2:3], labels = colnames(abcde), col = "chocolate", 
     font = 4, cex = 0.9)

# Assume that the sediment vector has been read into a character 
# vector (e.g., in Windows, copy the vector from the Excel file, 
# without the column label, and then type in the command:
#   sediment <- scan("clipboard", what="character")   

# convert sediment (consisting of "C","G" and "S") to factor object
sediment <- as.factor(sediment)

# compute mean positions in biplot and plot them
sediment.means <- cbind(tapply(abcde_mds$points[,1], sediment, mean), 
                        tapply(abcde_mds$points[,2], sediment, mean))
text(sediment.means[,1], sediment.means[,2], labels = c("C","G","S"), 
     font = 4, cex = 0.8, col = "chocolate")

# convert sediment categories to dummy variables
clay01   <- sediment == "C"
gravel01 <- sediment == "G"
sand01   <- sediment == "S"

# compute logistic regression coefficients and plot them connected 
# to origin
sediment_coefs <- glm(as.numeric(clay01) ~ abcde_mds$points[,1] + 
    abcde_mds$points[,2], family = "binomial")$coefficients
sediment_coefs <- rbind(sediment_coefs, 
    glm(as.numeric(gravel01) ~ abcde_mds$points[,1] + 
        abcde_mds$points[,2], family = "binomial")$coefficients)
sediment_coefs <- rbind(sediment_coefs, 
    glm(as.numeric(sand01) ~ abcde_mds$points[,1] + 
        abcde_mds$points[,2], family = "binomial")$coefficients)
segments(0, 0, sediment_coefs[,2], sediment_coefs[,3], 
         col = "chocolate")
text(sediment_coefs[,2:3], labels = c("C","G","S"), col = "chocolate", 
     font = 4)

# Chapter 5: Reduced-dimension biplots

# Small example from Chapter 1, biplotted
# (Note that the last two singular values, theoretically zero, and 
# their respective singular vectors, can turn out slightly differently 
# in different versions of R)
Y <- matrix(c(8,5,-2,2,4,2,0,-3,3,6,2,3,3,-3,-6,-6,-4,1,-1,-2), 
            nrow = 5)
svd(Y)
$d
[1] 1.412505e+01 9.822577e+00 6.351831e-16 3.592426e-33

$u
           [,1]       [,2]        [,3]          [,4]
[1,] -0.6634255 -0.4574027 -0.59215653  2.640623e-35
[2,] -0.3641420 -0.4939878  0.78954203  2.167265e-34
[3,]  0.2668543 -0.3018716 -0.06579517 -9.128709e-01
[4,] -0.2668543  0.3018716  0.06579517 -1.825742e-01
[5,] -0.5337085  0.6037432  0.13159034 -3.651484e-01

$v
           [,1]       [,2]       [,3]       [,4]
[1,] -0.7313508 -0.2551980 -0.6276102 -0.0781372
[2,] -0.4339970  0.4600507  0.2264451  0.7407581
[3,]  0.1687853 -0.7971898  0.0556340  0.5769791
[4,]  0.4982812  0.2961685 -0.7427873  0.3350628
colnames(Y) <- c("A","B","C","D")
rowcoord <- svd(Y)$u %*% diag(sqrt(svd(Y)$d))
colcoord <- svd(Y)$v %*% diag(sqrt(svd(Y)$d))
plot(rbind(rowcoord,colcoord)[,1:2], type = "n", xlab = "", ylab = "", 
     asp = 1, cex.axis = 0.7)
abline(h = 0, v = 0, lty = "dotted")
text(rowcoord[,1:2], labels = 1:5, col = "forestgreen")
text(colcoord[,1:2], labels = colnames(Y), col = "chocolate", 
     font = 3)

# Chapter 6: Principal component biplots

# 'attributes' data set already in data frame MT_ratings (Chap. 4)

# centre the data by the column means
MT_means <- apply(MT_ratings, 2, mean)
MT_Y     <- sweep(MT_ratings, 2, MT_means)

# compute matrix with equal row and column weights and compute SVD
MT_Y   <- MT_Y/sqrt(nrow(MT_Y)*ncol(MT_Y))
MT_SVD <- svd(MT_Y)

# form biplot: row principal, column standard coordinates
MT_F <- sqrt(nrow(MT_Y))*MT_SVD$u%*%diag(MT_SVD$d)
MT_G <- sqrt(ncol(MT_Y))*MT_SVD$v
plot(rbind(MT_F[,1:2],MT_G[,1:2]), type = "n", asp = 1, 
     xlim = c(-3.6,2.3), xlab = "dim 1", ylab = "dim 2", 
     cex.axis = 0.7) 
text(MT_F[,1:2], labels = rownames(MT_ratings), col = "forestgreen", 
     cex = 0.8) 
arrows(0, 0, MT_G[,1], MT_G[,2], col = "chocolate", length = 0.1, 
       angle = 10)
text(c(1.07,1.3,1.07,1.35,1.2,1.4)*MT_G[,1], 
     c(1.07,1.07,1.05,1,1.16,1.1)*MT_G[,2], 
     labels = colnames(MT_ratings), col = "chocolate", font = 4, 
     cex = 0.8)

# covariance biplot: row standard, column principal coordinates
MT_F <- sqrt(nrow(MT_Y))*MT_SVD$u
MT_G <- sqrt(ncol(MT_Y))*MT_SVD$v%*%diag(MT_SVD$d)
plot(rbind(MT_F[,1:2],MT_G[,1:2]), type = "n", asp = 1, 
     xlim = c(-3.6, 2.3), xlab = "dim 1", ylab = "dim 2", 
     cex.axis = 0.7) 
text(MT_F[,1:2], labels = rownames(MT_ratings), col = "forestgreen", 
     cex = 0.8) 
arrows(0, 0, MT_G[,1], MT_G[,2], col = "chocolate", length = 0.1, 
       angle = 10)
text(c(1.07,1.2,1.07,1.25,1.07,1.3)*MT_G[,1], 
     c(1.07,1.07,1.04,1.02,1.16,1.07)*MT_G[,2], 
     labels = colnames(MT_ratings), col = "chocolate", font = 4, 
     cex = 0.8)

# scree plot
MT_percents <- 100*MT_SVD$d^2/sum(MT_SVD$d^2)
MT_percents <- MT_percents[seq(6,1)]
barplot(MT_percents, horiz = TRUE, cex.axis = 0.7)

# Chapter 7: Log-ratio biplots

# get the US Arrests data provided in R
data(USArrests)

# perform (weighted) log-ratio analysis on columns 1, 2 and 4
# rm and cm are the row and column margins; 
# mr and mc are the weighted means used in the double-centering
N    <- USArrests[,c(1,2,4)]
P    <- N / sum(N)
rm   <- apply(P, 1, sum)
cm   <- apply(P, 2, sum)
Y    <- as.matrix(log(P))
mc   <- t(Y) %*% as.vector(rm)
Y    <- Y - rep(1, nrow(P)) %*% t(mc)
mr   <- Y %*% as.vector(cm)
Y    <- Y - mr %*% t(rep(1,ncol(P)))
Z    <- diag(sqrt(rm)) %*% Y %*% diag(sqrt(cm))
svdZ <- svd(Z)

# compute form biplot coordinates from results of SVD
USA_F <- diag(1/sqrt(rm)) %*% svdZ$u[,1:2] %*% diag(svdZ$d[1:2])
USA_G <- diag(1/sqrt(cm)) %*% svdZ$v[,1:2]

# biplot - axes with different scales plotted individually
plot(rbind(USA_F, USA_G/20), xlim = c(-0.35,0.45), 
     ylim = c(-0.18,0.23), asp = 1, type = "n", xlab = "", ylab = "", 
     xaxt = "n", yaxt = "n", cex.axis = 0.7)
axis(1, cex.axis = 0.7, col.axis = "forestgreen", 
     col.ticks = "forestgreen")
axis(2, cex.axis = 0.7, col.axis = "forestgreen", 
     col.ticks = "forestgreen", at = seq(-0.2, 0.2, 0.2))
axis(3, cex.axis = 0.7, col.axis = "chocolate", 
     col.ticks = "chocolate", at = seq(-0.4, 0.4, 0.2), 
     labels = seq(-8,8,4))
axis(4, cex.axis = 0.7, col.axis = "chocolate", 
     col.ticks = "chocolate", at = seq(-0.2, 0.2, 0.2), 
     labels = seq(-4,4,4))
text(USA_F, labels = rownames(N), col = "forestgreen", font = 2, 
     cex = 0.8)
text(USA_G/20, labels = colnames(N), col = "chocolate", font = 4, 
     cex = 0.9)

# total variance (inertia): either sum of squares of Z matrix, or sum 
# of squared singular values
sum(Z*Z)
[1] 0.01790182
sum(svdZ$d^2)
[1] 0.01790182
# fish morphology example: assume data read into data frame fish:
#   -- first two columns are sex, habitat
#   -- remaining 26 columns are the morphometric data
fish.morph<-fish[,3:ncol(fish)]

# log-ratio analysis is performed in the same way as above
N    <- fish.morph
P    <- N/sum(N)
rm   <- apply(P, 1, sum)
cm   <- apply(P, 2, sum)
Y    <- as.matrix(log(P))
mc   <- t(Y) %*% as.vector(rm)
Y    <- Y - rep(1,nrow(P)) %*% t(mc)
mr   <- Y %*% as.vector(cm)
Y    <- Y - mr %*% t(rep(1,ncol(P)))
Z    <- diag(sqrt(rm)) %*% Y %*% diag(sqrt(cm))
svdZ <- svd(Z)

# total variance
sum(Z*Z)
[1] 0.001960883
# calculate sex-habitat interaction and set up labels for four 
# sex*habitat groups
fish.sexhab                   <- 2*(fish[,2]-1)+fish[,1]
fish.labels                   <- rep("fL", nrow(fish))
fish.labels[fish.sexhab=="2"] <- "mL"
fish.labels[fish.sexhab=="3"] <- "fP"
fish.labels[fish.sexhab=="4"] <- "mP"

# EXTRA: log-ratio biplot of morphological data 
# Exhibit 7.3 (first axis coordinates reversed)
fish_F     <- diag(1/sqrt(rm)) %*% svdZ$u %*% diag(svdZ$d)
fish_G     <- diag(1/sqrt(cm)) %*% svdZ$v
fish_F[,1] <- -fish_F[,1]
fish_G[,1] <- -fish_G[,1]
plot(rbind(fish_F,fish_G/50), type = "n", xlab = "", ylab = "", 
     asp = 1, xaxt = "n", yaxt = "n", cex.axis = 0.7)
axis(1, col.axis = "forestgreen", col.ticks = "forestgreen")
axis(2, col.axis = "forestgreen", col.ticks = "forestgreen")
axis(3, col.axis = "chocolate", col.ticks = "chocolate", 
     at = seq(-0.06, 0.06, 0.02), labels = seq(-3,3,1))
axis(4, col.axis = "chocolate", col.ticks = "chocolate", 
     at = seq(-0.04, 0.04, 0.02), labels = seq(-2, 2, 1))
text(fish_F, labels = fish.labels, col = "forestgreen", font = 2, 
     cex = 0.7)
text(fish_G/50, labels = colnames(fish.morph), col = "chocolate", 
     font = 4, cex=0.8)

# scatterplot of log-ratios in Exhibit 7.4
logFdlFal <- log(fish.morph[,"Fdl"] / fish.morph[,"Fal"])	
logFdwFal <- log(fish.morph[,"Fdw"] / fish.morph[,"Fal"])
plot(logFdlFal, logFdwFal, asp = 1, pch = 24, cex = 0.7, 
     xlab = "log(Fdl/Fal)", ylab = "log(Fdw/Fal)", col = "forestgreen")
abline(a = 0.0107, b = 0.707, lty = 2, col = "chocolate")

# predictions of Fdw in Exhibit 7.5 
Fdw_pred <- 1.0108*fish.morph[,"Fdl"]^0.707*fish.morph[,"Fal"]^0.293
plot(Fdw_pred, fish.morph[,"Fdw"], xlim = c(18,30), ylim = c(18,30), 
     pch = 24, cex = 0.7, xlab = "predicted Fdw", ylab = "actual Fdw", 
     col = "forestgreen")
abline(a = 0, b = 1, lty = 2, col = "chocolate")

# correlation between predicted and observed
cor(Fdw_pred, fish.morph[,"Fdw"])
[1] 0.7496034

# Chapter 8: Correspondence analysis biplots

# Data set 'smoking' is in the ca package (to be installed from CRAN 
# first, before loading)
# Its name is actually 'smoke' (Exhibit 8.1)
library(ca)
data(smoke)

# computation of CA
N    <- smoke
P    <- N/sum(N)
rm   <- apply(P, 1, sum)
cm   <- apply(P, 2, sum)
Dr   <- diag(rm)
Dc   <- diag(cm)
Z    <- diag(sqrt(1/rm))%*%(as.matrix(P)-rm%*%t(cm))%*%diag(sqrt(1/cm))
svdZ <- svd(Z)

# principal coordinates of rows, standard coordinates of columns
smoke_F <- diag(1/sqrt(rm))%*%svdZ$u %*%diag(svdZ$d)
smoke_G <- diag(1/sqrt(cm))%*%svdZ$v

# plotting (notice that by default the first two columns are used)
plot(rbind(smoke_F, smoke_G), type = "n", xlab = "Dim 1", 
     ylab = "Dim 2", asp = 1)
text(smoke_F, labels = rownames(smoke), col = "forestgreen", font = 2, 
     cex = 0.7)
text(smoke_G, labels = colnames(smoke), col = "chocolate", font = 2, 
     cex = 0.8)

# obtaining Exhibit 8.2 using ca package
plot(ca(smoke), map = "rowprincipal", 
     col = c("forestgreen","chocolate"))

# Summary of CA, with contributions to inertia
summary(ca(smoke))
Principal inertias (eigenvalues):

 dim    value      %   cum%   scree plot               
 1      0.074759  87.8  87.8  *************************
 2      0.010017  11.8  99.5  ***                      
 3      0.000414   0.5 100.0                           
        -------- -----                                 
 Total: 0.085190 100.0                                 


Rows:
    name   mass  qlt  inr    k=1 cor ctr    k=2 cor ctr  
1 |   SM |   57  893   31 |  -66  92   3 | -194 800 214 |
2 |   JM |   93  991  139 |  259 526  84 | -243 465 551 |
3 |   SE |  264 1000  450 | -381 999 512 |  -11   1   3 |
4 |   JE |  456 1000  308 |  233 942 331 |   58  58 152 |
5 |   SC |  130  999   71 | -201 865  70 |   79 133  81 |

Columns:
    name   mass  qlt  inr    k=1 cor ctr    k=2 cor ctr  
1 | none |  316 1000  577 | -393 994 654 |  -30   6  29 |
2 | lght |  233  984   83 |   99 327  31 |  141 657 463 |
3 | medm |  321  983  148 |  196 982 166 |    7   1   2 |
4 | hevy |  130  995  192 |  294 684 150 | -198 310 506 |

# suppose the 'benthos' data have been read into data frame benthos
# calculation of CA and the row contributions (ca package stores 
# standard coordinates in $rowcoord)
benthos.ca     <- ca(benthos)
benthos.F      <- benthos.ca$rowcoord %*% diag(benthos.ca$sv)
benthos.rowcon <- benthos.ca$rowmass * (benthos.F[,1]^2 + 
                    benthos.F[,2]^2) / sum(benthos.ca$sv[1:2]^2)

# set up a vector of species labels where those with less than 0.01 
# contribution are labelled "á"
benthos.names                        <- rownames(benthos)
benthos.names[benthos.rowcon < 0.01] <- "á"

# nonlinear transformation of contributions for plotting of labels
benthos.rowsize <- log(1+exp(1)*benthos.rowcon^0.3)
benthos.rowsize[benthos.rowcon<0.01] <- 1

# plot the asymmetric map with label size related to contribution 
# (Exhibit 8.3)
FF <- benthos.ca$rowcoord
GG <- benthos.ca$colcoord %*% diag(benthos.ca$sv)
plot(c(FF[,1],GG[,1]), c(FF[,2],GG[,2]), type = "n", xlab = "", 
     ylab = "", asp = 1)
text(FF[,1], FF[,2], labels = benthos.names, cex = benthos.rowsize, 
     col = "chocolate", font = 4)
text(GG[,1], GG[,2], labels = colnames(benthos), cex = 0.8, 
     col = "forestgreen", font = 2)

# The contribution biplot can be obtained using the map="colgreen" 
# option in the ca package where "green" refers to Greenacre, not a 
# colour! (The option "colgab" refers to Gabriel...)
# (remember that the columns are the sites, and the species the rows, 
# which is the tranpose of the usual format).

# EXTRA: First substitute the species labels, using blank labels for 
# low contributions
benthos.names[benthos.rowcon < 0.01] <- " "
benthos.ca$rownames <- benthos.names

# Contribution biplot of Exhibit 8.4
plot(benthos.ca, map = "colgreen", mass = c(1,0), 
     col = c("chocolate","forestgreen"), pch = c(17,24,16,1)) 

# draw lines to the species points
segments(0, 0, benthos.ca$rowcoord[,1]*sqrt(benthos.ca$rowmass), 
         benthos.ca$rowcoord[,2]*sqrt(benthos.ca$rowmass), 
         col = "chocolate")

# Chapter 9: Multiple correspondence analysis biplots I

# The original Spanish data (N=2471) is in the Excel file 
# women_Spain2002_original.xls
# Later (in Chapter 10) we will show how to create the concatenated 
# matrix in R.
# Here we use the prepared concatenated file for the reduced sample 
# (N=2107) which has no missing data and where categories H4 and H5 
# have already been combined -- this concatenated matrix (Exhibit 9.1) 
# is in the Excel file women_Spain2002_concat.xls
# So we assume that this concatenated table has been read into the data 
# frame 'women.concat'

# the symmetric CA map, which is the default plotting option of the 
# ca package
plot(ca(women.concat))

# to change the sign of the second axis (if required)
women.ca <- ca(women.concat)
women.ca$rowcoord[,2] <- -women.ca$rowcoord[,2]
women.ca$colcoord[,2] <- -women.ca$colcoord[,2]
plot(women.ca)

# to obtain Exhibit 9.2: plot symbols, compute principal coordinates, 
# add labels and inertias
# NOTE: after plotting, adjust the window to accommodate the points 
# and labels as in Exhibit 9.2. Repeat the plotting again, if 
# necessary, in that window to reposition labels
par(cex.axis = 0.7)
plot(women.ca, labels = 0, col = c("forestgreen", "chocolate"))
women.F <- women.ca$rowcoord %*% diag(women.ca$sv)
women.G <- women.ca$colcoord %*% diag(women.ca$sv)
text(women.F, labels = women.ca$rownames, pos = 4, offset = 0.3, 
     col = "forestgreen", cex = 0.8, font = 2)
text(women.G, labels = women.ca$colnames, pos = 4, offset = 0.3, 
     col = "chocolate", cex = 0.8, font = 4)
text(max(women.G[,1]), 0, "0.0571 (82.1%)", adj = c(0.6,-0.6), 
     cex = 0.7) 
text(0, max(women.G[,2]), "0.0030 (4.4%)", adj = c(-0.1,-3), 
     cex = 0.7)

# asymmetric biplot of Exhibit 9.3
# (again, if you need to adjust the window, repeat plotting in new 
# window to position labels better)
plot(women.ca, map = "rowprincipal", labels = c(0,2), 
     col = c("forestgreen", "chocolate"))

# EXTRA: same as above, but labelling categories in colour
plot(women.ca, map = "rowprincipal", labels = 0, 
     col = c("forestgreen", "chocolate"))
text(women.ca$colcoord, labels = women.ca$colnames, pos = 4, 
     offset = 0.3, col = "chocolate", cex = 0.8, font = 4)

# EXTRA: Exhibit 9.4 -- showing categories 'e5' and 'ma4' connected 
# to their 8 corresponding variables
# computation of coordinates of these demographic categories per 
# substantive variable (the slight complication here is that the last 
# variable has 4 categories, while all the others have 5)
women.e5 <- as.numeric(women.concat["e5",])
women.e5 <- women.e5 / sum(women.e5/8)
women.e5.coord <- apply(women.e5[36:39] * 
                  women.ca$colcoord[36:39,], 2, sum)
for(j in 7:1){
  women.e5.coord <- rbind(apply(women.e5[(5*(j-1)+1):(5*j)] * 
       women.ca$colcoord[(5*(j-1)+1):(5*j),], 2, sum), women.e5.coord)
  }
women.ma4 <- as.numeric(women.concat["ma4",])
women.ma4 <- women.ma4 / sum(women.ma4/8)
women.ma4.coord <- apply(women.ma4[36:39] * 
                         women.ca$colcoord[36:39,], 2, sum)
for(j in 7:1){
  women.ma4.coord <- rbind(apply(women.ma4[(5*(j-1)+1):(5*j)] * 
      women.ca$colcoord[(5*(j-1)+1):(5*j),], 2, sum), women.ma4.coord)
  }

# plot them with principal coordinates of rows and link to their 
# centroid
plot(rbind(women.F, women.e5.coord, women.ma4.coord), type = "n", 
     asp = 1, xlab = "", ylab = "")
points(women.F, pch = 19, col = "forestgreen")
text(women.F, labels = women.ca$rownames, pos = 4, offset = 0.3, 
     col = "forestgreen", cex = 0.9, font = 2)
points(rbind(women.e5.coord, women.ma4.coord), pch = 22, 
       col = "chocolate")
text(rbind(women.e5.coord, women.ma4.coord), labels=LETTERS[1:8], 
     pos = 4, offset = 0.3, col = "chocolate", cex = 0.9, font = 4)
rownames(women.F) <- women.ca$rownames
segments(women.F["e5",1], women.F["e5",2], women.e5.coord[,1], 
         women.e5.coord[,2], col = "chocolate")
segments(women.F["ma4",1], women.F["ma4",2], women.ma4.coord[,1], 
         women.ma4.coord[,2], col = "chocolate")
abline(h = 0, v = 0, lty = 3, col = "gray")

# Exhibit 9.6: contribution biplot using map="rowgreen" option
par(mar=c(3,2,1,1))         # (optional, reduces margins around plot)
plot(women.ca, map  ="rowgreen", mass = c(FALSE,TRUE), 
     col = c("forestgreen", "chocolate"))

# same again, but labelling points in colour after computing their 
# positions and also adding average age and gender groups, as shown in 
# Exhibit 9.6
plot(women.ca, map = "rowgreen", labels = 0, mass = c(FALSE,TRUE), 
     col = c("forestgreen", "chocolate"))
women.sex   <- c(rep("m", 6),rep("f", 6)) 
women.age   <- rep(c("a1","a2","a3","a4","a5","a6"), 2)
women.sex.F <- cbind(tapply(women.F[12:23,1], women.sex, mean), 
                     tapply(women.F[12:23,2], women.sex, mean))
women.age.F <- cbind(tapply(women.F[12:23,1], women.age, mean), 
                     tapply(women.F[12:23,2], women.age, mean))
points(rbind(women.sex.F, women.age.F), pch = 21, col = "forestgreen", 
       cex = 0.7)
text(rbind(women.sex.F, women.age.F), 
     labels = c("f","m","a1","a2","a3","a4","a5","a6"), pos = 4, 
     offset = 0.3, col = "forestgreen", cex = 0.8, font = 2)
text(women.F, labels = women.ca$rownames, pos = 4, offset = 0.3, 
     col = "forestgreen", cex = 0.8, font = 2)
text(sqrt(women.ca$colmass)*women.ca$colcoord, 
     labels = women.ca$colnames, pos = 4, offset = 0.3, 
     col = "chocolate", cex = 0.8, font = 4)

# Chapter 10: Multiple correspondence analysis biplots II

# EXTRA: We make rather a long detour here to show how to 
#    1. read in the original Spanish data (N=2471) and remove cases 
#       with missing data
#    2. compute the corresponding indicator matrix
#    3. compute the Burt matrix and extract the concatenated matrix 
#       used in Chap. 9

# First read in the original data (2471x12 matrix in 
# women_Spain2002_original.xls) into data frame 'women'. We count 
# missing data (=9) for each case, and then keep the rows that have 
# no missing data
dim(women)
[1] 2471   12

women.9 <- rep(0, nrow(women))
for(i in 1:nrow(women)) {
  for(j in 1:12) {
    if(women[i,j] == 9) women.9[i] <- women.9[i]+1
  }
}
women <- women[women.9 == 0,]

dim(women)
[1] 2107   12

# Compute the indicator matrix - there are 8x5 categories for the 
# substantive variables and 2+5+6+6 categories for the demographics, 
# totalling 59 
women.Z <- matrix(0, nrow = nrow(women), ncol = 59)

# We admit the following single statement is cryptic but it works!
# (remember that education categories start at 0)
women.Z[cbind(as.numeric(matrix(rep(1:nrow(women), 12), byrow = TRUE, 
          nrow = 12)), as.numeric(t(women[,1:12]) + 
          c(seq(0,35,5),40,42,48,53)))] <- 1
colnames(women.Z) <- c(sort(paste(c("A","B","C","D","E","F","G","H"),
        rep(1:5, 8), sep = "")),"m","f","ma","wi","di","se","si",
        paste("e", 0:5, sep = ""), paste("a", 1:6, sep = ""))
 
# column sums of women.Z are the frequencies of each category
apply(women.Z,2,sum)
  A1   A2   A3   A4   A5   B1   B2   B3   B4   B5   C1   C2   C3   C4   
 397  932   91  598   89  124  956  237  665  125  165  988  237  620 
  C5   D1   D2   D3   D4   D5   E1   E2   E3   E4   E5   F1   F2   F3
  97   92  763  298  741  213  130  749  249  747  232  506 1193  147
  F4   F5   G1   G2   G3   G4   G5   H1   H2   H3   H4   H5    m    f
 236   25   94  425  180  888  520 1105  917   52   31    2  972 1135
  ma   wi   di   se   si   e0   e1   e2   e3   e4   e5   a1   a2   a3
1169  172   45   69  652  240  522  580  431  166  168  334  453  374 
  a4   a5   a6 
 325  251  370 

# amalgamate the H4 and H5 columns (39 & 40) into column 39 and then 
# delete the 40th column
women.Z[,39] <- women.Z[,39] + women.Z[,40]
women.Z      <- women.Z[,-40]
colnames(women.Z)[39] <- "H4,5"

# Compute the 39x39 Burt matrix pf the substantive variables from the 
# indicator matrix...
women.Burt <- t(women.Z[,1:39]) %*% women.Z[,1:39]

# ... or the Burt matrix is a by-product of the mjca function in the 
# ca package
women.Burt <- mjca(women[,1:8])$Burt

# after which you have to combine the categories H4 and H5
women.Burt[,39] <- women.Burt[,39]+women.Burt[,40]
women.Burt[39,] <- women.Burt[39,]+women.Burt[40,]
women.Burt <- women.Burt[-40,-40]
rownames(women.Burt)[39] <- colnames(women.Burt)[39] <- "H4,5"

# for the concatenated table we need to first create the interactive 
# variable and add it to the original data then compute the full Burt 
# matrix, then extract the required part of the Burt matrix and 
# combine H4 and H5 again
women               <- cbind(women, 6*(women[,"g"]-1)+ women[,"a"])
colnames(women)[13] <- "ga"
women.concat        <- mjca(women)$Burt[40+c(3:13,20:31),1:40]
women.concat[,39]   <- women.concat[,39]+women.concat[,40]
women.concat        <- women.concat[,1:39]
colnames(women.concat)[39] <- "H4,5"

# The matrix is correct, now just the labels of the interactively 
# coded variable have to be redefined!
# Presently they are ga1, ga2, ..., ga12; we want to code the first 6 
# as male, and the last 6 as female
rownames(women.concat)[12:23] <- paste(c(rep("ma", 6), 
                                       rep("fa", 6)), 1:6, sep = "")

# After that detour, back to MCA!
# The total inertia of the Burt matrix (check that it's 39x39 first)
dim(women.Burt)
[1] 39 39

sum(ca(women.Burt)$sv^2)
[1] 0.677625

# Adjusted total ienrtia (average off-diagonal inertia)
(8/7)*(sum(ca(women.Burt)$sv^2)-(31/64))
[1] 0.2208571

# Before continuing, check that you have the correct indicator matrix, 
# it should be 2107 x 39 in this case.
# Above we have described either reading in the full indicator matrix 
# from an Excel file, or computing it.
# In either case you probably have to extract the first 40 columns 
# corresponding to the substantive variables, then do the combining of 
# H4 and H5 in the indicator matrix:
women.Z      <- women.Z[,1:40]
women.Z[,39] <- women.Z[,39] + women.Z[,40]
women.Z      <- women.Z[,-40]

dim(women.Z)
[1] 2107 39

# total inertia of the 2107x39 indicator matrix:
sum(ca(women.Z)$sv^2)
[1] 3.875

# Exhibit 10.3 (with vertical axis possibly inverted, depending on R 
# version). More code given here for category labelling
par(cex.axis = 0.7)
plot(ca(women.Z), map = "rowprincipal", labels = 0, 
     col = c("lightgreen","chocolate"), pch = c(149,1,17,24), 
     mass = c(FALSE,TRUE))
women.Z.ca <- ca(women.Z)
text(women.Z.ca$colcoord, labels = women.Z.ca$colnames, pos = 2, 
     offset = 0.3, col = "chocolate", cex = 0.8, font = 4)

# How many singular values are greater than 1/8:
which(ca(women.Burt)$sv > (1/8))
[1] 1 2 3 4 5 6 7 8 9

# The adjusted singular values for these first 9 axes
(8/7)*(ca(women.Burt)$sv[1:9]-(1/8))
[1] 0.342185682 0.232597737 0.124150258 0.114994566 0.034506798 
[6] 0.025747549 0.014892541 0.008973295 0.006810267

# The parts of (adjusted) inertia on the first 9 axes
(64 / 49) * (ca(women.Burt)$sv[1:9] - (1/8))^2 / 0.2208571
[1] 0.5301665244 0.2449624996 0.0697885037 0.0598746893 0.0053913554 
[6] 0.0030016525 0.0010042139 0.0003645797 0.0002099988

# Compute MCA using ca function and then replace the singular values 
# with the adjusted ones
women.Burt.ca <- ca(women.Burt)
women.Burt.ca$sv[1:9] <- (8/7) * (ca(women.Burt)$sv[1:9] - (1/8))

# EXTRA: Exhibit 10.4 - first axis coordinates inverted to agree with 
# Exhibit, although this depends on version of R
# In Exhibit 10.4, the row labels were made lower case externally
women.Burt.ca$rowcoord[,1] <- -women.Burt.ca$rowcoord[,1]
women.Burt.ca$colcoord[,1] <- -women.Burt.ca$colcoord[,1]

plot(women.Burt.ca, map = "rowprincipal", labels = 0, 
     pch = c(149,1,17,24), col = c("lightgreen","chocolate"), 
     mass = c(FALSE,TRUE))
text(women.Burt.ca$colcoord, labels = women.Burt.ca$colnames, 
     pos = 2, offset = 0.3, col = "chocolate", cex = 0.8, font = 4)
women.Burt.Fadj <- women.Burt.ca$colcoord[,1:9] %*% 
                     diag((8/7) * (ca(women.Burt)$sv[1:9] - (1/8)))
points(women.Burt.Fadj, pch = 19, col = "lightgreen", cex = 0.7)
text(women.Burt.Fadj, labels = rownames(women.Burt), pos = 4, 
     offset = 0.3, col = "forestgreen", cex = 0.7, font = 2)

# EXTRA: Exhibit 10.5 - we assume that the singular values of 
# women.Burt.ca have already been substituted, and that the row 
# principal coordinates women.Burt.Fadj computed, as above.
# The column coordinates are shrunk by the square roots of their 
# masses, and then plotted as before:
women.Burt.Gctr <- sqrt(women.Burt.ca$colmass) * 
                     women.Burt.ca$colcoord
women.Burt.ca$colcoord <- women.Burt.Gctr

plot(women.Burt.ca, map = "rowprincipal", labels = 0, 
     pch = c(149,1,17,24), col = c("lightgreen","chocolate"))
text(women.Burt.Gctr, labels = women.Burt.ca$colnames, pos = 2, 
     offset = 0.3, col = "chocolate", cex = 0.8, font = 4)
points(women.Burt.Fadj, pch = 19, col = "lightgreen", cex = 0.7)
text(women.Burt.Fadj, labels = women.Burt.ca$rownames, pos = 2, 
     offset = 0.3, col = "forestgreen", cex = 0.7, font = 2)

# Exhibit 10.6 - some of the demographic category labels have been 
# edited externally in the book's version
women.BurtS.ca <- ca(rbind(women.Burt, women.concat), suprow = 40:62)
women.BurtS.ca$rowcoord[,1] <- -women.BurtS.ca$rowcoord[,1]
women.BurtS.ca$colcoord[,1] <- -women.BurtS.ca$colcoord[,1]
women.BurtS.Gctr <- sqrt(women.BurtS.ca$colmass) * 
                      women.BurtS.ca$colcoord
women.BurtS.ca$colcoord <- women.BurtS.Gctr
women.BurtS.ca$sv[1:9]  <- (8/7) * (women.BurtS.ca$sv[1:9]-(1/8))
plot(women.BurtS.ca, map = "rowprincipal", what = c("none","all"), 
     col = "chocolate", labels = 0, pch = c(20,1,17,24))
text(women.BurtS.Gctr, labels = women.Burt.ca$colnames, pos = 2, 
     offset = 0.3, col = "chocolate", cex = 0.8, font = 4)
women.BurtS.Fsup <- women.BurtS.ca$rowcoord[40:62,] %*% 
                      diag(women.BurtS.ca$sv)
points(women.BurtS.Fsup, pch = 21, col = "forestgreen")
text(women.BurtS.Fsup, labels = women.BurtS.ca$rownames[40:62], 
     pos = 2, offset = 0.3, col = "forestgreen", cex = 0.7, font = 2)

# EXTRA: an easier way to do the MCA using the mjca function in the 
# ca package.
# This package takes in original data, so one has to convert the 
# original data first (e.g., combining H4 and H5)
# before reading into the function
# So let's start again and read in the original data, cut out the 
# missing data, combine H4 and H5, and creating the interactive 
# variable, and then use the package.
# Suppose 'women' contains the original data again (2471 x 12) 

dim(women)
[1] 2471   12

# Take out any rows with missing data 
women.9 <- rep(0, nrow(women))
for(i in 1:nrow(women)) {
  for(j in 1:12) {
    if(women[i,j] == 9) women.9[i] <- women.9[i] + 1
    }
  }
women <- women[women.9 == 0,]

dim(women)
[1] 2107   12

# Replace 5's for variable 'H' with 4's
women[women[,"H"] == 5, "H"] <- 4

# Create interactive variable
women <- cbind(women, 6 * (women[,"g"]-1) + women[,"a"])
colnames(women)[13] <- "ga"

# Now use mjca function directly
plot(mjca(women[,1:8]), what = c("none","all"))

plot(mjca(women[,1:8]), what = c("none","all"), map = "rowgreen")

# Notice that to get exactly Exhibit 10.6, we will have to work via 
# the Burt matrix
women.BURT <- mjca(women)$Burt

dim(women.BURT)
[1] 70 70

#  8 questions + gender + marital + education + age + age*gender
#      39      +   2    +    5    +     6     +  6  +    12      = 70
plot(ca(women.BURT[c(1:39, 42:52, 59:70), 1:39], suprow = 40:62), 
     what = c("passive","all"), map = "rowgreen")

# only problem here is the scaling factor for the rows, because the 
# inertias are not adjusted