Throughout this section we choose $n=10,000$. Our first normal croissant (category quantifications on the left, object scores on the right) has four variables, knots at the integers from -3 to +3, and all correlations equal to .75. We see somewhat more scatter or filling in the object score croissant, basically because object scores are convex combinations of category quantifications. ```{r discnorm1, fig.align = "center", echo = FALSE} par(mfrow = c(1, 2)) # 4 variables, r = .75 set.seed (12345) x1 <- discreteNormal(10000, 4, matrix (.75, 4, 4) + .25 * diag (4), repList (c(-3,-2,-1,0,1,2,3), 4)) b1 <- burtTable (x1) c1 <- b1$c g1 <- geigen (c1, diag (diag (c1)), symmetric = TRUE) y1 <- g1$vectors [, c(nrow(c1) - 1, nrow (c1) - 2)] plot (y1, col = "RED", xlab = "dimension 1", ylab = "dimension 2", main = "r = .75, m = 4") plot (b1$g %*% y1, col = "RED", xlab = "dimension 1", ylab = "dimension 2", main = "r = .75, m = 4") ```

If we decrease the correlation to .25 the croissant crumbles. A correlation of .25 is too close to independence, which means that not enough variation will be captured by the first two eigenvalues. Nevertheless if we discard the outliers the croissant is still there.

```{r discnorm2, fig.align = "center", echo = FALSE, cache = TRUE} par(mfrow = c(1, 2)) # 4 variables, r = .25 set.seed (12345) x1 <- discreteNormal(10000, 4, matrix (.25, 4, 4) + .75 * diag (4), repList (c(-3,-2,-1,0,1,2,3), 4)) b1 <- burtTable (x1) c1 <- b1$c g1 <- geigen (c1, diag (diag (c1)), symmetric = TRUE) y1 <- g1$vectors [, c(nrow(c1) - 1, nrow (c1) - 2)] plot (y1, col = "RED", xlab = "dimension 1", ylab = "dimension 2", main = "r = .25, m = 4") plot (b1$g %*% y1, col = "RED", xlab = "dimension 1", ylab = "dimension 2", main = "r = .25, m = 4") ```

In the third simulation we increase the number of variables to 25, with correlation .75. The croissants are tight and symmetric.

```{r discnorm3, fig.align = "center", echo = FALSE, cache = TRUE} par(mfrow = c(1, 2)) # 25 variables, r = .75 set.seed (12345) x1 <- discreteNormal(10000, 25, matrix (.75, 25, 25) + .25 * diag (25), repList (c(-3,-2,-1,0,1,2,3), 25)) b1 <- burtTable (x1) c1 <- b1$c g1 <- geigen (c1, diag (diag (c1)), symmetric = TRUE) y1 <- g1$vectors [, c(nrow(c1) - 1, nrow (c1) - 2)] plot (y1, col = "RED", xlab = "dimension 1", ylab = "dimension 2", main = "r = .75, m = 25") plot (b1$g %*% y1, col = "RED", xlab = "dimension 1", ylab = "dimension 2", main = "r = .75, m = 25") ```

Moreover for 25 variables we still see a clear croissant if the correlation is .25, although there obviously is more filling. Lower correlation means objects scores will be convex combinations of more category quantifications, and thus they will cluster more around the origin, which is the fattest part of the croissant.

```{r discnorm4, fig.align = "center", echo = FALSE} par(mfrow = c(1, 2)) # 25 variables, r = .25 set.seed (12345) x1 <- discreteNormal(10000, 25, matrix (.25, 25, 25) + .75 * diag (25), repList (c(-3,-2,-1,0,1,2,3), 25)) b1 <- burtTable (x1) c1 <- b1$c g1 <- geigen (c1, diag (diag (c1)), symmetric = TRUE) y1 <- g1$vectors [, c(nrow(c1) - 1, nrow (c1) - 2)] plot (y1, col = "RED", xlab = "dimension 1", ylab = "dimension 2", main = "r = .25, m = 25") plot (b1$g %*% y1, col = "RED", xlab = "dimension 1", ylab = "dimension 2", main = "r = .25, m = 25") ```

We can also study the influence of skewness by choosing the knots to be `c(0,1,2,3)`, keeping the correlaton between the four variables at .75. A truncated versions of the croissant appears.

```{r discnorm5, fig.align = "center", echo = FALSE, cache = TRUE} par(mfrow = c(1, 2)) # 4 variables, r = .75, skew set.seed (12345) x1 <- discreteNormal(10000, 4, matrix (.75, 4, 4) + .25 * diag (4), repList (c(0,1,2,3), 4)) b1 <- burtTable (x1) c1 <- b1$c g1 <- geigen (c1, diag (diag (c1)), symmetric = TRUE) y1 <- g1$vectors [, c(nrow(c1) - 1, nrow (c1) - 2)] plot (y1, col = "RED", xlab = "dimension 1", ylab = "dimension 2", main = "r = .75, m = 4, skew") plot (b1$g %*% y1, col = "RED", xlab = "dimension 1", ylab = "dimension 2", main = "r = .75, m = 4, skew") ```

Finally we can choose a correlation matrix which does not have all correlations the same. We use the correlation matrix ```{r pcor, echo = FALSE} print (r <- matrix (c(1.000,-0.027,0.729,0.008,-0.027,1.000,0.001,0.343,0.729,0.001,1.000,-0.027,0.008,0.343,-0.027,1.000), 4, 4)) ``` This has the effect of perturbing the croissant, because it becomes a mixture of different croissants for the different variables.

```{r discnorm6, fig.align = "center", echo = FALSE, cache = TRUE} par(mfrow = c(1, 2)) # 4 variables, r matrix r <- matrix (c(1.000,-0.027,0.729,0.008,-0.027,1.000,0.001,0.343,0.729,0.001,1.000,-0.027,0.008,0.343,-0.027,1.000), 4, 4) set.seed (12345) x1 <- discreteNormal(10000, 4, r, repList (c(-3,-2,-1,0,1,2,3), 4)) b1 <- burtTable (x1) c1 <- b1$c g1 <- geigen (c1, diag (diag (c1)), symmetric = TRUE) y1 <- g1$vectors [, c(nrow(c1) - 1, nrow (c1) - 2)] plot (y1, col = "RED", xlab = "dimension 1", ylab = "dimension 2", main = "r varies, m = 4") plot (b1$g %*% y1, col = "RED", xlab = "dimension 1", ylab = "dimension 2", main = "r varies, m = 4") par(mfrow = c(1, 1)) ```

