INTRODUCTION

On se propose d’apprendre à utiliser R pour décrire et mettre en relation des variables quantitatives. A titre d’exemple, on va essayer de mettre en relation l’attraction-répulsion des pays du Monde telle qu’elle est perçue par les géographes de l’UGI en 2008 avec le PIB/habitant des pays du Monde en 2005. L’hypothèse est qu’il existe une corrélation positive entre l’attraction des pays et leur richesse. Mais aussi que cette règle peut comporter des exceptions intéressantes :

Au cours de cette démarche on sera amené à passer en revue toute une série d’outils de statistique univariée et bivariée. Les traitements seront effectués à l’aide du logiciel R en ne mobilisant que les fonctions élémentaires du logiciel (aucun ajout de package). Le programme pourra servir de “modèle” pour des applications à d’autres données.

A. CHARGEMENT DES DONNEES ET RECODAGES

A.1 Choix du répertoire de travail

Les données se trouvent dans le dossier UGI20082. La première étape consiste à indiquer à R le chemin d’accès vers ce dossier à l’aide de l’instruction setwd(). Puis à vérifier le contenu du dossier avec list.files() :

setwd("/Users/claudegrasland1/Documents/cg/cours/M1_Mise_a_niveau/UGI2008")
list.files()
##  [1] "article_echogeo.pdf"        "carte_ASY.pdf"             
##  [3] "carte_VOL_ATTREP.pdf"       "carte_VOL.pdf"             
##  [5] "continents.pdf"             "continents.png"            
##  [7] "exo1.html"                  "exo1.pdf"                  
##  [9] "exo1.Rmd"                   "exo2.pdf"                  
## [11] "exo2.Rmd"                   "Figure_continent.pdf"      
## [13] "Questionnaire_UGI_2008.pdf" "regression_asy_idr.csv"    
## [15] "survey_UGI_2008.csv"        "survey_UGI_2008.xls"       
## [17] "tab_monde_2005.csv"         "tab_UGI_2008.csv"          
## [19] "world2015.dbf"              "world2015.prj"             
## [21] "world2015.shp"              "world2015.shx"

A.2 Chargement des données

Les données sont dans deux fichiers différents : tab_monde_2005.csv pour les données statistiques officielles et tab_UGI_2008.csv pour les données sur les pays attractifs et répulsifs issues de l’enquête. On peut les charger dans R de deux façons différentes :

  • soit à l’aide des menus déroulants File/Import Datasets / From Text ce qui est recommandé pour les débutants.
  • soit en écrivant des lignes de programmes (Cf. ci-dessous) ce qui est plus pratique si l’on veut pouvoir ensuite executer le programme sans se servir de la souris.

Les plus malins peuvent aussi utiliser la première méthode (menus déroulant) puis récupérer dans la console les lignes de code qui ont servi à charger les données et les mettre dans leur programme : pas besoin de se servir de la souris la prochaine fois !

tab_UGI_2008 <- read.csv("tab_UGI_2008.csv", sep=";")
tab_monde_2005 <- read.csv("tab_monde_2005.csv", sep=";")

Vous pouvez examiner le contenu des tableaux de deux façons différentes

  • en cliquant dessus dans la fenêtre Environnement pour qu’ils s’affichent
  • en demandant à imprimer le début ou la fin du fichier avec les instructions head() et tail() ce qui transmet les résultats dans la consiole
head(tab_UGI_2008,3)
##   ISO3        NAME POS1 POS2 POS3 POS4 POS5 NEG1 NEG2 NEG3 NEG4 NEG5
## 1  AFG Afghanistan    0    0    0    0    0    6    7    3    5    1
## 2  AGO      Angola    0    0    0    0    0    1    2    0    0    0
## 3  ALB     Albania    0    0    0    0    0    0    1    0    0    0
tail(tab_UGI_2008,3)
##     ISO3         NAME POS1 POS2 POS3 POS4 POS5 NEG1 NEG2 NEG3 NEG4 NEG5
## 190  ZAF South Africa    0    1    0    0    2    7    4    5    3    2
## 191  ZMB       Zambia    0    1    0    0    0    0    0    0    0    0
## 192  ZWE     Zimbabwe    1    0    0    0    0    5    0    2    2    0
head(tab_monde_2005,3)
##   CODE     COUNTRY            REG          SUBREG  SUP   POP   PIB
## 1  AFG Afghanistan Asia + Pacific      South Asia  652 24507  7335
## 2  AGO      Angola         Africa Southern Africa 1247 16618 28736
## 3  ALB     Albania         Europe  Central Europe   27  3111  8524
tail(tab_monde_2005,3)
##     CODE      COUNTRY    REG          SUBREG  SUP   POP    PIB
## 179  ZAF South Africa Africa Southern Africa 1215 48073 234663
## 180  ZMB       Zambia Africa Southern Africa  743 11738   6804
## 181  ZWE     Zimbabwe Africa Southern Africa  387 12475   3221

A.3 Collage des deux tableaux (jointure)

On veut créer un tableau général tab qui fusionne les deux autres fichiers et ne conserve que les pays présents dans les deux. En effet on remarque qu’il y a 181 lignes dans un fichier et 192 dans l’autre.

