creation of multiple correspondence analysis chapter

This commit is contained in:
Pierre-Edouard Portier 2023-01-02 21:50:15 +01:00
parent 602ea3fc88
commit 501cde884b
3 changed files with 1005 additions and 0 deletions

405
05_d_svd_mca.R Normal file
View File

@ -0,0 +1,405 @@
## -----------------------------------------------------------------------------
set.seed(1123)
source('05_d_svd_mca_code.R')
## ----05-d-mat-Z---------------------------------------------------------------
Z <- matrix( c(0,1,1,0,0,0,1,1,
1,0,0,1,1,1,0,0,
0,1,0,0,0,0,0,0,
1,0,1,1,0,0,1,0,
0,0,0,0,1,1,0,1,
0,0,1,0,0,0,1,0,
1,1,0,1,1,1,0,1),
nrow = 8, ncol = 7,
dimnames = list(
c("i1", "i2", "i3", "i4",
"i5", "i6", "i7", "i8"),
c("j1-1", "j1-2", "j2-1", "j2-2", "j2-3",
"j3-1", "j3-2")))
kbl(Z, caption = "Exemple jouet d'un tableau disjonctif complet",
booktabs = TRUE) %>%
kable_styling(latex_options = "striped")
## ----05-d-mat-C---------------------------------------------------------------
C = t(Z) %*% Z
kbl(C, caption = "Matrice de Burt pour l'exemple jouet",
booktabs = TRUE) %>%
kable_styling(latex_options = "striped")
## -----------------------------------------------------------------------------
dat <- read.csv(file="data/housing.csv", header=TRUE)
str(dat)
## -----------------------------------------------------------------------------
dat$ocean_proximity <- as.factor(dat$ocean_proximity)
summary(dat)
## -----------------------------------------------------------------------------
hist(dat$longitude)
## -----------------------------------------------------------------------------
cuts <- quantile(dat$longitude, probs = seq(0,1,1/4))
hist(dat$longitude)
abline(v=cuts, lwd=4)
## -----------------------------------------------------------------------------
cuts <- kcuts(x = dat$longitude, centers = 4)
cuts
## -----------------------------------------------------------------------------
hist(dat$longitude)
abline(v=cuts, lwd=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')
summary(dat$c_longitude)
## -----------------------------------------------------------------------------
cuts <- kcuts(x = dat$latitude, centers = 4)
cuts
## -----------------------------------------------------------------------------
hist(dat$latitude)
abline(v=cuts, lwd=3)
## -----------------------------------------------------------------------------
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')
summary(dat$c_latitude)
## -----------------------------------------------------------------------------
hist(dat$housing_median_age)
## -----------------------------------------------------------------------------
table(dat$housing_median_age[dat$housing_median_age>45])
## -----------------------------------------------------------------------------
nb_age_52 <- length(dat$housing_median_age[dat$housing_median_age == 52])
pc_age_52 <- round(100 * (nb_age_52 / dim(dat)[1]))
hist(dat$housing_median_age, breaks = 40)
## -----------------------------------------------------------------------------
quantile(dat$housing_median_age[dat$housing_median_age<52])
## -----------------------------------------------------------------------------
cuts <- c(min(dat$housing_median_age), 15, 25, 35, 51, 52)
hist(dat$housing_median_age, breaks = 40)
abline(v=cuts, lwd=3)
## -----------------------------------------------------------------------------
dat$c_age <- cut(x = dat$housing_median_age, unique(cuts), include.lowest = TRUE)
levels(dat$c_age) <- c('A<=15','A(15,25]','A(25,35]','A(35,51]', 'A=52')
summary(dat$c_age)
## -----------------------------------------------------------------------------
cuts <- quantile(dat$total_rooms)
hist(dat$total_rooms)
abline(v=cuts, lwd=3)
## -----------------------------------------------------------------------------
dat$rooms <- dat$total_rooms / dat$households
nb_rooms_gt_8 <- length(dat$rooms[dat$rooms>8])
pc_rooms_gt_8 <- round(100 * (nb_rooms_gt_8 / dim(dat)[1]))
quantile(dat$rooms)
## -----------------------------------------------------------------------------
cuts <- c(min(dat$rooms), 4, 6, 8, max(dat$rooms))
hist(log10(dat$rooms))
abline(v=log10(cuts), lwd=3)
## -----------------------------------------------------------------------------
dat$c_rooms <- cut(x = dat$rooms, unique(cuts), include.lowest = TRUE)
levels(dat$c_rooms) <- c('R<=4','R(4,6]','R(6,8]', 'R>8')
summary(dat$c_rooms)
## -----------------------------------------------------------------------------
dat$bedrooms <- dat$total_bedrooms / dat$households
quantile(dat$bedrooms, na.rm = TRUE)
## -----------------------------------------------------------------------------
cuts <- c(min(dat$bedrooms, na.rm = TRUE), 1.1, max(dat$bedrooms, na.rm = TRUE))
hist(log10(dat$bedrooms))
abline(v=log10(cuts), lwd=3)
## -----------------------------------------------------------------------------
dat$c_bedrooms <- cut(x = dat$bedrooms, unique(cuts), include.lowest = TRUE)
levels(dat$c_bedrooms) <- c('B<=1','B>1')
summary(dat$c_bedrooms)
## -----------------------------------------------------------------------------
dat$pop <- dat$population / dat$households
quantile(dat$pop, probs = seq(0,1,1/4))
## -----------------------------------------------------------------------------
cuts <- c(min(dat$pop), 2, 3, 4, max(dat$pop))
hist(log10(dat$pop))
abline(v=log10(cuts), lwd=3)
## -----------------------------------------------------------------------------
dat$c_pop <- cut(x = dat$pop, unique(cuts), include.lowest = TRUE)
levels(dat$c_pop) <- c('P<=2','P(2,3]', 'P(3,4]', 'P>4')
summary(dat$c_pop)
## -----------------------------------------------------------------------------
quantile(dat$households, probs = seq(0,1,1/4))
## -----------------------------------------------------------------------------
cuts <- c(min(dat$households), 300, 400, 600, max(dat$households))
hist(log10(dat$households))
abline(v=log10(cuts), lwd=3)
## -----------------------------------------------------------------------------
dat$c_households <- cut(x = dat$households, cuts, include.lowest = TRUE)
levels(dat$c_households) <- c('H<=3', 'H(3,4]', 'H(4,6]', 'H>6')
summary(dat$c_households)
## -----------------------------------------------------------------------------
cuts <- quantile(dat$median_income, probs = seq(0,1,1/4))
hist(dat$median_income)
abline(v=cuts, lwd=3)
## -----------------------------------------------------------------------------
table(dat$median_income[dat$median_income>14])
## -----------------------------------------------------------------------------
nb_income_gt_15 <- length(dat$median_income[dat$median_income > 15])
pc_income_gt_15 <- round(100 * (nb_income_gt_15 / dim(dat)[1]), digits = 2)
hist(dat$median_income, breaks = 40)
## -----------------------------------------------------------------------------
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('IL', 'IML', 'IMH', 'IH', 'I>15')
summary(dat$c_income)
## -----------------------------------------------------------------------------
hist(dat$median_house_value, breaks = 30)
## -----------------------------------------------------------------------------
loc_mhv_gt_50k <- dat$median_house_value > 500000
nb_house_value_gt_50k <- length(dat$median_house_value[loc_mhv_gt_50k])
pc_house_value_gt_50k <- round(100 * (nb_house_value_gt_50k / dim(dat)[1]),
digits = 2)
table(dat$median_house_value[dat$median_house_value>499000])
## -----------------------------------------------------------------------------
quantile(dat$median_house_value[dat$median_house_value < 500000])
## -----------------------------------------------------------------------------
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('V<=115', 'V(115,175]', 'V(175,250]', 'V(250,500]',
'V>500')
summary(dat$c_house_value)
## -----------------------------------------------------------------------------
hist(dat$median_house_value)
abline(v=cuts, lwd=3)
## -----------------------------------------------------------------------------
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')]
summary(dat)
## -----------------------------------------------------------------------------
c_rooms_sup <- ventilate(dat$c_rooms, "R>8")
dat$c_rooms[c_rooms_sup$sup_i] <- c_rooms_sup$smpl
## -----------------------------------------------------------------------------
c_income_sup <- ventilate(dat$c_income, "I>15")
dat$c_income[c_income_sup$sup_i] <- c_income_sup$smpl
ocean_proximity_sup <- ventilate(dat$ocean_proximity, "ISLAND")
dat$ocean_proximity[ocean_proximity_sup$sup_i] <- ocean_proximity_sup$smpl
c_bedrooms_sup <- ventilate(dat$c_bedrooms, "NA")
dat$c_bedrooms[c_bedrooms_sup$sup_i] <- c_bedrooms_sup$smpl
dat <- droplevels(dat)
summary(dat)
## -----------------------------------------------------------------------------
# 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),]
## -----------------------------------------------------------------------------
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)]
## -----------------------------------------------------------------------------
B <- t(Z_act) %*% Z_act
B[1:5, 1:5]
## -----------------------------------------------------------------------------
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)]
## -----------------------------------------------------------------------------
K <- J - Q
a <- sweep(dec$vectors, 1, sqrt(r), FUN = "/")
a <- a[,(1 : K)]
f <- a %*% diag(delt)
## -----------------------------------------------------------------------------
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
## -----------------------------------------------------------------------------
temp <- sweep(f^2, 1, r, FUN = "*")
sum_ctr <- apply(temp, 2, sum)
ctr <- sweep(temp, 2, sum_ctr, FUN = "/")
## -----------------------------------------------------------------------------
temp <- f^2
sum_cor <- apply(temp, 1, sum)
cor <- sweep(temp, 1, sum_cor, FUN="/")
## -----------------------------------------------------------------------------
# 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
## ----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)

