Randomization Inference en R: una mejor forma de calcular p-values en experimentos aleatorios

Paso 1: escoger una hipótesis nula “sharp”

En un test de hipótesis estándar tenemos una hipótesis nula que afirma algo respecto a un parámetro poblacional, por ejemplo, \(\beta_1=0\) donde \(\beta_1\) es el coeficiente asociado a la variable de tratamiento \(D\). Lo que esta hipótesis nula plantea es que el efecto promedio del tratamiento es 0, pero no dice nada acerca de los efectos del tratamiento para cada unidad (\(\delta_i\)).

En cambio, al aplicar RI utilizamos una hipótesis nula sharp, es decir, una nula que afirma algo sobre los efectos en cada una de las unidades.

La más usada y conocida es la sharp null de Fischer:

\[ \delta_i=0 \]

Que nos plantea que el tratamiento tiene efecto cero para todas las unidades analizadas.

Para efectos de simplicidad, en este post se usará la sharp null de Fischer, pero notar que podríamos perfectamente testear hipótesis alternativas tales como \(\delta_i=1\), \(\delta_i=2\), etc. Incluso podríamos testear hipótesis en las que diferentes grupos o unidades tengan valores de \(\delta_i\) distintos. Lo único que se requiere es tener una hipótesis que plantee algo respecto de cada uno de los \(\delta_i\).

¿Qué ganamos con usar una nula de este tipo? Pues que nos permite completar las columnas de Potential Outcomes en nuestro dataset.

Los datos que observamos normalmente tienen esta estructura:

library(tidyverse)
library(magrittr)
library(causaldata)

ri <- causaldata::ri %>% 
  mutate(id_unit = row_number(),
         y0 = as.numeric(y0),
         y1 = as.numeric(y1))

ri

## # A tibble: 8 x 6
##   name       d     y    y0    y1 id_unit
##   <chr>  <dbl> <dbl> <dbl> <dbl>   <int>
## 1 Andy       1    10    NA    10       1
## 2 Ben        1     5    NA     5       2
## 3 Chad       1    16    NA    16       3
## 4 Daniel     1     3    NA     3       4
## 5 Edith      0     5     5    NA       5
## 6 Frank      0     7     7    NA       6
## 7 George     0     8     8    NA       7
## 8 Hank       0    10    10    NA       8

Para todas las unidades observamos \(D_i\) e \(Y_i\), pero solo uno de los outcomes potenciales (\(Y_i^0\) o \(Y_i^1\)). Apalancándonos en la sharp null, podemos completar las columnas y0 e y1 en base a \(Y_i=Y_i^0=Y_i^1\)

ri_fischer_null <- ri %>% 
  mutate(y0 = y,
         y1 = y)

ri_fischer_null

## # A tibble: 8 x 6
##   name       d     y    y0    y1 id_unit
##   <chr>  <dbl> <dbl> <dbl> <dbl>   <int>
## 1 Andy       1    10    10    10       1
## 2 Ben        1     5     5     5       2
## 3 Chad       1    16    16    16       3
## 4 Daniel     1     3     3     3       4
## 5 Edith      0     5     5     5       5
## 6 Frank      0     7     7     7       6
## 7 George     0     8     8     8       7
## 8 Hank       0    10    10    10       8

Como estamos asumiendo que el tratamiento tiene un efecto igual a cero para todas las unidades, ambos outcomes potenciales son iguales al outcome observado.

Por supuesto, no estamos afirmando que esto sea así realmente, simplemente queremos descubrir qué tan improbable es haber observado las diferencias entre grupos que estamos observando bajo el supuesto de esta sharp null.

Paso 2: escoger un estadístico de test

Para realizar un test de hipótesis necesitamos no solo una hipótesis nula sino que también un estadístico. En un test de hipótesis tradicional, sería necesario que este estadístico cuente con un estimador para su varianza, pero al aplicar randomization inference podemos usar literalmente cualquier estadístico escalar que pueda obtenerse a partir de un vector de asignaciones de tratamiento (\(D\)) y de un vector de outcomes (\(Y\)).

