9  Estimation with 2X2 and staggered DiD

9.1 Objectives

Le document présente l’estimation en double différence des groupes de traitement et de contrôle pour la période 2008 à 2021. Le principe est de comparer l’indice de richesse des ménages de contrôle et de traitement à une période antérieure et à une période postérieure à la mise en place des aires protégées. La période de pré-traitement sera la période de 1997 à 2008, et la période de post-traitement sera la période 2008 et 2021. Le groupe de traitement sont les ménages vivant à moins de 10km des aires protégées créées entre 2008 et 2021. Le groupe de contrôle sont les ménages vivant à plus de 10km des aires protégées créées jusqu’en 2021. Les ménages vivant à moins de 10km d’une aire protégée créée avant 2008 seront exclus de l’analyse.

Dans nos estimations par DID, nous avons utilisé des erreurs standards robustes aux corrélations intra-grappes. Cette estimation permet de gérer la dépendance des observations au sein des grappes. Le but est d’ajuster les erreurs standards en tenant compte des corrélations intra-grappes présentes au sein des groupes par l’échantillonnage de la variable de résultat Abadie (2021). Cet ajustement se calcule sur la base des résidus du modèle et d’une estimation de variance intra-cluster. On calcule les résidus à partir du modèle de régression puis on utilise ces résidus pour construire une matrice variance-covariance robuste. La matrice va ajuster la corrélation intra-grappe en tenant compte la structure corporelle des erreurs.

Reminder, we have the following hypothesis:

• H1-overall effect on livelihoods: PA reduce the living standards of nearby households by restricting access to natural resources, with insufficient compensation or benefits.

• H2-effect on inequalities: PA exacerbate economic inequalities, as better-off or better-connected individuals capture most of the benefits (tourism jobs, development projects).

• H3-heterogeneity: Impacts vary across PA depending on governance, community participation, and management approaches.

9.2 Overall effect on livelihoods

9.2.1 2X2 DiD

In our PAP, we stated:

“The difference-in-difference principle is to compare the wealth index of control and treatment households before and after the establishment of PA. Our treatment began in 2011, so we use 2008 DHS data for the pre-treatment and 2021 DHS data for the post-treatment. This method relies on the parallel trends assumption, meaning that, in the absence of treatment, treatment and control groups would have experienced similar changes over time. To verify this assumption, we include two pre-treatment, points (DHS 1997 and DHS 2008 data) and assess trends between these years. […] The difference-in-difference method relies on the assumption of parallel trends. To validate this assumption, we will use as a reference, among the rural households surveyed in 1997, those living in an area located in or within of a PA created between 2008 and 2021 (placebo treatment group) and those matched to them using the method described above (placebo control group). We will graphically represent the evolution of the treatment and control groups in 1997, 2008 and 2021 for visual confirmation. To confirm the validity of the parallel trends assumption, we will conduct a placebo test 1997-2008 outcomes, as well as matching and control variables, while defining the treatment variable for post-2008 PA. If the estimated effect between 1997 and 2008 is null or statistically insignificant, this supports the parallel trends assumption.”

library(tidyverse)
library(haven)
library(fixest)
library(didimputation)
library(broom)
library(ggplot2)

# Chargement des données
d97 <- read_rds("data/derived/data_matched_1997.rds") %>%
  rename(spei_wc_n_2 = spei_wc_1995,
         spei_wc_n_1 = spei_wc_1996,
         spei_wc_n   = spei_wc_1997) %>%
  mutate(hv219 = zap_labels(hv219), # hhh sex (1/2)
         hv220 = zap_labels(hv220)) # hhh age (num)

d08 <- read_rds("data/derived/data_matched_2008.rds") %>%
  rename(spei_wc_n_2 = spei_wc_2006,
         spei_wc_n_1 = spei_wc_2007,
         spei_wc_n   = spei_wc_2008) %>%
  mutate(hv219 = zap_labels(hv219),
         hv220 = zap_labels(hv220))

d11 <- read_rds("data/derived/data_matched_2011.rds") %>%
  rename(spei_wc_n_2 = spei_wc_2009,
         spei_wc_n_1 = spei_wc_2010,
         spei_wc_n   = spei_wc_2011) %>%
  mutate(hv219 = zap_labels(hv219),
         hv220 = zap_labels(hv220))

d13 <- read_rds("data/derived/data_matched_2013.rds") %>%
  rename(spei_wc_n_2 = spei_wc_2011,
         spei_wc_n_1 = spei_wc_2012,
         spei_wc_n   = spei_wc_2013) %>%
  mutate(hv219 = zap_labels(hv219),
         hv220 = zap_labels(hv220))

