Simulación de modelos probabilísticos
Teresa Ortiz
Distribuciones multivariadas
Hasta ahora hemos estudiado distribuciones univariadas, sin embargo, es común que un modelo probabilístico involucre más de una variable aleatoria por lo que es necesario el concepto de distribuciones de probabilidad multivariadas.
La distribución conjunta sobre un conjunto de variables aleatorias \(\{X_1,...,X_n\}\), que denotamos \(p(x_1,...,x_n)\), asigna probabilidades a todos los eventos determinados por el conjunto de variables aleatorias.
En el caso discreto bivariado, dado las variables aleatorias discretas \(X\) y \(Y\), definimos la función de densidad conjunta como \(f(x,y)=P(X=x, Y=y)\).
Ejemplo. Consideremos una distribución sobre la población de departamentos en renta de Hong Kong, el espacio de resultados es el conjunto de todos los departamentos en la población. En muchas ocasiones buscamos resolver preguntas que involucran más de una variable aleatoria, en este ejemplo nos interesan:
Renta mensual: toma los valores baja (≤1k), media ((1k,5k]), media alta ((5k,12k]) y alta (>12k).
Tipo de departamento: toma 3 valores, púlico, privado u otros.
La distribución conjunta de variables aleatorias discretas se puede representar por medio de tablas.
| Renta/Tipo | público | privado | otros |
|---|---|---|---|
| baja | 0.17 | 0.01 | 0.02 |
| media | 0.44 | 0.03 | 0.01 |
| media alta | 0.09 | 0.07 | 0.01 |
| alta | 0 | 0.14 | 0.10 |
En el caso continuo bivariado, decimos que la función \(p(x,y)\) es una función de densidad de probabilidad para las variables aleatorias \((X,Y)\) si:
\(p(x,y) \geq 0\) para toda \((x,y)\).
\(\int_{-\infty}^{\infty}p(x,y)dxdy=1\).
- Para cualquier conjunto \(A \subset \mathbb{R} \times \mathbb{R}\), \(P((X,Y) \in A) = \int\int_A p(x,y)dxdy\).
Ejemplo. Sean \((X,Y)\) uniformes en el cuadrado unitario, entonces \[ p(x,y) = \left\{ \begin{array}{lr} 1 & 0\leq x \leq 1,0\leq y \leq 1\\ 0 & e.o.c. \end{array} \right. \]
Para encontrar \(P(X < 1/2, Y<1/2)\), esto es la probailidad del evento \(A=\{X<1/2, Y<1/2\}\). La integral de \(p\) sobre este subconjunto corresponde, en este caso, a calcular el área del conjunto \(A\) que es igual a 1/4.
De la distribución conjunta \(p(x_1,...,x_n)\) podemos obtener la distribución de únciamente una variable aleatoria \(X_j\), donde \(X_j \in \{X_1,...,X_n\}\), la llamamos la distribución marginal de \(X_j\).
Sea \(\{X_1,...,X_n\}\) un conjunto de variables aleatorias con distribución conjunta \(p(x_1,...,x_n)\), la distribución marginal de \(X_j\) (\(j \in \{1,...,n\}\)) se define como, \[p_{X_j}(x_j) = \sum_{x_1,...,x_{j-1},x_{j+1},...,x_n}p(x_1,...,x_n)\mbox{ en el caso discreto,}\] \[p_{X_j}(x_j) = \int_{x_1,...,x_{j-1},x_{j+1},...,x_n}p(x_1,...,x_n)dx_1,...,dx_n\mbox{ en el caso continuo}\]
Ejemplo. Retomando el problema de los departamentos, ¿Cuál es la probabilidad de que un departamento elegido al azar sea público?
Probabilidad condicional
Sean \(A\), \(B\) dos eventos entonces, con \(P(B)>0\), entonces la probabilidad condicional de \(A\) dado \(B\) es
\[P(A|B)=\frac{P(AB)}{P(B)}\]
Ejemplo. ¿Cuál es la probabilidad de que un departamento privado tenga renta baja?
¿Cómo se compara con la probabilidad de que la renta sea baja (desconozco el tipo de departamento)?
La noción de probabilidad condicional se extiende a distribuciones condicionales:
Sean \(X\), \(Y\) dos variables aleatorias con función de densidad conjunta \(p(x,y)\), entonces la función de densidad condicional de \(X\) dado \(Y=y\), para toda \(y\) tal que \(p_Y(y) > 0\), se define como \[p_{X\vert Y}(x\vert y) = \frac{p(x, y)}{p_Y(y).}\]
Ejemplo. ¿Cuál es la distribución condicional de renta dado renta privado ?
Para obtener toda la distribución condicional calculamos los dos casos restantes (renta media, media alta y alta).
Vale la pena destacar que una distribución condicional es una distribución de probabilidad. En el ejemplo anterior, notemos que cada renglón de la tabla probabilidades suman uno, son no negativas y menores que uno.
Probabilidad Total
Sean \(E\), \(F\) dos eventos entonces, \[P(E) = P(E\vert F)P(F) + P(E\vert F^c)P(F^c).\] De manera más general, sean \(F_i\) \(i = 1,...,n\) eventos mutuamente excluyentes cuya unión es el espacio muestral, entonces \[P(E) = \sum_{i=1}^n P(E\vert F_i)P(F_i).\]
La identidad anterior nos dice que podemos expresar la probabilidad del evento \(E\) como un promedio ponderado de probabilidades condicionales.
Ejemplo. Supongamos que una aseguradora clasifica a la gente en tres grupos de acuerdo a su nivel de riesgo: bajo, medio y alto. De acuerdo a los registros, las probabilidades de incurrir en un accidente en un laspo de un año son 0.05, 0.15 y 0.30 respectivamente. Si el 20% de la población se clasifica en riesgo bajo, 50% en medio y 30% en alto, ¿qué proporción de la población tiene un accidente en un año dado?
Para variables aleatorias tenemos:
Sean \(X\), \(Y\) dos variables aleatorias, podemos expresar la distribución marginal de \(X\) como: \[p_X(x) = \sum_{y} p_{X \vert Y}(x\vert y)p_Y(y).\]
Supongamos que ruedo un dado, si observo un número par lanzo una moneda justa (la probabilidad de observar águila es la misma que la de observar sol), si el dado muestra un número impar lanzo una moneda sesgada en la que la probabilidad de observar águila es 0.9. Si observo sol, ¿Cuál es la probabilidad de que haya lanzado la moneda sesgada?
El ejercicio anterior introduce la noción de probabilidad inversa: inicialmente conozco la probabilidad de observar sol condicional a que la moneda es sesgada pero ahora me interesa conocer la probabilidad de que haya lanzado una moneda sesgada una vez que observé un sol en el volado.
Regla de Bayes
La regla de Bayes es una consecuencia de la definición de probabilidad condicional.
Sean \(F_i\) \(i = 1,...,n\) eventos mutuamente excluyentes cuya unión es el espacio muestral, entonces \[P(F_j\vert E) = \frac{P(E\vert F_j)P(F_j)}{\sum_{i=1}^n P(E\vert F_i)P(F_i)}\] esta identidad se conoce como la regla de Bayes.
Ejemplo. En el contexto del ejemplo de los seguros ahora nos hacemos la siguiente pregunta: si un asegurado tuvo accidentes en 2013, ¿cuál es la probabilidad de que clasifique en riesgo bajo?
Al igual que con probabilidad condicional, la Regla de Bayes tiene una definición análoga para variables aleatorias.
Sean \(X\), \(Y\) dos variables aleatorias, \[p_{X\vert Y}(x\vert y) = \frac{p_{Y\vert X}(y\vert x)p_X(x)}{p_Y(y)}.\]
Una consecuencia de la regla de Bayes es que cualquier distribución multivariada sobre \(n\) variables \(X_1,X_2,...X_n\) se puede expresar como:
\[p(x_1,x_2,...x_n) = p_{X_1}(x_1)p_{X_2\vert X_1}(x_2\vert x_1)p_{X_3\vert X_1X_2}(x_3\vert x_1x_2)···p_{X_n\vert X_1...X_{n-1}}(x_n\vert x_1...x_{n-1})\] esta igualdad se conoce como {egla de la cadena}.
Nótese que esta regla funciona para cualquier ordenamiento de las variables aleatorias.
Supongamos ahora que una compañía de seguros divide a la gente en dos clases: propensos a accidente (30% de las personas) y no propensos a accidente. En un año dado aquellos propensos a accidentes sufren un accidente con probabilidad 0.4, mientras que los del otro grupo sufren un accidente con probabilidad 0.2. ¿Cuál es la probabilidad de que un asegurado tenga un accidente en su segundo año condicional a que sufrió un accidente en el primer año?
Independencia
Los eventos \(E\), \(F\) son independientes sí y solo sí \[P(EF) = P(E)P(F)\]
De la definición de independencia se sigue que \(P(E\vert F) = P(E)\). Esto es, los eventos \(E\) y \(F\) son independientes si saber que uno de ellos ocurrió no afecta la probabilidad del otro. Utilizaremos la notación \(E\perp F\) que se lee “\(E\) es independiente de \(F\)”.
Dos variables aleatorias \(X\), \(Y\), son independientes sí y sólo sí \[p(x,y) = p_X(x)p_Y(y)\]
Más aún, \(X\) y \(Y\) son independientes sí y sólo sí \(p(x,y) \propto g(x)h(y)\), por lo que para demostrar independecia podemos omitir las constantes en la factorización de las densidades
Similar a la independencia en eventos, la independencia de variables aleatorias implica que \(p_{X\vert Y}(x\vert y) = p_X(x)\), esto es, \(Y = y\) no provee información sobre \(X\).
Ejemplo. Consideremos la función de densidad conjunta \(p(x,y) = \frac{1}{384} x^2y^4e^{-y-(x/2)}\), \(x>0\), \(y>0\), ¿\(X\) y \(Y\) son independientes?
Podemos definir \[ g(x) = \left\{ \begin{array}{lr} x^2e^{-x/2} & : x > 0\\ 0 & : x \le 0 \end{array} \right. \] y \[ h(y) = \left\{ \begin{array}{lr} y^4e^{-y} & : x > 0\\ 0 & : x \le 0 \end{array} \right. \] entonces \(p(x,y) \propto g(x)h(y)\), para toda \(x\), \(y\) \(\in \mathbb{R}\) y concluímos que \(X\) y \(Y\) son independientes.
**Ejemplo.*. Si la densidad conjunta de \(X\) y \(Y\) está dada por: \[ p(x, y) = \left\{ \begin{array}{lr} 2 & : 0 < x < y, 0 < y < 1\\ 0 & : e.o.c. \end{array} \right. \]
¿\(X\) y \(Y\) son independientes?
Ejercicio. Recordando el ejemplo de departamentos en Hong Kong, veamos si Renta y Tipo son independientes, para esto comparemos \(p(renta|tipo)\) y \(p(renta)\).
Independencia condicional
La independencia de eventos o variables aleatorias es poco común en la práctica, más frecuente es el caso en que dos eventos son independientes dado un tercer evento.
Ejemplo. En una competencia de velocidad, cada atleta se somete a dos pruebas de dopaje que buscan detectar si el deportista ingirió una substania prohibida. La prueba A consiste en un examen de sangre y la prueba B en un exámen de orina, cada prueba se realiza en un laboratorio distinto y no hay intercambio de información entre los laboratorios. Es razonable pensar que los resultados de los dos exámenes no son independientes. Ahora, supongamos que sabemos que el atleta consumió la substancia prohibida, en este caso podemos argumentar que conocer el resultado de la prueba A no cambia la probabilidad de que el atleta salga positivo en la prueba B. Decimos que el resultado de la prueba B es condicionalmente independiente del resultado de la prueba A dado que el atleta consumió la substancia.
Sean \(A\), \(B\) y \(C\), tres eventos decimos que \(A\) es independiente de \(B\) condicional a \(C\) (\(A \perp B \vert C\)) si, \[ P(A,B\vert C) = P(A\vert C)P(B\vert C)\]
Similar al caso de independencia, \(A\) y \(B\) son condicionalmente independientes dado \(C\) sí y solo sí \(P(A \vert B,C) = P(A \vert C)\), esto es, una vez que conocemos el valor de \(C\), \(B\) no proporciona información adicional sobre \(A\).
Ejemplo. Retomemos el ejercicio de asegurados. En la solución de este ejercicio utilizamos que \(P(A_2|AA_1) = 0.4\) y que \(P(A_2|A^cA_1) = 0.2\), al establecer esa igualdad estamos asumiendo que \(A_2\) (el asegurado tiene un accidente en el año 2) y \(A_1\) (el asegurado tiene un accidente en el año 1) son eventos condicionalmente independientes dado \(A\) (el asegurado es propenso a accidentes): \(P(A_2|AA_1) = P(A_2|A) = 0.4\) y \(P(A_2|A^cA_1) = P(A_2|A^c) = 0.2\).
En el caso de variables aleatorias definimos independencia condicional como sigue.
Sean \(X\), \(Y\) y \(Z\), tres variables aleatorias decimos que \(X\) es independiente de \(Y\) condicional a \(Z\) (\(X \perp Y \vert Z\)) si y sólo sí, \[p(x,y\vert z) = p_{X\vert Z}(x\vert z)p_{Y\vert Z}(y\vert z).\]
Y tenemos que \(X\) es independiente de \(Y\) condicional a \(Z\) sí y sólo sí, \(p(x,y,z) \propto g(x,z)h(y,z)\).
Ejemplo. Supongamos que ruedo un dado, si observo un número par realizo dos lanzamientos de una moneda justa (la probabilidad de observar águila es la misma que la de observar sol), si el dado muestra un número impar realizo dos lanzamientos de una moneda sesgada en la que la probabilidad de observar águila es 0.9. Denotemos por \(Z\) la variable aleatoria asociada a la selección de la moneda, \(X_1\) la correspondiente al primer lanzamiento y \(X_2\) la correspondiente al segundo. Entonces, \(X_1\) y \(X_2\) no son independientes, sin embargo, son condicionalmente independientes (\(X_1 \perp X_2 \vert Z\)), puesto que una vez que se que moneda voy a lanzar el resultado del primer lanzamiento no aporta información adicional para el segundo lanzamiento. Calcularemos la distribución conjunta y la distribución condicional de \(X_2\) dado \(X_1\).
La distribución conjunta esta determinada por la siguiente tabla:
| Z | X1 | X2 | P(Z,X1,X2) |
|---|---|---|---|
| justa | a | a | 0.125 |
| justa | a | s | 0.125 |
| justa | s | a | 0.125 |
| justa | s | s | 0.125 |
| ses | a | a | 0.405 |
| ses | a | s | 0.045 |
| ses | s | a | 0.005 |
| ses | s | s | 0.045 |
La distribución condicional \(p(X_2|X_1)\) es,
| X1/X2 | a | s | . |
|---|---|---|---|
| a | 0.757 | 0.243 | 1 |
| s | 0.567 | 0.433 | 1 |
y la distribución condicional \(p(X_2|X_1,Z)=p(X_2|Z)\) es,
| X1/X2 | a | s | . |
|---|---|---|---|
| a | 0.5 | 0.5 | 1 |
| s | 0.9 | 0.1 | 1 |
En este punto es claro que \(X \perp Y \vert Z\) no implica \(X \perp Y\), pues como vimos en el ejemplo de las monedas \(X_1 \perp X_2 \vert Z\) pero \(X_1 \not \perp X_2\). Más aún, \(X \perp Y\) tampoco implica \(X \perp Y \vert Z\).
La independencia condicional tiene importantes consecuencias, por ejemplo, si \(X\) es independiente de \(Y\) dado \(Z\) entonces, \[p(x,y,z) = p_Z(z)p_{X\vert Z}(x\vert z)p_{Y\vert Z}(y\vert z).\]
Esta expresión de la densidad conjunta es similar a la que obtendríamos usando la regla de la cadena; sin embargo, el número de parámetros necesarios bajo esta representación es menor lo que facilita la estimación.
Modelos gráficos
El objetivo de esta sección es la simulación de modelos, una manera conveniente de simular de un modelo probabilístico es a partir del modelo gráfico asociado. Un modelo gráfico representa todas las cantidades involucradas en el modelo mediante nodos de una gráfica dirigida, el modelo representa el supuesto que dados los nodos padres \(padres(v)\) cada nodo es independiente del resto de los nodos a excepción de sus descendientes.
Los nodos en las gráficas se clasifican en 3 tipos:
Constantes fijas por el diseño del estudio, siempre son nodos sin padres.
Estocásticos son variables a los que se les asigna una distribución.
Determinísticos son funciones lógicas de otros nodos.
Los supuestos de independencia condicional que representa la gráfica implican que la distribución conjunta de todas las cantidades V tiene una factorización en términos de la distribución condicional \(p(v|padres(v))\) de tal manera que: \[p(V) = \prod p(v|padres(v))\]
Veamos como usar las gráficas para simular de modelos probabilísticos. Los siguientes ejemplos están escritos con base en el libro Data Analysis using Regression and Multilevel Hierarchical Models de Gelman y Hill.
Ejemplo de simulación discreta predictiva
La probabilidad de que un bebé sea niño o niña es 48.8% y 51.2% respectivamente. Supongamos que hay 400 nacimientos en un hospital en un año dado. ¿Cuántas niñas nacerán?
Comencemos viendo el modelo gráfico asociado.
La gráfica superior muestra todas las variables relevantes en el problema, y las dependencias entre ellas. En este caso \(n\) es una constante que representa el número de nacimientos, (\(n=400\)), \(p=48.8\) es la probabilidad de que un nacimiento resulte en niña y \(k \sim Binomial(p, n)\). Debido a que el número de éxitos (nacimientos que resultan en niña) depende de la tasa p y el número de experimentos n, los nodos que representan a éstas dos últimas variables están dirigidos al nodo que representa k.
Una vez que tenemos la gráfica es fácil simular del modelo:
library(ggplot2)
library(dplyr)
library(arm)
#> Warning: package 'lme4' was built under R version 3.4.2
library(tidyr)
set.seed(918739837)
#> Warning in set.seed(918739837): '.Random.seed' is not an integer vector but
#> of type 'NULL', so ignored
n_ninas <- rbinom(1, 400, 0.488)esto nos muestra algo que podría ocurrir en 400 nacimientos. Ahora, para tener una noción de la distribución simulamos el proceso 1000 veces:
sims_ninas <- rerun(1000, rbinom(1, 400, 0.488)) %>% flatten_dbl()
mean(sims_ninas)
#> [1] 196
sd(sims_ninas)
#> [1] 10.1
ggplot() + geom_histogram(aes(x = sims_ninas), binwidth = 3)
El histograma de arriba representa la distribución de probabilidad para el número de niñas y toma en cuenta la incertidumbre.
Podemos agregar complejidad al modelo, por ejemplo con probabilidad 1/125 un nacimiento resulta en gemelos fraternales, y para cada uno de los bebés hay una posibilidad de aproximadamente 49.5% de ser niña. Además la probabilidad de gemelos idénticos es de 1/300 y estos a su vez resultan en niñas en aproximadamente 40.5% de los casos.
Podemos simular 400 nacimientos bajo este modelo como sigue:
tipo_nacimiento <- sample(c("unico", "fraternal", "identicos"),
size = 400, replace = TRUE, prob = c(1 - 1 / 125 - 1 / 1300, 1 / 125, 1 / 300))
n_unico <- sum(tipo_nacimiento == "unico") # número de nacimientos únicos
n_fraternal <- sum(tipo_nacimiento == "fraternal")
n_identicos <- 400 - n_unico - n_fraternal
n_ninas <- rbinom(1, n_unico, 0.488) +
rbinom(1, 2 * n_fraternal, 0.495) + # en cada nacimiento hay 2 bebés
2 * rbinom(1, n_identicos, 0.495)
n_ninas
#> [1] 183Repetimos la simulación 1000 veces para aproximar la distribución de número de niñas en 400 nacimientos.
modelo2 <- function(){
tipo_nacimiento <- sample(c("unico", "fraternal", "identicos"),
size = 400, replace = TRUE, prob = c(1 - 1 / 125 - 1 / 1300, 1 / 125, 1 / 300))
# número de nacimientos de cada tipo
n_unico <- sum(tipo_nacimiento == "unico") # número de nacimientos únicos
n_fraternal <- sum(tipo_nacimiento == "fraternal")
n_identicos <- 400 - n_unico - n_fraternal
# simulamos para cada tipo de nacimiento
n_ninas <- rbinom(1, n_unico, 0.488) +
rbinom(1, 2 * n_fraternal, 0.495) + # en cada nacimiento hay 2 bebés
2 * rbinom(1, n_identicos, 0.495)
n_ninas
}
sims_ninas_2 <- rerun(1000, modelo2()) %>% flatten_dbl()
mean(sims_ninas_2)
#> [1] 197
ggplot() + geom_histogram(aes(x = sims_ninas_2), binwidth = 4)
Ejemplo de simulación continua predictiva
El 52% de los adultos en EUA son mujeres y el 48% hombres, las estaturas de los hombres se distribuyen aproximadamente normal con media 175 cm y desviación estándar de 7.37 cm, en el caso de las mujeres la distribución es aproximadamente normal con media 161.80 cm y desviación estándar de 6.86 cm. Supongamos que seleccionamos 10 adultos al azar. ¿qué podemos decir del promedio de estatura?
sexo <- rbinom(10, 1, 0.52)
altura <- rnorm(sexo, mean = 161.8 * (sexo == 1) + 175 * (sexo == 0),
sd = 6.86 * (sexo == 1) + 7.37 * (sexo == 0))
mean(altura)
#> [1] 170Simulamos la distribución de la altura promedio:
mediaAltura <- function(){
sexo <- rbinom(10, 1, 0.52)
altura <- rnorm(sexo, mean = 161.8 * (sexo == 1) + 175 * (sexo == 0),
sd = 6.86 * (sexo == 1) + 7.37 * (sexo == 0))
}
sims_alturas <- rerun(1000, mediaAltura())
media_alturas <- sims_alturas %>% map_dbl(mean)
mean(media_alturas)
#> [1] 168
sd(media_alturas)
#> [1] 3.11
ggplot() + geom_histogram(aes(x = media_alturas), binwidth = 1.2)
¿Y que podemos decir de la altura máxima?
alt_max <- sims_alturas %>% map_dbl(max)
qplot(alt_max, geom = "histogram", binwidth = 1.5)
Supongamos que una compañía cambia la tecnología usada para producir una cámara, un estudio estima que el ahorro en la producción es de $5 por unidad con un error estándar de $4. Más aún, una proyección estima que el tamaño del mercado (esto es, el número de cámaras que se venderá) es de 40,000 con un error estándar de 10,000. Suponiendo que las dos fuentes de incertidumbre son independientes, usa simulación para estimar el total de dinero que ahorrará la compañía, calcula un intervalo de confianza.
Ejemplo de simulación de un modelo de regresión
En regresión utilizamos simulación para capturar tanto la incertidumbre en la predicción (término de error en el modelo) como la incertidumbre en la inferencia (errores estándar de los coeficientes e incertidumbre del error residual).
Comenzamos con un ejemplo en el que simulamos únicamente incertidumbre en la predicción.
Supongamos que el puntaje de un niño de tres años en una prueba cognitiva esta relacionado con las características de la madre, el siguiente modelo resume la diferencia en los puntajes promedio de los niños cuyas madres se graduaron de preparatoria y los que no.
\[y_i= \beta_0 + \beta_1 X_{i1} + \epsilon_i\]
donde \(y_i\) es el puntaje del i-ésimo niño, \(X_{i1}\) es una variable binaria que indica si la madre se graduó de preparatoria (codificado como 1) o no (codificado como 0) y \(\epsilon_i\) son los error aleatorios, estos son independientes con distribución normal \(\epsilon_i \sim N(0, \sigma^2)\).
Ahora consideremos el problema de simular el puntaje de 50 niños 30 con madres que terminaron la preparatoria y 20 cuyas madres no terminaron. Los coeficientes que usaremos son:
\[\beta_0 = 78\] \[\beta_1 = 12\] \[\sigma = 20\]
El modelo gráfico asociado sería como sigue:
vector_mu <- c(rep(78 + 12, 30), rep(78, 20)) # beta_0 + beta_1 X
y <- rnorm(50, vector_mu, 20)
sims_y <- rerun(2000, rnorm(50, vector_mu, 20))Podemos calcular la media y su intervalo de confianza:
medias <- sims_y %>% map_dbl(mean)
quantile(medias, c(0.025, 0.975))
#> 2.5% 97.5%
#> 79.5 90.9
qplot(medias, geom = "histogram", binwidth = 1.5)
Supongamos ahora que nos interesa incorporar que tenemos incertidumbre en los coeficientes de regresión, y expresamos nuestra incertidumbre a través de distribuciones de probabilidad, ¿cómo sería el modelo gráfico asociado?
Primero suponemos que \(\sigma^2\) tiene una distribución centrada en 202, proporcional a una distribución \(\chi^2\) con 432 grados de libertad.
\[ \begin{eqnarray*} \begin{pmatrix}\beta_{0}\\ \beta_{1} \end{pmatrix} & \sim & N\left[\left(\begin{array}{c} 77\\ 12 \end{array}\right), \sigma^2 \left(\begin{array}{cc} 0.01 & -0.01\\ -0.01 & 0.01 \end{array}\right)\right] \end{eqnarray*} \]
Finalmente, simulemos del modelo incorporando tanto la incertidumbre correpondiente a la predicción como la incertidumbre en los coeficientes de regresión.
Simula \(\sigma=20\sqrt{(432)/X}\) donde \(X\) es una generación de una distribución \(\chi^2\) con \(432\) grados de libertad.
Dado \(\sigma\) (obtenido del paso anterior), simula \(\beta\) de una distribución normal multivariada con media \((77,12)\) y matriz de covarianzas \(\sigma^2 V\).
Simula \(y\) el vector de observaciones usando los parámetros de 1 y 2.
simula_parametros <- function(){
# empezamos simulando sigma
sigma <- 20 * sqrt((432) / rchisq(1, 432))
# la usamos para simular betas
beta <- MASS::mvrnorm(1, mu = c(78, 12),
# Sigma = sigma ^ 2 * matrix(c(4.2, -4.2, -4.2, 5.4), nrow = 2))
Sigma = sigma ^ 2 * matrix(c(0.011, -0.011, -0.011, 0.013), nrow = 2))
# Simulamos las observaciones
list(sigma = sigma, beta = beta)
}
sims_parametros <- rerun(2000, simula_parametros())
# simulamos los puntajes
simula_puntajes <- function(beta, sigma, n_hs = 30, n_nhs = 20){
vector_mu <- c(rep(beta[1] + beta[2], n_hs), rep(beta[1], n_nhs)) # beta_0 + beta_1 X
obs = rnorm(50, vector_mu, sigma)
}
sims_puntajes <- map(sims_parametros, ~simula_puntajes(beta = .[["beta"]], sigma = .[["sigma"]]))
medias_incert <- sims_puntajes %>% map_dbl(mean)
quantile(medias_incert, c(0.025, 0.975))
#> 2.5% 97.5%
#> 79.8 91.0
qplot(medias_incert, geom = "histogram", binwidth = 1)
Los parametros se obtuvieron de ajustar el modelo de regresión lineal.
library(foreign)
kids_iq <- read.dta("data/kidiq.dta")
lm_kid <- lm(kid_score ~ mom_hs, kids_iq)
summary(lm_kid)
#>
#> Call:
#> lm(formula = kid_score ~ mom_hs, data = kids_iq)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -57.55 -13.32 2.68 14.68 58.45
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 77.55 2.06 37.67 <2e-16 ***
#> mom_hs 11.77 2.32 5.07 6e-07 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 19.9 on 432 degrees of freedom
#> Multiple R-squared: 0.0561, Adjusted R-squared: 0.0539
#> F-statistic: 25.7 on 1 and 432 DF, p-value: 5.96e-07
V <- vcov(lm_kid) / 20^2summary(lm_kid)
#>
#> Call:
#> lm(formula = kid_score ~ mom_hs, data = kids_iq)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -57.55 -13.32 2.68 14.68 58.45
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 77.55 2.06 37.67 <2e-16 ***
#> mom_hs 11.77 2.32 5.07 6e-07 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 19.9 on 432 degrees of freedom
#> Multiple R-squared: 0.0561, Adjusted R-squared: 0.0539
#> F-statistic: 25.7 on 1 and 432 DF, p-value: 5.96e-07Podemos usar simulación para calcular intervalos de confianza para \(\beta_0\) y \(\beta_1\),
sims_parametros %>% map_dbl(~(.$beta[1])) %>% sd()
#> [1] 2.06
sims_parametros %>% map_dbl(~(.$beta[2])) %>% sd()
#> [1] 2.27No parece que valga la pena el esfuerzo cuando podemos calcular los intervalos analíticamente, sin embargo con simulación podemos responder fácilmente otras preguntas, por ejemplo, la pregunta inicial: ¿cuál es la media esperada para un conjunto de 50 niños, 30 con madres que hicieron preparatoria y 20 que no? es fácil de responder con simulación.
Podríamos usar predict() para calcular el estimador puntual de la media en el examen para los niños:
pred_mi_pob <- predict(lm_kid, newdata = data.frame(mom_hs = c(rep(1, 30),
rep(0, 20))), se.fit = TRUE)
mean(pred_mi_pob$fit)
#> [1] 84.6¿Cómo calculas el error estándar? En este caso se puede resolver pues es una combinación lineal de los coeficientes del modelo, pero en muchos casos nuestro objetivo es más complicado que coeficientes o combinaciones lineales de estos.
Simulación de predicciones no lineales
Veamos un ejemplo de las elecciones en el congreso de EUA. Tenemos un modelo que usaremos para predecir la elección de 1990 basados en la de 1988.
Explicación del problema. EUA está dividido en 435 distritos congresionales, definimos la variable de interés \(y_i\) con \(i=1,...,n\), como la participación del partido Demócrata en el distrito \(i\) en \(1988\). La participación se calcula como el porcentaje de los votos correspondientes a los demócratas del total de votos que recibieron los demócratas y republicanos, esto es, se excluyen los votos a otros partidos.
El modelo del que simularemos se construyó usando datos de 1986 y 1988.
# Los datos están almacenados en 3 archivos correspondientes al año
paths <- dir("data/congress", full.names = TRUE)
paths <- set_names(paths, basename(paths))
# Leemos los datos y recodificamos variables
congress <- map_dfr(paths, read.table, quote = "\"",
stringsAsFactors = FALSE, .id = "year") %>%
mutate(id = rep(1:435, 3),
year = parse_number(year), # año de la elección
incumbency = ifelse(V3 == -9, NA, V3), # codificar NAs
dem_share = V4 / (V4 + V5)) %>% # participación demócrata
dplyr::select(id, year, incumbency, dem_share)
glimpse(congress)
#> Observations: 1,305
#> Variables: 4
#> $ id <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, ...
#> $ year <dbl> 1986, 1986, 1986, 1986, 1986, 1986, 1986, 1986, 198...
#> $ incumbency <int> 1, 1, 1, -1, 1, -1, 0, -1, -1, 1, 1, 1, 1, 1, 1, 0,...
#> $ dem_share <dbl> 0.745, 0.674, 0.696, 0.465, 0.391, 0.358, 0.549, 0....
# datos en forma horizontal
congress_share <- spread(dplyr::select(congress, -incumbency), year, dem_share)
colnames(congress_share) <- c("id", "share_86", "share_88", "share_90")
congress_inc <- spread(dplyr::select(congress, -dem_share), year, incumbency)
colnames(congress_inc) <- c("id", "inc_86", "inc_88", "inc_90")
# quitamos NAs
congress_h <- na.omit(left_join(congress_share, congress_inc))
#> Joining, by = "id"
glimpse(congress_h)
#> Observations: 431
#> Variables: 7
#> $ id <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16...
#> $ share_86 <dbl> 0.745, 0.674, 0.696, 0.465, 0.391, 0.358, 0.549, 0.22...
#> $ share_88 <dbl> 0.772, 0.636, 0.665, 0.274, 0.264, 0.334, 0.632, 0.33...
#> $ share_90 <dbl> 0.714, 0.597, 0.521, 0.234, 0.477, 0.256, 0.602, 0.49...
#> $ inc_86 <int> 1, 1, 1, -1, 1, -1, 0, -1, -1, 1, 1, 1, 1, 1, 1, 0, 1...
#> $ inc_88 <int> 1, 1, 1, -1, -1, -1, 1, -1, -1, 1, 1, 1, 1, 1, 1, 1, ...
#> $ inc_90 <int> 1, 1, 0, -1, 0, -1, 0, -1, -1, 1, 1, 1, 1, 1, 1, 1, 1...
ggplot(congress_h, aes(x = share_86, y = share_88, color = factor(inc_86))) +
geom_abline(color = "darkgray") +
geom_point()
# quitamos las elecciones que no se compitieron
congress_h <- filter(congress_h, share_88 != 1 & share_88 != 0)
fit_88 <- lm(share_88 ~ share_86 + inc_88, data = congress_h)
summary(fit_88)
#>
#> Call:
#> lm(formula = share_88 ~ share_86 + inc_88, data = congress_h)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -0.2416 -0.0428 -0.0059 0.0523 0.2263
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 0.3004 0.0151 19.9 <2e-16 ***
#> share_86 0.3858 0.0278 13.9 <2e-16 ***
#> inc_88 0.1043 0.0069 15.1 <2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 0.0774 on 351 degrees of freedom
#> Multiple R-squared: 0.846, Adjusted R-squared: 0.845
#> F-statistic: 965 on 2 and 351 DF, p-value: <2e-16Simulemos del modelo:
# Matriz X
X <- cbind(1, congress_h$share_88, congress_h$inc_90)
simula_modelo <- function(){
sigma <- 0.08 * sqrt((351) / rchisq(1, 351))
beta <- MASS::mvrnorm(1, mu = c(0.30, 0.39, 0.10),
Sigma = sigma ^ 2 * matrix(c(0.04, -0.07, 0.01, -0.07, 0.13, -0.02, 0.01,
-0.02, 0.01), nrow = 3))
mu <- X %*% beta
data_frame(id = 1:358, dem_share = rnorm(358, mu, sigma))
}
sims_congress <- rerun(2000, simula_modelo())
sims_congress <- sims_congress %>% bind_rows() %>% mutate(sim = rep(1:2000, each = 358))
ggplot(sims_congress, aes(x = reorder(id, dem_share), y = dem_share)) +
geom_boxplot() +
geom_hline(color = "red", yintercept = 0.5)
Podemos preguntarnos cuántas elecciones ganaron los demócratas en 1990: \(\sum I(\grave{y} > 0.5)\)
sims_congress %>%
group_by(sim) %>%
mutate(wins = sum(dem_share > 0.5)) %>%
ungroup() %>%
summarise(mean_wins = mean(wins),
median_wins = median(wins),
sd_wins = sd(wins))Veamos lo que ocurrió realmente
sum(congress_h$share_90 > 0.5)Podemos usar la función sim() del paquete arm para simular de modelos lineales y lineales generalizados.
Inferencia de gráficas
¿Qué observas en las siguientes gráficas?
El ímpetu por concluir (rage-to-conclude bias, Tufte) nos hace ver patrones en datos donde no existen dichos patrones. Esto puede conllevar a inferencias prematuras, el análisis estadístico busca alinear la realidad en la evidencia con la inferencia que se realiza a partir de dicha evidencia.
Ahora, hay 2 componentes en la inferencia: evaluar si existe una diferencia y que tan grande es ésta. En esta sección estudiaremos la primera. En particular buscamos responder la pregunta, ¿Lo que veo en una gráfica de la muestra refleja a la población entera? Las ideas que veremos, así como algunos ejemplos, se tomaron del artículo Graphical Inference for Infovis de Wickham, et. al.
Comenzemos con un ejemplo, consideremos un conjunto de datos que registra estaturas de hombres y mujeres. Supongamos que nos interesa describir de manera simple los datos, independientemente de si se trata de un hombre o una mujer.
head(singer_g)
#> gender height
#> 1 F 163
#> 2 F 157
#> 3 F 168
#> 4 F 165
#> 5 F 152
#> 6 F 155
ggplot(singer_g, aes(x = gender, y = height)) +
geom_jitter(position = position_jitter(width = 0.1, height = 1))
Suponemos que la estatura es una medición que se distribuye aproximadamente normal con media 171 cm y desviación estándar 10 cm. ¿Es razonable esta descripción?
mean(singer_g$height)
#> [1] 171
sd(singer_g$height)
#> [1] 9.71Una manera de probar que tan buena es esta descripción es considerando qué es lo que veríamos si el modelo es el que acabamos de mencionar, para esto hacemos 19 simulaciones bajo el modelo \(N(\mu, \sigma^2)\) y las comparamos con los datos observados. ¿Captura este modelo las características observadas?
library(nullabor)
sing_null <- lineup(null_dist('height', dist = 'normal',
params = list(mean = 171, sd = 10)), n = 20, singer_g)
#> decrypt("WOXC 1IZI Tx t4nTZT4x K")
ggplot(sing_null, aes(x = gender, y = height)) +
facet_wrap(~ .sample) +
geom_jitter(position = position_jitter(width = 0.1, height = 1),
size = 0.8, alpha = 0.5)
En la figura anterior, 19 de las gráficas corresponden a datos aleatorios (nulos) y una gráfica corresponde a datos reales. Si se distingue la gráfica verdadera de las nulas entonces hay evidencia de que es distinto. En realidad, el poder distinguirlo provee evidencia estadística rigurosa de que hay diferencia.
Antes de proseguir recordemos los conceptos de prueba de hipótesis:
Hipótesis nula (\(H_0\)): hipótesis que se desea contrastar, describe la conducta default del fenómeno de interés.
Hipótesis alternativa (\(H_1\)).
Estadística de prueba: es una estadística en en base a la cuál tomamos la decisión de rechazar o no rechazar. Se calcula considerando la hipótesis nula como verdadera.
4.. Valor-p: Nivel de significacncia alcanzado, el nivel de significancia más pequeño al cual los datos observados indican que la hipótesis nula debe ser rechazada.
| Escenario | \(H_0\) verdadera | \(H_0\) Falsa |
|---|---|---|
| Rechazar \(H_0\) | Error Tipo 1 (\(\alpha\)) | Decisión correcta |
| No rechazar \(H_0\) | Decisión correcta | Error tipo 2 (\(\beta\)) |
Veamos una prueba clásica usando el ejemplo de estaturas. \[H_0:\mu_m = \mu_h\] Esto es, la hipótesis nula corresponde a que la estura media de hombres y mujeres es la misma. La estadistica de prueba es: \[Z=\frac{\bar{X_m}-\bar{X_h}}{\hat{\sigma}\sqrt{1/n_1+1/n_2}}\] la prueba se basa en una distribución \(t\) con \(n_1 + n_2 -1\) grados de libertad.
t.test(x = singer_g$height[singer_g$gender == "F"],
y = singer_g$height[singer_g$gender == "M"], var.equal = TRUE)
#>
#> Two Sample t-test
#>
#> data: singer_g$height[singer_g$gender == "F"] and singer_g$height[singer_g$gender == "M"]
#> t = -20, df = 200, p-value <2e-16
#> alternative hypothesis: true difference in means is not equal to 0
#> 95 percent confidence interval:
#> -16.0 -12.6
#> sample estimates:
#> mean of x mean of y
#> 164 179
# veamos como se ven los "inocentes"
nulos <- data.frame(t = rt(10000, 233))
ggplot(nulos, aes(x = t)) +
geom_histogram(color = "darkgray", fill = "darkgray") +
geom_vline(xintercept = c(qt(0.025, 233), qt(0.975, 233)), color = "red",
alpha = 0.5)
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Los principios de pruebas de hipótesis son los mismos para pruebas visuales, a excepción de dos aspectos: la estadística de prueba y el mecanismo para medir similitud. La estadística de prueba ahora es una gráfica de los datos, y en lugar de una diferencia matemática usamos el ojo humano.
En la prueba gráfica, los verdaderos datos están escondidos entre 19 gráficas de datos nulos, donde un dato nulo es una muestra de la distribución bajo la hipótesis nula. Si es posible identificar los datos, hay evidencia indicando que estos son distintos a los datos nulos.
Este ejemplo muestra la analogía entre una prueba tradicional (numérica) y una prueba visual. Las pruebas estadísticas tradicionales se han estudiado por un largo tiempo y funcionan bien cuando los datos se comportan bien: siguiendo una distribución conocida en un escenario relativamente simple. Sin embargo, la estadística tradicional no cubre todas las complejidades que surgen cuando se exploran datos y la ventaja de las pruebas visuales es que se pueden usar en escenarios de análisis complejos donde no existen pruebas numéricas.
Protocolo de pruebas visuales:
Genera n-1 decoys (datos que siguen la hipótesis nula)
Grafica los decoy + los datos reales, donde los datos están posicionados de manera aleatoria.
Muestra la gráfica a un observador imparcial.
¿Pueden distinguir los datos? Si es el caso, hay evidencia de diferencia verdadera (p-value = 1/n).
Volvamos al ejemplo de las estaturas, proponemos el siguiente modelo: la estatura es aproximadamente normal con media 179 para hombres y 164 para mujeres, la desviación estándar es de 6.5 en ambos casos.
singer_g %>%
group_by(gender) %>%
summarise(round(mean(height)), round(sd(height), 1))
#> # A tibble: 2 x 3
#> gender `round(mean(height))` `round(sd(height), 1)`
#> <chr> <dbl> <dbl>
#> 1 F 164 6.3
#> 2 M 179 6.9¿Es razonable esta descripción?
# centramos los datos con la media correspondiente
singer_c <- singer_g %>%
group_by(gender) %>%
mutate(height_c = height - mean(height))
set.seed(26832)
sing_null_c <- lineup(null_dist('height_c', dist = 'normal',
params = list(mean = 0, sd = sd(singer_c$height_c))), n = 20, singer_c)
#> decrypt("WOXC 1IZI Tx t4nTZT4x N")
ggplot(sing_null_c, aes(x = gender, y = height_c)) +
facet_wrap(~ .sample) +
geom_jitter(position = position_jitter(width = 0.1, height = 1),
size = 0.8, alpha = 0.5)
Diagramas de dispersión
Un diagrama de dispersión muestra la relación entre dos variables continuas y responde a la pregunta: ¿existe una relación entre \(x\) y \(y\)? Una posible hipótesis nula es que no hay relación entre las variables. Supongamos que queremos usar preubas visuales para esta hipótesis, la función null_permute del paquete nullabor recibe el nombre de una variable de los datos y la salida de la funcción consiste en la variable permutada para obtener los datos nulos.
En este ejercicio usarás los datos diamonds, toma una muestra de tamaño 5000 (sin reemplazo) y procede como se indica:
Usa la función null_permute para crear una nueva base de datos con la variable depth permutada en los datos nulos.
Usa la función line_up (como se usó en el caso de las estaturas) y genera una gráfica con 20 páneles, donde grafiques depth en el eje horizontal y carat en el eje vertical tienes evidencia para afirmar que existe una relación entre depth y carat.
En muchos casos el supuesto de independencia es demasiado fuerte, es claro que las variables están relacionadas y queremos estudiar una relación particular. Por ejemplo, podemos pensar que los intentos de encestar a tres puntos en el basquetbol siguen una distribución cuadrática en el espacio: mientras el ángulo entre el jugador y la canasta aumenta el jugador se acerca más para asegurar el éxito. La siguiente figura prueba esta hipótesis usando datos de los tres punteros de los Lakers en la temporada 2008/2009, los datos se obtuvieron de Basketball Geek.
paths <- dir("data/2008-2009", pattern = "LAL", full.names = TRUE)
basket <- purrr::set_names(paths, paths) %>%
map_df(read_csv)
basket_LA <- basket %>%
filter(team == "LAL", type == "3pt", !is.na(x), !is.na(y)) %>% # datos Lakers
mutate(
x = x + runif(length(x), -0.5, 0.5),
y = y + runif(length(y), -0.5, 0.5),
r = sqrt((x - 25) ^ 2 + y ^ 2), # distancia a canasta
angle = atan2(y, x - 25)) %>% # ángulo
filter(r > 20 & r < 39) %>% # lanzamientos en el rango típico
dplyr::select(x, y, r, angle)
# guardar datos
write.table(basket_LA, file = "data/basket_LA.csv", sep = ",")glimpse(basket_LA)
#> Observations: 1,410
#> Variables: 4
#> $ x <dbl> 1.137, 41.175, 45.635, 37.477, 43.171, 7.752, 15.801, 0....
#> $ y <dbl> 5.96, 26.14, 17.24, 30.15, 24.28, 27.82, 28.73, 5.39, 29...
#> $ r <dbl> 24.6, 30.7, 26.9, 32.6, 30.3, 32.7, 30.2, 24.6, 33.6, 32...
#> $ angle <dbl> 2.897, 1.017, 0.696, 1.178, 0.928, 2.126, 1.881, 2.921, ...
basket_null <- lineup(null_lm(r ~ poly(angle, 2)), basket_LA, n = 10)
#> decrypt("WOXC 1IZI Tx t4nTZT4x p")
ggplot(basket_null, aes(x = angle * 180 / pi, y = r)) +
geom_point(alpha = 0.5, size = 0.8) +
scale_x_continuous("Angle (degrees)",
breaks = c(0, 45, 90, 135, 180),
limits = c(0, 180)) +
facet_wrap(~ .sample, nrow = 2)
Los datos reales están escondidos entre un conjunto de datos nulos que siguen la hipótesis de una relación cuadrática, los conjuntos nulos se crean ajustando el modelo, produciendo predicciones y residuales y sumando los residuales rotados a las predicciones.
Hasta ahora hemos usado las funciones del paquete nullabor, repitamos este ejemplo escribiendo las instrucciones nosotros mismos, en este caso, generamos los residuales como observaciones de una normal con los parámetros estimados.
fit <- lm(r ~ poly(angle, 2), data = basket_LA)
n <- 10
# matriz de covariables
X <- model.matrix(r ~ poly(angle, 2), basket_LA)
n_obs <- nrow(X)
# 9 simulaciones del modelo, cada una del mismo tamaño que los datos
sims <- rep(X %*% coef(fit), n - 1) + rnorm((n - 1) * n_obs, mean = 0, sd = 1)
# creo una base de datos que incluya datos originales
basket_LA_3 <- data.frame(
angle = rep(basket_LA$angle, n), # el valor de x (ángulo) para cada panel
r = c(basket_LA$r, sims), # distancias simuladas y observadas
id = rep(sample(1:n, size = n), each = n_obs)) # id aleatorio para esconder los datos
ggplot(basket_LA_3, aes(x = angle * 180 / pi, y = r)) +
geom_point(alpha = 0.5, size = 0.8) +
scale_x_continuous("Angle (degrees)",
breaks = c(0, 45, 90, 135, 180),
limits = c(0, 180)) +
facet_wrap(~ id, nrow = 2)
Otra alternativa es graficar los residuales, de la especificación del modelo esperamos que los residuales se distribuyan \(N(0, 1)\) por lo que los datos nulos consisten simplemente en muestras de una normal estándar.
basket_LA$resid <- resid(fit)
basket_null_2 <- lineup(
null_dist("resid", "normal", list(mean = 0, sd = 1)),
basket_LA, n = 10)
#> decrypt("WOXC 1IZI Tx t4nTZT4x N")
ggplot(basket_null_2, aes(x = angle * 180 / pi, y = resid)) +
geom_point(alpha = 0.6, size = 0.8) +
scale_x_continuous("Angle (degrees)",
breaks = c(0, 45, 90, 135, 180),
limits = c(0, 180)) +
facet_wrap(~ .sample, nrow = 2)
El paquete nullabor incluye métodos para generar distintas pruebas visuales: generar datos nulos de una distribución particular, datos nulos con residuales de un modelo y permutando una variable. Es fácil extender las ideas de pruebas visuales a otros modelos, por ejemplo, podemos aplicar el protocolo a un modelo de splines que veremos más adelante.
library(fda)
#> Loading required package: splines
#>
#> Attaching package: 'fda'
#> The following object is masked from 'package:graphics':
#>
#> matplot
knots <- quantile(x) # nudos
base <- create.bspline.basis(
norder = 4, # polinomios cúbicos
breaks = knots # nodos en los cuartiles de x
) # base de splines
H <- eval.basis(x, base) # y = H*beta + epsilon
beta_hat <- as.vector(solve(t(H) %*% H) %*% t(H) %*% toy$y)
mu <- function(x, betas){
as.numeric(betas %*% t(eval.basis(x, base)))
}
mu_hat <- as.numeric(beta_hat %*% t(H))
sigma_hat <- sqrt(1 / n * sum((toy$y - mu_hat) ^ 2))
# simulamos 9 conjuntos de datos de acuerdo a nuestro modelo.
y_sim <- rep(mu_hat, 9) + rnorm(n * 9, sd = sigma_hat)
codigo <- sample(1:10, 10) # para poder identificar los datos
toy_null <- data.frame(x = rep(x, 10), y = c(y, y_sim),
id = rep(codigo, each = n))
ggplot(toy_null, aes(x = x, y = y)) +
geom_point(size = 0.8) +
facet_wrap(~ id, nrow = 2)
Las pruebas de hipótesis visuales, no son la única herramienta que se debe usar para evaluar un modelo, por lo que no debemos descartar un modelo basados únicamente en una prueba visual. Sin embargo, las pruebas visuales y en general graficar los modelos ajustados, si son una herramienta útil que nos puede ayudar a comprender las implicaciones de un modelo y las fallas del mismo.
Cálculos de poder estadístico
Cuando se esta diseñando un estudio se determina la precisión en las inferencias que se desea, y esto determina el tamaño de muestra que se tomará. Usualmente se fija uno de los siguientes dos objetivos:
Se determina el error estándar de un parámetro o cantidad de interés. Por ejemplo en encuestas electorales es típico reportar los resultados de esta encuesta más menos 3 puntos porcentuales tienen un nivel del 95% de confianza, cúantas personas se debe entrevistar para lograr esto?
Se determina la probabilidad de que un estadístico determinado sea estadísticamente significativo. Por ejemplo cuando se hacen ensayos clínicos se determina un tamaño de muestra para que con probabilidad de x% se detecte una diferencia clinicamente relevante con el nuevo tratamiento (si es que este es efectivo).
En cualquiera de estos dos escenarios se necesita hacer supuestos para poder calcular el tamaño de muestra.
Tamaño de muestra para un error estándar determinado
Supongamos que queremos estimar el porcentaje de la población que apoya la pena de muerte. Sospechamos que la proporción es 60%, imaginemos que queremos una precisión (error estándar) de a lo más 0.05, o 5 puntos porcentuales. Bajo muestreo aleatorio simple, para una muestra de tamaño \(n\), el error estándar de la proporción \(p\) es \[\sqrt{p(1-p)/n}\] Sustituyendo nuestra expectativa \(p = 0.60\) llegamos a que el error estándar sería \(0.49/\sqrt{n}\), de tal manera que si queremos \(se(p) \le 0.05\) necesitamos \(n>0.96\), en el caso de proporciones es fácil determinar el tamaño de muestra de manera conservadora pues bastas con suponer \(p = 0.5\).
se_fun_n <- function(n, p) sqrt(p * (1-p) / n)
xy <- data.frame(x = 20:220, y = seq(0, 1, 0.005))
ggplot(xy, aes(x = x, y = y)) +
stat_function(fun = se_fun_n, args = list(p = 0.7), aes(color = "p=0.7")) +
stat_function(fun = se_fun_n, args = list(p = 0.9), aes(color = "p=0.9")) +
stat_function(fun = se_fun_n, args = list(p = 0.5), aes(color = "p=0.5")) +
labs(x = "n", y = "se", color = "") +
geom_segment(x = 20, xend = 100, y = 0.05, yend = 0.05, color = "red",
alpha = 0.3, linetype = "longdash") +
geom_segment(x = 100, xend = 100, y = 0.05, yend = 0, color = "red",
alpha = 0.3, linetype = "longdash")
Cómo calcularíamos el tamaño de muestra simulando? Este caso es trivial calcular pero conforme se aumenta complejidad en el diseño de la muestra o en el estadístico de interés también se complica encontrar una solución analítica.
sim_p_hat <- function(n, p, n_sims = 1000){
sim_muestra <- rbinom(n_sims, size = n, prob = p)
se_p_hat <- sd(sim_muestra / n)
data_frame(n = n, se_p_hat = se_p_hat, p = p)
}
sims_.7 <- map_df(seq(20, 220, 10), ~sim_p_hat(n = ., p = 0.7))
sims_.5 <- map_df(seq(20, 220, 10), ~sim_p_hat(n = ., p = 0.5))
sims_.6 <- map_df(seq(20, 220, 10), ~sim_p_hat(n = ., p = 0.6))
sims_.9 <- map_df(seq(20, 220, 10), ~sim_p_hat(n = ., p = 0.9))
sims <- bind_rows(sims_.7, sims_.5, sims_.6, sims_.9)
ggplot(sims, aes(x = n, y = se_p_hat, color = factor(p), group = p)) +
geom_smooth(se = FALSE, size = 0.5) +
labs(x = "n", y = "se", color = "") +
geom_segment(x = 20, xend = 100, y = 0.05, yend = 0.05, color = "red",
alpha = 0.3, linetype = "longdash") +
geom_segment(x = 100, xend = 100, y = 0.05, yend = 0, color = "red",
alpha = 0.3, linetype = "longdash")
#> `geom_smooth()` using method = 'loess' and formula 'y ~ x'
Tamaño de muestra determinado para obtener significancia estadística con una probabilidad determinada
Supongamos que nuestro objetivo es demostrar que más de la mitad de la población apoya la pena de muerte, esto es \(p>0.5\), nuevamente tenemos la hipótesis que el verdadero valor es \(p=0.6\).
Un prueba de potencia típica tiene un poder de 80%, es decir nos gustaría seleccionar \(n\) tal que 80% de los intervalos construidos con 95% de confianza no incluyan 0.5. Para encontrar la \(n\) tal que el 80% de las estimaciones estén al menos, 1.96 errores estándar por encima de 0.5 necesitamos que:
\[0.5 + 1.96 se = 0.6 - 0.84 se\]
Sustituyendo \(se = 0.5/\sqrt(n)\) obtenemos \(n=196\)
Y simulando sería
sim_potencia <- function(n, p, n_sims = 1000){
sim_muestra <- rbinom(n_sims, size = n, prob = p)
se_p_hat <- sd(sim_muestra / n)
acepta <- (sim_muestra / n - 1.96 * se_p_hat) > 0.5
data_frame(n = n, potencia = mean(acepta))
}
sims <- map_df(c(2, 10, 50, 80, 100, 150, 200, 300, 500), ~sim_potencia(n = ., p = 0.6))
ggplot(sims) +
geom_line(aes(x = n, y = potencia)) +
geom_hline(yintercept = 0.8, color = "red", alpha = 0.5)
Veamos un ejemplo más interesante, tenemos medidas del sistema inmune (porcentaje de CD4 en transformado con raíz cuadrada) de niños VIH positivos a lo largo de un periodo de 2 años. Las series de tiempo se ajustan de manera razonable con un modelo de intercepto y pendiente variable:
\[y_i \sim N(\alpha_{j[i]} + \beta_{j[i]}t_i, \sigma^2_y)\] donde \(i\) indexa las mediciones tomadas al tiempo \(i\) en el individuo \(j[i]\).
library(lme4)
allvar <- read.csv("data/allvar.csv")
cd4 <- allvar %>%
filter(treatmnt == 1, !is.na(CD4PCT), baseage > 1, baseage < 5) %>%
mutate(
y = sqrt(CD4PCT),
person = newpid,
time = visage - baseage
)
ggplot(cd4, aes(x = time, y = y, group = person)) +
geom_line(alpha = 0.5)
Veamos un ajuste usando la función lmer() del paquete lme4.
fit_cd4 <- lmer(formula = y ~ time + (1 + time | person), cd4)
fit_cd4
#> Linear mixed model fit by REML ['lmerMod']
#> Formula: y ~ time + (1 + time | person)
#> Data: cd4
#> REML criterion at convergence: 1096
#> Random effects:
#> Groups Name Std.Dev. Corr
#> person (Intercept) 1.329
#> time 0.680 0.15
#> Residual 0.748
#> Number of obs: 369, groups: person, 83
#> Fixed Effects:
#> (Intercept) time
#> 4.846 -0.468Notamos que las tendencias sobre el tiempo \(\beta\) tienen un promedio estimado en -0.5 con desviación estándar de 0.7, es decir, estimamos que la mayoría de los niños tienen niveles de CD4 decrecientes, pero no todos.
Usaremos estos resultados para hacer calculos de potencia para una nueva prueba que busca medir el efecto del consumo de zinc en la dieta. Quisiéramos que el estudio fuera suficientemente grande para que con probabilidad de al menos 80% la media del efecto del tratamiento sea significativo con un nivel de confianza del 95%.
Necesitamos hacer supuestos del efecto del tratamiento y del resto de los parámetros que caracterizan el estudio. El análisis de arriba muestra que en los niños VIH positivos que no recibieron zinc los niveles de CD4 caían en promedio 0.5 al año. Suponemos que con el zinc reduciremos la caída a cero.
$$ y_i N({j[i]} + {j[i]}t_i, ^2_y) \
\[\begin{eqnarray*} \begin{pmatrix}\alpha_{j}\\ \beta_{j} \end{pmatrix} & \sim & N\left[\left(\begin{array}{c} \gamma_0^{\alpha}\\ \gamma_0^{\beta}+\gamma_1^{\beta}z_j \end{array}\right), \left(\begin{array}{cc} \sigma^2_{\alpha} & \rho \sigma_{\alpha}\sigma_{\beta}\\ \rho \sigma_{\alpha}\sigma_{\beta} & \sigma^2_{\beta} \end{array}\right)\right] \end{eqnarray*}\]$$ donde
\[ z_j = \left\{ \begin{array}{lr} 1 & \text{si el }j \text{-ésimo niño recibió trataiento}\\ 0 & e.o.c \end{array} \right. \]
El tratamiento \(z_j\) afecta la pendiente \(\beta_j\) más no el intercepto \(\alpha_j\) pues el tratamiento no puede afectar en el tiempo cero. Usando los datos del ajuste de arriba tenemos que para el grupo control la pendiente será
\(\gamma_0^{\beta} = -0.5\) y el efecto del tratamiento \(\gamma_1^{\beta} = 0.5\), el resto de los parámetros los especificamos de acuerdo al ajuste de arriba. Por simplicidad fijaremos la correlación \(\rho\) en cero.
El siguiente paso es determinar el diseño del modelo, suponemos que dividiremos a \(J\) niños VIH positivos en dos grupos del mismo tamaño, \(J/2\) de ellos recibirán el cuidado usual y \(J/2\) recibirán suplementos de zinc. Más aún suponemos que se medirá el porcentaje de CD4 cada 2 meses durante un año.
Usaremos simulación para determinar el tamaño de muestra \(J\) que se requiere para tener una potencia de \(80%\) si el verdadero efecto es 0.5, ¿cuál es el modelo gráfico asociado?
cd4_sim <- function (J, K, mu.a.true = 4.8, g.0.true = -0.5, g.1.true = 0.5,
sigma.y.true = 0.7, sigma.a.true = 1.3, sigma.b.true = 0.7){
time <- rep(seq(0, 1, length = K), J) # K mediciones en el año
person <- rep (1:J, each = K) # ids
treatment <- sample(rep(0:1, J/2))
treatment1 <- treatment[person]
# parámetros a nivel persona
a.true <- rnorm(J, mu.a.true, sigma.a.true)
b.true <- rnorm(J, g.0.true + g.1.true * treatment, sigma.b.true)
y <- rnorm (J * K, a.true[person] + b.true[person] * time, sigma.y.true)
return (data.frame(y, time, person, treatment1))
}
# 2. Function to fit the model and calculate the power
cd4_power <- function (J, K){
fake <- cd4_sim(J, K)
lme_power <- lmer (y ~ time + time:treatment1 + (1 + time | person), data = fake)
theta_hat <- fixef(lme_power)["time:treatment1"]
theta_se <- summary(lme_power)$coefficients["time:treatment1", "Std. Error"]
theta_hat - 2 * theta_se > 0
}
cd4_rep <- function(n_sims, J, K){
rerun(n_sims, cd4_power(J, K = 7)) %>% flatten_dbl() %>% mean()
}
potencias <- map_df(c(4, 16, 60, 100, 150, 200, 225, 250, 300, 400),
~data_frame(n = ., p = cd4_rep(n_sims = 500, J = ., K = 7)))
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
#> $checkConv, : unable to evaluate scaled gradient
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
#> $checkConv, : Model failed to converge: degenerate Hessian with 1 negative
#> eigenvalues
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
#> $checkConv, : unable to evaluate scaled gradient
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
#> $checkConv, : Model failed to converge: degenerate Hessian with 1 negative
#> eigenvalues
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
#> $checkConv, : unable to evaluate scaled gradient
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
#> $checkConv, : Model failed to converge: degenerate Hessian with 1 negative
#> eigenvalues
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
#> $checkConv, : unable to evaluate scaled gradient
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
#> $checkConv, : Model failed to converge: degenerate Hessian with 1 negative
#> eigenvalues
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
#> $checkConv, : unable to evaluate scaled gradient
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
#> $checkConv, : Model failed to converge: degenerate Hessian with 1 negative
#> eigenvalues
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
#> $checkConv, : unable to evaluate scaled gradient
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
#> $checkConv, : Model failed to converge: degenerate Hessian with 1 negative
#> eigenvalues
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
#> $checkConv, : unable to evaluate scaled gradient
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
#> $checkConv, : Model failed to converge: degenerate Hessian with 1 negative
#> eigenvalues
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
#> $checkConv, : unable to evaluate scaled gradient
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
#> $checkConv, : Model failed to converge: degenerate Hessian with 1 negative
#> eigenvalues
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
#> $checkConv, : unable to evaluate scaled gradient
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
#> $checkConv, : Model failed to converge: degenerate Hessian with 1 negative
#> eigenvalues
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
#> $checkConv, : unable to evaluate scaled gradient
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
#> $checkConv, : Model failed to converge: degenerate Hessian with 1 negative
#> eigenvalues
ggplot(potencias, aes(x = n, y = p)) +
geom_hline(yintercept = 0.8, color = "red", alpha = 0.5) +
geom_line()
Notemos que la función cd4_rep() regresa la proporción de las simulaciones en las que el resultado es estadísticamente significarivo, est es, la potencia calculada con simulaicón, para un estudio con \(J\) niños medidos en \(K\) intervalos igualmente espaciados.
Notemos que en el límite, cuando \(J \to 0\) el poder es 0.025, esto es, con una muestra suficientemente chica el efecto del estimador es básicamente aleatorio y por tanto en 2.5% de los casos el estimador está 2 desviaciones por encima de cero.
Una ventaja de usar simulación para calcular potencia es que nos permite flexibilidad, por ejemplo, es fácil calcular para más escenarios:
¿qué ocurriría si solo puedo medir 3 veces al año?
Se sabe que es común que algunos participantes abandonen el estudio, o no asistan a todas las mediciones, con simulación es fácil incorporar faltantes.
potencias_2 <- map_df(c(4, 16, 60, 100, 150, 200, 225, 250, 300, 400),
~data_frame(n = ., p = cd4_rep(n_sims = 200, J = ., K = 3)))
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
#> $checkConv, : unable to evaluate scaled gradient
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
#> $checkConv, : Model failed to converge: degenerate Hessian with 1 negative
#> eigenvalues
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
#> $checkConv, : unable to evaluate scaled gradient
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
#> $checkConv, : Model failed to converge: degenerate Hessian with 1 negative
#> eigenvalues
ggplot(potencias, aes(x = n, y = p)) +
geom_hline(yintercept = 0.8, color = "red", alpha = 0.5) +
geom_line() +
geom_line(data = potencias_2)
¿Para qué simular?
Alternativa para presentar inferencias: en lugar de presentar un estimador puntual y/o intervalo de confianza podemos analizar simulaciones del modelo.
Inferencia predictiva: es fácil usar simulación para calcular errores estándar, o intervalos de confianza, resulta particularmente útil cuando estamos estimando cantidades que no son coeficientes de un modelo o transformaciones lineales de coeficientes.
Simulación para revisar el ajuste de un modelo. Podemos simular datos del modelo ajustado y compararlos con los datos verdaderos.
Referencias
Wickham, H., Cook, D., Hofmann, H. and Buja, A. (2010) Graphical Inference for Infovis
Hofmann, H., Follett, L., Majumder, M. and Cook, D. (2012) Graphical Tests for Power Comparison of Competing Designs.
Tufte, E. Making better inferences from statistical graphics.