Propensity Score

Explicación y ejemplo de replicación

1. Introducción

1.2 El problema del “Criterio de puerta trasera”

En modelos gráficos acíclicos existe la idea del “criterio de la puerta trasera”, que consiste en la mediación de una variable exógena entre una variable dependiente, por ejemplo la mortalidad (Y) y una independiente, como por ejemplo si la persona fuma (X) como se muestra en la siguiente figura:.

En este esquema DAG vemos que la edad (A) genera una “puerta trasera” en términos de que A afecta tanto directamente a la mortalidad pero también afecta la mortalidad dado si fuma. Por lo tanto, si tenemos una estrategia para solventar el problema de la puerta trasera, entonces tenemos capacidad para usar una estrategia para identificar un efecto causal

La subclasificación es un método que satisface el criterio de puerta trasera ponderando las diferencias de las medias a partir de “stratas” de valores específicos.La subclasificación ajustará las diferencias en medias de manera equivalente a su strata contrafactual, logrando un balance entre ambas condiciones.

Para ello es importante lograr el supuesto de independiencia condicional. Esto implica que:

\[(Y^1, Y^0) \perp D | X\]

\[E [ Y^1 | D = 1,X] = E [ Y^1 | D = 0,X] \]

\[ E[Y^0 | D = 1, X] = E[Y^0 | D = 0,X]\]

Para abordar estos temas y “cerrar” los problemas de puerta trasera existen diversos métodos, sin embargo nos centraremos en el propensity score.

1.3 Propensity score

Creado por Donald Rubin a mediados de los 70’s, el propensity score es similar al matching con vecino más cercano.

El propensity score no ha llegado a muchos campos, a diferencia de las regresiones discontinúas, ya que se sospecha de que realmente se logre el CIA en data sets de este tipo. Sin embargo no hay razones teóricas para desechar el procedimiento diseñado para estimar efectos causales. Sólo conocimientos previos y familiaridad profunda con los detalles institucionales puede decirte cuál es la estrategia apropiada, y en la medida en que el criterio de puerta trasera pueda ser cumplido, entonces los métodos de matching pueden ser apropiados.

Como señala el texto de Cunningham: “El propensity score toma estas covarianzas necesarias, estima un modelo de máxima verosimilitud de la probabilidad condicional de tratamiento (usualmente un logit o probit de manera de asegurarse que los valores ajustados están rondando el 0 y el 1), y usa los valores predichos de la estimacion para colapsar dichas covariables en una escala única llamada propensity score”

Pero hay algunos matices: “La idea del propensity score es comparar unidades que, basados en observables, tienen probabilidades similares de ser ubicado dentro del grupo tratamiento aunque esas unidades difieran en cuanto a la asignación real del tratamiento”

“Cuando dos unidades tienen la misma posibilidad de ser ubicados en tratamiento o control, entonces”En la medida que las unidades A y B tienen el mismo propensity score de 0.6, pero uno es tratamiento y otro no, y la CIA se mantiene en los datos, entonces las diferencias entre los resultados observados es atribuíble al tratamiento”

Existe otro supuesto importante para esto que es el common support assumption; esto simplemente implica que hayan unidades en tratamiento y control a través del propensity score estimado. En este caso el supuesto se cumple porque tenemos tanto en control como tratamiento unidades con 0.6 de probabilidad.

2. Replicación investigación de NSW Job Trainning program

La National Supported Work Demostration(NSW) Job-training program es un ṕrograma diseñado para ayudar a trabajadores desaventajados que tienen pocas habilidades laborales para moverse en el mercado laboral, dándoles experiencia y apoyo.

Se tienen dos grupos al azar: el grupo tratamiento que tendrá todos los beneficios del NSWJTP y el tratamiento control son personas que siguen buscando trabajo por cuenta propia. Quienes participan del grupo de tratamiento:

• Tienen trabajo garantizado entre 9 a 18 meses, según su situación
• Se dividen en grupos de tres o cinco participantes para trabajar juntos
• Se reunen con un consejero de la NSW para discutir problemas con el programa y su performance.
• Se les paga por su trabajo, aunque menos que en el mercado laboral convencional; sin embargo puede aumentar según méritos.

