1 Introduction

1.1 Installation

Renvoyer vers une vidéo

1.2 Bases du language R

Renvoyer vers un cours existant

1.3 Importer des données

#Depuis un fichier .txt
data <- read.table(file = "path/file.txt", header = TRUE / FALSE, sep = .. , dec = .. )

#Depuis un fichier .csv
data <- read.csv(file = "path/file.csv", header = TRUE / FALSE, sep = .. , dec = .. )

#Depuis un fichier .xlsx
data <- readxl::read_xlsx(path = "path/file.xlsx", sheet = .. , col_names = TRUE / FALSE)

Rmq

  • Si les observations sont séparées par TAB, utiliser sep="\t"
  • Si il y a une observation par case, utiliser sep=";"

1.4 Dataframe

1.4.1 Définition

Les données importées dans R se retrouvent dans un objet particulier appelé Dataframe. C’est comme son nom l’indique un tableau de données dans lequel les lignes correspondent aux individus et les colonnes aux variables.

1.4.2 Préparation

Supposons que nous disposions d’un dataframe df. Avant de commencer une quelconque analyse il faut le préparer.

Noms des variables
S’ils ne sont pas valides (ex: un nom de variable en 2 mots) ou bien simplement s’ils ne sont pas très clairs il est possible de renommer soi-même les variables.

colnames(df) <- c("Var1", "Var2", .., "VarM") # M variables dans df

Encodage des variables qualitatives
Celles-ci doivent être encodées commes des facteurs à plusieurs niveaux. Un niveau correspond à une modalité. Par défaut les variables qualitatives sont encodées comme numériques (si les modalités sont représentées par des chiffres) ou bien comme chaine de caractère.

data$VarQuali <- as.factor(data$VarQuali) # pour une variable qualitative nommée VarQuali

Données manquantes
La façon la plus simple de gérer cela est d’omettre chaque ligne de df où des données sont manquantes.

df <- na.omit(df)

1.4.3 Manipulation

Pour manipuler les données contenues dans un dataframe, la package dplyr est fortement recommandé car les fonctions sont très simples.

#Filtrer les individus
dplyr::filter(df, ..)

#Sélectioner des variables
dplyr::select(df, ..)

#Ajouter une nouvelle variable
dplyr::mutate(df, ..)

#Transformer une variable existante
dplyr::transmute(df, ..)

1.4.4 Création

Une autre façon que l’importation de données extérieures pour créer un dataframe est tout simplement de le créer à la main.

#Exemple avec un dataframe df composé de 3 variables nommées Var1, Var2, Var3
df <- data.frame(Var1 = ..,
                 Var2 = ..,  # A gauche de l'égalité le nom de chaque variable et à droite un vecteur.    
                 Var3 = ..)  # /!\ Les vecteurs à droite doivent tous être de même taille.

1.5 Rmarkdown

Renvoyer vers un cours existant

1.6 Packages

Installer un package

install.packages("NomPackage") #Cela ne se fait qu'une seule fois sur votre machine.

Charger un package pour utiliser ses fonctions

library(NomPackage) #Cela se fait à chaque nouvelle session de R

Appeler une fonction
Supposons que nous voulions utiliser la fonction filter() du package dplyr.

#1) Le package est installé et on a chargé en début de session le package dplyr avec "library()"
filter(..)

#2) Le package est juste installé mais pas chargé
dplyr::filter(..)

Liste des packages utiles dans ce cours

PACKAGE
readxl Importer des fichiers Excel
rmarkdown Créer des rapports en HTML depuis Rstudio
pander Créer des tableaux plus jolis dans vos rapports HTML
dplyr Gérer plus facilement les dataframes
car Tests statistiques et anova non balancée
Hmisc Test de corrélation 2 à 2
emmeans Test d’hypothèses et contrasts
Visreg Visualiser une régression
FactoMineR Analyse en composantes principales
factoextra Clustering
ggplot2 1 Graphiques plus beaux et plus complexes

1.7 Références internet utiles

Statistiques descriptives

  • …..
  • …..
  • …..

Inférence statistique

  • …..
  • …..
  • …..

Modèles linéaires

  • [Régression linéaire]
  • [ANOVA]

Modèles linéaires généralisés

  • [Régression logistique]
  • [Régression de Poisson]

Analyses mutivariées

2 Statistique descrptive

Que faut-il faire pour décrire numériquement et / ou visualiser tel type de variable?

2.1 Numérique

Tendance centrale

  • Moyenne ~> mean()
  • Médiane ~> median()

Variabilité

  • Minimum / maximum ~> min() / max()
  • Range ~> range()
  • Variance ~> var()
  • Ecart-type ~> sd()

Quantiles

  • Quantile ~> quantile()
quantile(df$varQuanti, p=0.25) #Quantile d'ordre 0.25
  • Ecarte interquartile ~> IQR()

Autres

  • Résumé statistique ~> summary()
summary(df$varQuanti)
  • Table de fréquence ou contingence ~> table()
  • Table de proportion ~> prop.table()
table(df$varQuali) #Fréquence
prop.table( table(df$varQuali) ) #Proportion
table(df$varQuali_1, df$varQuali_2) #Contingence
  • Coefficient / matrice de corrélation ~> cor()
