8 février 2021

Plan du cours

  • tableau d’effectifs et fréquence
    • calculer un effectif
    • graphiques univariés
  • mesures et tests d’association
    • tableau de contingence
    • comparaison de deux moyennes
    • corrélation linéaire
  • modèles statistiques
    • Analyse de variance à un facteur
    • Régression linéaire
    • Régression logistique

Présentation du jeu de données

Il s’agit d’un fichier sur la santée mentale en prison

smp <- read.csv("smp1.csv", header = TRUE, sep = ";", encoding = "UTF-8", na.strings = "NA")
str(smp)
## 'data.frame':    799 obs. of  6 variables:
##  $ age      : int  31 49 50 47 23 34 24 52 42 45 ...
##  $ prof     : chr  "autre" NA "prof.intermediaire" "ouvrier" ...
##  $ dep.cons : int  0 0 0 0 1 0 1 0 1 0 ...
##  $ scz.cons : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ grav.cons: int  1 2 2 1 2 1 5 1 5 5 ...
##  $ n.enfant : int  2 7 2 0 1 3 5 2 1 2 ...
head(smp)
##   age               prof dep.cons scz.cons grav.cons n.enfant
## 1  31              autre        0        0         1        2
## 2  49               <NA>        0        0         2        7
## 3  50 prof.intermediaire        0        0         2        2
## 4  47            ouvrier        0        0         1        0
## 5  23        sans emploi        1        0         2        1
## 6  34            ouvrier        0        0         1        3

Interroger les données

  • Variables quantitatives :
table(smp$n.enfant)
## 
##   0   1   2   3   4   5   6   7   8   9  10  11  13 
## 214 220 125 101  55  31   7   7   7   2   2   1   1
table(smp$n.enfant > 5)
## 
## FALSE  TRUE 
##   746    27

Interroger les données

  • Variables qualitatives
table(smp$prof)
## 
##        agriculteur            artisan              autre              cadre 
##                  6                 90                 31                 24 
##            employe            ouvrier prof.intermediaire        sans emploi 
##                135                227                 58                222
table(smp$prof == "agriculteur")
## 
## FALSE  TRUE 
##   787     6

Interroger les données

  • Variables qualitatives
which(smp$prof == "agriculteur") # indique les numéros d'observations
## [1]  15 312 384 391 439 442
smp$age[which(smp$prof == "agriculteur")] # indique l'âge des personnes dont la profession est agriculteur
## [1] 64 42 37 36 35 79
subset(smp, prof == "agriculteur", age) # selectionne des observations remplissant une condition
##     age
## 15   64
## 312  42
## 384  37
## 391  36
## 439  35
## 442  79
subset(smp, prof == "agriculteur", c(age, prof, n.enfant))
##     age        prof n.enfant
## 15   64 agriculteur        3
## 312  42 agriculteur        3
## 384  37 agriculteur        2
## 391  36 agriculteur        3
## 439  35 agriculteur        0
## 442  79 agriculteur        5
subset(smp, prof == "agriculteur" & n.enfant > 2, c(age, prof, n.enfant)) 
##     age        prof n.enfant
## 15   64 agriculteur        3
## 312  42 agriculteur        3
## 391  36 agriculteur        3
## 442  79 agriculteur        5
# les observations doivent remplir deux conditions

Tableaux d’effectifs et fréquences relatives

Puisque très peu de détenus ont plus de 5 enfants, il est possible de les regrouper :

table(smp$n.enfant > 5)
## 
## FALSE  TRUE 
##   746    27
smp$n.enfant.cat <- factor(smp$n.enfant) # pour créer une nouvelle variable factorielle
levels(smp$n.enfant.cat)[6:13] <- "5+" # pour transformer les niveaux 6 à 13 en 5 et plus
nlevels(smp$n.enfant.cat) # pour voir le nombre de niveaux dans la variable
## [1] 6
tab <- table(smp$n.enfant.cat)
tab
## 
##   0   1   2   3   4  5+ 
## 214 220 125 101  55  58

Tableaux d’effectifs et fréquences relatives

