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 una población de estudiantes que están tomando un determinado curso, el espacio de resultados es el conjunto de todos los estudiantes en la población. En muchas ocasiones buscamos resolver preguntas que involucran más de una variable aleatoria, consideremos una variable aleatoria X que mide la aptitud de los estudiantes y toma 2 valores: baja (0) ó alta (1) y una variable aleatoria Y relacionada con la calificación de los estudiante y toma tres valores: baja (0), media (1) ó alta (2). La distribución conjunta de las variables aleatorias discretas se puede representar por medio de tablas.

Y/X X = 0 X=1
Y=0 0.35 0.03
Y=1 0.28 0.09
Y=2 0.07 0.18

Entonces, \(P(X=0, Y=0) = 0.25\).

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:

  1. \(p(x,y) \geq 0\) para toda \((x,y)\).

  2. \(\int_{-\infty}^{\infty}p(x,y)dxdy=1\).

  3. 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.


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}\]


En nuestro ejemplo discreto la distribución marginal de Habilidad \(X\) se obtiene sumando sobre los tres posibles valores de Calificación \(Y\) y corresponde al márgen inferior de la siguiente tabla. Por su parte, la distribución marginal de Calificación corresponde al márgen derecho.

Y/X X = 0 X=1 .
Y=0 0.35 0.03 0.38
Y=1 0.28 0.09 0.37
Y=2 0.07 0.18 0.25
. 0.7 0.3 1

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 Habilidad dado Calificación = baja (\(Y = 0\))?

. X=0 X = 1
Y=0 0.92 0.08

Para obtener toda la distribución condicional calculamos los dos casos restantes : Y=1, Y=2:

. X=0 X = 1
Y=0 0.92 0.08
Y=1 0.76 0.24
Y=2 0.28 0.72

Vale la pena destacar que una distribución condicional es una distribución de probabilidad. En el ejemplo anterior notamos que para cada renglón de la tabla las probabilidades suman uno, son no negativas y menores que uno.

Análogo a la regla de probabilidad total que estudiamos para eventos, tenemos la siguiente definición para variables aleatorias:


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) \] \[p_X(x) = \int p_{X \vert Y}(x\vert y)p_Y(y)dy. \]


Verifiquemos la igualdad usando el ejemplo de habilidad y calificación.

Para variables aleatorias tenemos también una definición análoga para la regla de Bayes:

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 regla de la cadena.

Nótese que esta regla funciona para cualquier ordenamiento de las variables aleatorias.

Independencia

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

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.

Ejercicio. 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?

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\).

Ejercicio. Recordando el ejemplo de Habilidad y Calificación, compara \(p_{X\vert Y}(x\vert y)\) con \(p_X(x)\), para determinar si las variables aleatorias X y Y no son independientes.

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. Por ejemplo, supongamos que nos interesa la probabilidad de que un alumno sea admitido a Stanford o a MIT. Es razonable pensar que estos eventos no son independientes: si sabemos que un alumno fue admitido en Stanford nuestra estimación de la probabilidad de que sea admitido en MIT aumenta. Ahora, supongamos que ambas universidades basan su decisión de aceptación/rechazo únicamente en el promedio (promedio alto, medio o bajo) y sabemos que para el estudiante que analizamos promedio es alto. En este caso, podemos argumentar que conocer que el estudiante fue admitido en Stanford no cambia la probabilidad de que el estudiante sea admitido en MIT pues el evento promedio alto contiene toda la información relevante para conocer su oportunidad de ser admitido en MIT. En este caso, decimos que ser admitido en MIT es condicionalmente independiente de ser admitido en Stanford dado que el promedio es alto.

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\).

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).\]

De manera similar al caso de independencia, \(X\) es independiente de \(Y\) condicional a \(Z\) sí y sólo sí, \(p(x,y,z) \propto g(x,z)h(y,z)\). Adicionalmente, la definición de independencia condicional es equivalente a \(p_{X,Y \vert Z}(x, y\vert z) = p_{X\vert Z}(x\vert z)p_{Y\vert Z}(y\vert z)\). Esto es, una vez que conocemos \(Z\), si observamos \(Y\) ésta no provee información adicional sobre \(X\).

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:

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?

## 
## Attaching package: 'igraph'
## 
## The following objects are masked from 'package:stats':
## 
##     decompose, spectrum
## 
## The following object is masked from 'package:base':
## 
##     union

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(plyr)
library(dplyr)
library(arm)
library(tidyr)
set.seed(918739837)
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 <- rdply(1000, rbinom(1, 400, 0.488))
mean(sims_ninas$V1)
## [1] 195.773
ggplot(sims_ninas, aes(x = V1)) + geom_histogram(binwidth = 2)

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] 183

Repetimos 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_ninas2 <- rdply(1000, modelo2)
mean(sims_ninas2$V1)
## [1] 197.45
ggplot(sims_ninas2, aes(x = V1)) + geom_histogram(binwidth = 2)

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

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] 170.1001

Simulamos 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))
  mean(altura)
}
sims_alturas <- rdply(1000, mediaAltura)
mean(sims_alturas$V1)
## [1] 168.0834
ggplot(sims_alturas, aes(x = V1)) + geom_histogram(binwidth = 1)

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.

Ejemplo de simulación de un modelo de regresió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\]

vector_mu <- c(rep(78 + 12, 30), rep(78, 20)) # beta_0 + beta_1 X
y <- rnorm(50, vector_mu, 20)
sims_y <- rdply(2000, rnorm(50, vector_mu, 20))