df_quanti <- dplyr::select(df, varQuanti, varQuanti1, varQuanti2)

cor(df$varQuanti_1, df$varQuanti_2) #Coefficient
cor(df_quanti) #Matrice
  • Informations par modalités ~> tapply()
tapply(df$varQuanti, df$varQuali, FUN = mean) #ex: moyenne par modalités

2.2 Graphique

  • Scatter plot ~> plot()
plot(varQuanti1 ~ varQuanti2, data = df, main = "Scatter plot")

  • Matrice de scatter plots ~> pairs()
df_quanti <- dplyr::select(df, varQuanti, varQuanti1, varQuanti2)
pairs(df_quanti)

  • Histogramme ~> hist()
hist(df$varQuanti, freq = TRUE , main = "Histogramme de fréquence d'une variable discrète", xlab = "Var quanti discrète")

  • Bar plot ~> barplot()
barplot(table(df$varQuali2), main = "Bar plot d'une variable qualitative")

barplot(df$varQuanti1, main = "Bar plot d'une variable quantitative", xlab = "Var quantitative")

barplot(table(df$varQuali1,df$varQuali2), beside = TRUE, legend = levels(df$varQuali1), main ="Barplot par modalités")

  • Box plot ~> boxplot()
boxplot(df$varQuanti, main= "Boxplot d'une variable quantitative", ylab = "Var quanti")

boxplot(varQuanti ~ varQuali, data = df, main= "Boxplot d'une variable quantitative par modalités",  ylab = "Var quanti")

boxplot(varQuanti ~ varQuali1*varQuali2, data = df, main= "Boxplot d'une variable quantitative par modalités croisées",  ylab = "Var quanti")

  • Mosaic plot ~> mosaicplot()
mosaicplot(varQuali1 ~ varQuali2, data = df, main = "Mosaic plot de deux variables qualitatives")

3 Inference statistique

Données utilisées pour illustrer certaines fonctions: chocolat.xlsx

3.1 Tests sur les variances

#Comparaison de 2 variances - H0: var1 = var2 / H1: var1 != var2
car::leveneTest(choco$Cacao, group = choco$Marque)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)  
## group  1  6.1977 0.0228 *
##       18                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

3.2 Tests sur les moyennes

#Test sur 1 moyenne - H0: mu = 5 / H1: mu != 5
t.test(choco$Cacao, mu = 5, alternative = "two.sided")
## 
##  One Sample t-test
## 
## data:  choco$Cacao
## t = 24, df = 19, p-value = 1.131e-15
## alternative hypothesis: true mean is not equal to 5
## 95 percent confidence interval:
##  6.059293 6.261707
## sample estimates:
## mean of x 
##    6.1605
##Comparaison de 2 moyennes / var égale ou non - H0: mu1 = mu2 / H1: mu1 != mu2
t.test(x = choco$Cacao[choco$Marque=="Cabeau"], y = choco$Cacao[choco$Marque=="Gallet"], alternative = "two.sided", var.equal = FALSE)
## 
##  Welch Two Sample t-test
## 
## data:  choco$Cacao[choco$Marque == "Cabeau"] and choco$Cacao[choco$Marque == "Gallet"]
## t = -2.2137, df = 12.087, p-value = 0.04682
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.386775949 -0.003224051
## sample estimates:
## mean of x mean of y 
##     6.063     6.258
##Comparaison de 2 moyennes / obs pairées
t.test(x = .., y = .., alternative = "two.sided", paired = TRUE)

3.3 Tests sur les proportions

#Test sur 1 proportion - H0: p = 0.25 / H1: p!= 0.25 (Succès = etre de la marque Cabeau)
succes <- sum(choco$Marque=="Cabeau")
total <- length(choco$Marque)
prop.test(x = succes, n = total, p= 0.25, alternative = "two.sided")
## 
##  1-sample proportions test with continuity correction
## 
## data:  succes out of total, null probability 0.25
## X-squared = 5.4, df = 1, p-value = 0.02014
## alternative hypothesis: true p is not equal to 0.25
## 95 percent confidence interval:
##  0.2785367 0.7214633
## sample estimates:
##   p 
## 0.5
#Test sur plusieurs proportions - H0: p1 = 0.25, p2 = 0.75 / H1: p1 != 0.25 ou p2 !=0.75
chisq.test(x = table(choco$Marque), p = c(0.25, 0.75))
## 
##  Chi-squared test for given probabilities
## 
## data:  table(choco$Marque)
## X-squared = 6.6667, df = 1, p-value = 0.009823
#Comparer deux proportions - H0: p1 = p2 / H1: p1 != p2
M <- matrix( c(NbSuccesGrp1, NbSuccesGrp2, NbEchecGrp1, NbEchecGrp2), ncol=2)
fisher.test(M, alternative = "two.sided")

3.4 Tests de corrélation / indépendance

#Test de corrélation
cor.test(x = .., y = .., alternative = "two.sided") #x,y quantitatives
#Test d'indépendance
chisq.test(x =.., y = ..) 
#Tests de corrélation 2 à 2
res <- Hmisc::rcorr(df_quanti)
                # df_quanti = df sans les variables qualitatives

~> res$r = matrice de corrélation

~> res$p = matrice des p-valeurs associées

