managing clusters of individuals in mca

05_c_svd_ca.Rmd

@@ -277,7 +277,7 @@ Notons $\mathbf{F}$ (respectivement, $\mathbf{G}$) une matrice dont les lignes s
  & \mathbf{D_c^{-1}} \mathbf{B}\boldsymbol\Delta \\
 \end{align*}
 
-## Formules de transition {#05_c_transition_formula}
+## Formules de transition {#c-05-c-transition-formula}
 \begin{align*}
  & \mathbf{P} - \mathbf{r}\mathbf{c}^T = \mathbf{A}\boldsymbol\Delta\mathbf{B}^T \\
 = \{& \} \\

05_d_svd_mca.Rmd

@@ -117,7 +117,7 @@ cor <- sweep(temp, 1, sum_cor, FUN="/")
 
 ### Individus supplémentaires
 
-Le passage par la matrice de Burt ne permet pas de trouver directement les coordonnées principales des individus (seulement celles des variables). Cependant, nous pouvons a posteriori associer à chaque individu ses coordonnées principales. Pour ce faire, il faut utiliser les formules de transition (voir section \@ref(05_c_transition_formula)) qui permettent d'exprimer les coordonnées principales d'un profil ligne supplémentaire (noté $f^*_{is,k}$) comme barycentre des coordonnées standards des profils colonnes pondérés par le profil ligne supplémentaire. Nous notons $b^*_{is,j}$ la représentation de la $j$-ème composante de l'individu supplémentaire $is$ comme ligne supplémentaire du tableau disjonctif complet. La $j$-ème composante du profil ligne de cet individu supplémentaire est : $\frac{b^*_{is,j}}{b^*_{is,+}}$ (où $b^*_{is,+}$ signifie la somme des éléments de la ligne supplémentaire $\mathbf{b^*_{is}}$ au tableau disjonctif complet).
+Le passage par la matrice de Burt ne permet pas de trouver directement les coordonnées principales des individus (seulement celles des variables). Cependant, nous pouvons a posteriori associer à chaque individu ses coordonnées principales. Pour ce faire, il faut utiliser les formules de transition (voir section \@ref(c-05-c-transition-formula)) qui permettent d'exprimer les coordonnées principales d'un profil ligne supplémentaire (noté $f^*_{is,k}$) comme barycentre des coordonnées standards des profils colonnes pondérés par le profil ligne supplémentaire. Nous notons $b^*_{is,j}$ la représentation de la $j$-ème composante de l'individu supplémentaire $is$ comme ligne supplémentaire du tableau disjonctif complet. La $j$-ème composante du profil ligne de cet individu supplémentaire est : $\frac{b^*_{is,j}}{b^*_{is,+}}$ (où $b^*_{is,+}$ signifie la somme des éléments de la ligne supplémentaire $\mathbf{b^*_{is}}$ au tableau disjonctif complet).
 
 Sur notre exemple jouet de la Table \@ref(tab:05-d-mat-Z), considérons le premier individu : $\mathbf{b^*_1} = [0,1,0,1,0,0,1]^T$. Nous le notons $\mathbf{b^*_1}$ car nous faisons comme si nous l'ajoutions comme ligne supplémentaire au tableau de Burt $\mathbf{B}$ de la Table \@ref(tab:05-d-mat-B). D'après les formules de transition, nous avons :
 \[f^*_{1,k} = \sum_{i=1}^J \frac{b^*_{1,i}}{b^*_{1,+}} a_{i,k}\]
@@ -130,7 +130,7 @@ IS <- sweep(Z, 1, apply(Z, 1, sum), FUN = "/")
 FIS <- IS %*% a
 ```
 
-Affichons une carte des deux premiers axes factoriels (voir Figure \@ref(fig:05-d-map-toy-1-2)).
+Affichons une carte des deux premiers axes factoriels avec des points aux lieux des individus (voir Figure \@ref(fig:05-d-map-toy-1-2)).
 ```{r 05-d-map-toy-1-2, fig.width = 6, fig.cap = "Carte selon les axes factoriels 1 (x) et 2 (y)"}
 plot(f[,1], f[,2], type = "n",
      xlab="", ylab="", asp = 1, xaxt = "n", yaxt = "n")
