-
Notifications
You must be signed in to change notification settings - Fork 0
/
05-assignation_aleatoire.qmd
188 lines (149 loc) · 8.76 KB
/
05-assignation_aleatoire.qmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
---
output: html_document
editor_options:
chunk_output_type: console
---
# Méthode randomisée
```{=html}
<iframe src="TVM_RCT.pdf" width="100%" height="400px"></iframe>
```
[Cliquer ici pour télécharger la présentation](TVM_RCT.pdf).
> **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 -@sec-deforestation]. 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.
```{r}
# 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")
```
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.
```{r}
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")
```
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} [@hlavac2022].
```{r}
#| output: false
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)
```
```{r}
#| output: asis
#| echo: false
cat(rct_out1)
```
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.
```{r}
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")
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.
```{r}
#| output: false
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)
```
```{r}
#| output: asis
#| echo: false
cat(rct_out2)
```
Apparemment, le traitement reste significatif une fois que l'on contrôle pour la distance aux villes.