El NSWJTP al ser un experimento aleatoriozado, basta con el cálculo del ATE para medir sus efectos:

\[\frac{1}{N_T}\sum_{D_i =1} Y_i - \frac{1}{N_C} \sum_{D_i =0} Y_i \approx E[Y^1 - y^0] \approx ATE \]

Según lo indicado en el texto de Cunningham, el grupo tratado tuvo un aumento entre $900 a $1800 dólares en salario por sobre el grupo de control. Pero ahora se quiere ahondar en otros datos observacionales que tienen otro orden, con datos de corte transversal.

Lalonde (1986) realiza una investigación usando distintos métodos de inferencia causal que se ocupaban en la época y todos salieron muy diferentes a lo que apareció de manera experimental, incluso en signos negativos. ¿Por qué ocurre eso? por el sesgo de selección:

\[E[Y^0 | D = 1] \neq E [Y^0 | D = 0]\]

Esto quiere decir que las personas que participaron del NSWJTP tenían salarios menores que la población general de la muestra de corte transversal. Pero además de esto una violación del principio de independencia sería que las covariables de los grupos de control y tratamiento están desbalanceadas una de otra. Como se ve en la tabla 5.12 del libro:

2.1 Análisis de replicación

Average Treatment Effect

Dehejia y Wahba (1999) Querían re-analizar los datos usando propensity score para observar si había una mejora en la estimación, así como mostrar el valor de la técnica como un mecanismo de diagnóstico. Primero chequea con una estimación de maxima verosimilitud el principio de common support, haciendo un histograma de los propensity score entre tratamiento y control. Al constatarlo se ve que el solapamiento era tan malo que tuvieron que eliminar 12611 observaciones, dado que el grupo de control tenía scoring que el grupo de tratamiento no tenía. Esto resume la idea de que el sesgo de selección de las observables donde se concentran en tratamiento o control son relevantes a la hora de hacer inferencia causal.

El resultado final fue que, dependiendo de la cantidad de covariables exógenas integradas a la regresión, los trabajadores aumentaban desde 1672 a 1794US en quienes participaron del tratamiento,

Obtengo un resumen sumario de quienes tuvieron tratamiento:

nsw_dw %>% 
  filter(treat == 1) %>% 
  summary(re78)

##    data_id              treat        age             educ      
##  Length:185         Min.   :1   Min.   :17.00   Min.   : 4.00  
##  Class :character   1st Qu.:1   1st Qu.:20.00   1st Qu.: 9.00  
##  Mode  :character   Median :1   Median :25.00   Median :11.00  
##                     Mean   :1   Mean   :25.82   Mean   :10.35  
##                     3rd Qu.:1   3rd Qu.:29.00   3rd Qu.:12.00  
##                     Max.   :1   Max.   :48.00   Max.   :16.00  
##      black             hisp              marr           nodegree     
##  Min.   :0.0000   Min.   :0.00000   Min.   :0.0000   Min.   :0.0000  
##  1st Qu.:1.0000   1st Qu.:0.00000   1st Qu.:0.0000   1st Qu.:0.0000  
##  Median :1.0000   Median :0.00000   Median :0.0000   Median :1.0000  
##  Mean   :0.8432   Mean   :0.05946   Mean   :0.1892   Mean   :0.7081  
##  3rd Qu.:1.0000   3rd Qu.:0.00000   3rd Qu.:0.0000   3rd Qu.:1.0000  
##  Max.   :1.0000   Max.   :1.00000   Max.   :1.0000   Max.   :1.0000  
##       re74            re75            re78        
##  Min.   :    0   Min.   :    0   Min.   :    0.0  
##  1st Qu.:    0   1st Qu.:    0   1st Qu.:  485.2  
##  Median :    0   Median :    0   Median : 4232.3  
##  Mean   : 2096   Mean   : 1532   Mean   : 6349.1  
##  3rd Qu.: 1291   3rd Qu.: 1817   3rd Qu.: 9643.0  
##  Max.   :35040   Max.   :25142   Max.   :60307.9

