[INT2] : Enseigner avec R

GEO UNIV’R Tunisie 2024

Auteur·rice

Claude Grasland

Date de publication

2024-05-13

INTRODUCTION

On se propoe de passer en revue l’ensemble du programme de la semaine à travers un exercice pédagogique associant statistique, cartographie, analyse spatiale, modélisation … Cet exercice correspond typiquement au travail qu’on pourrait donner à des étudiants de licence 2 ou 3 ayant suivi des cours de statistique, cartographie et SIG.

On limite volontairement au minimum le nombre de packages R utilisés, en se limitant ici à deux : readxl et sf. Ceux-ci sont utilisés pour importer/exporter les données statistiques et les données géographiques.

On propose de réaliser une analyse de la distribution du taux d’équipement en ordinateur dans les délégations d’un gouvernorat en 2004 et 2014. Le programme est rédigé de telle sorte qu’on puisse facilement passer d’un gouvernorat à l’autre.

Code

# Importe les données
library(readxl)
data<-read_xls(path = "data/RP_Tunisie/data/don_del.xls",
              sheet = "data")

choix_gou <- "Sfax"

A chaque étape, on peut comparer les résultats obtenus en se servant de programmes R et ceux obtenus en se servant des logiciels utilisant la souris comme Excel, Xlstat, Magrit, QGIS, ArCGIS, Philcarto , …

(1) DONNEES STATISTIQUES

Télécharger le jeu de données

RP_Tunisie

Acquisition

Consigne : Après avoir chargé le tableau excel tun_del_2004.xls, sélectionnez les délégations du gouvernorat de Sfax et construisez un tableau décrivant le nombre ménages et leur équipement en ordinateur en 2004 et 2014
Résultat :

Code

# Importe les données
library(readxl)
data<-read_xls(path = "data/RP_Tunisie/data/don_del.xls",
              sheet = "data")

# Sélectionne les lignes
don <- data[data$gou_nom==choix_gou,]

# simplifie le code
don$code<-substr(don$del_code,7,8)

# Calcule le taux d'équipement
don$equip_2004 <- 100*don$ordin_2004/don$menag_2004
don$equip_2014 <- 100*don$ordin_2014/don$menag_2014

# Sélectionne les colonnes
don <- don[,c("code","del_nom_fr","del_nom_ar","gou_cap","menag_2004","menag_2014", "ordin_2004", "ordin_2014")]

# renomme
names(don)<-c("code","nomfr","nomar","cap", "men04","men14","equ04","equ14")

# Affiche
library(knitr)
kable(don, digits=0,
      main = "Tableau 1 : Ménages et équipements")

code	nomfr	nomar	cap	men04	men14	equ04	equ14
SM	Sfax Medina	صفاقس المدينة	1	28794	29107	3772	14494
SO	Sfax Ouest	صفاقس الغربية	0	26653	30371	3385	14841
SZ	Sakiet Ezzit	ساقية الزيت	0	17541	23096	2193	10978
SD	Sakiet Eddair	ساقية الدائر	0	23633	30277	2056	12418
SS	Sfax Sud	صفاقس الجنوبية	0	24215	31866	2664	14057
TI	Tina	طينة	0	10239	15577	307	4576
AG	Agareb	عقارب	0	7689	9743	177	1543
DJ	Jebeniana	جبنيانة	0	9502	11823	171	1948
AM	EL Amra	العامرة	0	6181	7501	80	966
HA	Hencha	الحنشة	0	8831	11385	97	1087
MC	Menzel Chaker	منزل شاكر	0	7084	8387	35	477
GH	El Ghraiba	الغريبة	0	3029	3390	21	269
BA	Bir Ali Ben khlifa	بئر علي بن خليفة	0	9598	11230	58	994
SK	Skhira	الصخيرة	0	5137	6807	72	936
MA	Mahres	المحرس	0	6506	8311	247	2245
KE	Kerkenah	قرقنة	0	3933	3869	102	1314

Transformation

Consigne : Ajoutez deux colonnes décrivant le taux d’équipement en % des ménages en 2004 et 2014

Code

don$pct04<-100*don$equ04/don$men04
don$pct14<-100*don$equ14/don$men14
kable(don, 
      main = "Tableau 2 : taux d'équipement",
      digits=c(0,0,0,0,0,0,0,0,1,1))

