```{r points, echo = FALSE, align = "center", fig.dim = c(10,10), cache = TRUE} par(pty="s") layout(matrix(c(1,0,2,4,3,0), 2, 3)) hGraphPlot(small, x, y) hGraphAll(small, x, y) ```

For this arbitrary initial configuration the loss, i.e. the sum of squares of the line lengths, or the sum of the squred distances between objects and the categories they fall in, is equal to `r hLoss (small, x, y)`. #Reciprocal Averaging In Homogeneity Analysis we minimize loss by what is known as *reciprocal averaging* or *alternating least squares*. We alternate two substeps. The first substep improves the category quantifications for a given set of object scores, the second substep improves and normalizes the object scores for a given set of category quantifications, namely those we have just computed in the first substep. Taken together these two substeps are an *iteration*. So each iteration starts with object scores and category quantifications and uses its two substeps to improve both. Each of the two substeps decreases the loss functions, i.e. the total squared length of the lines in the graph plots. The two substeps are both very simple. Let's look at the first one. We compute optimal category quantifications for given object scores by taking the averages (or *centroids*) of the objects scores in each of the categories. The corresponding graph plots are

```{r updateY, echo = FALSE, align = "center", fig.dim = c(10,10), cache = TRUE} y<-hUpdateY(small, x) par(pty="s") layout(matrix(1:3, 1, 3)) hGraphPlot(small, x, y) ```

and the loss has decreased to `r hLoss(small, x, y)`. Note that we have not improved the object scores yet, so they are still their initial configuration, equally spaced on a circle. Also note, in variable 2 for instance, that category quantifications coincide with object scores, and thus contribute zero to the loss, if the object is the only observation in the category. In addition, because category quantifications are averages of objects points, they are in the convex hull of the object points, which means in this figure that they are within the circle. Averaging objects points makes the category quantifications move closer to the origin. The second substep improves the object scores, while keeping the category quantifications in the locations we have just computed in the first substep. The second substaep has itself two substeps, say $2A$ and $2B$. In the substep $2A$ the score of an object for given category quantifications is computed as the average or centroid of the $m$ category quantifications the object is in.

```{r updateXa, echo = FALSE, align = "center", fig.dim = c(10,10), cache = TRUE} z<-hUpdateX(small, y)$z par(pty="s") layout(matrix(1:3, 1, 3)) hGraphPlot(small, z, y) ```

The loss function is down all the way to `r hLoss(small, z, y)`. This is not a proper loss value, however, because the object scores are no longer centered, standardized, and uncorrelated, and that was a Homogeneity Analysis requirement. Substep 1 shrinks the object scores towards the origin by averaging, substep 2A takes the resulting category quantifications and shrinks them more by even more averaging. Thus in substep $2B$ we have to renormalize the object scores such that they are centered, standardized, and uncorrelated. This gives

```{r updateXb, echo = FALSE, align = "center", fig.dim = c(10,10), cache = TRUE} x<-hUpdateX(small, y)$x par(pty="s") layout(matrix(1:3, 1, 3)) hGraphPlot(small, x, y) ```

Loss, which is now the proper loss for a normalized configuration, has decreased to `r hLoss(small, x, y)`. Now that we have new category quantifications and new suitably normalized object scores we can start the next iteration, and again improve both in two substeps. Ultimately, after a certain number of iterations, there is no change any more from one iteration to another, and we have reached the optimal solution. In other words, there is convergence, and our Homogeneity Analysis is finished. Note that the renormalization in step 2B is necessary, because without it both object scores and category quantifications would become smaller and smaller, and converge to the origin. Of course the origin does have loss zero, but it is never a proper description of the data. The optimal graph plots, after the iterations have converged, are

```{r optimum, echo = FALSE, align = "center", fig.dim = c(10,10), cache = TRUE} h <- hOptimum (small) par(pty="s") layout(matrix(c(1,0,2,4,3,0), 2, 3)) hGraphPlot(small, h$x, h$y) hGraphAll(small, h$x, h$y) ```

