Inferência Causal

Valor-p de Rosenbaum

Prof. Carlos Trucíos
ctrucios@unicamp.br

Instituto de Matemática, Estatística e Computação Científica (IMECC),
Universidade Estadual de Campinas (UNICAMP).

Introdução

Rosenbaum (1987) introduziu uma técnica de análise de sensibilidade para estudos obervacionais de tipo Matching.
O método funciona para matching geral, mas por simplicade desenvolveremos o método para matching 1-1.
O método é utilizado quando estamos interessados em testar a hipótese nula forte.

Sensibilidade e Matching

Considere um estudo observacional do tipo matching, em que o índice \((i, j)\) representa a unidade \(j\) no par \(i\) (\(i = 1, \cdots, n\) e \(j = 1, 2\)).

Com matching perfeito, as unidades \((i, 1)\) e \((i, 2)\) tem o mesmo valor de \(X_i\).

Consideremos a versão estendida do propensity score, isto é \[e_{ij} = P(Z_{ij} = 1 | X_i, Y_{ij}(1), Y_{ij}(0))\]

Sensibilidade e Matching

Assumindo ignorabilidade, \(e_{ij} = P(Z_{ij} = 1 | X_i, Y_{ij}(1), Y_{ij}(0)) = P(Z_{ij} = 1 | X_i)\).
- Assim, \(e_{i1} = e_{i2}\) e \(\pi_{i1} = 1/2\).
- Isto implica que, condicionado em \(X\) e \(\mathbb{S}\), o mecanismo de atribuição de tratamento é o mesmo do MPE e podemos utilizar o que aprendimos anteriormente.
Contudo, sem assumirmos ignorabilidade, \(e_{ij}\) é tambem uma função dos resultados potenciais.
- Rousenbaum (1987) propõe uma análise de sensibilidade limitando (por cima e por baixo) a razão de chances.

Sensibilidade e Matching

Se denotarmos \(o_{ij} = e_{ij}/(1 - e_{ij})\), a razão de chances é dada por \[\dfrac{o_{i1}}{o_{i2}} = \dfrac{e_{i1}/ (1 - e_{i1})}{e_{i2} / (1 - e_{i12})} = \dfrac{\dfrac{P(Z_{i1} = 1 |X_i, Y_{i1}(1), Y_{i1}(0))}{P(Z_{i1} = 0 |X_i, Y_{i1}(1), Y_{i1}(0))}}{\dfrac{P(Z_{i2} = 1 |X_i, Y_{i2}(1), Y_{i2}(0))}{P(Z_{i2} = 0 |X_i, Y_{i2}(1), Y_{i2}(0))}}.\]

Sensibilidade e Matching

Suposição de Rosenbaum

Para algum \(\Gamma \geq 1\) (pre-especificado), a razão de chances é limitada por \[o_{i1}/o_{i2} \leq \Gamma \quad e \quad o_{i2} / o_{i1} \leq \Gamma.\] Equivalentemente, \[\dfrac{1}{1 + \Gamma} \leq \pi_{i1} \leq \dfrac{\Gamma}{1 + \Gamma},\] em que \(\pi_{i1} = o_{i1} / (o_{i1} + o_{i2})\)

Se \(\Gamma = 1\), \(\pi_{ij} = 1/2\) e temos o MPE já estudado.

Rosenbaum

Consideremos a hipósete nula forte, ou seja, \[H_{0F}: Y_{ij}(1) = Y_{ij}(0), \quad i = 1, \cdots, n \text{ e } j = 1, 2.\]

Consideremos a seguinte classe de estatísticas de teste: \[T = \displaystyle \sum_{i = 1}^n S_i q_i,\] em que \(S_i = \mathbb{I}(\hat{\tau}_i > 0)\), \(\hat{\tau}_i = (2Z_{i1} - 1)(Y_{i1} - Y_{i2})\) e \(q_i \geq 0\) uma função de \((|\hat{\tau}_1|, \cdots, |\hat{\tau}_n|).\)

Rosenbaum

\[T = \displaystyle \sum_{i = 1}^n S_i q_i.\]

Alguns casos particulares são:

Estatística de sinais (\(T = \displaystyle \sum_{i = 1}^n S_i\))
Estatística T pareada (\(T = \displaystyle \sum_{i = 1}^n S_i |\hat{\tau}_i|\))
Estatística de Wilcoxon (\(T = \displaystyle \sum_{i = 1}^n S_i R_i\))

Rosenbaum

Qual a distribuição de \(T\) (sob \(H_{0F}\) e a suposição de Rosenbaum)

Más notícias

A distribuição de \(T\) é bastante complicada.

