library(here)
library(readxl)
library(dplyr)
library(MatchIt)
library(cobalt)
library(rbounds)
library(ggplot2)
library(sandwich)
library(lmtest)
library(car)
library(corrplot)
library(sf)
library(stringr)
library(broom)
library(kableExtra)
I’m filtering out localities that are not in forested biomes and, for this analysis, using the localities with Ejidos (no matter the area of the thiessen polygon covered by Ejido polygon) as the treated units and rural localities (withou ejido or comunidades) as controls
I matched on urbanization levels, agricultural area, area under PA, total population, slope (I decided to use slope instead of altitute because this is usually more common), distance to 1M people cities, and latitude and longitude. Also I did an exact matching at the biome level only. I allowed replacement for observations and trimmed observations with too high or too low propensity score to enhance common support. I matched 1 control to 1 treated unit
#propensity score matching----
##2010-2020----
matchit(formula = treat_ejido ~
perc_urban_10 +
perc_agriculture_10 +
perc_area_PA_fed_strict_10 +
perc_area_PA_fed_sustainable_10 +
perc_area_PA_otras_strict_10 +
perc_area_PA_otras_sustainable_10 +
pop_tot_10 +
slope +
dist_1m +
lat_dec +
long_dec,
data = table_matching_forest,
method = "nearest",
distance = "glm",
link = "logit",
#distance.options = list(), #this argument is used to supply the arguments to the distance function. I'm using a glm with logit to calculat propability propensity scores
estimand = "ATT",
exact = ~ biome,
mahvars = NULL, #I don't want to calculate mahalanobis distance (and create calipers) for specific variables
antiexact = NULL, #I don't need anti-exact matching for this
discard = "both", #I'm not discarding (trimming) observations in control and treatment groups to enhance 'common support'
reestimate = TRUE, #as not I'm trimming the ps observations, I will reestimate the ps
s.weights = NULL, #I don't have any prior expectations to put weight in the observations
replace = TRUE,
m.order = NULL, #only important if you are doing the matching without replacement
caliper = NULL, #there is no prior belief that our matching should be more restrict in one or other variable, thus I ignored the argument
#std.caliper = TRUE,
ratio = 1,
verbose = TRUE,
include.obj = FALSE) -> m_ejido_forest_psm20
match.data(m_ejido_forest_psm20) -> m_data_ejido_forest_psm20
##
## Call:
## matchit(formula = treat_ejido ~ perc_urban_10 + perc_agriculture_10 +
## perc_area_PA_fed_strict_10 + perc_area_PA_fed_sustainable_10 +
## perc_area_PA_otras_strict_10 + perc_area_PA_otras_sustainable_10 +
## pop_tot_10 + slope + dist_1m + lat_dec + long_dec, data = table_matching_forest,
## method = "nearest", distance = "glm", link = "logit", estimand = "ATT",
## exact = ~biome, mahvars = NULL, antiexact = NULL, discard = "both",
## reestimate = TRUE, s.weights = NULL, replace = TRUE, m.order = NULL,
## caliper = NULL, ratio = 1, verbose = TRUE, include.obj = FALSE)
##
## Summary of Balance for All Data:
## Means Treated Means Control
## distance 0.6946 0.6601
## perc_urban_10 1.9457 1.7877
## perc_agriculture_10 58.6675 58.9155
## perc_area_PA_fed_strict_10 0.3356 0.5083
## perc_area_PA_fed_sustainable_10 3.6514 5.2726
## perc_area_PA_otras_strict_10 0.1528 0.1030
## perc_area_PA_otras_sustainable_10 0.8998 0.6454
## pop_tot_10 151.5919 93.2923
## slope 9.9084 11.2187
## dist_1m 351.2062 316.6243
## lat_dec 20.6028 19.9140
## long_dec -99.3119 -98.5592
## biomeMangroves 0.0139 0.0233
## biomeMediterranean Forests, Woodlands & Scrub 0.0105 0.0159
## biomeTropical & Subtropical Coniferous Forests 0.2652 0.2812
## biomeTropical & Subtropical Dry Broadleaf Forests 0.3922 0.3122
## biomeTropical & Subtropical Moist Broadleaf Forests 0.3182 0.3674
## Std. Mean Diff. Var. Ratio
## distance 0.4131 1.0925
## perc_urban_10 0.0160 0.9832
## perc_agriculture_10 -0.0065 0.9239
## perc_area_PA_fed_strict_10 -0.0328 0.6611
## perc_area_PA_fed_sustainable_10 -0.0896 0.6612
## perc_area_PA_otras_strict_10 0.0151 1.5008
## perc_area_PA_otras_sustainable_10 0.0300 1.3173
## pop_tot_10 0.1348 1.4477
## slope -0.1885 0.8252
## dist_1m 0.1342 1.2062
## lat_dec 0.1954 1.2749
## long_dec -0.1281 1.2801
## biomeMangroves -0.0801 .
## biomeMediterranean Forests, Woodlands & Scrub -0.0531 .
## biomeTropical & Subtropical Coniferous Forests -0.0363 .
## biomeTropical & Subtropical Dry Broadleaf Forests 0.1639 .
## biomeTropical & Subtropical Moist Broadleaf Forests -0.1056 .
## eCDF Mean eCDF Max
## distance 0.1147 0.1727
## perc_urban_10 0.0144 0.0435
## perc_agriculture_10 0.0227 0.0545
## perc_area_PA_fed_strict_10 0.0020 0.0035
## perc_area_PA_fed_sustainable_10 0.0142 0.0159
## perc_area_PA_otras_strict_10 0.0006 0.0010
## perc_area_PA_otras_sustainable_10 0.0033 0.0062
## pop_tot_10 0.0233 0.0953
## slope 0.0471 0.0789
## dist_1m 0.0256 0.0824
## lat_dec 0.0471 0.1076
## long_dec 0.0479 0.1169
## biomeMangroves 0.0094 0.0094
## biomeMediterranean Forests, Woodlands & Scrub 0.0054 0.0054
## biomeTropical & Subtropical Coniferous Forests 0.0160 0.0160
## biomeTropical & Subtropical Dry Broadleaf Forests 0.0800 0.0800
## biomeTropical & Subtropical Moist Broadleaf Forests 0.0492 0.0492
##
## Summary of Balance for Matched Data:
## Means Treated Means Control
## distance 0.6946 0.6946
## perc_urban_10 1.9445 1.8430
## perc_agriculture_10 58.6679 58.2491
## perc_area_PA_fed_strict_10 0.3356 0.3746
## perc_area_PA_fed_sustainable_10 3.6515 3.4190
## perc_area_PA_otras_strict_10 0.1528 0.1450
## perc_area_PA_otras_sustainable_10 0.8999 1.0493
## pop_tot_10 150.6054 136.3017
## slope 9.9085 9.7640
## dist_1m 351.2059 359.9375
## lat_dec 20.6028 20.5509
## long_dec -99.3120 -99.2604
## biomeMangroves 0.0139 0.0139
## biomeMediterranean Forests, Woodlands & Scrub 0.0105 0.0105
## biomeTropical & Subtropical Coniferous Forests 0.2652 0.2652
## biomeTropical & Subtropical Dry Broadleaf Forests 0.3922 0.3922
## biomeTropical & Subtropical Moist Broadleaf Forests 0.3182 0.3182
## Std. Mean Diff. Var. Ratio
## distance 0.0000 1.0000
## perc_urban_10 0.0103 1.0647
## perc_agriculture_10 0.0109 0.9154
## perc_area_PA_fed_strict_10 -0.0074 0.9158
## perc_area_PA_fed_sustainable_10 0.0128 1.0464
## perc_area_PA_otras_strict_10 0.0024 1.0196
## perc_area_PA_otras_sustainable_10 -0.0176 0.7861
## pop_tot_10 0.0331 0.8202
## slope 0.0208 1.0205
## dist_1m -0.0339 0.9194
## lat_dec 0.0147 1.0320
## long_dec -0.0088 0.9971
## biomeMangroves -0.0000 .
## biomeMediterranean Forests, Woodlands & Scrub -0.0000 .
## biomeTropical & Subtropical Coniferous Forests -0.0000 .
## biomeTropical & Subtropical Dry Broadleaf Forests -0.0000 .
## biomeTropical & Subtropical Moist Broadleaf Forests -0.0000 .
## eCDF Mean eCDF Max
## distance 0.0000 0.0004
## perc_urban_10 0.0088 0.0298
## perc_agriculture_10 0.0222 0.0510
## perc_area_PA_fed_strict_10 0.0006 0.0016
## perc_area_PA_fed_sustainable_10 0.0023 0.0040
## perc_area_PA_otras_strict_10 0.0002 0.0005
## perc_area_PA_otras_sustainable_10 0.0015 0.0025
## pop_tot_10 0.0093 0.0809
## slope 0.0064 0.0176
## dist_1m 0.0097 0.0389
## lat_dec 0.0194 0.0600
## long_dec 0.0158 0.0479
## biomeMangroves 0.0000 0.0000
## biomeMediterranean Forests, Woodlands & Scrub 0.0000 0.0000
## biomeTropical & Subtropical Coniferous Forests 0.0000 0.0000
## biomeTropical & Subtropical Dry Broadleaf Forests 0.0000 0.0000
## biomeTropical & Subtropical Moist Broadleaf Forests 0.0000 0.0000
## Std. Pair Dist.
## distance 0.0005
## perc_urban_10 0.3520
## perc_agriculture_10 0.9851
## perc_area_PA_fed_strict_10 0.1319
## perc_area_PA_fed_sustainable_10 0.3352
## perc_area_PA_otras_strict_10 0.0884
## perc_area_PA_otras_sustainable_10 0.2232
## pop_tot_10 0.3787
## slope 0.7576
## dist_1m 0.8611
## lat_dec 0.7031
## long_dec 0.6912
## biomeMangroves 0.0000
## biomeMediterranean Forests, Woodlands & Scrub 0.0000
## biomeTropical & Subtropical Coniferous Forests 0.0000
## biomeTropical & Subtropical Dry Broadleaf Forests 0.0000
## biomeTropical & Subtropical Moist Broadleaf Forests 0.0000
##
## Sample Sizes:
## Control Treated
## All 40107. 86699
## Matched (ESS) 16079.14 86697
## Matched 28954. 86697
## Unmatched 11147. 0
## Discarded 6. 2
##
## Rosenbaum Sensitivity Test for Wilcoxon Signed Rank P-Value
##
## Unconfounded estimate .... 1
##
## Gamma Lower bound Upper bound
## 1 1 1
## 2 1 1
## 3 1 1
##
## Note: Gamma is Odds of Differential Assignment To
## Treatment Due to Unobserved Factors
##
I estimated the ATT without any regionalization, but included biome as a variable to control. The biome reference level is Mangrove
| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | -81.3614 | 2.6461 | -30.7480 | 0.0000 |
| treat_ejido | 1.4272 | 0.1959 | 7.2836 | 0.0000 |
| netMig_10_20 | 0.0000 | 0.0000 | 1.0828 | 0.2789 |
| biomeMediterranean Forests, Woodlands & Scrub | -35.9693 | 1.0086 | -35.6614 | 0.0000 |
| biomeTropical & Subtropical Coniferous Forests | 4.9179 | 0.6998 | 7.0273 | 0.0000 |
| biomeTropical & Subtropical Dry Broadleaf Forests | -2.2091 | 0.6897 | -3.2032 | 0.0014 |
| biomeTropical & Subtropical Moist Broadleaf Forests | -1.7533 | 0.6821 | -2.5704 | 0.0102 |
| perc_urban_10 | -0.4259 | 0.0057 | -74.6298 | 0.0000 |
| perc_agriculture_10 | -0.4332 | 0.0033 | -132.9308 | 0.0000 |
| perc_area_PA_fed_strict_10 | 0.0295 | 0.0150 | 1.9716 | 0.0487 |
| perc_area_PA_fed_sustainable_10 | 0.0373 | 0.0045 | 8.2926 | 0.0000 |
| perc_area_PA_otras_strict_10 | 0.0350 | 0.0167 | 2.0904 | 0.0366 |
| perc_area_PA_otras_sustainable_10 | -0.0438 | 0.0121 | -3.6126 | 0.0003 |
| pop_tot_10 | -0.0022 | 0.0002 | -10.0571 | 0.0000 |
| slope | 1.1654 | 0.0164 | 71.0681 | 0.0000 |
| dist_1m | 0.0113 | 0.0004 | 25.9341 | 0.0000 |
| lat_dec | -0.7054 | 0.0514 | -13.7297 | 0.0000 |
| long_dec | -1.3441 | 0.0346 | -38.8253 | 0.0000 |
| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 17.1466 | 0.0834 | 205.7098 | 0 |
| perc_area_communal | 0.4954 | 0.0021 | 232.6230 | 0 |
I guess what we want to do mostly is to show that the effect of ejidos on forest cover (I will do the analysis on poverty afterwards) is different in different parts of the country. That regionalization we discussed last time. The thing is, using any kind of regionalization we discuss so far (i.e. state, biome, economic regions…) has a too coarse spatial resolution and much better would be to show the variation at the locality (thiessen polygon) level. The approach I did before was wrong, and when I fixed it, I got one single value as you expected before. Any ideas on how I can show this variation, if that is possible at all, would be very helpful