Computo la media de Salario del grupo tratamiento y la alojo en una nueva variable llamada nsw_dw$y1:

mean1 <- nsw_dw %>% 
  filter(treat == 1) %>% 
  pull(re78) %>% 
  mean()

nsw_dw$y1 <- mean1

Realizo un análisis sumario del grupo de control:

nsw_dw %>% 
  filter(treat == 0) %>% 
  summary(re78)

##    data_id              treat        age             educ      
##  Length:260         Min.   :0   Min.   :17.00   Min.   : 3.00  
##  Class :character   1st Qu.:0   1st Qu.:19.00   1st Qu.: 9.00  
##  Mode  :character   Median :0   Median :24.00   Median :10.00  
##                     Mean   :0   Mean   :25.05   Mean   :10.09  
##                     3rd Qu.:0   3rd Qu.:28.00   3rd Qu.:11.00  
##                     Max.   :0   Max.   :55.00   Max.   :14.00  
##      black             hisp             marr           nodegree     
##  Min.   :0.0000   Min.   :0.0000   Min.   :0.0000   Min.   :0.0000  
##  1st Qu.:1.0000   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:1.0000  
##  Median :1.0000   Median :0.0000   Median :0.0000   Median :1.0000  
##  Mean   :0.8269   Mean   :0.1077   Mean   :0.1538   Mean   :0.8346  
##  3rd Qu.:1.0000   3rd Qu.:0.0000   3rd Qu.:0.0000   3rd Qu.:1.0000  
##  Max.   :1.0000   Max.   :1.0000   Max.   :1.0000   Max.   :1.0000  
##       re74              re75              re78             y1      
##  Min.   :    0.0   Min.   :    0.0   Min.   :    0   Min.   :6349  
##  1st Qu.:    0.0   1st Qu.:    0.0   1st Qu.:    0   1st Qu.:6349  
##  Median :    0.0   Median :    0.0   Median : 3139   Median :6349  
##  Mean   : 2107.0   Mean   : 1266.9   Mean   : 4555   Mean   :6349  
##  3rd Qu.:  139.4   3rd Qu.:  650.1   3rd Qu.: 7288   3rd Qu.:6349  
##  Max.   :39570.7   Max.   :23032.0   Max.   :39484   Max.   :6349

calculo la media del grupo de control y la guardo en una nueva variable llamada nsw_sw$y0:

mean0 <- nsw_dw %>% 
  filter(treat == 0) %>% 
  pull(re78) %>% 
  mean()

nsw_dw$y0 <- mean0

Y finalmente calculo la diferencia entre las medias:

ate <- unique(nsw_dw$y1 - nsw_dw$y0)

nsw_dw <- nsw_dw %>% 
  filter(treat == 1) %>% 
  select(-y1, -y0)

2.2 Otros casos:

En este segundo caso se utilizó como grupo de control una muuestra aleatoria del Current Population Survey. Esto implica que hay un altísimo sesgo de selección que será controlado con los propensity scores

Se realiza el modelo logit para estimar y obtener los valores ajustados:

# estimating
logit_nsw <- glm(treat ~ age + agesq + agecube + educ + educsq + 
                   marr + nodegree + black + hisp + re74 + re75 + u74 +
                   u75 + interaction1, family = binomial(link = "logit"), 
                 data = nsw_dw_cpscontrol)

nsw_dw_cpscontrol <- nsw_dw_cpscontrol %>% 
  mutate(pscore = logit_nsw$fitted.values)

Se calculan las medias del propensity score en el grupo de control y en el tratamiento:

# mean pscore 
pscore_control <- nsw_dw_cpscontrol %>% 
  filter(treat == 0) %>% 
  pull(pscore) %>% 
  mean()

pscore_contro_q <-nsw_dw_cpscontrol %>% 
  filter(treat == 0) %>% 
 pull(pscore) %>% 
  mean()


pscore_treated <- nsw_dw_cpscontrol %>% 
  filter(treat == 1) %>% 
  pull(pscore) %>% 
  mean()