@@ -155,7 +155,8 @@ Zsup <- matrix( c(0,0,1,1,1,1,1,0,
                  c("i1", "i2", "i3", "i4",
                    "i5", "i6", "i7", "i8"),
                  c("js-1", "js-2", "js-3")))
-kbl(Zsup, caption = "Exemple du codage disjonctif d'une variable supplémentaire sur un exemple jouet",
+kbl(Zsup, caption = "Exemple du codage disjonctif d'une variable supplémentaire \
+sur un exemple jouet",
     booktabs = TRUE) %>%
   kable_styling(latex_options = "striped")
 ```
@@ -485,7 +486,7 @@ summary(dat)
 
 ### Variables supplémentaires
 
-Nous considérons `c_house_value` comme variable supplémentaire (c'est-à-dire quelle ne doit pas influencer le calcul des facteurs).
+Nous considérons `c_house_value` comme variable supplémentaire. C'est-à-dire, comme nous l'avons vu au début de ce chapitre, qu'elle ne doit pas influencer le calcul des facteurs.
 
 ### Synthèse des transformations opérées sur le jeu de données
 
@@ -507,10 +508,66 @@ cam <- mca(dat) # cam pour correspondence analysis model
 
 Nous générons la Figure \@ref(fig:05-d-map-1-2) qui est une carte du plan factoriel 1-2.
 
-```{r 05-d-map-1-2, fig.width = 6, fig.cap = "Carte selon les facteurs 1 (x) et 2 (y)"}
+```{r 05-d-map-1-2, fig.width = 7, fig.cap = "Carte selon les facteurs 1 (x) et 2 (y)"}
 plot(cam)
 ```
 
+Observons les corrélations des modalités supplémentaires (i.e., valeurs des habitations) avec les axes factoriels (voir Table \@ref(tab:05-d-mat-cor-sup)).
+
+```{r 05-d-mat-cor-sup}
+kbl( round_preserve_sum(1000 * cam$sjcor)[,1:7],
+     caption = "Corrélations des modalités supplémentaires avec les axes factoriels",
+     booktabs = TRUE ) %>%
+  kable_styling(latex_options = "striped")
+```
+
+Créons une carte des deux premiers axes factoriels en ajoutant la représentation des modalités supplémentaires (voir Figure \@ref(fig:05-d-map-1-2-jsup))
+
+```{r 05-d-map-1-2-jsup, fig.width = 7, fig.cap = "Carte selon les facteurs 1 (x) et 2 (y) avec modalités supplémentaires"}
+plotjsup.mca(cam)
+```
+
+### Clustering des individus
+
+Calculons les corrélations des clusters d'individus avec les différents axes factoriels (voir Tables \@ref(tab:05-d-clust-cor-1-33), \@ref(tab:05-d-clust-cor-34-66) et \@ref(tab:05-d-clust-cor-67-100)).
+
+```{r}
+clstcor <- clstcor.mca(cam)
+```
+
+```{r 05-d-clust-cor-1-33}
+kbl( clstcor[1:33,],
+     caption = "Corrélations des clusters d'individus avec les axes factoriels (1/3)",
+     booktabs = TRUE ) %>%
+  kable_styling(latex_options = "striped")
+```
+
+```{r 05-d-clust-cor-34-66}
+kbl( clstcor[34:66,],
+     caption = "Corrélations des clusters d'individus avec les axes factoriels (2/3)",
+     booktabs = TRUE ) %>%
+  kable_styling(latex_options = "striped")
+```
+
+```{r 05-d-clust-cor-67-100}
+kbl( clstcor[67:100,],
+     caption = "Corrélations des clusters d'individus avec les axes factoriels (3/3)",
+     booktabs = TRUE ) %>%
+  kable_styling(latex_options = "striped")
+```
+
+Créons une carte des deux premiers axes factoriels en ajoutant la représentation des clusters d'individus (voir Figure \@ref(fig:05-d-map-1-2-tot))
+
+```{r 05-d-map-1-2-tot, fig.width = 7, fig.cap = "Carte selon les facteurs 1 (x) et 2 (y) avec modalités supplémentaires et clusters d'individus"}
+plotisupjsup.mca(cam)
+```
+
+Observons par exemple un diagramme en bâtons de la variable catégorielle supplémentaire `c_house_value` pour les individus du cluster 8 (voir Figure \@ref(fig:05-d-bar-clst-8-c-house-value)).
+
+```{r 05-d-bar-clst-8-c-house-value, fig.width = 7, fig.cap = "Barplot de house value pour le cluster 8"}
+barplotClst.mca(cam, dat, 8, 'c_house_value')
+```
+
 ## Annexe code source
 
 ```{r, code=readLines("05_d_svd_mca_code.R"), eval=FALSE}

05_d_svd_mca_code.R

@@ -1,7 +1,25 @@
+# https://www.r-bloggers.com/2016/07/round-values-while-preserve-their-rounded-sum-in-r/
+# E.G.
+# > sum(c(0.333, 0.333, 0.334))
+# [1] 1
+# > sum(round(c(0.333, 0.333, 0.334), 2))
+# [1] 0.99
+# > sum(round_preserve_sum(c(0.333, 0.333, 0.334), 2))
+# [1] 1.00
+round_preserve_sum <- function(x, digits = 0) {
+  up <- 10 ^ digits
+  x <- x * up
+  y <- floor(x)
+  indices <- tail(order(x-y), round(sum(x)) - sum(y))
+  y[indices] <- y[indices] + 1
+  y / up
+}
+
 # `kcuts` cuts variable `x` into `centers` categories with k-means.
 kcuts <-
 function(x, centers)
 {
+  set.seed(1123)
   km <- kmeans(x = x, centers = centers)
   cuts <- unlist(lapply(order(km$centers), function(clustId) {
     min(x[km$cluster == clustId]) }))
@@ -136,7 +154,7 @@ function()
 }
 
 mca <-
-function(dat)
+function(dat, nclst = 100)
 {
   lev_n <- unlist(lapply(dat$dat, nlevels))
   n <- cumsum(lev_n)
@@ -211,12 +229,80 @@ function(dat)
   sum_cor <- apply(temp, 1, sum)
   sjcor <- sweep(temp, 1, sum_cor, FUN="/")
 
+  # clusters of individuals
+  # Keep enough singular values to capture at least 90% of variance
+  nsv <- which(cumsum(sv / sum(sv)) > 0.9)[1]
+  nstart <- 25
+  set.seed(1123)
+  clsti <- kmeans(fsi[,1:nsv], nclst, nstart)
+  # Principal coordinates of clusters' centers
+  fclst <- clsti$centers
+
   r <- list(f=f, ctr=ctr, cor=cor, r=r, sv=sv,
-            fsi=fsi, sicor=sicor, fsj=fsj, sjcor=sjcor)
+            fsi=fsi, sicor=sicor, fsj=fsj, sjcor=sjcor,
+            clsti=clsti)
   class(r) <- "mca"
   return(r)
 }
 
+# compute clusters' correlations with factors
+clstcor.mca <-
+function(o)
+{
+  if(class(o) != "mca") {
+    warning("Object is not of class 'mca'")
+    return(invisible(NULL))
+  }
+
+  nclst <- length(o$clsti$size)
+
+  # Correlation of clusters and factorial axes
+  temp <- o$clsti$centers^2
+  sum_cor <- apply(temp, 1, sum)
+  clstcor <- sweep(temp, 1, sum_cor, FUN="/")
+
+  nfact <- 3
+
+  # select first nfact most correlated factors for each clst
+  selMostCorFact <- function(i) {
+    temp1 <- (sort(clstcor[i,], decreasing = TRUE))[1:nfact]
+    temp2 <- round_preserve_sum(1000 * temp1)
+    temp3 <- strtoi(sub('.', '', names(temp1)))
+    rbind(temp3, temp2)
+  }
+
+  tblClstCor <- t(sapply(1:nclst, selMostCorFact))
+  tblClstCor <- cbind(tblClstCor, o$clsti$size)
+  rwithinss <- o$clsti$withinss / o$clsti$size # withinss relative to the cluster size
+  clstqlty <- round_preserve_sum(1000 * rwithinss / sum(rwithinss))
+  tblClstCor <- cbind(tblClstCor, clstqlty)
+  rownames(tblClstCor) <- paste0('clst-', 1:nclst)
+  colnames(tblClstCor) <- c( rep(c("F", "COR"), nfact), "SIZE", "QLTY" )
+
+  tblClstCor[order(tblClstCor[,1], tblClstCor[,3], tblClstCor[,5]),]
+}
+
+textSize.mca <-
+function(o, d1 = NULL, d2 = NULL)
+{
+  if(class(o) != "mca") {
+    warning("Object is not of class 'mca'")
+    return(invisible(NULL))
+  }
+  if(is.null(d1)) d1<-1
+  if(is.null(d2)) d2<-2
+
+  # Part de l'inertie du plan factoriel d1-d2 expliquée par chaque profil
+  cont <- o$r * (o$f[,d1]^2 + o$f[,d2]^2) / (o$sv[d1]^2 + o$sv[d2]^2)
+  names <- rownames(o$f)
+  names[cont < 0.01] <- "."
+  optimPar <- nonlinearFontSize.mca(cont)
+  sizes <- log(1 + exp(optimPar[1]) * cont^optimPar[2])
+  sizes[cont < 0.01] <- 1
+  r <- list(names=names, sizes=sizes)
+  return(r)
+}
+
 plot.mca <-
 function(o, d1 = NULL, d2 = NULL)
 {
@@ -228,15 +314,105 @@ function(o, d1 = NULL, d2 = NULL)
   if(is.null(d1)) d1<-1
   if(is.null(d2)) d2<-2
 
-  # Part de l'inertie du plan factoriel d1-d2 expliquée par chaque profil
+  ns <- textSize.mca(o, d1, d2) # names and sizes
-  rowcon <- o$r * (o$f[,d1]^2 + o$f[,d2]^2) / (o$sv[d1]^2 + o$sv[d2]^2)
-  rnames <- rownames(o$f)
-  rnames[rowcon < 0.01] <- "."
-  rsize <- log(1 + exp(1) * rowcon^0.3)
-  rsize[rowcon < 0.01] <- 1
 
   plot(o$f[,d1], o$f[,d2], type = "n",
      xlab="", ylab="", asp = 1, xaxt = "n", yaxt = "n")
-  text(o$f[,d1], o$f[,d2], rnames, adj = 0, cex = rsize)
+  text(o$f[,d1], o$f[,d2], ns$names, adj = 0, cex = ns$sizes, col = 'blue', font = 2)
   points(0, 0, pch = 3)
 }
+
+plotjsup.mca <-
+function(o, d1 = NULL, d2 = NULL)
+{
+  if(class(o) != "mca") {
+    warning("Object is not of class 'mca'")
+    return(invisible(NULL))
+  }
+  if(is.null(d1)) d1<-1
+  if(is.null(d2)) d2<-2
+
+  ns <- textSize.mca(o, d1, d2) # names and sizes
+
+  plot(c(o$f[,d1], o$fsj[,d1]), c(o$f[,d2], o$fsj[,d2]), type = "n",
+     xlab="", ylab="", asp = 1, xaxt = "n", yaxt = "n")
+  text(o$f[,d1], o$f[,d2], ns$names, adj = 0, cex = ns$sizes,
+       col = 'blue', font = 2)
+  text(o$fsj[,d1], o$fsj[,d2], rownames(o$fsj), adj = 0,
+       font = 4, cex = 0.7, col = 'red')
+  points(0, 0, pch = 3)
+}
+
+plotisupjsup.mca <-
+function(o, d1 = NULL, d2 = NULL)
+{
+  if(class(o) != "mca") {
+    warning("Object is not of class 'mca'")
+    return(invisible(NULL))
+  }
+  if(is.null(d1)) d1<-1
+  if(is.null(d2)) d2<-2
+
+  nsm <- textSize.mca(o, d1, d2) # names and sizes for modalities
+  nsi <- textSizeClst.mca(o, d1, d2) # names and sizes for clusters of individuals
+
+  plot(c(o$f[,d1], o$fsj[,d1], o$clsti$centers[,d1]),
+       c(o$f[,d2], o$fsj[,d2], o$clsti$centers[,d2]),
+       type = "n", xlab="", ylab="", asp = 1, xaxt = "n", yaxt = "n")
+  text(o$f[,d1], o$f[,d2], nsm$names, adj = 0, cex = nsm$sizes,
+       col = 'blue', font = 2)
+  text(o$fsj[,d1], o$fsj[,d2], rownames(o$fsj), adj = 0, font = 4,
+       cex = 0.7, col = 'red')
+  text(o$clsti$centers[,d1], o$clsti$centers[,d2],
+       nsi$names, adj = 0, cex = nsi$sizes, font = 3)
+  points(0, 0, pch = 3)
+}
+
+textSizeClst.mca <-
+function(o, d1 = NULL, d2 = NULL)
+{
+  if(class(o) != "mca") {
+    warning("Object is not of class 'mca'")
+    return(invisible(NULL))
+  }
+  if(is.null(d1)) d1<-1
+  if(is.null(d2)) d2<-2
+
+  # correlation of each cluster center profile with factors d1+d2
+  temp <- o$clsti$centers^2
+  sum_cor <- apply(temp, 1, sum)
+  clstcor <- sweep(temp, 1, sum_cor, FUN="/")
+  clstcor <- rowSums(clstcor[,c(d1,d2)])
+  names <- names(clstcor)
+  names[clstcor < 0.01] <- "x"
+  optimPar <- nonlinearFontSize.mca(clstcor)
+  sizes <- log(1 + exp(optimPar[1]) * clstcor^optimPar[2])
+  sizes[clstcor < 0.01] <- 1
+  r <- list(names=names, sizes=sizes)
+  return(r)
+}
+
+nonlinearFontSize.mca <-
+function(dat, minFontCex = 0.3, maxFontCex = 1)
+{
+  # discover a non-linear transform of data that maps to font sizes
+  # the smallest font factor should be cex ~= minFontCex
+  # the largest font factor should be cex ~= maxFontCex
+  optim.err <- function(par) {
+    err <- (minFontCex - log(1 + exp(par[1]) * min(dat)^par[2]))^2
+    err <- err + (maxFontCex - log(1 + exp(par[1]) * max(dat)^par[2]))^2
+    return(err)
+  }
+  optim.res <- optim(par = c(1,0.3), fn = optim.err)
+  return(optim.res$par)
+}
+
+# barplot of categorical variable `cat`
+# for individuals of cluster `clstnb`
+# `o` is the correspondence analysis model
+# `dat` is the dataset
+barplotClst.mca <-
+function(o, dat, clstnb, cat)
+{
+  barplot(table(dat$dat[which(o$clsti$cluster == clstnb), cat]))
+}