8 Méthodes d’appariement

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

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) # Pour des graphs qui s'affichent bien sur une page web
library(units) # pour passer de m2 à km2

# 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_m2 = as.numeric(st_area(.)),
         surface_ha = surface_m2 / 10000,
         surface_km2 = surface_m2 / 1000000,
         couv_foret_2000 = `treecover_area_2000-01-01_treecover_ha` / surface_ha * 100,
         altitude = `elevation_2000-02-01_elevation_mean_m`,
         pente = `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_annuel_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(pente))

summary(Vahatra_defor_noNA)

Vahatra_defor_noNA %>%
  st_drop_geometry() %>%
  group_by(Groupe, traitement) %>%
  summarise(effectif = n())
 
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")

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_km2 + 
                       couv_foret_2000 + 
                       altitude +
                       pente + 
                       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_km2	-0.0002
	(0.0002)

couv_foret_2000	0.001
	(0.005)

altitude	-0.0004
	(0.0004)

pente	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 = "logit"),
                           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_km2, couv_foret_2000, altitude, pente, 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 = "logit"))

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_km2          650.9389      665.3977         -0.0184     0.5595
couv_foret_2000       71.9537       62.1290          0.3384     0.6778
altitude             589.4707      518.9928          0.1459     0.6982
pente                 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_km2        0.0993   0.2517
couv_foret_2000    0.0915   0.2333
altitude           0.0943   0.2294
pente              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_km2          675.2528      534.0017          0.1799     0.9542
couv_foret_2000       70.1568       72.2206         -0.0711     0.8134
altitude             563.7866      449.9603          0.2356     0.9791
pente                 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_km2        0.1049   0.2564          1.0223
couv_foret_2000    0.0634   0.2051          0.8333
altitude           0.0843   0.2308          1.0357
pente              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_annuel_deforestation_2000_2014 ~   
                          traitement +
                          surface_km2 + 
                          couv_foret_2000 + 
                          altitude +
                          pente + 
                          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_annuel_deforestation_2000_2014

traitement	-0.278^*
	(0.158)

surface_km2	0.0003^***
	(0.0001)

couv_foret_2000	-0.001
	(0.003)

altitude	0.0001
	(0.0002)

pente	-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

Retirez du score de propension et du calcul la variable de surface des aires protégées et réestimez le résultat.

8.5.3 Exercice avancé

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

8.6 Application sur les données en mailles

L’exercice ci-dessus consiste à réaliser un appariement sur des données d’une taille relativement grande. Cet exercice comporte des limites, car on dispose d’un nombre limité d’observations à comparer.

Une approche plus appropriée consiste à utiliser le maillage du territoire effectué dans la section Chapter 7 pour comparer des parcelles d’aires passées sous protection pendant la période d’étude (2000-2023) à des zones d’ayant jamais été protégées.

Les données ne peuvent pas contenir de données manquantes sur les variables d’appariement, donc on les écarte.

Code

library(tidyverse)
library(mapme.biodiversity)
library(MatchIt)
library(stargazer)
library(sf)
library(cobalt)
library(tmap)
library(htmltools) # Pour avoir de plus jolies tables
library(geodata) # Pour les frontières de Madagascar

grille_mada_summary_AP <- read_rds("data/grille_mada_summary_AP.rds") %>%
  # on passe le taux de déforestation en annuel
  mutate(taux_annuel_deforestation_2000_2023 = 
           (1 - (1 - taux_deforestation_2000_2023 / 100)^(1/22)) * 100)


# On référence le nom des variables qui vont servir à l'analyse
variables_analyse <- c("assetid",
                       "traitement",
                       "couv_foret_2000",
                       "pop_km2",
                       "altitude",
                       "pente",
                       "dist_ville",
                       "taux_annuel_deforestation_2000_2023")

# On renomme le ficher 'df' (dataframe) : plus concis dans les commandes ensuite
df <- grille_mada_summary_AP %>%
  # On supprime toutes les lignes pour lesquelles au moins 1 valeur variable 
  # est manquante parmi les variables d'analyse
  mutate(traitement = position_ap == "Intérieur") %>% 
  drop_na(any_of(variables_analyse))

On analyse maintenant le score de propension.

Code

# Get propensity scores
glm_out <- glm(traitement ~ couv_foret_2000 + 
                           pop_km2 + 
                           altitude + 
                           pente + 
                           dist_ville,  
               family = binomial(link = "logit"),
               data = df)

cellmatch_out1 <- stargazer(glm_out,
                            summary = TRUE,
                            type = "html",
                            title = "Score de propension")

**Score de propension**

	Dependent variable:

	traitement

couv_foret_2000	0.025^***
	(0.0005)

pop_km2	-0.047^***
	(0.001)

altitude	0.001^***
	(0.00004)

pente	-0.075^***
	(0.003)

dist_ville	0.001^***
	(0.0001)

Constant	-2.820^***
	(0.034)


Observations	58,509
Log Likelihood	-18,035.800
Akaike Inf. Crit.	36,083.610

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

Exercice : interpréter le résultat du score de propension.

Code

m_out <- matchit(traitement ~ couv_foret_2000 + 
                           pop_km2 + 
                           altitude + 
                           pente + 
                           dist_ville,
                 data = df,
                 method = "nearest",
                 replace = TRUE,
                 distance = "glm", 
                 discard = "both", # common support: drop units from both groups 
                 link = "probit")

print(m_out)

A `matchit` object
 - method: 1:1 nearest neighbor matching with replacement
 - distance: Propensity score [common support]

             - estimated with probit regression
 - common support: units from both groups dropped
 - number of obs.: 58509 (original), 12666 (matched)
 - target estimand: ATT
 - covariates: couv_foret_2000, pop_km2, altitude, pente, dist_ville

Code

# print(summary(m_out, un = FALSE))
bal_table <- bal.tab(m_out, un = TRUE)
print(bal_table)

Balance Measures
                    Type Diff.Un Diff.Adj
distance        Distance  0.9619   0.0000
couv_foret_2000  Contin.  0.7784  -0.0014
pop_km2          Contin. -1.2824  -0.0784
altitude         Contin.  0.3385  -0.0258
pente            Contin.  0.1807   0.0346
dist_ville       Contin.  0.4795  -0.0274

Sample sizes
                      Control Treated
All                  51295.      7214
Matched (ESS)         3904.59    7199
Matched (Unweighted)  5467.      7199
Unmatched            45293.         0
Discarded              535.        15

Code

grille_matched <- match.data(m_out) %>%
  st_sf()

if (!file.exists("data/grille_matched.rds")) {
  write_rds(grille_matched, "data/grille_matched.rds")
}

contour_mada <- gadm("Madagascar", level = 0, path = "data") %>% st_as_sf()
# On visualise les données appareillées
tm_shape(contour_mada) +
  tm_borders() +
  tm_shape(grille_matched) +
  tm_fill(
    fill = "traitement",
    fill.scale = tm_scale(values = "brewer.set1", labels = c("Contrôle", "Traitement")),
    fill.legend = tm_legend(title = "Groupes post-appariement")
  ) +
  tm_layout(legend.outside = TRUE)

Exercice: Réaliser des tests d’équilibre

Code

summary(m_out)


Call:
matchit(formula = traitement ~ couv_foret_2000 + pop_km2 + altitude + 
    pente + dist_ville, data = df, method = "nearest", distance = "glm", 
    link = "probit", discard = "both", replace = TRUE)

Summary of Balance for All Data:
                Means Treated Means Control Std. Mean Diff. Var. Ratio
distance               0.2507        0.1059          0.9619     2.2797
couv_foret_2000       75.9817       50.1663          0.7784     0.7665
pop_km2               11.3610       27.9297         -1.2824     0.0178
altitude             617.3927      433.9686          0.3385     1.5847
pente                 10.7017        9.4236          0.1807     1.4981
dist_ville           333.0405      198.4916          0.4795     2.3677
                eCDF Mean eCDF Max
distance           0.2858   0.4437
couv_foret_2000    0.2238   0.3156
pop_km2            0.1280   0.2337
altitude           0.0859   0.1973
pente              0.0713   0.1081
dist_ville         0.1461   0.2047

Summary of Balance for Matched Data:
                Means Treated Means Control Std. Mean Diff. Var. Ratio
distance               0.2497        0.2497          0.0000     0.9999
couv_foret_2000       75.9338       75.9786         -0.0014     0.9966
pop_km2               11.3641       12.3766         -0.0784     0.8900
altitude             614.4045      628.3593         -0.0258     1.4263
pente                 10.6975       10.4530          0.0346     1.6646
dist_ville           331.8206      339.4963         -0.0274     1.0271
                eCDF Mean eCDF Max Std. Pair Dist.
distance           0.0000   0.0015          0.0004
couv_foret_2000    0.0256   0.1043          0.5019
pop_km2            0.0213   0.0440          0.7552
altitude           0.0659   0.1064          0.8127
pente              0.0646   0.1181          0.9690
dist_ville         0.0217   0.0658          0.8226

Sample Sizes:
               Control Treated
All           51295.      7214
Matched (ESS)  3904.59    7199
Matched        5467.      7199
Unmatched     45293.         0
Discarded       535.        15

Code

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

Code

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

Code

bal.plot(m_out, var.name = "couv_foret_2000", which = "both")

On réalise la régression.

Code

modele <- lm(formula = taux_annuel_deforestation_2000_2023 ~
               traitement +
               couv_foret_2000 + 
               pop_km2 + 
               altitude + 
               pente + 
               dist_ville,
             data = grille_matched,
             weights = weights)
cellmatch_out2 <- stargazer(modele, type = "html")


	Dependent variable:

	taux_annuel_deforestation_2000_2023

traitement	-0.783^***
	(0.050)

couv_foret_2000	0.013^***
	(0.001)

pop_km2	0.002
	(0.002)

altitude	-0.001^***
	(0.0001)

pente	-0.019^***
	(0.005)

dist_ville	-0.001^***
	(0.0001)

Constant	2.145^***
	(0.072)


Observations	12,666
R²	0.072
Adjusted R²	0.071
Residual Std. Error	2.801 (df = 12659)
F Statistic	162.580^*** (df = 6; 12659)

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

Exercice : interprétez le résultat