09-mas-hipotesis.Rmd

# Más de pruebas de hipótesis e intervalos

```{r setup, include=FALSE, message=FALSE}
library(tidyverse)
library(patchwork)
knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning=FALSE, fig.align = 'center', fig.height=3, out.width = '95%', cache = TRUE)
comma <- function(x) format(x, digits = 2, big.mark = ",")
theme_set(theme_minimal())
```

En esta sección veremos enfoques más clásicos para analizar una prueba de
hipótesis. En particular veremos situaciones donde podemos hacer algunos
supuestos teóricos acerca de la distribución de las poblaciones. Esta es una
sección complementaria para entender prácticas estadísticas usuales: recuerda
que discutimos antes que hacer estimación por intervalos generalmente es más
útil que hacer pruebas de hipótesis, y adicionalmente, tenemos también la
técnica de pruebas de permutaciones que podemos aplicar en muchos de los casos
que discutiremos a continuación.

El enfoque básico es el mismo que cuando vimos pruebas de permutaciones:
calculamos una estadística de prueba de los datos y luego, con una distribución
de referencia (asociada a la hipótesis nula), calculamos un valor-$p$. Si el
valor-$p$ es chico, entonces los resultados observados no pueden explicarse
fácilmente por variación muestral, y rechazamos la hipótesis nula.

Con esta idea básica, y *supuestos distribucionales acerca de las poblaciones*,
podemos construir pruebas que requieren menos cómputo. La desventaja es que hay
que checar con cuidado los supuestos distribucionales que hagamos. Si los
supuestos son incorrectos, las valores-$p$ no tienen mucho sentido y son
difíciles de interpretar.

Para esta sección seguiremos más a @Wasserman (Capítulo 10), pero puedes
revisar también @Chihara (Capítulo 8).

## Prueba de Wald {-}

Como hemos visto, existe normalidad asintótica en varios estimadores que
hemos considerado, como medias y proporciones muestrales.
También vimos que estimadores de máxima verosimilitud cumplen muchas
veces un teorema central del límite.

Así que supongamos que tenemos una estadística $\hat{\theta}_n$ que estima
$\theta$ y es asintóticamente insesgada y normal. Denotamos por
$\hat{\textsf{ee}}_n$ una estimación de su error estándar ---hay distintas maneras
de hacerlo: por ejemplo, con simulación *(bootstrap)*, o por medios analíticos
*(teoría)*. Recuerda que el error estándar de una estadística es *la desviación
estándar de su distribución de muestreo*.

Si nos interesa probar la hipótesis de que $\theta = 125$, por ejemplo,
y $\hat{\theta}_n$ es aproximadamente normal, entonces podemos construir una
distribución de referencia aproximada como sigue:

- Si la nula es cierta, entonces la distribución de muestreo de $\hat{\theta}_n$
es aproximadamente $\mathsf{N}(125, \hat{\textsf{ee}}_n)$.
- Esto implica que la siguiente estadística $\hat W_n$ es aproximadamente normal
estándar bajo la nula:

$$\hat W_n = \frac{\hat{\theta}_n - 125}{\hat{\textsf{ee}}_n} \sim \mathsf{N}(0,1)\,.$$
Por lo que valores muy fuera de $[-2,2]$, por ejemplo, dan evidencia en contra
de la hipótesis nula. Como $\hat W_n$ **no depende de ningún parámetro
desconocido**, podemos usarla como distribución de referencia para comparar el
valor de $W$ que obtuvimos en la muestra.

Si observamos para nuestra muestra un valor $\hat W_n = \hat w_n$ entonces, el
valor-$p$ (dos colas) de esta prueba es, aproximadamente,


$$\mathsf{valor-}p \approx P(|Z| > |\hat w_n|) = 2(1 - \Phi(|\hat w_n|))\,.$$
donde $Z\sim \mathsf{N}(0,1)$ y $\Phi$ es su función de distribución acumulada.


**Ejemplo: media muestral**.  La media nacional de las escuelas de enlace está alrededor
de 454 (matemáticas en 6o. grado). Tomamos una muestra de 180
escuelas del Estado de México, y queremos saber si la media obtenida
es consistente o no con la media nacional. Ya que estamos usando como
estimador una media de una muestra iid, podemos estimar el error estándar
de la media con

$$\hat{\textsf{ee}}_n = \frac{s}{\sqrt{n}}\,.$$
Obtenemos:

```{r}
set.seed(29)
muestra_edomex <- read_csv("data/enlace.csv") %>%
  filter(estado == "ESTADO DE MEXICO") %>%
  sample_n(180)
resumen <- muestra_edomex %>%
  summarise(media = mean(mate_6), s = sd(mate_6), n = n()) %>%
  mutate(ee = s / sqrt(n))
resumen
```

La hipótesis nula es que la media poblacional del Estado de México
es igual a 454. Calculamos el valor-$p$ usando la prueba de Wald:

```{r}
dif <- (resumen %>% pull(media)) - 454
ee <- resumen %>% pull(ee)
w <- dif / ee
p <- 2 * (1 - pnorm(abs(w)))
p
```
y vemos que esta muestra es consistente con la media nacional. No tenemos
evidencia en contra de que la media del estado de México es muy similar
a la nacional.

```{block, type="ejercicio"}
Repite esta prueba con una muestra de Chiapas. ¿Qué resultado obtienes?
```


Tenemos entonces:

```{block2, type='mathblock'}

**Prueba de Wald.** Consideramos probar la hipótesis nula $H_0: \theta =
\theta_0$ contra la alternativa $H_1: \theta \neq \theta_0$. Suponemos que
$\hat{\theta}_n$ es asintóticamente normal e insesgada, de modo que bajo la
hipótesis nula
$$\hat W_n = \frac{\hat{\theta}_n - \theta_0}{\hat{\textsf{ee}}_n} \overset{.}{\sim} \mathsf{N}(0,1)\,.$$
Entonces el valor-$p$ de la **prueba de Wald** para esta hipótesis nula es
$$\mathsf{valor-}p \approx P(|Z| > |\hat w_n|) = 2(1 - \Phi(|\hat w_n|))\,.$$

```


**Ejemplo.** Podemos hacer la prueba de Wald para proporciones con el
estimador usual $\hat{p}_n$ que estima una proporción poblacional $p$.
En este caso, utilizamos la estimación usual del error estándar de $\hat{p}_n$,
que está dada por
$$\hat{\textsf{ee}}_n = \sqrt{\frac{\hat{p}_n(1-\hat{p}_n)}{n}}\,.$$
Supongamos por ejemplo que en nuestros datos observamos que en $n=80$
muestras independientes, tenemos $S_n=47$ éxitos. ¿Es esto consistente
con la hipótesis nula $p = 0.5$?

Calcuamos primero:

```{r}
p_hat <- 47 / 80
ee <- sqrt(p_hat * (1 - p_hat) / 80)
```

y la estadística $\hat W_n$ de prueba es:

```{r}
w <- (p_hat - 0.5) / ee
w
```

Calculamos su valor-$p$:

```{r}
valor_p <- 2 * (1 - pnorm(abs(w)))
valor_p
```

Y vemos que en este caso tenemos evidencia baja de que la proporción poblacional
es distinta de 0.5.



## Observación: pruebas $t$ y práctica estadística {-}

Con más supuestos distribucionales podemos hacer otros tipos de pruebas donde
no requerimos hacer supuestos asintóticos. Por ejemplo, si suponemos
que la muestra obtenida $X_1,\ldots, X_n$ proviene de una distribución
normal $\mathsf{N}(\mu, \sigma)$ (cosa que es **necesario** verificar), entonces
es posible demostrar que la estadística

$$T_n = \frac{\bar{X}_n - \mu}{S_n / \sqrt{n}}\,,$$
tiene una distribución exacta que es $t$ de Student con $n-1$ grados de
libertad, y no depende de otros parámetros, de manera que podemos usarla como
distribución de referencia y podemos calcular valores-$p$ exactos (revisa la
sección 8.1 de @Chihara).

La diferencia con usar una prueba de Wald está en que aquí consideramos también
la variablidad del error estándar estimado, lo que correctamente sugiere que
esperamos variaciones proporcionalmente más grandes en $T_n$ comparado con lo que
sucede si no consideramos esta variación (como en la prueba de Wald). Sin
embargo:

- Si la muestra $n$ es grande, la distribución $t$ de Student con
$n-1$ grados de libertad es muy similar a la normal estándar, de manera que la
aproximación de Wald es apropiada.
- Cuando la muestra $n$ es chica, es difícil validar el supuesto de normalidad, a menos
que tengamos alguna información adicional acerca de la distribución poblacional.
- La prueba tiene cierta robustez a desviaciones de normalidad de las observaciones,
pero si el sesgo es muy grande, por ejemplo, el supuesto es incorrecto y da
valores $p$ distorsionados.

