code library for mca

This commit is contained in:
Pierre-Edouard Portier 2023-01-12 23:04:50 +01:00
parent 75ea33ee42
commit e0bf925b55
2 changed files with 257 additions and 151 deletions

View File

@ -141,6 +141,43 @@ points(FIS[,1], FIS[,2], pch = 20)
### Variables supplémentaires
Pour gérer des variables supplémentaires, il suffit de faire la tabulation croisée des nouvelles variables avec les variables actives (c'est-à-dire celles qui ont été utilisées pour calculer les axes factoriels). Notons $\mathbf{Z^*}$ le tableau disjonctif complet des variables supplémentaires (voir Table \@ref(tab:05-d-mat-Z-sup)). Alors, $\mathbf{B^*} \triangleq \mathbf{Z^{*T}}\mathbf{Z}$ produit de nouvelles lignes de la matrice de Burt qui correspondent aux variables supplémentaires. Par les formules de transition, les coordonnées principales de la $s$-ème variable supplémentaire sont :
\[f^*_{s,k} = \sum_{i=1}^J \frac{b^*_{s,i}}{b^*_{s,+}} a_{i,k}\]
Sur notre exemple jouet, mettons que nous ajoutions une variable supplémentaire $js$ à trois modalités dont le codage disjonctif complet est donné par la Table \@ref(tab:05-d-mat-Z-sup).
```{r 05-d-mat-Z-sup}
Zsup <- matrix( c(0,0,1,1,1,1,1,0,
1,0,0,0,0,0,0,0,
0,1,0,0,0,0,0,1),
nrow = 8, ncol = 3,
dimnames = list(
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",
booktabs = TRUE) %>%
kable_styling(latex_options = "striped")
```
Calculons les coordonnées principales de ces modalités supplémentaires.
```{r}
# calcul des profils des variables supplémentaires
JS <- t(Zsup) %*% Z
JS <- sweep(JS, 1, apply(JS, 1, sum), FUN = "/")
# calcul des coordonnées principales des variables supplémentaires
FJS <- JS %*% a
```
Affichons une carte des deux premiers axes factoriels avec les variables supplémentaires (voir Figure \@ref(fig:05-d-map-toy-1-2-vsup)).
```{r 05-d-map-toy-1-2-vsup, fig.width = 6, fig.cap = "Carte selon les axes factoriels 1 (x) et 2 (y)"}
plot(c(f[,1], FJS[,1]), c(f[,2], FJS[,2]), type = "n",
xlab="", ylab="", asp = 1, xaxt = "n", yaxt = "n")
text(f[,1], f[,2], rownames(f), adj = 0)
text(FJS[,1], FJS[,2], colnames(Zsup), adj = 0)
points(0, 0, pch = 3)
```
## Jeu de données
`housing` est un jeu de données célèbre aux nombreuses vertues pédagogiques^[https://www.kaggle.com/datasets/harrywang/housing]. Il permet d'expérimenter sur un problème de régression réaliste, viz. prédire la valeur médiane d'une maison en fonction des caractéristiques de son quartier. Nous allons faire une analyse exploratoire de ce jeu de données par analyse des correspondances.
@ -449,167 +486,29 @@ 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).
### Synthèse des transformations opérées sur le jeu de données
Nous reprenons dans une fonction `makeDatasetHousing` (voir code source en annexe de ce chapitre) toutes les transformations opérées sur le jeu de données.
```{r}
# les catégories supplémentaires doivent être en dernières positions
# du tableau de données
sup_ind <- which(names(dat) == "c_house_value")
dat_act <- dat[,-sup_ind]
dat_sup <- dat[,sup_ind]
I <- dim(dat_act)[1]
Q <- dim(dat_act)[2]
dat[c(1:5, I),]
dat <- makeDatasetHousing()
```
### Tableau disjonctif complet
### Analyse des correspondances du tableau disjonctif complet
Nous construisons le tableau disjonctif complet.
```{r}
lev_n <- unlist(lapply(dat, nlevels))
n <- cumsum(lev_n)
J_t <- sum(lev_n)
Q_t <- dim(dat)[2]
Z <- matrix(0, nrow = I, ncol = J_t)
numdat <- lapply(dat, as.numeric)
offset <- c(0, n[-length(n)])
for (i in 1:Q_t)
Z[1:I + (I * (offset[i] + numdat[[i]] - 1))] <- 1
cn <- rep(names(dat), lev_n)
ln <- unlist(lapply(dat, levels))
dimnames(Z)[[1]] <- as.character(1:I)
dimnames(Z)[[2]] <- paste(cn, ln, sep = "")
Z_sup_min <- n[sup_ind[1] - 1] + 1
Z_sup_max <- n[sup_ind[length(sup_ind)]]
Z_sup_ind <- Z_sup_min : Z_sup_max
Z_act <- Z[,-Z_sup_ind]
J <- dim(Z_act)[2]
# A titre d'illustration, quelques lignes et colonnes de Z_act
Z_act[c(1:5, I), c(1,2,J)]
```
### Matrice de Burt
Nous construisons la matrice de Burt des modalités actives.
```{r}
B <- t(Z_act) %*% Z_act
B[1:5, 1:5]
```
### Décomposition de l'inertie des profils
Nous calculons la décomposition de l'inertie des profils à partir de la matrice de Burt.
```{r}
P <- B / sum(B)
r <- apply(P, 2, sum)
rr <- r %*% t(r)
S <- (P - rr) / sqrt(rr)
dec <- eigen(S)
# les Q dernières valeurs propres sont nécessairement nulles.
delt <- dec$values[1 : (J-Q)]
```
### Coordonnées standard et principales
Nous calculons les coordonnées principales des profils des modalités.
La fonction `mca` (voir code source en annexe de ce chapitre) reprend la méthode introduite en début de chapitre pour calculer l'analyse des correpondances du tableau disjonctif complet.
```{r}
K <- J - Q
a <- sweep(dec$vectors, 1, sqrt(r), FUN = "/")
a <- a[,(1 : K)]
f <- a %*% diag(delt)
```
### Labels
Nous associons des labels aux coordonnées principales et standard.
```{r}
lbl_dic <- c(
"O:<1H" = "ocean_proximity<1H OCEAN",
"O:INL" = "ocean_proximityINLAND",
"O:NB" = "ocean_proximityNEAR BAY",
"O:NO" = "ocean_proximityNEAR OCEAN",
"LO:W" = "c_longitudeLO-W",
"LO:MW" = "c_longitudeLO-MW",
"LO:ME" = "c_longitudeLO-ME",
"LO:E" = "c_longitudeLO-E",
"LA:S" = "c_latitudeLA-S",
"LA:MS" = "c_latitudeLA-MS",
"LA:MN" = "c_latitudeLA-MN",
"LA:N" = "c_latitudeLA-N",
"AG:15]" = "c_ageA<=15",
"AG:25]" = "c_ageA(15,25]",
"AG:35]" = "c_ageA(25,35]",
"AG:51]" = "c_ageA(35,51]",
"AG:52" = "c_ageA=52",
"RO:4]" = "c_roomsR<=4",
"RO:6]" = "c_roomsR(4,6]",
"RO:8]" = "c_roomsR(6,8]",
"BE:1]" = "c_bedroomsB<=1",
"BE:>1" = "c_bedroomsB>1",
"PO:2]" = "c_popP<=2",
"PO:3]" = "c_popP(2,3]",
"PO:4]" = "c_popP(3,4]",
"PO:>4" = "c_popP>4",
"HH:3]" = "c_householdsH<=3",
"HH:4]" = "c_householdsH(3,4]",
"HH:6]" = "c_householdsH(4,6]",
"HH:>6" = "c_householdsH>6",
"IC:L" = "c_incomeIL",
"IC:ML" = "c_incomeIML",
"IC:MH" = "c_incomeIMH",
"IC:H" = "c_incomeIH",
"HV:A" = "c_house_valueV<=115",
"HV:B" = "c_house_valueV(115,175]",
"HV:C" = "c_house_valueV(175,250]",
"HV:D" = "c_house_valueV(250,500]",
"HV:E" = "c_house_valueV>500"
)
lbl_act_dic <- lbl_dic[1:J]
fac_names <- paste("F", paste(1 : K), sep = "")
rownames(a) <- names(lbl_act_dic)
colnames(a) <- fac_names
rownames(f) <- names(lbl_act_dic)
colnames(f) <- fac_names
```
### Contributions et corrélations
Nous calculons $CTR_{i,k}$, la contribution de la modalité $i$ au facteur $k$.
```{r}
temp <- sweep(f^2, 1, r, FUN = "*")
sum_ctr <- apply(temp, 2, sum)
ctr <- sweep(temp, 2, sum_ctr, FUN = "/")
```
Nous calculons $COR_{i,k}$, la corrélation de la modalité $i$ avec le facteur $k$.
```{r}
temp <- f^2
sum_cor <- apply(temp, 1, sum)
cor <- sweep(temp, 1, sum_cor, FUN="/")
cam <- mca(dat) # cam pour correspondence analysis model
```
### Cartes factorielles
Nous générons la Figure \@ref(fig:05-d-map-1-2) qui est une carte du plan factoriel 1-2.
```{r}
# Parts de l'inertie du plan factoriel 1-2 expliquée par chaque profil
rowcon <- r * (f[,1]^2 + f[,2]^2) / sum(delt[1:2]^2)
rnames <- rownames(f)
rnames[rowcon < 0.01] <- "."
rsize <- log(1 + exp(1) * rowcon^0.3)
rsize[rowcon < 0.01] <- 1
```
```{r 05-d-map-1-2, fig.width = 6, fig.cap = "Carte selon les facteurs 1 (x) et 2 (y)"}
plot(f[,1], f[,2], type = "n",
xlab="", ylab="", asp = 1, xaxt = "n", yaxt = "n")
text(f[,1], f[,2], rnames,
adj = 0, cex = rsize)
points(0, 0, pch = 3)
plot(cam)
```
## Annexe code source

View File

@ -1,5 +1,6 @@
# `kcuts` cuts variable `x` into `centers` categories with k-means.
kcuts <- function(x, centers)
kcuts <-
function(x, centers)
{
km <- kmeans(x = x, centers = centers)
cuts <- unlist(lapply(order(km$centers), function(clustId) {
@ -7,9 +8,10 @@ kcuts <- function(x, centers)
cuts <- c(cuts, max(x))
}
# ventilate modality `mod` of categorical var `cat` of dataframe `dat`
# ventilate modality `mod` of categorical var `cat`
# into the remaining modalities according to their relative frequencies.
ventilate <- function(cat, mod)
ventilate <-
function(cat, mod)
{
if (mod == "NA") isna <- TRUE else isna <- FALSE
tab <- table(cat)
@ -33,3 +35,208 @@ ventilate <- function(cat, mod)
smpl = smpl)
return(r)
}
# pre-processing for exploratory analytics of the housing dataset
makeDatasetHousing <-
function()
{
# chargement du jeu de données
dat <- read.csv(file="data/housing.csv", header=TRUE)
# transformation de ocean_proximity en facteur
dat$ocean_proximity <- as.factor(dat$ocean_proximity)
levels(dat$ocean_proximity)[levels(dat$ocean_proximity)=="<1H OCEAN"] <- "O:<1H"
levels(dat$ocean_proximity)[levels(dat$ocean_proximity)=="ISLAND"] <- "O:ISL"
levels(dat$ocean_proximity)[levels(dat$ocean_proximity)=="INLAND"] <- "O:INL"
levels(dat$ocean_proximity)[levels(dat$ocean_proximity)=="NEAR BAY"] <- "O:NB"
levels(dat$ocean_proximity)[levels(dat$ocean_proximity)=="NEAR OCEAN"] <- "O:NO"
# quantification de longitude
cuts <- kcuts(x = dat$longitude, centers = 4)
dat$c_longitude <- cut(x = dat$longitude, unique(cuts), include.lowest = TRUE)
levels(dat$c_longitude) <- c('LO:W', 'LO:MW', 'LO:ME', 'LO:E')
# quantification de latitude
cuts <- kcuts(x = dat$latitude, centers = 4)
dat$c_latitude <- cut(x = dat$latitude, unique(cuts), include.lowest = TRUE)
levels(dat$c_latitude) <- c('LA:S','LA:MS','LA:MN','LA:N')
# quantification de housing_median_age
cuts <- c(min(dat$housing_median_age), 15, 25, 35, 51, 52)
dat$c_age <- cut(x = dat$housing_median_age, unique(cuts), include.lowest = TRUE)
levels(dat$c_age) <- c('AG:15]','AG:25]','AG:35]','AG:51]', 'AG:52')
# création et quantification de rooms
dat$rooms <- dat$total_rooms / dat$households
cuts <- c(min(dat$rooms), 4, 6, 8, max(dat$rooms))
dat$c_rooms <- cut(x = dat$rooms, unique(cuts), include.lowest = TRUE)
levels(dat$c_rooms) <- c('RO:4]','RO:6]','RO:8]', 'RO:>8')
# création et quantification de bedrooms
dat$bedrooms <- dat$total_bedrooms / dat$households
cuts <- c(min(dat$bedrooms, na.rm = TRUE), 1.1, max(dat$bedrooms, na.rm = TRUE))
dat$c_bedrooms <- cut(x = dat$bedrooms, unique(cuts), include.lowest = TRUE)
levels(dat$c_bedrooms) <- c('BE:1]','BE:>1')
# création et quantification de pop
dat$pop <- dat$population / dat$households
cuts <- c(min(dat$pop), 2, 3, 4, max(dat$pop))
dat$c_pop <- cut(x = dat$pop, unique(cuts), include.lowest = TRUE)
levels(dat$c_pop) <- c('PO:2]','PO:3]', 'PO:4]', 'PO:>4')
# quantification de households
cuts <- c(min(dat$households), 300, 400, 600, max(dat$households))
dat$c_households <- cut(x = dat$households, cuts, include.lowest = TRUE)
levels(dat$c_households) <- c('HH:3]', 'HH:4]', 'HH:6]', 'HH:>6')
# quantification de median_income
cuts <- quantile(dat$median_income, probs = seq(0,1,1/4))
cuts <- c(cuts[1:length(cuts)-1], 15, max(dat$median_income))
dat$c_income <- cut(x = dat$median_income, cuts, include.lowest = TRUE)
levels(dat$c_income) <- c('IC:L', 'IC:ML', 'IC:MH', 'IC:H', 'IC:>15')
# quantification de median_house_value
cuts <- c(min(dat$median_house_value), 115000, 175000, 250000, 500000,
max(dat$median_house_value))
dat$c_house_value <- cut(x = dat$median_house_value, cuts, include.lowest = TRUE)
levels(dat$c_house_value) <- c('HV:A', 'HV:B', 'HV:C', 'HV:D', 'HV:E')
# création du jeu de données quantifié
dat.all <- dat
dat <- dat.all[c('ocean_proximity', 'c_longitude', 'c_latitude', 'c_age',
'c_rooms', 'c_bedrooms', 'c_pop', 'c_households', 'c_income',
'c_house_value')]
vent <- list()
# ventilation de la modalité RO:>8 de c_rooms
vent$c_rooms <- ventilate(dat$c_rooms, "RO:>8")
dat$c_rooms[vent$c_rooms$sup_i] <- vent$c_rooms$smpl
# ventilation de la modalité IC:>15 de c_income
vent$c_income <- ventilate(dat$c_income, "IC:>15")
dat$c_income[vent$c_income$sup_i] <- vent$c_income$smpl
# ventilation de la modalité O:ISL de ocean_proximity
vent$ocean_proximity <- ventilate(dat$ocean_proximity, "O:ISL")
dat$ocean_proximity[vent$ocean_proximity$sup_i] <- vent$ocean_proximity$smpl
# ventilation des valeurs manquantes de c_bedrooms
vent$c_bedrooms <- ventilate(dat$c_bedrooms, "NA")
dat$c_bedrooms[vent$c_bedrooms$sup_i] <- vent$c_bedrooms$smpl
# suppression des modalités vides après ventilation
dat <- droplevels(dat)
# positionnement de c_house_value en variable supplémentaire
supind <- which(names(dat) == "c_house_value")
r <- list( vent=vent, dat=dat, supind=supind )
return(r)
}
mca <-
function(dat)
{
lev_n <- unlist(lapply(dat$dat, nlevels))
n <- cumsum(lev_n)
J_t <- sum(lev_n)
Q_t <- dim(dat$dat)[2]
I <- dim(dat$dat)[1]
# Indicator matrix
Z <- matrix(0, nrow = I, ncol = J_t)
numdat <- lapply(dat$dat, as.numeric)
offset <- c(0, n[-length(n)])
for (i in 1:Q_t)
Z[1:I + (I * (offset[i] + numdat[[i]] - 1))] <- 1
dimnames(Z)[[1]] <- as.character(1:I)
dimnames(Z)[[2]] <- as.character(unlist(lapply(dat$dat, levels)))
Z_sup_min <- n[dat$supind[1] - 1] + 1
Z_sup_max <- n[dat$supind[length(dat$supind)]]
Z_sup_ind <- Z_sup_min : Z_sup_max
Z_act <- Z[,-Z_sup_ind]
J <- dim(Z_act)[2]
Q <- dim(dat$dat)[2] - length(dat$supind)
# Burt matrix
B <- t(Z_act) %*% Z_act
# Decomposition of inertia
P <- B / sum(B)
r <- apply(P, 2, sum)
rr <- r %*% t(r)
S <- (P - rr) / sqrt(rr)
dec <- eigen(S)
# Last Q eigen value must be zero
K <- J - Q
sv <- dec$values[1 : K]
# standard coordinates
a <- sweep(dec$vectors, 1, sqrt(r), FUN = "/")
a <- a[,(1 : K)]
# principal coordinates
f <- a %*% diag(sv)
fac_names <- paste("F", paste(1 : K), sep = "")
rownames(a) <- colnames(Z_act)
colnames(a) <- fac_names
rownames(f) <- colnames(Z_act)
colnames(f) <- fac_names
# contributions
temp <- sweep(f^2, 1, r, FUN = "*")
sum_ctr <- apply(temp, 2, sum)
ctr <- sweep(temp, 2, sum_ctr, FUN = "/")
# correlations
temp <- f^2
sum_cor <- apply(temp, 1, sum)
cor <- sweep(temp, 1, sum_cor, FUN="/")
# profiles of "supplementary" individuals
si <- sweep(Z_act, 1, apply(Z_act, 1, sum), FUN = "/")
fsi <- si %*% a
temp <- fsi^2
sum_cor <- apply(temp, 1, sum)
sicor <- sweep(temp, 1, sum_cor, FUN="/")
# profiles of "supplementary" variables
Z_sup <- Z[,Z_sup_ind]
sj <- t(Z_sup) %*% Z_act
sj <- sweep(sj, 1, apply(sj, 1, sum), FUN = "/")
fsj <- sj %*% a
temp <- fsj^2
sum_cor <- apply(temp, 1, sum)
sjcor <- sweep(temp, 1, sum_cor, FUN="/")
r <- list(f=f, ctr=ctr, cor=cor, r=r, sv=sv,
fsi=fsi, sicor=sicor, fsj=fsj, sjcor=sjcor)
class(r) <- "mca"
return(r)
}
plot.mca <-
function(o, d1 = NULL, d2 = NULL)
{
if(class(o) != "mca") {
warning("Object is not of class 'mca'")
UseMethod("plot")
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
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)
points(0, 0, pch = 3)
}