En stockant le résultat dans un nouvel objet, il est alors possible de faire de nouvelles opérations dessus

sum(tab) # pour avoir la somme d'enfants
## [1] 773
tab / sum(tab) # pour avoir la fréquence relative d'enfants
## 
##          0          1          2          3          4         5+ 
## 0.27684347 0.28460543 0.16170763 0.13065977 0.07115136 0.07503234
prop.table(tab) # fonction dédiée
## 
##          0          1          2          3          4         5+ 
## 0.27684347 0.28460543 0.16170763 0.13065977 0.07115136 0.07503234
prop.table(tab) * 100 # pour l'avoir en pourcentage)
## 
##         0         1         2         3         4        5+ 
## 27.684347 28.460543 16.170763 13.065977  7.115136  7.503234
round(prop.table(tab) * 100, 2) # pour l'arrondir à deux décimales près
## 
##     0     1     2     3     4    5+ 
## 27.68 28.46 16.17 13.07  7.12  7.50

Tableaux d’effectifs et fréquences relatives

Ces opérations fonctionnent aussi avec des variables qualitatives

tab_q <- table(smp$prof)
tab_q
## 
##        agriculteur            artisan              autre              cadre 
##                  6                 90                 31                 24 
##            employe            ouvrier prof.intermediaire        sans emploi 
##                135                227                 58                222
sum(tab_q) # somme de détenus qui ont un travail
## [1] 793
table(smp$prof, useNA = "always") # pour utiliser les données manquantes
## 
##        agriculteur            artisan              autre              cadre 
##                  6                 90                 31                 24 
##            employe            ouvrier prof.intermediaire        sans emploi 
##                135                227                 58                222 
##               <NA> 
##                  6
round(prop.table(tab_q) * 100, 1)
## 
##        agriculteur            artisan              autre              cadre 
##                0.8               11.3                3.9                3.0 
##            employe            ouvrier prof.intermediaire        sans emploi 
##               17.0               28.6                7.3               28.0
round(prop.table(table(smp$prof, exclude = "sans emploi")) * 100)
## 
##        agriculteur            artisan              autre              cadre 
##                  1                 16                  5                  4 
##            employe            ouvrier prof.intermediaire               <NA> 
##                 23                 39                 10                  1

Tableaux de contingence

En croisant des modalités de deux variables qualitatives, la fonction table() permet aussi de faire des tableaux de contingence. La première variable mentionnée sera indiquée en lignes dans le tableau généré.

tab_c <- table(smp$dep.cons, smp$scz.cons)
tab_c
##    
##       0   1
##   0 451  31
##   1 282  35

Tableaux de contingence

Dans le cas d’un tableau à deux entrées, les fréquences relatives affichées sont calculées par rapport à l’effectif total.

round(prop.table(tab_c)*100)
##    
##      0  1
##   0 56  4
##   1 35  4

visualisation des données

Pour une meilleure visualisation du résultat, il est possible d’en faire un graphique

plot(tab)

visualisation des données

On peut aussi utiliser une fonction spécifique pour faire des graphiques en barre

barplot(tab)

visualisation des données

Il est possible d’améliorer le rendu esthétique du graphique :

barplot(prop.table(tab) * 100, ylab = "proportion", xlab = "Nb d'enfants", col = "steelblue2", 
        border = NA, ylim = c(0, 30), las = 1, main = "Proportion d'enfants par détenus")

Attention : ce type de commande ne fonctionne pas avec une variable, mais un tableau d’effectifs ou de fréquence.

visualisation des données

La fonction barplot() peut aussi s’utiliser sur un tableau de contingence en rajoutant l’option beside = TRUE. Il est aussi possible d’y ajouter une légende pour indiquer les différentes modalités.

barplot(tab_c, beside = TRUE, legend = c("schizo = Non", "schizo = Oui"),
        xlab = " Dépression", ylab = "Effectifs", col = c("palegreen3", "salmon3"))

visualisation des données

Il existe de nombreux autres types de visualisation possibles :

Le camembert pour représenter la distribution d’une variable aléatoire catégorielle

