1 Introducción

Algunos datos de algunas personas están perdidos (missing).

Si algunos datos están perdidos no es claro cómo hacer la estimación de los parámetros del modelo.

Es una mala práctica ajustar el modelo:

Los datos faltantes introducen sesgos en un análisis de los datos completos dependiendo de la estructura causal del proceso de perdida de datos.

Perdidos completamente al azar (MCAR, missing completely at random)

La falta de datos es independiente de los datos observados y no observados.

No existen diferencias sistemáticas entre los participantes con datos faltantes y aquellos con datos completos.

Los datos completos se pueden considerar como una muestra aleatoria de la población de interés.

Perdidos al azar (MCAR, missing at random)

La falta de datos está relacionada con los datos observados, pero no con los no observados.

Por ejemplo, en un estudio de depresión, la probabilidad de completar la encuesta puede estar relacionada con el sexo (variable de control completamente observada), pero no con la gravedad de la depresión (variable de interés).

Un manejo adecuado de los factores conocidos puede conducir a resultados imparciales.

Perdidos no al azar (MCAR, missing not at random)

La falta de datos está relacionada con los datos no observados (i.e., con eventos o factores que no son medidos por el investigador).

Continuando con el ejemplo anterior, los participantes con depresión severa tienen más probabilidades de negarse a completar la encuesta sobre la gravedad de la depresión (variable de interés).

Por lo general no es posible llevar a cabo resultados imparciales dado que no se observan factores asociados con la perdida de datos.

2 Modelo Normal Multivariado Semi-Conjugado

La distribución posterior \(p(\boldsymbol{\theta},\mathbf{\Sigma}\mid\mathbf{Y})\) depende directamente de la distribución muestral \(p(\mathbf{Y}\mid\boldsymbol{\theta},\mathbf{\Sigma}) = \prod_{i=1}^n p(\boldsymbol{y}_i\mid\boldsymbol{\theta},\mathbf{\Sigma})\), con \(\mathbf{Y}=[\boldsymbol{y}_1^\textsf{T},\ldots,\boldsymbol{y}_n^\textsf{T}]^\textsf{T}\), la cual no se puede calcular si faltan observaciones en algunas de las \(\boldsymbol{y}_i\).

3 Estimación con datos faltantes al azar

Sea \(\boldsymbol{o}_i = (o_{i,1},\ldots,o_{i,p})\) un vector con entradas binarias tales que \[ o_{i,j} = \begin{cases} 1, & \text{si $y_{i,j}$ es un dato observado;} \\ 0, & \text{si $y_{i,j}$ es un dato faltante,} \end{cases} \] para \(i=1,\ldots,n\) y \(j = 1,\ldots,p\). Además, \(\mathbf{O}= [\boldsymbol{o}_1^\textsf{T},\ldots,\boldsymbol{o}_n^\textsf{T}]^\textsf{T}\).

La matriz \(\mathbf{Y}\) consiste de dos partes, \(\mathbf{Y}= (\mathbf{Y}_{\text{obs}},\mathbf{Y}_{\text{miss}})\), donde \[ \mathbf{Y}_{\text{obs}} = \{y_{i,j}: o_{i,j} = 1\} \qquad\text{y}\qquad \mathbf{Y}_{\text{miss}} = \{y_{i,j}: o_{i,j} = 0\}\,. \]

Similarmente, \(\boldsymbol{y}_i = (\boldsymbol{y}_{i,\text{obs}},\boldsymbol{y}_{i,\text{miss}})\), para \(i=1,\ldots,n\).

La información observada del individuo \(i\) es \(\boldsymbol{o}_i\) y \(\boldsymbol{y}_{i,\text{obs}}\).

Bajo el supuesto de datos faltantes al azar:

3.1 Distribución posterior