d16 <- read_rds("data/derived/data_matched_2016.rds") %>%
  rename(spei_wc_n_2 = spei_wc_2014,
         spei_wc_n_1 = spei_wc_2015,
         spei_wc_n   = spei_wc_2016) %>%
  mutate(hv219 = zap_labels(hv219),
         hv220 = zap_labels(hv220))

d21 <- read_rds("data/derived/data_matched_2021.rds") %>%
  rename(spei_wc_n_2 = spei_wc_2019,
         spei_wc_n_1 = spei_wc_2020,
         spei_wc_n   = spei_wc_2021) %>%
  mutate(hv219 = zap_labels(hv219),
         hv220 = zap_labels(hv220))

# Préparation des données
dat <- bind_rows(d97, d08, d11, d13, d16, d21) %>%
  filter(GROUP %in% c("Treatment","Control")) %>%
  mutate(
    hv219   = factor(hv219, levels = c(1,2), labels = c("Homme","Femme")), # sexe (cat.)
    hv220   = as.numeric(hv220),                                           # âge
    treat   = as.integer(GROUP == "Treatment"),
    w_svy   = hv005 / 1e6,
    w_all   = w_svy * weights, # poids d'enquête × poids de matching (si 'weights' existe)
    id      = row_number(),
    # Map des années de statut -> première année d'observation post (treatment_phase)
    treatment_phase = case_when(
      STATUS_YR == 2010 ~ 2011,
      STATUS_YR == 2012 ~ 2013,
      STATUS_YR == 2015 ~ 2016,
      STATUS_YR == 2017 ~ 2021,
      is.na(STATUS_YR)  ~ 0,
      TRUE               ~ STATUS_YR
    )
  )

# Outcome h
yvar <- "wealth_centile_rural_weighted"

# fixest DID 2×2: placebo 1997–2008, traitement 2008–2021 ------------------

# Placebo 1997–2008
pre <- dat %>%
  filter(DHSYEAR %in% c(1997, 2008)) %>%
  mutate(post = as.integer(DHSYEAR == 2008),
         treat_post = treat * post)

f_pre <- as.formula(paste(
  yvar, "~ treat + post + treat_post +",
  "spei_wc_n_2 + spei_wc_n_1 + spei_wc_n + hv219 + hv220"
))

m_pre <- feols(f_pre, data = pre, weights = ~ w_all, cluster = ~ hv001)

# Traitement 2008–2021
main <- dat %>%
  filter(DHSYEAR %in% c(2008, 2021)) %>%
  mutate(post = as.integer(DHSYEAR == 2021),
         treat_post = treat * post)

f_main <- f_pre  # même formule

m_main <- feols(f_main, data = main, weights = ~ w_all, cluster = ~ hv001)

etable(m_pre, m_main, headers = c("Placebo 97–08", "Traitement 08–21"))

                                        m_pre                        m_main
                              Placebo 97–08            Traitement 08–21
Dependent Var.: wealth_centile_rural_weighted wealth_centile_rural_weighted
                                                                           
Constant                     44.96*** (4.313)              37.89*** (2.715)
treat                          -1.797 (4.156)                 1.484 (3.454)
post                           -6.394 (6.043)              22.55*** (4.543)
treat_post                      3.736 (5.392)                -2.751 (4.300)
spei_wc_n_2                     2.114 (3.907)               -5.482. (3.310)
spei_wc_n_1                    -1.055 (2.639)             -8.214*** (2.203)
spei_wc_n                       2.411 (4.317)              17.72*** (2.512)
hv219Femme                   -1.735. (0.9367)            -3.268*** (0.7580)
hv220                      0.0966*** (0.0255)            0.0921*** (0.0182)
_______________ _____________________________ _____________________________
S.E.: Clustered                     by: hv001                     by: hv001
Observations                            6,831                        11,686
R2                                    0.00655                       0.08010
Adj. R2                               0.00539                       0.07948
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Extraction compacte des deux effets DID
did_row <- function(model, year_post, vc = ~ hv001, term = "treat_post"){
  summary(model, vcov = vc) %>% broom::tidy() %>%
    filter(term == !!term) %>%
    transmute(year = year_post, estimate, se = std.error)
}
did_df <- bind_rows(
  did_row(m_pre, 2008),
  did_row(m_main, 2021)
) %>%
  mutate(period = factor(ifelse(year == 2008, "1997–2008", "2008–2021"),
                         levels = c("1997–2008", "2008–2021")),
         lo = estimate - 1.96*se,
         hi = estimate + 1.96*se)

