8 Méthodes d’appariement

8.1 La question de la comparabilité des groupes

On a vu dans le chapitre précédent que les comparaisons simples réalisées entre les premières et les dernières aires à avoir été formellement protégées pose problème.

On va maintenant chercher à renforcer la comparabilité entre le groupe de traitement et le groupe de contrôle en réalisant un appariemment (cf. diapos de présentation).

On va utiliser le package {MatchIt} : ne pas hésiter à se référer à la documentation du package.

On va commencer par réaliser quelques ajustements, car {MatchIt} requiert qu’aucune valeur des variables mobilisées ne soit manquante. On va donc retirer les observations comportant des NA.

Code

library(tidyverse) # Simplifie la manipulation de données
library(lubridate) # Simplifie les opérations sur des dates
library(sf) # Pour traiter les données spatiales
library(MatchIt) # Pour réaliser les appariements.
library(cobalt) # Pour les tests d'équilibre sur l'appariement
library(gt) # Pour faire de jolies tables
library(stargazer) # Pour présenter les résultats de régressions
library(mapme.biodiversity)
library(htmltools)

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

# Charger les données
AP_Vahatra <- read_rds("data/AP_Vahatra_mapme.rds") %>%
  portfolio_wide() %>%
  mutate(Groupe = ifelse(year(date_creation) < 2015, "Traitement", "Contrôle"))

# Préparer les données sans valeurs manquantes
Vahatra_defor_noNA <- AP_Vahatra %>%
  mutate(surface_ha = as.numeric(st_area(AP_Vahatra)) / 10000, 
         couv_foret_2000 = `treecover_area_2000-01-01_treecover_ha` / surface_ha * 100,
         altitude = `elevation_2000-02-01_elevation_mean_m`,
         indice_accidente = `slope_2000-02-01_slope_mean_degrees`,
         dist_ville = `traveltime_2015-01-01_5k_110mio_traveltime_mean_minutes`,
         traitement = ifelse(year(date_creation) < 2015, 1, 0),
         taux_deforestation_2000_2014 = 
           -((`treecover_area_2014-01-01_treecover_ha` / 
              `treecover_area_2000-01-01_treecover_ha`)^(1/14) - 1) * 100) %>%
  filter(!is.na(couv_foret_2000), !is.na(dist_ville), !is.na(altitude), !is.na(indice_accidente))

summary(Vahatra_defor_noNA)

Vahatra_defor_noNA %>%
  st_drop_geometry() %>%
  group_by(Groupe, traitement) %>%
  summarise(effectif = n())

Code

Vahatra_defor_noNA %>%
  st_drop_geometry() %>%
  group_by(Groupe) %>%
  summarize(`Nombre d'aires protégées` = n()) %>%
  gt() %>%
  tab_header("Observations par groupe avant appariemment") %>%
  tab_source_note("Source : Association Vahatra et données GFC")

Observations par groupe avant appariemment
Groupe	Nombre d'aires protégées
Contrôle	53
Traitement	44
Source : Association Vahatra et données GFC

8.2 Mesure de la propension

Pour commencer, nous allons spécifier le modèle probit qui estime dans quelle mesure la propension pour une aire d’avoir été protégée avant 2015 dépend de sa taille, de son taux de couverture forestière en 2000, de son altitude, de son caractère accidenté et de sa distance d’une ville d’au moins 5000 habitants.

Cette spécification peut se représenter selon l’équation suivante qui représente un modèle probit. Un modèle probit, tout comme le logit, est un modèle de choix binaire.

Code

# Spécification du modèle probit
pscor <- traitement ~  surface_ha + 
                       couv_foret_2000 + 
                       altitude +
                       indice_accidente + 
                       dist_ville

On va maintenant réaliser une régression pour connaître l’influence de ces facteurs dans la désignation des aires comme protégées.

Code

# Régression probit
reg_select <- glm(formula = pscor,
                  family = binomial(link = "probit"),
                  data = Vahatra_defor_noNA)

match_out1 <- stargazer(reg_select, type = "html")


	Dependent variable:

	traitement

surface_ha	-0.00000
	(0.00000)

couv_foret_2000	0.001
	(0.005)

altitude	-0.0004
	(0.0004)

indice_accidente	0.005
	(0.030)

dist_ville	0.003^***
	(0.001)

Constant	-0.547
	(0.343)


Observations	97
Log Likelihood	-60.920
Akaike Inf. Crit.	133.839

Note:	^p<0.1; ^p<0.05; ^**p<0.01

Exercice : analysez ce résultat. Quels facteurs sont corrélés avec la désignation précoce comme aire protégée ?

8.3 Appariement sur score de propension

On va maintenant utiliser ce modèle pour comparer les aires protégées traitées en premier par rapport à celles traitées plus récemment.

Code

# Calcul du matching
def_00_14_match <- matchit(formula = pscor,
                           family = binomial(link = "probit"),
                           method = "nearest",
                           discard = "both",
                           replace = FALSE,
                           distance = "glm",
                           data = Vahatra_defor_noNA)

print(def_00_14_match)