3.5 QQ-plot

#Tester la normalité d'une distribution
qqnorm(y = choco$Cacao)
qqline(y = choco$Cacao, col = "red")

4 Modeles lineaires

Objectif

Expliquer une variable réponse quantitative en fonction de variables explicatives quantitatives et / ou qualitatives + termes d’interaction.

4.1 Modelisation

#Régression linéaire simple - X quantitative
mod <- lm( Y ~ X, data = ..) 
#Régression linéaire multiple - X1,X2 quantitatives
mod <- lm( Y ~ X1 + X2, data = ..) 
#ANOVA 1 - X qualitative
mod <- lm( Y ~ -1 + X, data = ..)
#ANOVA 2 - X1,X2 qualitatives
mod <- lm( Y ~ -1 + X1 * X2, data = ..)
#Modèle linéaire général - X1 quantitative, X2 qualitative
mod <- lm( Y ~ X1 * X2, data = ..)

Formules dans R

Le premier argument de la fonction lm() est un objet de type formule.

  • Y ~ X = modéliser Y en fonction de X avec intercept.
  • Y ~ -1 + X = modéliser Y en fonction de X sans intercept.
  • Y ~ X1 + X2 = modéliser Y en fonction de X1 et X2 avec intercept.
  • Y ~ X1 * X2 = modéliser Y en fonction de X1, X2 et de l’interaction X1:X2 avec intercept = Y ~ X1 + X2 + X1:X2.
  • Y ~ f(X) = modéliser Y en fonction de f(X), une fonction quelconque de X, avec intercept.

Tout cela peut bien entendu être mixé pour créer un modèle bien spécifique.
Par exemple pour une régression quadratique sans intercept on utilisera Y ~ -1 + X + X^2.

data = ..

Pour le second argument, il est fortement conseillé de travailler avec un dataframe dans lequel se trouve toutes nos variables d’intéret. Supposons que les variables de l’objet formule se trouvent dans le dataframe df alors on notera data = df.

4.2 Fonctions utiles sur un objet de type mod <- lm(..)

Modèle summary(mod)
Signification de coefficients
+ Mise à jour du modèle		 drop1(mod, ..)
+ update(mod, ..)
Coefficients
+ IC		 coef(mod)
+ confint(mod)
Résidus / Valeurs ajustées residuals(mod) / fitted(mod)
Matrice Variance-Covariance vcov(mod)
Analyse de la variance anova(mod) TYPE I
car::Anova(mod, ..) TYPE II ou III
Graphes des diagnostics plot(mod)
Visualiser la modélisation visreg::visreg(mod, ..)
Prédictions predict(mod, ..)
Inférence car::leveneTest(mod)
emmeans::lsmeans(mod,..)
emmeans::contrast(..)

4.3 Exemple 1: regression lineaire simple

#Import et préparation des données
calibration <- readxl::read_xlsx(path = "data/calibration2.xlsx", col_names = TRUE)
#Visualisation des données
plot(SURFACE ~ ETHANOL, data = calibration, main = "Surface sous le chromatogramme \n en fonction de la concentration en éthanol")

#Modélisation
mod <- lm(SURFACE ~ ETHANOL, data = calibration)
#Modèle
summary(mod)
## 
## Call:
## lm(formula = SURFACE ~ ETHANOL, data = calibration)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.0679 -1.2971 -0.8394  0.1453  4.4014 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -1.3422     1.2603  -1.065    0.318    
## ETHANOL       8.2347     0.8858   9.296 1.46e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.307 on 8 degrees of freedom
## Multiple R-squared:  0.9153, Adjusted R-squared:  0.9047 
## F-statistic: 86.41 on 1 and 8 DF,  p-value: 1.46e-05
#Vérifier les hypothèses du modèle linéaire
par(mfrow=c(1,2))
plot(mod)

#Visualiser la régression
visreg::visreg(mod)

#Prédire avec intervalle de prédiction pour de nouvelles données ETHANOL = 0.25 / 2.8 / 5 
new <- data.frame(ETHANOL = c(0.25,2.8,5))
pred <- predict(mod, new, interval = "prediction")
data.frame(new, Fit = pred[,1], Lwr = pred[,2], Upr = pred[,3])
##   ETHANOL        Fit       Lwr       Upr
## 1    0.25  0.7164505 -5.165253  6.598154
## 2    2.80 21.7148585 15.206244 28.223473
## 3    5.00 39.8311321 30.204653 49.457611

4.4 Exemple 2: regression lineaire multiple

#Import et préparation des données
cubage <- readxl::read_xlsx(path = "data/Volume_arbres.xlsx", col_names = TRUE)
#Visualisation des données
pairs(lnVolume ~ lnDiam + lnHauteur, data = cubage, main = "Matrice des scatter plots \n entre Log(Volume), Log(Diametre) et Log(Hauteur)")