A partir de la información observada, se quiere obtener la distribución posterior de las cantidades desconocidas y faltantes: \[ \begin{align*} p(\boldsymbol{\theta},\mathbf{\Sigma},\mathbf{Y}_{\text{miss}}\mid\mathbf{Y}_{\text{obs}},\mathbf{O}) &\propto p(\mathbf{Y}_{\text{obs}},\mathbf{O}\mid \boldsymbol{\theta},\mathbf{\Sigma},\mathbf{Y}_{\text{miss}})\,\, p(\boldsymbol{\theta},\mathbf{\Sigma},\mathbf{Y}_{\text{miss}})\\ &= p(\mathbf{O}\mid \boldsymbol{\theta},\mathbf{\Sigma},\mathbf{Y}_{\text{miss}},\mathbf{Y}_{\text{obs}})\,\, p(\mathbf{Y}_{\text{obs}}\mid \boldsymbol{\theta},\mathbf{\Sigma},\mathbf{Y}_{\text{miss}})\,\, p(\mathbf{Y}_{\text{miss}}\mid\boldsymbol{\theta},\mathbf{\Sigma})\,\, p(\boldsymbol{\theta},\mathbf{\Sigma}) \\ &= p(\mathbf{O}\mid \mathbf{Y}_{\text{obs}})\,\, p(\mathbf{Y}_{\text{obs}}\mid \mathbf{Y}_{\text{miss}}, \boldsymbol{\theta},\mathbf{\Sigma})\,\, p(\mathbf{Y}_{\text{miss}}\mid\boldsymbol{\theta},\mathbf{\Sigma})\,\, p(\boldsymbol{\theta})\,\, p(\mathbf{\Sigma}) \\ &\propto p(\mathbf{Y}_{\text{obs}}\mid \mathbf{Y}_{\text{miss}}, \boldsymbol{\theta},\mathbf{\Sigma})\,\, p(\mathbf{Y}_{\text{miss}}\mid\boldsymbol{\theta},\mathbf{\Sigma})\,\, p(\boldsymbol{\theta})\,\,p(\mathbf{\Sigma})\,. \end{align*} \]

Para obtener las distribuciones condicionales completas, se observa que \[ p(\mathbf{Y}\mid\boldsymbol{\theta},\mathbf{\Sigma}) = p(\mathbf{Y}_{\text{obs}},\mathbf{Y}_{\text{miss}}\mid\boldsymbol{\theta},\mathbf{\Sigma}) = p(\mathbf{Y}_{\text{obs}}\mid \mathbf{Y}_{\text{miss}}, \boldsymbol{\theta},\mathbf{\Sigma})\,\, p(\mathbf{Y}_{\text{miss}}\mid\boldsymbol{\theta},\mathbf{\Sigma}) \] y \[ \begin{align*} p(\mathbf{Y}_{\text{miss}}\mid \mathbf{Y}_{\text{obs}}, \boldsymbol{\theta},\mathbf{\Sigma}) &\propto p(\mathbf{Y}_{\text{obs}},\mathbf{Y}_{\text{miss}}\mid\boldsymbol{\theta},\mathbf{\Sigma}) \\ &= \prod_{i=1}^n p (\boldsymbol{y}_{i,\text{obs}},\boldsymbol{y}_{i,\text{miss}}\mid \boldsymbol{\theta},\mathbf{\Sigma}) \\ &\propto \prod_{i=1}^n p (\boldsymbol{y}_{i,\text{miss}}\mid \boldsymbol{y}_{i,\text{obs}}, \boldsymbol{\theta},\mathbf{\Sigma})\,. \end{align*} \]

3.2 Muestreador de Gibbs

A partir del estado actual \((\boldsymbol{\theta}^{(b)},\mathbf{\Sigma}^{(b)},\mathbf{Y}_{\text{miss}}^{(b)})\), se genera un nuevo estado \((\boldsymbol{\theta}^{(b+1)},\mathbf{\Sigma}^{(b+1)},\mathbf{Y}_{\text{miss}}^{(b+1)})\) como sigue:

  1. Simular \(\mathbf{Y}_{\text{miss}}^{(b+1)} \sim p(\mathbf{Y}_{\text{miss}}\mid \boldsymbol{\theta}^{(b)}, \mathbf{\Sigma}_{\text{miss}}^{(b)}, \mathbf{Y}_{\text{obs}})\).
  2. Simular \(\boldsymbol{\theta}^{(b+1)} \sim p(\boldsymbol{\theta}\mid \mathbf{\Sigma}^{(b)}, \mathbf{Y}_{\text{miss}}^{(b+1)}, \mathbf{Y}_{\text{obs}})\).
  3. Simular \(\mathbf{\Sigma}^{(b+1)} \sim p(\mathbf{\Sigma}\mid \boldsymbol{\theta}^{(b+1)}, \mathbf{Y}_{\text{miss}}^{(b+1)}, \mathbf{Y}_{\text{obs}})\).