The minimum loss for these data is `r hLoss(small, h$x, h$y)`. Because the object scores are in deviations from the mean, and the category quantifications are weighted means of object scores it follows that category quantifications are in deviations from the weighted mean, with weights equal to the marginals of the variable. Thus category quantifications for each variable are distributed around the origin. In Homogeneity Analysis the graph plots for individual variables are often called *star plots*, because the optimal category quantification is the centroid of the scores of the objects in the category. Thus it is somewhere in the middle of a bunch of objects, which are connected to it by lines. Thus the subset of the graph for each category is a star graph, and the corresponding plot for the variable with these stars is a star plot. The top three plots in figure `r figure_nums("graph_optim", display = "n")` are examples of such star plots. One could formulate the objective of Homeogeneity Analysis as finding normalized object scores in such a way that the stars (over all categories of all variables) are as small as possible. Or, in yet another formulation, we want to maximize the between-category variation and minimize the within-category variation. Variables with a small star, which is necessarily close to the origin, have poor discrimination power. The average object score for each category is about the same. In general, categories with a large number of observations will have an average close to the average of all observations, and thus they will be close to the origin. And, conversely, categories with a small number of observations will tend to be relatively far from the origin. #Specifics ##Passive, Supplementary, and Constraining Variables Now, going back to the Durent Bend data, we do not have the values of all 5888 sherds on the three variables. The sherds are aggregated over various site/depth combinations and the original data cannot be recovered from the aggregated table. But the framework of Homogeneity Analysis can still be applied by using *equality restrictions* (@vanbuuren_deleeuw_A_92). The only thing added is that we require that sherds in the same site/depth get the same object score. Or, geometrically, all sherds in the same site/depth are mapped into the same point in the joint plot. The Homogeneity Analysis loss function is still minimized in two steps. The first step, updating the category quantifications, is still the same as in Homogeneity Analysis without equality restrictions. The second step, which updates the object scores, now has three substeps instead of two. In substep $2A$ we compute the average of the category quantifications of the categories the sherd is in. In substep $2B$ we replace these tentative object scores for the sherds by the site/depth averages, and in substep $2C$ we normalize the object scores, making them centered, standardized, and uncorrelated. The graph plots on unconstrained Homogeneity Analysis must now be replaced by plots of *valued graphs*. For any variable each site/depth point is now connected to all category points, and the edge connecting the object and category point has a value equal to the number of sherds in the category. ```{r galo, cache = TRUE, echo = FALSE} data(galo) galo[,2]<-as.factor(galo[,2]) galo[,5]<-as.factor(galo[,5]) k <- as.matrix(galo[,1:4]) h <- homals (k) x <- h$objectscores ``` As an example we use the GALO data, which has been used innumerable times before as a Homogeneity Analysis example, and is in the `Gifi` package in R (@mair_deleeuw_17). The GALO data can be used to show the difference between working with aggregated data (over sherds in the same site or students in the same school) and the raw data, which are actually unavailable for Durant Bend. Here is the description of the GALO data in the help file of the package. ``` galo Gifi R Documentation GALO dataset Description The objects (individuals) are 1290 school children in the sixth grade of elementary school in the city of Groningen (Netherlands) in 1959. Usage galo Format Data frame with the five variables Gender, IQ, Advice, SES and School. IQ (original range 60 to 144) has been categorized into 9 ordered categories and the schools are enumerated from 1 to 37. SES: LoWC = Lower white collar; MidWC = Middle white collar; Prof = Professional, Managers; Shop = Shopkeepers; Skil = Schooled labor; Unsk = Unskilled labor. Advice: Agr = Agricultural; Ext = Extended primary education; Gen = General; Grls = Secondary school for girls; Man = Manual, including housekeeping; None = No further education; Uni = Pre-University. References Peschar, J.L. (1975). School, Milieu, Beroep. Groningen: Tjeek Willink. ``` Note that IQ is measured by the GIT (Groningen Intelligence Test) and that Advice refers to the sixth grade teachers advice abut the most appropriate form of secondary educations for the students. We first ignore the school variable, and only analyze the four variables Gender, IQ, Advice, and SES. Separate joint plots for the four variables, with both object scores and category quatifications, are in figure `r figure_nums("joint_galo", display = "n")`. We do not make star plots (by drawing the lines from the object points to the category points they are in) in this case, because 1290 lines in a plot just create a big black blob. The joint plots show a curved one-dimensional solution with good students on the left and poor students on the right. Such curved solutions, sometimes called *horseshoes*, are a familiar outcome of Homogeneity Analysis when there is a dominant single dimension explaining the results (in this case student achievement). Both IQ and Advice differentiate students well (mainly because teachers rely on IQ scores in their advice), which means they will have the smallest stars. Girls tend to be better students than boys, and SES mainly contrasts the two extremes categories PROF and UNSK.