Boas notícias

Felizmente apenas precisamos da distribuição no caso pior, o caso que fornece uma distribuição de \(T\) que levará ao maior p-valor sob a suposição de Rosenbaum

Rosenbaum

No caso pior, \[S_i \sim Bernoulli(\Gamma / (1 + \Gamma)).\]

Assim, os primeiros momentos de \(T\) são dados por:

\[\mathbb{E}(T) = \dfrac{\Gamma}{1 + \Gamma} \displaystyle \sum_{i = 1}^n q_i \quad e \quad \mathbb{V}(T) = \dfrac{\Gamma}{(1 + \Gamma)^2} \displaystyle \sum_{i = 1}^n q_i^2.\]

Pelo TCL,

\[\dfrac{T -\mathbb{E}(T)}{\sqrt{\mathbb{V}(T)}} \sim N(0, 1).\]

Exemplos

Exemplo 1

Lalonde (1986) está interessado no efeito causal de um programa de treinamento sob o salario. Utilizaremos os dados referentes ao estudo observacional (cps1re74.csv).

Code

library(Matching)
library(dplyr)
dados <- read.table("datasets/cps1re74.csv", head = TRUE)
glimpse(dados)

Rows: 16,177
Columns: 10
$ treat    <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
$ age      <int> 37, 22, 30, 27, 33, 22, 23, 32, 22, 33, 19, 21, 18, 27, 17, 1…
$ educ     <int> 11, 9, 12, 11, 8, 9, 12, 11, 16, 12, 9, 13, 8, 10, 7, 10, 13,…
$ black    <int> 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
$ hispan   <int> 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ married  <int> 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0…
$ nodegree <int> 1, 1, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 1, 0, 0, 1…
$ re74     <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ re75     <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ re78     <dbl> 9930.0460, 3595.8940, 24909.4500, 7506.1460, 289.7899, 4056.4…

Exemplos

       ytreated ycontrol
[1,]  9930.0460        0
[2,]  9930.0460        0
[3,]  3595.8940     7810
[4,] 24909.4500     1900
[5,]  7506.1460        0
[6,]   289.7899      649

Exemplos

\[\text{Utilizaremos } \quad T = \displaystyle \sum_{i = 1}^n S_i |\hat{\tau}_i|.\]

Se \(\Gamma =1\), é possível simular a distribuição de \(T\) e obsermos o p-valor
Se \(\Gamma > 1\) utilizamos o caso pior e obtemos o o-valor

Code

library(sensitivitymw)
Gamma <- seq(1, 1.4, 0.001)
pvalor <- Gamma
for (i in 1:length(pvalor)) {
  pvalor[i] <- senmw(data_match, gamma = Gamma[i], method = "t")$pval
}

Exemplos

Code

plot(Gamma, pvalor, type = "l")

A medida que \(\Gamma\) aumenta, o p-valor cresce, sendo \(\Gamma = 1.232\) o valor máximo com que ainda podemos rejeitar \(H_0\) para um nível de significância de 0.05.

Code

Gamma[tail(which(pvalor < 0.05), 1)]

[1] 1.232

Exemplos

Code

gama <- Gamma[25]
tau_i <- data_match[, "ytreated"] - data_match[, "ycontrol"]
qi <- abs(tau_i)
mu <- gama / (1 + gama) * sum(qi)
sigma <- sqrt(gama / (1 + gama)^2 * sum(qi^2))
t <-  sum((tau_i >0) * qi)
1 - pnorm((t - mu)/sigma)

[1] 0.003898971

Code

pvalor[25]

[1] 0.003898971

Exemplos

Exemplo 2

O conjunto de dados erpcp do pacote sensitivitymw contém 39 observações matched por pares considerando as covariáveis idade e fumante.

Code

data(erpcp)
glimpse(erpcp)

Rows: 39
Columns: 2
$ welder  <dbl> 0.899, 1.233, 1.733, 3.156, 1.749, 0.431, 1.430, 1.283, 1.137,…
$ control <dbl> 0.751, 0.875, 0.161, 0.630, 1.462, 0.702, 0.936, 0.567, 0.726,…

Exemplos

Code

Gamma <- seq(1, 5, 0.005)
pvalor <- Gamma
for (i in 1:length(Gamma)) {
  pvalor[i] <- senmw(erpcp, gamma = Gamma[i], method = "t")$pval
}
Gamma[tail(which(pvalor < 0.05), 1)]

[1] 3.8

Exemplos

Code

plot(Gamma, pvalor, type = "l")

Referências

Peng Ding (2023). A First Course in Causal Inference. Capítulo 19.