Introduction

Un ormeau est connu sous plusieurs appellations: ormier, ormet de mer, haliotis, abalone, etc. Dans le langage courant, il désigne essentiellement une espèce de nos côtes françaises: Haliotis tuberculata. Ce mollusque marin, qui ressemble à une grosse moule épaisse, se nourrit d’algues qu’il râpe (coquillage brouteur) sur les rochers. Les ormeaux sont photofuges et leur élevage en ferme aquacole est l’halioticulture.

Ainsi, le but de ce projet est premièrement de prédire l’age des ormeaux en fonction des variables physiologiques étant donné que ’âge de l’ormeau est déterminé en coupant la coquille à travers le cône, en la colorant et en comptant le nombre d’anneaux au microscope, une tâche redondante qui prend du temps.

Dans un second temps, il nous est demandé de creer une regle de classification permettant de discriminer les males des femelles.

Importation des données

Abalone = read.csv("http://math.univ-lille1.fr/~preda/GIS4/PRJ/P11/abalone.data", header=FALSE, sep=",",dec=".")  
colnames(Abalone) <- c("Sex","Lenght","Diameter","Height","Whole weight","Shucked weight","Viscera weight","Shell weight","Rings")
Abalone$Sex<-as.factor(Abalone$Sex)

Chargement des packages

library(corrplot)
library(Hmisc)
library(dplyr)
library(rpart)
library(rpart.plot)
library(MASS)
library(car)
library(ggplot2)
library(gridExtra)
library(FactoMineR)

Statistiques univariées descriptives

dim(Abalone)
## [1] 4177    9
summary(Abalone)
##  Sex          Lenght         Diameter          Height        Whole weight   
##  F:1307   Min.   :0.075   Min.   :0.0550   Min.   :0.0000   Min.   :0.0020  
##  I:1342   1st Qu.:0.450   1st Qu.:0.3500   1st Qu.:0.1150   1st Qu.:0.4415  
##  M:1528   Median :0.545   Median :0.4250   Median :0.1400   Median :0.7995  
##           Mean   :0.524   Mean   :0.4079   Mean   :0.1395   Mean   :0.8287  
##           3rd Qu.:0.615   3rd Qu.:0.4800   3rd Qu.:0.1650   3rd Qu.:1.1530  
##           Max.   :0.815   Max.   :0.6500   Max.   :1.1300   Max.   :2.8255  
##  Shucked weight   Viscera weight    Shell weight        Rings       
##  Min.   :0.0010   Min.   :0.0005   Min.   :0.0015   Min.   : 1.000  
##  1st Qu.:0.1860   1st Qu.:0.0935   1st Qu.:0.1300   1st Qu.: 8.000  
##  Median :0.3360   Median :0.1710   Median :0.2340   Median : 9.000  
##  Mean   :0.3594   Mean   :0.1806   Mean   :0.2388   Mean   : 9.934  
##  3rd Qu.:0.5020   3rd Qu.:0.2530   3rd Qu.:0.3290   3rd Qu.:11.000  
##  Max.   :1.4880   Max.   :0.7600   Max.   :1.0050   Max.   :29.000

Une première analyse de ce jeu de données nous montre qu’on retrouve 4177 observations pour 9 variables :

Sex
Length
Diameter
Height
Whole weight
Shucked weight Viscera weight
Shell weight
Rings

On peut ainsi réaliser quelques analyses élémentaires à l’aide du summary. On remarque ainsi que pour la variable Sex on retrouve 1307 individus femelles, 1528 males et 1342 enfants. On peut aussi commenter la variable Rings, on remarque ainsi que celle ci est comprise entre 1 et 29, avec une moyenne étant égale à 9,934.

Une autre information que l’on peut tirer de ce summary est qu’il n’y aucune données manquantes dans ce jeu de données.

mcor <- cor(Abalone[,2:7])
corrplot(mcor,type="upper",order="hclust",tl.col = "black",tl.srt = 45)