#Modélisation
mod <- lm(lnVolume ~ lnDiam + lnHauteur, data = cubage)
#Modèle
summary(mod)
## 
## Call:
## lm(formula = lnVolume ~ lnDiam + lnHauteur, data = cubage)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.168561 -0.048488  0.002431  0.063637  0.129223 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -6.63162    0.79979  -8.292 5.06e-09 ***
## lnDiam       1.98265    0.07501  26.432  < 2e-16 ***
## lnHauteur    1.11712    0.20444   5.464 7.81e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.08139 on 28 degrees of freedom
## Multiple R-squared:  0.9777, Adjusted R-squared:  0.9761 
## F-statistic: 613.2 on 2 and 28 DF,  p-value: < 2.2e-16
#Vérifier les hypothèses du modèle linéaire
par(mfrow=c(1,2))
plot(mod)

#Visualiser la régression
visreg::visreg(mod)

#Prédire avec intervalle de prédiction pour de nouvelles données (lnDiam, lnHauteur) = (2.1, 4.28) / (2.8, 4.15) / (3, 4.4)
new <- data.frame(lnDiam = c(2.1,2.8,3), lnHauteur = c(4.28, 4.15, 4.4))
pred <- predict(mod, new, interval = "prediction")
data.frame(new, Fit = pred[,1], Lwr = pred[,2], Upr = pred[,3])
##   lnDiam lnHauteur      Fit      Lwr      Upr
## 1    2.1      4.28 2.313236 2.132834 2.493637
## 2    2.8      4.15 3.555864 3.359593 3.752136
## 3    3.0      4.40 4.231675 4.052660 4.410690

4.5 Exemple 3: ANOVA 1

#Import et préparation des données
diet <- readxl::read_xlsx(path = "data/diet.xlsx", col_names = TRUE)
colnames(diet) <- c("Produit", "Dif_poids", "Produit_T")
#AVANT TOUTE CHOSE: vérifier l'encodage des variables qualitatives utilisées
diet$Produit_T <- as.factor(diet$Produit_T)
#Visualisation des données
boxplot(Dif_poids ~ Produit_T, data = diet, xlab = "Produit_T", ylab= "Dif_poids", 
        main = "Différence de poids \n en fonction du produit utilisé")

#Modéliser
mod <- lm(Dif_poids ~ Produit_T, data = diet)
#Modèle
summary(mod)
## 
## Call:
## lm(formula = Dif_poids ~ Produit_T, data = diet)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.7833 -3.9521  0.3083  3.4229  9.6167 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -7.117      1.449  -4.912 1.29e-05 ***
## Produit_TB     9.058      2.049   4.421 6.34e-05 ***
## Produit_TC     3.200      2.049   1.562 0.125470    
## Produit_TP     8.258      2.049   4.031 0.000217 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.018 on 44 degrees of freedom
## Multiple R-squared:  0.3744, Adjusted R-squared:  0.3317 
## F-statistic: 8.777 on 3 and 44 DF,  p-value: 0.0001126
#Analyse de la variance
anova(mod)
## Analysis of Variance Table
## 
## Response: Dif_poids
##           Df  Sum Sq Mean Sq F value    Pr(>F)    
## Produit_T  3  663.12 221.040  8.7769 0.0001126 ***
## Residuals 44 1108.11  25.184                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Visualiser l'ANOVA
visreg::visreg(mod)

#Vérifier les hypothèses du modèle linéaire
par(mfrow=c(1,2))
plot(mod)

#Test d'homogénéité des variances
car::leveneTest(mod)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  3  1.1866 0.3258
##       44
#Estimations et comparaisons 2 à 2 des moyennes - Produit_T
lsm <- emmeans::lsmeans(mod, pairwise ~ Produit_T)
lsm
## $lsmeans
##  Produit_T lsmean   SE df lower.CL upper.CL
##  A          -7.12 1.45 44  -10.036   -4.197
##  B           1.94 1.45 44   -0.978    4.861
##  C          -3.92 1.45 44   -6.836   -0.997
##  P           1.14 1.45 44   -1.778    4.061
## 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast estimate   SE df t.ratio p.value
##  A - B       -9.06 2.05 44 -4.421  0.0004 
##  A - C       -3.20 2.05 44 -1.562  0.4105 
##  A - P       -8.26 2.05 44 -4.031  0.0012 
##  B - C        5.86 2.05 44  2.859  0.0316 
##  B - P        0.80 2.05 44  0.390  0.9795 
##  C - P       -5.06 2.05 44 -2.469  0.0791 
## 
## P value adjustment: tukey method for comparing a family of 4 estimates
plot(lsm)

#Contrasts
hyp <- list("Hyp1: A = B" = c(1,-1,0,0),                    # Ici vecteur des coefficients de l'hypothèse
            "Hyp2: Traitement = Controle" = c(1,1,1,-3))    # de l'hypothèse correspond à l'ordre  (A, B, C, P)
emmeans::contrast(lsm$lsmeans, hyp)
##  contrast                    estimate   SE df t.ratio p.value
##  Hyp1: A = B                    -9.06 2.05 44 -4.421  0.0001 
##  Hyp2: Traitement = Controle   -12.52 5.02 44 -2.494  0.0165

4.6 Exemple 4: ANOVA 2

#Import et préparation des données
mais <- readxl::read_xlsx(path = "data/mais_anova2.xlsx", col_names = TRUE)
#AVANT TOUTE CHOSE: vérifier l'encodage des variables qualitatives utilisées
mais$Graines <- as.factor(mais$Graines)
mais$Engrais <- as.factor(mais$Engrais)
#Visualisation des données
boxplot(Production ~ Graines*Engrais, data = mais, xlab = "Modalités croisées de Graines et Engrais", ylab= "Production", 
        main = "Production de mais \n en fonction du type de graine et d'engrais")