code	nomfr	nomar	cap	men04	men14	equ04	equ14	pct04	pct14
SM	Sfax Medina	صفاقس المدينة	1	28794	29107	3772	14494	13.1	49.8
SO	Sfax Ouest	صفاقس الغربية	0	26653	30371	3385	14841	12.7	48.9
SZ	Sakiet Ezzit	ساقية الزيت	0	17541	23096	2193	10978	12.5	47.5
SD	Sakiet Eddair	ساقية الدائر	0	23633	30277	2056	12418	8.7	41.0
SS	Sfax Sud	صفاقس الجنوبية	0	24215	31866	2664	14057	11.0	44.1
TI	Tina	طينة	0	10239	15577	307	4576	3.0	29.4
AG	Agareb	عقارب	0	7689	9743	177	1543	2.3	15.8
DJ	Jebeniana	جبنيانة	0	9502	11823	171	1948	1.8	16.5
AM	EL Amra	العامرة	0	6181	7501	80	966	1.3	12.9
HA	Hencha	الحنشة	0	8831	11385	97	1087	1.1	9.5
MC	Menzel Chaker	منزل شاكر	0	7084	8387	35	477	0.5	5.7
GH	El Ghraiba	الغريبة	0	3029	3390	21	269	0.7	7.9
BA	Bir Ali Ben khlifa	بئر علي بن خليفة	0	9598	11230	58	994	0.6	8.9
SK	Skhira	الصخيرة	0	5137	6807	72	936	1.4	13.8
MA	Mahres	المحرس	0	6506	8311	247	2245	3.8	27.0
KE	Kerkenah	قرقنة	0	3933	3869	102	1314	2.6	34.0

(2) STATISTIQUE UNIVARIEE

Paramètres principaux

Consigne : Etudiez l’évolution du taux d’équipement des ménages en ordinateur par délégation en 2004 et 2014 en vous servant de paramètres principaux (valeurs centrales, paramètres de dispersion). Puis établissez deux histogrammes permettant de visualiser l’évolution.

Code

# sélectionne les variables
sel <- don[,c("pct04","pct14")]

# Tableau standard
quant<-apply(sel,2,quantile)
moy<-apply(sel,2,mean)
ect<-apply(sel,2,sd)
cv<-100*ect/moy
tab<-rbind(quant,moy,ect,cv)
row.names(tab) <-c("Minimum","Q1","Médiane","Q3","Maximum","Moyenne","Ecart-type", "C.V. (%)")

kable(tab, caption="Paramètres principaux", digits =1,
      col.names = c("Situation en 2004","Situation en 2014"),
      )

Paramètres principaux
	Situation en 2004	Situation en 2014
Minimum	0.5	5.7
Q1	1.2	12.0
Médiane	2.4	21.7
Q3	9.3	41.8
Maximum	13.1	49.8
Moyenne	4.8	25.8
Ecart-type	4.9	16.3
C.V. (%)	101.6	63.3

Histogrammes

Consigne : Etablissez deux histogrammes permettant de visualiser la forme de la distribution du taux d’équipement et son évolution entre 2004 et 2014.

Code

par(mfrow=c(2,1))

mintot <-min(c(sel$pct04, sel$pct14))
maxtot <-max(c(sel$pct04, sel$pct14))

# Histogramme
hist(sel$pct04,
     breaks=quantile(sel$pct04),
     xlim=c(mintot,maxtot),
     col="gray80",
     main= "Situation en 2004",
     xlab = "Taux d'équipement (%)",
     ylab = "Fréquence moyenne")
rug(sel$pct04, col="black", lwd=2)
lines(density(sel$pct04), lty=3,lwd=2)

hist(sel$pct14,
    breaks=quantile(sel$pct14),
     xlim=c(mintot,maxtot),
     col="gray80",
     main= "Situation en 2014",
     xlab = "Taux d'équipement (%)",
     ylab = "Fréquence moyenne")
rug(sel$pct14, col="black", lwd=2)
lines(density(sel$pct14), lty=3,lwd=2)

(3) DONNEES GEOMETRIQUES

Acquision

Consigne : Après avoir chargé le shapefile Tunisie2014_del.shp, extraire les délégations correspondant à votre gouvernorat et afficher le fonds de carte avec le code des unités.