ggplot(did_df, aes(x = period, y = estimate)) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_errorbar(aes(ymin = lo, ymax = hi), width = .2) +
  geom_point(size = 3) +
  labs(x = NULL, y = "Effet DID sur le centile de richesse (pondéré)",
       title = "DID 2×2 avec IC clusterisés (hv001) - buffer 10 km")

Le test placebo donne une estimation avec un coefficient faible (+1.8 centile) et non significatif: pas de différence détectée avant intervention dans les tendances du groupe traitement et du groupe contrôle. L’estimation du traitement présente une plus grande magnitude (-4.1 centiles), mais il reste non significatif. Les ménages vivant près des AP semblent connaître, un niveau de bien-être légèrement inférieur à celui des ménages contrôle, mais l’effet n’est pas significatif. A ce stade, l’analyse ne confirme pas l’hypothèse 1 (que les AP réduise le niveau de vie) : le sens de la relation est négatif, mais on ne peut pas rejeter l’hypothèse nulle. L’intervalle de confiance couvre encore largement 0, ce qui signifie qu’on a une précision insuffisante pour détecter un effet. L’effet est donc soit nul, soit inférieur à notre capacité de détection. Dans le PAP, on avait estimé que l’effet devrait être d’au moins 7.5 centiles pour pouvoir être détecté.

9.2.2 Staggered diff-in-diff

In our PAP, we stated:

“The MIS data can therefore be mobilized for this study. However, the periods covered by MIS only allow for comparisons beyond a simple “before” and “after” framework, which complicates the analysis and necessitates methods with limited methodological consensus. Since 2020, the state of the art in difference-in-difference methods has been challenged by the recognition that classical two-way fixed effects methods can yield spurious results when the effects of an intervention are heterogeneous and assessed over multiple periods. Since then, a dozen alternative approaches have been developed, but econometricians continue to debate their reliability (Roth et al. 2023). To avoid undermining the credibility of our results, we will initially limit ourselves to two study periods (DHS 2008 and 2021). Only if this approach lacks sufficient statistical power will we consider incorporating additional data. […] Additionally, if we incorporate data from MIS 1997, MIS 2011, and MIS 2013, we could apply an estimator suitable for staggered adoption and multiple study periods accounting for potential heterogeneous treatment effects (Borusyak, Jaravel, and Spiess 2024). This will only be performed if power is insufficient with two periods of data.”

Pour l’instant, on n’arrive pas à appliquer l’approche de Borusyak, Jaravel, and Spiess (2024) avec des données en coupes répétées (pas du panel) avec le package R didimputation. En théorie cette méthode peut s’appliquer à des données “repeated cross sections” et il semblerait qu’un package STATA le fasse. Mais nous ne disposons pas de STATA pour tester. On va donc appliquer la méthode de Gardner, qui en principe est très proche (avec des intervalles de confiance un peu plus large).

Nota bene : Afin de stabiliser les estimations de l’event-study, nous avons appliqué un binning des années relatives au traitement, regroupant toutes les observations au-delà de ±5 années dans des catégories communes. Ce choix est motivé par l’espacement irrégulier des enquêtes disponibles, qui implique que certaines valeurs de temps relatif sont représentées par très peu d’observations et produisent des coefficients instables. Comme le soulignent Borusyak, Jaravel et Spiess (2024, p. 22) : “In practice, it is common to bin distant leads and lags into a single category, both to improve statistical precision and to avoid presenting very noisy and uninformative coefficients.” Cette approche permet ainsi de réduire le bruit aux extrêmes et de se concentrer sur l’interprétation des dynamiques principales autour du traitement.s

# did2s (Gardner) : statique + event-study--------------------------------
library(did2s)

dat3 <- dat %>%
  mutate(
    treat_on = as.integer(treatment_phase > 0 & DHSYEAR >= treatment_phase),
    rel_year = if_else(treatment_phase > 0, DHSYEAR - treatment_phase, Inf),
    # Binning prudent pour stabilité (-5..5)
    rel_year_binned = pmax(pmin(rel_year, 5), -5)
  )

# -- Statique (traitement "on/off")
did2s_static <- did2s(
  data        = dat3,
  yname       = yvar,
  first_stage = ~ spei_wc_n_2 + spei_wc_n_1 + spei_wc_n + hv219 + hv220 | DHSYEAR,
  second_stage= ~ i(treat_on, ref = FALSE),
  treatment   = "treat_on",
  cluster_var = "hv001",
  weights     = "w_all"
)
etable(did2s_static, headers = "did2s statique")

                                 did2s_static
                               did2s statique
Dependent Var.: wealth_centile_rural_weighted
                                             