Le collage de deux fichiers est une opération risquée qui peut entraîner beaucoup d’erreurs lorsqu’on le fait “à la main” avec une souris dans un tableur. La solution la plus sûre est de faire une jointure en repérant une colonne commune aux deux fichiers

Pour effectuer la jointure, on va utiliser ici les codes pays en 3 caractères. Il porte des noms différents dans les deux fichiers (“CODE” et “ISO3”) mais il s’agit de la même chose. Il suffit donc d’indiquer à R quelles colonnes utiliser dans chaque tableau pour faire la jointure.

tab<-merge(tab_monde_2005,tab_UGI_2008,by.x="CODE",by.y="ISO3")
head(tab,3)
##   CODE     COUNTRY            REG          SUBREG  SUP   POP   PIB
## 1  AFG Afghanistan Asia + Pacific      South Asia  652 24507  7335
## 2  AGO      Angola         Africa Southern Africa 1247 16618 28736
## 3  ALB     Albania         Europe  Central Europe   27  3111  8524
##          NAME POS1 POS2 POS3 POS4 POS5 NEG1 NEG2 NEG3 NEG4 NEG5
## 1 Afghanistan    0    0    0    0    0    6    7    3    5    1
## 2      Angola    0    0    0    0    0    1    2    0    0    0
## 3     Albania    0    0    0    0    0    0    1    0    0    0

A.4 Types de variables

Examinons les types actuels des variables de notre tableau à l’aide de l’instruction str()

str(tab)
## 'data.frame':    181 obs. of  18 variables:
##  $ CODE   : Factor w/ 181 levels "AFG","AGO","ALB",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ COUNTRY: Factor w/ 181 levels "Afghanistan",..: 1 4 2 170 5 6 7 8 9 26 ...
##  $ REG    : Factor w/ 6 levels "Africa","Asia + Pacific",..: 2 1 3 6 4 3 2 3 3 1 ...
##  $ SUBREG : Factor w/ 21 levels "Arabian Peninsula",..: 15 18 6 1 14 8 2 20 8 7 ...
##  $ SUP    : int  652 1247 27 84 2737 28 7682 83 83 26 ...
##  $ POP    : int  24507 16618 3111 4089 38732 3065 20395 8232 8453 7378 ...
##  $ PIB    : int  7335 28736 8524 103878 176985 4948 704832 302248 10963 780 ...
##  $ NAME   : Factor w/ 192 levels "Afghanistan",..: 1 4 2 181 6 7 8 9 10 27 ...
##  $ POS1   : int  0 0 0 0 1 0 7 2 0 0 ...
##  $ POS2   : int  0 0 0 0 0 1 9 2 0 0 ...
##  $ POS3   : int  0 0 0 1 0 0 4 2 0 0 ...
##  $ POS4   : int  0 0 0 4 1 0 10 2 0 0 ...
##  $ POS5   : int  0 0 0 0 1 0 4 2 0 0 ...
##  $ NEG1   : int  6 1 0 2 0 0 0 0 0 0 ...
##  $ NEG2   : int  7 2 1 2 0 0 0 1 0 0 ...
##  $ NEG3   : int  3 0 0 0 0 0 1 0 0 0 ...
##  $ NEG4   : int  5 0 0 2 1 0 0 0 0 0 ...
##  $ NEG5   : int  1 0 0 1 0 0 0 0 0 0 ...

A.3 Résumé général du tableau

Pour résumer notre tableau, on va utiliser la fonction summary() qui donne un résumé rapide de chacune des variables :

summary(tab)
##       CODE            COUNTRY                           REG    
##  AFG    :  1   Afghanistan:  1   Africa                   :53  
##  AGO    :  1   Albania    :  1   Asia + Pacific           :39  
##  ALB    :  1   Algeria    :  1   Europe                   :45  
##  ARE    :  1   Angola     :  1   Latin America + Caribbean:30  
##  ARG    :  1   Argentina  :  1   North America            : 2  
##  ARM    :  1   Armenia    :  1   West Asia                :12  
##  (Other):175   (Other)    :175                                 
##              SUBREG        SUP                POP         
##  Western Europe :20   Min.   :    0.00   Min.   :    102  
##  Central Europe :18   1st Qu.:   30.75   1st Qu.:   2618  
##  Western Africa :16   Median :  165.50   Median :   7868  
##  South America  :12   Mean   :  723.97   Mean   :  35934  
##  South East Asia:11   3rd Qu.:  571.50   3rd Qu.:  23613  
##  Southern Africa:11   Max.   :16378.00   Max.   :1312250  
##  (Other)        :93   NA's   :1                           
##       PIB                    NAME          POS1             POS2        
##  Min.   :      68   Afghanistan:  1   Min.   : 0.000   Min.   : 0.0000  
##  1st Qu.:    3935   Albania    :  1   1st Qu.: 0.000   1st Qu.: 0.0000  
##  Median :   13926   Algeria    :  1   Median : 0.000   Median : 0.0000  
##  Mean   :  243941   Angola     :  1   Mean   : 1.006   Mean   : 0.9724  
##  3rd Qu.:  102474   Argentina  :  1   3rd Qu.: 0.000   3rd Qu.: 0.0000  
##  Max.   :12506277   Armenia    :  1   Max.   :42.000   Max.   :27.0000  
##                     (Other)    :175                                     
##       POS3              POS4              POS5              NEG1        
##  Min.   : 0.0000   Min.   : 0.0000   Min.   : 0.0000   Min.   : 0.0000  
##  1st Qu.: 0.0000   1st Qu.: 0.0000   1st Qu.: 0.0000   1st Qu.: 0.0000  
##  Median : 0.0000   Median : 0.0000   Median : 0.0000   Median : 0.0000  
##  Mean   : 0.8785   Mean   : 0.7459   Mean   : 0.5856   Mean   : 0.8564  
##  3rd Qu.: 0.0000   3rd Qu.: 0.0000   3rd Qu.: 0.0000   3rd Qu.: 1.0000  
##  Max.   :16.0000   Max.   :16.0000   Max.   :11.0000   Max.   :21.0000  
##                                                                         
##       NEG2              NEG3             NEG4            NEG5        
##  Min.   : 0.0000   Min.   :0.0000   Min.   :0.000   Min.   : 0.0000  
##  1st Qu.: 0.0000   1st Qu.:0.0000   1st Qu.:0.000   1st Qu.: 0.0000  
##  Median : 0.0000   Median :0.0000   Median :0.000   Median : 0.0000  
##  Mean   : 0.7901   Mean   :0.6409   Mean   :0.558   Mean   : 0.4475  
##  3rd Qu.: 1.0000   3rd Qu.:1.0000   3rd Qu.:1.000   3rd Qu.: 1.0000  
##  Max.   :11.0000   Max.   :9.0000   Max.   :8.000   Max.   :10.0000  
##

