5  Expérimentations par assignation aléatoire

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

5.1 Mise en pratique

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.

5.2 Préparation des données

Pour commencer, on charge les librairies requises (après les avoir installées si nécessaire).

# 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
library(htmltools)
library(mapme.biodiversity)
library(units)
library(lmtest)      # Tests statistiques sur des modèles linéaires (waldtest, coeftest)
library(sandwich)    # Calcul d'erreurs standards robustes (vcovHC pour heteroskedasticité)



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

# On charge AP_Vahatra
AP_Vahatra <- read_rds("data/AP_Vahatra_mapme.rds")

5.3 Analyse des équilibres initiaux

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 2000.

# Définir les variables dans les bonnes unités
AP_RCT <- AP_Vahatra %>%
  portfolio_wide() %>% 
  rename(`Altitude moyenne (m)` = `elevation_2000-02-01_elevation_mean_m`,
         `Pente moyenne (degres)` = `slope_2000-02-01_slope_mean_degrees`,
         `Temps de trajet (minutes)` =
           `traveltime_2015-01-01_5k_110mio_traveltime_mean_minutes`,
         `Couvert forestier en 2000 (ha)` = 
           `treecover_area_2000-01-01_treecover_ha`,  
         `Couvert forestier en 2014 (ha)` = 
           `treecover_area_2014-01-01_treecover_ha`) %>%
  mutate(Groupe = ifelse(year(date_creation) < 2015, "Traitement", "Contrôle"),
         Traitement = ifelse(Groupe == "Traitement", 1, 0),
         `Surface (ha)` = as.numeric(st_area(AP_Vahatra)) / 10000,
         `Couvert forestier en 2000 (%)` = 
           (`Couvert forestier en 2000 (ha)` / `Surface (ha)`)*100)


# Définir la fonction pour les étoiles de significativité
stars <- function(p) {
  if (is.na(p)) return("")
  if (p < 0.01) return("***")
  else if (p < 0.05) return("**")
  else if (p < 0.10) return("*")
  else return("")
}


# Liste des covariables à inclure dans le tableau d’équilibre
covars <- c("Surface (ha)", "Couvert forestier en 2000 (ha)", "Couvert forestier en 2000 (%)", "Altitude moyenne (m)", "Pente moyenne (degres)", "Temps de trajet (minutes)")

# Création du tableau d’équilibre
balance_data <- map_dfr(covars, function(var) {
    
      # Moyennes et écart-types par groupe
      stats <- AP_RCT %>%
        group_by(Groupe) %>%
        summarise(
          mean = mean(.data[[var]], na.rm = TRUE),
          sd = sd(.data[[var]], na.rm = TRUE),
          n = sum(!is.na(.data[[var]]))
        ) %>%
        ungroup()
    
      # Extraire les statistiques pour les deux groupes
      treat <- stats %>% filter(Groupe == "Traitement")
      control <- stats %>% filter(Groupe == "Contrôle")
      
      # Différence brute
      diff <- treat$mean - control$mean
      
      # Estimation de l’erreur standard via régression avec erreurs robustes
      mod <- lm(as.formula(paste0("`", var, "` ~ Groupe")), data = AP_RCT)
      coef_test <- coeftest(mod, vcov. = vcovHC(mod, type = "HC1"))
      se_diff <- coef_test["GroupeTraitement", "Std. Error"]
        p_val <- coef_test["GroupeTraitement", "Pr(>|t|)"]
        etoiles <- stars(p_val)
      
      # Calcul de la différence standardisée
      pooled_sd <- sqrt(0.5 * (treat$sd^2 + control$sd^2))
      std_diff <- diff / pooled_sd
    
  # Construction de la ligne du tableau
      tibble(
        Variable = var,
        Moy_Traitement = 
          str_glue("{round(treat$mean, 2)}\n({round(treat$sd, 2)})"),
        Moy_Contrôle = 
          str_glue("{round(control$mean, 2)}\n({round(control$sd, 2)})"),
        Diff = str_glue("{round(diff, 2)}{etoiles}\n[{round(se_diff, 2)}]"),
        Diff_std = str_glue("{round(std_diff, 2)}")
      )
}) %>% bind_rows()
  
  # Test de significativité conjointe des covariables (p-value du F-test)
  f_test_mod <- lm(as.formula(paste("Traitement ~", paste0("`", covars, "`", collapse = " + "))), data = AP_RCT)
  f_test <- waldtest(f_test_mod, vcov = vcovHC(f_test_mod, type = "HC1"))
  p_val_joint <- signif(f_test$`Pr(>F)`[2], 3)
  
  # Ajouter une ligne au tableau avec la p-value du F-test
  balance_data <- bind_rows(
    balance_data,
    tibble(
      Variable = "P-value F-test conjoint",
      Moy_Traitement = "", Moy_Contrôle = "",
      Diff = sprintf("p = %.3f", p_val_joint),
      Diff_std = ""
    )
  )