Enfin ce que l’on peut observer ici est que les variables sont toutes très corrélées ce qui risque de nous poser des problèmes plus tard lors de la prédiction de l’âge et de la classification du sexe.

Prédire l’age des ormeaux en fonction des variables physiologiques

library(dplyr)
AbaloneAge = Abalone %>%  mutate(Age=case_when(
  Rings %in% 1:5 ~ "jeunes",
  Rings %in% 6:13 ~ "adultes",
  Rings %in% 14:29 ~ "vieux"
))

Créeons ainsi la variable Age à l’aide de la fonction mutate et de la variable Rings. On sait qu’il suffit de rajouter 1.5 à la variable Rings afin d’avoir l’Age. On décide ici de rendre la variable qualitative afin de réaliser un arbre de décision qui nous permettra de pérdire l’âge des ormeaux en fonction des autres variables. A l’aide de recherche biologiques et du summary on décide donc de répartir l’Age de la facon suivante : * Les individus ayant entre 1 et 5 anneaux seront “jeunes” * Ceux entre 6 et 13 anneaux seront des “adultes” * Enfin, les ormeaux ayant entre 14 et 30 anneaux seront des “vieux”

AbaloneAge = Abalone %>%  mutate(Age=case_when(
  Rings %in% 1:5 ~ "jeunes",
  Rings %in% 6:13 ~ "adultes",
  Rings %in% 14:29 ~ "vieux"
))

Voici quelques graphiques présentant la longueur, le diamètre, … des abalones en fonction de leur âge afin d’avoir une premiere idée de la prédiction possible de l’âge.

# boite à moustaches
par(mfrow = c(4, 2), mar = c(2, 2, 2, 0) + .1)
for (i in 2:8) {
  boxplot(AbaloneAge[,i] ~ AbaloneAge$Age, main = names(Abalone)[i], ylab=names(Abalone)[i], xlab="Age")
}

Tout d’abord, nous remarquons que pour chaque boîte à moustache, celle des jeunes a des valeurs plus basse que celles des aldultes et des vieux. Par exemple, 50% des abalones jeunes ont une longueur inférieure à 0.255, alors que chez 50% des abalones adultes ont une longueur inférieure à 0.54.
De plus, entre les abalones adultes et vieux, nous observons une différence de distribution. Prenons l’exemple précedent, 50% des abalones vieux ont une longueur inférieure à 0.59 alors que cette valeur est à 0.54 pour les abalones adultes.
Ainsi, cette diférence de distribution se remmarque pour les variables suivantes : la longueur, le diamètre, la taille (beaucoup de valeurs extrêmes chez les adultes), la taille du trou, la taille des écailles, la taille des viscères et la taille du coquillage.

set.seed(1710)
ormeaux <- AbaloneAge #Création d'un dataset d'apprentissage et d'un dataset de validation
nb_lignes <- floor((nrow(ormeaux)*0.75)) #Nombre de lignes de l'échantillon d'apprentissage : %
ormeaux <- ormeaux[sample(nrow(ormeaux)), ] #Ajout de numéros de lignes
ormeaux.train <- ormeaux[1:nb_lignes, ] #Echantillon d'apprentissage
ormeaux.test <- ormeaux[(nb_lignes+1):nrow(ormeaux), ] #Echantillon de test

Un premier piège à éviter est donc d’évaluer la qualité de notre modèle final à l’aide des mêmes données qui ont servi pour l’entraînement. En effet, le modèle est complètement optimisé pour les données à l’aide desquelles il a été créé. L’erreur sera précisément minimum sur ces données. Alors que l’erreur sera toujours plus élevée sur des données que le modèle n’aura jamais vues !

Pour minimiser ce problème, la meilleure approche est de séparer dès le départ notre jeu de données en deux parties distinctes :

Le training set, qui va nous permettre d’entraîner notre modèle et sera utilisé par l’algorithme d’apprentissage.