treat_on = 1                    1.720 (2.156)
_______________ _____________________________
S.E.: Clustered                     by: hv001
Observations                           21,923
R2                                    0.00071
Adj. R2                               0.00071
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# -- Event-study (avec binning -5..5, ref = -1 et never-treated)
did2s_es <- did2s(
  data        = dat3,
  yname       = yvar,
  first_stage = ~ spei_wc_n_2 + spei_wc_n_1 + spei_wc_n + hv219 + hv220 | DHSYEAR,
  second_stage= ~ i(rel_year_binned, ref = c(-1, Inf)),
  treatment   = "treat_on",
  cluster_var = "hv001",
  weights     = "w_all"
)
etable(did2s_es, headers = "did2s event-study (-5..5)")

                                          did2s_es
                         did2s event-study (-5..5)
Dependent Var.:      wealth_centile_rural_weighted
                                                  
rel_year_binned = -5                0.7757 (1.305)
rel_year_binned = -3                -2.535 (2.164)
rel_year_binned = 0                  3.944 (3.579)
rel_year_binned = 2                -0.3823 (4.957)
rel_year_binned = 3               13.95*** (3.322)
rel_year_binned = 5                0.0375 (0.7943)
____________________ _____________________________
S.E.: Clustered                          by: hv001
Observations                                21,923
R2                                         0.00378
Adj. R2                                    0.00355
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Plot ES did2s
plot_did2s <- broom::tidy(did2s_es, conf.int = TRUE) %>%
  filter(grepl("^rel_year_binned::", term)) %>%
  mutate(k = as.numeric(sub("^rel_year_binned::", "", term))) %>%
  arrange(k)

ggplot(plot_did2s, aes(k, estimate)) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_point() +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = .2) +
  labs(x = "Années relatives au 1er traitement (binnées -5..5)",
       y = "Effet DID sur le centile de richesse (pondéré)",
       title = "Event-study (did2s) - buffer de 10 km") + 
  scale_x_continuous(breaks = c(-5,-4,-3,-2,-1,0,1,2,3,4,5))

# Test joint de pré-tendances (tous les leads k < 0)
# On construit un motif regex pour k = -5..-1 présents dans le modèle
leads_present <- plot_did2s$k[plot_did2s$k < 0]
if(length(leads_present) > 0){
  keep_regex <- paste0("^rel_year_binned::(", paste(leads_present, collapse="|"), ")$")
  print(fixest::wald(did2s_es, keep = keep_regex))
}

Wald test, H0: joint nullity of rel_year_binned::-5 and rel_year_binned::-3
 stat = 0.866281, p-value = 0.420527, on 2 and 21,917 DoF, VCOV: Corrected Clustered (hv001).$stat
[1] 0.8662812

$p
[1] 0.4205269

$df1
[1] 2

$df2
[1] 21917

$vcov
[1] "Corrected Clustered (hv001)"

L’estimation en “statique” donne un coefficient de ~1 centile (0.99), avec une erreur type de 2.1, ce qui indique un résultat non significatif. L’analyse en série temporelle avec binning ±5 ans. Etaye l’hypothèse de tendances parallèle. Les estimations en post semblent apparaître, mais ne semblent pas persister. Cela abonde dans le sens d’un impact faible et non-significatif. De possibles avancées à court-terme (~3 ans après la création de l’AP), mais pas soutenues dans le temps.

9.2.3 Quantile treatment effects

Le quantile treatment effect mesure l’effet du traitement à différents points de la distribution de la variable de résultat, le wealth index. Il montre comment le traitement affecte les ménages les plus pauvres (quantiles bas), la classe moyenne et les ménages les plus aisés (quantiles élevés).

# Testing Quantile treatment effects --------------------------------------
library(qte)

# Avec CiC -------------------

## Traitement 2008 -> 2021
set.seed(123)
dat_2per <- dat %>%
  filter(DHSYEAR %in% c(2008, 2021)) %>%
  mutate(treat = as.integer(GROUP == "Treatment"))

cic_res <- suppressWarnings(CiC(
  formla = wealth_centile_rural_weighted ~ treat,
  t = 2021, tmin1 = 2008, tname = "DHSYEAR",
  data = dat_2per,
  panel = FALSE, # repeated cross-sections
  se = TRUE, iters = 200, # bootstrap
  probs = seq(0.05, 0.95, 0.05),
  xformla = ~ spei_wc_n_2 + spei_wc_n_1 + spei_wc_n + hv219 + hv220
))

summary(cic_res)

Quantile Treatment Effect:
        
