Introduccion

Este documento es un pequeño tutorial para trabajar con Propensity Score Matching, o pareamiento por puntajes de propensión. Esta técnica permite dar solución al problema de sesgo de selección. Es decir, ayuda a tomar algún experimento natural y realiza comparaciones entre variables con características similares para igualar el problema al de un verdadero experimento.

En este documento haremos análisis con una base de datos creada con el propósito de averiguar si una campaña publicitaria tiene efecto en las ventas de un producto.

El primer paso es colocar a R en la carpeta adecuada, que contiene nuestro archivo. El ejemplo está en la base de datos en linea. Descarga la base de datos y colócala en tu computadora en una carpeta que puedas indentificar. En R, teclea el siguiente código:

setwd("C:Direccion/de/La/Carpeta/de/Trabajo")

Debes de cambiar lo que está entre comillas por la dirección en la que se ubica el archivo de la base de datos. Puedes observar un video que muestre cómo identificar esta dirección de la capeta de trabajo.

Cargar los datos

La forma más directa de cargar los datos ya que estamos en la carpeta correcta es usar una función como read.csv(), cuyo primer argumento es el nombre de la base de datos. En nuestro caso, este nombre es la función mtk_camp.csv, de “Marketing Campaign”.

df <- read.csv("mkt_camp.csv")

Es posible agregar otros argumentos a la función. Por ejemplo, si las primeras \(n\) filas de la base de datos, es posible usar read.csv(base_de_datos, skip = n). Para más información sobre el funcionamiento de esta función, sólo debes teclear help(read.csv). Ahora verificaremos la base de datos. La función dim() nos regresa la dimensión de la base de datos y la función head() nos muestra las primeras filas de la base de datos.

dim(df)
## [1] 1000    4
head(df)
##   X...Edad Ingreso resp_camp compra
## 1       27   21840         0      0
## 2       26   44638         0      0
## 3       50   33023         0      0
## 4       23   28877         0      0
## 5       33   28636         1      1
## 6       29   21561         1      1

La base de datos tiene 1000 observaciones de 4 variables. Las variables incluyen Edad e Ingresos, resp_camp es una variable dicotómica que indica si el individuo respondió a la campaña, en cuyo caso resp_camp[i] = 0 para cualquier \(i\in[1,n]\), donde \(n\) es el número de observaciones, en este caso 1000. La variable compro también es dicotómica, e indica si el individuo compró nuestro producto o no.

Note que el nombre de la edad viene como X...Edad. Esto es por que hay un espacio en el nombre de la variable. Para cambiarlo a únicamente Edad

names(df)[1] <- 'Edad'
head(df)
##   Edad Ingreso resp_camp compra
## 1   27   21840         0      0
## 2   26   44638         0      0
## 3   50   33023         0      0
## 4   23   28877         0      0
## 5   33   28636         1      1
## 6   29   21561         1      1

Tratamiento vs Control

En el corazón de esta técnica está la idea de algún tratamiento. Puede ser un tratamiento médico, o se puede tratar de los beneficiarios de un programa social. En este caso, el tratamiento es la respuesta a nuestra publicidad. Vamos a separar la base de datos en el grupo de tratamiento y en el grupo de control

Tratamiento <- df[df$resp_camp == 1,]
colMeans(Tratamiento)
##        Edad     Ingreso   resp_camp      compra 
##    33.99529 24242.64706     1.00000     0.68000

Si seleccionas df$resp_camp únicamente y das Control+Enter, verás la columna de respuestas a la campaña. De hecho, para saber cuántas personas respondieron a la campaña, sólo es necesario ejecutar sum(df$resp_camp) y para saber cuánto es esto en proporción del total de observaciones, puedes usar mean(df$resp_camp). La primera linea entonces genera el subconjunto de aquellas observaciones en las cuales la comparación df$resp_camp == 1 es un valor verdadero (TRUE). La función colMeans() regresa el valor promedio por columna. En el caso de la edad, por ejemplo, la edad promedio de aquellos que respondieron a nuestro anuncio es de 43.65 y tienen un ingreso promedio de $29,207.83. El 68% de quienes respondieron a la campaña compraron nuestro producto.

Ahora veamos a aquellos que no respondieron a la campaña:

Control <- df[df$resp_camp == 0,]
colMeans(Control)
##          Edad       Ingreso     resp_camp        compra 
##    34.0991304 24643.7200000     0.0000000     0.1095652

Como podemos observar, cerca de un 10.95% de los individuos que no respondieron a nuestra campaña compraron nuestro producto. Puede ser que hayan entrado por otros medios o que hayan respondido a una publicidad pasada y decidido regresar.

Podemos hacer una estimación del efecto de las variables que tenemos disponibles en las compras por medio de una regresión lineal. Considere el modelo \[ C = \beta_0+\beta_1R+\beta_2+\beta_3I+\varepsilon \qquad \varepsilon\sim\mathcal{N}(0,1) \] Donde \(C\) es el nivel proporcional de compras con respecto al total de observaciones, \(R\) es la proporción que responde a la campaña publicitaria, \(E\) es la edad y el nivel de ingreso está determinado por \(I\). Esta regresión se realiza en R con

modelo1 <- lm(compra ~ resp_camp + Edad + Ingreso, data = df)

summary(modelo1)
## 
## Call:
## lm(formula = compra ~ resp_camp + Edad + Ingreso, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.6804 -0.1097 -0.1095  0.3199  0.8907 
## 
## Coefficients:
##                  Estimate    Std. Error t value Pr(>|t|)    
## (Intercept)  0.1096137815  0.0542011974   2.022   0.0434 *  
## resp_camp    0.5704382873  0.0247111118  23.084   <2e-16 ***
## Edad        -0.0000095209  0.0012851921  -0.007   0.9941    
## Ingreso      0.0000000112  0.0000010808   0.010   0.9917    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3862 on 996 degrees of freedom
## Multiple R-squared:  0.3486, Adjusted R-squared:  0.3467 
## F-statistic: 177.7 on 3 and 996 DF,  p-value: < 2.2e-16

En este caso particular, dado que ambos grupos son muy similares y están disribuidos de manera alreatoria, no es posible discernir ningún patrón particular. De hecho, la edad y el ingreso no son significativas, lo cual es posible observar en el valor P que es muy alto, es decir, la probabilidad de que los resultados se hayan obtenido por casualidad es alta, debida a que existe mucho ruido en esta información. Este ruido existe porque las variables de Edad e Ingreso están distribuidas de manera uniforme. En cambio, el efecto de la respuesta a la campaña es positivo y estadísticamente significativo. Podemos guardar el efecto de la campaña en las compras en una variable que llamaremos efecto.

efecto <- modelo1$coefficients[2]
efecto
## resp_camp 
## 0.5704383

Puntajes de Propensión

Haremos una regresión logística para obtener puntajes de propensión, es decir, queremos clasificar a nuestros sujetos en base a que tan propensos son a entrar en nuestra campaña, dados su edad e ingresos. Considere el modelo \[ P(R = 1|E,I) = \frac{e^{\beta_0+\beta_1E+\beta_2I}}{1+e^{\beta_0+\beta_1E+\beta_2I}} \] Donde \(e\) es la base del logaritmo natural. Este modelo lo podemos ejecutar en R con el siguiente código:

modelo_pp <- glm(resp_camp ~ Edad + Ingreso, family = binomial("logit"), data = df)
summary(modelo_pp)
## 
## Call:
## glm(formula = resp_camp ~ Edad + Ingreso, family = binomial("logit"), 
##     data = df)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -1.085  -1.054  -1.034   1.303   1.342  
## 
## Coefficients:
##                 Estimate   Std. Error z value Pr(>|z|)
## (Intercept) -0.183246855  0.277562659  -0.660    0.509
## Edad        -0.001227284  0.006732608  -0.182    0.855
## Ingreso     -0.000003160  0.000005662  -0.558    0.577
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1363.7  on 999  degrees of freedom
## Residual deviance: 1363.4  on 997  degrees of freedom
## AIC: 1369.4
## 
## Number of Fisher Scoring iterations: 4