Le testing set, qui permet de mesurer l’erreur du modèle final sur des données qu’il n’a jamais vues. On va simplement passer ces données comme s’il s’agissait de données que l’on n’a encore jamais rencontrées (comme cela va se passer ensuite en pratique pour prédire de nouvelles données) et mesurer la performance de notre modèle sur ces données.

C’est donc pour cela qu’on decide de faire deux échantillons ormeaux.train et ormeaux.test.

m1 <- rpart(
  formula = ormeaux.train$Age ~.,
  data    = ormeaux.train,
)

rpart.plot(m1,extra = 1)

Voici notre arbre, celui ci se lit de cette façon :

Premierement,ce fichier d’apprentissage contient 2642+142+366= 3132 ormeaux.

Ensuite, pour chaque feuille, R indique la classe prédite, le nombre d’individus de la classe prédite à gauche et le nombre d’individus des autres classes à droite. Par exemple, si on prend la première feuille à gauche qui correspond aux individus ayant un Shell weight <0.32, le chiffre à gauche correspond aux ormeaux adultes celui du milieu aux jeunes et celui de droite aux ormeaux vieux, la feuille contient 2582 ormeaux adultes, 40 ormeaux jeunes et 366 ormeaux vieux. Et cette sous-population est classée dans « Adulte ».

La variable la plus discriminante au sens de l’indice de Gini est Viscera weight(Poids des viscéres), soit ici pour une valeur de Viscera weight étant égale à 0.025.

A partir de cette variable discriminante le fichier peut être partagé en deux sous-fichiers qui séparent au mieux la variable-cible (Age). On trouve alors à gauche (Viscera Weight >0.025 VRAI) 2589 Adultes + 40 Jeunes + 366 Vieux, et à droite (Viscera Weight >0.025 FAUX), 35 Adultes, 102 Jeunes et 0 Vieux.

Au bout du compte, le fichier initial est “classé” en 7 classes (les 7 feuilles finales de l’arbre de décision). Les fréquences de Adultes/Jeunes/Vieux dans ces feuilles peuvent être utilisées comme des estimations de probabilités à priori par un module de prédiction.

Maintenant, nous devons valider notre fichier d’apprentissage.

Prediction <-predict(m1, ormeaux.test, type = 'class')
mc <- table(ormeaux.test$Age, Prediction)
mc
##          Prediction
##           adultes jeunes vieux
##   adultes     834     14    26
##   jeunes       14     33     0
##   vieux        88      0    36

On remarque que la matrice de confusion classent assez bien les différents ormeaux. Regardons ainsi l’erreur de classement.

erreur.classement<-1.0-(mc[1,1]+mc[2,2]+mc[3,3])/sum(mc)
erreur.classement
## [1] 0.1358852

L’erreur de classement est 13% ce qui est assez significatif. On peut en conclure que cette arbre de décision semble assez pertinent.

Réalisons maintenant une regression logistique sur cette variable Age afin de voir si cela coincide avec notre arbre de décision. Comme précisé au dessus nous sommes ici dans le cas de multicolinéarité. Ainsi, nous pouvons réaliser une régression logistique a l’aide de stepAIC afin de remedier à cela.

AbaloneAge$Age = as.factor(AbaloneAge$Age)
logit=glm(Age ~ ., family = "binomial",data=AbaloneAge)
modeleAIC=stepAIC(logit,direction = "both")
summary(modeleAIC)
## 
## Call:
## glm(formula = Age ~ Lenght + `Whole weight` + `Shucked weight` + 
##     `Viscera weight` + `Shell weight`, family = "binomial", data = AbaloneAge)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.3857  -0.5159  -0.3377  -0.1898   4.1790  
## 
## Coefficients:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         4.283      0.336  12.748  < 2e-16 ***
## Lenght            -19.317      1.085 -17.796  < 2e-16 ***
## `Whole weight`      9.983      1.038   9.615  < 2e-16 ***
## `Shucked weight`  -18.236      1.357 -13.434  < 2e-16 ***
## `Viscera weight`   -3.707      1.891  -1.961   0.0499 *  
## `Shell weight`     10.845      1.631   6.651 2.91e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 3708.2  on 4176  degrees of freedom
## Residual deviance: 2720.1  on 4171  degrees of freedom
## AIC: 2732.1
## 
## Number of Fisher Scoring iterations: 6

