Problemset 02-AGE3802-1

Variable	Definición
countyid	Número identificador de la comuna
pre_acres	% de la superficie comunal bajo conservación en 1999
income	Mediana del ingreso a nivel de hogar para el año 1999
popdensity	Densidad poblacional de la comuna en 1999
consprog	Variable tipo dummy indicando si la comuna participó en el programa (0 indica no-participante)
post_acres	% de la superficie comunal bajo conservación en 2002

1. Muestre si la participación en el programa de conservación fue asignada aleatoriamente.

Respuesta: Para que exista aleatoriedad: \[T \perp \{Y(0),Y(1)\}\] Efectuamos test de medias con el indicador t-student para saber si hay medias diferentes entre tratados y contrroles, por cada variable:

#Abro el dataset
landuse=readxl::read_xlsx("data/input/landuse.xlsx") 

#Quito notación científica
options(scipen=9999) 

#Aplico t-student test para ver si hay aleatorización

#Test de medias en las hectáreas antes del tratamiento entre quienes se aplicó y no el programa.
t.test(landuse$pre_acres[landuse$consprog==1],landuse$pre_acres[landuse$consprog==0])

## 
##  Welch Two Sample t-test
## 
## data:  landuse$pre_acres[landuse$consprog == 1] and landuse$pre_acres[landuse$consprog == 0]
## t = -12.572, df = 439.46, p-value < 0.00000000000000022
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.05287259 -0.03857653
## sample estimates:
## mean of x mean of y 
## 0.1558641 0.2015886

#Test de medias en los ingresos entre quienes se aplicó y no el programa.
t.test(landuse$income[landuse$consprog==1],landuse$income[landuse$consprog==0])

## 
##  Welch Two Sample t-test
## 
## data:  landuse$income[landuse$consprog == 1] and landuse$income[landuse$consprog == 0]
## t = -7.7269, df = 430.89, p-value = 0.00000000000007816
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -2710.033 -1610.922
## sample estimates:
## mean of x mean of y 
##  10020.22  12180.69

#Test de medias en la densidad poblacional entre quienes se aplicó y no el programa.
t.test(landuse$popdensity[landuse$consprog==1],landuse$popdensity[landuse$consprog==0])

## 
##  Welch Two Sample t-test
## 
## data:  landuse$popdensity[landuse$consprog == 1] and landuse$popdensity[landuse$consprog == 0]
## t = 12.634, df = 366.76, p-value < 0.00000000000000022
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  31.61011 43.26427
## sample estimates:
## mean of x mean of y 
##  138.7773  101.3401

Todas las variables presentan diferencias signicantes entre los grupos tratados y no. Por lo tanto el experimento no fue aleatorizado.

2. Muestre si la participación en el programa de conservación fue tan buena como si hubiera sido realizada en forma aleatoria condicional sobre el ingreso y la densidad poblacional.

Viendo si el tratamiento es independiente de los controles \[T \perp \{income, popdensity \} \]

summary(lm(consprog ~income + popdensity, data=landuse))

## 
## Call:
## lm(formula = consprog ~ income + popdensity, data = landuse)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.77275 -0.32735 -0.04402  0.30170  1.02580 
## 
## Coefficients:
##                 Estimate   Std. Error t value             Pr(>|t|)    
## (Intercept)  0.120068369  0.091251629   1.316                0.189    
## income      -0.000042277  0.000005549  -7.619    0.000000000000131 ***
## popdensity   0.006519923  0.000498287  13.085 < 0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4006 on 497 degrees of freedom
## Multiple R-squared:  0.3354, Adjusted R-squared:  0.3328 
## F-statistic: 125.4 on 2 and 497 DF,  p-value: < 0.00000000000000022

Nos podemos dar cuenta que no son independientes.Pero deberíamos comprobar que con tratados y controles pre_acres no difere.

Sin controles

summary(lm(pre_acres ~consprog, data=landuse))

## 
## Call:
## lm(formula = pre_acres ~ consprog, data = landuse)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.109352 -0.025725  0.002483  0.026247  0.133926 
## 
## Coefficients:
##              Estimate Std. Error t value            Pr(>|t|)    
## (Intercept)  0.201589   0.002321   86.86 <0.0000000000000002 ***
## consprog    -0.045725   0.003670  -12.46 <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0402 on 498 degrees of freedom
## Multiple R-squared:  0.2377, Adjusted R-squared:  0.2361 
## F-statistic: 155.3 on 1 and 498 DF,  p-value: < 0.00000000000000022

Con controles

summary(lm(pre_acres ~consprog+ income + popdensity, data=landuse))