# Affichage final
balance_data %>%
  gt() %>%
  tab_header(
    title = "Tableau de contrôle d'équilibre initial entre groupes",
    subtitle = "(exercice : \"comme si\" c'était un RCT)"
  ) %>%
  cols_label(
    Variable = "Variable",
    Moy_Traitement = "Traitement\n(SD)",
    Moy_Contrôle = "Contrôle\n(SD)",
    Diff = "Différence\n[SE]",
    Diff_std = "Différence\nstandardisée"
  ) %>%
  tab_source_note("Source : Association Vahatra et données GFC (*: p < 0,1 ; **: p < 0,05 ; ***: p < 0,01)")

Tableau de contrôle d'équilibre initial entre groupes
(exercice : "comme si" c'était un RCT)
Variable	Traitement (SD)	Contrôle (SD)	Différence [SE]	Différence standardisée
Surface (ha)	63728 (78166.99)	66539.77 (104980.3)	-2811.77 [18541.73]	-0.03
Couvert forestier en 2000 (ha)	45507.19 (69945.04)	38974.21 (73268.24)	6532.98 [14574.56]	0.09
Couvert forestier en 2000 (%)	71.95 (29.03)	62.13 (35.26)	9.82 [6.53]	0.3
Altitude moyenne (m)	578.9 (482.84)	518.99 (578.2)	59.9 [107.19]	0.11
Pente moyenne (degres)	10.93 (6.61)	9.59 (6.45)	1.34 [1.32]	0.21
Temps de trajet (minutes)	255.12 (216.23)	147.08 (138.44)	108.04*** [37.41]	0.6
P-value F-test conjoint			p = 0.004
Source : Association Vahatra et données GFC (*: p < 0,1 ; **: p < 0,05 ; ***: p < 0,01)

Exercice : Analysez le résultat de cette table.

5.4 Différences de déforestation “après intervention”

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

# Calculer le taux de déforestation pour la période 2000-2014

AP_RCT <- AP_RCT %>%
  mutate(taux_deforestation_2000_2014 = 
           -((`Couvert forestier en 2014 (ha)` - 
              `Couvert forestier en 2000 (ha)`) /
           `Couvert forestier en 2000 (ha)`)* 100)

# Comparer les taux de déforestation entre les groupes
comparaison_deforestation <- AP_RCT %>%
  st_drop_geometry() %>%
  group_by(Groupe) %>%
  summarise(`Taux de déforestation 2000-2014 (%)` = 
              mean(taux_deforestation_2000_2014, na.rm = TRUE)) %>%
  mutate(across(!Groupe, ~round(., 2)))

# Ce qui suit est une série d'opérations pour formater le rendu en tableau
comparaison_deforestation %>%
  t() %>%
  as.data.frame() %>%
  tibble::rownames_to_column() %>%
  `colnames<-` (filter(., row_number() == 1)) %>%
  filter(row_number() != 1) %>%
  gt() %>%
  tab_header(title = "Moyennes des taux de déforestation",
             subtitle = "(exercice : \"comme si\" c'était une RCT)") %>%
  tab_source_note("Source : Association Vahatra et données GFC")

Moyennes des taux de déforestation
(exercice : "comme si" c'était une RCT)
Groupe	Contrôle	Traitement
Taux de déforestation 2000-2014 (%)	8.76	4.56
Source : Association Vahatra et données GFC

Exercice : Commentez le résultat de cette table.

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).

# On exécute une régression pour la période 2000-2014
def_2000_2014 <- lm(taux_deforestation_2000_2014 ~ Groupe, data = AP_RCT)

# On consolide les résultats des régressions dans une table qu'on formate
# avec le package stargazer
rct_out1 <- stargazer(def_2000_2014, type = "html",
          title = "Impact de la conservation sur la perte de couvert forestier",
          notes = "Données : Association Vahatra et données GFC")  

Impact de la conservation sur la perte de couvert forestier

	Dependent variable:
	taux_deforestation_2000_2014
GroupeTraitement	-4.201**
	(1.777)
Constant	8.757***
	(1.197)
Observations	97
R2	0.056
Adjusted R2	0.046
Residual Std. Error	8.714 (df = 95)
F Statistic	5.588** (df = 1; 95)
Note:	*p<0.1; **p<0.05; ***p<0.01
	Données : Association Vahatra et données GFC

Exercice : Analysez le résultat de cette table. Qu’est-ce qu’elle suggère ?

Rappel sur les régressions

Application finale :

• Régressez à nouveau le taux de déforestation entre 2000 et 2014 sur la variable de traitement en ajoutant le taux de déforestation en 2000 comme variable de contrôle (vous ajoutez un deuxième régresseur en plus de ‘Groupe’).

• Utilisez le même modèle pour estimer l’effet du traitement sur le pourcentage de couvert forestier en 2023.

Hlavac, Marek. 2022. Stargazer: Well-Formatted Regression and Summary Statistics Tables. https://CRAN.R-project.org/package=stargazer.