#Modéliser
mod <- lm(Production ~ Graines*Engrais, data = mais)
#Intercation significative? 
drop1(mod, test="Chisq") 
## Single term deletions
## 
## Model:
## Production ~ Graines * Engrais
##                 Df Sum of Sq    RSS    AIC  Pr(>Chi)    
## <none>                        38.50 31.685              
## Graines:Engrais  4    243.67 282.17 59.538 3.102e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Si oui, on passe à la suite
#Si non, 'mod <- update(mod, . ~ . - Graines::Engrais)'
#Modèle
summary(mod)
## 
## Call:
## lm(formula = Production ~ Graines * Engrais, data = mais)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##   -2.5   -1.0    0.0    1.0    2.5 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        112.000      1.462  76.582 5.59e-14 ***
## GrainesB           -15.000      2.068  -7.252 4.80e-05 ***
## GrainesC           -16.500      2.068  -7.978 2.26e-05 ***
## Engrais2           -14.500      2.068  -7.011 6.25e-05 ***
## Engrais3            -8.500      2.068  -4.110  0.00264 ** 
## GrainesB:Engrais2    2.000      2.925   0.684  0.51134    
## GrainesC:Engrais2    5.500      2.925   1.880  0.09275 .  
## GrainesB:Engrais3    7.500      2.925   2.564  0.03048 *  
## GrainesC:Engrais3   21.000      2.925   7.180 5.20e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.068 on 9 degrees of freedom
## Multiple R-squared:  0.9713, Adjusted R-squared:  0.9458 
## F-statistic:  38.1 on 8 and 9 DF,  p-value: 4.767e-06
#Données balancées? 
table(mais$Graines, mais$Engrais)
##    
##     1 2 3
##   A 2 2 2
##   B 2 2 2
##   C 2 2 2
#Analyse de la variance

#Si balancées, on utilise 'anova(mod)' pour analyser la variance
#Si non balancées, on utilise 'car::Anova(mod, type=3)' pour analyser les variance

anova(mod)
## Analysis of Variance Table
## 
## Response: Production
##                 Df Sum Sq Mean Sq F value    Pr(>F)    
## Graines          2 432.33 216.167  50.532 1.278e-05 ***
## Engrais          2 628.00 314.000  73.403 2.676e-06 ***
## Graines:Engrais  4 243.67  60.917  14.240 0.0006255 ***
## Residuals        9  38.50   4.278                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Visualiser l'ANOVA
visreg::visreg(mod, "Graines", by = "Engrais", main = "Visualisation de l'ANOVA de Production ~ Graines \n par modalités d'Engrais")

#Vérifier les hypothèses du modèle linéaire
par(mfrow=c(1,2))
plot(mod)

#Test d'homogénéité des variances
car::leveneTest(mod)
## Warning in anova.lm(lm(resp ~ group)): ANOVA F-tests on an essentially
## perfect fit are unreliable
## Levene's Test for Homogeneity of Variance (center = median)
##       Df    F value    Pr(>F)    
## group  8 7.5689e+30 < 2.2e-16 ***
##        9                         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Estimations et comparaisons 2 à 2 des moyennes - Graines
lsmG <- emmeans::lsmeans(mod, pairwise ~ Graines)
## NOTE: Results may be misleading due to involvement in interactions
lsmG
## $lsmeans
##  Graines lsmean    SE df lower.CL upper.CL
##  A        104.3 0.844  9    102.4    106.2
##  B         92.5 0.844  9     90.6     94.4
##  C         96.7 0.844  9     94.8     98.6
## 
## Results are averaged over the levels of: Engrais 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast estimate   SE df t.ratio p.value
##  A - B       11.83 1.19  9  9.910  <.0001 
##  A - C        7.67 1.19  9  6.420  0.0003 
##  B - C       -4.17 1.19  9 -3.489  0.0170 
## 
## Results are averaged over the levels of: Engrais 
## P value adjustment: tukey method for comparing a family of 3 estimates
plot(lsmG)

#Estimations et comparaisons 2 à 2 des moyennes - Engrais
lsmE <- emmeans::lsmeans(mod, pairwise ~ Engrais)
## NOTE: Results may be misleading due to involvement in interactions
lsmE
## $lsmeans
##  Engrais lsmean    SE df lower.CL upper.CL
##  1        101.5 0.844  9     99.6    103.4
##  2         89.5 0.844  9     87.6     91.4
##  3        102.5 0.844  9    100.6    104.4
## 
## Results are averaged over the levels of: Graines 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast estimate   SE df t.ratio p.value
##  1 - 2          12 1.19  9  10.049 <.0001 
##  1 - 3          -1 1.19  9  -0.837 0.6905 
##  2 - 3         -13 1.19  9 -10.887 <.0001 
## 
## Results are averaged over the levels of: Graines 
## P value adjustment: tukey method for comparing a family of 3 estimates
plot(lsmE)