Code

# Chargement du package spatial features
library(sf)

# Importation du fonds de carte complet
map<-st_read("data/RP_Tunisie/geom/map_del.geojson", quiet=T)

# Selection d'un gouvernorat
map <- map[map$gou_nom == choix_gou,]

# simplifie le code
map$code <- substr(map$del_code,7,8)

# ne conserve que le code,la capitale et la géométrie
map <- map[,c("code","gou_cap","geometry")]


# Affichage du fonds de carte
par(mar=c(0,0,3,0))
plot(map$geometry, 
     col="gray90",
     main = "Code des unités spatiales de la zone d'étude")

# Ajout du code des unités spatiales
coo<-st_coordinates(st_centroid(map))
text(coo, map$code, cex=0.5,col="black",)

Transformation

Consigne : ajoutez une colonne correspondant à la distance en km au chef-lieu de gouvernorat et faites en une cartographie en prenant comme bornes de classes 0, 5, 10, 20, 40, 80, 160 km.

Code

cap<-map[map$gou_cap==1,]
map$dist<-as.numeric(st_distance(st_centroid(map),st_centroid(cap)))/1000
plot(map["dist"], main="Distance au chef-lieu (en km)",
     breaks=c(0,5, 10,20,40,80, 160),
     pal=c("gray10", "gray30","gray50","gray70","gray90", "gray100"),)

(4) CARTOGRAPHIE THEMATIQUE

Cartes de stock

Consigne : Réalisez deux cartes de stock décrivant le nombre de ménages équipés en ordinateur en 2004 et 2014. Vous utiliserez la même échelle de taille pour rendre les deux cartes comparables.

Code

library(mapsf)
map<-map[,c("code","geometry")]
map_don <- merge(map, don, by="code")
maxequ<-max(don$equ04,don$equ14)

par(mfrow=c(1,2))
mf_map(map_don$geometry, col="white")
mf_map(map_don, type="prop", var="equ04",
       val_max = maxequ, inches=0.1, col="gray20", 
       leg_title = "Nb. de ménages équipés",)
mf_layout(title="2004",frame = T, credits = "Source : INS Tunisie")

mf_map(map_don$geometry, col="white")
mf_map(map_don, type="prop", var="equ14",
       val_max = maxequ, inches=0.1, col="gray20",
       leg_title = "Nb. de ménages équipés")
mf_layout(title="2014",frame = T, credits = "Source : INS Tunisie")

Cartes de ratio (choroplèthes)

Consigne : Réalisez deux cartes de taux décrivant le nombre de ménages équipés en ordinateur en 2004 et 2014. Pour les rendre comparables vous utiliserez dans chaque carte une partition en quintiles (5 classes d’effectifs égaux)

Code

library(mapsf)
map_don <- merge(map, don, by="code")
maxequ<-max(don$equ04,don$equ14)

par(mfrow=c(1,2))
mf_map(map_don, type="choro", var="pct04",
       breaks = "quantile",nbreaks = 5, pal ="Grays",
       leg_title = "% ménages équipés",leg_val_rnd = 1)
mf_layout(title="2004",frame = T, credits = "Source : INS Tunisie")

mf_map(map_don, type="choro", var="pct14",
       breaks = "quantile",nbreaks = 5, pal ="Grays",
       leg_title = "% ménages équipés",leg_val_rnd = 1)
mf_layout(title="2014",frame = T, credits = "Source : INS Tunisie")

(5) STATISTIQUES BIVARIEES

Nuage de points

Consigne : Tracez un nuage de point montrant l’évolution de l’indicateur entre les deux dates.

Code

# prépration de l'analyse
code<-don$code
nomfr<-don$nomfr
nomar<-don$nomar
X<-don$pct04
Y<-don$pct14
tab<-data.frame(code,nomfr, nomar,X,Y)

# Diagramme
plot(tab$X,tab$Y, 
     pch=20,
     cex=0.8,
     col="red",
     main = "Evolution du taux d'équipement",
     xlab="tx. equipement. 2004",
     ylab ="tx. equipement 2014")
text(tab$X,tab$Y,tab$code, 
     pos=2,
     cex=0.5,
     col="blue")