Puedes ver [aquí](https://en.wikipedia.org/wiki/Student%27s_t-distribution), o el apéndice
B.11 de @Chihara para ver descripciones de la distribución $t$ y cómo se compara
con una normal estándar dependiendo de los grados de libertad.

En muchas ocasiones, en la práctica es común no checar supuestos y saltar directamente
a hacer pruebas $t$, lo cual no es muy seguro.
Si tenemos duda de esos supuestos, podemos hacer
pruebas gráficas o de permutaciones, si son apropiadas.


## Prueba de Wald para dos medias o proporciones {-}

Cuando tenemos dos muestras extraidas de manera independiente de dos poblaciones
distintas, y queremos ver si la hipótesis de medias poblacionales
iguales es consistente con los datos, podemos usar también una prueba de Wald.

Sea $\bar{X}_{n_1}$ y $\bar{X}_{n_2}$ las medias muestrales correspondientes. Si
la hipótesis de normalidad aplica para ambas distribuciones muestrales
(normalidad asintótica), el estimador
$$\hat{\delta}_n = \bar{X}_{n_1} - \bar{X}_{n_2}$$\,,
es aproximadamente normal con media $\mathsf{N}(\mu_1 - \mu_2, \textsf{ee}),$
donde $\mu_1$ y $\mu_2$ son las medias poblacionales correspondientes; $n = n_1
+ n_2;$ y el error estándar de $\hat{\delta}$ es la raíz de la suma de los
cuadrados de los errores estándar de $\bar{X}$ y $\bar{Y}$:

$$ \textsf{ee} = \sqrt{\textsf{ee}_1^2 + \textsf{ee}_{2}^2}\,.$$
Se sigue entonces que:
$$\textsf{ee} =\sqrt{\frac{\sigma_1^2}{n_1}+\frac{\sigma_2^2}{n_2}}\,,$$
(Nota: usa probabilidad para explicar por qué es cierto esto). De esto
se deduce que bajo la hipótesis nula de igualdad de medias $\mu_1 = \mu_2$,
tenemos que la estadística de Wald

$$\hat W_n = \frac{\hat{\delta}_n - 0}{\sqrt{\frac{s_1^2}{n_1}+\frac{s_2^2}{n_2}} } \sim \mathsf{N}(0,1)\,,$$
es aproximamente normal estándar. Procedemos entonces a calcular el valor
$p$ usando la función de distribución acumulada de la normal estándar.


En el caso particular de proporciones, podemos simplificar, como hicimos arriba,
a
$$W = \frac{\hat{p}_1 - \hat{p}_{n_2}}{\sqrt{\frac{\hat{p}_1(1-\hat{p}_1)}{n_1}+\frac{\hat{p}_2(1-\hat{p}_2)}{n_2}} } \sim \mathsf{N}(0,1)\,.$$

```{block, type = 'ejercicio'}
Haz una prueba comparando las medias en enlace de la Ciudad de México vs
Estado de México. ¿Hay evidencia de que tienen distintas medias?
```


**Ejemplo (@Wasserman).** Supongamos tenemos dos conjuntos de prueba para
evaluar algoritmos de predicción, de tamaños $n_1=100$ y $n_2=250$
respectivamente, tenemos dos algoritmos para generar predicciones de clase
(digamos positivo y negativo). Usaremos el primer conjunto para evaluar el
algoritmo 1 y el segundo para evaluar el algoritmo 2. El algoritmo 1 corre en 1
hora, y el algoritmo 2 tarda 24 horas.

Supón que obtenemos que la tasa de clasificación correcta del primer algoritmo
es $\hat{p}_{n_1} = 0.85$, y la tasa del segundo es de $\hat{p}_{n_2} = 0.91$.
¿Estos datos son consistentes con la hipótesis de que los algoritmos tienen
desempeño muy similar? Es decir, queremos probar la hipótesis $p_1 = p_2$.

Calculamos la estadística de Wald:

```{r}
n_1 <- 100
n_2 <- 250
p_hat_1 <- 0.86
p_hat_2 <- 0.90
ee <- sqrt(p_hat_1 * (1 - p_hat_1) / n_1 + p_hat_2 * (1 - p_hat_2) / n_2)
delta = p_hat_1 - p_hat_2
w  <- delta / ee
w
```

que da un valor-$p$ de:

```{r}
2 * (1 - pnorm(abs(w)))
```

Y vemos que valor-$p$ es grande, de forma que los datos son consistentes
con la hipótesis de que los algoritmos tienen desempeño similar. ¿Cómo tomaríamos
nuestra decisión final? Si la diferencia entre 1 hora y 24 horas no es muy
importante, entonces preferíamos usar el algoritmo 2. Sin embargo, si el costo
de 24 horas es más alto que 1 hora de corrida, los datos no tienen indicios fuertes
de que vayamos a perder en desempeño, y podriamos seleccionar el algoritmo 1.


## Datos pareados {-}

Las pruebas que acabamos de ver para comparar medias requieren
poblaciones independientes. Si las dos muestras están pareadas (es decir,
son dos mediciones en una misma muestra), podemos
tomar considerar las diferencias $D_i = X_i - Y_i$ y utilizar la prueba para una
sola muestra con la media $\bar{D}_n$. Esta es una prueba de Wald pareada.


**Ejemplo (@Wasserman).** Ahora supongamos que utilizamos la misma
muestra de tamaño $n=300$ para probar los dos algoritmos. En este caso,
no debemos hacer la prueba para medias de muestras independientes. Sin embargo,
esto podemos ponerlo en términos de una prueba para una sola muestra.

Tenemos las observaciones $X_1,\ldots, X_n$ y $Y_1,\dots, Y_n$, donde
$X_i=1$ si el algoritmo 1 clasifica correctamente, y 0 en otro caso. Igualmente,
$Y_i=1$ si el algoritmo 2 clasifica correctamente, y 0 en otro caso. Definimos

$$D_i= X_i - Y_i\,.$$
Y $D_1,\ldots, D_n$ es una muestra iid. Ahora observemos que la media $\bar{D}_n$
tiene valor esperado $p_1 - p_2$, donde $p_1$ y $p_2$ son las tasas de correctos
del algoritmo 1 y del algoritmo 2 respectivamente. Podemos hacer una prueba
de Wald como al principio de la sección:

$$\hat W_n = \frac{\bar{D}_n - 0}{{\textsf{ee}}}\,.$$
Y notemos que el error estándar **no** se calcula como en el ejemplo anterior. Podríamos
usar bootstrap para estimarlo, pero en este caso podemos usar el estimador usual

$$\hat{\textsf{ee}}_n = \frac{s_n}{\sqrt{n}}\,,$$
donde
$$s_n = \frac{1}{n}\sum_{i=1}^n (D_i - \bar{D}_n)^2\,,$$
y nótese que necesitamos las decisiones indiviudales de cada algoritmo para
cada caso, en contraste al ejemplo anterior de muestras independientes donde
los errores estándar se calculaban de manera independiente. Esto tiene sentido,
pues la variablidad de $\bar{D}_n$ depende de cómo están correlacionados
los aciertos de los dos algoritmos.

Supongamos por ejemplo que los datos que obtenemos son:

```{r, echo = FALSE, include = FALSE}
set.seed(31521)
datos_clasif <- tibble(x = rep(1, 135), y = rep(1, 135)) %>%
  bind_rows(tibble(x = rep(0, 65), y = rep(1, 65))) %>%
  bind_rows(tibble(x = rep(1, 40), y = rep(0, 40))) %>%
  bind_rows(tibble(x = rep(1 ,100), y = rep(1, 100))) %>%
  sample_n(300) %>%
  mutate(caso = as.character(1: 300)) %>%
  select(caso, x, y)
```

```{r}
datos_clasif %>% head
```

Como explicamos arriba, nos interesa la diferencia. Calculamos $d$:

```{r}
datos_clasif <- datos_clasif %>%
  mutate(d = x - y)
datos_clasif %>% head
```
Y ahora calculamos la media de $d$ (y tasa de correctos de cada clasificador:)

```{r}
medias_tbl <-
  datos_clasif %>% summarise(across(where(is.numeric), mean, .names = "{col}_hat"))
d_hat <- pull(medias_tbl, d_hat)
medias_tbl
```
Ahora necesitamos calcular el error estándar. Como explicamos arriba, hacemos

```{r}
ee <- datos_clasif %>%
  mutate(d_hat = mean(d)) %>%
  mutate(dif_2 = (d - d_hat)) %>%
  summarise(ee = sd(dif_2) / sqrt(n())) %>%
  pull(ee)
ee
```

Y ahora podemos calcular la estadística $W$ y el valor-$p$ correspondiente:

```{r}
w <- d_hat / ee
valor_p <- 2 * (1 - pnorm(abs(w)))
c(w = w, valor_p = valor_p) %>% round(3)
```

Y vemos que tenemos evidencia considerable de que el desempeño no es el mismo:
el algoritmo 2 parece ser mejor.

```{block, type='ejercicio'}
¿Qué pasaría si incorrectamente usaras la prueba de dos muestras para este
ejemplo? ¿Qué cosa cambia en la fórmula de la estadística de Wald?
```

## Pruebas de cociente de verosimilitud {-}

Otra técnica clásica para hacer pruebas de hipótesis es el de cociente
de verosimilitudes. Con esta técnica podemos hacer pruebas que involucren
varios parámetros, y podemos contrastar hipótesis nulas contra alternativas
especificas.

Para aplicar este tipo de pruebas es necesario hacer supuestos distribucionales
(modelos probabilísticos), pues estas pruebas se basan en la función de
verosimilitud $\mathcal{L}(\theta; x_1,\ldots, x_n)$.

**Ejemplo**. Supongamos que tenemos la hipótesis nula de que una moneda
es justa ($p =0.05$ de sol). En 120 tiros de la moneda (que suponemos
independientes), observamos 75 soles.
Recordemos la función de log-verosimilitud para el modelo binomial
(ignorando constantes que no dependen de $p$) es
$$\ell(p) = 75 \log(p) + (120 - 75)\log(1-p) \,.$$

- Primero calculamos el estimador de máxima verosimilitud de $p$, que
es $\hat{p}_n = 75/120 = 0.625$. Evaluamos la verosimilitud

$$\ell(\hat{p}_n) = \ell(0.625) = 75\log(0.625) + 45\log(0.375) = -79.388\,.$$
- Ahora evaluamos la verosimlitud según la hipótesis nula, donde asumimos
que $p = 0.5$:

$$\ell(p_0) = \ell(0.5) = 75\log(0.5) + 45\log(0.5) = -83.177$$
- Finalmente, contrastamos estos dos números con una estadística que denotamos
con $\lambda_n$:

$$\lambda_n = 2\left[\ell(\hat{p}_n) - \ell(p_0)\right] = 2[\ell(0.625)- \ell(0.5)] = 2(3.79)=7.58\,.$$

- A $\lambda_n$ se le llama la **estadística de cociente de verosimilitud**.
Tomamos la diferencia de log verosimilitudes, que es los mismo que tomar el
logaritmo del **cociente** de verosimilitudes, y de ahí el nombre de la prueba.

- Nótese que cuando este número $\lambda_n$ es muy grande, esto implica que la hipótesis
nula es menos creíble, o menos consistente con los datos, pues la nula tiene mucho
menos verosimilitud de lo que los datos indican. Por otro lado, cuando este valor
es cercano a 0, entonces tenemos menos evidencia en contra de la
hipótesis nula. Esto se explica en la siguiente gráfica:

```{r}
log_verosim <- function(p){
  75 * log(p) + (120 - 75) * log(1 - p)
}
verosim_tbl <- tibble(p = seq(0.4, 0.7, 0.01)) %>%
  mutate(log_verosim = log_verosim(p))
ggplot(verosim_tbl, aes(x = p, y = log_verosim)) +
  geom_line() +
  geom_segment(x = 75/120, xend = 75/120, y = -130, yend = log_verosim(75/120), colour = "red") +
  geom_segment(x = 0.5, xend = 0.5, y = -130, yend = log_verosim(0.5), colour = "gray") +
  geom_errorbar(aes(x = 0.5,  ymin = log_verosim(0.5), ymax = log_verosim(75/120)),
               colour = "orange", width = 0.05) +
  annotate("text", x = 0.48, y = -81.5, label = "3.79") +
  annotate("text", x = 0.515, y = -91, label ="nula", colour = "gray20") +
  annotate("text", x = 0.665, y = -91, label ="max verosímil", colour = "red") +
  labs(subtitle = expression(paste(lambda, "=2(3.79)=7.58")))


```


Este método puede generalizarse para que no solo aplique a hipótesis
nulas donde $\theta = \theta_0$, sino en general, $\theta \in \Theta_0$.
Por ejemplo, podemos construir pruebas para $\theta < 0.4$.

```{block, type='mathblock'}
**Definición**. Consideramos la hipótesis nula $\theta= \theta_0$.
La **estadística del cociente de verosimilitudes** está dada por:

$$\lambda_n = 2\log\left( \frac{\max_{\theta \in \Theta}\, \mathcal{L}(\theta)}{\mathcal{L}(\theta_0)}        \right ) = 2\log\left(  \frac{\mathcal{L}(\hat{\theta}_n)}{\mathcal{L}(\theta_0)}  \right)\,.$$

donde $\hat{\theta}_n$ es el estimador de máxima verosimilitud.
```

Para construir una prueba asociada, como siempre, necesitamos una distribución de
referencia. Esto podemos hacerlo con simulación, o usando resultados asintóticos.


### Distribución de referencia para pruebas de cocientes {-}

Para nuestro ejemplo anterior, podemos simular datos bajo la hipótesis nula,
y ver cómo se distribuye la estadística $\lambda_n$:

**Ejemplo.** Simulamos bajo la hipótesis nula como sigue:

```{r}
n_volados <- 120 # número de volados
simulados_nula <- rbinom(4000, n_volados, p = 0.5)
lambda <- function(n, x, p_0 = 0.5){
  # estimador de max verosim
  p_mv <- x / n
  # log verosimilitud bajo mv
  log_p_mv <- x * log(p_mv) + (n - x) * log(1 - p_mv)
  # log verosimllitud bajo nula
  log_p_nula <- x * log(p_0) + (n - x) * log(1 - p_0)
  lambda <- 2*(log_p_mv - log_p_nula)
  lambda
}
lambda_obs <- lambda(n_volados, 75, 0.5)
sims_tbl <- tibble(sim_x = simulados_nula) %>%
  mutate(lambda = map_dbl(sim_x, ~ lambda(n_volados, .x, p_0 = 0.5)))
ggplot(sims_tbl, aes(x = lambda)) +
  geom_histogram(binwidth = 0.7) +
  geom_vline(xintercept = 2.92, colour = "red")
```
Con esta aproximación a la distribución de referencia podemos calcular
el valor-$p$ en nuestro ejemplo anterior:

```{r}
valor_p <- mean(sims_tbl$lambda >= lambda_obs)
valor_p
```
y observamos que tenemos evidencia fuerte en contra de la hipótesis nula:
la moneda no está balanceada.

**Ejemplo**. Este ejemplo es un poco artificial, pero lo usamos
para entender mejor las pruebas de cocientes
de verosimlitud. Supongamos que tenemos una muestra de $\mathsf{N}(\mu, 1)$, y queremos
probar si $\mu = 8$. Asumimos que el supuesto de normalidad y desviación
estándar iugal a 1 se cumplen.

```{r}
set.seed(3341)
n_muestra <- 100
muestra_1 <- rnorm(n_muestra, 7.9, 1)
```



```{r}
crear_log_p <- function(x){
  # crear log verosim para dos muestras normales independientes.
  log_p <- function(params){
    mu <- params[1]
    log_vero <- dnorm(x, mean = mu, sd = 1, log = TRUE) %>% sum
    log_vero
  }
}
lambda_calc <- function(muestra, crear_log_p){
  log_p <- crear_log_p(muestra)
  res <- optim(c(0), log_p, control = list(fnscale = -1))
  lambda_mv <- log_p(res$par)
  lambda_nula <- log_p(8.0)
  lambda <- 2 * (lambda_mv - lambda_nula)
  lambda
}
lambda <- lambda_calc(muestra_1, crear_log_p)
lambda
```

Ahora construimos con simulación la distribución de referencia usando simulaciones
bajo la nula

```{r, cache = TRUE}
sims_nula <- map(1:10000, ~ rnorm(n_muestra, 8, 1))
lambda_nula_sim <- map_dbl(sims_nula, ~ lambda_calc(.x, crear_log_p))
tibble(lambda = lambda_nula_sim) %>%
  ggplot(aes(x = lambda)) + geom_histogram() +
  geom_vline(xintercept = lambda, colour = "red")
```
```{r}
valor_p <- mean(lambda_nula_sim >= lambda)
valor_p
```

Estos datos muestran consistencia con la hipótesis $\mu = 8$.

**Discusión**: Nota en los dos ejemplos anteriores la similitud entre
las distribuciones de referencia. En ambos casos, estas distribuciones
resultan ser aproximadamente **$\chi$-cuadrada con 1 grado de libertad** (ji-cuadrada). Podemos
checar para el último ejemplo:

```{r}
teorica <- tibble(x = seq(0.1, 10, 0.01)) %>%
  mutate(f_chi_1 = dchisq(x, df = 1))
tibble(lambda = lambda_nula_sim) %>%
  ggplot() + geom_histogram(aes(x = lambda, y = ..density..), binwidth = 0.1) +
  geom_line(data = teorica, aes(x = x, y = f_chi_1), colour = "red")
```

O mejor, con una gráfica de cuantiles de las simulaciones vs
la téorica:

```{r}
tibble(lambda = lambda_nula_sim) %>%
  ggplot(aes(sample = lambda)) +
  geom_qq(distribution = stats::qchisq, dparams = list(df = 1)) +
  geom_qq_line(distribution = stats::qchisq, dparams = list(df = 1))
```

Este resultado asintótico no es trivial, y se usa comúnmente para calcular
valores $p$. Discutiremos más este punto más adelante.


## Otro tipo de pruebas {-}

Con cocientes de verosimlitud podemos diseñar pruebas para contrastar
condiciones que sólo un subconjunto de parámetros cumple.


**Ejemplo.** Supongamos que queremos hacer una prueba de
igualdad de medias $\mu_1 = \mu_2$ para dos poblaciones normales
$\mathsf{N}(\mu_1, \sigma_1)$ y $\mathsf{N}(\mu_2, \sigma_2)$, donde extraemos las
muestras de manera independiente, y no conocemos
las desviaciones estándar. Obtenemos dos muestras (que supondremos
provienen de distribuciones normales, pues ese es nuestro supuesto)

```{r}
set.seed(223)
muestra_1 <- rnorm(80, 0.8, 0.2)
muestra_2 <- rnorm(120, 0.8, 0.4)
```

Necesitamos: 1) calcular el valor de la estadística $\lambda_n$ de cociente de
verosimilitudes, 2) Calcular la distribución de referencia para $\lambda_n$ bajo
la hipótesis nula y finalmente 3) Ver qué tan extremo es el valor obtenido de
$\lambda_n$ en relación a la distribución de referencia.