Podemos calcular la media y su intervalo de confianza:

head(sims_y[, 1:5])
##   .n        V1        V2        V3        V4
## 1  1  82.47508 101.64550  52.78054 104.18743
## 2  2 121.44829  70.26634  96.90540  64.90629
## 3  3 111.06877 102.09530 102.44267  91.46075
## 4  4  63.64287  83.31220  86.10466  92.42929
## 5  5  59.66619 105.37822  87.54980  51.27433
## 6  6  83.74825  92.94726 108.71808 114.01212
medias <- apply(sims_y[, -1], 1, mean)
quantile(medias, c(0.025, 0.975))
##     2.5%    97.5% 
## 79.47659 90.89766
hist(medias)

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:

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.1 & -0.1\\ -0.1 & 0.1 \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.

  1. Simula \(\sigma=20\sqrt{(432)/X}\) donde \(X\) es una generación de una distribución \(\Chi^2\) con \(432\) grados de libertad.

  2. 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\).

  3. Simula \(y\) el vector de observaciones usando los parámetros de 1 y 2.

simsIter <- function(){
  # empezamos simulando sigma
  sigma <- 20 * sqrt((432) / rchisq(1, 432))
  # la usamos para simular betas
  beta <- mvrnorm(1, mu = c(78, 12), 
    Sigma = sigma ^ 2 * matrix(c(0.1, -0.1, -0.1, 0.1), nrow = 2))
  vector_mu <- c(rep(beta[1] + beta[2], 30), rep(beta[1], 20)) # beta_0 + beta_1 X
  # Simulamos las observaciones
  mean(rnorm(50, vector_mu, sigma))
}

sims_y2 <- rdply(2000, simsIter)

hist(sims_y2$V1)

Y podemos graficar la distribución

ggplot(sims_y2, aes(x = V1)) +
  geom_histogram(binwidth = 3)

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ócrata 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)
names(paths) <- basename(paths)

# Leemos los datos y recodificamos variables
congress <- tbl_df(ldply(paths, read.table, quote = "\"", 
  stringsAsFactors = FALSE)) %>%
  mutate(id = rep(1:435, 3),
         year = extract_numeric(.id),          # 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)

congress
## Source: local data frame [1,305 x 4]
## 
##    id year incumbency dem_share
## 1   1 1986          1 0.7450362
## 2   2 1986          1 0.6738455
## 3   3 1986          1 0.6964566
## 4   4 1986         -1 0.4645901
## 5   5 1986          1 0.3910945
## 6   6 1986         -1 0.3582454
## 7   7 1986          0 0.5485980
## 8   8 1986         -1 0.2267028
## 9   9 1986         -1 0.2218160
## 10 10 1986          1 0.6593966
## .. ..  ...        ...       ...
# 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(join(congress_share, congress_inc))
## Joining by: id
head(congress_h)
##   id  share_86  share_88  share_90 inc_86 inc_88 inc_90
## 1  1 0.7450362 0.7724427 0.7140294      1      1      1
## 2  2 0.6738455 0.6361816 0.5970501      1      1      1
## 3  3 0.6964566 0.6649283 0.5210433      1      1      0
## 4  4 0.4645901 0.2738342 0.2343770     -1     -1     -1
## 5  5 0.3910945 0.2636131 0.4774393      1     -1      0
## 6  6 0.3582454 0.3341927 0.2562970     -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)
display(fit_88)
## lm(formula = share_88 ~ share_86 + inc_88, data = congress_h)
##             coef.est coef.se
## (Intercept) 0.30     0.02   
## share_86    0.39     0.03   
## inc_88      0.10     0.01   
## ---
## n = 354, k = 3
## residual sd = 0.08, R-Squared = 0.85
summary(fit_88)
## 
## Call:
## lm(formula = share_88 ~ share_86 + inc_88, data = congress_h)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.241640 -0.042841 -0.005898  0.052343  0.226273 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.300413   0.015107   19.89   <2e-16 ***
## share_86    0.385833   0.027847   13.86   <2e-16 ***
## inc_88      0.104263   0.006904   15.10   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.07738 on 351 degrees of freedom
## Multiple R-squared:  0.8461, Adjusted R-squared:  0.8452 
## F-statistic: 964.7 on 2 and 351 DF,  p-value: < 2.2e-16

Simulemos del modelo

# Matriz X
X <- cbind(1, congress_h$share_88, congress_h$inc_90)
simsIter <- function(){
  sigma <- 0.08 * sqrt((351) / rchisq(1, 351))
  beta <- 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
  rnorm(358, mu, sigma)
}

sims_congress <- rdply(2000, simsIter)
sims_long <- gather(sims_congress, district, share, V1:V358)

ggplot(sims_long, aes(x = reorder(district, share), y = share)) +
  geom_hline(color = "red", yintercept = 0.5) +
  geom_boxplot()

Podemos preguntarnos cuántas elecciones ganaron los demócratas en 1990: \(\sum I(\grave{y} > 0.5)\)

sims_long %>%
  group_by(.n) %>%
  mutate(wins = sum(share > 0.5)) %>%
  ungroup() %>%
  summarise(mean_wins = mean(wins),
            median_wins = median(wins), 
            sd_wins = sd(wins))
## Source: local data frame [1 x 3]
## 
##   mean_wins median_wins  sd_wins
## 1   201.977         202 3.893519

Veamos lo que ocurrió realmente

sum(congress_h$share_90 > 0.5)
## [1] 207