In real data we have a combination of all these effects. We will have low correlations and/or skewness for some of the variables. Thus croissants can be asymmetric and have a lot of scatter. And, of course, there may be deviations from underlying normality as well, which can easily destroy the croissant. #Real Data: GALO The GALO example has served us well for almost 40 years. The objects (individuals) are 1290 school children in the sixth grade of elementary school in the city of Groningen (Netherlands) in 1959. There are four variables * Gender: M/F. * IQ: The original range (60 to 144) has been categorized into 9 ordered categories. * SES: + LoWC = Lower white collar; + MidWC = Middle white collar; + Prof = Professional, Managers; + Shop = Shopkeepers; + Skil = Schooled labor; + Unsk = Unskilled labor. * Teacher's Advice Secondary Education: + Agr = Agricultural; + Ext = Extended primary education; + Gen = General; + Grls = Secondary school for girls; + Man = Manual, including housekeeping; + None = No further education; + Uni = Pre-University. ##MCA ```{r galodata, echo = FALSE, cache = TRUE} data(galo, package="homals") x <- galo[,1:4] tgalo <- burtTable (x, rep(-1, 4)) cgalo <- tgalo$c ogalo <- tgalo$ord ggalo <- geigen (cgalo / 4, diag (diag (cgalo)), symmetric = TRUE) y <- ggalo$vectors[,c (nrow (cgalo) - 1, nrow (cgalo) - 2)] s <- ggalo$vectors[1, nrow (cgalo)] y1 <- y [1:ogalo[1], ] / s y2 <- y [ogalo[1] + (1:ogalo[2]), ] / s y3 <- y [ogalo[1] + ogalo[2] + (1:ogalo[3]), ] / s y4 <- y [ogalo[1] + ogalo[2] + ogalo[3] + (1:ogalo[4]), ] / s xo <- (tgalo$g %*% y) / 4 ``` We use the first two dimensions, which have eigenvalues `r ggalo$values[nrow (cgalo) - 1]` and `r ggalo$values[nrow (cgalo) - 2]`. The category quantifications from MCA for each variable are in `r figure_nums("galocat")`. Small croissants are visible throughout.

```{r galocat, echo = FALSE, fig.align= "center"} par(mfrow = c(2, 2)) plot (y1, type = "n", main = "GENDER", xlab = "dimension 1", ylab = "dimension 2", xlim = c(-4,4), ylim = c(-4,4)) text (y1, as.character(unique(galo[,1])), col = "RED") plot (y2, type = "n", main = "IQ", xlab = "dimension 1", ylab = "dimension 2", xlim = c(-4,4), ylim = c(-4,4)) text (y2, as.character(unique(galo[,2])), col = "RED") plot (y3, type = "n", main = "ADVICE", xlab = "dimension 1", ylab = "dimension 2", xlim = c(-4,4), ylim = c(-4,4)) text (y3, as.character(unique(galo[,3])), col = "RED") plot (y4, type = "n", main = "SES", xlab = "dimension 1", ylab = "dimension 2", xlim = c(-4,4), ylim = c(-4,4)) text (y4, as.character(unique(galo[,4])), col = "RED") par(mfrow = c(1, 1)) ```

The croissant is somewhat clearer if we put the category quantification of all variables in a single plot.

```{r galotcat, echo = FALSE, fig.align = "center"} plot(y, col="RED", main = "GALO Category Quantifications", xlab = "dimension 1", ylab = "dimension 2") ```

We plot the object scores in a somewhat nonstandard way, by putting unnormalized object scores (centroids of category quantifications) within the hull of the category quantifications. This illustrates that, at least to some extent, the shape of the object score plot is determined by the shape of the category quantification plot. We see an obvious croissant for these data.

```{r galojoint, fig.align = "center", echo = FALSE} plot(y, xlab = "dimension 1", ylab = "dimension 2", type = "n") i <- chull (y) points (y[i,], col = "RED") lines(y[c(i, i[1]), ], col = "BLUE", lwd = 3) points (xo, col = "RED") ```

##Two-step MCA: PRIMALS One of the standard ways to get rid of MCA croissants is to impose constraints on the category quantifications. In his canonical MCA paper @guttman_41 already argued that only the first dimension of an MCA should be used in linear multivariate analysis, because the remaining dimensions are components of non-linearity. The MCA does not decompose variance, as PCA does in classical multivariate analysis, but it decomposes chi-square, and thus its interpretation should be completely different. There are various ways to incorporate this insight. The simplest one is to use the first MCA dimension to quantify the variables, and to use these quantified variables in a standard linear technique such as PCA or regression. The combination of MCA quantification and PCA is sometimes called PRIMALS [@meulman_gifi_81]. ```{r primals, echo = FALSE} yd <- directSum (list (y1[, 1, drop = FALSE], y2[, 1, drop = FALSE], y3[, 1, drop = FALSE], y4[, 1, drop = FALSE])) rr <- crossprod (yd, cgalo %*% yd) rownames (rr) <- colnames (rr) <- c ("SEX", "IQ", "ADV", "SES") dr <- sqrt (diag (rr)) r <- rr / outer (dr, dr) e <- eigen (r) ev <- e$values ee <- e$vectors ef <- cbind (-ee[,1] * sqrt (ev[1]), ee[,2] * sqrt (ev[2])) ``` The correlation matrix of the quantified variables is ```{r cormat, echo = FALSE} r ``` with eigenvalues (divided by 4) ```{r coreval, echo = FALSE} ev / 4 ``` and eigenvectors ```{r corevec, echo = FALSE} ee ``` Using the first two eigenvectors we can make a PCA plot.

```{r galopca, fig.align = "center", echo = FALSE} plot (ef, xlim=c(-1,1), ylim=c(-1,1), type = "n", main = "PCA Plot", asp =1, xlab = "dimension 1", ylab = "dimension 2") text (ef, c("SEX", "IQ", "ADV", "SES")) arrows (0, 0, ef[1, 1], ef[1, 2], col = "RED", lwd = 2) arrows (0, 0, ef[2, 1], ef[2, 2], col = "RED", lwd = 2) arrows (0, 0, ef[3, 1], ef[3, 2], col = "RED", lwd = 2) arrows (0, 0, ef[4, 1], ef[4, 2], col = "RED", lwd = 2) w<-seq(0, 2*pi, length=100) lines (sin (w), cos (w), col = "BLUE", lwd = 2) abline (h = 0) abline (v = 0) ```

To show how this analysis is related to MCA, we first normalize the category quantifications of variable $j$ to $y_j'D_jy_j=1$. We then create the first column of $Y$ by multiplying the $y_j$ by the corresponding entries of the first eigenvector of the correlation matrix. And do the same for the second column of $Y$. We still have $Y'DY=I$, while $\mathbf{tr}\ Y'CY$ is equal to the sum of the first two eigenvalues of the correlation matrix, i.e. `r (ev[1]+ev[2])/4`. For the MCA the corresponding sum is `r sum (ggalo$values[c(nrow (cgalo) - 1, nrow (cgalo) - 2)])`. ```{r primsol, echo = FALSE} dd <- 1 / sqrt (diag (crossprod(yd, diag(cgalo) * yd))) yd <- t (dd * t(yd)) yy <- yd %*% ee[,1:2] ```