Uno de los estadísticos más simples que podemos usar es la diferencia simple de medias (SDO) entre ambos grupos.

# Funcion que define el estadístico: recibe como input un dataframe con 
# las columnas `y` (outcome) y `d` (asignacion de tratamiento), y retorna 
# un valor escalar
sdo <- function(data) {
  data %>%
    summarise(te1 = mean(y[d == 1]),
              te0 = mean(y[d == 0]),
              sdo = te1 - te0) %>%
    pull(sdo)
}

sdo(ri_fischer_null)

## [1] 1

Sin embargo, nada impide usar otros estadísticos de test. Esta libertad es de hecho una de las principales ventajas de aplicar Randomization Inference.

Algunos ejemplos de estadísticos alternativos que podríamos usar son:

• Estadísticos de quantiles, por ejemplo la diferencia de medianas entre grupos (o cualquier otro percentil). Útiles si tenemos outliers u observaciones de alto leverage.

• Estadísticos de ranking, que convierten la variable de outcomes (\(Y\)) en valores de ranking (1 para el outcome más bajo, 2 para el segundo más bajo, etc) y luego computan una métrica en base a esos rankings (por ejemplo, diferencia de medias o de medianas).

• Estadístico KS (Kolmogorov-Smirnov), que permite identificar diferencias en las distribuciones de outcomes mediante medir la distancia máxima entre ambas funciones de distribución acumuladas. Este estadístico nos permitiría, por ejemplo, detectar el efecto de un tratamiento que no afecte la media de \(Y\), pero sí su varianza o dispersión.

Para mantener la simplicidad, a continuación se usará el SDO como estadístico de test, pero más abajo mostraré ejemplos de código utilizando el estadístico KS.

Paso 3: simular distintas asignaciones de tratamiento y obtener la distribución del estadístico

El siguiente consiste en obtener la distribución de valores que podría tomar el estadístico del test bajo la sharp null.

Para esto:

1. Realizaremos permutaciones en el vector \(D\) mediante el mismo proceso aleatorio que se utilizó para obtener las asignaciones de tratamiento originales en primer lugar. Por ejemplo, si para asignar el tratamiento se usó un proceso equivalente a lanzar una moneda al aire para cada observación (\(B(n, 0.5)\)), entonces cada una de las permutaciones deben generarse con el mismo proceso.³

   En nuestro ejemplo, el proceso consiste en sortear 4 asignaciones de tratamiento para un total de 8 observaciones, lo cual produce 70 permutaciones posibles en total. El siguiente código obtiene tales permutaciones:

   perms <- t(combn(ri_fischer_null$id_unit, 4)) %>% asplit(1) 
   
   perms_df <- 
     tibble(treated_units = perms) %>% 
     transmute(id_perm = row_number(),
               treated_units = map(treated_units, unlist))
   
   ri_permuted <- 
     crossing(perms_df, ri_fischer_null) %>% 
     mutate(d = map2_dbl(id_unit, treated_units, ~.x %in% .y))
   
   ri_permuted

   ## # A tibble: 560 x 8
   ##    id_perm treated_units name       d     y    y0    y1 id_unit
   ##      <int> <list>        <chr>  <dbl> <dbl> <dbl> <dbl>   <int>
   ##  1       1 <int [4]>     Andy       1    10    10    10       1
   ##  2       1 <int [4]>     Ben        1     5     5     5       2
   ##  3       1 <int [4]>     Chad       1    16    16    16       3
   ##  4       1 <int [4]>     Daniel     1     3     3     3       4
   ##  5       1 <int [4]>     Edith      0     5     5     5       5
   ##  6       1 <int [4]>     Frank      0     7     7     7       6
   ##  7       1 <int [4]>     George     0     8     8     8       7
   ##  8       1 <int [4]>     Hank       0    10    10    10       8
   ##  9       2 <int [4]>     Andy       1    10    10    10       1
   ## 10       2 <int [4]>     Ben        1     5     5     5       2
   ## # ... with 550 more rows

   ri_permuted %>% 
     pull(id_perm) %>% 
     n_distinct()

   ## [1] 70