tau QTE Std. Error
0.05     -1.43    1.88
0.1  -9.21    2.62
0.15    -15.35    2.60
0.2 -15.55    2.33
0.25    -18.34    2.08
0.3 -18.61    2.23
0.35    -20.81    2.25
0.4 -21.16    2.29
0.45    -23.54    2.52
0.5 -23.94    2.34
0.55    -23.33    2.26
0.6 -22.56    2.17
0.65    -21.85    2.02
0.7 -19.46    1.97
0.75    -17.27    1.88
0.8 -14.72    1.77
0.85    -11.86    1.70
0.9  -8.05    1.71
0.95     -5.44    1.58

Average Treatment Effect:   -15.76
     Std. Error:        1.44

ggqte(cic_res) + labs(x="Quantiles", y="QTET", title="CiC QTET: 2008-2021")

## placebo

## Placebo: 1997 -> 2008
dat_placebo <- dat %>%
  filter(DHSYEAR %in% c(1997, 2008)) %>%
  mutate(treat = as.integer(GROUP == "Treatment"))

# (Optional) sanity check:
# with(dat_placebo, table(DHSYEAR, treat))

cic_pre <-  suppressWarnings(CiC(
  formla = wealth_centile_rural_simple ~ treat,
  t = 2008, tmin1 = 1997, tname = "DHSYEAR",
  data = dat_placebo,
  panel = FALSE, # repeated cross-sections
  se = TRUE, iters = 100, # bootstrap
  probs = seq(0.05, 0.95, 0.05),
  xformla = ~ spei_wc_n_2 + spei_wc_n_1 + spei_wc_n + hv219 + hv220
))

summary(cic_pre)

Quantile Treatment Effect:
        
tau QTE Std. Error
0.05    14.50    3.44
0.1 14.84    2.64
0.15    16.88    2.61
0.2 20.19    3.12
0.25    19.74    3.81
0.3 19.54    4.35
0.35    18.68    4.24
0.4 17.07    4.12
0.45    16.17    4.04
0.5 15.51    5.06
0.55    10.45    6.03
0.6  5.05    6.27
0.65     0.10    6.71
0.7 -4.42    6.07
0.75    -8.24    6.01
0.8 -7.50    5.64
0.85    -8.14    5.18
0.9 -6.39    4.50
0.95    -8.58    5.92

Average Treatment Effect:   6.69
     Std. Error:        3.37

ggqte(cic_pre) +
  labs(x = "Quantiles", y = "QTET",
       title = "Placebo CiC QTET: 1997-2008 - buffer 10 km")

Pour l’estimation principale (2008-2021), l’ATE est fortement négatif et très significatif (-17.9 centiles!). Les effets négatifs se retrouvent sur l’ensemble de la distribution, mais sont particulièrement marqués pour les centiles les plus faibles. Cela conforterait fortement l’hypothèse 2 (les AP renforcent les inégalités), mais c’est surprenant que ce soit à ce point marqué sur l’ensemble de la distribution alors que le DiD sur le wealth centile n’était pas significatif (Hypothèse 1).

Pour le test placebo, la courbe n’est pas exactement plate : il existe donc des différences significatives entre les groupes. Les ménages pauvres des zones traitées avaient déjà des niveaux de richesse différents de ceux des zones controles. Ces différences suggèrent une hétérogénéité entre les groupes.

9.3 Effect on inequalities

9.3.1 Staggered DiD

yvar <- "zscore_wealth"

# did2s (Gardner) : statique + event-study--------------------------------
# -- Statique (traitement "on/off")
did2s_static <- did2s(
  data        = dat3,
  yname       = yvar,
  first_stage = ~ spei_wc_n_2 + spei_wc_n_1 + spei_wc_n + hv219 + hv220 | DHSYEAR,
  second_stage= ~ i(treat_on, ref = FALSE),
  treatment   = "treat_on",
  cluster_var = "hv001",
  weights     = "w_all"
)
etable(did2s_static, headers = "did2s statique")

                   did2s_static
                 did2s statique
Dependent Var.:   zscore_wealth
                               
treat_on = 1    0.0015 (0.0057)
_______________ _______________
S.E.: Clustered       by: hv001
Observations             21,923
R2                     -4.43e-5
Adj. R2                -4.43e-5
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# -- Event-study (avec binning -5..5, ref = -1 et never-treated)
did2s_es <- did2s(
  data        = dat3,
  yname       = yvar,
  first_stage = ~ spei_wc_n_2 + spei_wc_n_1 + spei_wc_n + hv219 + hv220 | DHSYEAR,
  second_stage= ~ i(rel_year_binned, ref = c(-1, Inf)),
  treatment   = "treat_on",
  cluster_var = "hv001",
  weights     = "w_all"
)
etable(did2s_es, headers = "did2s event-study (-5..5)")

                                      did2s_es
                     did2s event-study (-5..5)