pie(table(smp$prof))

visualisation des données

L’histogramme pour représenter la distribution d’une variable aléatoire quantitative continue

hist(smp$age, col = "grey", main = "Répartition des détenus \n en fonction de leur age", xlab = "age")

visualisation des données

Il est aussi possible d’utiliser des courbes de densité pour représenter une distribution en fonction d’une autre variable. Elles presentent l’avantage de pouvoir être superposées dans le même graphique.

plot(density(smp$age[smp$dep.cons == 1], na.rm = TRUE), lwd = 2, 
     col = "orange", xlim = c(19,83), main = "")
lines(density(smp$age[smp$dep.cons == 0], na.rm = TRUE), lwd = 2, 
      col = "deepskyblue4")

visualisation des données

La boîte à moustache pour représenter de façon plus synthétique la distribution d’une variable aléatoire quatitative

boxplot(smp$age, xlab = "age")

visualisation des données

Les boîtes à moustaches sont réellement utiles pour représenter graphiquement la distribution d’une variable quantitative en fonction de sous-groupes

library(RColorBrewer)
coul <- RColorBrewer::brewer.pal(6, "PuRd") # pour créer un vecteur de couleurs

visualisation des données

boxplot(age ~ n.enfant.cat, data = smp, frame = FALSE, xlab = "nb d'enfants",
        col = coul)

visualisation des données

Les graphiques cartésiens (ou graphique en x et y) permettent de représenter la distribution conjointe de deux variables aléatoires quantitatives

plot(smp$age, smp$n.enfant)

visualisation des données

plot(jitter(smp$age), jitter(smp$n.enfant))

La librairie ggplot2

C’est une librairie dédiée à la visualisation qui permet de réaliser des graphiques plus perfectionnés et au design plus moderne.

library(ggplot2)
qplot(age, data = smp, geom = "histogram", fill = I("orange"), colour = I("black"))

# qplot étant la fonction la plus basique de la librairie

La librairie ggplot2

La librairie ggplot2

On peut aussi faire des nuages de points

ggplot(smp, aes(x = age, y = n.enfant.cat)) +
  geom_point() # la fonction ggplot permet plus de customisation
## Warning: Removed 2 rows containing missing values (geom_point).

La librairie ggplot2

Ou des graphiques en barre

ggplot(smp, aes(prof)) +
  geom_bar()

La librairie ggplot2

Il est aussi possible de faire des comparaisons d’une variable en fonction de deux autres

ggplot(smp, aes(prof)) +
  geom_bar() +
  facet_grid(dep.cons ~ n.enfant.cat) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = rel(0.8)))

La librairie ggplot2

ggplot2 permet aussi de faire des graphiques plus complexes, et parce qu’ils sont extrêmement customisables, de les mettre en forme pour l’impression, directement dans R.

