Analyse graphique avec R et STATA
Un article de IMSP - Formation continue.
Il a souvent été affirmé qu'une image vaut mieux qu'un long discours. Certains aspects que nous avons vu dans analyse préliminaire avec R et STATA nous apparaissent beaucoup plus évidents et parlants lorsqu'ils sont présentés sous forme graphique. Nous en abordons ici les rudiments.
Tableau synoptique:
| opération | STATA | R |
|---|---|---|
| configuration | version 8.0 | library(sma) options(scipen=1000) options(digits=3) |
| ouvrir le fichier de données | use test001 | load("test001.Rdata") attach('test001') |
| histogramme | histogram p1ph1,bin(17) graph display Graph,margin(tiny) scheme(s1mono) | hist(test001$X.P1PH1, breaks=17, col="light blue") |
| diagramme quantile-quantile ("q-q plot"): comparaison à une distribution normale | qnorm p1ph1, rlopts(lcolor(red)) graph display Graph, margin(tiny) scheme(s1color) | qqnorm(test001$X.P1PH1) qqline(test001$X.P1PH1, col="red",lwd=3) |
| diagramme quantile-quantile ("q-q plot"): comparaison entre deux distributions | qqplot p1ph1 p1ph2, rlopts(lcolor(red)) graph display Graph, margin(tiny) scheme(s1color) | qqplot(test001$X.P1PH2,test001$X.P1PH1) plot.qqline(test001$X.P1PH2,test001$X.P1PH1, col="red",lwd=3) |
| stripchart | N/A | stripchart(test001$X.P1PH1,method="jitter", xlim=c(20,55)) |
| corrélation et régression (graphique) | twoway (scatter p1ph2 p1ph1) (lfit p1ph2 p1ph1) (lfit p1ph1 p1ph2) | plot(test001$X.P1PH2 ~ test001$X.P1PH1) abline(lm(test001$X.P1PH2 ~ test001$X.P1PH1),col="red") abline(lm(test001$X.P1PH1 ~ test001$X.P1PH2),col="blue") |
| corrélation et régression (coefficient, test et intervalles de confiance) | correlate p1ph1 p1ph2 ou pwcorr p1ph1 p1ph2, sig obs | cor(test001$X.P1PH1,test001$X.P1PH2, use="complete") cor.test(test001$X.P1PH1,test001$X.P1PH2, alternative="two.sided", method="pearson", conf.level= 0.98) |
| diagramme de Tukey ("boîte à moustaches") | graph hbox p1ph2 p1ph1 | boxplot(test001$X.P1PH1,test001$X.P1PH2, horizontal = TRUE, notch = TRUE, col="light blue") |
Sommaire |
Configuration préalable
STATA
Version 8.0 use test001
R
Nous supposons que les extensions spécifiées dans analyse préliminaire avec R et STATA ont été installées. Dans ce cas il suffit de charger en mémoire celles qui nous seront nécessaires, et de configurer quelques paramètres d'affichage par défaut:
library(sma) options(scipen=1000) options(digits=3)
chargeons le fichier test001.Rdata produit dans la session précédente:
load("test001.Rdata")
objects()
réponse:
[1] "test001"
avec la fonction attach() nous incorporons l'objet test001 comme objet par défaut: nous n'aurons en principe plus besoin de préciser le nom de l'objet dans le nom des variables: test001$X.P1PH1 peut alors être identifiée comme test001$X.P1PH1 ou X.P1PH1 Si cette fonction s'avère pratique, il faut rester prudent en particulier si on travailole avec plusieurs objets qui utilisent des noms de variables similaires.
attach('test001')
Histogramme
L'histogramme est un graphe qui représente la fréquence d'ocurrence absolue ou relative d'une variable divisée en classes d'amplitudes égales. Il nous renseigne de façon immédiate sur les caractéristiques de distribution de cette variable.
STATA
histogram p1ph1, bin(17) frequency
Stata appelle cet objet "Graph", nous pouvons ensuite en modifier les couleurs:
graph display Graph, margin(tiny) scheme(s1mono)
on remarquera que par défaut STATA fait des histogrammes de densité au lieu d'histogrammes de fréquence:
histogram p1ph1, bin(17) graph display Graph, margin(tiny) scheme(s1mono)
ajoutons une distribution normale pour comparaison:
histogram p1ph1, bin(17) normal normopts(lpattern(dash) lcolor(red) lwidth(medthick))
et additionnellement une courbe de densité de la distribution empirique:
histogram p1ph1, bin(17) normal normopts(lpattern(dash) lcolor(red) lwidth(medthick)) kdensity
R
Nous utilisons la fonction hist(). Consultons le manuel:
start.help() help(hist)
Nous souhaitons visualiser l'histogramme de la variable X.P1PH1 du tableau test001:
hist(test001$X.P1PH1, breaks=17, col="light blue")
nous obtenons un histogramme de fréquences:
Nous pouvons aussi calculer un histogramme de probablilités. Observer le changement dans l'axe des ordonnées. L'aire totale de l'histogramme égale 1.
hist(test001$X.P1PH1, breaks=17, col="light blue", probability = TRUE)
Il devient alors possible de comparer la distribution de notre variable à celle d'une loi normale. Nous devons pour cela calculer une distribution normale avec la même moyenne et déviation standard que nos données expérimentales:
f <- function(x) {
dnorm(x,
mean=mean(test001$X.P1PH1),
sd=sd(test001$X.P1PH1),)
}
puis nous l'appliquons au graphique généré au préalable:
curve(f,add=T,col="red",lwd=3,lty=2)
avec le résultat suivant:
Nous pouvons alors surimposer la fonction de densité des données d'observation:
lines(density(test001$X.P1PH1), col='red', lwd=3)
Cette représentation a l'avantage d'être intuitive et didactique, mais la méthode de calcul décrite ici fonctionne mal lorsqu'il y a des données manquantes (essayez de reproduire le résultat avec la variable X.P1PH2 par exemple). L'interprétation laisse beaucoup à la subjectivité. Le recours aux diagrammes quantile-quantile est recommandé dans ce cas.
Diagramme quantile-quantile (q-q plot)
Dans la pratique, lorsqu'on veut représenter la comparaison d'une fonction de répartition empirique avec une distribution normale, on utilise le diagrame quantile-quantile ou "q-q plot": "diagramme représentant la fonction de répartition empirique, gradué selon les fractiles de la loi normale réduite". Plus la répartition empirique se rapproche de la répartition d'une loi normale, plus les points apparaissent alignés le long d'une droite appelée droite de Henry.
Malgré leur attrait, les méthodes graphiques ne nous fournissent pas de critères véritablement objectifs pour juger de la normalité d'une distribution de données expérimentales[1]. De ce fait, il est important de savoir mettre en oeuvre les tests de normalité chaque fois que la question se pose.
STATA
La fonction qnorm nous permet de comparer une distribution empirique à une distribution selon la loi normale:
qnorm p1ph1, rlopts(lcolor(red)) graph display Graph, margin(tiny) scheme(s1color)
Le diagramme quantile-quantile est utile lorsqu'on cherche à répondre à des questions telles que:
- est-ce que les échantillons proviennent de populations avec une distribution semblable?
- est-ce que les échantillons ont des profils de distribution similaires?
- est-ce que les variables périphériques ont un comportement similaire?
Comparons par exemple les distributions du premier et du deuxième hématocrites p1ph1 et p1ph2
qqplot p1ph1 p1ph2, rlopts(lcolor(red)) graph display Graph, margin(tiny) scheme(s1color)
R
qqnorm(test001$X.P1PH1) qqline(test001$X.P1PH1,col="red", lwd=3)
Le diagramme quantile-quantile est utile lorsqu'on cherche à répondre à des questions telles que:
- est-ce que les échantillons proviennent de populations avec une distribution semblable?
- est-ce que les échantillons ont des profils de distribution similaires?
- est-ce que les variables périphériques ont un comportement similaire?
Comparons par exemple les distributions du premier et du deuxième hématocrites X.P1PH1 et X.P1PH2. Nous utilisons pour cela la fonction qqplot(). La fonction plot.qqline() qui nous permet de dessiner la ligne de référence est fournie par l'extension sma:
qqplot(test001$X.P1PH2,test001$X.P1PH1) plot.qqline(test001$X.P1PH2,test001$X.P1PH1,col="red",lwd=3)
En l'ocurrence, les deux distributions sont légèrement différentes.
par(mfrow=c(2,1)) ; segmente en deux lignes le dispositif graphique hist(test001$X.P1PH1, breaks=17, col="light blue", probability = TRUE,
xlim = c(20,50), main = "") hist(test001$X.P1PH2, breaks=17, col="light blue", probability = TRUE,
xlim = c(20,50), main = "") par(mfrow=c(1,1)) ; retourne le dispositif graphique à la configuration antérieure
Nous devrons tester si cette différence est significative.
Mais surtout l'existence d'une tracé de densité bimodal associé à la non normalité des données doit nous faire suspecter un mélange de populations avec des caractéristiques différentes. Voici ce qu'on observe si on sépare les enfants de 1 à 6 ans (en noir) des enfants de 6 à 12 ans (en rouge):
(pour la méthode utilisée pour élaborer ce graphique voir ici)
"Stripcharts"
STATA
N/A
R
stripchart(test001$X.P1PH2,method="jitter",xlim=c(20,55)) stripchart(test001$X.P1PH1,method="jitter",xlim=c(20,55))
très moyennement utile dans ce cas:
Graphe de fréquences cumulées
Ce type de graphe est parfois utilisé pour inspecter la distribution des données d'une variable, tout comme le graphe quantile-quantile.
STATA
cumul p1ph1, generate(cp1ph1) line cp1ph1 p1ph1, sort
R
Nous allons reprendre une fonction ad hoc de Vincent Zoonekynd pour générer un graphe de effectifs cumulés. Copier ce texte et le coller dans la libne de commandes de R:
cumulated.frequencies <- function (x, main="") {
x.name <- deparse(substitute(x))
n <- length(x)
plot( 1:n ~ sort(x),
xlab = x.name,
ylab = 'Fréquences cumulées',
main = main
)
}
nous voulons le graphe de effectifs cumulés de la variable X.P1PH1:
cumulated.frequencies(X.P1PH1, main = "p1ph1")
Corrélation et régression
STATA
distribution de p1ph2 par rapport à p1ph1 avec les deux droites de régression.
twoway (scatter p1ph2 p1ph1) (lfit p1ph2 p1ph1) (lfit p1ph1 p1ph2)
examinons le coefficient de corrélation:
correlate p1ph1 p1ph2
résultat:
(obs=705)
| p1ph1 p1ph2
-------------+------------------
p1ph1 | 1.0000
p1ph2 | 0.9419 1.0000
autre possibilité, avec niveau de signifiance et décompte des observations appariées:
pwcorr p1ph1 p1ph2, sig obs
résultat:
| p1ph1 p1ph2
-------------+------------------
p1ph1 | 1.0000
|
| 754
|
p1ph2 | 0.9419 1.0000
| 0.0000
| 705 705
|
R
plot(test001$X.P1PH2 ~ test001$X.P1PH1)
et les deux droites de régression:
abline(lm(test001$X.P1PH2 ~ test001$X.P1PH1),col="red") abline(lm(test001$X.P1PH1 ~ test001$X.P1PH2),col="blue")
Les deux séries d'observations sont hautement corrélées.
Vérifions le coefficient de corrélation:
cor(test001$X.P1PH1,test001$X.P1PH2, use="complete") [1] 0.942
et avec le test de corrélation:
cor.test(test001$X.P1PH1,test001$X.P1PH2, alternative="two.sided",
method="pearson", conf.level= 0.98)
réponse:
Pearson's product-moment correlation
data: test001$X.P1PH1 and test001$X.P1PH2
t = 74.3, df = 703, p-value < 0.00000000000000022
alternative hypothesis: true correlation is not equal to 0
98 percent confidence interval:
0.931 0.951
sample estimates:
cor
0.942
Diagramme de Tukey ("boxplot", "boîte à moustaches")
STATA
graph hbox p1ph2 p1ph1
R
boxplot(test001$X.P1PH1,test001$X.P1PH2, horizontal = TRUE, notch = TRUE,
col="light blue")
Exercices
- générez un histogramme, un diagramme de Tukey et un diagramme quantile-quantile pour la variable X.P1PTEMP (ou p1ptemp). Que remarquez-vous?
Références
Vincent Zoonekind: Statistics with R
Statistical graphics and experimental data Nikolic-Doric et coll. 2006
STATA Tutorial. Carolina Population Center
Introduction to Stata Graphics with associated Stata code and output (U. Michigan)
Dartmouth Social Science Computing: Frequently asked questions on Stata
Lectures recommandées
Park HM (2004) Testing normality in SAS, STATA, and SPSS
Pour continuer
|
|