B. STATISTIQUE UNIVARIEE : RESUMER LA VARIABLE INDEPENDANTE (PIB/hab)

On va décrire dans cette sections quelques fonctions simples permettant de résumer une variable quantitative et de la représenter sous forme graphique. Pour pouvoir facilement refaire les calculs ensuite on commence par créer une variable X qui correspond à la variable quantitative de notre choix :

X<-tab$PIB*1000/tab$POP

B.1 calcul des paramètres principaux

Des petites fonctions simples permettent de calculer les paramètres statistiques utiles pour résumer une distribution :

# moyenne
mean(X)
## [1] 8126.963
# ecart-type
sd(X)
## [1] 13037.71
# minimum
min(X)
## [1] 105.7197
# maximum 
max(X)
## [1] 64108.52
# mediane
median(X)
## [1] 2451.957

Il existe une fonction sépcifique pour les quantiles, permettant de spécifier les fréquences cumulées que l’on souhaite.

#premier quartile (Q1)
Q1<-quantile(X,0.25)
Q1
##      25% 
## 653.3628
#Deuxième quartile (Q2 = Mediane)
Q2<-quantile(X,0.50)
Q2
##      50% 
## 2451.957
#troisième quartile (Q3)
Q3<-quantile(X,0.75)
Q3
##      75% 
## 7319.379
#Intervalle interquartile (Q3-Q1)
Q3-Q1
##      75% 
## 6666.017
# Exemple de calcul des déciles (D1..D9)
myfreq<-c(0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1)
quantile(X,myfreq)
##         0%        10%        20%        30%        40%        50% 
##   105.7197   299.3022   541.3066   871.5084  1331.8691  2451.9573 
##        60%        70%        80%        90%       100% 
##  3262.0773  5171.2939 11636.3636 33190.0000 64108.5221

B.2 Réalisation d’un histogramme

Solution la plus simple …

# histogramme simple
hist(X)

Le même avec quelques paramètres en plus :

hist(X,
     main="PIB par habitant des pays du Monde (2005)",
     ylab = "Nb. de pays",
     xlab= "en $/hab",
     col="gray",
     sub="source : UNEP",
     breaks = 12)

B.3 Réalisation d’une boîte à moustache

La boîte à moustache ou boxplot est pratique pour étudier la forme d’une distribution et repérer l’existence de valeurs expceptionnelles (outliers).

boxplot(X)

La même en version horizontale et un peu améliorée :

boxplot(X,
     main="PIB par habitant des pays du Monde (2005)",
     horizontal = T,
     xlab = "en $/habitant",
     col="gray",
     sub="source : UNEP")

C. STATISTIQUE UNIVARIEE : RESUMER LA VARIABLE DEPENDANTE (ASY)

C.1 Création de nouvelles variables

On souhaite résumer l’avis des géographes à l’aide d’un indicateur synthétique. Pour cela on va calculer d’abord l’attraction d’un pays (ATT = somme des avis positifs) et sa répulsion (REP = somme des avis négatifs). Puis on calcule le volume total de réponses (VOL = ATT+REP) et le solde des avis positifs moins les avis négatifs (SOL=ATT-REP). Et enfin on va calculer notre indice d’attraction-répulsion en divisant le solde par le volume ATTREP=SOL/VOL).

Cela donne une valeur comprise entre -1 (tous les avis sont négatifs) et +1 (tous les avis sont positifs) avec comme point de référence 0 (autant d’avis positifs que négatifs). Même si le calcul paraît compliqué, cela donne une mesure qui est facile à interpréter et permet de repérer tout de suite les pays les plus attirants ou les plus repoussants aux yeux des 191 géographes qui ont répondu à l’enquête.

