0. Configuración de datos

Empezaremos configurando el entorno de R para responder las preguntas pertinentes.

#Cargar tidyverse
library(tidyverse)

#Cargar datos
data=read.csv("earthfriends.csv")

#Pivotear la tabla para hacer figuras y análisis
data_longer=pivot_longer(data,cols=c(pre_pollution,post_pollution),
                         names_to = "period",
                         values_to = "pollution")

#Cambiar los valores del periodo
#data_longer=mutate(data_longer, period = ifelse(period == "pre_pollution", 0, 1))

I. Explorando los datos de polución

data_longer$name=factor(data_longer$period,
       levels = c("pre_pollution", "post_pollution"),
       labels = c("Pre intervención", "Post intervención"))

ggplot(data_longer, aes(x=name,y=pollution,col=as.factor(protocol)))+
  geom_boxplot()+
  labs(color = "Protocolo",
       x = "Periodo",
       y = "Niveles de polución",
       title= "Box-plotiveles de polución")

ggplot(data_longer, aes(x=name,y=age,col=as.factor(protocol)))+
  geom_boxplot()+
  labs(color = "Protocolo",
       x = "Periodo",
       y = "edad",
       title = "Box-plot edad de las empresas")

Qué proporción cambió la contaminación después del tratamiento respecto al nivel inicial del grupo

library(knitr)

#Diferencias pre y post pollution
data$diff=data$post_pollution-data$pre_pollution

#Calculando la proporción de cambio
tratamiento=mean(data$diff[data$assignment==1&data$protocol==1])/mean(data$pre_pollution[data$assignment==1])
#Protocolo 1
protocolo1=mean(data$diff[data$assignment==1&data$protocol==1])/mean(data$pre_pollution[data$assignment==1&data$protocol==1])
#Protocolo 2
protocolo2=mean(data$diff[data$assignment==1&data$protocol==2])/mean(data$pre_pollution[data$assignment==1&data$protocol==2])
#Sin tratamiento Aumenta
control=mean(data$diff[data$assignment==0])/mean(data$pre_pollution[data$assignment==0])

prepospollution= data.frame(Etapa=c("Con tratamiento","Sin tratamiento", "Protocolo 1", "Protocolo 2"),
                            Cambios=c(tratamiento,control,protocolo1,protocolo2)
                            )


kable(prepospollution)
Etapa Cambios
Con tratamiento -0.3462160
Sin tratamiento 0.4713735
Protocolo 1 -0.2858306
Protocolo 2 -0.1459551

II. Desarrollo de tarea

1. Viendo si existen diferencias entre los grupos en términos de sus niveles de polución previos al experimento, edad y número de empleados.

Para esto paso usaremos t-test o un test de student

# PREGUNTA 1 ------------------
## PRE POLLUTION -------------
#T-test pre pollution
t.test(data$pre_pollution[data$assignment==0],data$pre_pollution[data$assignment==1])
## 
##  Welch Two Sample t-test
## 
## data:  data$pre_pollution[data$assignment == 0] and data$pre_pollution[data$assignment == 1]
## t = -1.3773, df = 38.785, p-value = 0.1763
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -2968.431   563.771
## sample estimates:
## mean of x mean of y 
##  3753.975  4956.305
## EDAD -------------
#T-test pre edad
t.test(data$age[data$assignment==0],data$age[data$assignment==1])
## 
##  Welch Two Sample t-test
## 
## data:  data$age[data$assignment == 0] and data$age[data$assignment == 1]
## t = -1.3773, df = 38.785, p-value = 0.1763
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -3.265274  0.620148
## sample estimates:
## mean of x mean of y 
##  5.129372  6.451936
## EMPLOYEES -------------
#T-test empleados
t.test(data$employees[data$assignment==0],data$employees[data$assignment==1])
## 
##  Welch Two Sample t-test
## 
## data:  data$employees[data$assignment == 0] and data$employees[data$assignment == 1]
## t = -1.1006, df = 40.965, p-value = 0.2775
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -58.05709  17.09895
## sample estimates:
## mean of x mean of y 
##  162.5126  182.9917

Ninguna de las variables parece ser significante al comparar sus medias. O sea, no hay diferencias significativas. Por otro lado, debemos considerar que los resultados son descriptivos y muestran cambios promedio dentro de cada grupo pero no permiten determinar causalidad entre las variables usadas.