```{r primcat, echo = FALSE, fig.align = "center"} plot(yy, col = "RED", xlab = "dimension 1", ylab = "dimension 2") lines (yy[1:ogalo[1], ], col = "BLUE") lines (yy[ogalo[1] + 1:ogalo[2], ], col = "BLUE") lines (yy[ogalo[1] + ogalo[2] + 1:ogalo[3], ], col = "BLUE") lines (yy[ogalo[1] + ogalo[2] + ogalo[3] + 1:ogalo[4], ], col = "BLUE") ```

The corresponding object scores are the left singular vectors of the matrix of quantified variables. There is no croissant in sight any more.

```{r primobj, echo = FALSE, fig.align = "center"} plot ((tgalo$g %*% yy) / 4, col = "RED", xlab = "dimension 1", ylab = "dimension 2") ```

##Nonlinear PCA: PRINCALS The PRIMALS approach is a two-step approach, in the sense that we first maximize by choosing the largest eigenvalue of the Burt table and then we maximize by choosing the largest eigenvalues of the induced correlation matrix. Thus there are two different criteria involved, and this could be seen as a disadvantage. In the PRINCALS approach we choose quantifications of the variables that maximize the sum of the first p eigenvalues of the induced correlation matrix. For $p=1$ this is PRIMALS, but for $p>1$ the two approaches are different. We perform the necessary computations, in the GALO case for $p=2$, using the `homals` package [@deleeuw_mair_A_09a]. Alternatively, the `aspect` package [@mair_deleeuw_A_10] could be used. ```{r galohomals, echo = FALSE, cache = TRUE} library (homals) h <- homals (x, rank = 1, eps = 1e-10) hc <- h$catscores y <- rbind(hc[[1]],hc[[2]],hc[[3]],hc[[4]]) ``` ```{r princals, echo = FALSE} yd <- directSum (list (hc[[1]][, 1, drop = FALSE], hc[[2]][, 1, drop = FALSE], hc[[3]][, 1, drop = FALSE], hc[[4]][, 1, drop = FALSE])) rr <- crossprod (yd, cgalo %*% yd) rownames (rr) <- colnames (rr) <- c ("SEX", "IQ", "ADV", "SES") dr <- sqrt (diag (rr)) r <- rr / outer (dr, dr) e <- eigen (r) ev2 <- e$values ee <- e$vectors ef <- cbind (-ee[,1] * sqrt (ev2[1]), ee[,2] * sqrt (ev2[2])) ``` The correlation matrix of the quantified variables is ```{r cormatr, echo = FALSE} r ``` with eigenvalues (divided by 4) ```{r corevalr, echo = FALSE} ev2 / 4 ``` and eigenvectors ```{r corevecr, echo = FALSE} ee ``` The largest eigenvalue is now `r ev2[1]/4`, while for PRIMALS it was the larger value `r ev[1]/4`. But the sum of the two largest eigenvalues for PRINCALS is `r (ev2[1]+ev2[2])/4`, while for PRIMALS it is the smaller value `r (ev[1]+ev[2])/4`. Using the first two eigenvectors of the induced correlation matrix we can make a PCA plot.

```{r galoprin, fig.align = "center", echo = FALSE} plot (ef, xlim=c(-1,1), ylim=c(-1,1), type = "n", main = "PCA Plot", asp =1, xlab = "dimension 1", ylab = "dimension 2") text (ef, c("SEX", "IQ", "ADV", "SES")) arrows (0, 0, ef[1, 1], ef[1, 2], col = "RED", lwd = 2) arrows (0, 0, ef[2, 1], ef[2, 2], col = "RED", lwd = 2) arrows (0, 0, ef[3, 1], ef[3, 2], col = "RED", lwd = 2) arrows (0, 0, ef[4, 1], ef[4, 2], col = "RED", lwd = 2) w<-seq(0, 2*pi, length=100) lines (sin (w), cos (w), col = "BLUE", lwd = 2) abline (h = 0) abline (v = 0) ```

In the same way as for PRIMALS we can make the PRINCALS category quantifications and object scores.

```{r princat, fig.align = "center", echo = FALSE} plot(y, col = "RED", xlab = "dimension 1", ylab = "dimension 2") lines(hc[[1]], col = "BLUE") lines(hc[[2]], col = "BLUE") lines(hc[[3]], col = "BLUE") lines(hc[[4]], col = "BLUE") ```

```{r prinobj, fig.align = "center", echo = FALSE} plot((tgalo$g %*% y) / 4, col = "RED", xlab = "dimension 1", ylab = "dimension 2") ```

