6  Méthode randomisée

Cliquer ici pour télécharger la présentation.

ATTENTION : Il va de soi que les AP malgaches n’ont à aucun moment été assignées aléatoirement. Lors de cette séquence, on fait “comme si”, pour montrer la manière dont les données sont analysées quand il y a eu assignation aléatoire. On verra en fin de session les limites d’une telle approche et dans les suivantes des manières de construire des contrefactuels plus vraisemblables pour un sujet comme celui-ci.

L’analyse est effectuée partir des données préparées dans le Chapitre 3. On commence par vérifier s’il existe des déséquilibres flagrants entre les aires qui ont été protégées avant 2015 et celles qui ont été protégées en 2015, en matière de surface totale ou de part couverte par des forêts en 1996.

Code 
# On charge les librairies utiles pour cette analyse
library(tidyverse) # Facilite la manipulation de données
library(gt) # Aide à formater de jolis tableaux de rendu
library(broom) # Aide à formater les rendus de régressions
library(stargazer) # idem
library(sf) # Pour les données spatiales
library(lubridate) # Pour gérer des dates

# Désactiver les notations scientifiques
options(scipen =999)

load("data/ch3_AP_Vahatra.rds")

rct_AP_Mada <- AP_Vahatra %>%
  st_drop_geometry() %>%
  rename(`Déforestation 1996-2006 (%)` = 
           `Forest loss (ha) between 1996-2006 (percent loss)`,
         `Déforestation 2006-2016 (%)` = 
           `Forest loss (ha) between 2006-2016 (percent loss)`,
         `Surface (ha)` = hectares) %>%
  mutate(Groupe = ifelse(year(date_creation) < 2015, "Traitement", "Controle"),
         `Couvert forestier en 1996 (%)` = `Forest cover (ha) in 1996` / 
                                              `Surface (ha)` * 100,
         `Déforestation 1996-2016 (%)` = 
           (`Forest loss (ha) between 1996-2006 (absolute loss)` + 
             `Forest loss (ha) between 2006-2016 (absolute loss)`) /
           `Forest cover (ha) in 1996` * 100)

# On fait une série de tests de comparaison de moyenne
t_tests <- rct_AP_Mada %>% 
  # On applique aux variables de déforestation, couvert en 96 et taille
  summarise(across(ends_with("(%)") | ends_with("(ha)"),# toutes finissent ainsi
                   ~ t.test(.[Groupe == "Controle"], # on applique un t.test
                            .[Groupe == "Traitement"])$p.value)) %>%
  mutate(Groupe = "t-test")

equilibre_avant <- rct_AP_Mada %>%
  group_by(Groupe) %>%
  summarise(`Nombre d'aires` = n(),
            `Sans forêt` = sum(is.na(`Couvert forestier en 1996 (%)`)), 
            `Surface (ha)` = mean(`Surface (ha)`),
            `Couvert forestier en 1996 (%)` = 
              mean(`Couvert forestier en 1996 (%)`, na.rm = TRUE)) %>%
  bind_rows(t_tests) %>% # On colle tous les t-tests 
  mutate(across(!Groupe, ~round(., 2))) %>%# arrondit tout sauf colonne "Groupe"
  select(-starts_with("Déforestation"))# On ne garde que les t-tests de baseline

# Ce qui suit est une série d'opération pour formater le rendu en tableau
equilibre_avant %>%
  t() %>% # On transpose lignes <=> colonnes
  as.data.frame() %>% # La transposition a altéré le format, on remet en tableau
  tibble::rownames_to_column() %>% # On met le nom des lignes en 1° colonne
  # "Truc pour renommer avec le contenu de la première ligne
  `colnames<-` (filter(., row_number() == 1)) %>% 
  filter(row_number() != 1)%>% # Enlève la 1° ligne qui est maintenant en entête
  gt() %>%
  tab_header(title = "Equilibre des variables avant intervention",
             subtitle = "(exercice : \"comme si\" c'était une RCT)") %>%
  tab_source_note("Source : Association Vahatra et Carvalho et al. 2018")