```{r}
crear_log_p <- function(x, y){
  # crear log verosim para dos muestras normales independientes.
  log_p <- function(params){
    mu_1 <- params[1]
    mu_2 <- params[2]
    sigma_1 <- params[3]
    sigma_2 <- params[4]
    log_vero_1 <- dnorm(x, mean = mu_1, sd = sigma_1, log = TRUE) %>% sum
    log_vero_2 <- dnorm(y, mean = mu_2, sd = sigma_2, log = TRUE) %>% sum
    log_vero <- log_vero_1 + log_vero_2 #se suman por independiencia
    log_vero
  }
}
log_p <- crear_log_p(muestra_1, muestra_2)
```

```{r}
crear_log_p_nula <- function(x, y){
  log_p <- function(params){
    # misma media
    mu <- params[1]
    sigma_1 <- params[2]
    sigma_2 <- params[3]
    log_vero_1 <- dnorm(x, mean = mu, sd = sigma_1, log = TRUE) %>% sum
    log_vero_2 <- dnorm(y, mean = mu, sd = sigma_2, log = TRUE) %>% sum
    log_vero <- log_vero_1 + log_vero_2 #se suman por independiencia
    log_vero
  }
}
log_p_nula <- crear_log_p_nula(muestra_1, muestra_2)
```