#Estimations et comparaisons 2 à 2 des moyennes - modalités croisées
lsm <- emmeans::lsmeans(mod, pairwise ~ Graines*Engrais)
lsm
## $lsmeans
##  Graines Engrais lsmean   SE df lower.CL upper.CL
##  A       1        112.0 1.46  9    108.7    115.3
##  B       1         97.0 1.46  9     93.7    100.3
##  C       1         95.5 1.46  9     92.2     98.8
##  A       2         97.5 1.46  9     94.2    100.8
##  B       2         84.5 1.46  9     81.2     87.8
##  C       2         86.5 1.46  9     83.2     89.8
##  A       3        103.5 1.46  9    100.2    106.8
##  B       3         96.0 1.46  9     92.7     99.3
##  C       3        108.0 1.46  9    104.7    111.3
## 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast  estimate   SE df t.ratio p.value
##  A 1 - B 1     15.0 2.07  9   7.252 0.0009 
##  A 1 - C 1     16.5 2.07  9   7.978 0.0004 
##  A 1 - A 2     14.5 2.07  9   7.011 0.0012 
##  A 1 - B 2     27.5 2.07  9  13.296 <.0001 
##  A 1 - C 2     25.5 2.07  9  12.329 <.0001 
##  A 1 - A 3      8.5 2.07  9   4.110 0.0406 
##  A 1 - B 3     16.0 2.07  9   7.736 0.0006 
##  A 1 - C 3      4.0 2.07  9   1.934 0.6113 
##  B 1 - C 1      1.5 2.07  9   0.725 0.9967 
##  B 1 - A 2     -0.5 2.07  9  -0.242 1.0000 
##  B 1 - B 2     12.5 2.07  9   6.044 0.0035 
##  B 1 - C 2     10.5 2.07  9   5.077 0.0114 
##  B 1 - A 3     -6.5 2.07  9  -3.143 0.1502 
##  B 1 - B 3      1.0 2.07  9   0.483 0.9998 
##  B 1 - C 3    -11.0 2.07  9  -5.318 0.0084 
##  C 1 - A 2     -2.0 2.07  9  -0.967 0.9804 
##  C 1 - B 2     11.0 2.07  9   5.318 0.0084 
##  C 1 - C 2      9.0 2.07  9   4.351 0.0294 
##  C 1 - A 3     -8.0 2.07  9  -3.868 0.0563 
##  C 1 - B 3     -0.5 2.07  9  -0.242 1.0000 
##  C 1 - C 3    -12.5 2.07  9  -6.044 0.0035 
##  A 2 - B 2     13.0 2.07  9   6.285 0.0026 
##  A 2 - C 2     11.0 2.07  9   5.318 0.0084 
##  A 2 - A 3     -6.0 2.07  9  -2.901 0.2061 
##  A 2 - B 3      1.5 2.07  9   0.725 0.9967 
##  A 2 - C 3    -10.5 2.07  9  -5.077 0.0114 
##  B 2 - C 2     -2.0 2.07  9  -0.967 0.9804 
##  B 2 - A 3    -19.0 2.07  9  -9.186 0.0001 
##  B 2 - B 3    -11.5 2.07  9  -5.560 0.0062 
##  B 2 - C 3    -23.5 2.07  9 -11.362 <.0001 
##  C 2 - A 3    -17.0 2.07  9  -8.219 0.0003 
##  C 2 - B 3     -9.5 2.07  9  -4.593 0.0213 
##  C 2 - C 3    -21.5 2.07  9 -10.395 0.0001 
##  A 3 - B 3      7.5 2.07  9   3.626 0.0783 
##  A 3 - C 3     -4.5 2.07  9  -2.176 0.4853 
##  B 3 - C 3    -12.0 2.07  9  -5.802 0.0047 
## 
## P value adjustment: tukey method for comparing a family of 9 estimates
plot(lsm)

4.7 Exemple 5: modele lineaire general

#Import et préparation des données
digestion <- readxl::read_xlsx(path = "data/DigestionEnz09.xlsx", col_names = TRUE)
#AVANT TOUTE CHOSE: vérifier l'encodage des variables qualitatives utilisées
digestion$Standard <- as.factor(digestion$Standard)
#Visualisation des données
plot(DigestionLOG ~ PH, data = digestion, col = levels(Standard), main = "Log(Digestion) en fonction du PH \n par modalités de Standard")
legend("topleft", legend = levels(digestion$Standard), pch = c(1,1,1), col = levels(digestion$Standard), cex=0.5)