On doit toutefois tenir compte du fait que l’indice n’a pas véritablement de sens pour les pays peu cités. Aussi on fixe un seuil d’au moins 8 réponses pour choisir les pays qui seront analysés.

tab$ATT<-tab$POS1+tab$POS2+tab$POS3+tab$POS4+tab$POS5
tab$REP<-tab$NEG1+tab$NEG2+tab$NEG3+tab$NEG4+tab$NEG5
tab$VOL<-tab$ATT+tab$REP
tab$SOL<-tab$ATT-tab$REP
tab$ATTREP<-tab$SOL/tab$VOL
tab2<-tab[tab$VOL>=8,]

C.2 Résumé de la variable

On reprend le même programme que pour le PIB/hab. , mais en changeant juste le début du programme où l’on définit la variable X puis en modifiant les titres et sous-titres des graphiques

X<-tab2$ATTREP

# moyenne
mean(X)
## [1] 0.01950831
# ecart-type
sd(X)
## [1] 0.778677
# minimum
min(X)
## [1] -1
# maximum 
max(X)
## [1] 1
# mediane
median(X)
## [1] 0.2039474
#premier quartile (Q1)
Q1<-quantile(X,0.25)
Q1
##        25% 
## -0.8346154
#Deuxième quartile (Q2 = Mediane)
Q2<-quantile(X,0.50)
Q2
##       50% 
## 0.2039474
#troisième quartile (Q3)
Q3<-quantile(X,0.75)
Q3
##       75% 
## 0.7745098
#Intervalle interquartile (Q3-Q1)
Q3-Q1
##      75% 
## 1.609125
# Exemple de calcul des déciles (D1..D9)
myfreq<-c(0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1)
quantile(X,myfreq)
##         0%        10%        20%        30%        40%        50% 
## -1.0000000 -0.9945946 -0.9110000 -0.7250000 -0.3888889  0.2039474 
##        60%        70%        80%        90%       100% 
##  0.5413333  0.7326923  0.8303030  0.9269565  1.0000000
hist(X,
     main="Indice d'attraction-répulsion des pays du Monde (2008)",
     ylab = "Nb. de pays",
     xlab= "Asymétrie = (Att-Rep)/(Att+Rep)",
     col="gray",
     sub="source : Enquête UGI 2008",
     breaks = 10)

On voit apparaître une superbe distribution bimodale avec beaucoup de pays très répulsifs (proches de -1) et beaucoup de pays très répulsifs (proches de +1) mais très peu de pays mélangeant les avis positifs et négatifs à parts égales (proches de 0).

D. STATISTIQUE BIVARIEE : CORRELATION ET REGRESSION ENTRE DEUX VARIABLES

On commence par choisir deux variables dont l’une sera la variable explicative (= variable indépendante) et l’autre la variable à explique (= variable dépendante). On leur donne respectivement les noms X et Y pour pouvoir facilement utiliser le programme sur d’autres variables par la suite. On crée un petit tableau mytab comportant ces deux variables ainsi que leur code et leur nom :

CODE<-tab2$CODE
NOM<-tab2$COUNTRY
X<-tab2$PIB*1000/tab2$POP
Y<-tab2$ATTREP
mytab<-data.frame(CODE,NOM,X,Y, stringsAsFactors = F)

D.1 Création d’un graphique croissant X et Y

On utilise la fonction plot(X,Y) pour tracer le graphique

plot(X,Y)

On améliore le graphique en ajoutant des titres ainsi que le nom des gouvernorats avec la fondtion text(X,Y,NOM)

plot(X,Y,
     main = "Relation entre richesse et attraction-répulsion des pays ",
     xlab = "Indice d'attraction-répulsion (2008)",
     ylab = "PIB/habitant (2005)",
     sub = "source : Nations Unies 2005 & Enquête UGI 2008")
     text(X,Y,NOM, cex=0.6,pos=4,col="blue")
     abline(h=mean(Y), lty=2)
     abline(v=mean(X),lty=2)

D.2 Mesure du coefficient de corrélation et test de significativité

On peut calculer le coefficient de corrélation linéaire de bravais-Pearson avec la fonction cor() :

cor(X,Y)
## [1] 0.7747704

Mais il est préférable d’utiliser la fonction cor.test() pour obtenir également la significativité de la relation en fonction du nombre de degrés de liberté (nombre d’éléménts moins nombre de paramètres connus) :

cor.test(X,Y)
## 
##  Pearson's product-moment correlation
## 
## data:  X and Y
## t = 7.7503, df = 40, p-value = 1.729e-09
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.6158598 0.8731055
## sample estimates:
##       cor 
## 0.7747704

On peut voir ici que la corrélation de + 0.82 pour 22 degrés de liberté (df=22) est très significative avec une p-value égale à 1.203e-06 = 0.00000123. On peut affirmer qu’il y a moins d’une chance sur 1000 (p< 0.001) que la relation observée soit l’effet du hasard. Il y a bien un lien entre niveau de développement (IDR2011) et perception positive ou négative (ASY) des gouvernorats.

D.3 Calcul de l’équation de la droite de régression

On utilise la fonction lm() (linear model) en stockant tous les résultats dans un objet dont on choisit le nom. Par exemple, toto. Puis on examine les paramètres du modèle à l’aide de l’instruction summary()