```{r galo_sex, fig.align = "center", echo = FALSE, fig.dim = c(15,15), cache = TRUE} par(mfrow=c(2,2), pty = "s") g <- ifelse(outer(galo[,1], levels(galo[,1]), "=="), 1, 0) y <- crossprod (g, x) / colSums(g) plot(x, col = "RED", main = "Gender") text (y, levels(galo[,1]), col = "BLUE", cex = 3) g<-ifelse(outer(galo[,2], levels(galo[,2]), "=="), 1, 0) y <- crossprod (g, x) / colSums(g) plot(x, col = "RED", main = "IQ") text (y, levels(galo[,2]), col = "BLUE", cex = 3) g<-ifelse(outer(galo[,3], levels(galo[,3]), "=="), 1, 0) y <- crossprod (g, x) / colSums(g) plot(x, col = "RED", main = "Advice") text (y, levels(galo[,3]), col = "BLUE", cex = 3) g<-ifelse(outer(galo[,4], levels(galo[,4]), "=="), 1, 0) y <- crossprod (g, x) / colSums(g) plot(x, col = "RED", main = "SES") text (y, levels(galo[,4]), col = "BLUE", cex = 3) ```

Our first analysis does not use the School variable at all. In the terminology of Gifi School is a *passive variable*. There are a number of different ways in which we can incorporate School. The first is the obvious one: just repeat the Homogeneity Analysis over all five variables and include School. The joint plots of the first four variables are in figure `r figure_nums("joint_galo_school", display = "n")`. ```{r galo_all, cache = TRUE, echo = FALSE, cache = TRUE} k <- as.matrix(galo) h <- homals (k) z <- h$objectscores ```

```{r galo_with_school, fig.align = "center", echo = FALSE, fig.dim = c(15,15), cache = TRUE} par(mfrow=c(2,2), pty = "s") g <- ifelse(outer(galo[,1], levels(galo[,1]), "=="), 1, 0) y <- crossprod (g, z) / colSums(g) plot(z, col = "RED", main = "Gender") text (y, levels(galo[,1]), col = "BLUE", cex = 3) g<-ifelse(outer(galo[,2], levels(galo[,2]), "=="), 1, 0) y <- crossprod (g, z) / colSums(g) plot(z, col = "RED", main = "IQ") text (y, levels(galo[,2]), col = "BLUE", cex = 3) g<-ifelse(outer(galo[,3], levels(galo[,3]), "=="), 1, 0) y <- crossprod (g, z) / colSums(g) plot(z, col = "RED", main = "Advice") text (y, levels(galo[,3]), col = "BLUE", cex = 3) g<-ifelse(outer(galo[,4], levels(galo[,4]), "=="), 1, 0) y <- crossprod (g, z) / colSums(g) plot(z, col = "RED", main = "SES") text (y, levels(galo[,4]), col = "BLUE", cex = 3) ```

The horseshoe pattern is still there, but it is less pronounced, mainly because of several outlying students at the bottom of the plot more or less defining the vertical dimension. This is due to including the School variable. The joint plot for the School variable is in figure `r figure_nums("joint_galo_school_sep", display = "n")`. We see the outliers are in school 25, a small school with 11 low-IQ students, possibly some type of special education school. The figure also shows the star for school 25. The actual data for the eleven students are ```{r} galo[galo[,5]=="25",] ```

```{r galo_school_sep, fig.align = "center", echo = FALSE, fig.dim = c(10,10), cache = TRUE} par(pty = "s") g<-ifelse(outer(galo[,5], levels(galo[,5]), "=="), 1, 0) y <- crossprod (g, z) / colSums(g) plot(z, col = "RED", main = "SCHOOL") text (y, levels(galo[,5]), col = "BLUE", cex = 3) for (i in which(galo[,5] == "25")) lines (rbind(z[i,], y[25,])) ```