2. Para cada permutación, simulamos los valores de \(Y\) en base a la hipótesis nula escogida previamente. Para la sharp null de Fischer literalmente no hay que hacer nada, ya que al plantear esta que \(Y_i=Y_i^0=Y_i^1\), el vector \(Y\) se mantiene constante en cada permutación de \(D\) (sin embargo, si nuestra sharp null fuera algo como \(\delta_i=1\) entonces sí habría que cambiar los valores de \(Y\) con cada permutación).

3. Calcularemos el valor del estadístico del test (el SDO en nuestro caso) para cada una de estas permutaciones y guardamos su valor.

   perms_stats <- 
     ri_permuted %>%
     group_by(id_perm) %>%
     summarise(te1 = mean(y[d == 1]),
               te0 = mean(y[d == 0]),
               sdo = te1 - te0)
   
   perms_stats %>% 
     select(id_perm, sdo)

   ## # A tibble: 70 x 2
   ##    id_perm   sdo
   ##      <int> <dbl>
   ##  1       1   1  
   ##  2       2   2  
   ##  3       3   3  
   ##  4       4   3.5
   ##  5       5   4.5
   ##  6       6  -4.5
   ##  7       7  -3.5
   ##  8       8  -3  
   ##  9       9  -2  
   ## 10      10  -2.5
   ## # ... with 60 more rows

4. Listo! Los valores del estadístico obtenidos en (3) nos indican la distribución de este bajo la sharp null.

   ggplot(perms_stats, aes(sdo)) +
     geom_histogram() +
     labs(x = "Estadístico del test (SDO)")

   

¿Cuántas permutaciones se deben realizar? Lo ideal sería agotar todas permutaciones posibles para conseguir una representación exacta de la distribución del estadístico de test bajo la sharp null. Sin embargo, esto solo es factible en datasets muy pequeños: incluso con un dataset de 2000 observaciones ya se vuelve computacionalmente prohibitivo obtener todas las combinaciones posibles del vector \(D\).

La alternativa en ese caso es conformarnos con hacer las permutaciones suficientes para obtener una aproximación razonable a la distribución del estadístico⁴.

Paso 4: comparar el estadístico “real” con la distribución simulada para obtener el p-value

Luego de haber obtenido la distribución (exacta o aproximada) del estadístico bajo la sharp null, procedemos a ubicar el estadístico “verdadero” (aquel correspondiente los valores reales de \(D\) e \(Y\)) dentro de esa distribución.

Notar que si queremos realizar un test de dos colas (lo más usual) debemos aplicar valor absoluto a la distribución de estadísticos (que es de hecho lo que haremos en este ejemplo).

true_statistic <- perms_stats %>% 
  filter(id_perm == 1) %>% 
  pull(sdo)

ggplot(perms_stats, 
       # Aplicamos valor absoluto porque haremos un test de 2 colas
       aes(abs(sdo))) +
  geom_histogram() +
  labs(x = "Valor absoluto del estadístico del test (SDO)",
       y = "count") +
  geom_vline(xintercept = true_statistic, colour = "red", size = 2) +
  annotate(geom = "label",
           x = true_statistic,
           y = 10,
           label = "Estadístico 'verdadero'",
           colour = "red") +
  annotate("segment", x = true_statistic+1, xend = true_statistic+3, y = 9, yend =9,
           colour = "purple", size = 2, arrow = arrow()) +
  annotate(geom = "label",
           x = true_statistic+2,
           y = 10,
           label = "Estadísticos más extremos",
           colour = "purple")

Habiendo hecho esto procedemos a calcular el p-value (que es lo que hemos estado buscando desde un comienzo).

Este p-value se calcula en dos pasos:

1. Se crea un ranking descendiente de los estadísticos de test obtenidos en las permutaciones (el con mayor valor tiene ranking 1 y el con menor valor tiene ranking N). Para los estadísticos “empatados” en un mismo valor, su ranking corresponde el máximo de los rankings originales de los estadísticos que comparten tal valor (por ejemplo, si tres estadísticos son iguales y tienen rankings 2, 3, y 4, todos ellos quedan con ranking igual a 4).