#Modéliser
mod <- lm(DigestionLOG ~ PH * Standard, data = digestion)
#Intercation significative? 
drop1(mod, test = "Chisq")
## Single term deletions
## 
## Model:
## DigestionLOG ~ PH * Standard
##             Df Sum of Sq       RSS     AIC  Pr(>Chi)    
## <none>                   0.0076721 -127.69              
## PH:Standard  2  0.023196 0.0308682 -106.63 3.619e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Si oui, on passe à la suite
#Si non, 'mod <- update(mod, . ~ . - Graines::Engrais)'
#Modèle
summary(mod)
## 
## Call:
## lm(formula = DigestionLOG ~ PH * Standard, data = digestion)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.034167 -0.017938  0.001167  0.019708  0.031250 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   0.20767    0.07655   2.713 0.018862 *  
## PH            0.05050    0.01264   3.994 0.001779 ** 
## Standard2     0.34983    0.10827   3.231 0.007202 ** 
## Standard3    -0.43667    0.10827  -4.033 0.001659 ** 
## PH:Standard2  0.00200    0.01788   0.112 0.912783    
## PH:Standard3  0.09425    0.01788   5.271 0.000197 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.02529 on 12 degrees of freedom
## Multiple R-squared:  0.9851, Adjusted R-squared:  0.9789 
## F-statistic: 159.1 on 5 and 12 DF,  p-value: 1.548e-10
#Analyse de la variance
anova(mod)
## Analysis of Variance Table
## 
## Response: DigestionLOG
##             Df  Sum Sq  Mean Sq F value    Pr(>F)    
## PH           1 0.08184 0.081840 128.007 9.292e-08 ***
## Standard     2 0.40362 0.201810 315.654 4.213e-11 ***
## PH:Standard  2 0.02320 0.011598  18.141 0.0002357 ***
## Residuals   12 0.00767 0.000639                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Visualiser le modèle
visreg::visreg(mod, "PH", by = "Standard")

#Vérifier les hypothèses du modèle linéaire
par(mfrow=c(1,2))
plot(mod)

#Estimations et comparaisons 2 à 2 des moyennes - Standard
lsmS <- emmeans::lsmeans(mod, pairwise ~ Standard)
## NOTE: Results may be misleading due to involvement in interactions
lsmS
## $lsmeans
##  Standard lsmean     SE df lower.CL upper.CL
##  1         0.511 0.0103 12    0.488    0.533
##  2         0.873 0.0103 12    0.850    0.895
##  3         0.639 0.0103 12    0.617    0.662
## 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast estimate     SE df t.ratio p.value
##  1 - 2      -0.362 0.0146 12 -24.786 <.0001 
##  1 - 3      -0.129 0.0146 12  -8.825 <.0001 
##  2 - 3       0.233 0.0146 12  15.961 <.0001 
## 
## P value adjustment: tukey method for comparing a family of 3 estimates
plot(lsmS)

#Contrasts


#CH05 GLM - page 15: je ne vois pas comment obtenir les constrasts "pente1 - pente2" "pente1 - pente3"  "pente2 - pente3"

5 Modeles lineaires generalises

Objectif

Regression logistique ~> Expliquer une variable réponse binaire en fonction de variables explicatives quantitatives et / ou qualitatives + termes d’interaction.

Regression de Poisson ~> Expliquer une variable de comptage en fonction de variables explicatives quantitatives et / ou qualitatives + termes d’interaction.

5.1 Modelisation

#Régression logistique - X1 quantitatives ou qualitative, X2 quantitatives ou qualitative
mod <- glm( Y ~ X1 * X2, data = .., family = "binomial") 

#Régression de Poisson - X1 quantitatives ou qualitative, X2 quantitatives ou qualitative
mod <- glm( Y ~ X1 * X2, data = .., family = "poisson")

5.2 Fonctions utiles sur un objet de type mod <- glm()

Les mêmes fonctions que pour un objet de type mod <- lm() peuvent être utilisées sur un objet de type mod <- glm(). Ces fonctions ont la même utilité. Je vous renvoie donc vers la section 3.1.

5.3 Exemple 1: Régression logistique

5.4 Exemple 2: Régression de Poisson

6 Analyses multivariees

Objectif

  • ACP
    ~> Réduire le nombre de variables en transformant des variables corrélées en nouvelles variables décorrélées afin de mieux visualier les données.

  • Clustering
    ~> Diviser un ensemble de données en différents groupes.
    ~> Nombre de groupes fixé ou non à l’avance.
    ~> Les observations d’un même groupe sont supposées partager des caractéristiques communes.

6.1 Anaylse en composantes principales (ACP)

res <- FactoMineR::PCA(..)

Remarque
La fonction PCA() du package FactoMineR standarise automatiquement les données. Il n’est donc pas nécessaire de le faire à la main.


Valeurs propres / variances
Extraire factoextra::get_eig(res)
Visualiser factoextra::fviz_screeplot(res, ..)
Variables
Informations var <- factoextra::get_pca_var(res)
~> Coordonnées var$coord
~> Corrélation var$cor
~> Cos2 var$cos2
~> Contributions var$contrib
Cercle de corrélation factoextra::fviz_pca_var(res, ..)
Visualiser les contributions factoextra::fviz_contrib(res, ..)
Individus
Informations ind <- factoextra::get_pca_ind(res)
~> Coordonnées head(ind$coord)
~> Cos2 head(ind$cos2)
~> Contributions head(ind$contrib)
Carte des individus factoextra::fviz_pca_ind(res, ..)
Autres
Biplot variables-individus factoextra::fviz_pca_biplot(res, ..)

6.2 Exemple 1: ACP

#Visualiser les données
pairs(data[,-5], col = data$Species , main="Scatter plots des variables centrées réduites prises 2 à 2", oma=c(3,3,3,12), pch=20)
par(xpd = TRUE)
legend("bottomright", title = "Species", legend = levels(data$Species), col = c(1,2,3), cex = 0.8, pch = c(20,20,20))