Selon le stepAIC, le modele ayant un AIC le plus faible est celui basé sur les variables suivantes : Age ~ Lenght + Whole weight + Shucked weight + Viscera weight + Shell weight. Les variables Length,Whole weight, Shucked Weight, et Shell Weight sont toutes très significatives, de même Viscera weight est aussi une variable qui explique bien l’Age. Ces variables sont les mêmes que celles présentes dans notre arbre de décision on en déduit donc que les résultats de l’arbre sont certainements bons et qu’il prédit bien la variable Age.

Régle de classification permettant de discriminer les males des femelles

Réalisons premièrement quelques graphiques afin de voir s’il existe un lien entre la variable Sex et les autres variables graphiquement.

grid.arrange(SL,SD,SH,SWW,SSW,SSA,SexShellW,SSR, ncol=3, nrow = 3)

On remarque à l’aide de ces graphiques qu’il n’existe pas a priori de réelles différences entre les femelles et les males au vu des différentes densités. Essayons donc de réaliser une régression logistique afin d’avoir une meilleure idée.

Enlevons premièrement la variable Infant de la variable Sex étant donné qu’elle est inutile a notre futur règle de classification.

AbaloneSex <- subset(Abalone, Abalone$Sex != "I")
AbaloneSex$Sex=as.factor(AbaloneSex$Sex)

Afin d’établir une règle de classification nous pouvons réaliser une ACP, de plus celles-ci est utile lors d’un cas de multicolinéarité comme-celui ci.

library(FactoMineR)
res.PCA<-PCA(AbaloneSex,quali.sup=c(1),graph=FALSE)
res.PCA$eig
##         eigenvalue percentage of variance cumulative percentage of variance
## comp 1 6.211521177             77.6440147                          77.64401
## comp 2 0.940240249             11.7530031                          89.39702
## comp 3 0.413594420              5.1699303                          94.56695
## comp 4 0.184595328              2.3074416                          96.87439
## comp 5 0.126046742              1.5755843                          98.44997
## comp 6 0.093853119              1.1731640                          99.62314
## comp 7 0.020843840              0.2605480                          99.88369
## comp 8 0.009305125              0.1163141                         100.00000

Les 2 premiers axes de l’ analyse expriment 89.4% de l’inertie totale du jeu de données ; cela signifie que 89.4% de la variabilité totale du nuage des individus (ou des variables) est représentée dans ce plan. C’est un pourcentage élevé, et le premier plan représente donc bien la variabilité contenue dans une très large part du jeu de données actif.

Le premier axe est prépondérant : il explique a lui seul 77.64% de la variabilité totale des données. Il convient de noter que dans un tel cas, la variabilité liée aux autres composantes peut être dénuée de sens, en dépit d’un pourcentage élevé.

Une estimation du nombre pertinent d’axes à interpréter suggère de restreindre l’analyse à la description du premier axe. Cette observation suggère que seul cet axe est porteur d’une véritable information. En conséquence, la description de l’analyse sera restreinte à ce seul axe.

La dimension 1 oppose des individus caractérisés par une coordonnée fortement positive sur l’axe (à droite du graphe) à des individus caractérisés par une coordonnée fortement négative sur l’axe (à gauche du graphe).

Le groupe 1 (caractérisés par une coordonnée positive sur l’axe) partage :

de fortes valeurs pour les variables Rings, Height, Shell.weight, Diameter, Lenght et Whole.weight (de la plus extrême à la moins extrême).

Le groupe 2 (caractérisés par une coordonnée positive sur l’axe) partage :