## 
## Call:
## lm(formula = pre_acres ~ consprog + income + popdensity, data = landuse)
## 
## Residuals:
##              Min               1Q           Median               3Q 
## -0.0000000195157 -0.0000000041659 -0.0000000001527  0.0000000039940 
##              Max 
##  0.0000000246308 
## 
## Coefficients:
##                         Estimate           Std. Error       t value
## (Intercept)  0.09999999966026010  0.00000000144197112  69349516.042
## consprog     0.00000000004588635  0.00000000070759152         0.065
## income       0.00001250000001812  0.00000000000009251 135122280.905
## popdensity  -0.00049999999477121  0.00000000000911419 -54859496.687
##                        Pr(>|t|)    
## (Intercept) <0.0000000000000002 ***
## consprog                  0.948    
## income      <0.0000000000000002 ***
## popdensity  <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.000000006319 on 496 degrees of freedom
## Multiple R-squared:      1,  Adjusted R-squared:      1 
## F-statistic: 8.813e+15 on 3 and 496 DF,  p-value: < 0.00000000000000022

La regresión anterior nos muestra que la diferencia entre tratamiento y controles se vuelve no significativo en la cantidad de acres antes del tratamiento, si consideramos otras variables explicativas como popdensity y income.

Ahora veamos el caso en los outcome de post: Caso sin controles

summary(lm(post_acres~consprog, data=landuse))

## 
## Call:
## lm(formula = post_acres ~ consprog, data = landuse)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.199839 -0.040789  0.000843  0.044065  0.180753 
## 
## Coefficients:
##             Estimate Std. Error t value            Pr(>|t|)    
## (Intercept) 0.200701   0.003628  55.317 <0.0000000000000002 ***
## consprog    0.004317   0.005737   0.753               0.452    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.06284 on 498 degrees of freedom
## Multiple R-squared:  0.001136,   Adjusted R-squared:  -0.0008697 
## F-statistic: 0.5664 on 1 and 498 DF,  p-value: 0.452

Estimador del tratamiento: 0.004317 (No significante)
Error residual estándar: 0.06284

Caso con controles

summary(lm(post_acres~consprog+income+popdensity, data=landuse))

## 
## Call:
## lm(formula = post_acres ~ consprog + income + popdensity, data = landuse)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.14712 -0.03013  0.00008  0.03348  0.17049 
## 
## Coefficients:
##                  Estimate    Std. Error t value             Pr(>|t|)    
## (Intercept)  0.0900095984  0.0112966463   7.968  0.00000000000001116 ***
## consprog     0.0456409756  0.0055433920   8.233  0.00000000000000162 ***
## income       0.0000123441  0.0000007247  17.033 < 0.0000000000000002 ***
## popdensity  -0.0003914407  0.0000714021  -5.482  0.00000006698954154 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0495 on 496 degrees of freedom
## Multiple R-squared:  0.3827, Adjusted R-squared:  0.3789 
## F-statistic: 102.5 on 3 and 496 DF,  p-value: < 0.00000000000000022

Estimador del tratamiento: 0.0456409756 (Muy significante)
Error residual estándar: 0.0495

Respuesta Las dos últimas regresiones muestran que sin considerar los controles, el tratamiento no es significativo para post_acres, pero al incluir estas variables explicativas la participación en el programa se comporta como si fuera asignada aleatoriamente.

3. Use regresión lineal para obtener un estimador válido del ATE y su error estándar. ¿Cuál es la estrategia de identificación que permite darle una interpretación causal al parámetro estimado y que está asociado a la variable de participación en el programa?

\[Y(0),Y(1)\perp consprog|(income, popdensity) \]

#Si incluye pre_acres: Resultado raro, mejor ignorar
#summary(lm(post_acres ~ consprog + pre_acres + income + popdensity, data = landuse))
summary(lm(post_acres ~ consprog + income + popdensity, data = landuse))

## 
## Call:
## lm(formula = post_acres ~ consprog + income + popdensity, data = landuse)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.14712 -0.03013  0.00008  0.03348  0.17049 
## 
## Coefficients:
##                  Estimate    Std. Error t value             Pr(>|t|)    
## (Intercept)  0.0900095984  0.0112966463   7.968  0.00000000000001116 ***
## consprog     0.0456409756  0.0055433920   8.233  0.00000000000000162 ***
## income       0.0000123441  0.0000007247  17.033 < 0.0000000000000002 ***
## popdensity  -0.0003914407  0.0000714021  -5.482  0.00000006698954154 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0495 on 496 degrees of freedom
## Multiple R-squared:  0.3827, Adjusted R-squared:  0.3789 
## F-statistic: 102.5 on 3 and 496 DF,  p-value: < 0.00000000000000022