#Analyse en composantes principales
res <- FactoMineR::PCA(data, graph=FALSE, quali.sup=5)
#Valeurs propres / variances
factoextra::get_eig(res)
##       eigenvalue variance.percent cumulative.variance.percent
## Dim.1 2.91849782       72.9624454                    72.96245
## Dim.2 0.91403047       22.8507618                    95.81321
## Dim.3 0.14675688        3.6689219                    99.48213
## Dim.4 0.02071484        0.5178709                   100.00000
#Valeurs propres / variances - Visualiser
factoextra::fviz_screeplot(res, addlabels = TRUE)

#Variables
var <- factoextra::get_pca_var(res) 
var
## Principal Component Analysis Results for variables
##  ===================================================
##   Name       Description                                    
## 1 "$coord"   "Coordinates for the variables"                
## 2 "$cor"     "Correlations between variables and dimensions"
## 3 "$cos2"    "Cos2 for the variables"                       
## 4 "$contrib" "contributions of the variables"
#Variables - Cercle de corrélation
factoextra::fviz_pca_var(res, axes = c(1, 2), col.var = "contrib")

#Variables
ind <- factoextra::get_pca_ind(res) 
ind
## Principal Component Analysis Results for individuals
##  ===================================================
##   Name       Description                       
## 1 "$coord"   "Coordinates for the individuals" 
## 2 "$cos2"    "Cos2 for the individuals"        
## 3 "$contrib" "contributions of the individuals"
#Variables - Carte des individus
factoextra::fviz_pca_ind(res, axes = c(1, 2), geom.ind = "point", col.ind = data$Species, addEllipses = TRUE, legend.title = "Species")

#Biplot variables - individus
factoextra::fviz_pca_biplot(res, repel = TRUE, col.var = "black", col.ind = "red", label = "var")

6.3 Clustering

#Clustering hiérarchique = CAH
res <- factoextra::hcut(x = .., k = .., hc_method = "single"/"complete"/"average"/"ward.D"/"ward.D2", hc_metric = "euclidean", stand = TRUE)

#Clustering non hiérarchique = K-Means
res <- kmeans(x = .., centers = .., nstart = 20)
CAH
Clustering res
Infos sur le clustering ~> Groupes res$cluster
~> Nb de groupes res$nbclust
~> Taille des groupes res$size
Visualiser le clustering factoextra::fviz_cluster(res)
Comparer le clustering avec les vrais catégories
si elles sont connues		 resPCA <- FactoMineR::PCA(..)
factoextra::fviz_pca_ind(resPCA, ..)
K-Means
Clustering res
Infos sur le clustering ~> Groupes res$cluster
~> Centres res$centers
~> Somme totale des carrés res$totss
~> Somme des carrés Within res$withinss
~> Somme totale des carrés Within res$tot.withinss
~> Somme des carrés Between res$betweenss
Visualiser le clustering factoextra::fviz_cluster(res, ..)
Comparer le clustering avec les vrais catégories
si elles sont connues		 resPCA <- FactoMineR::PCA(..)
factoextra::fviz_pca_ind(resPCA, ..)

6.4 Exemple 2: CAH

#Si on n'a aucune idée en combien de cluster diviser les données
factoextra::fviz_nbclust(data[,-5], factoextra::hcut, method = "wss")

#Clustering
res <- factoextra::hcut(data[,-5], k = 3, hc_method = "average", hc_metric = "euclidean", stand = TRUE)
#Dendrogramme
factoextra::fviz_dend(res, rect = TRUE, cex = 0.5)

#Groupes
grp <- as.factor(res$cluster)
grp
##  [1] 1 1 1 1 1 1 2 3 3 2 2 2 2 2 2 2 2 2
## Levels: 1 2 3
#Visualiser les groupes
factoextra::fviz_cluster(res)

#Comparaison clustering et vraies catégories - plots
resPCA <- FactoMineR::PCA(data, graph=FALSE, quali.sup=5)
factoextra::fviz_pca_ind(resPCA, col.ind = data$Species, addEllipses = TRUE, legend.title = "Species", title="Vraie classification")

factoextra::fviz_pca_ind(resPCA, col.ind = grp, addEllipses = TRUE, legend.title = "Cluster", title="Classification après CAH")
## Too few points to calculate an ellipse

6.5 Exemple 3: K-means

#Si on n'a aucune idée en combien de cluster diviser les données
factoextra::fviz_nbclust(data[,-5], kmeans, method = "wss")

#Clustering
res <- kmeans(data[,-5], centers = 3, nstart = 20)
#Groupes
grp <- as.factor(res$cluster)
grp
##  [1] 2 2 2 2 2 2 1 1 1 1 1 1 3 3 1 3 3 3
## Levels: 1 2 3
#Visualiser les groupes
factoextra::fviz_cluster(res, data[,-5])

#Comparaison clustering et vraies catégories
resPCA <- FactoMineR::PCA(data, graph=FALSE, quali.sup=5)
factoextra::fviz_pca_ind(resPCA, col.ind = data$Species, addEllipses = TRUE, legend.title = "Species", title="Vraie classification")

factoextra::fviz_pca_ind(resPCA, col.ind = grp, addEllipses = TRUE, legend.title = "Cluster", title="Classification après K-Means")

  1. ggplot2 ne sera pas utilisé ici mais il vaut le coup d’être mentionné et de s’y intéresser.↩︎