##Simultaneous Bilinearizing MCA A more complete and more revealing analysis is possible by using the theory in @deleeuw_C_82, see also @bekker_deleeuw_C_88 and @deleeuw_E_15g. We study an alternative way of diagonalizing the matrix $E=D^{-\frac12}CD^{-\frac12}$. Suppose we can find square orthonormal $K_j$ such that $K_j'K_j=I$ and $K_jE_{j\ell}K_\ell$ is diagonal for all $j,\ell=1,\cdots,m$, say equal to the diagonal matrices $\Gamma_{j\ell}$. The condition is that if two submatrices of $E$ have a subscript in common, then then also have the corresponding set of singular vectors in common. @deleeuw_A_88a calls this *simultaneous bilinearizability*. Since that name is a bit verbose, we call the condition *SBL*. We know SBL is possible for the continuous multinormal, with the $K_j$ the Hermite-Chebyshev polynomials. It is also possible if all variables are binary, and it is possible if $m=2$. Those are some important gauges, and we expect SBL to be true approximately in many practical situations. Define the direct sum $K:=K_1\oplus\cdots\oplus K_m$, and $\Gamma:=K'EK$. Under SBL all $\Gamma_{j\ell}$ are diagonal. We can then find a permutation $P$ of rows and columns such that $R:=P'\Gamma P$ is the direct sum of a number of correlation matrices. If we start with $m$ variables with $k$ categories each, then $R=R_1\oplus\cdots\oplus R_k$, where each $R_s$ is of order $m$. The eigenvectors of $R$ are the direct sum $L:=L_1\oplus\cdots\oplus L_k$, and under SBL we have $$ L'RL=L'P'\Gamma PL=L'P'K'EKPL=m\Lambda. $$ Thus under SBL the eigenvectors of $E$ are given by $Z=KPL$ and the eigenvalues of $E$ are also the eigenvalues of $\Gamma$ and also the eigenvalues of $R$. Thus for each eigenvalue of $E$ we can say which is the $R_s$ it belongs to. If the variables have a different number of categories then under SBL we have $k$ matrices $R_s$, where $k$ is now the maximum of the $k_j$. The $R_s$ are generally of different orders, because variable $j$ is represented in $R_s$ if and only if $k_j\geq s$. Thus a binary variable only occurs in $R_1$ and $R_2$. In all cases the matrix $R_1$ can be chosen to be the $m\times m$ matrix with all elements equal to one. For all variables binary the matrix $R_2$ is the matrix of phi-coefficients (point correlations). Now SBL will not be precisely be satisfied for empirical data, except in the trivial cases $k_j\equiv 2$ or $m=2$. To create a suitable dat analysis techniue we choose the $K_j$ in such a way that the sum of squares of the off-diagonal elements of the $\Gamma_{j\ell}$ is minimized, or, equivalently, the sum of squares of the diagonal elements is maximized. The details and the code are in @deleeuw_R_08a. For this paper the approximate diagonalization procedure is in the function `threeStepApprox()` which uses `kplSVD()` to find the $K_j$ and `kplPerm()` to find the permutation $P$. Because the SBL theory is somewhat more complicated than other MCA theory, we analyze a simple normal example in detail. ```{r normex} set.seed (12345) x <- discreteNormal(1000, 4, matrix (.75, 4, 4) + .25 * diag (4), repList (c(-2,2), 4)) b <- burtTable (x) d <- 1 / sqrt (diag (b$c)) e <- b$c * outer (d, d) h <- threeStepApprox (e, rep(3, 4)) ``` The Burt Table is ```{r normex_b, echo = FALSE} print (c <- b$c, digits = 2) ``` The normalized Burt Table $E$ is ```{r normex_e, echo = FALSE} ppp <- apply (e, 1, function (x) cat (formatC (x, digits = 4, width = 7, format = "f"), "\n")) ``` and its eigenvalues are ```{r normex_v, echo = FALSE} cat (formatC (eigen(e)$values, digits = 4, width = 7, format = "f")) ``` After application of `kplSVD()` to $E$ we have ```{r normex_a, echo = FALSE} ppp <- apply (h$bval, 1, function (x) cat (formatC (x, digits = 4, width = 7, format = "f"), "\n")) ``` Permuting rows and columns gives the matrix $R=R_1\oplus R_2\oplus R_3$ ```{r normex_c, echo = FALSE} ppp <- apply (h$cval, 1, function (x) cat (formatC (x, digits = 4, width = 7, format = "f"), "\n")) ``` And we can now diagonalize the blocks. ```{r normex_d, echo = FALSE} ppp <- apply (h$dval, 1, function (x) cat (formatC (x, digits = 4, width = 7, format = "f"), "\n")) ``` and we can sort the diagonal elements ```{r normex_w, echo = FALSE} cat (formatC (sort (diag (h$dval), decreasing = TRUE), digits = 4, width = 7, format = "f")) ``` to find they are really close to the eigenvalues of $E$. The correlations between the eigenvector of $E$ and the SBL approximations are ```{r normex_x, echo = FALSE} ppp <- apply (h$cc, 1, function (x) cat (formatC (x, digits = 4, width = 7, format = "f"), "\n")) ``` Thus the eight non-trivial eigenvectors of $E$ correspond, repectively, with first of $R_2$, first of $R_3$, second of $R_3$, second of $R_2$, third of $R_2$, third of $R_3$, fourth of $R_2$, fourth of $R_3$. In terms of the polynomials underlying the category quantifications this means linear, quadratic, quadratic, linear, linear, quadratic, linear, quadratic. ```{r 3step, echo = FALSE, cache = TRUE} d <- 1 / sqrt (diag (cgalo)) e <- cgalo * outer (d, d) h <- threeStepApprox (e, ogalo) ``` Back to GALO. For GALO the eigenvalues of $E/4$ are ```{r egalo1, echo = FALSE} cat (formatC (h$eval/4, digits = 4, width = 10, format = "f")) ``` The direct sum of the eigenvalues of the $R_s/m$, which have orders `r h$gg`, is ```{r egalo2, echo = FALSE} hg <- h$gg isum <- 0 for (i in 1:length(hg)) { cat (formatC (h$aval[isum + 1:hg[i]]/4, digits = 4, width = 10, format = "f"), "\n") isum <- isum + hg[i] } ``` and if we sort these we find ```{r egalo3, echo = FALSE} cat (formatC (sort(h$aval, decreasing = TRUE)/4, digits = 4, width = 10, format = "f")) ``` Note that the orders of the matrices $R_s$ in this case are `r h$gg`. Clearly the approximate eigenvalues from the KPL diagonalization and the actual eigenvalues are very similar, especially the larger ones and the smaller ones. The correlation between actual eigenvectors and approximate eigenvectors is $Q:=Z'KPL$. Let us look at the first three non-trivial eigenvectors of $E$, which correspond with rows 2, 3, and 4 of $Q$. Regress them on the non-trivial columns of $KPL$, which are all columns except the first four. The percentage of variance explained by the colums of $KPL$ are the squares of the elements of rows 2, 3, and 4 of $Q$. ```{r correle, echo = FALSE} cat (formatC (h$cc[2,-(1:4)]^2, digits = 2, format ="f")) cat (formatC (h$cc[3,-(1:4)]^2, digits = 2, format ="f")) cat (formatC (h$cc[4,-(1:4)]^2, digits = 2, format ="f")) ``` Thus we see the first non-trivial eigenvector of $E$ corresponds with the first eigenvector of $R_2$, which the next eigenvector of $E$ corresponds with the first eigenvector of $R_3$. The third non-trivial eigenvector of $E$ is most closely related to the second eigenvector of $R_2$, but there are also contributions of the first eigenvectors of $R_4$ and $R_5$. There is a somewhat more rigorous method to match the eigenvectors in $Z$ and $KPL$. The Hadamard product $Q*Q$ is doubly stochastic. Finding the permutation matrix closest in the least suares sense to $Q*Q$ is a linear assignment problem that can be solved by using `lp.assign()` from the `lpSolve` package [@lpsolve_15]. We indicate for each row the column where the permutation matrix has a one. ```{r assign, echo = FALSE} library(lpSolve) q <- lp.assign(h$cc^2,direction="max")$solution apply (q, 1, which.max) ``` Again this indicates that the first three non-trivial eigenvectors of $E$ correspond with the dominant eigenvalue of $R_2$, the dominant eigenvalue of $R_3$, and the second eigenvalue of $R_2$. If the eigenvectors in $K$ correspond with the orthogonal polynomials, as they do in the multinormal case, then the first two nontrivial eigenvectors produce a croissant, because they come from the linear and quadratic transformations, respectively. But if we select the first and the third non-trivial eigenvector of $E$ to plot, then they both come from the linear quantifications and they do not produce a croissant. We first plot $Y=D^{-\frac12}KPL$, using the columns corresponding to the first two non-trivial eigenvectors of $E$. This shows the croissant.

```{r 3galocatm, fig.align= "center", echo = FALSE} yy <- d * h$avec [, c(5, 9)] plot (yy[, ], col ="RED", xlab = "dimension 1", ylab = "dimension 2") yv <- yy[1:ogalo[1], ] ov <- order (yv[, 1]) lines (yv[ov, ], col = "BLUE") yv <- yy[ogalo[1] + 1:ogalo[2], ] ov <- order (yv[, 1]) lines (yv[ov, ], col = "BLUE") yv <- yy[ogalo[1] + ogalo[2] + 1:ogalo[3], ] ov <- order (yv[, 1]) lines (yv[ov, ], col = "BLUE") yv <- yy[ogalo[1] + ogalo[2] + ogalo[3] + 1:ogalo[4], ] ov <- order (yv[, 1]) lines (yv[ov, ], col = "BLUE") ```