Ahora tenemos el problema de que no conocemos las sigma. Estas deben ser
estimadas para después calcular el cociente de verosimilitud:

```{r}
res <- optim(c(0, 0, 1, 1), log_p, method = "Nelder-Mead",
             control = list(fnscale = -1))
res$convergence
est_mv <- res$par
names(est_mv) <- c("mu_1", "mu_2", "sigma_1", "sigma_2")
est_mv
```
Y tenemos
```{r}
lambda_1 <- log_p(est_mv)
lambda_1
```
Ahora calculamos el máximo bajo el supuesto de la hipótesis nula:

```{r}
res <- optim(c(0, 1, 1), log_p_nula, method = "Nelder-Mead",
             control = list(fnscale = -1))
res$convergence
est_mv_nula <- res$par
names(est_mv) <- c("mu", "sigma_1", "sigma_2")
est_mv_nula
```

y evaluamos

```{r}
lambda_2 <- log_p_nula(est_mv_nula)
lambda_2
```

Finalmente, nuestra estadística $\lambda_n$ es
```{r}
lambda <- 2 * (lambda_1 - lambda_2)
lambda
```

Y ahora necesitamos calcular un valor-$p$. El problema que tenemos en este
punto es que bajo la hipótesis nula no están determinados todos los parámetros,
así que no podemos simular de manera simple muestras para obtener la
distribución de referencia. Podemos sin embargo usar bootstrap paramétrico
usando los estimadores de máxima verosimilitud bajo la nula

```{r, cache = TRUE}
simular_boot <- function(n_1, n_2, est_mv_nula){
  x <- rnorm(n_1, est_mv_nula[1], est_mv_nula[2])
  y <- rnorm(n_2, est_mv_nula[1], est_mv_nula[3])
  list(x = x, y = y)
}
lambda_nula_sim <- function(est_mv_nula){
  muestras <- simular_boot(80, 120, est_mv_nula)
  x <- muestras$x
  y <- muestras$y
  log_p <- crear_log_p(x, y)
  log_p_nula <- crear_log_p_nula(x, y)
  est_1 <- optim(c(0,0,1,1), log_p, control = list(fnscale = -1))
  est_2 <- optim(c(0,1,1), log_p_nula, control = list(fnscale = -1))
  lambda <- 2*(log_p(est_1$par) - log_p_nula(est_2$par))
  lambda
}
lambda_sim <- map_dbl(1:2000, ~ lambda_nula_sim(est_mv_nula = est_mv_nula))
```


Y graficamos la distribución de referencia junto con el valor de $\lambda_n$
que obtuvimos:

```{r}
tibble(lambda = lambda_sim) %>%
  ggplot(aes(x = lambda)) +
  geom_histogram() +
  geom_vline(xintercept = lambda, colour = "red")
```
Y claramente los datos son consistentes con medias iguales. El valor-$p$ es
```{r}
mean(lambda_sim > lambda)
```


Verificamos
una vez más que la distribución de referencia es cercana a una $\chi$-cuadrada
con un grado de libertad.

```{r}
tibble(lambda = lambda_sim) %>%
  ggplot(aes(sample = lambda)) +
  geom_qq(distribution = stats::qchisq, dparams = list(df = 1)) +
  geom_qq_line(distribution = stats::qchisq, dparams = list(df = 1))

```

Esta es la definición generalizada de las pruebas de cociente de verosimilitudes

```{block, type='mathblock'}
**Definición**. Consideramos la hipótesis nula $\theta \in \Theta_0$.
La **estadística del cociente de verosimilitudes** está dada por:

$$\lambda_n = 2\log\left( \frac{{\max}_{\theta \in \Theta} \,\, \mathcal{L}(\theta)}{\max_{\theta \in \Theta_0} \, \mathcal{L}(\theta)}        \right ) = 2\log\left(  \frac{ \mathcal{L}(\hat{\theta}_n)}{\mathcal{L}(\hat{\theta}_0)}  \right)\,,$$

donde $\hat{\theta}_n$ es el estimador de máxima verosimilitud de $\theta$ y
$\hat{\theta}_0$ es el estimador de máxima verosimilitud de $\theta$ cuando
restringimos a que $\theta \in \Theta_0$.
```

En nuestro ejemplo anterior, el espacio $\Theta_0$ era
$\{ (\mu,\mu,\sigma_1, \sigma_2)\}$, que es un subconjunto de
$\{ (\mu_1,\mu_2,\sigma_1, \sigma_2)\}$. Nótese que el espacio $\Theta_0$
tiene tres parámetros libres, mientras que el espacio total tiene 4.

Aunque podemos usar el bootstrap paramétrico para construir distribuciones
de referencia para estas pruebas y calcular un valor-$p$, el siguiente
teorema, cuya demostración no es trivial, explica las observaciones que
hicimos arriba. Este teorema enuncia la estrategia del enfoque clásico,
que utiliza una aproximación asintótica.

```{block2, type='mathblock'}
**Valores-p para pruebas de cocientes de verosimilitud**. Supongamos
que $\theta = (\theta_1,\theta_2, \ldots, \theta_p)$.  Sea
$$\Theta_0 = \{\theta : \theta_1 = a_1, \theta_2 = a_2, \dots, \theta_q = a_q \}\,,$$
  es decir la hipótesis $\theta \in \Theta_0$ es que los primeros $q$ parámetros
de $\theta$ estan fijos en algún valor. Los otros parámetros no se consideran en esta
prueba.

Si $\lambda_n$ es la estadística de cociente de verosimilitudes de esta prueba, entonces,
bajo la nula $\theta \in \Theta_0$ tenemos que la distribución
de $\lambda_n$ es asintóticamente $\chi$-cuadrada con $q$ grados de libertad, denotada por $\chi^2_q$.

El valor-$p$ para esta prueba es
$$P(\chi^2_{q} > \lambda)\,.$$
```