perms_ranked <- perms_stats %>%
  # Igual que antes, se aplica valor absoluto para tener un test de dos colas
  mutate(abs_sdo = abs(sdo)) %>% 
  select(id_perm, abs_sdo) %>% 
  arrange(desc(abs_sdo)) %>% 
  mutate(rank = row_number(desc(abs_sdo))) %>% 
  group_by(abs_sdo) %>% 
  mutate(new_rank = max(rank))

perms_ranked

## # A tibble: 70 x 4
## # Groups:   abs_sdo [11]
##    id_perm abs_sdo  rank new_rank
##      <int>   <dbl> <int>    <int>
##  1      25     6       1        2
##  2      46     6       2        2
##  3      24     5.5     3        4
##  4      47     5.5     4        4
##  5       5     4.5     5       12
##  6       6     4.5     6       12
##  7      22     4.5     7       12
##  8      23     4.5     8       12
##  9      48     4.5     9       12
## 10      49     4.5    10       12
## # ... with 60 more rows

2. Se calcula el ratio entre el ranking del estadístico verdadero y el total de permutaciones obtenidas (N).

n <- nrow(perms_ranked)

p_value <- 
  perms_ranked %>% 
  filter(id_perm == 1) %>%
  pull(new_rank)/n

p_value

## [1] 0.8571429

En nuestro caso, el p-value obtenido es de 0.86, lo que implica que no podemos rechazar la sharp null con un 5% de significancia (ni con ningún nivel de significancia convencional 😅). Dicho de otra forma, el valor de SDO observado no es un valor improbable de observar si asumimos que el tratamiento no tuvo efecto en ninguna unidad y que la única fuente de variabilidad en el estadístico es la asignación aleatoria de \(D\).

Bonus: usando el estadístico KS para medir diferencias en distribuciones de outcomes

Una de las ventajas de RI mencionadas al comienzo de este post era la posibilidad de utilizar estadísticos alternativos, siendo el único requisito sobre estos el que correspondan a valores escalares obtenidos a partir de \(D\) e \(Y\).

Uno de los estadísticos “alternativos” más interesantes (en mi opinión) es el estadístico KS (Kolmogorov-Smirnov) que mide la diferencia máxima entre dos funciones de distribución acumuladas empíricas (en este caso, entre la de las unidades tratadas \(\hat{F_T}\) y la de las unidades de control \(\hat{F_C}\)):

\[ T_{KS}=\max|\hat{F_T}(Y_i) - \hat{F_C}(Y_i)| \]

Lo cual se puede apreciar visualmente en el siguiente gráfico de ejemplo, donde el valor del estadístico corresponde al largo de la línea roja punteada.

Ejemplo de estadístico KS con 2 funciones de distribución acumuladas. Fuente: StackExchange

¿Por qué querríamos usar el estadístico KS en vez de algo más familiar, como la diferencia de medias?

Porque puede haber casos en que el tratamiento afecte la distribución de los outcomes, pero sin generar cambios importantes en la media, por ejemplo aumentando la dispersión, o generando una distribución bimodal. He aquí un ejemplo (tomando del Mixtape) con datos simulados:

tb <- tibble(
  d = c(rep(0, 20), rep(1, 20)),
  y = c(0.22, -0.87, -2.39, -1.79, 0.37, -1.54, 
        1.28, -0.31, -0.74, 1.72, 
        0.38, -0.17, -0.62, -1.10, 0.30, 
        0.15, 2.30, 0.19, -0.50, -0.9,
        -5.13, -2.19, 2.43, -3.83, 0.5, 
        -3.25, 4.32, 1.63, 5.18, -0.43, 
        7.11, 4.87, -3.10, -5.81, 3.76, 
        6.31, 2.58, 0.07, 5.76, 3.50)
)

kdensity_d1 <- tb %>%
  filter(d == 1) %>% 
  pull(y)