3.3 Distribuciones condicionales completas

  • La distribución condicional de \(\boldsymbol{y}_{i,[b]}\mid\boldsymbol{y}_{i,[a]},\boldsymbol{\theta},\mathbf{\Sigma}\) es \(\boldsymbol{y}_{i,[b]}\mid\text{resto} \sim \textsf{N}(\boldsymbol{\theta}_{b\mid a},\mathbf{\Sigma}_{b\mid a})\) donde \[ \boldsymbol{\theta}_{b\mid a} = \boldsymbol{\theta}_{[b]} - \mathbf{\Sigma}_{[b,a]}\mathbf{\Sigma}_{[a,a]}^{-1}(\boldsymbol{y}_{i,[a]} - \boldsymbol{\theta}_{[a]}) \qquad\text{y}\qquad \mathbf{\Sigma}_{b\mid a} = \mathbf{\Sigma}_{[b,b]} - \mathbf{\Sigma}_{[b,a]}\mathbf{\Sigma}_{[a,a]}^{-1}\mathbf{\Sigma}_{[a,b]}\, \] donde \(b\) y \(a\) almacenan los subíndices de \(\boldsymbol{y}_i\) que son \(\text{miss}\) y \(\text{obs}\), respectivamente.

  • La distribución condicional completa de \(\boldsymbol{\theta}\) es \(\boldsymbol{\theta}\mid\text{resto}\sim\textsf{N}(\boldsymbol{\mu}_n,\mathbf{\Lambda}_n)\), donde \[ \boldsymbol{\mu}_n = \left( \mathbf{\Lambda}_0^{-1} + n\mathbf{\Sigma}^{-1} \right)^{-1} \left( \mathbf{\Lambda}_0^{-1}\boldsymbol{\mu}_0 + n\mathbf{\Sigma}^{-1}\bar{\boldsymbol{y}} \right) \qquad\text{y}\qquad \mathbf{\Lambda}_n = \left( \mathbf{\Lambda}_0^{-1} + n\mathbf{\Sigma}^{-1} \right)^{-1}\,. \]

  • La distribución condicional completa de \(\mathbf{\Sigma}\) es \(\mathbf{\Sigma}\mid\text{resto}\sim\textsf{WI}(\nu_n,\mathbf{S}_n^{-1})\), donde \[ \nu_n = \nu_0 + n \qquad\text{y}\qquad \mathbf{S}_n = \mathbf{S}_0 + \sum_{i=1}^n (\boldsymbol{y}_i - \boldsymbol{\theta})(\boldsymbol{y}_i - \boldsymbol{\theta})^{\textsf{T}} = \mathbf{S}_0 + (n-1)\mathbf{S}+ n(\bar{\boldsymbol y} - \boldsymbol{\theta})(\bar{\boldsymbol y} - \boldsymbol{\theta})^{\textsf{T}} \,. \]

4 Ejemplo: Diabetes

A una población de mujeres que tenían al menos 21 años, de herencia India Pima y que vivían cerca de Phoenix, Arizona, se les realizó una prueba de diabetes según los criterios de la Organización Mundial de la Salud.

Los datos fueron recopilados por el Instituto Nacional de Diabetes y Enfermedades Digestivas y Renales de EE. UU.

Datos disponibles en el paquete MASS de R.

4.1 Datos

# datos
# Y0 -> base de datos completa
# Y  -> base de datos con datos perdidos al azar 
# O  -> indicadora de datos observados
suppressMessages(suppressWarnings(library(MASS))) 
data(Pima.tr)
Y0 <- Pima.tr[,2:5]
(n <- dim(Y0)[1])
## [1] 200
(p <- dim(Y0)[2])
## [1] 4
# datos perdidos al azar
Y <- Y0
set.seed(1)
O <- matrix(data = rbinom(n = n*p, size = 1, prob = 0.9), nrow = n, ncol = p)
Y[O == 0] <- NA
# encabezado
head(Y)
##   glu bp skin  bmi
## 1  86 68   28 30.2
## 2 195 70   33   NA
## 3  77 82   NA 35.8
## 4  NA 76   43 47.9
## 5 107 60   NA   NA
## 6  97 76   27   NA
# vector de medias muestral
round(colMeans(Y, na.rm = T), 2)
##    glu     bp   skin    bmi 
## 123.58  70.84  29.10  32.21
round(var(Y0, na.rm = T), 2)
##          glu     bp   skin   bmi
## glu  1002.81  97.93  80.79 42.08
## bp     97.93 131.78  35.66 16.81
## skin   80.79  35.66 137.47 47.37
## bmi    42.08  16.81  47.37 37.58
# visualización
par(mfrow = c(p,p), mar = 1.75*c(1,1,0.5,0.5), mgp = c(1.75,0.75,0))
for(j1 in 1:p) {
  for(j2 in 1:p) {
    if (j1 == j2) {
      hist(Y[,j1], col = "gray90", border = "gray90", main = "")
      mtext(colnames(Y)[j1], side=3, line = -0.1, cex = 0.7)
    }
    if (j1 != j2) {
      plot(x = Y[,j1], y = Y[,j2], xlab = "", ylab = "", pch = 16, cex = 0.7)
    }
  }
}

