Démarche statistique
remi.mahmoud@agrocampus-ouest.fr
https://demarche-stat-lesson.netlify.app
Objectifs
- Aborder les problèmes courants en analyse de données
- Argumenter le choix de procédures d’analyses
- Mettre en oeuvre une démarche d’analyse de données avec
![]()
- Interpréter & restituer les résultats d’une analyse
On attend:
Attention
Réflexion
Participation
Evaluation:
2 CC (50%) + 1 projet (50%)
Les statistiques sont (utilisées) partout
- En cela elles sont un outil politique, traduisant une certaine vision du monde et à manipuler avec précaution.
Coeur de la statistique
- Avoir les bonnes données:
- observations d’un phénomène
- sont-elles représentatives ?
Ex. où renforcer l’avion ?
Les résumer, les visualiser
Se poser des questions:
- Ex. pollution aux algues vertes: quelle variable utiliser ? Surface ? Masse ? Où/Quand/Comment faire des recueils ?
- Cause / Conséquence (variété de blé / Rendement ; utilisation de phyto / diversité entomologique ; complexité du paysage)
- Variables explicatives vs variable(s) réponse(s)
- Tous les effets ont-ils été pris en compte ?
Ex. paradoxe de Simpson (Vidéo complète sur le sujet de Science Etonnante)
Outil: R (moteur) / Rstudio (carrosserie)
- Gratuit/libre/multi-plateforme
- Nombreux (+++) domaines d’applications (ex. ce cours a été réalisé sous R)
- BEAUCOUP d’aide dispo
Installation:
- Télécharger et installer le moteur
![]()
\(\rightarrow\) https://cran.r-project.org/
- Télécharger et installer la carrosserie
![]()
On peut aussi dessiner avec R (mais c’est une autre histoire)
Quelques définitions (1/2)
Population: ensemble d’entités objet de l’investigation statistique
Individu: élément de la population d’étude
Echantillon: ensemble des individus pour lesquels des valeurs ont été observées pour les variables de l’étude
Variable: descripteur ou caractère des individus de la population d’étude
Inférence: décider pour une population à partir des données observées de l’échantillon
Décrire les données à l’aide d’indicateurs
Soit \(x_1, x_2, ..., x_n\) une série de \(n\) valeurs d’une variable X. On peut les décrire en utilisant des indicateurs de position et de dispersion (que vous connaissez probablement)
- moyenne: \(\color{blue}{\bar{x} = \frac{1}{n} \sum x_i = \frac{x_1 + x_2 + ... + x_n}{n}}\)
- médiane: \(\color{blue}{q_{0.5}(x) = }\) valeur telle que 50% des \(x_i\) ont une valeur inférieure et 50% une valeur supérieure
- 1er quartile: \(\color{blue}{q_{0.25}(x) = }\) valeur telle que 25% des \(x_i\) ont une valeur inférieure et 75% une valeur supérieure
- 3ème quartile: \(\color{blue}{q_{0.75}(x) = }\) valeur telle que 75% des \(x_i\) ont une valeur inférieure et 25% une valeur supérieure
- quantile \(\alpha\): \(\color{blue}{q_\alpha(x) = }\) valeur telle que \(100 \times \alpha\) % des \(x_i\) ont une valeur inférieure et \(100 - 100 \times \alpha\)% une valeur supérieure
- variance: \(\color{red}{s^2(x) = \hat{\sigma}^2(x) = \frac{1}{n-1} \sum (x_i - \bar{x})^2}\)
- écart-type: \(\color{red}{s(x) = \sqrt{s^2(x)}}\)
La visualisation: pourquoi ?
- Quels points communs entre ces 4 jeux de données ?
Visualisation: quels choix ?
- La visualisation dépend de:
La nature des données (nature des variables, nb d’individus)
La problématique: qu’est-ce que je veux montrer ?
le site data-to-viz.com permet de trouver de bonnes idées
- La visualisation permet de:
- comprendre / explorer / vérifier les données
- suggérer des analyses
- faire passer des idées
Les graphiques sont ce qu’on retient de nombreux rapports / analyses / articles scientifiques.
Sur 
, le package ggplot2 permet de faire des visualisations allant du plus simple au plus compliqué.
GGplot - Quick tuto
GGplot - Exemple: distribution d’une variable quantitative
GGplot - Exemple: distribution d’une variable quantitative
GGplot - Exemple: distribution d’une variable quantitative
GGplot - Exemple: distribution d’une variable quantitative
GGplot - Exemple: distribution d’une variable quantitative
GGplot - Exemple: distribution d’une variable quantitative
Moche ? Pensez à theme_bw() (black & white)
Distribution d’une variable qualitative (1/2)
Le camembert c’est bon, mais seulement en fromage. SINON CA PUE TROP.
Quelle classe est la plus représentée ?
Distribution d’une variable qualitative (2/2)
Alternative: graphiques en barre
Code
Code
library(dplyr)
ozone %>%
summarise(prop_vent = n()/nrow(ozone), .by = vent) %>% # On résume le jeu de donnée par direction de vent, on calcule le nb de chaque catégorie sur la longueur totale du jeu de données
mutate(vent = reorder(vent, prop_vent)) %>% # trier le facteur pour ordonner les catégories
ggplot(aes(x = vent, y= prop_vent)) +
labs(y = "%") +
geom_col(alpha = .5)+ # 50% de transparence
scale_y_continuous(labels=scales::percent) +
geom_text(aes(label=scales::percent(prop_vent)), vjust = -1.5)+
theme_bw()
Effet d’une / deux variable(s) sur une variable quanti
Effet d’une / deux variable(s) quali sur une variable quanti (1/2)
ggplot(ozone, aes(x = vent,
y= maxO3,
fill=vent,
col=vent)) + # dans aes: ce qui DEPEND du jeu de données
geom_boxplot(outlier.shape=NA,alpha=0.4) + # alpha : transparence
geom_point(position = position_jitter(width = 0.05, height= 0)) + # ajout
# points, avec perturbation horizontale pour faciliter visu
theme_bw() +
labs(title ="Effet du vent sur le max d'ozone")
Effet de deux variable(s) quali sur une variable quanti (2/2)
Quid d’une interaction potentielle vent pluie ?
Code
library(dplyr)
ozone %>%
summarise(moy_O3 = mean(maxO3), .by = c(vent, pluie)) %>%
# MOYENNER LES VALEURS DE MAXO3 PAR VENT ET PLUIE
#ET PASSER LE RESULTAT DANS GGPLOT VIA %>%
ggplot(aes(x = vent, y= moy_O3, group = pluie, color =pluie)) +
geom_line() +
geom_point() +
theme_bw() +
labs(title ="Interaction pluie:vent sur le max03",
y = "Moyenne de maxO3", x = "Direction du vent")
Des visualisations vers les tests statistiques
Les visus précédentes nous ont permis d’intuiter certaines tendances.
Peut-on les généraliser ? Rappel:
C’est toute la question des tests statistiques.
Rappel sur la notion de distribution statistique
Définition :
Une distribution statistique décrit comment les valeurs d’une variable se répartissent dans un ensemble de données.Types de variables :
- Quantitatives : mesurables (ex : taille des plantes).
- Qualitatives : catégories (ex : variétés de cultures).
Quelques paramètres descriptifs d’une distribution
- Tendance centrale :
- Moyenne : la “valeur attendue”.
- Médiane : la valeur au milieu de la distribution.
- Mode : la valeur la plus fréquente.
- Dispersion :
- Variance / écarte-type
- IQR / étendue
Illustration de la règle 68-95-99
Le TCL, c’est PUISSANT
- Ce que dit le TCL, c’est que peu importe la distribution de base, les moyennes empiriques d’échantillons issues cette distribution auront une distribution normale centrée sur l’espérance de la distribution !
Une simulation rapide pour visualiser les choses
- Ex distribution \(\chi_2(3)\)
- On connaît la vraie moyenne de cette loi: 3
Mais dans la vraie vie, on aurait accès qu’à un échantillon de cette loi, par ex. de taille \(n = 50\). Regardons comment évolue la distribution des moyennes empiriques en fonction de la taille de l’échantillon.
Dans notre cas, n = 112, donc \(T\) est censée suivre une loi de Student à 111 ddl, sous \(H_0\):
One Sample t-test
data: ozone$maxO3
t = -3.6406, df = 111, p-value = 0.0004148
alternative hypothesis: true mean is not equal to 100
95 percent confidence interval:
85.02578 95.58136
sample estimates:
mean of x
90.30357 [1] 0.0004147533Pour une loi de Student à n-1 ddl, -3.64 est une valeur peu probable donc on rejette \(H_0\), au seuil de 5%.
Ici, on a testé \(H_0\) vs H1 : \(\mu \neq 100\) (alternative = “two.sided”) mais on aurait pu testeer H1: \(\mu < 100\) (alternative = “less”) ou H1: \(\mu > 100\) (alternative = “greater”)
Comparaison de 2 moyennes (1/2)
Question: le temps (sec/pluvieux) a-t-il un effet sur le max d’ozone ?
Réflexe: visu !
Passer de l’échantillon à n’importe quel jour induit de l’incertitude. Mais, il est plus facile de conclure qu’il y a une différence entre les deux moyennes dans la population si:
- Les moyennes sont très différentes ( vs )
- La variabilité du maximum d’ozone est faible au sein des jours pluvieux et au sein des jours secs
- Il y a beaucoup de données
Test de comparaison de 2 variances
- \(H_0\) hypothèse nulle, \(H_0:\sigma_1^2=\sigma_2^2\) vs \(H_1\) hypothèse alternative, \(H1: \sigma_1^2\ne\sigma_2^2\)
\(\Longleftrightarrow~~H_0:\frac{\sigma_1^2}{\sigma_2^2}=1\) contre \(H_1:\frac{\sigma_1^2}{\sigma_2^2}\ne 1\)
- Statistique de test: \(F=\frac{S_1^2}{S_2^2}\)
- Sous \(H_0\), F suit une loi de Fisher à \(n_1 - 1\) et \(n_2 - 1\) ddl
Erreur et puissance d’un test
Décision d’un test prise pour la population (“il y a plus d’ozone les jours sec”) alors qu’on n’observe qu’un échantillon
\(\Rightarrow\) Décision incertaine, et deux erreur possibles
- Erreur de première espèce: rejeter \(H_0\) alors que celle-ci est vraie: faux positif
- Erreur de deuxième espèce : ne pas rejeter \(H_0\) alors que celle-ci est fausse: faux négatif
En image ( qu’est-ce que H0 ici ?):
un test est dit puissant si son erreur de 2ème espèce est petite. La puissance d’un test est la probabilité de rejeter \(H_0\) et d’avoir raison
Test de l’effet d’un facteur
Par ex. la direction du vent a-t-elle un effet sur le maximum d’ozone ?
Variable réponse quantitative, notée Y: maxO3
Variable explicative qualitative à I modalités (groupes): direction du vent
C’est le cadre de l’ANOVA (ANalysis Of VAriance)
maxO3 T9 T12 T15 Ne9 Ne12 Ne15 Vx9 Vx12 Vx15 maxO3v vent pluie
87 15.6 18.5 18.4 4 4 8 0.6946 -1.7101 -0.6946 84 Nord Sec
82 17.0 18.4 17.7 5 5 7 -4.3301 -4.0000 -3.0000 87 Nord Sec
92 15.3 17.6 19.5 2 5 4 2.9544 1.8794 0.5209 82 Est Sec
114 16.2 19.7 22.5 1 1 0 0.9848 0.3473 -0.1736 92 Nord Sec
94 17.4 20.5 20.4 8 8 7 -0.5000 -2.9544 -4.3301 114 Ouest Sec
80 17.7 19.8 18.3 6 6 7 -5.6382 -5.0000 -6.0000 94 Ouest PluieRéflexe …
Visu!
Question de base modifiée en : la moyenne du maximum d’ozone est-elle la même pour chaque direction du vent ?
Deux définitions du modèle
Si \(Y_{ij}\) est la valeur de la réponse du j\(^{eme}\) individu (\(j = 1, . . . , n_i\)) du i\(^{eme}\) groupe (\(i = 1, . . . , I\)):
\[\left\{ \begin{array}{l} \forall i,j~~~~Y_{ij}=\mu_i+ \varepsilon_{ij}\\ \forall i,j~~~~{\cal L}(\varepsilon_{ij})={\cal N}(0,\sigma^2)\\ \forall (i,j)\neq (i',j')~~cov(\varepsilon_{ij},\varepsilon_{i'j'})=0\\ \end{array}\right.\]
\(I\) paramètres: les \(\mu_i\)
\[\left\{ \begin{array}{l} \forall i,j~~~~Y_{ij}=\mu+\alpha_i+ \varepsilon_{ij}\\ \forall i,j~~~~{\cal L}(\varepsilon_{ij})={\cal N}(0,\sigma^2)\\ \forall (i,j)\neq (i',j')~~cov(\varepsilon_{ij},\varepsilon_{i'j'})=0\\ \end{array}\right.\]
\(I + 1\) paramètres: l’effet moyen \(\mu\) et les \(I\) coefficients \(\alpha_i\) (effet du niveau i)
On doit avoir \(\forall~i,~~\mu+\alpha_i=\mu_i\) besoin de contraintes
Décomposition de la variabilité
Equation d’analyse de la variance
| \(\displaystyle{\underbrace{\sum_{i,j}(Y_{ij}-Y_{\bullet \bullet})^2}_{SC_T}}\) | \(= \displaystyle{\underbrace{\sum_{i}(Y_{i\bullet}-Y_{\bullet \bullet})^2}_{SC_F}}\) | \(+ \displaystyle{\underbrace{\sum_{i,j}(Y_{ij}-Y_{i\bullet})^2}_{SC_R}}\) | |
| Variabilité | totale | modèle | résiduelle |
| ddl | \(n-1\) | \(I-1\) | \(n-I\) |
- Vous pouvez commencez à comprendre d’où vient le nom ANalysis Of VAriance
Indicateur de liaison : rapport de corrélation
\(\eta^2 = \frac{SC_{\mbox{modèle}}}{SC_{\mbox{total}}} = 1-\frac{SC_{\mbox{résiduelle}}}{SC_{\mbox{total}}}\)
Expliqué à mon grand-père (ou presque): la proportion de variabilité du phénomène que mon modèle explique
Propriétés:
- \(0 \leq \eta^2 \leq 1\) (obvious)
- \(\eta^2 = 0 \Leftrightarrow SC_{\mbox{modèle}}=0\)
- \(\eta^2=1 \Leftrightarrow SC_{\mbox{modèle}}=SC_{\mbox{total}}\)
Deux rapports de corrélation différents pour des \(Y_{i\bullet}\) identiques.
Inférence: test global
Rappel de l’objectif: y a-t-il un effet du facteur sur Y ?
\(\rightarrow\) La variabilité de Y est-elle expliquée par le facteur groupe ? Ou bien peut-on considérer que les données proviennt d’une même loi \(\mathcal{N}(\mu, \sigma^2)\) ?
Hypothèses:
\[H_0 :~\forall i,\ \mu_i=\mu\] \[H_1 :~\exists i\ /\ \mu_i\ne \mu\]
\(\Leftrightarrow\)
\[H_0 :~\forall i,\ \alpha_i=0\] \[H_1 :~\exists i\ /\ \alpha_i\ne 0\]
\[\mbox{On a : }\mathbb{E}\left(\frac{SC_{mod}}{I-1}\right)=\sigma^2+\frac{1}{I-1}\displaystyle{\sum_{i=1}^I n_i\alpha_i^2} \ \ \ \mathbb{E}\left(\frac{SC_{R}}{n-I}\right)=\sigma^2\]
Statistique de test : \(F_{obs}=\displaystyle{\frac{SC_{mod}/(I-1)}{SC_{R}/(n-I)}}\)
Loi de \(F_{obs}\) sous \(H_0\) : \({\mathcal L}(F_{obs})={\mathcal F}_{n-I}^{I-1}\)
Table d’analyse de variance: décomposition de variabilité et test
- Ce qui est souvent reporté dans le cadre d’une ANOVA
| Variabilité | Somme Carrés | ddl | Carré moyen | \(F_{obs}\) |
|---|---|---|---|---|
| Facteur | \(\displaystyle{\sum_in_i(Y_{i\bullet}-Y_{\bullet\bullet})^2}\) | \(I-1\) | \(\displaystyle{\frac{SC_F}{I-1}}\) | \(\displaystyle{\frac{CM_F}{CM_R}}\) |
| Résiduelle | \(\displaystyle{\sum_{i,j}(Y_{ij}-Y_{i\bullet})^2}\) | \(n-I\) | \(\displaystyle{\frac{SC_R}{n-I}}\) | |
| Totale | \(\displaystyle{\sum_{i,j}(Y_{ij}-Y_{\bullet \bullet})^2}\) | \(n-1\) |
Call:
LinearModel(formula = maxO3 ~ vent, data = ozone)
Residual standard error: 27.32 on 108 degrees of freedom
Multiple R-squared: 0.08602
F-statistic: 3.388 on 3 and 108 DF, p-value: 0.02074
AIC = 744.8 BIC = 755.7
Ftest
SS df MS F value Pr(>F)
vent 7586 3 2528.69 3.3881 0.02074
Residuals 80606 108 746.35
Ttest
Estimate Std. Error t value Pr(>|t|)
(Intercept) 94.7382 3.0535 31.0265 < 2e-16
vent - Est 10.8618 6.8294 1.5904 0.11466
vent - Nord -8.6092 4.6219 -1.8627 0.06522
vent - Ouest -10.0382 4.0972 -2.4500 0.01589
vent - Sud 7.7856 5.2052 1.4957 0.13764
Test de comparaison 2 à 2
Autre stratégie : comparer toutes les paires de moyennes
Pb : on effectue beaucoup de tests \(\Rightarrow\) risque de multiplier les erreurs en rejetant des hypothèses \(H_0\).
\(\Rightarrow\) correction des tests: modifier le seuil \(\alpha =5\%\) et prendre \(\alpha=\frac{5\%}{(\mbox{nb tests)}}\)
\(\Rightarrow\) les tests sont peu puissants
library(FactoMineR)
res <- LinearModel(maxO3 ~vent, data = ozone)
meansComp(res, ~vent, adjust ="Bonferroni", graph = FALSE)$adjMean
vent emmean SE df lower.CL upper.CL
Est 105.6 8.64 108 83.7 127.5
Nord 86.1 4.91 108 73.7 98.6
Ouest 84.7 3.86 108 74.9 94.5
Sud 102.5 5.96 108 87.4 117.7
Confidence level used: 0.95
Conf-level adjustment: bonferroni method for 4 estimates
$groupComp
$groupComp$Letters
Ouest Nord Sud Est
"a" "a" "a" "a"
$groupComp$LetterMatrix
a
Ouest TRUE
Nord TRUE
Sud TRUE
Est TRUE
attr(,"class")
[1] "meansComp"
Analyse des résidus du modèle
Rappel du modèle:
\[\left\{ \begin{array}{l} \forall i,j~~~~Y_{ij}=\mu_i+ \varepsilon_{ij}\\ \forall i,j~~~~{\cal L}(\varepsilon_{ij})={\cal N}(0,\sigma^2)\\ \forall (i,j)\neq (i',j')~~cov(\varepsilon_{ij},\varepsilon_{i'j'})=0\\ \end{array}\right.\]
On a des hypothèses sur les résidus du modèle
Les résidus, c’est ce qui n’est pas expliqué par le modèle. Donc pour les estimer, il faut ajuster le modèle
On peut ensuite vérifier les hypothèses du modèle
Liaison linéaire
Questions
- Influence de la température sur le max d’ozone ?
- A quel max d’ozone peut-on s’attendre s’il fait 19.5°C à 9h ?
Objectifs
Etudier qualitativement et quantitativement la dépendance d’une variable réponse quantitative Y en fonction d’une variable quantitative x
La variable x permet elle d’expliquer la variabilité de la variable Y?
Prédire Y à partir de x
Corrélation, causalité etc. (2/2)
Ex. pour rigoler:
“Corrélation n’implique pas causalité”
Ok, mais cette affirmation ne doit pas servir à réfuter toute association constatée dans des données.
Quelques valeurs de corrélations
Estimation des paramètres
\[SCER = \sum_{i=1}^n (y_i - \hat{y}_i)^2 = \sum_{i=1}^n (y_i - \beta_0 - \beta_1 x_i)^2\]
Paramètres d’espérance:
- \(\hat{\beta_1} = \frac{\sum_{i=1}^n (x_i-\bar x)(Y_i-\bar{Y})}{\sum_{i=1}^n (x_i - \bar{x})^2}\)
- \(\hat{\beta_0} = \bar{Y} - \hat{\beta_1}\bar{x}\)
Variance résiduelle:
\[ \hat{\sigma}^2 = \frac{1}{n-2} \sum_{i=1}^n (Y_i-\hat{Y}_i)^2 ~~~~~~ddl_{\mbox{résiduelle}}=n-2~~~~~~\mathbb{E}(\hat\sigma^2)=\sigma^2 \]
Décomposition de la variabilité
\[\begin{array}{ccccc} \sum_{i=1}^{n} ( Y_{i} - \bar{Y} )^{2} & = & \sum_{i=1}^{n} ( \hat{Y}_{i} - \bar{Y} )^{2} & + & \sum_{i=1}^{n} ( Y_{i} - \hat{Y}_{i} )^{2} \\ SCT & = & SCM & + & SCR \\ n - 1 & = & 1 & + & n - 2 \end{array}\]
Prédiction et intervalle de confiance
Pour une valeur de \(x_0\) particulière, on peut maintenant prédire \(Y\):
\[\hat{Y}_0 = \hat{\beta}_0 + \hat{\beta}_1 x_0\]
Prédiction de la valeur moyenne de Y pour un \(x_0\) particulier
\(\mathbb{E}(\hat{Y}_0|x_0) \sim \mathcal{N}(Y_0, \sigma \sqrt{\frac{1}{n} + \frac{(x_0 - \bar{x})^2}{\sum_{i=1}^n (x_i - \bar{x})^2)}})\)
Différent de:
Prédiction d’une nouvelle valeur de Y pour un \(x_0\) donné (notez le 1 supplémentaire dans la variance).
\(\hat{Y}_0 \sim \mathcal{N}(Y_0, \sigma \sqrt{1+\frac{1}{n} + \frac{(x_0 - \bar{x})^2}{\sum_{i=1}^n (x_i - \bar{x})^2)}})\)
Code
data_predict <- bind_cols(data.frame(T9 = ozone$T9), predict(lm(maxO3 ~ T9, data= ozone), interval = "predict"))
ozone %>%
ggplot(aes(x = T9, y = maxO3)) +
geom_point() +
geom_smooth(method ="lm") +# affiche intervalle de confiance par défaut
geom_ribbon(data = data_predict, aes(x = T9, ymax = upr, ymin = lwr, y = fit), alpha = .2, fill = "lightblue3") +
theme_bw()+
labs(title = "Intervalle de confiance du comportement moyen") +
theme(text = element_text(size = 15))
Validité du modèle
Vérification des hypothèses de départ:
- Normalité des résidus, stabilité de la variance, indépendance des résidus
Une variable quanti peut avoir un effet différent selon différentes modalités d’une quali
Cela s’écrit ainsi:
\[MaxO3 \sim vent + T9 + \mbox{vent:T9}\]
\[MaxO3_{ij} \sim \mu+\left\{ \begin{array}{l}
\alpha_{1} ~\mbox{si vent d'est} \\
\alpha_{2} ~\mbox{si vent du nord} \\
\alpha_{3} ~\mbox{si vent d'ouest} \\
\alpha_{4} ~\mbox{si vent du sud} \\
\end{array}
\right\} ~~+~~ \left(\beta+\left\{ \begin{array}{l}
\gamma_{1} ~\mbox{si vent d'est} \\
\gamma_{2} ~\mbox{si vent du nord} \\
\gamma_{3} ~\mbox{si vent d'ouest} \\
\gamma_{4} ~\mbox{si vent du sud} \\
\end{array}
\right\}\right)\times T9{ij}+ alea_{ij}\] \[\left\{
\begin{array}{l}
\forall i,j~~~~Y_{ij}=\mu+\alpha_i+(\beta+\gamma_i)\times x_{ij} + \varepsilon_{ij}\\
\forall i,j~~~~\varepsilon_{ij} ~~\mbox{i.i.d.}~~,~~\mathbb{E}(\varepsilon_{ij})=0,~~\mathbb{V}(\varepsilon_{ij})=\sigma^2\\
\forall i,j~~~~cov(\varepsilon_{ij},\varepsilon_{i'j'})=0\\
\end{array}\right.\]
Interaction entre deux variables qualitatives
Définition courante: réaction réciproque de deux phénomènes l’un sur l’autre
Définition statistique : l’effet d’un facteur sur Y diffère selon les modalités de l’autre facteur
library(dplyr)
ozone %>%
group_by(vent, pluie) %>%
summarize(MOY = mean(maxO3)) %>%
ggplot(aes(x=vent, y=MOY, col=pluie, group=pluie)) +
geom_line() +
geom_point() +
labs("Interaction pluie:vent sur maxO3", x = "Direction du vent", y = "Maximum d'ozone") +
theme_bw()+
theme(text = element_text(size = 15))
Cela s’écrit ainsi:
\[MaxO3 \sim vent + pluie + \mbox{vent:pluie}\] \[MaxO3_{ijk} \sim \mu+\left\{ \begin{array}{l}
\alpha_{1} ~\mbox{si vent d'est} \\
\alpha_{2} ~\mbox{si vent du nord} \\
\alpha_{3} ~\mbox{si vent d'ouest} \\
\alpha_{4} ~\mbox{si vent du sud} \\
\end{array}
\right\} + \left\{ \begin{array}{l}
\beta_{1} ~\mbox{si pluie} \\
\beta_{2} ~\mbox{si sec} \\
\end{array}
\right\} +
\left\{ \begin{array}{l}
\alpha\beta_{11} ~\mbox{si vent d'est ET pluie} \\
\alpha\beta_{12} ~\mbox{si vent d'est ET sec} \\
\alpha\beta_{21} ~\mbox{si vent du nord ET pluie} \\
\cdots
\alpha\beta_{42} ~\mbox{si vent du sud ET sec} \\
\end{array}
\right\} + alea_{ijk}\] \[\left\{
\begin{array}{l}
\forall i,j,k~~~~Y_{ijk}=\mu+\alpha_i+\beta_j +\alpha\beta_{ij} +\varepsilon_{ijk}\\
\forall i,j,k~~~~\varepsilon_{ijk} ~~\mbox{i.i.d.}~~,~~\mathbb{E}(\varepsilon_{ijk})=0,~~\mathbb{V}(\varepsilon_{ijk})=\sigma^2\\
\forall i,j,k~~~~cov(\varepsilon_{ijk},\varepsilon_{i'j'k'})=0\\
\end{array}\right.\]
Quatres types d’effet possibles
Code
theme_set(theme_bw())
ozone %>%
ggplot(aes(x = T9, y = maxO3)) +
geom_point() +
geom_smooth(method ="lm", se= FALSE) +
labs(title = "Quanti sur quanti")
ozone %>%
ggplot(aes(x = vent, y = maxO3, color = vent)) +
geom_point(position= position_jitter(width =.05, height= 0)) +
geom_boxplot(alpha= .3, outlier.shape = NA) +
labs(title = "Quali sur quanti") +
theme_bw()
Code
ozone %>%
ggplot(aes(x = T9, y = maxO3, color = vent))+
geom_point() +
geom_smooth(method = "lm", se = FALSE) +
labs(title = "Quali / quanti sur quanti") +
theme_bw()
ozone %>%
group_by(vent, pluie) %>%
summarize(MOY = mean(maxO3)) %>%
ggplot(aes(x=vent, y=MOY, col=pluie, group=pluie)) +
geom_line() +
geom_point() +
labs(title = "Quali / quali sur quanti") +
theme_bw()
Sélection de variables
Sélection de modèle = trouver un compromis entre un modèle qui s’ajuste bien aux données et qui n’a pas “trop” de paramètres (Rappel: )
- Souvent utilisé: sélection du modèle qui minimise l’AIC / BIC: compromis entre maximisation de la vraisemblance \(L\) (à quel point le modèle s’ajuste bien aux données) et le nb de paramètres
\(AIC = 2p - 2ln(\hat{L})\)
\(BIC = pln(n) - 2ln(\hat{L})\)
A votre avis, quelle est la différence entre les 2 critères ( !)
Plusieurs stratégies
- Construction exhaustive de tous les sous-modèles (long et même impossible si trop de variables)
- Méthode descendante (backward): construire le modéle complet; supprimer la variable explicative la moins intéressante et reconstruire le modèle sans cette variable; itérer jusqu’à ce que toutes les variables explicatives soient intéressantes
- Méthode ascendante (forward): partir du modéle avec la variable la plus intéressante; ajouter la variable qui, connaissant les autres variables du modèle, apporte le plus d’information complémentaire; itérer jusqu’é ce qu’aucune variable n’apporte d’information intéressante
- Méthode stepwise: compromis entre les 2 méthodes ci-dessus
Comparaison de moyennes ajustées
Moyennes ajustées = \(\hat{\mu} + \hat{\alpha}_i\): Permet de s’affranchir de l’effet des autres variables
Comparaisons possibles 2 à 2, moyennant une correction de Bonferroni par ex.
$adjMean
Produit emmean SE df lower.CL upper.CL
CoteOr_bio70 4.50 0.264 118 3.75 5.25
CoteOr_bio85 2.42 0.275 118 1.63 3.20
Equador72 3.88 0.264 118 3.12 4.63
Lindt70 4.79 0.264 118 4.04 5.55
Lindt78 4.29 0.264 118 3.54 5.05
Lindt85 2.58 0.264 118 1.83 3.34
Repere_bio74 3.88 0.275 118 3.09 4.66
Repere_bio85 2.46 0.264 118 1.70 3.21
Villars72 5.54 0.264 118 4.79 6.30
Villars85 3.17 0.264 118 2.41 3.92
Results are averaged over the levels of: Juge
Confidence level used: 0.95
Conf-level adjustment: bonferroni method for 10 estimates
$groupComp
CoteOr_bio85 Repere_bio85 Lindt85 Villars85 Equador72 Repere_bio74
"a" "a" "a" "ab" "bc" "bc"
Lindt78 CoteOr_bio70 Lindt70 Villars72
"bc" "cd" "cd" "d"
attr(,"class")
[1] "meansComp"Analyse des résidus
Les résidus sont la part du phénomène non-expliquée par le modèle
On peut vérifier les hypothèses initiales du modèle
Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 3 1.5589 0.2036
108 - Test de Shapiro-Wilk de normalité des Résidus (Rq : la non-normalité n’est pas un pb tant que la distribution est symétrique)
Shapiro-Wilk normality test
data: res
W = 0.97117, p-value = 0.01587Liste de messages d’erreurs communs
- Erreur dans file(file, “rt”) : impossible d’ouvrir la connexion:
- Nom du fichier mal spécifié
- Chemin mal spécifié (vérifiez le projet dans lequel vous êtes en haut à droite)
Erreur dans file(file, “rt”) : impossible d’ouvrir la connexion vers ‘https://r-stat-sc-donnees.github.io/ozone.tx’
- objet xxx introuvable \(\Rightarrow\) objet non chargé (vérifiez votre environnement)
- impossible de trouver la fonction xxx \(\Rightarrow\) Package non chargé / fonction mal orthographiée (R esT SeNSiBle A La CAsSe !)
Error in LinearModel(maxO3 ~ pluie, data = ozone): impossible de trouver la fonction "LinearModel"Error in Linearmodel(maxO3 ~ pluie, data = ozone): impossible de trouver la fonction "Linearmodel"