toto<-lm(Y~X)
summary(toto)
## 
## Call:
## lm(formula = Y ~ X)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.87721 -0.38318 -0.06988  0.31726  1.34445 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -6.381e-01  1.145e-01  -5.572 1.89e-06 ***
## X            3.277e-05  4.229e-06   7.750 1.73e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4984 on 40 degrees of freedom
## Multiple R-squared:  0.6003, Adjusted R-squared:  0.5903 
## F-statistic: 60.07 on 1 and 40 DF,  p-value: 1.729e-09

On voit que l’équation de la droite de régression est ici (Y = 0.000033*X - 0.63) avec un coefficient de détermination de 60% (Multiple R-Squared). Le modèle explique donc environ 2/3 de la variance des niveaux d’attraction-répulsion. Il montre également que chaque fois que le PIB augmente de 10000 $/hab., l’attractivité augmente de 0.33

on peut refaire le graphique précédent en lui ajoutant la droite de régression avec l’instruction abline()

plot(X,Y,
     main = "Relation entre dévelopement et perception des gouvernorats",
     xlab = "Indicateur de développement régional (2011)",
     ylab = "Indice d'attraction-répulsion (2017)",
     sub = "source = MDRP & enquête PICS sur 44 étudiants du NO de Tunisie")

text(X,Y,NOM, cex=0.6,pos=4,col="blue")

     abline(h=mean(Y), lty=2)
     abline(v=mean(X),lty=2)
     abline(toto,col="red",lwd=2)

On voit que certains points sont juste sur la droite ce qui correspond à un estimation parfaite du modèle (ex. Ireland). D’autres points sont au-dessus de la droite et sont plus attrractifs que prévus par le modèle (ex. Syrie; Maroc, Grèce, Liban, Finlande, …). D’autres enfin sont situés sous la droite et sont moins attrractifs que prévis (ex. israël, Arabie Saoudite, USA, Danemark, …). Ces résidus positifs ou négatifs correspondent à des facteurs d’attraction ou de répulsion autres que le PIB/habitant. Il est intéressant de les calculer et de les stocker dans un tableau pour chercher des explications à ces écarts au modèle de prévision.

D.4 Calcul des valeurs estimées et des résidus

Il est intéressant de stocker les valeurs estimées et les résidus pour pouvoir ensuite les analyser ou les cartographier. On va donc les extraire du modèle et les ajouter au petit tableau de données mytab. Puis on trie celui-ci par ordre croissant de résidus.

mytab$Y_ESTIM<-toto$fitted.values
mytab$Y_RESID<-toto$residuals
mytab<-mytab[order(mytab$Y_RESID),]

les 5 plus forts résidus négatifs

head(mytab,5)
##    CODE                NOM         X          Y    Y_ESTIM    Y_RESID
## 23  ISR             Israel 18163.628 -0.9200000 -0.0427902 -0.8772098
## 11  DNK            Denmark 46925.041  0.2500000  0.8998653 -0.6498653
## 35  SAU       Saudi Arabia 13130.479 -0.8461538 -0.2077517 -0.6384021
## 39  USA                USA 40727.109  0.1304348  0.6967280 -0.5662933
## 34  RUS Russian Federation  5203.492 -0.9459459 -0.4675587 -0.4783872

Le plus fort résidu négatif concerne Israël. Au vu de sa richesse (18163 $/hab), on s’attend à ce que ce pays ait un niveau d’attraction-répulsion moyen (-0.05). Mais l’indice observé est clairement négatif (-0.92). Donc on conclue que le pays est moins attractif que ce que laisserait prévoir sa richesse (-0.87). Il en va de même pour le Danemark, l’Arabie Saoudite, les USA ou la Russie …

les 10 plus forts résidus négatifs

tail(mytab,5)
##    CODE     NOM         X         Y     Y_ESTIM   Y_RESID
## 14  ESP   Spain 25699.489 0.8596491  0.20419770 0.6554514
## 26  LBN Lebanon  5249.143 0.2727273 -0.46606256 0.7387898
## 18  GRC  Greece 19979.573 0.7647059  0.01672742 0.7479785
## 28  MAR Morocco  1682.243 0.4666667 -0.58296772 1.0496344
## 37  SYR   Syria  1331.869 0.7500000 -0.59445123 1.3444512

Le plus fort résidu négatif concerne … la Syrie ! Nous sommes en 2008, avant le début de la crise économique, des Printemps Arabes et de la guerre de Syrie. A l’époque, beaucoup de géographes tunisiens présents à l’UGI ont une opinion très positive de ce pays qui apparaît stable, laïque et chargé d’histoire. Au vu de sa richesse (1331 $/hab), on s’attend à ce que ce pays ait un niveau d’attraction-répulsion clairement négatif (-0.59). Mais l’indice observé est franchement positif (+0.75). Donc on conclue que le pays est - en 2008 - beaucoup plus attractif que ce que laisserait prévoir sa richesse (+1.34). Il en va de même pour le Liban, le Grèce, le Maroc ou l’Espagne.

On exporte le tableau pour pouvoir l’étudier, le cartographier ou l’insérer dans un texte.

write.table(mytab,
            "regression_asy_idr.csv",
            sep=";",
            dec=",",
            row.names = F)

E. BONUS CARTOGRAPHIQUE …