Respuesta: El coeficiente del ATE es \(0.0456409756\) y su error estándar \(0.0055433920\) y lo podemos interpretar que participar en el programa de conservación aumenta en promedio, 4.6 puntos porcentuales, condicionado en las variables de ingreso y densidad. Como ya se ha mencionado en las preguntas anteriores, estamos asumiendo que la variable tratamiento consprogno está correlacionada con el error de la ecuación, o sea:

\[ \mathbb{E}[\epsilon|consprog,income, popdensity]=0\]

4. Use teffect en STATA ó matching en R para estimar ATE y ATT usando matching por covariables con el vecino más cercano y los 4 vecinos más cercanos. ¿De qué forma sus respuestas se comparan con las otras estimaciones obtenidas anteriormente?

library(Matching)

## Warning: package 'Matching' was built under R version 4.3.3

## Loading required package: MASS

## ## 
## ##  Matching (Version 4.10-15, Build Date: 2024-10-14)
## ##  See https://www.jsekhon.com for additional documentation.
## ##  Please cite software as:
## ##   Jasjeet S. Sekhon. 2011. ``Multivariate and Propensity Score Matching
## ##   Software with Automated Balance Optimization: The Matching package for R.''
## ##   Journal of Statistical Software, 42(7): 1-52. 
## ##

# Preparar datos
Y <- landuse$post_acres
Tr <- landuse$consprog
X <- cbind(landuse$income,landuse$popdensity)

# Matching 1-neighbour

# ATT 1-neighbour
att_1 <- Match(Y = Y,
               Tr = Tr,
               X = X,
               M = 1)   # 1 vecino más cercano

summary(att_1)

## 
## Estimate...  0.03727 
## AI SE......  0.010537 
## T-stat.....  3.537 
## p.val......  0.00040476 
## 
## Original number of observations..............  500 
## Original number of treated obs...............  200 
## Matched number of observations...............  200 
## Matched number of observations  (unweighted).  200

# ATE 1-neighbour
ate_1 <- Match(Y = Y,
               Tr = Tr,
               X = X,
               M = 1,
               estimand = "ATE")

summary(ate_1)

## 
## Estimate...  0.041902 
## AI SE......  0.0073977 
## T-stat.....  5.6642 
## p.val......  0.000000014771 
## 
## Original number of observations..............  500 
## Original number of treated obs...............  200 
## Matched number of observations...............  500 
## Matched number of observations  (unweighted).  500

# Matching 4-neighbours

# ATT 4-neighbours
att_4 <- Match(Y = Y,
               Tr = Tr,
               X = X,
               M = 4)   # 1 vecino más cercano

summary(att_4)

## 
## Estimate...  0.034662 
## AI SE......  0.0072371 
## T-stat.....  4.7896 
## p.val......  0.0000016715 
## 
## Original number of observations..............  500 
## Original number of treated obs...............  200 
## Matched number of observations...............  200 
## Matched number of observations  (unweighted).  800

# ATE 4-neighbours
ate_4 <- Match(Y = Y,
               Tr = Tr,
               X = X,
               M = 4,
               estimand = "ATE")

summary(ate_4)

## 
## Estimate...  0.040193 
## AI SE......  0.0061535 
## T-stat.....  6.5318 
## p.val......  0.000000000064996 
## 
## Original number of observations..............  500 
## Original number of treated obs...............  200 
## Matched number of observations...............  500 
## Matched number of observations  (unweighted).  2000

Tabla acumulada

RESULTADOS	1 VECINO	4 VECINOS	OLS CONTROL
ATT	0.03727	0.034662	-
ATE	0.041902	0.040193	0.0456409756

Matching corrige y crea controles comparables. Los resultados también muestran que el matching corrigió a la baja, o sea, el estimador ahora tiene un valor más bajo en el efecto del tratamiento. Además parece ser más consistente con OLS al usar matching para el ATE con 1 vecino. Por último, podemos ver la validez y robustés del tratamiento ya que los coeficientes asociados al tratamiento consprog son cercanos en cuantoa sus valores.

5. Estime el ATE y ATT usando propensity score matching y siguiendo el siguiente procedimiento:

a. Use una función probit para estimar la participación en el programa de conservación como una función que depende del ingreso y la densidad poblacional (esto corresponde al propensity score).

modelo1=glm(consprog~income + popdensity,family=binomial(link="probit"),data=landuse)
summary(modelo1)