If we choose the first and second eigenvector of $R_2$, i.e. the solution corresponding to the first and third non-trivial eigenvectors of $E$, there is no croissant. In fact, the solution is similar to the one given by PRIMALS, consisting of a line through the origin for each variable, with the category quantifications on the line.

```{r 3galocata, fig.align= "center", echo = FALSE} yy <- d * h$avec [, c(5, 6)] plot (yy[, ], col ="RED", xlab = "dimension 1", ylab = "dimension 2") lines (yy[1:ogalo[1], ], col = "BLUE") lines (yy[ogalo[1] + 1:ogalo[2], ], col = "BLUE") lines (yy[ogalo[1] + ogalo[2] + 1:ogalo[3], ], col = "BLUE") lines (yy[ogalo[1] + ogalo[2] + ogalo[3] + 1:ogalo[4], ], col = "BLUE") ```

We can compare this to the solution we get by actaully plotting the first and third non-trivial eigenvector of $E$. There is no croissant in this case either.

```{r 3galocate, fig.align= "center", echo = FALSE} yy <- d * h$evec [, c(2, 4)] plot (yy[, ], col ="RED", xlab = "dimension 1", ylab = "dimension 3") yv <- yy[1:ogalo[1], ] ov <- order (yv[, 1]) lines (yv[ov, ], col = "BLUE") yv <- yy[ogalo[1] + 1:ogalo[2], ] ov <- order (yv[, 1]) lines (yv[ov, ], col = "BLUE") yv <- yy[ogalo[1] + ogalo[2] + 1:ogalo[3], ] ov <- order (yv[, 1]) lines (yv[ov, ], col = "BLUE") yv <- yy[ogalo[1] + ogalo[2] + ogalo[3] + 1:ogalo[4], ] ov <- order (yv[, 1]) lines (yv[ov, ], col = "BLUE") ```

#Real Data: Norway The data is register data from Norway, with information on a high number (n = 159539) of individuals. The analysis is part of a social mobility study conducted by Arild Danielsen, Department of History, Sociology and Innovation, Faculty of Economics and Social Sciences, University College of Southeast Norway. See @danielsen_15. The data have four variables with data on class-fractions (based on occupational categories) of the individualâ€™s two parents (in 2008) and two grand-fathers (in 1980). ```{r read_data, echo = FALSE} x <- read.table ("Four\040variables\040for\040MCA\040-\04001.12.15.dat", header = TRUE, sep = "\t") row.names (x) <- as.character (x[, 1]) x <- x[, -1] ``` Thus we have information for each of the `r nrow(x)` individuals about the membership of Mother, Father, Father's Father, and Mother's Father in the following categories. 1. Cultural elite 2. Professional elite 3. Economic elite 4. Cultural upper-middle 5. Professional upper-middle 6. Economic upper-middle 7. Cultural lower-middle 8. Professional lower-middle 9. Economic lower-middle 10. Skilled workers 11. Unskilled and semi-skilled workers 12. Benefits (family / unempl. / social) 13. Not possible to classify 14. Missing observations ```{r missing, echo = FALSE} RECODE_MISSING <- TRUE REMOVE_MISSING <- TRUE if (RECODE_MISSING) { for (j in 1:4) { x[, j] <- ifelse (is.na (x[, j]), 0, x[, j]) } } if (REMOVE_MISSING) { x <- x[apply(x,1,function (z) all((z > 0) & (z < 12))), ] } ``` In order to simplify our analysis we decided to eliminate all individuals with scores in categories 12-14. This reduces the number of observations to `r nrow(x)`. For the actual sociological analysis this would probably be a wasteful and unwise decision, but for our croissant illustration it makes sense. The marginals of the remaining eleven categories are as follows. ```{r indicators, echo = FALSE} g1 <- ifelse(outer(x[, 1], sort (unique (x[, 1])),"=="), 1L, 0L) g2 <- ifelse(outer(x[, 2], sort (unique (x[, 2])),"=="), 1L, 0L) g3 <- ifelse(outer(x[, 3], sort (unique (x[, 3])),"=="), 1L, 0L) g4 <- ifelse(outer(x[, 4], sort (unique (x[, 4])),"=="), 1L, 0L) gg <- cbind (g1, g2, g3, g4) rbind (colSums (g1), colSums (g2), colSums (g3), colSums (g4)) ``` ##MCA ```{r mca, echo = FALSE, cache = TRUE} library (geigen) burt <- crossprod (cbind(g1, g2, g3, g4)) gburt <- geigen (burt / 4, diag( diag(burt)), symmetric = TRUE) lbd <- rev (gburt$values)[1:41] y <- gburt$vectors[,c (nrow (burt) - 1, nrow (burt) - 2)] y[, 1] <- y[, 1] * sqrt (lbd [2]) y[, 2] <- y[, 2] * sqrt (lbd [3]) y <- y / gburt$vectors[1, nrow (burt)] y1 <- y [1:ncol (g1), ] y2 <- y [ncol (g1) + (1:ncol (g2)), ] y3 <- y [ncol (g1) + ncol (g2) + (1:ncol (g3)), ] y4 <- y [ncol (g1) + ncol (g2) + ncol (g3) + (1:ncol (g4)), ] ``` The first two eigenvalues of an MCA on these data are `r lbd[2]` and `r lbd[3]`. We plot the category quantifications of the four variables. Calling the shape of the plot a croissant is a stretch, it is better to call it a *wedge*.

```{r quantall, echo = FALSE, fig.align = "center"} plot(y, xlab = "dimension 1", ylab = "dimension 2", col = "RED") ```

The first nine categories have factorial structure, combining the three possibilities Cultural, Professional, Economic with the three Elite, Upper Middle, and Lower Middle. We use this factorial structure to connect the nine categories by lines. The line "Cultural"" connects points 1, 4, 7, the line "Elite" connects points 1, 2, 3, and so on. In total there are six lines. Categories 10 and 11, skilled and unskilled labour, are not connected. We'll leave the interpretation to people who know what they are talking about.