On peut facilement réaliser des cartes à l’aide des fichiers exportés par R, par exemple en utilisant MAGRIT (Cf. cours de mise à niveau en cartographie). Mais les étudiants qui prennent l’habitude d’utiliser R pourront aussi essayer de faire leurs cartes directement dans R en utilisant le package cartography(). Deux exemples rapides :

E.1 Chargement du fonds de carte

Nous souhaitons charger un fonds de carte au format shapefile (.shp) qui se trouve dans le dossier. On commence par examiner la liste des fonds de carte présentes dans ce dossier :

list.files()
##  [1] "article_echogeo.pdf"        "carte_ASY.pdf"             
##  [3] "carte_VOL_ATTREP.pdf"       "carte_VOL.pdf"             
##  [5] "continents.pdf"             "continents.png"            
##  [7] "exo1.html"                  "exo1.pdf"                  
##  [9] "exo1.Rmd"                   "exo2_files"                
## [11] "exo2.pdf"                   "exo2.Rmd"                  
## [13] "Figure_continent.pdf"       "Questionnaire_UGI_2008.pdf"
## [15] "regression_asy_idr.csv"     "survey_UGI_2008.csv"       
## [17] "survey_UGI_2008.xls"        "tab_monde_2005.csv"        
## [19] "tab_UGI_2008.csv"           "world2015.dbf"             
## [21] "world2015.prj"              "world2015.shp"             
## [23] "world2015.shx"

Le fonds de carte qui nous intéresse est le fonds de carte world2015.shp. Pour le charger dans R, nous avons besoin d’utiliser un package spécifique. Nous allons choisir ici le package sf qui a remplacé sp

library(sf)
## Linking to GEOS 3.6.1, GDAL 2.1.3, PROJ 4.9.3
map<-st_read("world2015.shp",stringsAsFactors = F)
## Reading layer `world2015' from data source `/Users/claudegrasland1/Documents/cg/cours/M1_Mise_a_niveau/UGI2008/world2015.shp' using driver `ESRI Shapefile'
## Simple feature collection with 192 features and 4 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -179.9585 ymin: -55.99725 xmax: 179.9959 ymax: 83.61347
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs

On peut examiner le contenu de cet objet avec str()

str(map)
## Classes 'sf' and 'data.frame':   192 obs. of  5 variables:
##  $ ISO2    : chr  "AE" "AF" "AG" "AL" ...
##  $ ISO3    : chr  "ARE" "AFG" "ATG" "ALB" ...
##  $ name    : chr  "United Arab Emirates" "Afghanistan" "Antigua and Barbuda" "Albania" ...
##  $ ISO3_SO : chr  "ARE" "AFG" "ATG" "ALB" ...
##  $ geometry:sfc_MULTIPOLYGON of length 192; first list element: List of 1
##   ..$ :List of 1
##   .. ..$ : num [1:12, 1:2] 56.3 56.4 55.8 55.8 55.2 ...
##   ..- attr(*, "class")= chr  "XY" "MULTIPOLYGON" "sfg"
##  - attr(*, "sf_column")= chr "geometry"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA
##   ..- attr(*, "names")= chr  "ISO2" "ISO3" "name" "ISO3_SO"

On voit qu’il y a une liste de colonnes comme dans un data.frame mais aussi une dernière colonne geometry qui contient la liste des points composant chaque polygone

On peut visualiser facilement ce fonds de carte avec l’instruction plot() appliquée à la colonne contenant la variable geometry

plot(map$geometry)

E.2 Jointure entre données et géométries

Il s’agit d’une jointure à l’aide de la fonction merge() en précisant le nom des colonnes servant à la jointure dans le tableau de données et dans la géométrie :

