library(tidyverse) # ggplot(), %>%, mutate(), and friends
library(ggdag) # Make DAGs
library(dagitty) # Do DAG logic with R
library(broom) # Convert models to data frames
library(MatchIt) # Match things
library(modelsummary) # Make side-by-side regression tables
set.seed(1234) # Make all random draws reproduciblePareo y Probabilidad Inversa Ponderada
Para este ejercicio utilizaremos la base de datos mosquito_nets.csv disponible aquí:
Antecedentes del programa
Los investigadores están interesados por saber si el uso de mosquiteros reduce el riesgo de contraer malaria. Recabaron datos de 1752 hogares en un país no identificado e incluyeron variables relacionadas con factores ambientales, salud individual y las características del hogar. Los datos no son experimentales : los investigadores no controlan quién usa mosquiteros, y cada hogar decide si solicita mosquiteros gratuitos o los compra, así como si los usa si los tiene.
El archivo CSV contiene las siguientes columnas:
Riesgo de malaria (
malaria_risk): Probabilidad de que alguien en el hogar se infecte con malaria. Se mide en una escala de 0 a 100, donde los valores más altos indican un mayor riesgo.Mosquitero (
netynet_num): Variable binaria que indica si en el hogar se utilizaban mosquiteros.Elegible para el programa (
eligible): Una variable binaria que indica si el hogar es elegible para el programa de red gratuita.Ingresos (
income): Ingresos mensuales del hogar, en dólares estadounidenses.Temperaturas nocturnas (
temperature): La temperatura media nocturna, en grados Celsius.Salud (
health): Autopercepción de la salud en el hogar. Medida en una escala de 0 a 100, donde los valores más altos indican mejor salud.Número de personas en el hogar (
household): Número de personas que viven en el hogar.Resistencia a los insecticidas (
resistance): Algunas cepas de mosquitos son más resistentes a los insecticidas y, por lo tanto, presentan un mayor riesgo de infectar a las personas con malaria. Esto se mide en una escala de 0 a 100, donde los valores más altos indican una mayor resistencia.
Nuestro objetivo
Nuestro objetivo en este ejemplo es estimar el efecto causal del uso de mosquiteros sobre el riesgo de malaria utilizando únicamente datos observacionales. No se trata de un ensayo clínico aleatorio (ECA), por lo que podría parecer un poco impreciso afirmar causalidad. Sin embargo, si podemos dibujar un DAG correcto y ajustar los nodos correctos, podemos aislar la relación entre el uso de mosquiteros y la malaria y hablar de causalidad.
Dado que se trata de datos simulados, conocemos el verdadero efecto causal del programa de red, ya que lo incorporé a los datos. El verdadero efecto promedio del tratamiento (EPT) es -10. El uso de un mosquitero reduce el riesgo de malaria en 10 puntos, en promedio.
¡Veamos si podemos encontrar ese efecto de 10 puntos!
Carguemos datos y paquetes
# Load the data.
# It'd be a good idea to click on the "nets" object in the Environment panel in
# RStudio to see what the data looks like after you load it
nets <- read_csv("data/mosquito_nets.csv")Rows: 1752 Columns: 10
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
dbl (8): id, net_num, malaria_risk, income, health, household, temperature, ...
lgl (2): net, eligible
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.DAG y conjuntos de ajuste
Antes de ejecutar cualquier modelo, necesitamos determinar qué debemos ajustar. Para ello contruímos el siguiente modelo causal para mostrar el proceso de generación de datos que subyace a la relación entre el uso de mosquiteros y el riesgo de malaria:
mosquito_dag <- dagify(
malaria_risk ~ net + income + health + temperature + resistance,
net ~ income + health + temperature + eligible + household,
eligible ~ income + household,
health ~ income,
exposure = "net",
outcome = "malaria_risk",
coords = list(x = c(malaria_risk = 7, net = 3, income = 4, health = 5,
temperature = 6, resistance = 8.5, eligible = 2, household = 1),
y = c(malaria_risk = 2, net = 2, income = 3, health = 1,
temperature = 3, resistance = 2, eligible = 3, household = 2)),
labels = c(malaria_risk = "Risk of malaria", net = "Mosquito net", income = "Income",
health = "Health", temperature = "Nighttime temperatures",
resistance = "Insecticide resistance",
eligible = "Eligible for program", household = "Number in household")
)
ggdag_status(mosquito_dag, use_labels = "label", text = FALSE) +
guides(fill = "none", color = "none") + # Disable the legend
theme_dag()
Siguiendo la lógica del cálculo do , podemos encontrar todos los nodos que confunden la relación entre el uso de la red y el riesgo de malaria, ya que estos nodos abren caminos de puerta trasera y distorsionan el efecto causal que nos interesa. Podemos hacerlo gráficamente, buscando cualquier nodo que apunte tanto al uso de la red como al riesgo de malaria, o podemos usar R:
adjustmentSets(mosquito_dag){ health, income, temperature }Con base en las relaciones entre todos los nodos del DAG, el ajuste por salud, ingresos y temperatura es suficiente para cerrar todas las puertas traseras e identificar la relación entre el uso de la red y el riesgo de malaria.
Estimación ingenua de correlación no causal
nets %>%
group_by(net) %>%
summarize(number = n(),
avg = mean(malaria_risk))# A tibble: 2 × 3
net number avg
<lgl> <int> <dbl>
1 FALSE 1071 41.9
2 TRUE 681 25.6O podemos usar la siguiente regresión:
model_wrong <- lm(malaria_risk ~ net, data = nets)
tidy(model_wrong)# A tibble: 2 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 41.9 0.405 104. 0
2 netTRUE -16.3 0.649 -25.1 2.25e-119Según esta estimación, el uso de mosquitero se asocia con una reducción de 16,33 puntos porcentuales en el riesgo de malaria, en promedio. Sin embargo, legalmente no podemos hablar de esto como un efecto causal, ya que existen variables de confusión que deben considerarse.
Pareo
Podemos utilizar técnicas de correspondencia para emparejar observaciones similares y hacer el supuesto de no confusión: si vemos dos observaciones que son prácticamente idénticas y una utilizó una red y la otra no, la elección de utilizar una red fue aleatoria.
Dado que sabemos por el DAG que los ingresos, las temperaturas nocturnas y la salud ayudan a causar tanto el uso de mosquiteros como el riesgo de malaria (¡y confunden esa relación!), intentaremos encontrar observaciones con valores similares de ingresos, temperaturas y salud que usaron y no usaron mosquiteros.
Podemos usar la matchit()función del paquete {MatchIt} de R para emparejar puntos según la distancia de Mahalanobis. Hay muchas otras opciones disponibles; consulta la [documentación del paquete]((http://gking.harvard.edu/matchit) para obtener más información.
Podemos incluir la replace = TRUEopción para que los puntos que ya han sido emparejados puedan volver a ser emparejados (es decir, no estamos forzando una coincidencia uno a uno; en su lugar, tenemos una coincidencia uno a muchos).
Paso 1: Pre-procesamiento
matched_data <- matchit(net ~ income + temperature + health,
data = nets,
method = "nearest",
distance = "mahalanobis",
replace = TRUE)
summary(matched_data)
Call:
matchit(formula = net ~ income + temperature + health, data = nets,
method = "nearest", distance = "mahalanobis", replace = TRUE)
Summary of Balance for All Data:
Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
income 955.1938 872.7526 0.4089 1.3633 0.1044
temperature 23.3809 24.0880 -0.1685 1.0846 0.0424
health 54.9090 48.0570 0.3619 1.2083 0.0714
eCDF Max
income 0.1983
temperature 0.0972
health 0.1683
Summary of Balance for Matched Data:
Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
income 955.1938 950.1836 0.0248 1.1061 0.0093
temperature 23.3809 23.4066 -0.0061 1.0335 0.0067
health 54.9090 54.5433 0.0193 1.0714 0.0073
eCDF Max Std. Pair Dist.
income 0.0485 0.1202
temperature 0.0264 0.1196
health 0.0264 0.1239
Sample Sizes:
Control Treated
All 1071. 681
Matched (ESS) 320.94 681
Matched 436. 681
Unmatched 635. 0
Discarded 0. 0Aquí podemos ver que los 681 usuarios de la red se emparejaron con no usuarios de apariencia similar (436 de ellos). 635 personas no se emparejaron y serán descartadas. Si tienes curiosidad, puede ver qué filas tratadas se emparejaron con qué filas de control ejecutando matched_data$match.matrix.
Podemos crear un nuevo marco de datos de esas coincidencias con match.data():
matched_data_for_real <- match.data(matched_data)Paso 2: Estimacion
Ahora que los datos se han emparejado, deberían funcionar mejor para el modelado. Además, dado que utilizamos ingresos, temperaturas y salud en el proceso de emparejamiento, hemos ajustado esos nodos DAG y hemos cerrado esas puertas traseras, por lo que nuestro modelo puede ser bastante simple:
model_matched <- lm(malaria_risk ~ net,
data = matched_data_for_real)
tidy(model_matched)# A tibble: 2 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 38.3 0.599 64.0 0
2 netTRUE -12.7 0.768 -16.6 2.28e-55La disminución de 12,73 puntos que se observa aquí es mejor que la estimación ingenua, pero no refleja el verdadero efecto causal de 10 puntos (que incorporé en los datos). Quizás se deba a que las coincidencias no son excelentes, o quizás a que descartamos demasiados datos. Hay una serie de diagnósticos que puedes consultar para comprobar la precisión de la coincidencia (consulta la documentación de MatchIt para ver ejemplos).
En realidad, la causa más probable de la estimación incorrecta es un desequilibrio en los datos. Dado que establecimos replace = TRUE, no realizamos un emparejamiento 1:1; las observaciones sin tratamiento se emparejaron con más de una observación tratada. Como resultado, las observaciones con emparejamiento múltiple se contabilizan en exceso y tienen demasiada importancia en el modelo. Afortunadamente, matchit()nos proporciona una columna llamada weightsque nos permite reducir la escala de las observaciones con sobreemparejamiento al ejecutar el modelo. Es importante destacar que estas ponderaciones no tienen nada que ver con la inferencia causal, las puertas traseras ni la ponderación de probabilidad inversa; su único propósito es ayudar a reducir la escala del desequilibrio que surge del sobreemparejamiento. Si se utiliza replace = FALSEy aplica el emparejamiento 1:1, toda la weightscolumna será solo 1.
Podemos incorporar esos pesos al modelo y obtener una estimación más precisa:
model_matched_wts <- lm(malaria_risk ~ net,
data = matched_data_for_real,
weights = weights)
tidy(model_matched_wts)# A tibble: 2 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 36.1 0.597 60.4 0
2 netTRUE -10.5 0.765 -13.7 1.17e-39Ponderación de probabilidad inversa
Una posible desventaja de la comparación es que generalmente hay que descartar una parte considerable de los datos: todo lo que no tenga comparación no se incluye en los datos coincidentes finales.
Un enfoque alternativo al emparejamiento es asignar a cada observación cierta probabilidad de recibir tratamiento y luego ponderar cada observación por su probabilidad inversa: las observaciones que se predice que recibirán tratamiento y luego no lo reciben, o las observaciones que se predice que no recibirán tratamiento y luego lo reciben, recibirán más peso que las observaciones que reciben o no reciben tratamiento según lo previsto.
La generación de estos pesos de probabilidad inversa requiere un proceso de dos pasos: (1) primero generamos los puntajes de propensión, o la probabilidad de recibir tratamiento, y luego (2) utilizamos una fórmula especial para convertir dichos puntajes de propensión en pesos. Una vez que tenemos los pesos de probabilidad inversa, podemos incorporarlos a nuestro modelo de regresión.
Curso intensivo simplificado de regresión logística
Hay muchas maneras de generar puntajes de propensión (como regresión logística, regresión probit e incluso técnicas de aprendizaje automático como bosques aleatorios y redes neuronales), pero la regresión logística es probablemente el método más común.
Los detalles técnicos completos de la regresión logística están más allá del alcance de esta clase, pero si tienes curiosidad, deberías consultar [este tutorial muy accesible]()https://uc-r.github.io/logistic_regression .
Lo único que necesitas saber es que la variable de resultado en los modelos de regresión logística debe ser binaria, y las variables explicativas que incluya en el modelo ayudan a explicar la variación en la probabilidad de su resultado binario. La Y (o resultado) en la regresión logística es un cociente logarítmico de probabilidades, que obliga a que el resultado del modelo esté en un rango de 0 a 1:
\[ \log \frac{p_y}{p_{1-y}} = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \dots + \beta_n x_n + \epsilon \]
Así se ve visualmente. Dado que el uso de la red es un resultado binario, existen líneas de observaciones en 0 y 1 a lo largo del eje y. La línea curva azul muestra el resultado de un modelo de regresión logística: las personas con bajos ingresos tienen una baja probabilidad de usar la red, mientras que las personas con altos ingresos tienen una probabilidad mucho mayor.
También incluí una línea roja que muestra los resultados de un modelo MCO clásico lm(). Sigue la línea azul bastante bien durante un tiempo, pero predice probabilidades negativas si se usan valores de ingresos más bajos, como menos de 250. Por extrañas razones históricas y matemáticas, muchos economistas prefieren usar MCO con resultados binarios (incluso tienen un nombre elegante: modelos de probabilidad lineal (MPL)), pero yo prefiero la regresión logística, ya que no genera probabilidades mayores del 100% ni menores del 0%. (PERO NO TE QUEJES NUNCA DE LOS MPL EN INTERNET. ¡Provocarás peleas entre economistas y otros científicos sociales! 🤣)
ggplot(nets, aes(x = income, y = net_num)) +
geom_point(alpha = 0.1) +
geom_smooth(method = "lm", color = "red", size = 0.5) +
geom_smooth(method = "glm", method.args = list(family = binomial(link = "logit"))) +
labs(x = "Income", y = "Probability of using a net")
Los coeficientes de un modelo de regresión logística se interpretan de forma diferente a la habitual (¡y sus interpretaciones pueden ser controvertidas!). Aquí tienes un ejemplo del modelo del gráfico anterior:
# Notice how we use glm() instead of lm(). The "g" stands for "generalized"
# linear model. We have to specify a family in any glm() model. You can
# technically run a regular OLS model (like you do with lm()) if you use
# glm(y ~ x1 + x2, family = gaussian(link = "identity")), but people rarely do that.
#
# To use logistic regression, you have to specify a binomial/logit family like so:
# family = binomial(link = "logit")
model_logit <- glm(net ~ income + temperature + health,
data = nets,
family = binomial(link = "logit"))
tidy(model_logit)# A tibble: 4 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) -1.32 0.376 -3.50 0.000464
2 income 0.00209 0.000421 4.95 0.000000727
3 temperature -0.0589 0.0125 -4.70 0.00000264
4 health 0.00688 0.00430 1.60 0.109 Los coeficientes aquí no son números normales; se llaman “logaritmos de probabilidades” y representan el cambio en los logaritmos de probabilidades al aumentar las variables explicativas. Por ejemplo, aquí se muestran los logaritmos de probabilidades de usar un aumento neto de 0,00688 por cada punto de aumento en la puntuación de salud. Pero ¿qué significan los logaritmos de probabilidades? Nadie lo sabe.
Puedes hacer que estos coeficientes sean un poco más interpretables al deslogarlos y crear algo llamado “odds ratio”. Estos coeficientes se registraron con un logaritmo natural, por lo que se deslogarizan elevandoa la potencia del coeficiente. La razón de probabilidades para la temperatura es, \(e^{-0.0589}\)o 0,94. Las razones de probabilidades se interpretan de forma ligeramente diferente a los coeficientes del modelo convencional. Todas las razones de probabilidades se centran en 1: valores superiores a 1 indican un aumento en la probabilidad del resultado, mientras que valores inferiores a 1 indican una disminución en la probabilidad del resultado.
Nuestra razón de probabilidades de temperatura nocturna es de 0,94, que es 0,06 por debajo de 1. Esto significa que, por cada punto de aumento en la temperatura nocturna, una persona tiene un 6 % menos de probabilidad de usar un mosquitero. Si el coeficiente fuera de 1,34, podríamos decir que tendría un 34 % más de probabilidad de usar un mosquitero; si fuera de 5,02, podríamos decir que tendría un 5 % más de probabilidad de usar un mosquitero; si fuera de 0,1, podríamos decir que tendría un 90 % menos de probabilidad de usar un mosquitero.
Puedes hacer que R exponencie los coeficientes automáticamente incluyéndolos exponentiate = TRUEen tidy():
tidy(model_logit, exponentiate = TRUE)# A tibble: 4 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 0.268 0.376 -3.50 0.000464
2 income 1.00 0.000421 4.95 0.000000727
3 temperature 0.943 0.0125 -4.70 0.00000264
4 health 1.01 0.00430 1.60 0.109 Pero, una vez más, ¡esto va más allá del alcance de esta clase! Solo recuerda que, al construir un modelo de regresión logística, estás usando variables explicativas para predecir la probabilidad de un resultado.
Una última explicación breve de por qué tenemos que usar logaritmos de probabilidades. Al trabajar con resultados binarios, trabajamos con probabilidades, y podemos crear algo llamado “odds” con probabilidades. Si hay un 70 % de probabilidad de que algo ocurra, hay un 30 % de probabilidad de que no ocurra. La razón entre esas dos probabilidades se llama “odds”:Las probabilidades suelen seguir una relación curva: a medida que se asciende a niveles más altos de la variable explicativa, las probabilidades aumentan más rápidamente. Sin embargo, si se registran estas probabilidades (), la relación se vuelve lineal, lo que significa que podemos usar la regresión lineal tradicional para calcular probabilidades. ¡Mágico!
Puede ver la relación entre las probabilidades logarítmicas y las razones de probabilidades en los primeros dos paneles aquí (estos son datos falsos donde X varía entre -5 y 5, e Y es 0 o 1; puede ver los puntos de datos en el panel final como puntos):
Los coeficientes de la regresión logística son probabilidades logarítmicas porque provienen del modelo que crea esa línea recta en el primer panel. Las probabilidades logarítmicas son imposibles de interpretar, por lo que podemos deslogaritmarlas (Los coeficientes de la regresión logística son probabilidades logarítmicas porque provienen del modelo que crea esa línea recta en el primer panel. Las probabilidades logarítmicas son imposibles de interpretar, por lo que podemos deslogaritmarlas para convertirlos en razones de probabilidades.
El panel inferior muestra las probabilidades predichas. Puedes realizar una transformación matemática más con las probabilidades (\(\frac{p}{1 - p}\)) para generar una probabilidad en lugar de probabilidades: \(\frac{\text{odds}}{1 + \text{odds}}\). Eso es lo que es una puntuación de propensión.para convertirlos en razones de probabilidades.
Paso 1: Generar puntuaciones de propensión
Con esa pequeña incursión en la regresión logística hecha, podemos construir un modelo para generar puntuaciones de propensión (o probabilidades predichas).
Al incluir variables en el modelo que genera puntuaciones de propensión, realizamos ajustes y cerramos las puertas traseras en el DAG, al igual que con el emparejamiento. Pero, a diferencia del emparejamiento, ¡no descartamos ningún dato! Simplemente, damos más importancia a algunas observaciones y menos a otras.
Primero construimos un modelo que predice el uso de la red en función de los ingresos, las temperaturas nocturnas y la salud (ya que esos nodos son nuestros factores de confusión del DAG):
model_net <- glm(net ~ income + temperature + health,
data = nets,
family = binomial(link = "logit"))
# We could look at these results if we wanted, but we don't need to for this class
# tidy(model_net, exponentiate = TRUE)Luego podemos ingresar los ingresos, las temperaturas y la salud para cada fila de nuestro conjunto de datos y generar una probabilidad prevista utilizando este modelo:
# augment_columns() handles the plugging in of values. You need to feed it the
# name of the model and the name of the dataset you want to add the predictions
# to. The type.predict = "response" argument makes it so the predictions are in
# the 0-1 scale. If you don't include that, you'll get predictions in an
# uninterpretable log odds scale.
net_probabilities <- augment_columns(model_net,
nets,
type.predict = "response") %>%
# The predictions are in a column named ".fitted", so we rename it here
rename(propensity = .fitted)
# Look at the first few rows of a few columns
net_probabilities %>%
select(id, net, income, temperature, health, propensity) %>%
head()# A tibble: 6 × 6
id net income temperature health propensity
<dbl> <lgl> <dbl> <dbl> <dbl> <dbl>
1 1 TRUE 781 21.1 56 0.367
2 2 FALSE 974 26.5 57 0.389
3 3 FALSE 502 25.6 15 0.158
4 4 TRUE 671 21.3 20 0.263
5 5 FALSE 728 19.2 17 0.308
6 6 FALSE 1050 25.3 48 0.429Los puntajes de propensión se muestran en la propensitycolumna. Es poco probable que algunas personas, como la persona 3, usen mosquiteros (solo un 15,8 % de probabilidad) dados sus niveles de ingresos, temperatura y salud. Otras, como la persona 6, tienen una mayor probabilidad (42,9 %), ya que sus ingresos y salud son más altos. ¡Genial!
A continuación, necesitamos convertir esos puntajes de propensión en ponderaciones de probabilidad inversas, lo que aumenta la importancia de las observaciones inusuales (es decir, personas que tenían una alta probabilidad de usar una red pero no lo hicieron, y viceversa). Para ello, seguimos esta ecuación:
\[ \frac{\text{Treatment}}{\text{Propensity}} + \frac{1 - \text{Treatment}}{1 - \text{Propensity}} \]
Esta ecuación generará ponderaciones que proporcionan el efecto promedio del tratamiento (EPT), pero existen otras versiones que permiten calcular el efecto promedio del tratamiento en los sujetos tratados (EPT), el efecto promedio del tratamiento en los sujetos de control (ETC) y muchos otros.
Usaremos mutate()para crear una columna para el peso de probabilidad inversa:
net_ipw <- net_probabilities %>%
mutate(ipw = (net_num / propensity) + ((1 - net_num) / (1 - propensity)))
# Look at the first few rows of a few columns
net_ipw %>%
select(id, net, income, temperature, health, propensity, ipw) %>%
head()# A tibble: 6 × 7
id net income temperature health propensity ipw
<dbl> <lgl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 1 TRUE 781 21.1 56 0.367 2.72
2 2 FALSE 974 26.5 57 0.389 1.64
3 3 FALSE 502 25.6 15 0.158 1.19
4 4 TRUE 671 21.3 20 0.263 3.81
5 5 FALSE 728 19.2 17 0.308 1.44
6 6 FALSE 1050 25.3 48 0.429 1.75Estas primeras filas tienen ponderaciones bastante bajas: quienes tenían baja probabilidad de usar redes no las usaron, mientras que quienes tenían alta probabilidad sí. ¡Pero observen a la persona 4! Solo tenía un 26 % de probabilidad de usar una red, ¡y la usó ! ¡Qué raro! Por lo tanto, tiene una ponderación de probabilidad inversa más alta (3,81).
Paso 2: Estimación
Ahora que hemos generado ponderaciones de probabilidad inversa basadas en nuestros factores de confusión, podemos ejecutar un modelo para determinar el efecto causal del uso de mosquiteros en el riesgo de malaria. Nuevamente, no necesitamos incluir ingresos, temperaturas ni salud en el modelo, ya que los usamos al crear las puntuaciones de propensión y las ponderaciones:
model_ipw <- lm(malaria_risk ~ net,
data = net_ipw,
weights = ipw)
tidy(model_ipw)# A tibble: 2 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 39.7 0.468 84.7 0
2 netTRUE -10.1 0.658 -15.4 3.21e-50Después de usar las ponderaciones de probabilidad inversas, encontramos un efecto causal de -10,13 puntos. Está un poquito lejos del valor real de 10, ¡pero no está nada mal!
Es importante verificar los valores de las ponderaciones de probabilidad inversa. A veces pueden ser demasiado altas, como si alguien tuviera ingresos de 0 y el nivel de salud más bajo posible, viviera en un campo nevado y aun así usara una red. Tener valores de IPW muy altos puede distorsionar la estimación. Para solucionar esto, podemos truncar las ponderaciones a un nivel inferior. No existe una regla general para una buena ponderación máxima; he visto usar 10 con frecuencia. En nuestros datos de redes antimosquitos, ninguna IPW es superior a 10 (el máximo es exactamente 10 con la persona 877), por lo que no debemos preocuparnos por el truncamiento.
Si quisieras truncar, harías algo como esto (aquí estamos truncando en 8 en lugar de 10, por lo que el truncamiento realmente hace algo):
net_ipw <- net_ipw %>%
# If the IPW is larger than 8, make it 8, otherwise use the current IPW
mutate(ipw_truncated = ifelse(ipw > 8, 8, ipw))
model_ipw_truncated <- lm(malaria_risk ~ net,
data = net_ipw,
weights = ipw_truncated)
tidy(model_ipw_truncated)# A tibble: 2 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 39.7 0.467 85.1 0
2 netTRUE -10.2 0.656 -15.5 4.32e-51Ahora el efecto causal es -10,19, que es ligeramente menor y probablemente menos preciso ya que realmente no tenemos ningún caso excepcional que altere nuestra estimación original de IPW.
Resultados de todos los modelos
Solo usamos datos observacionales para estimar los efectos causales. ¡Causalidad sin RCTs!
Comparemos todos los ATE que acabamos de calcular:
modelsummary(list("Ingenuo" = model_wrong,
"Ponderado" = model_matched, "Ponderado + pesos" = model_matched_wts,
"IPW" = model_ipw, "IPW truncado en 8" = model_ipw_truncated))| Ingenuo | Ponderado | Ponderado + pesos | IPW | IPW truncado en 8 | |
|---|---|---|---|---|---|
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 | |||||
| (Intercept) | 41.937*** | 38.335*** | 36.090*** | 39.679*** | 39.679*** |
| (0.405) | (0.599) | (0.597) | (0.468) | (0.467) | |
| netTRUE | -16.332*** | -12.730*** | -10.485*** | -10.131*** | -10.190*** |
| (0.649) | (0.768) | (0.765) | (0.658) | (0.656) | |
| Num.Obs. | 1752 | 1117 | 1117 | 1752 | 1752 |
| R2 | 0.265 | 0.198 | 0.144 | 0.119 | 0.121 |
| R2 Adj. | 0.265 | 0.197 | 0.144 | 0.119 | 0.121 |
¿Cuál es correcto? En este caso, dado que se trata de datos simulados falsos donde construí un efecto de 10 puntos, podemos ver cuál de estos modelos se acerca más: aquí, el modelo IPW no truncado gana. Pero no siempre será así. Tanto la coincidencia como el IPW funcionan bien para cerrar puertas traseras y ajustar factores de confusión. En la vida real, no se conoce el valor real, así que se recomienda probar varias maneras.