Se realizan histogramas para observar el propensity Score por cada condición:

# histogram
nsw_dw_cpscontrol %>% 
  filter(treat == 0) %>% 
  ggplot() +
  geom_histogram(aes(x = pscore))

nsw_dw_cpscontrol %>% 
  filter(treat == 1) %>% 
  ggplot() +
  geom_histogram(aes(x = pscore))

Si observamos los datos de la tabla, el propensity score del grupo tratamiento es de 0.43 y el del grupo control (CPS) de 0.007, lo que muestra que las probabilidades están mal distribuídas. esto se comprueba en los hiistogramas donde el grupo tratamiento aparece bastante bien distribído en cada percentil, sin embargo en el grupo control están muy cargadas a los valores 0, lo que se puede interpresar como que las características de los individuos en el grupo tratamiento son muy raras en la muestra de CPS. Por lo tanto si son tan sdistinto,s se verá reflejado en un propensity score muy distinto.

En resumen, esto indicaría que el principio de apoyo común no esté dado por un bajo solapamiento de concidencias.

Volviendo a la teoría,para obtener el propensity score es importante entender que se usan los coeficientes de la regresión logit para estimar la probabilidad condicional del tratamiento, asumiendo que las probabilidades están basadas en una distribución logística acumulada:

\[P(D = 1| X) = F(\beta_0 + \gamma Treat + \alpha X)\] donde \(F()= \frac{e}{1+e}\) y \(X\) es una covariable exógena que incluímos en el modelo.

Es importante también señalar que el propensity score es una probabilidad del tratamiento independiente del tratamiento actual asignado, por lo que la definición del propensity scotre es “la probabilidad de seleción condicional en las variables de confusión”:

\[p(X) = Pr(D = 1 | X)\]

Todo esto quiere decir que, para lograr independencia, asumiendo CIA, lo único que debemos hacer es condicionar al propensity score. Sumado a esto, el teoría del propensity score señala que sólo debes controlar por covariables que determinan la verosimilitud de que una unidad reciba el tratamiento.

Esto redunda que podemos estimar el ATE ponderando apropiadamente:

\[Pr(D = 1 | X, p(X)) = Pr(D = 1 | X)\] \[ = p(X)\] Esto significa que X ya no condiciona a D, sino que ahora es la \(p(X)\) quien hace mediar esa relación, por tanto D y X son independientes una de otra.

2.3 Weighting propensity score

Existen diversas formas de cómo los ponderadores de los propensity score son integrados a diseños, y la probabilidad inversa logra dar buenos resultados.

Una forma es considerando que cada propensity score individual es un ponderador del reusltado inidividual, que cuando se agrega tiene el potencial de calcular el ATE:

\[ \delta_{ATE} = E [Y^1 - Y^0]\] \[ = E [Y ·\frac{D - p(X)}{p(X)(1-p(X))}]\] ¸\[ = \frac{1}{Pd(D = 1)}· E[Y ·\frac{D - p(X)}{1 - p(X)}] \]

Para obtener el ATE y el ATT se realizan dos procedimientos: primero se estima el propensity score usando estimaciones logit o probit. En el segundo paso el investigador usa el scoring estimado para producir muestras de uno de los estimadores del ATE que pueden ser escritos como:

\[\hat{\delta}_{ATE} = \frac{1}{N}\sum_{i=1}^{N}Y_i · \frac{D_i - \hat{p}(X_i)}{\hat{p}(X_i)· (1 - \hat{p}(X_i))}\]

O para el ATT:

\[\hat{\delta}_{ATT} = \frac{1}{N_T} \sum_{i=1}^NY_i · \frac{D_i - \hat{p}(X_i)}{1 - \hat{p}(X_i)} \]

Para calcular la varianza de esta estimación se puede usar la técnica de bootstrapping, el cual nos permitirá tener la distribución de los parámetros y calcular la desviación estándar. Sin embargo podría tener algunos problemas con valores extremos, por lo que usar las probabilidades negativas también es una buena opción. En ella pondera dentro de cada grupo, el cual tiene a hacerlo más estable:

\[\hat{\delta}_{ATT} = \left [ \sum_{i=1}^1 \frac{Y_iD_i}{\hat{p}}] \right] / \left[\sum_{i=1}^{N} \frac{D_i}{\hat{p}} \right] - \left[\sum_{i=1}^N \frac{Y_i(1 - D_i)}{(1 - \hat{p})} \right] / \left[ \sum_{i=1}^N \frac{(1 - D_i)}{(1 -\hat{p})} \right]\]

Para realizar la estimación creamos las variables D1 Y d0, y luego ver las diferencias entre medias usando la variable de ingreso re78 considerando el propensity score:

#continuation
N <- nrow(nsw_dw_cpscontrol)
#- Manual with non-normalized weights using all data
nsw_dw_cpscontrol <- nsw_dw_cpscontrol %>% 
  mutate(d1 = treat/pscore,
         d0 = (1-treat)/(1-pscore))

s1 <- sum(nsw_dw_cpscontrol$d1)
s0 <- sum(nsw_dw_cpscontrol$d0)

nsw_dw_cpscontrol <- nsw_dw_cpscontrol %>% 
  mutate(y1 = treat * re78/pscore,
         y0 = (1-treat) * re78/(1-pscore),
         ht = y1 - y0)

mean(nsw_dw_cpscontrol$ht)

## [1] -11876.79

En el siguiente código se comparen los resultados de las estimaciones ponderadasnormalizadas vs las que no lo están. vemos que en la ponderación que no está normalizada el ingreso es de $2006.365 US mientras que la normalizada muestra una diferencia de $1806.73 US

#- Manual with normalized weights
nsw_dw_cpscontrol <- nsw_dw_cpscontrol %>% 
  mutate(y1 = (treat*re78/pscore)/(s1/N),
         y0 = ((1-treat)*re78/(1-pscore))/(s0/N),
         norm = y1 - y0)

nsw_dw_cpscontrol %>% 
  pull(ht) %>% 
  mean()

## [1] -11876.79

nsw_dw_cpscontrol %>% 
  pull(norm) %>% 
  mean()

## [1] -7238.14

Luego se realiza un corte de datos evitando los valores extremos. se seleccionan los elementos con valores mayores o iguales a 0.1 y menores o iguales a 0.9:

#-- trimming propensity score
nsw_dw_cpscontrol <- nsw_dw_cpscontrol %>% 
  select(-d1, -d0, -y1, -y0, -ht, -norm) %>% 
  filter(!(pscore >= 0.9)) %>% 
  filter(!(pscore <= 0.1))

N <- nrow(nsw_dw_cpscontrol)

Posteriormente, se vuelven a computar los cálculos con los valores cortados:

#- Manual with non-normalized weights using trimmed data
nsw_dw_cpscontrol <- nsw_dw_cpscontrol %>% 
  mutate(d1 = treat/pscore,
         d0 = (1-treat)/(1-pscore))

s1 <- sum(nsw_dw_cpscontrol$d1)
s0 <- sum(nsw_dw_cpscontrol$d0)

nsw_dw_cpscontrol <- nsw_dw_cpscontrol %>% 
  mutate(y1 = treat * re78/pscore,
         y0 = (1-treat) * re78/(1-pscore),
         ht = y1 - y0)

#- Manual with normalized weights with trimmed data
nsw_dw_cpscontrol <- nsw_dw_cpscontrol %>% 
  mutate(y1 = (treat*re78/pscore)/(s1/N),
         y0 = ((1-treat)*re78/(1-pscore))/(s0/N),
         norm = y1 - y0)

nsw_dw_cpscontrol %>% 
  pull(ht) %>% 
  mean()

## [1] 2006.365

nsw_dw_cpscontrol %>% 
  pull(norm) %>% 
  mean()

## [1] 1806.73

Con los valores cortadoos vemos que se ajustan mucho mejor a los datos experimentales originales, lo que conmfirmarìa esta estrategia de análisis causal.