There is another, and perhaps more interesting, way to incorporate School in our analysis, by using it as what is commonly known as a *supplementary variable*. Such a supplementary variable does not actively enter into the Homogeneity Analysis, but after the analysis of the remaining variables we can compute category quantifications of the supplementary variable as centroids of object scores in the categories. Thus we can make star plots for the passive variables that have not been used in the analysis. This is done in figure `r figure_nums("school_galo_sup", display = "n")`. The horseshoe of object scores does not change from the one in figure `r figure_nums("joint_galo", display = "n")`, but by not including School in the analysis we do not give school 25 the opportunity to dominate the second dimension. It is still true that the same schools (5, 24, 25, 37) perform poorly, and the same schools (4, 17, 31, 32) perform well, but generaly the category quantifications are more well-behaved.

```{r supple, fig.align = "center", echo = FALSE, fig.dim = c(10,10), cache = TRUE} par(pty = "s") h <-ifelse(outer(galo[,5], levels(galo[,5]), "=="), 1, 0) y <- crossprod(h, x) / colSums(h) plot (x, col = "RED", xlab = "", ylab = "", axes = FALSE) box() axis(1, labels = FALSE) axis(2, labels = FALSE) text(y, as.character(1:nrow(y)), col = "BLUE", cex = 2) for (i in which(galo[,5] == "25")) lines (rbind(x[i,], y[25,])) ```

We now repeat the analysis, requiring that students in the same school get the same object score. We treat school as a *constraining variable*, performing a Homogeneity Analysis with equality restrictions on the object scores (@vanbuuren_deleeuw_A_92). In terms of the joint plot we require the stars for the school variable to collapse into single points. Computationally this is easiest to do using the the R package `anacor` (@deleeuw_mair_A_09b). In this constrained Homogeneity Analysis each school gets an object score, and these scores are plotted in figure `r figure_nums("school_galo", display = "n")`. Not surprisingly school 25 is now even more of an outlier, but otherwise schools are dispersed pretty much in the same way as before. For the GALO example this constrained analysis throws away useful information and gives a result which is inferior to the supplementary variable approach. For the archeological data we do not ignore within-site information, because the data are aggregated over sherds in the same site to begin with.

```{r aggregate, fig.align = "center", echo = FALSE, fig.dim = c(10,10), cache = TRUE} galo[,5]<-as.factor(galo[,5]) h <- ifelse(outer(galo[,5], levels(galo[,5]), "=="), 1, 0) g1 <- ifelse(outer(galo[,1], levels(galo[,1]), "=="), 1, 0) g2 <- ifelse(outer(galo[,2], levels(galo[,2]), "=="), 1, 0) g3 <- ifelse(outer(galo[,3], levels(galo[,3]), "=="), 1, 0) g4 <- ifelse(outer(galo[,4], levels(galo[,4]), "=="), 1, 0) f <- crossprod (h, cbind(g1,g2,g3,g4)) m <- anacor (f, ndim = 3) xa <- m$row.scores d <- rowSums(f) e<-diag(crossprod(xa, d*xa)) xa<-xa %*% diag(1 / sqrt (e)) plot(xa, type = "n", xlab = "", ylab = "", axes = FALSE) box() axis(1, labels = FALSE) axis(2, labels = FALSE) text(xa, as.character(1:nrow(xa)), col = "RED", cex = 2) ```

##Binary Variables Besides aggregation, another proerty of the Durant Bend data is that the three variables describing the sherds are binary (CS/Plain, Dark/Light, Thin/Thick). This implies some special properties of the Homogeneity Analysis. We have seen that category quantifications are in deviations from the weighted mean, with the weights equal to the marginal frequencies of the variable. If a variable has only two categories, and our Homogeneity Analysis has two dimensions, that means that the two category quantifications for a variable are on a line through the origin. The direction of the line is determined by the marginals of the variables. What Homogeneity Analysis gives us is how far away from the origin the category quantifications are placed on the line to get the smallest stars. We have said very little so far about the number of dimensions we choose for our Homogeneity Analysis. The default is to choose two, because two-dimensional joint and graph plots are the easiest to look at. The maximum number of dimensions in Homogeneity Analysis, i.e. the number of dimensions that are needed to represent all variation in the data, is equal to the total number of categories minus the number of variables. In the GALO example (without School) that is $2+9+6+7-4=20$ but in the Durant Bend example it is $2+2+2-3=3$. Only three dimensions will capture all variation. #Analysis Durant Bend Data The Durant Bend analysis is an aggregated Homogeneity Analysis of three binary variables, requiring equal object scores for all sherds in the same site/depth. The joint plot is in figure `r figure_nums("durant_joint", display = "n")`. For a discussion of these results we refer to the companion paper by @nance_deleeuw_18. There is not much variation in the category quantifications of the three variables (the lines are rather short). In particular, the averages for light sherds and for dark sherds are very close, indicating not much discriminatiry power for that variable.