2. Viendo si la asignación a los nuevos dos protocolos fue efectivamente aleatoria. Si existen diferencias explique qué podría explicar esta situación.

# PREGUNTA 2 ----------
## X:PRE-POLLUTION & CONTROL v/s Y:PRE-POLLUTION & PROTOCOLO 1 ------------
### Resultado: P-value <0.05)
t.test(data$pre_pollution[data$protocol==0&data$assignment==0],data$pre_pollution[data$protocol==1&data$assignment==1])
## 
##  Welch Two Sample t-test
## 
## data:  data$pre_pollution[data$protocol == 0 & data$assignment == 0] and data$pre_pollution[data$protocol == 1 & data$assignment == 1]
## t = -2.3285, df = 37.79, p-value = 0.02534
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -4205.4458  -293.3791
## sample estimates:
## mean of x mean of y 
##  3753.975  6003.387
## X:PRE-POLLUTION CONTROL v/s Y:PRE-POLLUTION PROTOCOLO 2 ------------
### No significante
t.test(data$pre_pollution[data$protocol==0&data$assignment==0],data$pre_pollution[data$protocol==2&data$assignment==1])
## 
##  Welch Two Sample t-test
## 
## data:  data$pre_pollution[data$protocol == 0 & data$assignment == 0] and data$pre_pollution[data$protocol == 2 & data$assignment == 1]
## t = -0.15339, df = 37.983, p-value = 0.8789
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -2204.207  1893.712
## sample estimates:
## mean of x mean of y 
##  3753.975  3909.223
## X:AGE & CONTROL v/s Y:AGE & PROTOCOLO 1 ------------
#T-test age, protocolo 1 (P-value <0.05)
t.test(data$age[data$protocol==0&data$assignment==0],data$age[data$protocol==1&data$assignment==1])
## 
##  Welch Two Sample t-test
## 
## data:  data$age[data$protocol == 0 & data$assignment == 0] and data$age[data$protocol == 1 & data$assignment == 1]
## t = -2.3285, df = 37.79, p-value = 0.02534
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -4.6259904 -0.3227172
## sample estimates:
## mean of x mean of y 
##  5.129372  7.603726
## X:AGE & CONTROL v/s Y:AGE & PROTOCOLO 2 ------------
### No significante
t.test(data$age[data$protocol==0&data$assignment==0],data$age[data$protocol==2&data$assignment==1])
## 
##  Welch Two Sample t-test
## 
## data:  data$age[data$protocol == 0 & data$assignment == 0] and data$age[data$protocol == 2 & data$assignment == 1]
## t = -0.15339, df = 37.983, p-value = 0.8789
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -2.424628  2.083083
## sample estimates:
## mean of x mean of y 
##  5.129372  5.300145
## X:EMPLOYEES & CONTROL v/s Y:EMPLOYEES & PROTOCOLO 1 ------------
### No significante
t.test(data$employees[data$protocol==0&data$assignment==0],data$employees[data$protocol==1&data$assignment==1])
## 
##  Welch Two Sample t-test
## 
## data:  data$employees[data$protocol == 0 & data$assignment == 0] and data$employees[data$protocol == 1 & data$assignment == 1]
## t = -0.15392, df = 37.701, p-value = 0.8785
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -47.68715  40.94943
## sample estimates:
## mean of x mean of y 
##  162.5126  165.8815
## X:EMPLOYEES & CONTROL v/s Y:EMPLOYEES & PROTOCOLO 2 ------------
### No significante
t.test(data$employees[data$protocol==0&data$assignment==0],data$employees[data$protocol==2&data$assignment==1])
## 
##  Welch Two Sample t-test
## 
## data:  data$employees[data$protocol == 0 & data$assignment == 0] and data$employees[data$protocol == 2 & data$assignment == 1]
## t = -1.7679, df = 37.954, p-value = 0.08512
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -80.634650   5.456115
## sample estimates:
## mean of x mean of y 
##  162.5126  200.1019