map_don<-merge(map, tab,by.x="ISO3",by.y="CODE")
str(map_don)
## Classes 'sf' and 'data.frame':   181 obs. of  27 variables:
##  $ ISO3    : chr  "AFG" "AGO" "ALB" "ARE" ...
##  $ ISO2    : chr  "AF" "AO" "AL" "AE" ...
##  $ name    : chr  "Afghanistan" "Angola" "Albania" "United Arab Emirates" ...
##  $ ISO3_SO : chr  "AFG" "AGO" "ALB" "ARE" ...
##  $ COUNTRY : Factor w/ 181 levels "Afghanistan",..: 1 4 2 170 5 6 7 8 9 26 ...
##  $ REG     : Factor w/ 6 levels "Africa","Asia + Pacific",..: 2 1 3 6 4 3 2 3 3 1 ...
##  $ SUBREG  : Factor w/ 21 levels "Arabian Peninsula",..: 15 18 6 1 14 8 2 20 8 7 ...
##  $ SUP     : int  652 1247 27 84 2737 28 7682 83 83 26 ...
##  $ POP     : int  24507 16618 3111 4089 38732 3065 20395 8232 8453 7378 ...
##  $ PIB     : int  7335 28736 8524 103878 176985 4948 704832 302248 10963 780 ...
##  $ NAME    : Factor w/ 192 levels "Afghanistan",..: 1 4 2 181 6 7 8 9 10 27 ...
##  $ POS1    : int  0 0 0 0 1 0 7 2 0 0 ...
##  $ POS2    : int  0 0 0 0 0 1 9 2 0 0 ...
##  $ POS3    : int  0 0 0 1 0 0 4 2 0 0 ...
##  $ POS4    : int  0 0 0 4 1 0 10 2 0 0 ...
##  $ POS5    : int  0 0 0 0 1 0 4 2 0 0 ...
##  $ NEG1    : int  6 1 0 2 0 0 0 0 0 0 ...
##  $ NEG2    : int  7 2 1 2 0 0 0 1 0 0 ...
##  $ NEG3    : int  3 0 0 0 0 0 1 0 0 0 ...
##  $ NEG4    : int  5 0 0 2 1 0 0 0 0 0 ...
##  $ NEG5    : int  1 0 0 1 0 0 0 0 0 0 ...
##  $ ATT     : int  0 0 0 5 3 1 34 10 0 0 ...
##  $ REP     : int  22 3 1 7 1 0 1 1 0 0 ...
##  $ VOL     : int  22 3 1 12 4 1 35 11 0 0 ...
##  $ SOL     : int  -22 -3 -1 -2 2 1 33 9 0 0 ...
##  $ ATTREP  : num  -1 -1 -1 -0.167 0.5 ...
##  $ geometry:sfc_MULTIPOLYGON of length 181; first list element: List of 1
##   ..$ :List of 1
##   .. ..$ : num [1:36, 1:2] 74.9 74.5 72.6 71.2 71.6 ...
##   ..- attr(*, "class")= chr  "XY" "MULTIPOLYGON" "sfg"
##  - attr(*, "sf_column")= chr "geometry"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
##   ..- attr(*, "names")= chr  "ISO3" "ISO2" "name" "ISO3_SO" ...

Nous allons maintenant apprendre à réaliser quelques types simples de carte avec le package cartography développé par les ingénieurs de l’UMS RIATE. On charge le package avec la dernière version :

library(cartography)
library(lwgeom)
## Linking to liblwgeom 2.5.0dev r16016, GEOS 3.6.1, PROJ 4.9.3

E.3 Réalisation d’une carte de stock

On se propose de cartographier la variable VOL en nous servant de la fonction PropSymbolsLayer(). La seule difficulté est de bien repérer quelle est la colonne qui servira à faire la jointure dans le fichier de données (tab) et dans le fonds de carte (map).

plot(map_don$geometry)
propSymbolsLayer(map_don,
                 var = "VOL")               
## Warning in st_centroid.sfc(x = sf::st_geometry(x), of_largest_polygon
## = max(sf::st_is(sf::st_as_sf(x), : st_centroid does not give correct
## centroids for longitude/latitude data

Pour faire une carte plus jolie, il faut ajuster la taille des cercles avec l’instruction inches = qui donne la taille maximum de la carte. Et ajouter différents éléments de légende, comme dans l’exemple ci-dessous :

par(mar=c(2,2,2,2))
plot(map_don$geometry,col="lightyellow",lwd=0.2)
propSymbolsLayer(map_don,
                 var = "VOL",
                 inches=0.10,
                 legend.title.txt = "nb. de citations",
                 legend.pos = "topleft",
                 legend.style = "e",
                 legend.values.rnd = -2,
                 add=T)
## Warning in st_centroid.sfc(x = sf::st_geometry(x), of_largest_polygon
## = max(sf::st_is(sf::st_as_sf(x), : st_centroid does not give correct
## centroids for longitude/latitude data
layoutLayer(title = "PAYS LES PLUS CITES (positivement ou négativement)",
            col = "grey",
            scale=NULL,
            north = T,
            author = "C. Grasland, U. Paris Diderot, 2018",
            sources = "Source : Enquête UGI 2008")

Si on veut sauvegarder la carte au format vectoriel, on execute le même programme en ajoutant deux instructions en début et fin de programme :

pdf("carte_VOL.pdf",width = 8,height = 5)
par(mar=c(2,2,2,2))
plot(map_don$geometry,col="lightyellow",lwd=0.2)
propSymbolsLayer(map_don,
                 var = "VOL",
                 inches=0.10,
                 legend.title.txt = "nb. de citations",
                 legend.pos = "topleft",
                 legend.style = "e",
                 legend.values.rnd = -2)
## Warning in st_centroid.sfc(x = sf::st_geometry(x), of_largest_polygon
## = max(sf::st_is(sf::st_as_sf(x), : st_centroid does not give correct
## centroids for longitude/latitude data
layoutLayer(title = "PAYS LES PLUS CITES (positivement ou négativement)",
            col = "grey",
            scale=NULL,
            north = T,
            author = "C. Grasland, U. Paris Diderot, 2018",
            sources = "Source : Enquête UGI 2008")
dev.off()
## quartz_off_screen 
##                 2

E.4 Réalisation d’une carte de taux (choroplèthe)

On se propose de cartographier la variable ATTREP (indice d’attraction-répulsion) en nous servant de la fonction ChoroLayer(). On commence par appliquer les paramètres par défaut :

choroLayer(map_don,var = "ATTREP")

On va maintenant améliorer la carte en ajoutant toute une série de paramètres et en superposant un habillage composé des limites de gouvernorats.