**Observaciones**:

- Para hacer cálculos con la distribución $\chi$-cuadrada usamos rutinas numéricas
(por ejemplo la función `pchisq` en R).

- Nótese que $p$ es la dimensión del espacio $\Theta$ ($p$ parámetros),
y que $p-q$ es la dimensión del espacio $\Theta_0$ (pues $q$ parámetros están fijos),
de  modo que los grados de libertad son la dimensión de $\Theta$ menos la
dimensión de $\Theta_0$.

- En nuestro primer ejemplo (proporción de éxitos) solo teníamos un parámetro. El espacio
$\Theta_0$ es de dimensión 0, así que los grados de libertad son $1 = 1 - 0$
- En este último ejemplo donde probamos igualdad de medias, el espacio $\Theta$ tiene
dimensión 4, y el espacio $\Theta_0$ es de dimensión 3 (tres parámetros libres), por
lo tanto los grados de libertad son $1 = 4 -3$.


**Ejemplo** En nuestro ejemplo de prueba de igualdad de medias,
usaríamos
```{r}
pchisq(lambda, df =1, lower.tail = FALSE)
```
que es similar al que obtuvimos con la estrategia del bootstrap paramétrico.



## Errores tipo I y tipo II {-}

En algunas ocasiones, en lugar de solamente calcular un valor-$p$ queremos
tomar una decisión asociada a distintas hipótesis que consideramos posibles. Por
ejemplo, nuestra hipótesis nula podría ser

- Hipótesis nula $H_0$: Una medicina nueva que estamos probando no es efectiva
en reducir el colesterol en pacientes.

Y queremos contrastar con una alternativa:

- Hipótesis alternativa $H_A$: la medicina nueva reduce los niveles de colesterol
en los pacientes.

La decisión que está detrás de estas pruebas es: si no podemos rechazar la nula,
la medicina no sale al mercado. Si rechazamos la nula, entonces la medicina es
aprobada para salir al mercado.

Para diseñar esta prueba, procedemos como sigue:

1. Definimos cómo recolectar datos $X$ de interés. 
2. Definimos una estadística $T(X)$ de los datos.
3. Definimos una **región de rechazo** $R_0$ de valores tales que si $T(X)\in
R_0$, entonces rechazaremos la hipótesis nula (e implícitamente tomaríamos la
*decisión* asociada a la alternativa).

Ejecutamos la prueba observando datos $X=x$, calculando $T(x)$, y checando si
$T(x) \in R_0$. Si esto sucede entonces decimos que rechazamos la hipótesis
nula, y tomamos la *decisión* asociada a la alternativa.

**Ejemplo**. Si tenemos la hipótesis nula $p_1=0.5$ para una proporción, y al
alternativa es $p_1\neq 0.5$, podemos usar la estadística de Wald $T(x) =
\frac{\hat{p_1} - 0.5}{\hat{\textsf{ee}}_n}$. Podríamos definir la región de
rechazo como $R_0 = \{T(x) : |T(x)| > 3 \}$ (rechazamos si en valor absoluto la
estadística de Wald es mayor que 3).

Cuando diseñamos una prueba de este tipo, quisiéramos minimizar dos tipos de errores:

- Rechazar la hipótesis nula $H_0$ cuando es cierta: **Error tipo I**
- No rechazar la hipótesis nula $H_0$ cuando $H_0$ es falsa: **Error tipo II**

La gravedad de cada error depende del problema. En nuestro ejemplo de la medicina, por
ejemplo:

- Un error tipo II resultaría en una medicina efectiva que no sale al mercado, lo
que tiene consecuencias financieras (para la farmaceútica) y costos de oportunidad
en salud (para la población). Por otra parte,
- Un error tipo I resultaría en salir al mercado con una medicina que no es efectiva. Esto tiene
costos de oportunidad financieros que pueden ser grandes para la sociedad.
Todos estos costos dependen, por
ejempĺo, de qué tan grave es la enfermedad, qué tan costosa es la medicina, y así sucesivamente.
- En el enfoque más clásico, los errores tipo I y tipo II generalmente no se balancean
según su severidad o probabilidad. En lugar de eso, generalmente se establece
un límite para la probabilidad de cometer un error del tipo I (usualmente 5%, por una
tradición que no tiene mucho fundamento)

En vista de este ejemplo simple, y las observaciones de arriba:

- **Reducir una decisión compleja** a una prueba de hipótesis con resultados
binarios (rechazar o no) es generalmente **erróneo**.
- **Las pruebas de hipótesis se usan muchas veces incorrectamente cuando lo más
apropiado es usar estimación** por intervalos o algo similar que cuantifique la
incertidumbre de las estimaciones.