kdensity_d1 <- density(kdensity_d1)

kdensity_d0 <- tb %>%
  filter(d == 0) %>% 
  pull(y)
kdensity_d0 <- density(kdensity_d0)

kdensity_d0 <- tibble(x = kdensity_d0$x, y = kdensity_d0$y, d = 0)
kdensity_d1 <- tibble(x = kdensity_d1$x, y = kdensity_d1$y, d = 1)

kdensity <- full_join(kdensity_d1, kdensity_d0,
                      by = c("x", "y", "d"))
kdensity$d <- as_factor(kdensity$d)

ggplot(kdensity)+
  geom_point(size = 0.3, aes(x,y, color = d))+
  xlim(-7, 8)+
  scale_color_discrete(labels = c("Control", "Treatment"))

En el ejemplo de arriba, si hacemos randomization inference usando el SDO como estadístico no podemos descartar la sharp null de que el efecto del tratamiento es cero para todas las unidades, a pesar de que claramente las distribuciones son muy distintas.

set.seed(1989)

conduct_ri(
  test_function = sdo,
  assignment = "d",
  outcome = "y",
  declaration = declare_ra(N = 40, m = 20),
  sharp_hypothesis = 0,
  data = tb,
  sims = 10000
)

##                    term estimate two_tailed_p_value
## 1 Custom Test Statistic    1.415             0.1361

Veamos qué resultado obtenemos al usar el estadístico KS en vez del SDO:

# Declarando la función que obtiene el estadístico KS a partir de un dataframe
ks_statistic <- function(data) {
  
  control <- data[data$d == 0, ]$y
  treated <- data[data$d == 1, ]$y
  
  cdf1 <- ecdf(control) 
  cdf2 <- ecdf(treated)
  
  minMax <- seq(min(control, treated),
                max(control, treated),
                length.out=length(c(control, treated))) 
  
  x0 <-
    minMax[which(abs(cdf1(minMax) - cdf2(minMax)) == max(abs(cdf1(minMax) - cdf2(minMax))))]
  
  y0 <- cdf1(x0)
  y1 <- cdf2(x0) 
  
  diff <- unique(abs(y0 - y1))
  
  diff
  
}

# obteniendo el p-value con randomization inference
set.seed(1989)

ri_ks <-
  conduct_ri(
    test_function = ks_statistic,
    assignment = "d",
    outcome = "y",
    declaration = declare_ra(N = 40, m = 20),
    sharp_hypothesis = 0,
    data = tb,
    sims = 10000
  )

ri_ks

##                    term estimate two_tailed_p_value
## 1 Custom Test Statistic     0.45             0.0227

El estadístico KS sí permite rechazar la nula de efecto 0 a un nivel de significancia del 5%⁵.

En la siguiente figura podemos visualizar el estadístico KS correspondiente al vector \(D\) verdadero, junto con las CDF empíricas de los grupos de tratamiento y control.

tb2 <- tb %>%
group_by(d) %>%
arrange(y) %>%
mutate(rn = row_number()) %>%
ungroup()


cdf1 <- ecdf(tb2[tb2$d == 0, ]$y) 
cdf2 <- ecdf(tb2[tb2$d == 1, ]$y) 

minMax <- seq(min(tb2$y),
              max(tb2$y),
              length.out=length(tb2$y)) 

x0 <- 
  minMax[which(abs(cdf1(minMax) - cdf2(minMax)) == max(abs(cdf1(minMax) - cdf2(minMax))) )] 

y0 <- cdf1(x0) 
y1 <- cdf2(x0) 

ggplot(tb2) +
  geom_step(aes(x=y, y=rn, color=factor(d)),
            size = 1.5, stat="ecdf") +
  labs(y = "Densidad") +
  geom_segment(aes(x = x0[1], y = y0[1], xend = x0[1], yend = y1[1]),
               linetype = "dashed", color = "red", size = 1) +
  geom_point(aes(x = x0[1] , y= y0[1]), color="red", size=4) +
  geom_point(aes(x = x0[1] , y= y1[1]), color="red", size=4)