```{r quant0, echo = FALSE, fig.align = "center"} par(mfrow = c(2, 2)) plot (y1, type = "n", main = "Mother", xlab = "dimension 1", ylab = "dimension 2") text (y1, as.character (1 : ncol (g1)), col = "RED") lines(y1[c(1,4,7),],col="BLUE") lines(y1[c(2,5,8),],col="BLUE") lines(y1[c(3,6,9),],col="BLUE") lines(y1[c(7,8,9),],col="BLUE") lines(y1[c(4,5,6),],col="BLUE") lines(y1[c(1,2,3),],col="BLUE") plot (y2, type = "n", main = "Father", xlab = "dimension 1", ylab = "dimension 2") text (y2, as.character (1 : ncol (g2)), col = "RED") lines(y2[c(1,4,7),],col="BLUE") lines(y2[c(2,5,8),],col="BLUE") lines(y2[c(3,6,9),],col="BLUE") lines(y2[c(7,8,9),],col="BLUE") lines(y2[c(4,5,6),],col="BLUE") lines(y2[c(1,2,3),],col="BLUE") plot (y3, type = "n", main = "Father's Father", xlab = "dimension 1", ylab = "dimension 2") text (y3, as.character (1 : ncol (g3)), col = "RED") lines(y3[c(1,4,7),],col="BLUE") lines(y3[c(2,5,8),],col="BLUE") lines(y3[c(3,6,9),],col="BLUE") lines(y3[c(7,8,9),],col="BLUE") lines(y3[c(4,5,6),],col="BLUE") lines(y3[c(1,2,3),],col="BLUE") plot (y4, type = "n", main = "Mother's Father", xlab = "dimension 1", ylab = "dimension 2") text (y4, as.character (1 : ncol (g4)), col = "RED") lines(y4[c(1,4,7),],col="BLUE") lines(y4[c(2,5,8),],col="BLUE") lines(y4[c(3,6,9),],col="BLUE") lines(y4[c(7,8,9),],col="BLUE") lines(y4[c(4,5,6),],col="BLUE") lines(y4[c(1,2,3),],col="BLUE") par(mfrow = c(1, 1)) ```

We can replot the separate category quantifications, labeled by the category frequencies. Obviously the sharp point of the wedge corresponds with the categories with the highest frequencies, illustrating the general point that in MCA heavily populated categories tend to be close to the origin of the plot.

```{r frequencies, echo = FALSE, fig.align= "center"} par(mfrow = c(2, 2)) plot (y1, type = "n", main = "Mother", xlab = "dimension 1", ylab = "dimension 2", xlim = c(-2, 3.5)) text (y1, as.character (colSums (g1)), col = "RED") plot (y2, type = "n", main = "Father", xlab = "dimension 1", ylab = "dimension 2", xlim = c(-2, 3.5)) text (y2, as.character (colSums (g2)), col = "RED") plot (y3, type = "n", main = "Father's Father", xlab = "dimension 1", ylab = "dimension 2", xlim = c(-2, 3.5)) text (y3, as.character (colSums (g3)), col = "RED") plot (y4, type = "n", main = "Mother's Father", xlab = "dimension 1", ylab = "dimension 2", xlim = c(-2, 3.5)) text (y4, as.character (colSums (g4)), col = "RED") par(mfrow = c(1, 1)) ```

The `r nrow(x)` object scores form a very sharp and well-filled wedge.

```{r objscores, echo = FALSE, fig.align = "center", cache = TRUE} z <- - qr.Q (qr (gg %*% y)) plot(z, col="RED", main = "object scores", xlab = "dimension 1", ylab = "dimension 2") ```

##Additive Constraints ```{r biggie, echo = FALSE, cache = TRUE} d1 <- rbind(diag(3),diag(3),diag(3)) d2 <- cbind(c(1,1,1,0,0,0,0,0,0),c(0,0,0,1,1,1,0,0,0),c(0,0,0,0,0,0,1,1,1)) d12 <- cbind (d2,d1,rep(0,9),rep(0,9)) d22 <- cbind (matrix (0, 2, 6), diag (2)) dd <- rbind (d12, d22) dbig <- directSum(list(dd,dd,dd,dd)) dbig <- (svd(dbig)$u)[, 1:28] bbig <- crossprod (dbig, burt %*% dbig) ebig <- crossprod (dbig, diag (burt) * dbig) gbig <- geigen (bbig/4, ebig, symmetric = TRUE) lbig <- rev (gbig$values) ybig <- dbig %*% gbig$vectors[, 27:26] y1 <- ybig [1:ncol (g1), ] y2 <- ybig [ncol (g1) + (1:ncol (g2)), ] y3 <- ybig [ncol (g1) + ncol (g2) + (1:ncol (g3)), ] y4 <- ybig [ncol (g1) + ncol (g2) + ncol (g3) + (1:ncol (g1)), ] ``` In figure `r figure_nums("quant0", display = "cite")` we have used the factorial $3\times 3$ structure of the nine first categories to draw a grid in the plot. Instead, we can require that the nine points for each variable are on a regular grid, by using a design matrix that forces the "Cultural Elite" category point to be the sum of a "Cultural" point and an "Elite" point. And so on for all nine combinations. The grids for the four variables are generally different. The MCA under these constraints has first two eigenvalues `r lbig[2]` and `r lbig[3]`, which is only marginally smaller than the unconstrained eigenvalues `r lbd[2]` and `r lbd[3]`. The grid lines clearly display the structure within the wedge.

```{r bigcat, fig.align = "center", echo = FALSE} par(mfrow=c(2, 2)) plot (y1, type = "n", main = "Mother", xlab = "dimension 1", ylab = "dimension 2") text (y1, as.character (1 : ncol (g1)), col = "RED") lines(y1[c(1,4,7),],col="BLUE") lines(y1[c(2,5,8),],col="BLUE") lines(y1[c(3,6,9),],col="BLUE") lines(y1[c(7,8,9),],col="BLUE") lines(y1[c(4,5,6),],col="BLUE") lines(y1[c(1,2,3),],col="BLUE") plot (y2, type = "n", main = "Father", xlab = "dimension 1", ylab = "dimension 2") text (y2, as.character (1 : ncol (g2)), col = "RED") lines(y2[c(1,4,7),],col="BLUE") lines(y2[c(2,5,8),],col="BLUE") lines(y2[c(3,6,9),],col="BLUE") lines(y2[c(7,8,9),],col="BLUE") lines(y2[c(4,5,6),],col="BLUE") lines(y2[c(1,2,3),],col="BLUE") plot (y3, type = "n", main = "Father's Father", xlab = "dimension 1", ylab = "dimension 2") text (y3, as.character (1 : ncol (g3)), col = "RED") lines(y3[c(1,4,7),],col="BLUE") lines(y3[c(2,5,8),],col="BLUE") lines(y3[c(3,6,9),],col="BLUE") lines(y3[c(7,8,9),],col="BLUE") lines(y3[c(4,5,6),],col="BLUE") lines(y3[c(1,2,3),],col="BLUE") plot (y4, type = "n", main = "Mother's Father", xlab = "dimension 1", ylab = "dimension 2") text (y4, as.character (1 : ncol (g4)), col = "RED") lines(y4[c(1,4,7),],col="BLUE") lines(y4[c(2,5,8),],col="BLUE") lines(y4[c(3,6,9),],col="BLUE") lines(y4[c(7,8,9),],col="BLUE") lines(y4[c(4,5,6),],col="BLUE") lines(y4[c(1,2,3),],col="BLUE") par(mfrow=c(1, 1)) ```

It is somewhat surprising that the object scores show a pattern which is quite different from the unconstrained scores in `r figure_nums("objscores", display = "cite")`. The wedge is still there, but it is clearly partitioned into five or six parallel clouds of individuals.