Dependent Var.:                  zscore_wealth
                                              
rel_year_binned = -5          0.0089* (0.0039)
rel_year_binned = -3           0.0040 (0.0064)
rel_year_binned = 0           -0.0034 (0.0116)
rel_year_binned = 2         0.0221*** (0.0067)
rel_year_binned = 3           0.0302* (0.0140)
rel_year_binned = 5           -0.0024 (0.0021)
____________________        __________________
S.E.: Clustered                      by: hv001
Observations                            21,923
R2                                     4.97e-5
Adj. R2                               -0.00018
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Plot ES did2s
plot_did2s <- broom::tidy(did2s_es, conf.int = TRUE) %>%
  filter(grepl("^rel_year_binned::", term)) %>%
  mutate(k = as.numeric(sub("^rel_year_binned::", "", term))) %>%
  arrange(k)

ggplot(plot_did2s, aes(k, estimate)) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_point() +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = .2) +
  labs(x = "Années relatives au 1er traitement (binnées -5..5)",
       y = "Effet DID sur le centile de richesse (pondéré)",
       title = "Event-study (did2s) - buffer 10 km") + 
  scale_x_continuous(breaks = c(-5,-4,-3,-2,-1,0,1,2,3,4,5))

# Test joint de pré-tendances (tous les leads k < 0)
# On construit un motif regex pour k = -5..-1 présents dans le modèle
leads_present <- plot_did2s$k[plot_did2s$k < 0]
if(length(leads_present) > 0){
  keep_regex <- paste0("^rel_year_binned::(", paste(leads_present, collapse="|"), ")$")
  print(fixest::wald(did2s_es, keep = keep_regex))
}

Wald test, H0: joint nullity of rel_year_binned::-5 and rel_year_binned::-3
 stat = 2.75235, p-value = 0.0638, on 2 and 21,917 DoF, VCOV: Corrected Clustered (hv001).$stat
[1] 2.752355

$p
[1] 0.06379954

$df1
[1] 2

$df2
[1] 21917

$vcov
[1] "Corrected Clustered (hv001)"

L’hypothèse des tendances parallèles est vérifié car l’hypothèse nulle n’est pas rejeté (p-value = 0.06379954 > 0.05). Après la mise en place de l’aire protégée, un effet positif apparaît deux ans après l’intervention (rel_year = 2: 0.0221***), puis devient particulièrement marqué au bout de trois ans (rel_year = 3: 0.0302*). Toutefois, cet effet disparaît au bout de cinq ans, indiquant que l’effet n’est pas durable. Le R² et R² ajusté montrent que la part de variance expliquée par le traitement est minimale.

9.4 Heterogeneity

L’étude des effets hétérogènes vise à déterminer quelles catégories de ménages bénéficient davantage – ou moins – de la création des aires protégées, et sous quelles conditions. Pour cela, nous analysons l’influence de la gouvernance en distinguant les aires protégées dites « strictes » (statuts IUCN I-IV) et les aires « multi-usages » (statuts IUCN V-VI).

# DID 2x2-----------------------------------------------------
dat2 <- dat %>%
  mutate(
    treat_strict = as.integer(GROUP == "Treatment" & IUCN_CAT %in% c("Ia","Ib",
                                                                     "II","III",
                                                                     "IV")),
    treat_multi  = as.integer(GROUP == "Treatment" & IUCN_CAT %in% c("V","VI"))
  )


# Placebo ---------------------------------------------------
pre <- dat2 %>%
  filter(DHSYEAR %in% c(1997, 2008)) %>%
  mutate(post = as.integer(DHSYEAR == 2008))

# Strict
f_pre_strict <- as.formula(paste(
  yvar, "~ treat_strict + post + treat_strict:post +",
  "spei_wc_n_2 + spei_wc_n_1 + spei_wc_n + hv219 + hv220"
))

m_pre_strict <- feols(f_pre_strict, data = pre, weights = ~ w_all, cluster = ~ hv001)

# Multi
f_pre_multi <- as.formula(paste(
  yvar, "~ treat_multi + post + treat_multi:post +",
  "spei_wc_n_2 + spei_wc_n_1 + spei_wc_n + hv219 + hv220"
))

m_pre_multi <- feols(f_pre_multi, data = pre, weights = ~ w_all, cluster = ~ hv001)

# Main estimation ---------------------------

main <- dat2 %>%
  filter(DHSYEAR %in% c(2008, 2021)) %>%
  mutate(post = as.integer(DHSYEAR == 2021))