565
05_d_svd_mca.Rmd Normal file
View File

@ -0,0 +1,565 @@
# SVD et analyse exploratoire d'un jeu de données
```{r}
set.seed(1123)
source('05_d_svd_mca_code.R')
```
## Analyse des correspondances multiples
### Matrice de Burt
Soit un tableau $\mathbf{N}$ qui croise un ensemble d'individus $\{i\}_{i=1}^I$ et un ensemble de variables $\{j\}_{j=1}^J$. Chaque individu $\mathbf{n_i}$ est décrit par $J$ valeurs : $n_{i,1},\dots,n_{i,J}$.
\[
\mathbf{N} =
\begin{array}{c|cccc}
& j_1 & j_2 & \dots & j_J \\
\hline
\mathbf{i_1} & n_{1,1} & n_{1,2} & \dots & n_{1,J} \\
\mathbf{i_2} & n_{2,1} & n_{2,2} & \dots & n_{2,J} \\
\dots & \dots & \dots & \dots & \dots \\
\mathbf{i_I} & n_{I,1} & n_{I,2} & \dots & n_{I,J} \\
\end{array}
\]
Nous considérons la situation où chaque variable $j$ est catégorielle. Mettons que la variable $j_1$ ait 2 modalités, la variable $j_2$ ait 3 modalités,\dots, la variable $j_J$ ait 2 modalités. Nous constuisons le tableau binaire $\mathbf{Z}$, nommé \emph{tableau disjonctif complet}.
\[
\mathbf{Z} =
\begin{array}{c|cc|ccc|c|cc}
& \multicolumn{2}{c|}{j_1} & \multicolumn{3}{c|}{j_2} & \dots & \multicolumn{2}{c|}{j_J} \\
& 1 & 2 & 1 & 2 & 3 & \dots & 1 & 2 \\
\hline
\mathbf{i_1} & 0 & 1 & 1 & 0 & 0 & \dots & 1 & 0 \\
\mathbf{i_2} & 0 & 1 & 0 & 0 & 1 & \dots & 0 & 1 \\
\dots & \dots & \dots & \dots & \dots & \dots & \dots & \dots & \dots \\
\mathbf{i_I} & 1 & 0 & 1 & 0 & 0 & \dots & 1 & 0 \\
\end{array}
\]
Nous construisons la \emph{matrice de Burt} $\mathbf{C} = \mathbf{Z}^T\mathbf{Z}$. C'est une matrice symétrique. Les blocs diagonaux sont les histogrammes de chaque variable catégorielle. Les autres blocs sont les tables de contingence entre paires de variables catégorielles. Prenons un petit exemple numérique pour fixer les idées. Nous partons d'un tableau disjonctif complet $\mathbf{Z}$, voir Table \@ref(tab:05-d-mat-Z), puis nous calculons la matrice de Burt $\mathbf{C}$, voir Table \@ref(tab:05-d-mat-C).
```{r 05-d-mat-Z}
Z <- matrix( c(0,1,1,0,0,0,1,1,
1,0,0,1,1,1,0,0,
0,1,0,0,0,0,0,0,
1,0,1,1,0,0,1,0,
0,0,0,0,1,1,0,1,
0,0,1,0,0,0,1,0,
1,1,0,1,1,1,0,1),
nrow = 8, ncol = 7,
dimnames = list(
c("i1", "i2", "i3", "i4",
"i5", "i6", "i7", "i8"),
c("j1-1", "j1-2", "j2-1", "j2-2", "j2-3",
"j3-1", "j3-2")))
kbl(Z, caption = "Exemple jouet d'un tableau disjonctif complet",
booktabs = TRUE) %>%
kable_styling(latex_options = "striped")
```
```{r 05-d-mat-C}
C = t(Z) %*% Z
kbl(C, caption = "Matrice de Burt pour l'exemple jouet",
booktabs = TRUE) %>%
kable_styling(latex_options = "striped")
```
### Analyse des correspondances de la matrice de Burt
L'analyse des correspondances du tableau disjonctif complet $\mathbf{Z}$ ou l'analyse des correspondances de la matrice de Burt $\mathbf{C}$ sont très proches. Les deux analyses diffèrent en ceci que :
* L'analyse de $\mathbf{C}$ ne donne des résultats que pour les catégories (les individus n'apparaissent plus).
* Pour chaque facteur principal, l'inertie associée par l'analyse de $\mathbf{C}$ est le carré de celle associée par l'analyse de $\mathbf{Z}$.
* $\mathbf{C}$ étant une matrice symétrique définie positive, la décomposition en valeurs singulières est identique à une décomposition en valeurs propres.
Nous pouvons reprendre les résultats dérivés au chapitre sur l'analyse des correspondances pour les appliquer à l'analyse de la matrice de Burt.
* $n \triangleq \sum_i \sum_j c_{i,j}$
* $\mathbf{P} \triangleq \left[ p_{i,j} \right] \quad ; \quad p_{i,j} \triangleq c_{i,j} / n$
* $r_i \triangleq \sum_j p_{i,j}$
* $\mathbf{S} \triangleq \left[ s_{i,j} \right] \quad ; \quad s_{i,j} \triangleq \left( p_{i,j} - r_i r_j \right) / \sqrt{r_i r_j}$ (écarts à l'indépendance standardisés)
* $\mathbf{S} = \mathbf{V} \boldsymbol \Lambda \mathbf{V}^T$ avec $\mathbf{V}^T\mathbf{V}=\mathbf{I}$
* Coordonnées standard du facteur $k$ : $a_{i,k} = v_{i,k} / \sqrt{r_i}$
* Coordonnées principales du facteur $k$ : $f_{i,k} = a_{i,k} \lambda_k$
* Contribution de la modalité $i$ à l'inertie du facteur $k$ : $CTR_{i,k} \triangleq r_i f_{i,k}^2 / \sum_i r_i f_{i,k}^2$
* Corrélation de la modalité $i$ avec le facteur $k$ : $COR_{i,k} \triangleq f_{i,k}^2 / \sum_k f_{i,k}^2$
### Individus et variables 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.
## 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.
```{r}
dat <- read.csv(file="data/housing.csv", header=TRUE)
str(dat)
```
## Analyse individuelle des variables et discrétisation
Nous analysons chaque variable du jeu de données. Nous proposons également une discrétisation de chaque variable afin de mener dans un second une analyse des correspondances.
### Proximité de l'océan
La variable `ocean_proximity` est du type `chr`. C'est en fait une variable catégorielle avec des modalités comme "NEAR BAY", "NEAR OCEAN", etc. Nous la transformons explicitement en variable catégorielle, aussi appelée \emph{facteur} dans le langage R.
```{r}
dat$ocean_proximity <- as.factor(dat$ocean_proximity)
summary(dat)
```
### Longitude
La distribution des longitudes apparait comme étant bi-modale, avec un premier pic autour de -122°, près de la côte ouest, et un second pic autour de -118°, plus dans les terres.
```{r}
hist(dat$longitude)
```
Nous pouvons découper la longitude en, par exemple, 4 intervalles d'égales densités.
```{r}
cuts <- quantile(dat$longitude, probs = seq(0,1,1/4))
hist(dat$longitude)
abline(v=cuts, lwd=4)
```
Alternativement, nous pouvons adopter une approche par clustering pour découper la variable en intervalles.
```{r}
cuts <- kcuts(x = dat$longitude, centers = 4)
cuts
```
```{r}
hist(dat$longitude)
abline(v=cuts, lwd=4)
```
Nous proposons de discrétiser les longitudes en une nouvelle variable catégorielle `c_longitude` avec pour modalités `LO-W` (West), `LO-M` (Mid-West), `LO-ME` (Mid-East), `LO-E` (East).
```{r}
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')
summary(dat$c_longitude)
```
Analyser des données peut être vu comme initier un dialogue avec elles et les laisser nous guider vers un modèle. C'est un processus itératif. Par exemple, nous venons, en première intention, de tester un découpage des longitudes en quatre modalités obtenues par l'application d'un algorithme de classification par moyennes mobiles (k-means). Nous pourrions, plus tard, revenir sur ce choix.
### Latitude
Nous procédons similairement pour les latitudes. Nous les discrétisons en une nouvelle variable catégorielle nommée `c_latitude` avec pour modalités `LA-S` (South), `LA-MS` (Mid-South), `LA-MN` (Mid-North) and `LA-N` (North).
```{r}
cuts <- kcuts(x = dat$latitude, centers = 4)
cuts
```
```{r}
hist(dat$latitude)
abline(v=cuts, lwd=3)
```
```{r}
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')
summary(dat$c_latitude)
```
### Âge médian des habitations
L'histogramme de la variable `housing_median_age` parait unimodal.
```{r}
hist(dat$housing_median_age)
```
Les valeurs supérieures à 50 semblent étonnamment nombreuses. Nous comptons les valeurs différentes supérieures à 45 pour mieux comprendre ce phénomène.
```{r}
table(dat$housing_median_age[dat$housing_median_age>45])
```
Une limite semble avoir été fixée à la valeur 52. Elle apparait clairement sur un histogramme avec plus d'intervalles. Ce phénomène pourrait avoir un impact sur les analyses à venir.
```{r}
nb_age_52 <- length(dat$housing_median_age[dat$housing_median_age == 52])
pc_age_52 <- round(100 * (nb_age_52 / dim(dat)[1]))
hist(dat$housing_median_age, breaks = 40)
```
`r pc_age_52`% des observations ont un âge médian des habitations égal à 52. C'est sans doute suffisant pour qu'il soit intéressant de les regrouper dans une catégorie à part.
```{r}
quantile(dat$housing_median_age[dat$housing_median_age<52])
```
Nous sommes guidés par l'analyse des quantiles pour discrétiser cette variable en une nouvelle variable catégorielle nommée `c_age` avec des modalités dont les frontières correspondent à des nombres entiers faciles à saisir.
```{r}
cuts <- c(min(dat$housing_median_age), 15, 25, 35, 51, 52)
hist(dat$housing_median_age, breaks = 40)
abline(v=cuts, lwd=3)
```
```{r}
dat$c_age <- cut(x = dat$housing_median_age, unique(cuts), include.lowest = TRUE)
levels(dat$c_age) <- c('A<=15','A(15,25]','A(25,35]','A(35,51]', 'A=52')
summary(dat$c_age)
```
### Nombre de pièces
Sur l'histogramme de la variable `total_rooms`, nous observons des valeurs élevées peu fréquentes (l'histogramme est dit à longue queue). Par ailleurs, les valeurs, qui s'étendent de `r min(dat$total_rooms)` à `r format(max(dat$total_rooms), scientific = FALSE)`, sont difficiles à interpréter puisqu'elles dépendent du nombre de foyers.
```{r}
cuts <- quantile(dat$total_rooms)
hist(dat$total_rooms)
abline(v=cuts, lwd=3)
```
Nous créons une nouvelle variable `rooms` pour compter le nombre relatif de pièces par foyers.
```{r}
dat$rooms <- dat$total_rooms / dat$households
nb_rooms_gt_8 <- length(dat$rooms[dat$rooms>8])
pc_rooms_gt_8 <- round(100 * (nb_rooms_gt_8 / dim(dat)[1]))
quantile(dat$rooms)
```
Nous sommes guidés par l'analyse des quantiles pour discrétiser cette variable en une nouvelle variable catégorielle nommée `c_rooms` avec pour modalités : `R<=4` moins de 4 pièces, `R(4,6]` entre 4 et 6 pièces, `R(6,8]` entre 6 et 8 pièces, `R>8` plus de 8 pièces. Nous ajoutons cette dernière catégorie, `R>8`, pour prendre en compte la longue queue de la distribution : seulement `r pc_rooms_gt_8`% des observations appartiennent à cette catégorie.
Pour mieux visualiser les coupures, nous affichons un histogramme après une transformation logarithmique des nombres de pièces.
```{r}
cuts <- c(min(dat$rooms), 4, 6, 8, max(dat$rooms))
hist(log10(dat$rooms))
abline(v=log10(cuts), lwd=3)
```
```{r}
dat$c_rooms <- cut(x = dat$rooms, unique(cuts), include.lowest = TRUE)
levels(dat$c_rooms) <- c('R<=4','R(4,6]','R(6,8]', 'R>8')
summary(dat$c_rooms)
```
### Chambres
Comme pour les pièces, nous créons une nouvelle variable `bedrooms` pour compter le nombre relatif de chambres par foyers. Nous remarquons par ailleurs que cette variable possède des valeurs manquantes indiquées par le symbole `NA`. Nous devons demander à certaines fonctions du langage R d'ignorer les valeurs manquantes. Pour ce faire, nous fixons le paramètre `na.rm` à `TRUE`.
```{r}
dat$bedrooms <- dat$total_bedrooms / dat$households
quantile(dat$bedrooms, na.rm = TRUE)
```
Nous sommes guidés par l'analyse des quantiles pour discrétiser cette variable en une nouvelle variable catégorielle nommée `c_bedrooms` avec pour modalités : `B<=1` 1 chambre ou moins, `B>1` plus de 1 chambre .
```{r}
cuts <- c(min(dat$bedrooms, na.rm = TRUE), 1.1, max(dat$bedrooms, na.rm = TRUE))
hist(log10(dat$bedrooms))
abline(v=log10(cuts), lwd=3)
```
```{r}
dat$c_bedrooms <- cut(x = dat$bedrooms, unique(cuts), include.lowest = TRUE)
levels(dat$c_bedrooms) <- c('B<=1','B>1')
summary(dat$c_bedrooms)
```
### Population
Comme pour les variables comptant les pièces et les chambres, nous créons une nouvelle variable `pop` pour compter le nombre de personnes par foyers.
```{r}
dat$pop <- dat$population / dat$households
quantile(dat$pop, probs = seq(0,1,1/4))
```
Nous sommes guidés par l'analyse des quantiles pour discrétiser cette variable en une nouvelle variable catégorielle nommée `c_pop` avec pour modalités : `P<=2`, `P(2,3]`, `P(3,4]` et `P>4`.
```{r}
cuts <- c(min(dat$pop), 2, 3, 4, max(dat$pop))
hist(log10(dat$pop))
abline(v=log10(cuts), lwd=3)
```
```{r}
dat$c_pop <- cut(x = dat$pop, unique(cuts), include.lowest = TRUE)
levels(dat$c_pop) <- c('P<=2','P(2,3]', 'P(3,4]', 'P>4')
summary(dat$c_pop)
```
### Nombre de foyers
Après une analyse des quantiles, nous créons une variable catégorielle `c_households` avec pour modalités : moins de 300, entre 300 et 400, entre 400 et 600, plus de 600.
```{r}
quantile(dat$households, probs = seq(0,1,1/4))
```
```{r}
cuts <- c(min(dat$households), 300, 400, 600, max(dat$households))
hist(log10(dat$households))
abline(v=log10(cuts), lwd=3)
```
```{r}
dat$c_households <- cut(x = dat$households, cuts, include.lowest = TRUE)
levels(dat$c_households) <- c('H<=3', 'H(3,4]', 'H(4,6]', 'H>6')
summary(dat$c_households)
```
### Revenu médian
Nous remarquons un nombre étonnamment important de revenus supérieurs à 15.
```{r}
cuts <- quantile(dat$median_income, probs = seq(0,1,1/4))
hist(dat$median_income)
abline(v=cuts, lwd=3)
```
Nous pouvons compter le nombre de valeurs distinctes supérieures à 14 pour mieux mesurer ce phénomène.
```{r}
table(dat$median_income[dat$median_income>14])
```
Nous observons une limite, sans doute artificielle, pour les revenus plus grand que 15. Elle apparait clairement sur un histogramme aux intervalles plus fins. Ce phénomène pourrait impacter les analyses à venir.
```{r}
nb_income_gt_15 <- length(dat$median_income[dat$median_income > 15])
pc_income_gt_15 <- round(100 * (nb_income_gt_15 / dim(dat)[1]), digits = 2)
hist(dat$median_income, breaks = 40)
```
Seules `r pc_income_gt_15`% des observations correspondent à un revenu médian supérieur à 15. Nous les isolons cependant dans une catégorie à part, elle pourrait nous être utile plus tard.
```{r}
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('IL', 'IML', 'IMH', 'IH', 'I>15')
summary(dat$c_income)
```
### Prix médian des habitations
Les prix médians supérieurs à 50000 semblent avoir été ramenés à cette valeur limite.
```{r}
hist(dat$median_house_value, breaks = 30)
```
```{r}
loc_mhv_gt_50k <- dat$median_house_value > 500000
nb_house_value_gt_50k <- length(dat$median_house_value[loc_mhv_gt_50k])
pc_house_value_gt_50k <- round(100 * (nb_house_value_gt_50k / dim(dat)[1]),
digits = 2)
table(dat$median_house_value[dat$median_house_value>499000])
```
`r pc_house_value_gt_50k`% des observations ont subi cette simplification.
```{r}
quantile(dat$median_house_value[dat$median_house_value < 500000])
```
Après analyse des quantiles, nous proposons une discrétisation tout en isolant dans une catégorie à part les prix médians supérieurs à 50000.
```{r}
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('V<=115', 'V(115,175]', 'V(175,250]', 'V(250,500]',
'V>500')
summary(dat$c_house_value)
```
```{r}
hist(dat$median_house_value)
abline(v=cuts, lwd=3)
```
## Analyse des correspondances du jeu de données `housing`
### Ventilation des petites modalités
Nous commençons par ventiler les modalités avec de trop petits effectifs, c'est-à-dire que nous les répartissons aléatoirement parmi les autres modalités.
```{r}
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')]
summary(dat)
```
Nous ventilons la modalité `R>8` de la variable `c_rooms`.
```{r}
c_rooms_sup <- ventilate(dat$c_rooms, "R>8")
dat$c_rooms[c_rooms_sup$sup_i] <- c_rooms_sup$smpl
```
De même pour les modalités `I>15` de `c_income`, `ISLAND` de `ocean_proximity` et pour les valeurs manquantes de `c_bedrooms`.
```{r}
c_income_sup <- ventilate(dat$c_income, "I>15")
dat$c_income[c_income_sup$sup_i] <- c_income_sup$smpl
ocean_proximity_sup <- ventilate(dat$ocean_proximity, "ISLAND")
dat$ocean_proximity[ocean_proximity_sup$sup_i] <- ocean_proximity_sup$smpl
c_bedrooms_sup <- ventilate(dat$c_bedrooms, "NA")
dat$c_bedrooms[c_bedrooms_sup$sup_i] <- c_bedrooms_sup$smpl
dat <- droplevels(dat)
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).
```{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),]
```
### 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.
```{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="/")
```
### 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)
```
## Annexe code source
```{r, code=readLines("05_d_svd_mca_code.R"), eval=FALSE}
```

35
05_d_svd_mca_code.R Normal file
View File

@ -0,0 +1,35 @@
# `kcuts` cuts variable `x` into `centers` categories with k-means.
kcuts <- function(x, centers)
{
km <- kmeans(x = x, centers = centers)
cuts <- unlist(lapply(order(km$centers), function(clustId) {
min(x[km$cluster == clustId]) }))
cuts <- c(cuts, max(x))
}
# ventilate modality `mod` of categorical var `cat` of dataframe `dat`
# into the remaining modalities according to their relative frequencies.
ventilate <- function(cat, mod)
{
if (mod == "NA") isna <- TRUE else isna <- FALSE
tab <- table(cat)
if (isna)
{
sup_i <- which(is.na(cat))
sup_n <- length(sup_i)
act_mod <- tab
}
else
{
sup_i <- which(cat == mod)
sup_n <- tab[mod]
act_mod <- tab[! names(tab) %in% mod]
}
prob <- act_mod / sum(act_mod)
smpl <- sample(names(prob), size = sup_n, replace = TRUE, prob = prob)
r <- list(sup_mod = mod,
sup_n = sup_n,
sup_i = sup_i,
smpl = smpl)
return(r)
}