de fortes valeurs pour les variables Lenght, Diameter, Shucked.weight, Viscera.weight, Whole.weight, Height et Shell.weight (de la plus extrême à la moins extrême). de faibles valeurs pour la variable Rings.

Le groupe 3 (caractérisés par une coordonnée positive sur l’axe) partage :

de fortes valeurs pour les variables Whole.weight, Viscera.weight, Shucked.weight, Shell.weight, Lenght, Diameter, Height et Rings (de la plus extrême à la moins extrême).

Le groupe 4 (caractérisés par une coordonnées négative sur l’axe) partage :

de faibles valeurs pour les variables Diameter, Lenght, Height, Whole.weight, Shell.weight, Viscera.weight, Shucked.weight et Rings (de la plus extrême à la moins extrême).

Notons que les variables Lenght, Diameter et Whole weight sont extrêmement corrélées à cette dimension (corrélations respectives de 0.92, 0.92, 0.96). Ces variables pourraient donc résumer à elles seules la dimension 1.

Cependant, nous remarquons qu’il existe de la multicolinéarité entre les variables. Les résultats définis préalablement sont donc à nuancer. En effet, les deux premières composantes principales représentent, à elles seules, environ 89,4% de l’information principale. De plus, en dehors de la première valeur propre, les composantes principales ont des valeurs en dessous de 1.
Il y a bien de la multicolinéarité dans les variables.

Une première solution à la multicolinéarité très présente pour cette règle de classification contrairement à tout à l’heure, est de choisir un modèle avec le moins de variables possibles et un critère AIC faible.

summary(m_aic)
## 
## Call:
## glm(formula = Sex ~ Diameter + Height + `Shucked weight` + `Viscera weight`, 
##     family = "binomial", data = train)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.6630  -1.1923   0.8942   1.1323   1.5322  
## 
## Coefficients:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        2.2173     0.4009   5.531 3.19e-08 ***
## Diameter          -4.7727     1.3658  -3.494 0.000475 ***
## Height            -3.4786     2.1897  -1.589 0.112141    
## `Shucked weight`   2.6947     0.5144   5.239 1.62e-07 ***
## `Viscera weight`  -2.7651     1.1070  -2.498 0.012495 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 3137.3  on 2267  degrees of freedom
## Residual deviance: 3082.9  on 2263  degrees of freedom
## AIC: 3092.9
## 
## Number of Fisher Scoring iterations: 4

Ce sont les variables du diamètre, de la taille, de la taille des écailles et de la taille des viscères qui sont sélectionnées par le stepAIC. Toutes les variables sont significatives à 5% sauf la variable taille.
De plus, nous avons calculé le Root Mean Squares Error of Prediction (erreur quadratique moyenne). Nous obtenons pour ce modèle un RMSEP de 1.16 alors que pour celui du modèle m0 le RMSEP est de 1.17. Une différence aussi minimime reste cependant à nuancer. En effet, comme vu plus haut lors de l’analyse graphique et de l’analyse de l’ACP on remarque qu’il est très difficile de différencier les deux sexes. Une idée serait encore d’utiliser un arbre de décision afin de vérifier cela.

Arbre de décision de la variable Sexe

Nous allons maintenant construire l’arbre de decision

set.seed(1717)
m2 <- rpart(
  formula = ormeauxsex.train$Sex ~ .  ,
  data    = ormeauxsex.train,
)

rpart.plot(m2,extra = 1)

#Matrice de confusion

Prediction <-predict(m2, ormeauxsex.test, type = 'class')
mc <- table(ormeauxsex.test$Sex, Prediction)
mc
##    Prediction
##       F   I   M
##   F  93   0 222
##   I   0   0   0
##   M  83   0 311
erreur.classement<-1.0-(mc[1,1]+mc[2,2])/sum(mc)
erreur.classement
## [1] 0.8688293

On observe donc une erreur de classement étant égale à 86%. On peut en conclure grâce aux différentes analyses antérieurs qu’il n’existe pas de réelles variables physiologiques pouvant classer d’un côté les mâles et d’un autre les femelles.