m_main_strict <- feols(f_pre_strict, data = main, weights = ~ w_all, cluster = ~ hv001)
m_main_multi  <- feols(f_pre_multi, data = main, weights = ~ w_all, cluster = ~ hv001)


did_row <- function(model, year_post, term){
  summary(model) %>% broom::tidy() %>%
    filter(term == !!term) %>%
    transmute(year = year_post, estimate, se = std.error)
}

did_df <- bind_rows(
  did_row(m_pre_strict, 2008, "treat_strict:post") %>% mutate(type = "strict"),
  did_row(m_main_strict, 2021, "treat_strict:post") %>% mutate(type = "strict"),
  did_row(m_pre_multi, 2008, "treat_multi:post") %>% mutate(type = "multi"),
  did_row(m_main_multi, 2021, "treat_multi:post") %>% mutate(type = "multi")
) %>%
  mutate(
    period = factor(ifelse(year == 2008, "1997–2008", "2008–2021"),
                    levels = c("1997–2008","2008–2021")),
    lo = estimate - 1.96*se,
    hi = estimate + 1.96*se
  )


pd <- position_dodge(width = 0.2)

ggplot(did_df, aes(x = period, y = estimate, colour = type)) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_point(size = 3, position = pd) +
  geom_errorbar(aes(ymin = lo, ymax = hi), width = .2, position = pd) +
  labs(x = NULL, y = "Effet DID (centile de richesse)",
       title = "DID 2×2 par type d’AP (strict vs usage multiple)")

# did2s (Gardner) : statique + event-study--------------------------------
# Création du groupe IUCN 
dat3 <- dat3 %>%
  filter(!is.na(IUCN_CAT), !is.na(treat_on)) %>%
  mutate(
    IUCN_group = case_when(
    IUCN_CAT %in% c("Ia", "Ib", "II", "III", "IV") ~ "strict",
    IUCN_CAT %in% c("V", "VI") ~ "usage_multiple", 
    TRUE ~ NA_character_
  ),
  IUCN_group = factor(IUCN_group, levels = c("strict", "usage_multiple")),
  rel_year_binned = factor(rel_year_binned, levels = -5:5)
  ) %>%
  filter(!is.na(IUCN_group))

run_did2s_es <- function(data, yvar) {
  
  # DID Statique
  did2s_static <- did2s(
    data        = data,
    yname       = yvar,
    first_stage = ~ spei_wc_n_2 + spei_wc_n_1 + spei_wc_n + hv219 + hv220 | DHSYEAR,
    second_stage = ~ treat_on * IUCN_group,
    treatment   = "treat_on",
    cluster_var = "hv001",
    weights     = "w_all"
  )
  print(etable(did2s_static, headers = paste("did2s statique -", yvar)))
  
# DID event study
  
  did2s_es <- did2s(
    data = data,
    yname = yvar, 
    first_stage = ~ spei_wc_n_2 + spei_wc_n_1 + spei_wc_n + hv219 + hv220 | DHSYEAR,
    second_stage = ~ i(rel_year_binned, IUCN_group, ref = -1),
    treatment   = "treat_on",
    cluster_var = "hv001",
    weights     = "w_all"
  )
  
  print(etable(did2s_es, 
               headers = paste("did2s event-study (-5..5) par statut IUCN - buffer 10km", yvar)))
  
    # Extraction des coefficients
  tidy_es <- broom::tidy(did2s_es, conf.int = TRUE) %>%
    filter(grepl("^rel_year_binned::", term)) %>%
    mutate(
      year = as.numeric(stringr::str_extract(term, "(?<=::)-?[0-9]+")),
      
      # identifier strict vs usage_multiple
      group = case_when(
        grepl("usage_multiple", term, ignore.case = TRUE) ~ "usage_multiple",
        grepl("strict", term, ignore.case = TRUE) ~ "strict",
        TRUE ~ NA_character_
      ),
      outcome = yvar
    ) %>%
    filter(!is.na(group))

  #Plot
  plot_title <- paste0("Event-study (did2s) -", yvar, "\nHétérogénéité par statut IUCN - buffer 10km")
  
  g <- ggplot(tidy_es, aes(x = year, y = estimate, color = group)) +
    geom_hline(yintercept = 0, linetype = "dashed") +
    geom_point(size = 2) +
    geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2) +
    labs(
      x = "Années relatives au 1er traitement (binnées -5..5) - buffer 10 km",
      y = "Effet DID (centile de richesse)",
      title = plot_title,
      color = "Catégorie IUCN"
    ) +
    theme_minimal(base_size = 14)+ 
  scale_x_continuous(breaks = c(-5,-4,-3,-2,-1,0,1,2,3,4,5))
  
  print(g)
  
  # Test de pré-tendances
  leads <- tidy_es %>% filter(year < 0)
  
  if(nrow(leads) > 0){
    keep_regex <- paste0("^rel_year_binned::(", paste(unique(leads$year), collapse="|"), "):")
    print(wald(did2s_es, keep = keep_regex))
  }
  
  return(list(
    static = did2s_static,
    es     = did2s_es,
    coef   = tidy_es,
    plot   = g
  ))
}