4.2 Ajuste del modelo

# hiperparámetros
# alternativa: usar información auxiliar de la encuesta nacional de salud
mu0 <- colMeans(Y, na.rm = T)
L0  <- var(Y, na.rm = T)
nu0 <- p + 2
S0  <- L0
# valores iniciales
theta <- mu0 
Sigma <- S0
Yfull <- Y
for(j in 1:p)
  Yfull[is.na(Yfull[,j]), j] <- mean(Yfull[,j], na.rm = T)
# número de muestras
B <- 10000
ncat <- floor(B/10)
# almacenamiento
THETA <- SIGMA <- YMISS <- LP <- NULL
# cadena
iL0 <- solve(L0)
Lm0 <- solve(L0)%*%mu0
nun <- nu0 + n
for (b in 1:B) {
  # actualizar Y_miss
  for(i in 1:n) {
    miss <- O[i,] == 0
    if (any(miss)) {
      obs <- O[i,] == 1
      iSa <- solve(Sigma[obs,obs])
      m_j <- theta[miss] + Sigma[miss,obs]%*%iSa%*%t(Yfull[i,obs] - theta[obs])
      V_j <- Sigma[miss,miss] - Sigma[miss,obs]%*%iSa%*%Sigma[obs,miss]
      Yfull[i,miss] <- c(mvtnorm::rmvnorm(n = 1, mean = m_j, sigma = V_j))
    }
  }
  # actualizar estadísticos suficientes
  yb  <- colMeans(Yfull)
  SS  <- cov(Yfull)
  SSn <- S0 + (n-1)*SS
  # actualizar theta
  iSigma <- solve(Sigma)
  Ln     <- solve(iL0 + n*iSigma)
  theta  <- c(mvtnorm::rmvnorm(n = 1, mean = Ln%*%(Lm0 + n*(iSigma%*%yb)), sigma = Ln))
  # actualizar Sigma
  Sigma <- solve(rWishart::rWishart(n = 1, df = nun, Sigma = solve(SSn + n*((yb - theta)%*%t(yb - theta))))[,,1])     
  # log-verosimilitud
  LP[b] <- sum(apply(X = Yfull, MARGIN = 1, FUN = function(y) mvtnorm::dmvnorm(x = y, mean = theta, sigma = Sigma, log = T)))
  # almacenar
  THETA <- rbind(THETA, theta)
  SIGMA <- rbind(SIGMA, c(Sigma))
  YMISS <- rbind(YMISS, Yfull[O == 0])
  # progreso
  if (b%%ncat == 0) 
       cat(100*round(b/B, 1), "% completado ... \n", sep = "")
}

4.3 Convergencia

Cadena de la log-verosimilitud:

Coeficiente de Variación (%) de Monte Carlo:

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.013   0.114   0.244   0.249   0.421   0.494

4.4 Inferencia

Inferencia sobre el vector de medias:

colnames(THETA) <- paste0("theta", 1:p)
round(t(apply(X = THETA, MARGIN = 2, FUN = quantile, probs = c(0.025,0.5,0.975))), 2)
##          2.5%    50%  97.5%
## theta1 119.02 123.45 127.85
## theta2  69.49  71.06  72.68
## theta3  27.64  29.36  31.13
## theta4  31.30  32.17  33.03

Inferencia sobre la matriz de correlación:

COR <- array(data = NA, dim = c(p,p,B))
colnames(COR) <- rownames(COR) <- colnames(Y)
for(b in 1:B) {
  Sig <- matrix(data = SIGMA[b,], nrow = p, ncol = p)
  COR[,,b] <- Sig/sqrt(outer(diag(Sig), diag(Sig)))
}
round(apply(X = COR, MARGIN = c(1,2), FUN = mean), 2)
##       glu   bp skin  bmi
## glu  1.00 0.23 0.25 0.19
## bp   0.23 1.00 0.25 0.24
## skin 0.25 0.25 1.00 0.66
## bmi  0.19 0.24 0.66 1.00

Inferencia sobre las correlaciones \(\frac{\sigma_{i,j}}{\sigma_i\sigma_j}\) (izquierda) y los coeficientes de regresión de las predicciones \(\boldsymbol{\beta}_{b\mid a}^\textsf{T} = \mathbf{\Sigma}_{[b,a]}\mathbf{\Sigma}_{[a,a]}^{-1}\) (derecha):

4.5 Predicción

Referencias