Equilibre des variables avant intervention
(exercice : "comme si" c'était une RCT)
Groupe Controle Traitement t-test
Nombre d'aires 53 45 NA
Sans forêt 4 0 NA
Surface (ha) 66308.62 63513.89 0.88
Couvert forestier en 1996 (%) 54.31 62.24 0.14
Source : Association Vahatra et Carvalho et al. 2018

On a à première vue des déséquilibres limités “avant intervention”. En moyenne, les deux groupes sont assez proches en termes de surface et de couvert forestier et le test de Student ne permet pas de rejeter l’hypothèse nulle concernant une différence de moyenne sur ces critères.

On va maintenant s’intéresser aux différences de déforestation observées “après intervention” dans le groupe de traitement.

Code 
comparaison_apres <- rct_AP_Mada %>%
  group_by(Groupe) %>%
  summarise(across(starts_with("Déforestation"), ~ mean(., na.rm = TRUE))) %>%
  bind_rows(t_tests) %>% # On colle tous les t-tests 
  mutate(across(!Groupe, ~round(., 2))) %>%# arrondit tout sauf colonne "Groupe"
  select(Groupe, starts_with("Déforestation"))# On ne garde que les t-tests de baseline


# Même procédure que plus haut pour formater le rendu en tableau
comparaison_apres  %>%
  t() %>% # On transpose lignes <=> colonnes
  as.data.frame() %>% # La transposition a altéré le format, on remet en tableau
  tibble::rownames_to_column() %>% # On met le nom des lignes en 1° colonne
  # "Truc pour renommer avec le contenu de la première ligne
  `colnames<-` (filter(., row_number() == 1)) %>% 
  filter(row_number() != 1)%>% # Enlève la 1° ligne qui est maintenant en entête
  gt() %>%
  tab_header(title = "Moyennes après intervention",
             subtitle = "(exercice : \"comme si\" c'était une RCT)") %>%
  tab_source_note("Source : Association Vahatra et Carvalho et al. 2018")
Moyennes après intervention
(exercice : "comme si" c'était une RCT)
Groupe Controle Traitement t-test
Déforestation 1996-2006 (%) 8.21 4.12 0.07
Déforestation 2006-2016 (%) 8.50 3.91 0.01
Déforestation 1996-2016 (%) 15.90 7.87 0.00
Source : Association Vahatra et Carvalho et al. 2018

On peut également réaliser une régression simple, qu’on présente selon le format courant pour la littérature en économie grâce au package {stargazer} (Hlavac 2022).

Code 
def_96_06 <- lm(`Déforestation 1996-2006 (%)`  ~ Groupe, data = rct_AP_Mada)
def_06_16 <- lm(`Déforestation 2006-2016 (%)`  ~ Groupe, data = rct_AP_Mada)
def_96_16 <- lm(`Déforestation 1996-2016 (%)`  ~ Groupe, data = rct_AP_Mada)

rct_out1 <- stargazer(def_96_06, def_06_16, def_96_16, type = "html",
          title = "Impact de la conservation sur la perte de couvert forestier",
          notes = "Données : Association Vahatra et Carvalho et al. 2018") %>%
  str_replace_all("\\*", "\\\\*") 

# Dans un bloc plus bas et non affiché, on a le code suivant
# cat(rct_out1)
Impact de la conservation sur la perte de couvert forestier
Dependent variable:
Déforestation 1996-2006 (%) Déforestation 2006-2016 (%) Déforestation 1996-2016 (%)
(1) (2) (3)
GroupeTraitement -4.092* -4.593*** -8.031***
(2.308) (1.745) (2.831)
Constant 8.214*** 8.499*** 15.902***
(1.597) (1.208) (1.959)
Observations 94 94 94
R2 0.033 0.070 0.080
Adjusted R2 0.023 0.060 0.070
Residual Std. Error (df = 92) 11.179 8.454 13.713
F Statistic (df = 1; 92) 3.143* 6.924*** 8.046***
Note: *p<0.1; **p<0.05; ***p<0.01
Données : Association Vahatra et Carvalho et al. 2018

On analyse ensuite la relation aux variables topologiques (altitude, indice de terrain accidenté) et de temps de trajet à la ville la plus proche en 2015. Le seuil retenu ici pour considérer une localité comme une ville est qu’elle ait au moins 5000 habitants.

Code 
t_tests_autres <- rct_AP_Mada %>% 
  # On applique aux variables d'altitude, TRI et temps de trajet aux villes.
  summarise(across(indice_accidente:altitude,# toutes finissent ainsi
                   ~ t.test(.[Groupe == "Controle"], # on applique un t.test
                            .[Groupe == "Traitement"])$p.value)) %>% 
  mutate(Groupe = "t-test")


equilibre_autres <- rct_AP_Mada %>%
  group_by(Groupe) %>%
  summarise(`Nombre d'aires` = n(),
            indice_accidente = mean(indice_accidente, na.rm = TRUE), 
            dist_ville = mean(dist_ville, na.rm = TRUE),
            altitude = mean(altitude, na.rm = TRUE)) %>%
  bind_rows(t_tests_autres) %>% # On colle tous les t-tests 
  mutate(across(!Groupe, ~round(., 2))) %>%# arrondit tout sauf colonne "Groupe"
  select(-starts_with("Déforestation"))# On ne garde que les t-tests de baseline

equilibre_autres %>%
  gt() %>%
  tab_header(title = "Equilibre entre les groupes en matière topologique",
             subtitle = "(exercice : \"comme si\" c'était une RCT)") %>%
  tab_source_note("Source : Association Vahatra et Carvalho et al. 2018")
Equilibre entre les groupes en matière topologique
(exercice : "comme si" c'était une RCT)
Groupe Nombre d'aires indice_accidente dist_ville altitude
Controle 53 11.02 148.77 540.45
Traitement 45 12.08 233.31 538.63
t-test NA 0.53 0.03 0.99
Source : Association Vahatra et Carvalho et al. 2018
Code 
save(rct_AP_Mada, file = "data/rct_AP_Mada.rds")

Le temps de trajet aux villes est significativement distinct entre les deux groupes. Attention, car cette variable pose un problème d’endogénéité, car le jeu de données utilisé pour cela date de 2015, alors que notre période d’étude démarre en 1996. Or, il est possible que la présence d’aires protégées ait eu une incidence sur la construction ou l’amélioration des tronçons routiers à proximité (discussion en séance sur ce point). Il semble important de tenir compte de l’accessibilité géographiques des aires protégées, mais une donnée antérieure à 2015 serait préférable.

On essaye de limiter ce biais en ajoutant le temps de trajet à une ville comme variables de contrôle à notre régression.

Code 
def_96_16_controle <- lm(`Déforestation 1996-2016 (%)` ~ 
                           Groupe + dist_ville, 
                         data = rct_AP_Mada)

rct_out2 <- stargazer(def_96_16, def_96_16_controle, type = "html")  %>%
  str_replace_all("\\*", "\\\\*") 

# Dans un bloc plus bas et non affiché, on a le code suivant
# cat(rct_out1)
Dependent variable:
Déforestation 1996-2016 (%)
(1) (2)
GroupeTraitement -8.031*** -6.660**
(2.831) (2.958)
dist_ville -0.013
(0.008)
Constant 15.902*** 18.013***
(1.959) (2.341)
Observations 94 90
R2 0.080 0.103
Adjusted R2 0.070 0.082
Residual Std. Error 13.713 (df = 92) 13.669 (df = 87)
F Statistic 8.046*** (df = 1; 92) 4.989*** (df = 2; 87)
Note: *p<0.1; **p<0.05; ***p<0.01

Apparemment, le traitement reste significatif une fois que l’on contrôle pour la distance aux villes.