# ce graphique est stocké dans l'objet "ga"
ga  <- ggplot(p_eff2[[1]], aes(x = ID_SITE, y = CT2, color = site)) +
  # p_eff2 correspond aux données qui détaillent des types d'objets par sites
  
    # mise en forme des titres et des axes
  
  labs(size = "Effectifs", color = "", title = periode_names[[1]], subtitle =
         "Répartition des types par sites") +
  xlab("Sites") + 
  ylab("Types et sous-types d'objets") + 

  # ajout d'annotation pou faciliter la lecture
  
  annotate("rect", xmin = -Inf, xmax = Inf, ymin = first(o_list2[[1]][[1]]), 
           ymax = "DBA11000", fill = "coral1", alpha = 0.3) +
  annotate("text", x = "SIT0940", y = "BbG7a000", label = "Bagues", size = 
             rel(2), colour = "grey23", fontface = "bold") +
  annotate("rect", xmin = -Inf, xmax = Inf, ymin = "DBA11000", ymax = 
             "DBHH000T", fill = "darkslategray3", alpha = 0.3) +
  annotate("text", x = "SIT0940", y = "DBG3B00u", label = "Bracelets", size = 
             rel(2), colour = "grey23", fontface = "bold") +
  annotate("segment", x = -Inf, xen = Inf, y = "ECAA2000", yend = "ECAA2000", 
           colour = "gold2", 
           alpha = 0.3, size = 2) +
  annotate("text", x = "SIT0940", y = "ECAA2000", label = "Ceintures", size = 
             rel(2), colour = "grey23", fontface = "bold") +
  annotate("rect", xmin = -Inf, xmax = Inf, ymin = first(o_list2[[1]][[4]]), 
           ymax = last(o_list2[[1]][[4]]), fill = "navajowhite4", alpha = 0.3) +
  annotate("text", x = "SIT0940", y = "IMRâR000", label = "Monnaies", size = 
             rel(2), colour = "grey23", fontface = "bold") +
  annotate("segment", x = -Inf, xen = Inf, y = first(o_list2[[1]][[5]]), 
           yend = last(o_list2[[1]][[5]]), colour = "chartreuse2", 
           alpha = 0.3, size = 2) +
  annotate("text", x = "SIT0940", y = "KpA11000", label = "Pendentifs", size =
             rel(2), colour = "grey23", fontface = "bold") +
  annotate("rect", xmin = -Inf, xmax = Inf, ymin = "KpA11000", ymax = "LPZS400B",
           fill = "Khaki2", alpha = 0.3) +
  annotate("text", x = "SIT0940", y = "LPZA000A", label = "Perles", size = 
             rel(2), colour = "grey23", fontface = "bold") +
  
  # Détermine le type de graphique
  ## ici = nuage de point en fonction de deux variables non continues
  
  geom_count(alpha = .6) +
  
  # déterminer l'échelle et ses informations
  
  scale_size(breaks = c(1, 15, 30), range = c(1, 4), labels = 
               c("1", "<15", "15<30")) +
  scale_color_manual(breaks = c("1", "2"), values = c("grey28", "red3", "blue3"),
                     labels = c("Types présents \n sur une majorité \n de sites",
                                "Sites avec une \n majorité de types")) +
  scale_y_discrete(breaks = c("DBGcB004", "DBH3a00B", "LPH2000t"), labels =
                     c("Gebhard s. 38", "Haev. 3a", "Haev. 20")) +
  scale_x_discrete(breaks = c("SIT0001","SIT0086", "SIT0095", "SIT0175"), 
                   labels = c("Wederath","La Bure", "Bâle-Gasfabrik", 
                              "Kirchzarten")) +
  
  # mise en page du graphique et de la légende
  
  theme_minimal() +
  theme(axis.text = element_text(size = rel(0.6), face = 2, hjust = 1)) + 
  theme(axis.text.x = element_text(angle = 90)) + 
  theme(axis.title = element_text(size = rel(0.8), face = "bold")) + 
  theme(legend.box.background = element_rect(fill = "white", colour = "grey28")) +
  theme(legend.text = element_text(size = rel(0.6))) +
  theme(legend.title = element_text(size = rel(0.6))) +
  theme(legend.position = "top", legend.justification = "right") +
  guides(colour = guide_legend(override.aes = list(size = 3))) +
  theme(plot.title = element_text(size = rel(1.2), face = "bold")) +
  theme(plot.subtitle = element_text(size = rel(1)))

# pour ajouter du texte sur le graphique
ga <- ggdraw(add_sub(ga, label = "L. Scholtus", y = 20.8, x = 0.95, hjust = 0, 
                     size = 5))

# pour exporter le graphique en png et en définir les dimensions
ggsave(paste0("03Figures/periode/P1.png"), plot = ga, width = 22, height = 16, 
       units = "cm")

La librairie ggplot2

Spatialisation des données

Il est aussi possible d’utiliser R pour générer des cartes pour la visualiation des données, mais aussi dans le cadre d’analyses spatiales. Il existe plusieurs types d’objets spatiaux :