par(mar=c(2,2,2,2))
choroLayer(map_don,
                 var = "ATTREP",
           method = "equal",
           nclass = 8,
           col = carto.pal(pal1="orange.pal", n1=4, pal2 = "green.pal", n2 = 4),
           border = "grey40",
           lwd=0.4,
           add = FALSE,
           legend.pos = "bottomleft",
           legend.title.txt = "(POS-NEG)/(POS+NEG)",
           legend.values.rnd = 2)

layoutLayer(title = "ATTRACTION-REPULSION DES PAYS DU MONDE",
            col = "grey",
            scale=NULL,
            north = T,
            author = "C. Grasland, U. Paris Diderot, 2018",
            sources = "Source : Enquête UGI 2008")

Comme précédemment on peut exporter la figure au format pdf pour la retoucher ensuite. On peut même spécifier la longueur et la largeur de l’image que l’on souhaite produire :

pdf("carte_ASY.pdf",width = 8,height = 5)
par(mar=c(2,2,2,2))
choroLayer(map_don,
                 var = "ATTREP",
           method = "equal",
           nclass = 8,
           col = carto.pal(pal1="orange.pal", n1=4, pal2 = "green.pal", n2 = 4),
           border = "grey40",
           lwd=0.4,
           add = FALSE,
           legend.pos = "bottomleft",
           legend.title.txt = "(POS-NEG)/(POS+NEG)",
           legend.values.rnd = 2)

par(mar=c(2,2,2,2))
choroLayer(map_don,
                 var = "ATTREP",
           method = "equal",
           nclass = 8,
           col = carto.pal(pal1="orange.pal", n1=4, pal2 = "green.pal", n2 = 4),
           border = "grey40",
           lwd=0.4,
           add = FALSE,
           legend.pos = "bottomleft",
           legend.title.txt = "(POS-NEG)/(POS+NEG)",
           legend.values.rnd = 2)

layoutLayer(title = "ATTRACTION-REPULSION DES PAYS DU MONDE",
            col = "grey",
            scale=NULL,
            north = T,
            author = "C. Grasland, U. Paris Diderot, 2018",
            sources = "Source : Enquête UGI 2008")
dev.off()
## quartz_off_screen 
##                 2

E.5 Réalisation d’une carte stock + taux

On va utiliser le programme propSymbolsChoroLayer. Le programme la plus simple est :

plot(map_don$geometry)
propSymbolsChoroLayer(map_don,var="VOL",var2 = "ATTREP")
## Warning in st_centroid.sfc(x = sf::st_geometry(x), of_largest_polygon
## = max(sf::st_is(sf::st_as_sf(x), : st_centroid does not give correct
## centroids for longitude/latitude data

Mais on peut évidemment l’améliorer :

par(mar=c(2,2,2,2))
plot(map_don$geometry,col="lightgray",lwd=0.2)
propSymbolsChoroLayer(map_don,
           var = "VOL",
           inches = 0.08,
           var2 = "ATTREP",
           method = "equal",
           nclass = 8,
           col = carto.pal(pal1="orange.pal", n1=4, pal2 = "green.pal", n2 = 4),
           border = "grey40",
           lwd=0.4,
           legend.var.pos="bottomleft",
           legend.var.values.rnd = 0,
           legend.var.title.txt = "Nb. réponses", 
           legend.var.style = "e",
           legend.var2.pos = "topleft",
           legend.var2.title.txt = "Attraction-Repulsion",
           legend.var2.values.rnd = 2)
## Warning in st_centroid.sfc(x = sf::st_geometry(x), of_largest_polygon
## = max(sf::st_is(sf::st_as_sf(x), : st_centroid does not give correct
## centroids for longitude/latitude data
layoutLayer(title = "CONNAISSANCE ET ATTRACTION-REPULSION DES PAYS DU MONDE",
            col = "grey",
            scale = NULL,
            north = T,
            author = "C. Grasland, U. Paris Diderot, 2018",
            sources = "Source : Enquête UGI 2008")

?layoutLayer()

Comme d’habitude, on exporte le résultat en pdf :

pdf("carte_VOL_ATTREP.pdf",width = 8,height = 5)
par(mar=c(2,2,2,2))
plot(map_don$geometry,col="lightgray",lwd=0.2)
propSymbolsChoroLayer(map_don,
           var = "VOL",
           inches = 0.08,
           var2 = "ATTREP",
           method = "equal",
           nclass = 8,
           col = carto.pal(pal1="orange.pal", n1=4, pal2 = "green.pal", n2 = 4),
           border = "grey40",
           lwd=0.4,
           legend.var.pos="bottomleft",
           legend.var.values.rnd = 0,
           legend.var.title.txt = "Nb. réponses", 
           legend.var.style = "e",
           legend.var2.pos = "topleft",
           legend.var2.title.txt = "Attraction-Repulsion",
           legend.var2.values.rnd = 2)
## Warning in st_centroid.sfc(x = sf::st_geometry(x), of_largest_polygon
## = max(sf::st_is(sf::st_as_sf(x), : st_centroid does not give correct
## centroids for longitude/latitude data
layoutLayer(title = "CONNAISSANCE ET ATTRACTION-REPULSION DES PAYS DU MONDE",
            col = "grey",
            scale=NULL,
            north = T,
            author = "C. Grasland, U. Paris Diderot, 2018",
            sources = "Source : Enquête UGI 2008")
dev.off()
## quartz_off_screen 
##                 2