Hay dos diferencias significativas al comparar la edad del control con el protocolo 1. Por otro lado, la diferencia de los niveles de pre-polución entre el control y el protocolo 1 también son significativos, esto es reforzado en los boxplot de la sección I. Una explicación, puede ser que las empresas a las que se le aplicó el protocolo 1 sean porque no están correctamente balanceadas al azar respecto a agey pre_pollution, o sea, que en promedio el protocolo 1 se aplicó a las empresas con más años y que contaminaban más que el resto. Aparentemente, parece ser que este protocolo es el idóneo para las empresas que promediaban más años, porque una vez aplicado el tratamiento, diferencian un decrecimiento de casi el doble respecto a las que se les aplicó el protocolo 2.

3. Calcualando

\[\mathbb{E}(Y)\]

options(scipen = 999999)
weighted.mean(data$post_pollution)
## [1] 4383.195
t.test(data$post_pollution)
## 
##  One Sample t-test
## 
## data:  data$post_pollution
## t = 12.397, df = 59, p-value < 0.00000000000000022
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  3675.697 5090.694
## sample estimates:
## mean of x 
##  4383.195

\[\mathbb{E}(Y|T=1)\]

options(scipen = 999999)

weighted.mean(data$post_pollution[data$assignment==1])
## [1] 3813.044
t.test(data$post_pollution[data$assignment==1])
## 
##  One Sample t-test
## 
## data:  data$post_pollution[data$assignment == 1]
## t = 9.8981, df = 39, p-value = 0.000000000003418
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  3033.840 4592.247
## sample estimates:
## mean of x 
##  3813.044

\[\mathbb{E}(Y|T=0)\]

options(scipen = 999999)

weighted.mean(data$post_pollution[data$assignment==0])
## [1] 5523.499
t.test(data$post_pollution[data$assignment==0])
## 
##  One Sample t-test
## 
## data:  data$post_pollution[data$assignment == 0]
## t = 8.2092, df = 19, p-value = 0.000000114
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  4115.223 6931.776
## sample estimates:
## mean of x 
##  5523.499

Los intervalos de confianza (IC) al 95% indican el rango de valores posibles para el promedio poblacional de polución si repitiéramos el muestreo muchas veces bajo las mismas condiciones.

Los IC no significan que haya 95% de probabilidad de que el promedio esté en ese rango, sino que si construyéramos intervalos del 95% en muchas muestras hipotéticas, aproximadamente el 95% de esos intervalos contendrían el verdadero valor del promedio poblacional.

4.Usando un modelo de regresión lineal para estimar el impacto de los nuevos protocolos sobre niveles de polución.

summary(lm(post_pollution~assignment,data))
## 
## Call:
## lm(formula = post_pollution ~ assignment, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5465.8 -2287.5  -313.2  2065.3  4516.8 
## 
## Coefficients:
##             Estimate Std. Error t value          Pr(>|t|)    
## (Intercept)   5523.5      589.8   9.365 0.000000000000333 ***
## assignment   -1710.5      722.4  -2.368            0.0212 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2638 on 58 degrees of freedom
## Multiple R-squared:  0.08815,    Adjusted R-squared:  0.07242 
## F-statistic: 5.607 on 1 and 58 DF,  p-value: 0.02125

Según el \(\hat{\beta}\) estimado el nivel promedio de contaminación del grupo tratado (T=1) es menor que el del grupo control (T=0). Esta regresión aun es solamente descriptiva y no puede interpretarse aún como un efecto causal del tratamiento, porque no controla posibles diferencias iniciales entre los grupos ni por variables confusoras. Esto quiere decir, que la estimación puede estar sesgada por factores omitidos.

El coeficiente asociado a la variable de tratamiento es significativo y está representando lo siguiente:

\[ \hat{\beta}=\mathbb{E}[Y|T=1] - \mathbb{E}[Y|T=0]\]

5. Estimando el impacto de los nuevos protocolos sobre los niveles de polución controlando por edad de la planta y número de empleados.

A continuación, desarrollaremos el siguiente modelo controlando por algunas variables disponibles de interés.Donde \(T\) es la variable tratamiento y \(X\) las variables control.

\[ y=\beta_{1}X+\beta_{2}T+\epsilon\]

