Analyse préliminaire avec R et STATA
Un article de IMSP - Formation continue.
Librement inspiré par: "Statistics with Computing" (2005), Taylor M, Kirkwood B, Tholen A, Borchert MWong J, Horrill S. LSHTM, University of London.
Tableau synoptique:
| opération | STATA | R |
|---|---|---|
| ouvrir le logiciel | (icône "Intercooled Stata") | R |
| configuration | version 8.0 | install.packages('foreign','pastecs','class', 'e1071','BSDA','sma') library('foreign') library('pastecs') library('BSDA') library(sma) options(scipen=1000) options(digits=3) help.start() |
| consulter l'aide | help nom_de_commande | help(nom_de_commande) ou help.search(nom) |
| importer des données au format CSV | insheet using test001.csv, names | test001 <- read.csv('test001.csv', header = TRUE) |
| sommaire | describe | str(test001) |
| liste de données (n premières lignes) | list in 1/n | head(naiss,n) |
| statistiques descriptives | summarize | stat.desc(test001) ou bien: summary(test001) sd(test001, na.rm = TRUE) |
| moyenne, médiane, écart-type pour certaines colonnes | summ p1ph1 summ p1ph1,detail | mean(test001$X.P1PH1) median(test001$X.P1PH1) sd(test001$X.P1PH1) |
| quantiles, maxima, minima | summarize p1ph1, detail | stat.desc(test001$X.P1PH1) quantile(test001$X.P1PH1, probs=c(.01,.05,.1,.25,.5,.75,.9,.95,.99)) p1ph1.sorted <- sort(test001$X.P1PH1) head(p1ph1.sorted) tail(p1ph1.sorted) ou max(test001$X.P1PH1), min(test001$X.P1PH1) |
| enregistrer les données | save test001 | save(test001,file = "test001.Rdata") |
| sortir du programme | exit | q('no') |
Les notions utilisées dans ce chapitre pratique sont explicitées dans les mesures de tendance centrale et les mesures de dispersion.
Sommaire |
Configuration préalable
Stata
nécessaire pour Stata 9.0 lorsque les instructions sont données au format 8.0
version 8.0
R
Nous devons charger quelques extensions, une seule fois par installation. La routine vous demandera de choisir l'endroit de téléchargement et installera automatiquement les paquets.
install.packages('foreign')
install.packages('pastecs')
install.packages('class')
install.packages('e1071')
install.packages('BSDA')
install.packages('sma')
vous devez charger en mémoire ces extensions pour pouvoir les utiliser. Vous exécuterez donc les commandes suivantes une fois par session (typiquement dans le script de commandes):
library('foreign')
library('pastecs')
library('BSDA')
configurons quelques défauts pour l'affichage des décimales:
options(scipen=1000) options(digits=3)
Les données
Le format classique dans lequel la plupart des logiciels s'attendent à trouver les données consiste en un tableau où les colonnes représentent les variables observées et les lignes les cas ou observations.
Soit un fragment de table de données provenant de travaux de recherche sur le paludisme au Mali. Le fichier est au format csv (comma-separated values [1]), très courant pour l'échange de données. Vous serez probablement amenés à utiliser ce format pour échanger des données avec les tableurs (Openoffice, Excel...)
http://www.santepublique.org/fc/data/test001.csv
Enregistrez-le (bouton de droite) dans le répertoire où vous avez l'intention de travailler avec R et/ou Stata.
Les données sont prêtes pour l'importation. Nous traiterons séparément de la préparation des données, qui dans la pratique courante consomme passablement de temps et de ressources.
La documentation qui accompagne ce fichier nous apprend que les variables, soit les colonnes du tableau, représentent les observations suivantes:
| variable | description |
|---|---|
| quartier | quartier de résidence |
| ddn | date de naissance |
| sexe | sexe |
| p1ph1 | lecture du premier hématocrite |
| p1ph2 | lecture du deuxième hématocrite |
| p1ptemp | température corporelle |
| p1prate | taille de la rate (score de Ridley) |
| p1pge | goutte épaisse positive |
| p1pdenpf | densité de plasmodium falciparum |
| p1pgam | présence ou absence de gamétocytes |
dans ce premier exercice nous allons surtout nous intéresser aux variables numériques continues de cette table. Est-ce que vous voyez desquelles il s'agit?
Importer les données
STATA
vérifions que le fichier de données se trouve bien dans le répertoire de travail:
ls
réponse:
35.9k 11/16/06 18:46 test001.csv
puis importez-le avec la commande insheet. L'option "names" indique à STATA que la première ligne contient les noms des variables.
insheet using test001.csv, name
réponse:
(10 vars, 754 obs)
R
Vérifions que le fichier de données se trouve bien dans le répertoire de travail
dir()
réponse:
[1] "test001.csv"
On lit l'objet et on l'attribue à un objet test001 de type "data frame"
test001 <- read.csv('test001.csv')
R reconnaît automatiquement la première ligne comme contenant les noms des variables si elle contient une colonne de moins que les lignes suivantes. si ce n'était pas le cas, il faudrait le préciser:
test001 <- read.csv('test001.csv', header = TRUE)
Voir aussi: Importer des données dans R
liste de données
STATA
liste des trois premières lignes
list in 1/3
réponse:
+------------------------------------------------------------------+
1. | quartier | ddn | sexe | p1ph1 | p1ph2 | p1ptemp | p1prate |
| ZONDE | 01.01.1989 | M | 34 | 37 | 37 | 0 |
|------------------------------------------------------------------|
| p1pge | p1pdenpf | p1pgam |
| N | 0 | N |
+------------------------------------------------------------------+
+------------------------------------------------------------------+
2. | quartier | ddn | sexe | p1ph1 | p1ph2 | p1ptemp | p1prate |
| ZONDE | 01.01.1992 | M | 30 | 28 | 36.4 | 0 |
|------------------------------------------------------------------|
| p1pge | p1pdenpf | p1pgam |
| Y | 33400 | Y |
+------------------------------------------------------------------+
+------------------------------------------------------------------+
3. | quartier | ddn | sexe | p1ph1 | p1ph2 | p1ptemp | p1prate |
| ZONDE | 01.01.1989 | F | 35 | 35 | 37.1 | 0 |
|------------------------------------------------------------------|
| p1pge | p1pdenpf | p1pgam |
| Y | 150 | Y |
+------------------------------------------------------------------+
R
liste des trois premières lignes:
head(test001,3)
réponse:
X.QUARTIER X.DDN X.SEXE X.P1PH1 X.P1PH2 X.P1PTEMP X.P1PRATE X.P1PGE 1 ZONDE 01.01.1989 M 34 37 37.0 0 N 2 ZONDE 01.01.1992 M 30 28 36.4 0 Y 3 ZONDE 01.01.1989 F 35 35 37.1 0 Y X.P1PDENPF X.P1PGAM 1 0 N 2 33400 Y 3 150 Y
Sommaire des données
STATA
instruction:
describe
réponse:
Contains data
obs: 754
vars: 10
size: 31,668 (97.0% of memory free)
-------------------------------------------------------------------------------
storage display value
variable name type format label variable label
-------------------------------------------------------------------------------
quartier str14 %14s QUARTIER
ddn str10 %10s DDN
sexe str1 %9s SEXE
p1ph1 byte %8.0g P1PH1
p1ph2 byte %8.0g P1PH2
p1ptemp float %9.0g P1PTEMP
p1prate byte %8.0g P1PRATE
p1pge str1 %9s P1PGE
p1pdenpf long %12.0g P1PDENPF
p1pgam str1 %9s P1PGAM
-------------------------------------------------------------------------------
Sorted by:
Note: dataset has changed since last saved
On remarque que les dates ont été interprétées comme des chaînes de caractères. Nous traiterons cet aspect plus tard, pour l'instant ce sont les champs numériques qui nous intéressent: ils ont été correctement reconnus comme tels.
instruction:
summarize
réponse:
Variable | Obs Mean Std. Dev. Min Max
-------------+--------------------------------------------------------
quartier | 0
ddn | 0
sexe | 0
p1ph1 | 754 34.39655 4.088079 19 52
p1ph2 | 705 34.95177 4.11282 20 53
-------------+--------------------------------------------------------
p1ptemp | 744 36.97513 .6246507 35.6 47.2
p1prate | 744 .2674731 .5867226 0 3
p1pge | 0
p1pdenpf | 715 3551.119 27096.24 0 666000
p1pgam | 0
Remarquons une température corporelle (p1ptemp) maximale à 47.2, manifestement une erreur de saisie. Les valeurs des hématocrites (p1ph1 et p1ph2) semblent plausibles, la taille de la rate (p1prate) est codée de 0 à 3, ce qui est corrrect.
R
instruction:
str(test001)
réponse:
'data.frame': 754 obs. of 10 variables: $ X.QUARTIER: Factor w/ 12 levels "FONAH","GOUNYAHOUN",..: 9 9 9 9 9 9 9 9 9 9 ... $ X.DDN : Factor w/ 288 levels "01.01.1989","01.01.1990",..: 1 4 1 6 4 7 83 1 18 24 ... $ X.SEXE : Factor w/ 2 levels "F","M": 2 2 1 2 1 1 1 1 1 2 ... $ X.P1PH1 : int 34 30 35 45 36 37 35 35 35 25 ... $ X.P1PH2 : int 37 28 35 48 36 36 42 38 35 29 ... $ X.P1PTEMP : num 37 36.4 37.1 37.2 37 37.2 36.2 37.2 36.6 37.5 ... $ X.P1PRATE : int 0 0 0 0 0 0 0 0 0 0 ... $ X.P1PGE : Factor w/ 3 levels "","N","Y": 2 3 3 2 3 3 3 3 2 3 ... $ X.P1PDENPF: int 0 33400 150 0 250 100 400 100 0 9600 ... $ X.P1PGAM : Factor w/ 3 levels "","N","Y": 2 3 3 2 2 2 2 2 2 2 ...
instruction:
summary(test001);
réponse:
X.QUARTIER X.DDN X.SEXE X.P1PH1 X.P1PH2
ZONDE :215 15.05.1999: 15 F:384 Min. :19.0 Min. :20.0
FONAH :172 15.02.1998: 12 M:370 1st Qu.:32.0 1st Qu.:32.0
SEREKENI : 89 15.01.1997: 11 Median :34.5 Median :35.0
ZONDE koko: 79 15.01.1998: 11 Mean :34.4 Mean :35.0
GOUNYAHOUN: 64 15.02.1992: 11 3rd Qu.:37.0 3rd Qu.:38.0
ZIANH : 64 15.03.1996: 11 Max. :52.0 Max. :53.0
(Other) : 71 (Other) :683 NA's :49.0
X.P1PTEMP X.P1PRATE X.P1PGE X.P1PDENPF X.P1PGAM
Min. :35.6 Min. : 0.000 : 39 Min. : 0 : 39
1st Qu.:36.7 1st Qu.: 0.000 N:252 1st Qu.: 0 N:620
Median :37.0 Median : 0.000 Y:463 Median : 200 Y: 95
Mean :37.0 Mean : 0.267 Mean : 3551
3rd Qu.:37.2 3rd Qu.: 0.000 3rd Qu.: 1600
Max. :47.2 Max. : 3.000 Max. :666000
NA's :10.0 NA's :10.000 NA's : 39
et si à ce stade on veut aussi consulter les déviations standard des variables numériques:
sd(test001, na.rm = TRUE)
ce qui nous donne:
X.QUARTIER X.DDN X.SEXE X.P1PH1 X.P1PH2 X.P1PTEMP X.P1PRATE
NA NA NA 4.088 4.113 0.625 0.587
X.P1PGE X.P1PDENPF X.P1PGAM
NA 27096.244 NA
autre possibilité: la fonction stat.desc() fournie par le paquet pastecs:
stat.desc(test001)
réponse:
X.QUARTIER X.DDN X.SEXE X.P1PH1 X.P1PH2 X.P1PTEMP X.P1PRATE
nbr.val NA NA NA 754.000 705.000 744.0000 744.0000
nbr.null NA NA NA 0.000 0.000 0.0000 596.0000
nbr.na NA NA NA 0.000 49.000 10.0000 10.0000
min NA NA NA 19.000 20.000 35.6000 0.0000
max NA NA NA 52.000 53.000 47.2000 3.0000
range NA NA NA 33.000 33.000 11.6000 3.0000
sum NA NA NA 25935.000 24641.000 27509.5000 199.0000
median NA NA NA 34.500 35.000 37.0000 0.0000
mean NA NA NA 34.397 34.952 36.9751 0.2675
SE.mean NA NA NA 0.149 0.155 0.0229 0.0215
CI.mean NA NA NA 0.292 0.304 0.0450 0.0422
var NA NA NA 16.712 16.915 0.3902 0.3442
std.dev NA NA NA 4.088 4.113 0.6247 0.5867
coef.var NA NA NA 0.119 0.118 0.0169 2.1936
X.P1PGE X.P1PDENPF X.P1PGAM
nbr.val NA 715.00 NA
nbr.null NA 302.00 NA
nbr.na NA 39.00 NA
min NA 0.00 NA
max NA 666000.00 NA
range NA 666000.00 NA
sum NA 2539050.00 NA
median NA 200.00 NA
mean NA 3551.12 NA
SE.mean NA 1013.34 NA
CI.mean NA 1989.49 NA
var NA 734206439.57 NA
std.dev NA 27096.24 NA
coef.var NA 7.63 NA
exercice: tester la fonction EDA() fournie par l'extension BSDA. Qu'est-ce qu'elle apporte un plus?
EDA(test001$X.P1PH1)
Quantiles, maxima et minima
STATA
instruction:
summarize p1ph1, detail
réponse:
P1PH1
-------------------------------------------------------------
Percentiles Smallest
1% 25 19
5% 28 22
10% 30 23 Obs 754
25% 32 24 Sum of Wgt. 754
50% 34.5 Mean 34.39655
Largest Std. Dev. 4.088079
75% 37 47
90% 40 48 Variance 16.71239
95% 41 48 Skewness .2427224
99% 45 52 Kurtosis 3.676212
R
instruction:
stat.desc(test001$X.P1PH1)
réponse:
nbr.val nbr.null nbr.na min max range
754.000 0.000 0.000 19.000 52.000 33.000
sum median mean SE.mean CI.mean.0.95 var
25935.000 34.500 34.397 0.149 0.292 16.712
std.dev coef.var
4.088 0.119
et les quantiles:
quantile(test001$X.P1PH1, probs=c(.01,.05,.1,.25,.5,.75,.9,.95,.99))
réponse:
1% 5% 10% 25% 50% 75% 90% 95% 99% 25.0 28.0 30.0 32.0 34.5 37.0 40.0 41.0 45.0
classons les données de X.P1PH1 par ordre croissant:
p1ph1.sorted <- sort(test001$X.P1PH1)
nous pouvons consulter les valeurs minimales:
> head(p1ph1.sorted) [1] 19 22 23 24 24 24
et maximales:
> tail(p1ph1.sorted) [1] 45 45 47 48 48 52
Enregistrer vos données et quitter le programme
nous voulons enregistrer la table dans le format natif de chacun des logiciels.
STATA
save test001
réponse:
file test001.dta saved
vérifions:
dir
réponse:
35.9k 11/16/06 18:46 test001.csv <dir> 11/16/06 23:48 .. <dir> 11/16/06 23:48 . 29.7k 11/16/06 23:48 test001.dta
sortir de STATA:
exit
R
save(test001,file = "test001.Rdata")
vérifions:
dir()
réponse:
[1] "test001.csv" "test001.Rdata"
Sortir de R sans enregistrer l'environnement:
q('no')
Références
Vincent Zoonekynd publie un site de référence très complet, orienté vers les graphiques: Statistics with R
Lectures recommandées
Aragón et Enanoria (2006) Epidemiology with R chap. 1: Getting started with R
R class notes: entering data (UCLA)
R class notes: exploring data (UCLA)
Stata for a PC or Mac: brief information (Columbia Univ.)
Introduction to Stata, Hands-on training exercises. Harvard-MIT data center
STATA Quick introduction (New Mexico Univ.)
Exercices
- Calculez les quantiles pour les variables p1ph2 et p1pdenpf
- téléchargez le fichier http://www.santepublique.org/fc/data/athletes.csv et importez-le dans votre logiciel de statistiques. Vous trouverez une description des variables et du contexte ici
- Que pouvez-vous nous dire de ces données? établissez un tableau avec:
- nombre de variables
- type des variables : quelles sont les variables numériques continues? les facteurs?
- nombre d'observations
- valeurs maximale et minimale pour chacune des variables numériques continues
- moyenne pour chacune des variables numériques continues
- médiane pour chacune des variables numériques continues
- déviation-standard pour chacune des variables numériques continues
- enregistrez la table de données dans un fichier athletes au format natif de votre logiciel. Nous l'utiliserons dans d'autres exercices. Gardez un régistre des opérations ("log") que vous avez du effectuer pour mener à bien cette analyse. Est-ce que vous lui trouvez une utilité?
- Que pouvez-vous nous dire de ces données? établissez un tableau avec:
Pour continuer
|
|