```{r anacor, align = "center", echo = FALSE, fig.dim = c(10,10), cache = TRUE} par(pty="s") h <- anacor (sherds) x <- h$row.scores d <- rowSums(sherds) e<-diag(crossprod(x, d*x)) x<-x %*% diag(1 / sqrt (e)) y1 <- crossprod (sherds[, 1:2], x) / colSums(sherds[, 1:2]) y2 <- crossprod (sherds[, 3:4], x) / colSums(sherds[, 3:4]) y3 <- crossprod (sherds[, 5:6], x) / colSums(sherds[, 5:6]) mA <- 1.1 * max (x) mI <- 1.1 * min (x) plot(0, xlim=c(mI,mA), ylim = c (mI, mA), type = "n", xlab = "", ylab = "", axes = FALSE) box() axis(1, labels = FALSE) axis(2, labels = FALSE) text(x[1:5,],row.names(sherds)[1:5], col = "RED", cex = 2) text(x[10:13,], row.names(sherds)[10:13], col = "BLUE", cex = 2) text(x[20:24,], row.names(sherds)[20:24], col = "GREEN", cex = 2) text(x[14:16,], row.names(sherds)[14:16], col = "ORANGE", cex = 2) text(x[6:9,], row.names(sherds)[6:9], col = "MAGENTA", cex = 2) text(x[17:19,], row.names(sherds)[17:20], col = "BLACK", cex = 2) text(x[25:27,], row.names(sherds)[25:27], col = "PURPLE", cex = 2) lines(y1, lwd = 3) text(y1, row.names(y1), cex = 2) lines(y2, lwd = 3) text(y2, row.names(y2), cex = 2) lines(y3, lwd = 3) text(y3, row.names(y3), cex = 2) ```

We have to realize, of course, that there are only three dimensions available to describe our data. This makes it interesting to look at the three-dimensional solution, specifically at light vs dark sherds. We more or less expect each variable to define a dimension, indicating relatively low correlations between the three variables, and consequently not much difference between sites. The summary of the three dimensional Homogeneity Analysis from `anacor` is ```{r 3d, echo = FALSE, cache = TRUE} h <- anacor (sherds, ndim = 3, scaling = c ("standard", "Benzecri")) print(h) ``` The two three-dimensional scatterplots, one for categoiry quantifications and one for sites, are in figures `r figure_nums("durent_3d_col", display = "n")` and `r figure_nums("durent_3d_row", display = "n")`. We see that the third dimension indeed separates the light from the dark.

```{r 3dc, fig.align = "center", fig.dim = c(10,10), echo = FALSE, cache = TRUE} library(scatterplot3d) h <- anacor (sherds, ndim = 3, scaling = c ("standard", "Benzecri")) hc <- h$col.scores hr <- h$row.scores s<-scatterplot3d(hc, pch = 20, type = "n", cex.symbols = 3) s$points3d(rbind(hc[1,],hc[2,]), type = "l", lwd = 2, col = "BLUE") s$points3d(rbind(hc[3,],hc[4,]), type = "l", lwd = 2, col = "BLUE") s$points3d(rbind(hc[5,],hc[6,]), type = "l", lwd = 2, col = "BLUE") text(s$xyz.convert(hc), labels=row.names(hc), col = "RED", cex = 2) ```

```{r 3dr, fig.align = "center", fig.dim = c(10,10), echo = FALSE, cache = TRUE} s<-scatterplot3d(hr, color = "RED", bg = "RED", pch = 20, type = "h", cex.symbols = .1) text(s$xyz.convert(hr), labels=row.names(hr), col = "RED", cex = 1) ```

#References