```{r bigobj, fig.align = "center", echo = FALSE} z <- cbind(g1,g2,g3,g4) %*% ybig z[, 1] <- z[, 1] / sqrt (sum (z[, 1] ^ 2)) z[, 2] <- z[, 2] / sqrt (sum (z[, 2] ^ 2)) plot(z, col="RED", main = "object scores", xlab = "dimension 1", ylab = "dimension 2") ```

##Additive and Equality Constraints ```{r small, echo = FALSE, cache = TRUE} dsmall <- rbind(dd,dd,dd,dd) bsmall <- crossprod (dsmall, burt %*% dsmall) esmall <- crossprod (dsmall, diag (burt) * dsmall) gsmall <- geigen (bsmall/4, esmall, symmetric = TRUE) ysmall <- dsmall %*% gsmall$vectors[,c(7,6)] lsmall <- rev (gsmall$values) ysmall[, 2] <- - ysmall[, 2] y1 <- ysmall [1:ncol (g1), ] y2 <- ysmall [ncol (g1) + (1:ncol (g2)), ] y3 <- ysmall [ncol (g1) + ncol (g2) + (1:ncol (g3)), ] y4 <- ysmall [ncol (g1) + ncol (g2) + ncol (g3) + (1:ncol (g1)), ] ``` We can go a step further with constraining the solution by requiring additivity, as before, and in addition that all four grids are the same. The MCA under these constraints has first two eigenvalues `r lsmall[2]` and `r lsmall[3]`, which are somewhat smaller than the previous eigenvalues.

```{r smallcat, fig.align = "center", echo = FALSE} par(mfrow=c(2, 2)) plot (y1, type = "n", main = "Mother", xlab = "dimension 1", ylab = "dimension 2") text (y1, as.character (1 : ncol (g1)), col = "RED") lines(y1[c(1,4,7),],col="BLUE") lines(y1[c(2,5,8),],col="BLUE") lines(y1[c(3,6,9),],col="BLUE") lines(y1[c(7,8,9),],col="BLUE") lines(y1[c(4,5,6),],col="BLUE") lines(y1[c(1,2,3),],col="BLUE") plot (y2, type = "n", main = "Father", xlab = "dimension 1", ylab = "dimension 2") text (y2, as.character (1 : ncol (g2)), col = "RED") lines(y2[c(1,4,7),],col="BLUE") lines(y2[c(2,5,8),],col="BLUE") lines(y2[c(3,6,9),],col="BLUE") lines(y2[c(7,8,9),],col="BLUE") lines(y2[c(4,5,6),],col="BLUE") lines(y2[c(1,2,3),],col="BLUE") plot (y3, type = "n", main = "Father's Father", xlab = "dimension 1", ylab = "dimension 2") text (y3, as.character (1 : ncol (g3)), col = "RED") lines(y3[c(1,4,7),],col="BLUE") lines(y3[c(2,5,8),],col="BLUE") lines(y3[c(3,6,9),],col="BLUE") lines(y3[c(7,8,9),],col="BLUE") lines(y3[c(4,5,6),],col="BLUE") lines(y3[c(1,2,3),],col="BLUE") plot (y4, type = "n", main = "Mother's Father", xlab = "dimension 1", ylab = "dimension 2") text (y4, as.character (1 : ncol (g4)), col = "RED") lines(y4[c(1,4,7),],col="BLUE") lines(y4[c(2,5,8),],col="BLUE") lines(y4[c(3,6,9),],col="BLUE") lines(y4[c(7,8,9),],col="BLUE") lines(y4[c(4,5,6),],col="BLUE") lines(y4[c(1,2,3),],col="BLUE") par(mfrow=c(1, 1)) ```

The object scores show the same clouds within the wedge, but because many more individuals get the same scores there are fewer points and the clouds are far less dense.

```{r smallobj, fig.align = "center", echo = FALSE} z <- cbind(g1,g2,g3,g4) %*% ysmall z[, 1] <- z[, 1] / sqrt (sum (z[, 1] ^ 2)) z[, 2] <- z[, 2] / sqrt (sum (z[, 2] ^ 2)) plot(z, col="RED", main = "object scores", xlab = "dimension 1", ylab = "dimension 2") ```

##Simultaneous Bilinearizing MCA Using additivity does not really get rid of the wedge, but allows us to study it in more detail and take it apart. SBL theory gives us another tool to look at the wedge. ```{r 3stepN, echo = FALSE} d <- 1 / sqrt (diag (burt)) e <- burt * outer (d, d) h <- threeStepApprox (e, rep(11, 4)) ``` For Norway the eigenvalues of $E/4$ are ```{r enorway1, echo = FALSE} cat (formatC (h$eval/4, digits = 4, format = "f", width = 10)) ``` The direct sum of the eigenvalues of the $R_s/m$ is ```{r enorway2, echo = FALSE} isum <- 0 for (i in 1:11) { cat (formatC (h$aval[isum + 1:4]/4, digits = 4, format = "f", width = 10), "\n") isum <- isum + 4 } ``` Note that each row here adds up to one. Also note that the orders of the matrices $R_s$ in this case are `r h$gg`, because all variables have 11 categories. If we sort the eigenvalues of all $R_s/m$ we find ```{r enorway3, echo = FALSE} cat (formatC (sort(h$aval, decreasing = TRUE)/4, format = "f", width = 10, digits = 4)) ``` The two sets of sorted eigenvalues are plotted in `r figure_nums("noreigplot", display = "cite")`. It is pretty obvious they are close.

```{r noreigplot, fig.align = "center", echo = FALSE} par(pty="s") plot (h$eval/4, sort(h$aval, decreasing = TRUE)/4, col = "RED", xlab = "E", ylab = "R", asp = 1, xlim = c(0,1), ylim = c(0,1)) abline (0, 1, col = "BLUE") ```

The assignment problem for the eigenvectors gives the solution ```{r norassign, echo = FALSE} q <- lp.assign(h$cc^2,direction="max")$solution apply (q, 1, which.max) ``` Again we see the familiar pattern: eigenvalues of $E$ are the first eigenvalue of (trivial) $R_1$, followed by the first eigenvalue of $R_2$, of $R_3$, and of $R_4$. The remaining eigenvalue of $R_2$ are actually the smallest ones of the whole sequence, except for the trivial zeroes. This tells us that these data are highly one-dimensional. We see the familiar wedge again in `r figure_nums("3norwaycatm", display = "cite")`, which plots the dominant eigenvectors corresponding to $R_2$ and $R_3$.

```{r 3norwaycatm, fig.align= "center", echo = FALSE} yy <- d * h$avec [, c(5, 9)] plot (yy[, ], col ="RED", xlab = "dimension 1", ylab = "dimension 2") isum = 0 for (i in 1:4) { yv <- yy[isum + 1:11, ] ov <- order (yv[, 1]) lines (yv[ov, ], col = "BLUE") isum <- isum + 11 } ```