Analyse de la corrélation

Consigne : calculez les coefficients de corrélation de Pearson et Spearman et testez leur sgnificativité.

Code

cor.test(X,Y, method="pearson")


    Pearson's product-moment correlation

data:  X and Y
t = 10.264, df = 14, p-value = 6.756e-08
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.8306927 0.9791923
sample estimates:
      cor 
0.9395226

Code

cor.test(X,Y, method="spearman")


    Spearman's rank correlation rho

data:  X and Y
S = 12, p-value < 2.2e-16
alternative hypothesis: true rho is not equal to 0
sample estimates:
      rho 
0.9823529

Droite de régression

Consigne : calculez l’equation de la droite de régression et tracez- là sur le graphique.

Code

modreg <- lm(Y~X)
summary(modreg)


Call:
lm(formula = Y ~ X)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.5455 -2.6977 -1.7530  0.8751 15.1464 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  10.6836     2.0643   5.175 0.000141 ***
X             3.1357     0.3055  10.264 6.76e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 5.79 on 14 degrees of freedom
Multiple R-squared:  0.8827,    Adjusted R-squared:  0.8743 
F-statistic: 105.4 on 1 and 14 DF,  p-value: 6.756e-08

Code

plot(tab$X,tab$Y, 
     pch=20,
     cex=0.8,
     col="red",
     main = "Droite de régression",
     xlab="tx. equipement. 2004",
     ylab ="tx. equipement 2014")

abline(modreg,col="blue",lwd=1)

Analyse des résidus

Consigne : Calculez les valeurs théoriques prévus par le modèle de régression et les résidus. Affichez le tableau correspondant après l’avoir trié par ordre de résidus croissants.

Code

tab$Y_est <- modreg$fitted.values
tab$Y_res <- modreg$residuals
tab<-tab[order(tab$Y_res),]
kable(tab, digits=1)

	code	nomfr	nomar	X	Y	Y_est	Y_res
11	MC	Menzel Chaker	منزل شاكر	0.5	5.7	12.2	-6.5
12	GH	El Ghraiba	الغريبة	0.7	7.9	12.9	-4.9
10	HA	Hencha	الحنشة	1.1	9.5	14.1	-4.6
13	BA	Bir Ali Ben khlifa	بئر علي بن خليفة	0.6	8.9	12.6	-3.7
3	SZ	Sakiet Ezzit	ساقية الزيت	12.5	47.5	49.9	-2.4
7	AG	Agareb	عقارب	2.3	15.8	17.9	-2.1
1	SM	Sfax Medina	صفاقس المدينة	13.1	49.8	51.8	-2.0
9	AM	EL Amra	العامرة	1.3	12.9	14.7	-1.9
2	SO	Sfax Ouest	صفاقس الغربية	12.7	48.9	50.5	-1.6
14	SK	Skhira	الصخيرة	1.4	13.8	15.1	-1.3
5	SS	Sfax Sud	صفاقس الجنوبية	11.0	44.1	45.2	-1.1
8	DJ	Jebeniana	جبنيانة	1.8	16.5	16.3	0.1
4	SD	Sakiet Eddair	ساقية الدائر	8.7	41.0	38.0	3.1
15	MA	Mahres	المحرس	3.8	27.0	22.6	4.4
6	TI	Tina	طينة	3.0	29.4	20.1	9.3
16	KE	Kerkenah	قرقنة	2.6	34.0	18.8	15.1

Cartographie des résidus

Consigne : Cartographiez les résidus après les avoir standardisés.

Code

library(mapsf)

# Standardisation des résidus
tab$Y_res_std<-tab$Y_res/sd(tab$Y_res)

# Jointure avec la carte
map<-map[,c("code","geometry")]
map_reg <- merge(map, tab, by="code")

# Choix de la palette et des classes
library(RColorBrewer)
mypal<-brewer.pal(n = 6, name = "RdYlBu")
mybreaks = c(-10, -2,-1,0,1,2,10)

mf_map(map_reg, type="choro", var="Y_res_std",
       pal = mypal, breaks=mybreaks,
       leg_title = "Résidus standardisés",leg_val_rnd = 1)
mf_layout(title="Ecarts à la tendance 2004-2014",frame = T, credits = "Source : INS Tunisie")