Objet Explications Fonction
Spatial Points Contient des points spatialisés SpatialPoints()
Spatial Points Data Frame contient des points spatialisés et les données correspondantes SpatialPointsDataFrame()
Spatial Grid Data Frame Grille spatialisée contenant les données correspondantes (comme l’élévation) SpatialGridDataFrame()
Spatial Lines Liste de lignes spatialisées SpatialLines()
Spatial Lines Data Frame Lignes spatialisées avec des données correspondant à chaque ligne SpatialLinesDataFrame()
Raster Objet raster contenant des données(comme l’élévation) raster()

Spatialisation des données

Pour utiliser ces fonctions, il est nécessaire d’installer plusieurs librairies.

Librairie Proprietes PDF
raster Lire, écrire, manipuler, analyser et modéliser des grilles spatialisées. Peux gérer des fichiers très lourds https://cran.r-project.org/web/packages/raster/raster.pdf
rgdal Permet d’importer dans R des cartes rasterisées ou vectorisées et de les exporter https://cran.r-project.org/web/packages/rgdal/rgdal.pdf
spded Collection de fonction pour créer des objets spatiaux et pour les utiliser dans des analyses spatiales https://cran.r-project.org/web/packages/spdep/spdep.pdf
mapview Affiche un fond de carte OpenStreetMap https://cran.r-project.org/web/packages/mapview/mapview.pdf

Spatialisation des données

library(sp)
library(raster)
library(rgdal)
library(mapview)
site <- read.csv("SIT_Corpus.csv", header = TRUE, sep = ";", encoding = "UTF-8", na.strings = "NA")
str(site)
## 'data.frame':    1044 obs. of  7 variables:
##  $ ID_SITE  : chr  "SIT0001" "SIT0002" "SIT0003" "SIT0004" ...
##  $ Lieu.dit : chr  "Wederath-Belginum" "Cherain-Brisy" "Han-sur-Lesse" "Klosteräcker" ...
##  $ Commune  : chr  "Wederath" "Cherain" "Han-sur-Lesse" "Breisach-Hochstetten" ...
##  $ Region   : chr  "Rhénanie Palatinat" "Wallonie" "Wallonie" "Bade-Wurtemberg" ...
##  $ Pays     : chr  "Allemagne" "Belgique" "Belgique" "Allemagne" ...
##  $ Longitude: num  7.16 5.87 5.19 7.63 7.58 ...
##  $ Latitude : chr  "49.86499" "50.17996" "50.12674" "48.02291" ...

Spatialisation des données

La première étape consiste à transformer le jeu de données en un objet spatialisé

site_1 <- site # pour conserver données sources
# transformer les colonnes coordonnées en numérique
site$x <- as.numeric(as.character(site$Longitude))
site$y <- as.numeric(as.character(site$Latitude))
## Warning: NAs introduits lors de la conversion automatique
# enlever les lignes avec des données manquantes
site <- site[complete.cases(site$x), ]
site <- site[complete.cases(site$y), ]
# pojeter avec les données
sp::coordinates(site)=~x+y # indique où se trouve les coordonnées
sp::proj4string(site) <- sp::CRS("+proj=longlat +datum=WGS84") # système de projection
# modifier le système de projection
site <- sp::spTransform(site, sp::CRS("+proj=utm +zone=32 +datum=WGS84 +units=m"))

class(site)
## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"
str(site)
## Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
##   ..@ data       :'data.frame':  1043 obs. of  7 variables:
##   .. ..$ ID_SITE  : chr [1:1043] "SIT0001" "SIT0002" "SIT0003" "SIT0004" ...
##   .. ..$ Lieu.dit : chr [1:1043] "Wederath-Belginum" "Cherain-Brisy" "Han-sur-Lesse" "Klosteräcker" ...
##   .. ..$ Commune  : chr [1:1043] "Wederath" "Cherain" "Han-sur-Lesse" "Breisach-Hochstetten" ...
##   .. ..$ Region   : chr [1:1043] "Rhénanie Palatinat" "Wallonie" "Wallonie" "Bade-Wurtemberg" ...
##   .. ..$ Pays     : chr [1:1043] "Allemagne" "Belgique" "Belgique" "Allemagne" ...
##   .. ..$ Longitude: num [1:1043] 7.16 5.87 5.19 7.63 7.58 ...
##   .. ..$ Latitude : chr [1:1043] "49.86499" "50.17996" "50.12674" "48.02291" ...
##   ..@ coords.nrs : num(0) 
##   ..@ coords     : num [1:1043, 1:2] 368087 276199 227433 397550 394110 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:1043] "1" "2" "3" "4" ...
##   .. .. ..$ : chr [1:2] "x" "y"
##   ..@ bbox       : num [1:2, 1:2] 200397 5123841 530637 5588741
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:2] "x" "y"
##   .. .. ..$ : chr [1:2] "min" "max"
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. ..@ projargs: chr "+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs"
plot(site)