De nuevo, los efectos son pequeños y poco significativos. Esto es porque, como la edad y el ingreso están distribuidos de manera uniforme, no existe un patrón claro. Esto implica que no existe sesgo de selección en nuestro problema. Sin embargo, seguiremos, considerando que esto es un caso atípico y que en problemas posteriores nos encontraremos con casos en los que el sesgo de selección si está presente. De hecho, este tipo de resultado es el que deberíamos esperar idealmente si pudiéramos realizar un experimento en el que aplicamos el tratamiento a los sujetos seleccionados de manera aleatoria. Para observar esto, ejecute la siguiente línea:

df$Puntajes <- modelo_pp$fitted.values

hist(df$Puntajes[df$resp_camp == 1], main = "Puntajes con resp = 1")

hist(df$Puntajes[df$resp_camp == 0], main = "Puntajes con resp = 0")

Note que los puntajes están bien distribuidos y concentrados alrededor de 0.425, que es el puntaje medio.

Pareamiento de puntajes de propensión

Primero generamos un vector con los dos nombres de las variables con las que trabajaremos y cargaremos la biblioteca tableone para generar tablas estratificadas. Esta tabla la presentamos para mostrar, nuevamente, que ambos estratos (SI respondió a la campaña y NO respondió a la campaña) son muy similares. En el caso de que estos sean muy diferentes, lo podemos observar porque las diferencias medias estandarizadas en la columna SMD serían valores mucho más grandes.

xvars <- c("Edad","Ingreso")
library(tableone)

table1 <-CreateTableOne(vars = xvars, strata = "resp_camp", data = df, test = FALSE)
print(table1, smd = T)
##                      Stratified by resp_camp
##                       0                   1                   SMD   
##   n                        575                 425                  
##   Edad (mean (sd))       34.10 (9.62)        34.00 (9.38)      0.011
##   Ingreso (mean (sd)) 24643.72 (11555.95) 24242.65 (10979.83)  0.036

Nota que, si la biblioteca tableone no está descargada en R, la consola marcará un error. Para descargar la biblioteca, usa la función ìnstall.packages("tableone"). Al correr esta función R automáticamente descargará la biblioteca desde CRAN.

Para implementar el algoritmo de pareamiento se utiliza la biblioteca MatchIt. Asegúrate de tener descargada la biblioteca siguiendo los mismos pasos anteriormente descritos y corre el siguiente código:

library(MatchIt)
match1 <- matchit(modelo_pp, method = "nearest", data = df)
match1.data <- match.data(match1)

table_match1 <- CreateTableOne(vars = xvars,
                               strata = "resp_camp",
                               data = match1.data)
print(table_match1, smd = T)
##                      Stratified by resp_camp
##                       0                   1                   p      test
##   n                        425                 425                       
##   Edad (mean (sd))       34.12 (9.46)        34.00 (9.38)      0.847     
##   Ingreso (mean (sd)) 24329.01 (11394.09) 24242.65 (10979.83)  0.910     
##                      Stratified by resp_camp
##                       SMD   
##   n                         
##   Edad (mean (sd))     0.013
##   Ingreso (mean (sd))  0.008
plot(match1, type = "hist")

En la primera línea, se carga a R la biblioteca MatchIt. En la segunda se realiza un pareamiento por medio del modelo logit descrito anteriormente y con el algoritmo del vecino más cercano. La biblioteca MatchIt contiene la función match.data(), que presenta la información del pareamiento. Esta información se puede resumir con print(), y la función plot() muestra los histogramas respectivos del tratamiento y el control previos al pareamiento y con pareamiento. Nuevamente, dado que nuestros datos ya están distribuídos aleatoriamente, ambos histogramas se verán idénticos.

Finalmente, para realizar un comparativo podemos extraer el vector de los datos de compra en aquellos individuos que respondieron a la publicidad y los que no. La diferencia media será nuestro nivel de respuesta de la publicidad en las compras.

y_trat <- match1.data$compra[match1.data$resp_camp == 1]
y_ctr <- match1.data$compra[match1.data$resp == 0]

diferencia <- y_trat - y_ctr
mean(diferencia)
## [1] 0.5788235

Note que, en este caso, este nivel coincide con el nivel extraído en la regresión y almacenado en la variable efecto.

Ejercicio

Considere esta base de datos diferente. Realice el mismo ejercicio de Pareamiento. ¿Existe sesgo de selección? ¿Cuál es el efecto de la publicidad en este caso? ¿Conviene realizar la campaña?