A matchit object
 - method: 1:1 nearest neighbor matching without replacement
 - distance: Propensity score [common support]
             - estimated with logistic regression
 - common support: units from both groups dropped
 - number of obs.: 97 (original), 78 (matched)
 - target estimand: ATT
 - covariates: surface_ha, couv_foret_2000, altitude, indice_accidente, dist_ville

On peut maintenant observer les équilibres entre les groupes traités et contrôle avant et après l’appariement.

Code

summary(def_00_14_match)


Call:
matchit(formula = pscor, data = Vahatra_defor_noNA, method = "nearest", 
    distance = "glm", discard = "both", replace = FALSE, family = binomial(link = "probit"))

Summary of Balance for All Data:
                 Means Treated Means Control Std. Mean Diff. Var. Ratio
distance                0.5143        0.4032          0.5970     1.9343
surface_ha          65093.8897    66539.7728         -0.0184     0.5595
couv_foret_2000        71.9537       62.1290          0.3384     0.6778
altitude              589.4707      518.9928          0.1459     0.6982
indice_accidente       11.0161        9.5870          0.2145     1.0678
dist_ville            260.4379      147.0780          0.5255     2.4284
                 eCDF Mean eCDF Max
distance            0.1840   0.3049
surface_ha          0.0993   0.2517
couv_foret_2000     0.0915   0.2333
altitude            0.0943   0.2294
indice_accidente    0.0729   0.1758
dist_ville          0.1822   0.3263

Summary of Balance for Matched Data:
                 Means Treated Means Control Std. Mean Diff. Var. Ratio
distance                0.4666        0.4483          0.0983     1.1599
surface_ha          67525.2840    53400.1694          0.1799     0.9542
couv_foret_2000        70.1568       72.2206         -0.0711     0.8134
altitude              563.7866      449.9603          0.2356     0.9791
indice_accidente       10.6267        9.8470          0.1171     1.0674
dist_ville            204.5025      168.8781          0.1651     0.9053
                 eCDF Mean eCDF Max Std. Pair Dist.
distance            0.0404   0.1795          0.1149
surface_ha          0.1049   0.2564          1.0223
couv_foret_2000     0.0634   0.2051          0.8333
altitude            0.0843   0.2308          1.0357
indice_accidente    0.0550   0.1538          1.0324
dist_ville          0.0941   0.2308          0.3906

Sample Sizes:
          Control Treated
All            53      44
Matched        39      39
Unmatched       9       0
Discarded       5       5

Exercice : Étudiez les tables ci-dessus. Quel effet a eu l’appariement sur l’équilibre des variables entre le groupe de traitement et le groupe de contrôle ? Combien d’observations ont été écartées ?

On peut observer la distance entre groupe de traitement et de contrôle.

Code

plot(def_00_14_match, type = "jitter", interactive = FALSE)

On peut également représenter l’équilibre entre les variables avant et après traitement avec les graphiques suivants.

Code

bal.plot(def_00_14_match, var.name = "dist_ville", which = "both")

Exercice : Quel effet a eu l’appariement sur la variable de distance à la ville ? Les autres variables d’appariement produisent-elles un effet aussi visible ?

8.4 Estimation du résultat en contrôlant pour les variables d’appariement

Le modèle qu’on utilise pour estimer l’impact est très proche de celui exposé ci-dessus, à la différence que la variable de traitement passe dans la partie droite, et qu’elle est remplacée par la déforestation.

Code

# Spécification du modèle pour l'impact
estimp <- taux_deforestation_2000_2014 ~   
                          traitement +
                          surface_ha + 
                          couv_foret_2000 + 
                          altitude +
                          indice_accidente + 
                          dist_ville

On va donc réaliser une régression, en tenant compte des pondérations générées par l’algorithme d’appariement (variable “weight”).

Code

# On extrait les données de l'appariement
def_00_14_match_data <- match.data(def_00_14_match)

# Régression avec pondérations
def_00_14_match_est <- lm(formula = estimp,
                          data = def_00_14_match_data,
                          weights = weights)

# Présentation des résultats
match_out2 <- stargazer(def_00_14_match_est, type = "html")


	Dependent variable:

	taux_deforestation_2000_2014

traitement	-0.278^*
	(0.158)

surface_ha	0.00000^***
	(0.00000)

couv_foret_2000	-0.001
	(0.003)

altitude	0.0001
	(0.0002)

indice_accidente	-0.009
	(0.017)

dist_ville	-0.001
	(0.001)

Constant	0.795^***
	(0.228)


Observations	78
R²	0.156
Adjusted R²	0.084
Residual Std. Error	0.689 (df = 71)
F Statistic	2.184^* (df = 6; 71)

Note:	^p<0.1; ^p<0.05; ^**p<0.01

8.5 Exercices

8.5.1 Exercice simple

Analysez, interprétez et critiquez les résultats ci-dessus.

8.5.2 Exercice intermédiaire

Ajoutez des variables d’intérêt et modifiez les paramètres de la fonction de matching.

8.5.3 Exercice avancé

Réalisez une analyse analogue avec les données de feux. Rédigez une analyse interprétative.