summary(lm(post_pollution~assignment+age+employees,data))
## 
## Call:
## lm(formula = post_pollution ~ assignment + age + employees, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4937.8 -2504.5   -59.9  1984.0  4604.0 
## 
## Coefficients:
##               Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)  4827.9194  1108.0125   4.357 0.0000568 ***
## assignment  -1894.0121   739.8779  -2.560    0.0132 *  
## age           141.8294    98.0831   1.446    0.1537    
## employees      -0.1964     4.9597  -0.040    0.9686    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2636 on 56 degrees of freedom
## Multiple R-squared:  0.121,  Adjusted R-squared:  0.07388 
## F-statistic: 2.569 on 3 and 56 DF,  p-value: 0.06339

Ahora controlando por las variables disponibles (edad de la planta y número de empleados) podemos decir que estamos reduciendo el posible sesgo por variables omitidas que parecían sospechosas en la parte 2, como age. Por otro lado, el tratamiento sigue siendo significante cuando lo asociamos a menores niveles de polución.

Si bien los resultados se parecen al apartado anterior por lo que en cuanto a magnitud no serían distintas, pero si conceptualmente, ya que estamos controlando por las variables que parecían afectar a los resultados, como la edad de la planta y cantidad de empleados. Sin embargo, aun no podemos decir que hay efectos causales ya que podrían existir variables confusoras no observables que están afectando al resultado.

6. Preocupación que pueda tener sobre la validez interna y externa de este experimento aleatorio.

No se observan amenazas considerables sobre la validez interna del experimento ya que al comparar la asignación del tratamiento, esta no tiene diferencias significativas entre las variables que parecían ser significativamente diferentes en la sección 2, pero estas hacían referencia a los protocolos y no al balance entre tratamiento y controles.

t.test(data$pre_pollution ~ data$assignment)
## 
##  Welch Two Sample t-test
## 
## data:  data$pre_pollution by data$assignment
## t = -1.3773, df = 38.785, p-value = 0.1763
## alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
## 95 percent confidence interval:
##  -2968.431   563.771
## sample estimates:
## mean in group 0 mean in group 1 
##        3753.975        4956.305
t.test(data$age ~ data$assignment)
## 
##  Welch Two Sample t-test
## 
## data:  data$age by data$assignment
## t = -1.3773, df = 38.785, p-value = 0.1763
## alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
## 95 percent confidence interval:
##  -3.265274  0.620148
## sample estimates:
## mean in group 0 mean in group 1 
##        5.129372        6.451936

Lo anterior indica que los grupos de control (T=0) y tratamiento (T=1) están balanceados antes de la intervención. Aun así, para seguir mitigando el riesgo de que nuestro experimento no esté balanceado, además de usar las variables control como en la sección 5, podríamos recurrir a otros modelos como Diferencias en Diferencias (DiD) que permite ajustar por diferencias iniciales y capturar el cambio en contaminación antes y después del tratamiento.

data_longer=pivot_longer(data,cols=c(pre_pollution,post_pollution),
                         names_to = "period",
                         values_to = "pollution")

data_longer=mutate(data_longer, period = ifelse(period == "pre_pollution", 0, 1))


did_model <- lm(pollution ~ assignment * period + age + employees, data = data_longer)
summary(did_model)
## 
## Call:
## lm(formula = pollution ~ assignment * period + age + employees, 
##     data = data_longer)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6226.4 -1377.7   -61.3  1374.0  5631.8 
## 
## Coefficients:
##                     Estimate Std. Error t value           Pr(>|t|)    
## (Intercept)        1074.6521   773.5855   1.389            0.16748    
## assignment          509.3868   637.1414   0.799            0.42567    
## period             1769.5242   726.5558   2.435            0.01642 *  
## age                 525.4602    60.4587   8.691 0.0000000000000303 ***
## employees            -0.0982     3.0572  -0.032            0.97443    
## assignment:period -2912.7857   889.8456  -3.273            0.00141 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2298 on 114 degrees of freedom
## Multiple R-squared:  0.4322, Adjusted R-squared:  0.4073 
## F-statistic: 17.36 on 5 and 114 DF,  p-value: 0.0000000000009545

Aquí el coeficiente de interacción assignment:period (negativo y significativo) es el estimador DiD, y representa el efecto del tratamiento sobre la contaminación en el tiempo discreto ajustando por diferencias iniciales entre los grupos.

Por último, la preocupación en cuanto a la validez externa, es que los los resultados deben interpretarse con precaución ya que contamos con una muestra pequeña y con dos protocolos particulares, lo que no permitiría generalizar los resultados a otros sectores o regiones.