Spatialisation des données

Pour vérifier que la spatialisation n’a pas d’erreur, on affiche un fond de carte sous la couche de points.

mapview(site)

Spatialisation des données

Il est aussi possible d’importer un fond de carte existant et de le mettre en forme.

srtm <- readGDAL("fdc_vrai.tif") # pour ouvrir le fond de carte et le stocker dans un objet
## fdc_vrai.tif has GDAL driver GTiff 
## and has 12899 rows and 12115 columns
class(srtm)
## [1] "SpatialGridDataFrame"
## attr(,"package")
## [1] "sp"
file_limite <- "LimiteG"
limite <- function(file_limite){
  spdf_peuple <- readOGR("shape",file_limite)
  spdf_peuple <- spTransform(spdf_peuple, CRS("+proj=utm +zone=32 +datum=WGS84 +units=m"))
  return(spdf_peuple)}

carte <- function(srtm, spdf_peuple){
  bb <- spdf_peuple@bbox
  ext <- extent(spdf_peuple@bbox)
  ext <- bb
  raster_srtm <- raster(srtm)
  sgdf_z <- crop(raster_srtm, ext, snap='near', overwrite=TRUE)
  sgdf_limite <- as(sgdf_z, 'SpatialGridDataFrame')
  return(sgdf_limite)
}
spdf_peuple <- limite(file_limite)
## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\Lizzie\Documents\Archéologie\Cours donnes\UE_humanite_num\cours2\shape", layer: "LimiteG"
## with 1 features
## It has 1 fields
sgdf_limite <- carte(srtm, spdf_peuple)
remove(srtm)

plot(sgdf_limite)

Spatialisation des données

library(colorspace)
coul <- sequential_hcl(40, h=c(90,-20), c=c(40,5), l=c(99,5), cmax = 100)
# pour créer un vecteur de 40 couleurs séquentielles
plot(sgdf_limite, col = coul, what = "image") # pour appliquer les couleurs et enlever l'échelle

Spatialisation des données

Enfin, on peut ajouter les points sur le fond de carte et les mettre en forme.

library(GISTools)

{plot(sgdf_limite, col = coul, what = "image")
  map.scale(232827, 5598702,50000,"km", 2, 50)
  north.arrow(234678, 5567129, 2000, lab="N", cex.lab=0.9, tcol = "black", col = "black")
  lines(spdf_peuple, col = "#737373",lty = 4, lwd = 1.5)
  points(site, col = "red", pch = 1, lwd = 3)
  title(main = "Carte de répartition", adj = 1, line = -2)}

Analyses spatiales

L’analyse des voisins les plus proches consiste à trouver dans un ensemble de points quels sont les plus proches à partir d’un point donné.

library(spdep)
## Warning: package 'spdep' was built under R version 4.0.3
## Warning: package 'spData' was built under R version 4.0.3
anyDuplicated(site_1$Longitude)
## [1] 17
site2 <- remove.duplicates(site)
voisin <- tri2nb(coordinates(site2)) # crée une liste de voisins
coords <- as.matrix(coordinates(site2)) # extrait les coordonnées au sein d'une matrice

del <- nb2lines(voisin, coords = coords, proj4string = CRS("+proj=utm +zone=32 +datum=WGS84 +units=m"))

class(del)
## [1] "SpatialLinesDataFrame"
## attr(,"package")
## [1] "sp"
plot(del)