Consulta por ejemplo [el comunicado de la ASA acerca de p-values y pruebas de hipótesis](https://amstat.tandfonline.com/doi/full/10.1080/00031305.2016.1154108)

En el caso de la medicina, por ejemplo, realmente no nos interesa que
la medicina sea mejor que un placebo. Nos importa que tenga un efecto considerable en los pacientes.
Si estimamos este efecto, incluyendo incertidumbre, tenemos una mejor herramienta
para hacer análisis costo-beneficio y tomar la decisión más apropiada.

Como dijimos, típicamente se selecciona la región de rechazo de forma
que bajo la hipótesis nula la probabilidad de cometer un error tipo I está acotada.

```{block, type='mathblock'}
**Definición**. Supongamos que los datos $X_1,X_2,\ldots, X_n$ provienen de una
distribución $F_\theta$, donde no conocemos $\theta$. Supongamos que la
hipótesis nula es que $\theta = \theta_0$ (que llamamos una hipótesis simple).

La **función de potencia** de una prueba con región de rechazo $R_0$ se define
como la probabilidad de rechazar para cada posible valor del parámetro $\theta$
$$\beta(\theta) = \mathbb{P}_\theta (T(X)\in R_0)\,.$$

El **tamaño** de una prueba se define como el valor
$$\alpha := \beta(\theta_0)\,,$$
es decir, la probabilidad de rechazar la nula ($\theta = \theta_0$)
erróneamente.
```

**Observación**. Esto se generaliza para hipótesis compuestas, donde la nula
es que el parámetro $\theta$ está en un cierto conjunto $\Theta_0$. Por ejemplo,
una hipótesis nula puede ser $\theta < 0.5$. En este caso, $\alpha$ se define
como el valor más grande que $\beta(\theta)$ toma cuando $\theta$ está en $\Theta_0$,
es decir, la probabilidad de rechazo más grande cuando la hipótesis nula se cumple.

Decimos que una prueba tiene **nivel de significancia** de $\alpha$ si su tamaño
es menor o igual a $\alpha$.

Decimos que la **potencia de una prueba** es la probabilidad de, correctamente,
rechazar la hipótesis nula cuando la alterna es verdadera:
$$1-\beta := \beta(\theta_a) = \mathbb{P}_{\theta_a} (T(X) \in R_0)\,.$$
**Observación:** Sería deseable encontrar la prueba con mayor potencia bajo
$H_a$, entre todas las pruebas con tamaño $\alpha$. Esto no es trivial y no
siempre existe.

**Observación:** El valor $p$ es el menor tamaño con el que podemos rechazar $H_0$.

**Ejemplo** (@Chihara) Supongamos que las calificaciones de Enlace de alumnos en México
se distribuye aproximadamente como una normal con media 515 y desviación estándar
de 120. En una ciudad particular, se quiere decidir
si es neceario pedir fondos porque la media de la ciudad es más baja
que la nacional. Nuestra hipótesis nula es $H_0: \mu \geq 515$ y la
alternativa es $\mu < 515$, así que si rechazamos la nula se pedirían los fondos.

Supondremos que la distribución de calificaciones
en la ciudad es también aproximadamente normal con desviación estándar de 130. Se
plantea tomar una muestra de 100 alumnos, y rechazar si la media muestral $\bar{X}$ es
menor que 505. ¿Cuál es la probabilidad $\alpha$ de tener un error de tipo I?


La función de potencia es
$$\beta(\mu) = \mathbb{P}_\mu(\bar{X} < 505)\,.$$
Restando la media $\mu$ y estandarizando obtenemos
$$\beta(\mu) = \mathbb{P} \left (\frac{\bar{X} - \mu}{130/\sqrt{100}} < \frac{505 -\mu}{130/\sqrt{100}} \right )\,,$$
así que
$$\beta(\mu) = \Phi \left (\frac{505 -\mu}{130/\sqrt{100}}\right )\,,$$
donde $\Phi$ es la función acumulada de la normal estándar. La gráfica
de la función potencia es entonces

```{r}
potencia_tbl <- tibble(mu = seq(450, 550, 1)) %>%
  mutate(beta = pnorm((505 - mu)/13)) %>% # probabilidad de rechazar
  mutate(nula_verdadera = factor(mu >= 515)) # nula verdadera
ggplot(potencia_tbl, aes(x = mu, y = beta, colour = nula_verdadera)) +
  geom_line()
```

Es decir, si la media $\mu$ de la ciudad es muy baja, con mucha seguridad rechazamos. Si
es relativamente alta entonces no rechazamos. El tamaño de la prueba es el mayor
valor de probabilidad de rechazo que se obtiene sobre los valores $\mu\geq 515$ (la nula).
Podemos calcularlo analíticamente como sigue:

Si $\mu \geq 515$, entonces
$$\beta(\mu) \leq \beta(515) = \Phi\left (\frac{505 -515}{130/\sqrt{100}}\right ) = \Phi( - 10 / 13) = \Phi(-0.7692)\,,$$
que es igual a
```{r}
pnorm(-0.7692)
```
Y este es el tamaño de la prueba. En otras palabras: si la ciudad no está
por debajo de la media nacional, hay una probabilidad de 22% de que erróneamente
se pidan fondos (al rechazar $H_0$).

**Ejemplo** Supongamos que los que programan el presupuesto deciden que se requiere
tener una probabilidad de a lo más 5% de rechazar erróneamente la hipótesis nula (es decir,
pedir fondos cuando en realidad su media no está debajo de la nacional) para
poder recibir fondos. ¿Cuál es la región de rechazo que podríamos escoger?

En el caso anterior usamos la región $\bar{X}<505$. Si el tamaño de muestra
está fijo en $n=100$ (por presupuesto), entonces tenemos que escoger un punto
de corte más extremo. Si la región de rechazo es $\bar{X} < C)$ entonces
tenemos, siguiendo los cálculos anteriores, que
$$0.05 = \alpha = \Phi \left ( \frac{C -515}{130/\sqrt{100}}\right) = \Phi \left( \frac{C- 515}{13}   \right)\,.$$
Buscamos el cuantil 0.05 de la normal estándar, que es

```{r}
z_alpha <- qnorm(0.05)
z_alpha
```
Y entonces requerimos que

$$\frac{C- 515}{13} = -1.6448\,.$$
Despejando obtenemos

```{r}
C <- 13*z_alpha + 515
C
```

Así que podemos usar la región $\bar{X} < 493.5$, que es más *estricta*
que la anterior de $\bar{X} < 505$.

```{block, type='ejercicio'}
Considera la potencia de la prueba $\beta(\mu)$ que vimos arriba. Discute y
corre algunos ejemplos para contestar las siguientes preguntas:

- Recuerda la definición: ¿qué significa $\beta(\mu)$?
- ¿Qué pasa con la potencia cuando $\mu$ está más lejos de los valores de la hipótesis nula?
- ¿Qué pasa con la potencia cuando hay menos variabilidad en la población? ¿Y cuando
la muestra es más grande?
- ¿Qué pasa si hacemos más chico el nivel de significancia?

```


## Consideraciones prácticas {-}

Algunos recordatorios de lo que hemos visto:

- Rechazar la nula **no quiere decir que la nula es falsa**, ni que **encotramos un
"efecto"**. Un valor-$p$ chico tampoco
quiere decir que la nula es falsa. Lo que quiere decir es que la nula es poco consistente
con los datos que observamos, o que es muy poco probable que la nula produzca los datos
que observamos.

- Rechazar la nula (encontrar un efecto "significativo") **no quiere decir que
el efecto tiene importancia práctica**. Si la potencia es alta (por ejemplo
cuando el tamaño de muestra es grande), puede ser que la discrepancia de los
datos con la nula es despreciable, entonces para fines prácticos podríamos
trabajar bajo el supuesto de la nula. Por eso en general preferimos hacer
estimación que pruebas de hipótesis para entender o resumir los datos y tamaños
de las discrepancias.

- Adicionalmente, muchas de las hipótesis nulas que generalmente se utilizan se
pueden rechazar sin datos (por ejemplo, igualdad de proporciones en dos poblaciones
reales). Lo que importa es qué tan diferentes son, y qué tan bien podemos estimar
sus diferencias.

- En la literatura, muchas veces parece que "encontrar una cosa interesante" es
rechazar una hipótesis nulas con nivel 5% de significancia. Es **más importante
entender cómo se diseñó el estudio, cómo se recogieron los datos, cuáles fueron
las decisiones de análisis** que pasar el mítico nivel de 5%.

- Cuando la potencia es baja (por ejemplo porque el tamaño de muestra es muy
chico), tenemos que observar diferencias muy grandes para rechazar. Si probamos
algo **poco factible** (por ejemplo, que la vitamina $C$ aumenta la estatura a
los adultos), entonces **los rechazos generalmente se deben a variabilidad en la
muestra** (error tipo II).

- Cuando diseñamos y presentamos resultados de un estudio o análisis, es mejor
pensar en describir los datos y su variabilidad, y mostrar estimaciones
incluyendo fuentes de incertidumbre, en lugar de intentar resumir con un
valor-$p$ o con el resultado de una prueba de hipótesis.

## Pruebas múltiples {-}

En algunas ocasiones se hacen muchas pruebas para "filtrar" las cosas que son
interesantes y las que no. Por ejemplo, cuando comparamos miles de genes entre
dos muestras (la nula es que son similares). Si cada prueba se conduce a un
nivel $\alpha$, la probablilidad de tener al menos un rechazo falso (un error
tipo I) es considerablemente más alta que $\alpha$.