#Appendix: Code ```{r showcode, eval = FALSE, tidy = TRUE} discreteNormal <- function (n, m, r, knots) { x <- matrix (rnorm (n * m), n , m) x <- x %*% chol (r) for (j in 1:m) { x[,j] <- apply (outer (x[, j], c (knots[[j]], Inf), ">"), 1, which.min) } return (x) } repList <- function (x, m) { z <- vector ("list", m) return (lapply (z, function (i) x)) } burtTable <- function (x, degrees = rep (-1, ncol (x)), knots = NULL, center = FALSE, standardize = FALSE, orthonormalize = FALSE) { n <- nrow (x) m <- ncol (x) g <- matrix (0, n, 0) l <- rep (0, m) for (j in 1:m) { z <- x[,j] if (degrees[j] < 0) { h <- ifelse (outer (z, sort (unique (z)), "=="), 1, 0) } else { h <- bsplineBasis (z, degrees [j], knots [[j]]) } if (center) { h <- center (h)[, -1] } if (standardize) { h <- standardize (h) } if (orthonormalize) { h <- gs (h)$q } g <- cbind (g, h) l[j] <- ncol (h) } return (list(c = crossprod (g), g = g, ord = l)) } directSum <- function (x) { m <- length (x) nr <- sum (sapply (x, nrow)) nc <- sum (sapply (x, ncol)) z <- matrix (0, nr, nc) kr <- 0 kc <- 0 for (i in 1:m) { ir <- nrow (x[[i]]) ic <- ncol (x[[i]]) z[kr + (1:ir), kc + (1:ic)] <- x[[i]] kr <- kr + ir kc <- kc + ic } return (z) } kplPerm <- function (cc, k) { kl<-unlist (sapply (k, function (i) 1:i)) p <- ifelse (outer (1:sum (k), order (kl), "=="), 1, 0) return (list (pcp = t(p) %*% cc %*% p, perm = p, ord = as.vector (table (kl)))) } makeE <- function (cc, k) { dd <- mInvSqrt (blockSelect (cc, k)) return (dd %*% cc %*% dd) } mInvSqrt <- function (x) { ex <- eigen (x) ew <- abs (ex$values) ev <- ifelse (ew == 0, 0, 1 / sqrt (ew)) ey <- ex$vectors return (ey %*% (ev * t (ey))) } blockSelect <- function (cc, k) { l <- unlist (lapply (1:length (k), function(i) rep (i,k[i]))) return (cc * ifelse (outer (l, l, "=="), 1, 0)) } kplSVD <- function (e, k, eps = 1e-6, itmax = 500, verbose = TRUE, vectors = TRUE) { m <- length (k) sk <- sum (k) ll <- kk <- ww <- diag (sk) itel <- 1 ossq <- 0 klw <- 1 + cumsum (c (0, k))[1:m] kup <- cumsum (k) ind <- lapply (1:m, function (i) klw[i]:kup[i]) for (i in 1:m) kk[ind[[i]],ind[[i]]] <- t (svd (e[ind[[i]], ])$u) kek <- kk %*% e %*% t(kk) for (i in 1:m) for (j in 1:m) ww[ind[[i]],ind[[j]]] <- ifelse(outer(1:k[i],1:k[j],"=="),1,0) repeat { for (l in 1:m) { if (k[l] == 2) next() li <- ind[[l]] for (i in (klw[l] + 1):(kup[l]-1)) for (j in (i+1):kup[l]) { bi <- kek[i,-li] bj <- kek[j,-li] wi <- ww[i,-li] wj <- ww[j,-li] acc <- sum(wi*bi^2)+sum(wj*bj^2) acs <- sum((wi-wj)*bi*bj) ass <- sum(wi*bj^2)+sum(wj*bi^2) u <- eigen(matrix(c(acc,acs,acs,ass),2,2))$vectors[,1] c <- u[1] s <- u[2] kek[-li,i] <- kek[i,-li] <- c*bi+s*bj kek[-li,j] <- kek[j,-li] <- c*bj-s*bi if (vectors) { ki <- kk[i,li]; kj <- kk[j,li] kk[i,li] <- c*ki+s*kj kk[j,li] <- c*kj-s*ki } } } nssq <- sum (ww * kek ^ 2) - sum (diag (kek) ^ 2) if (verbose) cat("Iteration ",formatC(itel,digits=4),"ssq ",formatC(nssq,digits=10,width=15),"\n") if (((nssq - ossq) < eps) || (itel == itmax)) break() itel <- itel + 1 ossq <- nssq } return(list(kek = kek, kk = kk, itel = itel, ssq = nssq)) } inducedR <- function (c, y, k) { m <- length (k) l <- unlist (lapply (1:m, function(i) rep (i,k[i]))) g <- ifelse (outer (l, 1:m, "=="), 1, 0) s <- g * matrix (y, length(y), m) r <- crossprod (s, c %*% s) e <- abs (diag (r)) d <- ifelse (e == 0, 0, e) return (r / sqrt (outer (d, d))) } threeStepApprox <- function (e, k, eps = 1e-6, itmax = 500, verbose = FALSE, vectors = TRUE) { ee <- eigen (e) ev <- ee$vectors ea <- ee$values f <- kplSVD (e, k, eps = eps, itmax = itmax, verbose = verbose, vectors = TRUE) ef <- f$kek kf <- f$kk g <- kplPerm (ef, k) eg <- g$pcp kg <- crossprod (g$perm, kf) gg <- g$ord gl <- cumsum (gg) lg <- 1 + c (0, gl)[-length(gl) - 1] bb <- lapply (1:length(gl), function (i) eigen (eg[lg[i]:gl[i],lg[i]:gl[i]])$vectors) vc <- directSum (bb) eh <- crossprod (vc, eg %*% vc) kh <- crossprod (kg, vc) cc <- crossprod (ev, kh) mc <- apply (cc, 1, function (x) max (abs (x))) wc <- apply (cc, 1, function (x) which.max (abs (x))) yc <- apply (cc, 2, function (x) max (abs (x))) vc <- apply (cc, 2, function (x) which.max (abs (x))) return (list (eval = ea, evec = ev, aval = diag (eh), avec = kh, bval = ef, cval = eg, dval = eh, cc = cc, mc = mc, wc = wc, yc = yc, vc = vc, gg = gg)) } ``` #Appendix: NEWS 0.01 12/04/15 * First working version posted 0.02 12/06/15 * more runs of examples 0.03 12/07/15 * normal examples 1.00 12/08/15 * added text * added constrained examples * first published version 1.01 12/09/15 * GALO text * GALO single * GALO PCA 1.02 12/10/15 * subsections * removed Norway data from repository * short MCA introduction * added PRIMALS section to GALO * added PRINCALS section to GALO 1.03 12/12/15 * code for three-step MCA * GALO analysis by three-step MCA 1.04 12/13/15 * GALO analysis by three-step MCA * intro theory for simultaneous bilinearity 2.00 12/14/15 * GALO analysis with SBL completed * Norway analysis with SBL * That's it for now 2.01 12/16/15 * small detailed SBL example #References