res_wealth <- run_did2s_es(dat3, "wealth_centile_rural_weighted")

                                                                      did2s_static
                                    did2s statique - wealth_centile_rural_weighted
Dependent Var.:                                      wealth_centile_rural_weighted
                                                                                  
treat_on                                                           -0.5813 (9.508)
IUCN_groupstrict                                                     3.691 (3.061)
IUCN_groupusage_multiple                                            -1.348 (1.089)
treat_on x IUCN_groupusage_multiple                                 -1.556 (12.54)
___________________________________                  _____________________________
S.E.: Clustered                                                          by: hv001
Observations                                                                 4,889
R2                                                                         0.00745
Adj. R2                                                                    0.00684
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
                                                                                                                                did2s_es
                                                   did2s event-study (-5..5) par statut IUCN - buffer 10km wealth_centile_rural_weighted
Dependent Var.:                                                                                            wealth_centile_rural_weighted
                                                                                                                                        
rel_year_binned = -5 x IUCN_group = strict                                                                                 3.645 (2.687)
rel_year_binned = -5 x IUCN_group = usage_multiple                                                                        -1.591 (1.274)
rel_year_binned = -3 x IUCN_group = strict                                                                                 3.934 (5.255)
rel_year_binned = -3 x IUCN_group = usage_multiple                                                                        -1.384 (1.278)
rel_year_binned = 0 x IUCN_group = strict                                                                                -0.3570 (10.73)
rel_year_binned = 0 x IUCN_group = usage_multiple                                                                         9.464. (5.139)
rel_year_binned = 2 x IUCN_group = strict                                                                                  7.983 (5.681)
rel_year_binned = 2 x IUCN_group = usage_multiple                                                                      -10.96*** (2.898)
________________________________________                                                                   _____________________________
S.E.: Clustered                                                                                                                by: hv001
Observations                                                                                                                       4,889
R2                                                                                                                               0.01202
Adj. R2                                                                                                                          0.01060
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Wald test, H0: joint nullity of rel_year_binned::-5:IUCN_group::strict, rel_year_binned::-5:IUCN_group::usage_multiple, rel_year_binned::-3:IUCN_group::strict and rel_year_binned::-3:IUCN_group::usage_multiple
 stat = 0.756033, p-value = 0.553848, on 4 and 4,881 DoF, VCOV: Corrected Clustered (hv001).$stat
[1] 0.7560327

$p
[1] 0.5538482

$df1
[1] 4

$df2
[1] 4881

$vcov
[1] "Corrected Clustered (hv001)"

L’estimation du DID statique montre que, pour la période placebo, les aires protégées à usage multiple présentent un effet positif d’environ cinq centiles de richesse, tandis que les aires protégées strictes affichent un effet négatif. Dans les deux cas, les intervalles de confiance chevauchent zéro, indiquant que les effets ne sont pas significatifs. Pour la période de traitement (2008-2021), les effets estimés pour les aires protégées à usage multiple sont proches de zéro et légèrement négatif pour les aires protégées à usage stricte.

L’analyse en event study confirme dans l’ensemble l’absence d’effet significatif avant le traitement, qui valide l’hypothèse de tendances parallèles. A l’année 0, les estimations indiquent un effet moyen positif pour les aires protégées à usage multiple, tandis que les aires à usage strictes présentent une tendance plutôt négative. Deux ans après le traitement (rel_year_binned = 2), la situation s’inverse. Les ménages proches d’une aire protégée à usage multiple connaissent une diminution significative de leur centile de richesse (coefficient = -10.96*) mais ceux près des aires protégées à usage stricte connaissent une tendance positive, mais non significative (coefficient = 7.983).

Abadie, Alberto. 2021. “Using Synthetic Controls: Feasibility, Data Requirements, and Methodological Aspects.” Journal of Economic Literature 59 (2): 391425.

Borusyak, Kirill, Xavier Jaravel, and Jann Spiess. 2024. “Revisiting Event-Study Designs: Robust and Efficient Estimation.” The Review of Economic Studies, February, rdae007. https://doi.org/10.1093/restud/rdae007.