Por ejemplo, si repetimos una prueba de hipótesis con nivel $\alpha$ con
muestras independientes, la probabilidad de tener al menos un rechazo falso es
$1-(1-\alpha)^n$, que es muy cercano a uno si $n$ es grande (¿cómo derivas esta
fórmula?). Por ejemplo, si $\alpha = 0.05$ y $n = 100$, con más de 99%
probabilidad tendremos al menos un rechazo falso, o un "hallazgo" falso. Sin $n$
es muy grande, varios de los hallazgos que encontremos serán debidos a
variabilidad muestral.

Puedes ver en [@Wasserman], sección 10.7 métodos conservadores como corrección
de Bonferroni (sólo rechazar cuando el valor-$p$ es menor a $0.05/n$), o la
técnica más moderna de control de tasa de descubrimientos falsos (FDR).

Cuando estamos en una situación como esta (que es más retadora en cuanto a
análisis), sin embargo, sugerimos usar estimaciones que tomen cuenta todos los
datos con regularización apropiada: por ejemplo, en lugar de trabajar con cada
muestra por separado, intentamos construir un modelo para el proceso completo de
muestreo. Una posibilidad son modelos jerárquicos bayesianos. Ver por ejemplo
[@gelmansignif].

## Meta-análisis {-}

El meta-análisis es un estudio estadístico de los resultados de múltiples
estudios científicos. Usualmente se utilizan cuando diversos estudios
intentan responder a la misma pregunta, y cada estudio incluye un medida de los 
errores de estimación. Por ejemplo, el siguiente caso viene de @anxiety. 

Los autores han recopilado distintos estudios sobre la prevalencia ---la
proporción de personas en un estudio que sufren de una enfermedad--- de ansiedad
en los jóvenes durante la pandemia del COVID-19.

En total se tienen 26 estudios distintos. Donde cada uno reporta el número de 
jóvenes con diagnóstico positivo de ansiedad dentro de sus grupos de estudio. 
Por ejemplo, podemos ver una muestra de los datos. 

```{r, echo = FALSE}

datos <- tibble(
       publicacion = factor(seq(26, 1)), 
       casos   = c(278, 213, 252, 864, 3334, 67, 1899, 5136, 4240, 265, 40, 309, 
                   1941, 175, 403, 805, 857, 98, 362, 179, 1402, 248, 155, 122, 
                   103, 30), 
       total   = c(436, 384, 455, 1794, 7772, 168, 4805, 13002, 11180, 859, 156, 
                   1241, 7890, 721, 1784, 3613, 4342, 515, 2050, 1018, 8079, 
                   1550, 1025, 1036, 2078, 1352))

datos <- datos %>% mutate(
        estimate = casos/total, 
        se       = sqrt(estimate * (1 - estimate)/total), 
        lo.bound = estimate - 2 * se, 
        hi.bound = estimate + 2 * se
)

datos |> head()

```

```{block, type="ejercicio"}

¿Qué modelo paramétrico crees que fue utilizado para generar estos resúmenes?

```

Esta información la podríamos graficar como sigue (cada estudio se muestra como
una observación)

```{r}

datos %>% 
    ggplot(aes(x = estimate, y = publicacion)) +
    geom_linerange(aes(xmin = lo.bound, xmax = hi.bound)) + 
    geom_point() + 
  xlab("Estimador")+ ylab("Publicación (estudio)") + 
  ggtitle("Meta-análisis de prevalencia de ansiedad en jóvenes")

```
De lo cual puede resultar un poco confuso que todos los estudios tienen como
objetivo estimar la misma cantidad. Pero parece ser que hay resultados muy
diferentes entre todos estos estudios. Más aún, la precisión con la que 
se estima la proporción de jóvenes con ansiedad parece ser muy distinta entre 
estudios. 

Podríamos ajustar de manera ingenua un estimador utilizando como se muestra 
a continuación: 

```{r mean estimates, echo = FALSE}

datos %>% 
    summarise( casos = sum(casos), 
               total = sum(total)) %>% 
    mutate(
        estimate = casos/total, 
        se       = sqrt(estimate * (1 - estimate)/total), 
        lo.bound = estimate - 2 * se, 
        hi.bound = estimate + 2 * se, 
        width    = hi.bound - lo.bound
        )

```

```{block, type = "ejercicio"}

¿Consideras que es un buen estimador de acuerdo a lo observado en los estudios 
individuales?

```

Podríamos suponer que aunque las medias sean diferentes los errores estándar en
las estimaciones deberían de ser similares. Para esto podemos utilizar una
estimación de los errores estándar por grupo ([*pooled variance*](https://en.wikipedia.org/wiki/Pooled_variance)):

$$\hat \sigma^2 = \frac{(n_1 - 1) \hat s_1^2 + \cdots + (n_k - 1) \hat s_k^2}{n_1 + \cdots + n_k - k}\,,$$
donde $k$ denota el número de estudios, $\hat s_i$ la estimación del error 
estándar y $n_k$ el tamaño de muestra de cada estudio. 

Es claro, por la gráfica de arriba, que este supuesto en la varianza no es el
correcto. Pero, nos lleva a un estimador con un error estándar mas amplio, pues
nos ayuda a incoporar heterogeneidad en las medias de los grupos. Esto a su vez,
nos puede "proteger" al tener una mejor cuantificación del error en la
estimación.

```{r pooled estimates}

datos %>% 
    summarise(average  = mean(estimate), 
              se       = sqrt(sum((total - 1) * (se**2))/(sum(total) - n())), 
              lo.bound = average - 2 * se, 
              hi.bound = average + 2 * se, 
              width    = hi.bound - lo.bound)
    
```

```{block, type = "ejercicio"}

Utiliza *bootstrap* no paramétrico para estimar la distribución de muestreo de la
prevalencia de ansiedad en jóvenes durante la pandemia. ¿Cómo se compara 
tu estimador con los propuestos arriba? Al hacer el *bootstrap* ¿qué supuestos 
del método generador de datos estás incorporando?

```


Este último ejemplo es particular en algunos aspectos. El más importante es que
nos muestra que un estimador con menor incertidumbre no necesariamente es mejor
que uno con mayor. El proceso de estimación puede incorporar supuestos que hacen
transparente el proceso generador de datos, lo cual puede llevar a estimaciones
mas prudentes o conservadoras. Esto último, lo leemos como el *slogan* de
nuestro curso: "no es malo tener incertidumbre, el problema es **no**
cuantificarla".

```{block, type = "ejercicio"}

¿Qué características de nuestros estimadores no están siendo capturados en 
los ejercicios de arriba con respecto a la prevalencia de ansiedad en jóvenes 
durante la pandemia? ¿Encuentras una manera de subsanar esas deficiencias?

```


En la sección siguiente veremos una paradigma distinto para hacer inferencia
estadística. Posiblemente se te ocurrirá cómo atacar el problema del
meta-análisis con las técnicas a continuación.