## 
## Call:
## glm(formula = consprog ~ income + popdensity, family = binomial(link = "probit"), 
##     data = landuse)
## 
## Coefficients:
##                Estimate  Std. Error z value             Pr(>|z|)    
## (Intercept) -1.27945477  0.33310078  -3.841             0.000123 ***
## income      -0.00015702  0.00002256  -6.960      0.0000000000034 ***
## popdensity   0.02329757  0.00218599  10.658 < 0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 673.01  on 499  degrees of freedom
## Residual deviance: 472.63  on 497  degrees of freedom
## AIC: 478.63
## 
## Number of Fisher Scoring iterations: 5

Respuesta:El Propensity Score muestra que la participación al programa depende de las variables ingresos y densidad poblacional.

b. Use el comando predict en Stata o su análogo en R para generar las predicciones de las probabilidades de participación en el programa para cada observación.

landuse$pscore <- predict(modelo1, type = "response")

c. Use el comando teffect en STATA y efectúe el matching usando el propensity score con el vecino más cercano. ¿De qué forma sus respuestas se comparan con las otras estimaciones obtenidas anteriormente? ¿Cuál es ahora la estrategia de identificación que permite estimar efectos causales consistentes y no-sesgados?

# ATT con 1 vecino (estimador clásico)
att_1 <- Match(Y = landuse$post_acres,
               Tr = landuse$consprog,
               X = landuse$pscore,
               M = 1)

summary(att_1)

## 
## Estimate...  0.05373 
## AI SE......  0.014199 
## T-stat.....  3.784 
## p.val......  0.00015432 
## 
## Original number of observations..............  500 
## Original number of treated obs...............  200 
## Matched number of observations...............  200 
## Matched number of observations  (unweighted).  229

# ATT con 1 vecino (estimador clásico)
ate_1 <- Match(Y = landuse$post_acres,
               Tr = landuse$consprog,
               X = landuse$pscore,
               M = 1,
               estimand = "ATE",
               Weight = 2)

summary(ate_1)

## 
## Estimate...  0.047035 
## AI SE......  0.0091635 
## T-stat.....  5.1328 
## p.val......  0.00000028543 
## 
## Original number of observations..............  500 
## Original number of treated obs...............  200 
## Matched number of observations...............  500 
## Matched number of observations  (unweighted).  566

Tabla acumulada

RESULTADOS	1 VECINO	4 VECINOS	OLS CONTROL	1 VECINO PS
ATT	0.03727	0.034662	-	0.05373
ATE	0.041902	0.040193	0.0456409756	0.047035

Respuesta

Bueno, usamos tres métodos, matching tradicional, OLS, y Propensity Score y todos los coeficientes fueron cercanos respecto a la variable tratamiento consprog tal como muestra la tabla anterior, entonces podemos confirmar el ATE obtenido inicialmente.

Por otro lado, en este paso estamos usando la siguiente estrategia de identificación:

\[Y(1),Y(0)\perp consprog|pscore\] O sea, que podemos ignorar el tratamiento dado el propensity score que captura el efecto de las covariables sobre el tratamiento, tal como se ve en el probit de esta etapa, por lo que podemos decir que la asignación al programa fue aleatoria, y que pudimos encontrar un efecto significante del tratamiento.

6. Presente una representación gráfica de las relaciones causales que permiten justificar el uso de matching como método válido para estimar el efecto causal del programa de conservación.

Respuesta

En nuestro caso podriamos esperar que las covariables income y popdensity están afectando tanto al tratamiento consprog como a la variable dependiente post_acres, o sea, que el tratamiento no fue asignado aleatoriamente. A continuación, un ejemplo:

library(dagitty)

g <- dagitty("
dag {
  income -> consprog
  income -> post_acres
  popdensity -> consprog
  popdensity -> post_acres
  consprog -> post_acres
}
")

plot(g)

## Plot coordinates for graph not supplied! Generating coordinates, see ?coordinates for how to set your own.

Por ende, el matching sería una excelente estrategia de bloquear estos “backdoor paths” ya que controles que sean más homogeneos a donde se aplico al tratamiento, condicionando en dos variables significativas incomey popdensity

Por último, presento la distribución del Propensity Score obtenido del modelo probit que pone al tratamiento como variable dependiente. Hay un traslape entre las distribuciones con y sin tratamiento, lo que permite apoyar la idea de que podemos encontrar comunas no tratadas con valores similares.

library(ggplot2)

## Warning: package 'ggplot2' was built under R version 4.3.3

ggplot(landuse, aes(x = pscore, fill = factor(consprog))) +
  geom_density(alpha = 0.4) +
  labs(
    x = "Propensity score",
    y = "Densidad",
    fill = "Tratamiento",
    title = "Distribución del propensity score por estado de tratamiento"
  ) +
  theme_minimal()