09-Variables latentes y algoritmo EM
Teresa Ortiz
Preparamos paquetes.
library(plyr)
library(dplyr)
library(tidyr)Recordatorio de datos faltantes
El modelo para los datos completos está parametrizado con \(\theta\): \[p_{\theta}(y)\]
El proceso de censura lo incluimos mediante el vector \(I\) de indicadores de faltantes (\(I_{j}=1\) cuando el dato \(j\) está censurado faltante), entonces si el modelo para el mecanismo de faltantes está parametrizado con \(\psi\): \[p_{\psi} (I|y)\]
Entonces generamos los datos según (1) y luego censuramos observaciones según (2). Así, el modelo completo para nuestras observaciones es: \[p_{\theta,\psi} (y,I)=p_{\theta}(y)p_{\psi}(I|y).\]
Escribimos \(y=(y_{obs},y_{falta})\) para denotar por separado datos observados y datos faltantes. De tal manera que la verosimilitud para nuestros datos observados está dada por \[p_{\theta,\psi}(y_{obs},I),\] pues sabemos los valores de las variables observadas y qué datos faltan. Para hacer máxima verosimilitud calculamos esta conjunta. \[p_{\theta,\psi} (y_{obs},I)=\int p_{\theta,\psi} (y_{obs}, y_{falta},I)d y_{falta}\] \[=\int p_{\psi}(I|y_{obs},y_{falta})p_{\theta} (y_{obs}, y_{falta})d y_{falta}.\] Nótese que esta integral (o suma, dependiendo del caso), promedia los valores faltantes según el modelo de los datos \(p_{\theta}\) De la ecuación de arriba, vemos que en general tendríamos que modelar también el mecanismo de datos faltantes y estimar todos los parámetros. Esto es difícil no solamente porque hay más parámetros, sino porque en la práctica puede ser difícil proponer un modelo razonable para datos faltantes. Preferimos hacer un supuesto (MCAR, MAR).
Si se cumple MAR, entonces tenemos que \[p_{\psi}(I|y_{obs},y_{falta})=p_{\psi}(I|y_{obs}),\] y por tanto, \[p_{\theta,\psi} (y_{obs},I)= p_{\psi}(I|y_{obs})\int p_{\theta} (y_{obs}, y_{falta})dy_{falta}.\] notamos que los parámetros \(\psi\) y \(\theta\) están en factores separados y para encontrar los estimadores de máxima verosimilitud, no es necesario trabajar con \(\psi\) ni el mecanismo al azar. El mecanismo de faltantes es entonces ignorable.
Algoritmo de Esperanza-Maximización (EM)
Recordemos que el problema de maximización de verosimilitud se vuelve más difícil en la presencia de datos faltantes.
Ejemplo. Escogemos al azar una moneda y luego tiramos un volado con esa moneda (modelo Moneda -> Sol):
ej_1 <- data.frame(moneda = c('A', 'A', 'B', 'B', 'B', NA),
sol = c(1, 0, 0, 0, 1, 0))
ej_1## moneda sol
## 1 A 1
## 2 A 0
## 3 B 0
## 4 B 0
## 5 B 1
## 6 <NA> 0Si parametrizamos con \(\theta=P(moneda=A)\), \(\theta_A=P(sol|moneda=A)\) y \(\theta_B=P(sol|moneda=B)\), la log verosimilitud para los datos completos es \[\log (\theta\theta_A) + \log(\theta(1-\theta_A)) + 2\log((1-\theta)(1-\theta_B))+ \log((1-\theta)\theta_B)\] \[= 2\log\theta+3\log(1-\theta)+\log\theta_A+\log(1-\theta_A)+2\log(1-\theta_B)+\log\theta_B\] En este caso es fácil encontrar los estimadores de máxima verosimilitud. Ahora, para incluir la contribución del caso con faltante suponemos MAR y promediamos los valores faltantes (como indica la fórmula en (5), cambiando la integral por suma): \[p_{\theta}(x^{6}_{sol}=0 )=p_{\theta}(x^{6}_{sol}=0 |x^{6}_{moneda}=A) p_{\theta}(x^{6}_{moneda}=A) + p_{\theta}(x^{6}_{sol}=0 |x^{6}_{moneda}=B) p_{\theta}(x^{6}_{moneda}=B),\] y ahora buscamos maximizar \[p_{\theta}(y_{obs})=2\log\theta+3\log(1-\theta)+\log\theta_A+\log(1-\theta_A)+2\log(1-\theta_B)+\log\theta_B + \log((1-\theta_A)\theta + (1-\theta_B)(1-\theta)).\]
Notemos que esta expresión es más difícil de maximizar. El algoritmo EM da un algoritmo iterativo para encontrar máximos locales de la verosimilitud \(p_{\theta} (y_{obs})\). Este algoritmo sigue una estrategia de optimización motivada por la noción misma de datos faltantes y considerando la distribución condicional de los faltantes dado lo observado.
El algoritmo de esperanza-maximización (EM) es un proceso iterativo para maximizar verosimilitud en presencia de datos faltantes. Cada iteración consiste de dos pasos:
Calcular el valor esperado de la log-verosimilitud promediando sobre datos faltantes con una aproximación de la solución (\(\theta^{(t)}\)).
- Obtener una nueva aproximación maximizando la cantidad resultante en el paso previo.
Podemos ver que el algorimto EM que el paso de esperanza consiste en un soft assignment de cada valor faltante de acuerdo al modelo de los datos y los datos observados.
Denotemos por \(\theta^{(t)}\) el maximizador estimado en la iteración \(t\) (\(t=0,1,...\)). En lugar de tratar directamente con \(\log p_{\theta}(y_{obs})=\log p(y_{obs}\theta)\) el algoritmo EM trabaja sobre \(Q(\theta \vert \theta^{(t)})\) que definimos como la esperanza de la log-verosimilitud de los datos completos \(y=(y_{obs},y_{falta})\) condicional a los datos observados \(y_{obs}\):
\[Q(\theta \vert \theta^{(t)})=E\big[\log\mathcal{L}(\theta|y) \big|y_{obs}, \theta^{(t)}\big]\]
\[=E \big[\log p(y|\theta) \big |y_{obs}, \theta^{(t)}\big]\]
\[=\int_{y_{falta}} [\log p(y|\theta)] p(y_{falta}|y_{obs}, \theta^{(t)})dy_{falta}.\]
Es así que en el paso esperanza calculamos \[Q(\theta|\theta^{(t)})\] y en el paso maximización: \[\theta^{(t+1)} = argmax_{\theta}Q(\theta|\theta^{(t)})\]
Ejemplo. En nuestro ejemplo anterior, haríamos lo siguiente. Primero definimos la función que calcula el esperado de la log verosimilitud:
# logLik: verosimilitude de datos completos
logLik <- function(theta){
2 * log(theta[1]) + 3 * log(1 - theta[1]) + log(theta[2]) + log(1 - theta[2]) +
log(theta[3]) + 2 * log(1 - theta[3])
}
pasoEsperanza <- function(theta_ant){
function(theta){
# a = p(moneda=A|sol=1, theta_ant)
# b = p(moneda=B|sol=1, theta_ant)
a_0 <- theta_ant[1] * (1 - theta_ant[2])
b_0 <- (1 - theta_ant[1]) * (1 - theta_ant[3])
a <- a_0 / (a_0 + b_0)
b <- b_0 / (a_0 + b_0)
logLik(theta) + log(theta[1] * (1 - theta[2])) * a +
log((1 - theta[1]) * (1 - theta[3])) * b
}
}Ahora escogemos una solución inicial y calculamos el esperado
theta_0 <- c(0.1, 0.1, 0.1)
esperanza_0 <- pasoEsperanza(theta_0)Optimizamos este valor esperado:
max_1 <- optim(c(0.5, 0.5, 0.5), esperanza_0, control = list(fnscale = -1))
theta_1 <- max_1$par
theta_1## [1] 0.35 0.48 0.26esperanza_1 <- pasoEsperanza(theta_1)
max_2 <- optim(c(0.5, 0.5, 0.5), esperanza_1,control=list(fnscale=-1) )
theta_2 <- max_2$par
theta_2## [1] 0.38 0.44 0.27esperanza_2 <- pasoEsperanza(theta_2)
max_3 <- optim(c(0.5, 0.5, 0.5), esperanza_2,control=list(fnscale=-1) )
theta_3 <- max_3$par
theta_3## [1] 0.39 0.43 0.27Y vemos que recuperamos la misma solución que en el caso de maximizar usando la función optim.
Veamos porque funciona trabajar con \(Q(\theta|\theta^{(t)})\). Si escribimos: \[ p(y_{obs}|\theta)=\frac{p(y|\theta)}{p(y_{falta}|y_{obs},\theta)}.\]
Tomando logaritmos, \[\log{p}(y_{obs}|\theta)=\log{p}(y|\theta)-log{p}(y_{falta}|y_{obs}, \theta).\]
Y tomando el valor esperado con respecto la distribuión de \(y_{falta}|y_{obs}\) con parámetro \(\theta^{(t)}\) (que consideramos como una aproximación de los verdaderos parámetros), obtenemos:
\[E\big[\log p(y_{obs}|\theta)\big | y_{obs}, \theta^{(t)}\big]= E\big[\log p(y|\theta)\big | y_{obs}, \theta^{(t)}\big] - E\big[\log p(y_{falta}|\theta)\big | y_{obs}, \theta^{(t)}\big]\]
\[=\int_{y_{falta}}[\log p(y|\theta)] p(y_{falta}|y_{obs}, \theta^{(t)})dy_{falta} - \int_{y_{falta}}[\log{p}(y_{falta}| y_{obs},\theta)]p(y_{falta}|y_{obs},\theta^{(t)})dy_{falta}\] \[= Q(\theta|\theta^{(t)}) - H(\theta|\theta^{(t)})\]
La igualdad anterior es cierta para cualquier valor de \(\theta\), en particular, \(\theta=\theta^{(t)}\) por lo que: \[\log{p}(y_{obs}|\theta) -\log{p}(y_{obs}|\theta^{(t)})= Q(\theta|\theta^{(t)}) - Q(\theta^{(t)}|\theta^{(t)}) - [H(\theta|\theta^{(t)}) - H(\theta^{(t)}|\theta^{(t)})]\] utilizando la desigualdad de Jensen se puede demostrar que \(H(\theta|\theta^{(t)}) \leq H(\theta^{(t)}|\theta^{(t)})\), por lo que (en términos de la verosimilutd) \[\log\mathcal{L}(\theta|y_{obs}) -\log\mathcal{L}(\theta^{(t)}|t_{obs}) \geq Q(\theta|\theta^{(t)}) - Q(\theta^{(t)}|\theta^{(t)})\] y vemos que si incrementamos \(Q(\theta|\theta^{(t)})\) entonces también incrementamos \(\mathcal{L}(\theta|y_{obs})\).
Observaciones del algoritmo EM
El algoritmo EM es una generalización natural de la estiamación por máxima verosimilitud cuando hay presencia de datos faltantes.
El problema de maximización que ataca EM es más complejo que el problema de máxima verosimilitud con datos completos. En el caso usual \(\log p(x|\theta)\) tiene un único máximo global que muchas veces se puede encontrar de forma cerrada, por su parte la verosimilitud con datos faltantes tiene muchos máximos locales y no tiene solución cerrada.
El algoritmo EM reduce el problema de múltiples máximos a una secuencia de subproblemas (\(Q(\theta|\theta^{(t)}\)) cada uno con un único máximo global.
Estos subproblemas son tales que se garantiza convergencia de las soluciones \(\theta^{(1)}, \theta^{(2)}, ...\) a un máximo, pues a cada paso se incrementa monótonamente el valor de la log-verosimilitud de los datos observados.El algoritmo EM garantiza convergencia únicamente a un máximo local pues la secuencia de \(\theta^{(t)}\) depende del valor con el que se inicializa el proceso \(\theta^{(1)}\). Por tanto, es conveniente repetir el algoritmo con varios valores iniciales.
Una manera de entender el algoritmo EM es pensar que en cada paso de esperanza imputas los faltantes, sin embargo, en lugar de una imputación dura se realiza un soft assignment de cada valor faltante. El soft assignment se calcula obteniendo la probabilidad de cada alternativa para el valor faltante (usando el parámetro \(\theta^{(t)})\)) para crear una base de datos ponderada que considera todos los posibles escenarios de los faltantes. Al usar los valores esperados (o datos ponderados) en lugar de realizar una imputación dura el algoritmo toma en cuenta la confianza del modelo cada vez que completa los datos.
Variables latentes
Un caso importante de datos faltantes es cuando una variable está totalmente censurada. Esto puede suceder por dos razones:
Alguna variable claramente importante está totalmente censurada (por ejemplo, peso en un estudio de consumo de calorías).
Cuando buscamos añadir estructura a nuestro modelo para simplificar su estimación, interpretación o forma. Por ejemplo: hacer grupos de actitudes ante la comida y el ejercicio para explicar con una sola variable el consumo de comida chatarra (por ejemplo, análisis de factores).
En estos casos, estas variables se llaman variables latentes, pues consideramos que tienen un efecto que observamos a través de otras variables, pero no podemos observar directamente los valores de esta variable.
¿Cuál es el supuesto apropiado acerca de este tipo de valores censurados (variables latentes)?
- MCAR
- MAR
- MNAR
- Ninguno de estos
La siguiente tabla es una clasificación de los modelos de variable latente de acuerdo a la métrica de las variables latentes y observadas.
| Latentes/Observadas | Métricas | Categóricas |
|---|---|---|
| Métricas | Análisis de factores (FA) | Modelos de rasgos latentes (LTM) |
| Categóricas | Modelos de perfiles latentes (LPM) | Modelos de clases latentes (LCM) |
Modelos de perfiles latentes: Mezcla de normales
El ejemplo más clásico de variables latentes es el de mezcla de normales.
Ejemplo. Modelo de mezcla de dos normales. Consideremos los siguientes datos:
library(ggplot2)
set.seed(280572)
N <- 800
n <- sum(rbinom(N, 1, 0.6))
x_1 <- rnorm(N - n, 0, 0.5)
x_2 <- rnorm(n, 2, 0.3)
qplot(c(x_1, x_2))## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.
Estos datos tienen una estructura bimodal. Es poco apropiado modelar estos datos con un modelo normal \((\mu,\sigma^2)\).
Podemos entonces modelar pensando que los datos vienen de dos clases, cada una con una distribución normal pero con distintos parámetros. ¿Cómo ajustaríamos tal modelo?
Ahora, si vemos la mezcla Gaussiana desde la representación generativa, o formulación en variable latente, tenemos el modelo gráfico \(\Delta\) -> \(X\) donde \(\Delta\) es una indicadora de clase. En el caso del modelo de dos clases tenemos \(\delta \in \{0,1\}\) y sea \(P(\delta=1)=\pi\), escribimos la conjunta \[p(\delta, x)=\pi^{\delta}(1-\pi)^{1-\delta}[\delta\phi_{\theta_1}(x)+(1-\delta)\phi_{\theta_2}(x)]\]
y podemos verificar que la distribución marginal es una mezcla gaussiana: \[p(x)\sum_{\delta}p(x|\delta)p(\delta)\] \[=\phi_{\theta_1}(x) \pi + \phi_{\theta_2}(x)(1-\pi)\]
Ahora, si conocieramos la clase a la que pertenece cada observación (\(\delta^i\)) podríamos escribir la log-verosimilitud completa (sin censura) como \[\sum_{i=1}^N \log(\delta^i \phi_{\theta_1} (x^i)+ (1-\delta^i)\phi_{\theta_2}(x^i)) + \delta^i \log\pi + (1-\delta^i)\log(1-\pi).\]
Aquí, es fácil ver que la verosimilitud se separa en dos partes, una para \(\delta^i=1\) y otra para \(\delta^i=0\), y los estimadores de máxima verosimilitud son entonces:
\[\hat{\mu}_1=\frac{\sum_i\delta^i x^i}{\sum_i (\delta^i)}\] \[\hat{\mu}_2=\frac{\sum_i(1-\delta^i) x^i}{\sum_i (1-\delta^i)}\]
\[\hat{\sigma}_1^2=\frac{\sum_i\delta^i (x^i-\mu_1)^2}{\sum_i (\delta^i)}\] \[\hat{\sigma}_2^2=\frac{\sum_i(1-\delta^i) (x^i-\mu_2)^2}{\sum_i (1-\delta^i)},\]
y \(\hat{\pi}\) es la proporción de casos tipo 1 en los datos. Este problema es entonces trivial de resolver.
En el caso de variables latentes \(\delta^i\) están censuradas y tenemos que marginalizar con respecto a \(\delta^i\), resultando en:
\[\sum_{i=1}^N \log(\pi \phi_{\theta_1} (x^i)+ (1-\pi)\phi_{\theta_2}(x^i)).\]
donde \(\pi\) es la probabilidad de que la observación venga de la primera densidad. Este problema es más difícil pues tenemos tanto \(\pi\) como \(\theta_1\) y \(\theta_2\) dentro del logaritmo. Podemos resolver numéricamente como sigue:
crearLogLike <- function(x){
logLike <- function(theta){
pi <- exp(theta[1]) / (1 + exp(theta[1]))
mu_1 <- theta[2]
mu_2 <- theta[3]
sigma_1 <- exp(theta[4])
sigma_2 <- exp(theta[5])
sum(log(pi*dnorm(x, mu_1, sd=sigma_1)+(1-pi)*dnorm(x,mu_2,sd=sigma_2)))
}
logLike
}
func_1 <- crearLogLike(c(x_1,x_2))
system.time(salida <- optim(c(0.5,0,0,1,1), func_1, control=list(fnscale=-1)))## user system elapsed
## 0.033 0.001 0.035salida$convergence## [1] 0exp(salida$par[1]) / (1 + exp(salida$par[1]))## [1] 0.59salida$par[2:3]## [1] 1.99 0.03exp(salida$par[4:5])## [1] 0.30 0.52Y vemos que hemos podido recuperar los parámetros originales.
Ahora implementamos EM para resolver este problema. Empezamos con la log-verosimilitud para datos completos (que reescribimos de manera más conveniente): \[\sum_{i=1}^N \delta^i\log\phi_{\theta_1} (x^i)+ (1-\delta^i)\log\phi_{\theta_2}(x^i) + \delta^i \log\pi + (1-\delta^i)\log(1-\pi).\]
Tomamos valores iniciales para los parámetros \(\hat{\mu}_1,\hat{\mu}_2,\hat{\sigma}_1^2, \hat{\sigma}_2^2, \hat{\pi}\) y comenzamos con el paso Esperanza promediando sobre las variables aleatorias, que en este caso son las \(\delta^i\). Calculamos entonces \[\hat{\gamma}^i=E_{\hat{\theta}}(\delta^i|x^i)=P(\delta^i=1|x^i),\] y usamos bayes para expresar en términos de los parámetros: \[\hat{\gamma}^i= \frac{\hat{\pi}\phi_{\hat{\theta_1}}} {\hat{\pi}\phi_{\hat{\theta_1}}(x_i)+(1-\hat{\pi})\phi_{\hat{\theta_2}}(x_i)}\]
\(\hat{\gamma}^i\) se conocen como la responsabilidad del modelo 1 para explicar la i-ésima observación.
Utilizando estas asignaciones de los faltantes pasamos al paso Maximización, donde la función objetivo es: \[\sum_{i=1}^N \hat{\gamma}^i\log \phi_{\theta_1} (x^i)+ (1-\hat{\gamma}^i)\log\phi_{\theta_2}(x^i) + \hat{\gamma}^i \log\pi + (1-\hat{\gamma}^i)\log(1-\pi).\]
La actualización de \(\pi\) es fácil:
\[\hat{\pi}=\frac{1}{N}\sum_i{\gamma^i}.\]
y se puede ver sin mucha dificultad que
\[\hat{\mu}_1=\frac{\sum_i\hat{\gamma}^i x^i}{\sum_i \hat{\gamma}^i}\] \[\hat{\mu}_2=\frac{\sum_i(1-\hat{\gamma}^i) x^i}{\sum_i (1-\hat{\gamma}^i})\]
\[\hat{\sigma}_1^2=\frac{\sum_i\hat{\gamma}^i (x^i-\mu_1)^2}{\sum_i \hat{\gamma}^i}\] \[\hat{\sigma}_2^2=\frac{\sum_i(1-\hat{\gamma}^i) (x^i-\mu_2)^2}{\sum_i (1-\hat{\gamma}^i)},\]
Implementa EM para el ejemplo de mezcla de normales.
Mezclas gaussianas más general
Un caso más general es suponer que hay \(K\) posibles clases y las distribuciones de la mezcla son normal multivariadas. De tal manera que \[P(\delta_k^i=1)=\pi_k, \sum_{k=1}^K \pi_k =1\]
Entonces, la distribución de los datos observados es la mezcla Gaussiana:
\[p(x)=\sum_{k=1}^K \pi_k p_{\theta_k}(x)\]
donde \(p_{\theta_k}\) es normal multivariada con parámetros \(\theta_k=\{\mu_k, \Sigma_k\}\),
\[p_{\theta_k}(x)=\frac{1}{(2\pi|\Sigma_k|)^{1/2}}exp\big\{\frac{1}{2}(x-\mu_k)^T \Sigma_k^{-1}(x-\mu_k)\big\}\]
La estimación en el caso multivariado con \(K\) clases se realiza con el algoritmo EM de manera análoga al caso de dos clases y normal univariada; sin embargo, es posible restringir el tipo de mezcla a clases de distribuciones Gaussianas determinadas por la matriz de covarianzas.
Consideremos la descomposición de espectral de la \(k\)-ésima matriz de covarianzas:
\[\Sigma_k=\lambda_k D_k A_k D_k^T\]
\(A_k\) es una matriz diagonal con \(det(A_k)=1\) y entradas diagonales proporcionales a los eigenvalores de \(\Sigma_k\). \(A_k\) controla la forma (shape) del \(k\)-ésimo cluster.
\(\lambda_k\) es una constante igual al \(det(\Sigma_k)\) que controla el volumen del \(k\)-ésimo cluster.
\(D_k\) es la matriz ortonormal de eigenvectores y controla la orientación del \(k\)-ésimo cluster.
Utilizando la descomposición de arriba hay modelos de mezclas Gaussianas que se restringen a cluster con propiedades particulares, por ejemplo:
\(\Sigma_k=\lambda I\): clusters esféricos de igual volumen.
\(\Sigma_k=\lambda_k I\): clusters esféricos de distinto volumen.
\(\Sigma_k=\lambda_k A\): clusters elipsoidales de distinto volumen pero misma forma , orientaciones alineadas con el eje.
…
- \(\Sigma_k=\lambda_k D_k A_k D_k^T\) modelo sin restricciones.
Ahora veremos como ajustar estos modelos. En R hay varios paquetes para ajustar mezclas gaussianas, dos de ellos son mclust y mixtools. También está flexmix
Veamos un ejemplo usando mclust, en este paquete los modelos disponibles son
?mclustModelNames
Mezclas Univariadas
"E" = equal variance (one-dimensional)
"V" = variable variance (one-dimensional)
Mezclas multivariadas
"EII" = spherical, equal volume
"VII" = spherical, unequal volume
"EEI" = diagonal, equal volume and shape
"VEI" = diagonal, varying volume, equal shape
"EVI" = diagonal, equal volume, varying shape
"VVI" = diagonal, varying volume and shape
"EEE" = ellipsoidal, equal volume, shape, and orientation
"EEV" = ellipsoidal, equal volume and equal shape
"VEV" = ellipsoidal, equal shape
"VVV" = ellipsoidal, varying volume, shape, and orientationLa nomenclatura es: E=equal, V=varying, I=identity, y las letras están ordenadas 1ra=volumen, 2a=forma, 3ra=orientación.
En todos los modelos hay \(K-1\) parámetros para las probabilidades iniciales \(\pi_k\) mas \(Kp\) parámetros para las medias mas los parámetros de la matriz de covarianzas:
Los modelos con menos restricciones tienen más parámetros por estimar y por tanto necesitan de más datos para alcanzar la misma precisión.
Usaremos los datos wine.
library(mclust)
wine <- read.csv("wine.csv", stringsAsFactors=FALSE)[, -1]
w_mclust <- Mclust(wine)
summary(w_mclust)## ----------------------------------------------------
## Gaussian finite mixture model fitted by EM algorithm
## ----------------------------------------------------
##
## Mclust VEI (diagonal, equal shape) model with 6 components:
##
## log.likelihood n df BIC ICL
## -3207 178 101 -6938 -6953
##
## Clustering table:
## 1 2 3 4 5 6
## 40 18 24 27 28 41Podemos ver más detalles:
summary(w_mclust, parameters = TRUE)## ----------------------------------------------------
## Gaussian finite mixture model fitted by EM algorithm
## ----------------------------------------------------
##
## Mclust VEI (diagonal, equal shape) model with 6 components:
##
## log.likelihood n df BIC ICL
## -3207 178 101 -6938 -6953
##
## Clustering table:
## 1 2 3 4 5 6
## 40 18 24 27 28 41
##
## Mixing probabilities:
## 1 2 3 4 5 6
## 0.21 0.11 0.14 0.15 0.16 0.23
##
## Means:
## [,1] [,2] [,3] [,4] [,5] [,6]
## Alcohol.content 13.70 13.9 12.3 13.27 12.92 12.23
## Malic.acid 2.09 1.8 2.2 3.34 2.95 1.83
## As 2.42 2.5 2.4 2.46 2.32 2.23
## Alcalinity.in.ash 16.88 16.8 20.8 21.91 19.86 20.52
## Magnesium 105.73 106.5 105.0 100.83 96.84 89.14
## Total.phenols 2.73 3.1 2.7 1.77 1.63 2.07
## Falavanoids 2.81 3.3 2.7 0.86 0.81 1.91
## Nonflavanoid.phenols 0.28 0.3 0.3 0.44 0.46 0.38
## Proantocyanins 1.78 2.1 2.2 1.36 0.89 1.45
## Color.Intensity 4.98 6.7 3.4 9.11 4.77 2.86
## Hue 1.04 1.1 1.1 0.62 0.82 1.07
## OD280.OD315 3.27 2.9 3.0 1.64 1.77 2.84
## Proline 1026.23 1322.6 587.9 647.84 602.20 477.66
##
## Variances:
## [,,1]
## Alcohol.content Malic.acid As Alcalinity.in.ash
## Alcohol.content 0.19 0.00 0.000 0.0
## Malic.acid 0.00 0.65 0.000 0.0
## As 0.00 0.00 0.046 0.0
## Alcalinity.in.ash 0.00 0.00 0.000 5.4
## Magnesium 0.00 0.00 0.000 0.0
## Total.phenols 0.00 0.00 0.000 0.0
## Falavanoids 0.00 0.00 0.000 0.0
## Nonflavanoid.phenols 0.00 0.00 0.000 0.0
## Proantocyanins 0.00 0.00 0.000 0.0
## Color.Intensity 0.00 0.00 0.000 0.0
## Hue 0.00 0.00 0.000 0.0
## OD280.OD315 0.00 0.00 0.000 0.0
## Proline 0.00 0.00 0.000 0.0
## Magnesium Total.phenols Falavanoids
## Alcohol.content 0 0.000 0.00
## Malic.acid 0 0.000 0.00
## As 0 0.000 0.00
## Alcalinity.in.ash 0 0.000 0.00
## Magnesium 101 0.000 0.00
## Total.phenols 0 0.096 0.00
## Falavanoids 0 0.000 0.11
## Nonflavanoid.phenols 0 0.000 0.00
## Proantocyanins 0 0.000 0.00
## Color.Intensity 0 0.000 0.00
## Hue 0 0.000 0.00
## OD280.OD315 0 0.000 0.00
## Proline 0 0.000 0.00
## Nonflavanoid.phenols Proantocyanins Color.Intensity
## Alcohol.content 0.0000 0.00 0.00
## Malic.acid 0.0000 0.00 0.00
## As 0.0000 0.00 0.00
## Alcalinity.in.ash 0.0000 0.00 0.00
## Magnesium 0.0000 0.00 0.00
## Total.phenols 0.0000 0.00 0.00
## Falavanoids 0.0000 0.00 0.00
## Nonflavanoid.phenols 0.0077 0.00 0.00
## Proantocyanins 0.0000 0.11 0.00
## Color.Intensity 0.0000 0.00 0.81
## Hue 0.0000 0.00 0.00
## OD280.OD315 0.0000 0.00 0.00
## Proline 0.0000 0.00 0.00
## Hue OD280.OD315 Proline
## Alcohol.content 0.000 0.000 0
## Malic.acid 0.000 0.000 0
## As 0.000 0.000 0
## Alcalinity.in.ash 0.000 0.000 0
## Magnesium 0.000 0.000 0
## Total.phenols 0.000 0.000 0
## Falavanoids 0.000 0.000 0
## Nonflavanoid.phenols 0.000 0.000 0
## Proantocyanins 0.000 0.000 0
## Color.Intensity 0.000 0.000 0
## Hue 0.016 0.000 0
## OD280.OD315 0.000 0.085 0
## Proline 0.000 0.000 16386
## [,,2]
## Alcohol.content Malic.acid As Alcalinity.in.ash
## Alcohol.content 0.18 0.00 0.000 0.0
## Malic.acid 0.00 0.62 0.000 0.0
## As 0.00 0.00 0.044 0.0
## Alcalinity.in.ash 0.00 0.00 0.000 5.2
## Magnesium 0.00 0.00 0.000 0.0
## Total.phenols 0.00 0.00 0.000 0.0
## Falavanoids 0.00 0.00 0.000 0.0
## Nonflavanoid.phenols 0.00 0.00 0.000 0.0
## Proantocyanins 0.00 0.00 0.000 0.0
## Color.Intensity 0.00 0.00 0.000 0.0
## Hue 0.00 0.00 0.000 0.0
## OD280.OD315 0.00 0.00 0.000 0.0
## Proline 0.00 0.00 0.000 0.0
## Magnesium Total.phenols Falavanoids
## Alcohol.content 0 0.000 0.0
## Malic.acid 0 0.000 0.0
## As 0 0.000 0.0
## Alcalinity.in.ash 0 0.000 0.0
## Magnesium 97 0.000 0.0
## Total.phenols 0 0.091 0.0
## Falavanoids 0 0.000 0.1
## Nonflavanoid.phenols 0 0.000 0.0
## Proantocyanins 0 0.000 0.0
## Color.Intensity 0 0.000 0.0
## Hue 0 0.000 0.0
## OD280.OD315 0 0.000 0.0
## Proline 0 0.000 0.0
## Nonflavanoid.phenols Proantocyanins Color.Intensity
## Alcohol.content 0.0000 0.0 0.00
## Malic.acid 0.0000 0.0 0.00
## As 0.0000 0.0 0.00
## Alcalinity.in.ash 0.0000 0.0 0.00
## Magnesium 0.0000 0.0 0.00
## Total.phenols 0.0000 0.0 0.00
## Falavanoids 0.0000 0.0 0.00
## Nonflavanoid.phenols 0.0073 0.0 0.00
## Proantocyanins 0.0000 0.1 0.00
## Color.Intensity 0.0000 0.0 0.77
## Hue 0.0000 0.0 0.00
## OD280.OD315 0.0000 0.0 0.00
## Proline 0.0000 0.0 0.00
## Hue OD280.OD315 Proline
## Alcohol.content 0.000 0.000 0
## Malic.acid 0.000 0.000 0
## As 0.000 0.000 0
## Alcalinity.in.ash 0.000 0.000 0
## Magnesium 0.000 0.000 0
## Total.phenols 0.000 0.000 0
## Falavanoids 0.000 0.000 0
## Nonflavanoid.phenols 0.000 0.000 0
## Proantocyanins 0.000 0.000 0
## Color.Intensity 0.000 0.000 0
## Hue 0.015 0.000 0
## OD280.OD315 0.000 0.081 0
## Proline 0.000 0.000 15604
## [,,3]
## Alcohol.content Malic.acid As Alcalinity.in.ash
## Alcohol.content 0.53 0.0 0.00 0
## Malic.acid 0.00 1.8 0.00 0
## As 0.00 0.0 0.13 0
## Alcalinity.in.ash 0.00 0.0 0.00 15
## Magnesium 0.00 0.0 0.00 0
## Total.phenols 0.00 0.0 0.00 0
## Falavanoids 0.00 0.0 0.00 0
## Nonflavanoid.phenols 0.00 0.0 0.00 0
## Proantocyanins 0.00 0.0 0.00 0
## Color.Intensity 0.00 0.0 0.00 0
## Hue 0.00 0.0 0.00 0
## OD280.OD315 0.00 0.0 0.00 0
## Proline 0.00 0.0 0.00 0
## Magnesium Total.phenols Falavanoids
## Alcohol.content 0 0.00 0.0
## Malic.acid 0 0.00 0.0
## As 0 0.00 0.0
## Alcalinity.in.ash 0 0.00 0.0
## Magnesium 287 0.00 0.0
## Total.phenols 0 0.27 0.0
## Falavanoids 0 0.00 0.3
## Nonflavanoid.phenols 0 0.00 0.0
## Proantocyanins 0 0.00 0.0
## Color.Intensity 0 0.00 0.0
## Hue 0 0.00 0.0
## OD280.OD315 0 0.00 0.0
## Proline 0 0.00 0.0
## Nonflavanoid.phenols Proantocyanins Color.Intensity
## Alcohol.content 0.000 0.00 0.0
## Malic.acid 0.000 0.00 0.0
## As 0.000 0.00 0.0
## Alcalinity.in.ash 0.000 0.00 0.0
## Magnesium 0.000 0.00 0.0
## Total.phenols 0.000 0.00 0.0
## Falavanoids 0.000 0.00 0.0
## Nonflavanoid.phenols 0.022 0.00 0.0
## Proantocyanins 0.000 0.31 0.0
## Color.Intensity 0.000 0.00 2.3
## Hue 0.000 0.00 0.0
## OD280.OD315 0.000 0.00 0.0
## Proline 0.000 0.00 0.0
## Hue OD280.OD315 Proline
## Alcohol.content 0.000 0.00 0
## Malic.acid 0.000 0.00 0
## As 0.000 0.00 0
## Alcalinity.in.ash 0.000 0.00 0
## Magnesium 0.000 0.00 0
## Total.phenols 0.000 0.00 0
## Falavanoids 0.000 0.00 0
## Nonflavanoid.phenols 0.000 0.00 0
## Proantocyanins 0.000 0.00 0
## Color.Intensity 0.000 0.00 0
## Hue 0.044 0.00 0
## OD280.OD315 0.000 0.24 0
## Proline 0.000 0.00 46366
## [,,4]
## Alcohol.content Malic.acid As Alcalinity.in.ash
## Alcohol.content 0.23 0.00 0.000 0.0
## Malic.acid 0.00 0.78 0.000 0.0
## As 0.00 0.00 0.056 0.0
## Alcalinity.in.ash 0.00 0.00 0.000 6.6
## Magnesium 0.00 0.00 0.000 0.0
## Total.phenols 0.00 0.00 0.000 0.0
## Falavanoids 0.00 0.00 0.000 0.0
## Nonflavanoid.phenols 0.00 0.00 0.000 0.0
## Proantocyanins 0.00 0.00 0.000 0.0
## Color.Intensity 0.00 0.00 0.000 0.0
## Hue 0.00 0.00 0.000 0.0
## OD280.OD315 0.00 0.00 0.000 0.0
## Proline 0.00 0.00 0.000 0.0
## Magnesium Total.phenols Falavanoids
## Alcohol.content 0 0.00 0.00
## Malic.acid 0 0.00 0.00
## As 0 0.00 0.00
## Alcalinity.in.ash 0 0.00 0.00
## Magnesium 123 0.00 0.00
## Total.phenols 0 0.12 0.00
## Falavanoids 0 0.00 0.13
## Nonflavanoid.phenols 0 0.00 0.00
## Proantocyanins 0 0.00 0.00
## Color.Intensity 0 0.00 0.00
## Hue 0 0.00 0.00
## OD280.OD315 0 0.00 0.00
## Proline 0 0.00 0.00
## Nonflavanoid.phenols Proantocyanins Color.Intensity
## Alcohol.content 0.0000 0.00 0.00
## Malic.acid 0.0000 0.00 0.00
## As 0.0000 0.00 0.00
## Alcalinity.in.ash 0.0000 0.00 0.00
## Magnesium 0.0000 0.00 0.00
## Total.phenols 0.0000 0.00 0.00
## Falavanoids 0.0000 0.00 0.00
## Nonflavanoid.phenols 0.0093 0.00 0.00
## Proantocyanins 0.0000 0.13 0.00
## Color.Intensity 0.0000 0.00 0.99
## Hue 0.0000 0.00 0.00
## OD280.OD315 0.0000 0.00 0.00
## Proline 0.0000 0.00 0.00
## Hue OD280.OD315 Proline
## Alcohol.content 0.000 0.0 0
## Malic.acid 0.000 0.0 0
## As 0.000 0.0 0
## Alcalinity.in.ash 0.000 0.0 0
## Magnesium 0.000 0.0 0
## Total.phenols 0.000 0.0 0
## Falavanoids 0.000 0.0 0
## Nonflavanoid.phenols 0.000 0.0 0
## Proantocyanins 0.000 0.0 0
## Color.Intensity 0.000 0.0 0
## Hue 0.019 0.0 0
## OD280.OD315 0.000 0.1 0
## Proline 0.000 0.0 19895
## [,,5]
## Alcohol.content Malic.acid As Alcalinity.in.ash
## Alcohol.content 0.27 0.00 0.000 0.0
## Malic.acid 0.00 0.92 0.000 0.0
## As 0.00 0.00 0.066 0.0
## Alcalinity.in.ash 0.00 0.00 0.000 7.7
## Magnesium 0.00 0.00 0.000 0.0
## Total.phenols 0.00 0.00 0.000 0.0
## Falavanoids 0.00 0.00 0.000 0.0
## Nonflavanoid.phenols 0.00 0.00 0.000 0.0
## Proantocyanins 0.00 0.00 0.000 0.0
## Color.Intensity 0.00 0.00 0.000 0.0
## Hue 0.00 0.00 0.000 0.0
## OD280.OD315 0.00 0.00 0.000 0.0
## Proline 0.00 0.00 0.000 0.0
## Magnesium Total.phenols Falavanoids
## Alcohol.content 0 0.00 0.00
## Malic.acid 0 0.00 0.00
## As 0 0.00 0.00
## Alcalinity.in.ash 0 0.00 0.00
## Magnesium 145 0.00 0.00
## Total.phenols 0 0.14 0.00
## Falavanoids 0 0.00 0.15
## Nonflavanoid.phenols 0 0.00 0.00
## Proantocyanins 0 0.00 0.00
## Color.Intensity 0 0.00 0.00
## Hue 0 0.00 0.00
## OD280.OD315 0 0.00 0.00
## Proline 0 0.00 0.00
## Nonflavanoid.phenols Proantocyanins Color.Intensity
## Alcohol.content 0.000 0.00 0.0
## Malic.acid 0.000 0.00 0.0
## As 0.000 0.00 0.0
## Alcalinity.in.ash 0.000 0.00 0.0
## Magnesium 0.000 0.00 0.0
## Total.phenols 0.000 0.00 0.0
## Falavanoids 0.000 0.00 0.0
## Nonflavanoid.phenols 0.011 0.00 0.0
## Proantocyanins 0.000 0.16 0.0
## Color.Intensity 0.000 0.00 1.2
## Hue 0.000 0.00 0.0
## OD280.OD315 0.000 0.00 0.0
## Proline 0.000 0.00 0.0
## Hue OD280.OD315 Proline
## Alcohol.content 0.000 0.00 0
## Malic.acid 0.000 0.00 0
## As 0.000 0.00 0
## Alcalinity.in.ash 0.000 0.00 0
## Magnesium 0.000 0.00 0
## Total.phenols 0.000 0.00 0
## Falavanoids 0.000 0.00 0
## Nonflavanoid.phenols 0.000 0.00 0
## Proantocyanins 0.000 0.00 0
## Color.Intensity 0.000 0.00 0
## Hue 0.022 0.00 0
## OD280.OD315 0.000 0.12 0
## Proline 0.000 0.00 23351
## [,,6]
## Alcohol.content Malic.acid As Alcalinity.in.ash
## Alcohol.content 0.23 0.00 0.000 0.0
## Malic.acid 0.00 0.78 0.000 0.0
## As 0.00 0.00 0.056 0.0
## Alcalinity.in.ash 0.00 0.00 0.000 6.5
## Magnesium 0.00 0.00 0.000 0.0
## Total.phenols 0.00 0.00 0.000 0.0
## Falavanoids 0.00 0.00 0.000 0.0
## Nonflavanoid.phenols 0.00 0.00 0.000 0.0
## Proantocyanins 0.00 0.00 0.000 0.0
## Color.Intensity 0.00 0.00 0.000 0.0
## Hue 0.00 0.00 0.000 0.0
## OD280.OD315 0.00 0.00 0.000 0.0
## Proline 0.00 0.00 0.000 0.0
## Magnesium Total.phenols Falavanoids
## Alcohol.content 0 0.00 0.00
## Malic.acid 0 0.00 0.00
## As 0 0.00 0.00
## Alcalinity.in.ash 0 0.00 0.00
## Magnesium 123 0.00 0.00
## Total.phenols 0 0.12 0.00
## Falavanoids 0 0.00 0.13
## Nonflavanoid.phenols 0 0.00 0.00
## Proantocyanins 0 0.00 0.00
## Color.Intensity 0 0.00 0.00
## Hue 0 0.00 0.00
## OD280.OD315 0 0.00 0.00
## Proline 0 0.00 0.00
## Nonflavanoid.phenols Proantocyanins Color.Intensity
## Alcohol.content 0.0000 0.00 0.00
## Malic.acid 0.0000 0.00 0.00
## As 0.0000 0.00 0.00
## Alcalinity.in.ash 0.0000 0.00 0.00
## Magnesium 0.0000 0.00 0.00
## Total.phenols 0.0000 0.00 0.00
## Falavanoids 0.0000 0.00 0.00
## Nonflavanoid.phenols 0.0093 0.00 0.00
## Proantocyanins 0.0000 0.13 0.00
## Color.Intensity 0.0000 0.00 0.98
## Hue 0.0000 0.00 0.00
## OD280.OD315 0.0000 0.00 0.00
## Proline 0.0000 0.00 0.00
## Hue OD280.OD315 Proline
## Alcohol.content 0.000 0.0 0
## Malic.acid 0.000 0.0 0
## As 0.000 0.0 0
## Alcalinity.in.ash 0.000 0.0 0
## Magnesium 0.000 0.0 0
## Total.phenols 0.000 0.0 0
## Falavanoids 0.000 0.0 0
## Nonflavanoid.phenols 0.000 0.0 0
## Proantocyanins 0.000 0.0 0
## Color.Intensity 0.000 0.0 0
## Hue 0.019 0.0 0
## OD280.OD315 0.000 0.1 0
## Proline 0.000 0.0 19814Y podemos ver una tabla con el BIC para cada modelo, es importante tomar en cuenta que en mclust el BIC se define como:
\[BIC = 2\log l(x;\hat{\theta}) - d\log n\]
\(d\) es el número de parámetros y \(n\) el número de observaciones. Por tanto, buscamos maximizar el BIC.
round(w_mclust$BIC)## EII VII EEI VEI EVI VVI EEE EEV VEV VVV
## 1 -27318 -27318 -8161 -8161 -8161 -8161 -7201 -7201 -7201 -7201
## 2 -24478 -24297 -7534 -7553 -7440 -7496 -7138 -7284 -7240 -7207
## 3 -23210 -22679 -7125 -7065 -7034 -7003 -7061 -7409 -7456 -7498
## 4 -22037 -21722 -7079 -7014 -7108 -6978 -7011 -7687 -7713 -7902
## 5 -21552 -20791 -7063 -6972 -7082 -7046 -6975 -8182 -8189 -8304
## 6 -20776 -20013 -7059 -6938 -7063 -7009 -7090 -8578 -8514 -8650
## 7 -19736 -19594 -7060 -6990 -7141 -7143 -7083 -9012 -8774 NA
## 8 -19288 -19107 -7024 -6958 -7226 -7225 -7019 -9317 -9150 NA
## 9 -19057 -18667 -7027 -7000 -7334 -7342 -7065 -9618 -9415 NA
## attr(,"G")
## [1] 1 2 3 4 5 6 7 8 9
## attr(,"modelNames")
## [1] "EII" "VII" "EEI" "VEI" "EVI" "VVI" "EEE" "EEV" "VEV" "VVV"
## attr(,"oneD")
## [1] FALSEw_mclust## 'Mclust' model object:
## best model: diagonal, equal shape (VEI) with 6 componentsEn la tabla de arriba, los renglones indican el número de conglomerados.
El paquete incluye algunas gráficas auxiliares.
plot(w_mclust, what = "BIC")
# hay 3 plots más:
# plot(w_mclust)Observaciones
Las mezclas Gaussianas son sensibles a datos atípicos (outliers).
Puede haber inestabilidad si los componentes Gaussianos no están separados apropiadamente o si los datos son chicos.
A comparación de otros métodos de clustering el tener un modelo probabilístico explícito hace posible verificar supuestos y tener medidas de ajuste. Por ejemplo, podemos usar BIC para elegir el número de clusters.
Desventajas: los resultados dependen de la inicialización del algoritmo EM, puede ser demasiado flexible o inestable.
Las mezclas gaussianas se pueden extender a mezclas de modelos de regresión que pueden ser modelos lineales estándar o modelos lineales generalizados. Esto esta implementado en el paquete flexmix que utiliza el algoritmo EM.
Aplicación: Background Subtraction
Background Subtraction es una técnica de procesamiento de imágenes donde se buscan extraer los objetos en el frente de la imagen. Por ejemplo, puede ser de interés detectar autos, animales o personas de una imagen. Ejemplos:
Monitoreo de tráfico (contar vehiculos, detectar y seguir vehículos).
Reconocimiento de actividades humanas (correr, caminar, brincar).
Seguimiento de objetos (seguir la bola en el tenis).
Los modelos de mezclas gaussianas son populares para esta tarea: la idea básica es suponer que la intensidad de cada pixel en una imagen se puede modelar usando un modelo de mezclas Gaussianas. Después se determina que intensidades corresponden con mayor probabilidad al fondo seleccionando como fondo las gaussianas con menos varianza y mayor evidencia (\(\pi_j\)).
En OpencCV hay tres algoritmos implementados para sustraer fondo, estos usan variantes de mezclas gaussianas:
Análisis de clase latente
El análisis de clase latente es un tipo de modelo de mezclas finitas, se utiliza cuando observamos variables categóricas y nos interesa investigar si existe alguna fuente que explique patrones en las variables observadas, o identificar y caracterizar agrupaciones de observaciones o individuos. Algunos lugares donde se aplican estos modelos son: encuestas de opinión pública, exámenes psicométricos y análisis de comportamiento del consumidor.
Más formalmente, el análisis de clase latente busca estratificar una tabla de contingencias de variables observadas con una variable no observada que induzca independencias condicionales en las variables observadas. El modelo agrupa cada observación en una clase latente (un nivel de la variable).
Ejemplo. Supongamos que aplicamos una prueba a un grupo de estudiantes y nos interesa analizar si los resultados en la prueba se pueden explicar mediante una variable latente que mida la habilidad de cada estudiante, más aún, nos gustaría clasificar a los estudiantes en grupos de acuerdo a la habilidad observada. En este ejemplo, cada pregunta en la prueba es una variable aleatoria que denotamos \(Y_i\) y toma dos valores \(0\) si la respuesta es incorrecta y \(1\) en otro caso.
Ahora simulamos este ejemplo donde suponemos que se realiza la prueba a 1,500 estudiante, el examen tiene 4 secciones (mide 4 áreas de habilidad) que representan las variables observadas \(Y_1,..., Y_4\). La probabilidad de que un estudiante ascierte cada pregunta depende del grupo de habilidad al que pertenece, en la práctica esta es la variable no observada. Adicionalmente, en nuestro ejemplo la probabilidad de pertenecer a cada grupo de habilidad es distinta.
Entonces, simulamos como sigue:
- 1,500 alumnos pertenecientes a 3 grupos, la probabilidad de pertenecer a cada grupo es 0.25, 0.5 y 0.25:
set.seed(346713443)
N <- 1500
x_0 <- rmultinom(N, size = 1, prob = c(0.25, 0.50, 0.25))
obs <- tbl_df(data.frame(id = 1:N, gpo = c(1:3 %*% (x_0))))
obs## Source: local data frame [1,500 x 2]
##
## id gpo
## 1 1 3
## 2 2 2
## 3 3 2
## 4 4 1
## 5 5 2
## 6 6 3
## 7 7 1
## 8 8 1
## 9 9 1
## 10 10 3
## .. .. ...- La probabilidad de acertar en cada área de habilidad varía de acuerdo al grupo.
pi <- tbl_df(data.frame(gpo = rep(1:3, 4), prueba = rep(1:4, each = 3),
p = c(0.05, 0.45, 0.80, 0.10, 0.50, 0.90, 0.20, 0.80, 0.75, 0.15, 0.90, 0.45)))
pi## Source: local data frame [12 x 3]
##
## gpo prueba p
## 1 1 1 0.05
## 2 2 1 0.45
## 3 3 1 0.80
## 4 1 2 0.10
## 5 2 2 0.50
## 6 3 2 0.90
## 7 1 3 0.20
## 8 2 3 0.80
## 9 3 3 0.75
## 10 1 4 0.15
## 11 2 4 0.90
## 12 3 4 0.45\(p=4\) items, cada uno puede tomar dos valores, 0 si la respuesta es incorrecta y 1 en otro caso.
obs_items <- inner_join(obs, pi) %>%
group_by(id, gpo, prueba) %>%
mutate(Y = rbinom(1, 1, p)) %>%
select(-p) %>%
spread(prueba, Y)## Joining by: "gpo"colnames(obs_items)[3:6] <- paste("Y", 1:4, sep = "_")
obs_items## Source: local data frame [1,500 x 6]
##
## id gpo Y_1 Y_2 Y_3 Y_4
## 1 1 3 1 0 1 0
## 2 2 2 0 0 1 1
## 3 3 2 1 0 1 0
## 4 4 1 0 0 0 0
## 5 5 2 1 1 0 1
## 6 6 3 1 1 1 0
## 7 7 1 0 0 0 0
## 8 8 1 0 0 0 0
## 9 9 1 0 0 0 0
## 10 10 3 1 1 0 0
## .. .. ... ... ... ... ...Veamos la red bayesiana asociada a las variables observadas que generamos,
library(igraph)
library(bnlearn)
obs_items_f <- obs_items %>%
mutate_each(funs(factor)) %>%
data.frame
net_obs <- hc(obs_items_f[, 3:6], score='bic')
graphviz.plot(net_obs)
Ahora veamos que ocurre cuando añadimos la variable latente que similumaos
net_obs_cl <- hc(obs_items_f[, 2:6], score='bic')
graphviz.plot(net_obs_cl)
¿Cuántos parámetros requiere cada modelo?
- 15 y 11
- 16 y 12
- 15 y 10
- 15 y 8
- Ninguno de estos
Realizamos la estimación de la variable latente utilizando el paquete poLCA, más adelante explicaremos los algoritmos detrás de la función que usamos.
library(poLCA)
f <- cbind(Y_1, Y_2, Y_3, Y_4) ~ 1
mod_1 <- poLCA(f, (obs_items[, 3:6] + 1), nclass = 3, maxiter = 8000) ## Conditional item response (column) probabilities,
## by outcome variable, for each class (row)
##
## $Y_1
## Pr(1) Pr(2)
## class 1: 0.23 0.77
## class 2: 0.56 0.44
## class 3: 0.94 0.06
##
## $Y_2
## Pr(1) Pr(2)
## class 1: 0.098 0.90
## class 2: 0.529 0.47
## class 3: 0.863 0.14
##
## $Y_3
## Pr(1) Pr(2)
## class 1: 0.28 0.72
## class 2: 0.20 0.80
## class 3: 0.74 0.26
##
## $Y_4
## Pr(1) Pr(2)
## class 1: 0.58 0.419
## class 2: 0.06 0.940
## class 3: 0.91 0.095
##
## Estimated class population shares
## 0.25 0.51 0.24
##
## Predicted class memberships (by modal posterior prob.)
## 0.16 0.6 0.23
##
## =========================================================
## Fit for 3 latent classes:
## =========================================================
## number of observations: 1500
## number of estimated parameters: 14
## residual degrees of freedom: 1
## maximum log-likelihood: -3833
##
## AIC(3): 7693
## BIC(3): 7767
## G^2(3): 2 (Likelihood ratio/deviance statistic)
## X^2(3): 2 (Chi-square goodness of fit)
## obs_items_f$pred_class <- factor(mod_1$predclass, levels = c(1, 2, 3), labels = 1:3)
table(obs_items_f$gpo, obs_items_f$pred_class) ##
## 1 2 3
## 1 14 63 306
## 2 49 678 30
## 3 178 166 16round(prop.table(table(obs_items_f$gpo, obs_items_f$pred_class), 1) * 100)##
## 1 2 3
## 1 4 16 80
## 2 6 90 4
## 3 49 46 4¿Qué pasa si estimamos una nueva red bayesiana incluyendo la variable latente estimada?
net_obs <- hc(obs_items_f[, 3:7], score='bic')
BIC(net_obs, data = obs_items_f[, 3:7])## [1] -3923graphviz.plot(net_obs)
net_obs_2 <- hc(obs_items_f[, 3:7], score='bic', k = 8)
BIC(net_obs_2, data = obs_items_f[, 3:7])## [1] -4003graphviz.plot(net_obs_2)
Estimación del modelo de clase latente
Comencemos con notación:
Observamos \(Y=(Y_1,...Y_p)\) variables categóricas, cada una con \(K_j\) posibles valores (niveles de la variable)
\(y_{j,k}^i = 1\) si observamos el \(k\)-ésimo valor, en la \(j\)-ésima variable, para \(i\)-ésimo individuo, \(y_{j,k}^i = 0\) en otro caso.
\(\pi_{jrk}\) denota la probabilidad condicional, de que una observación en la clase (latente) \(r\) produzca una observación \(k\) en la \(j\)-ésima variable. Por tanto para cada variable y cada clase:
\[\sum_{k=1}^{K_j} \pi_{rjk} = 1\]
Sea \(\delta\) un vector de R dimensiones donde \(\delta_r \in \{0,1\}\) con \(r=1,...,R\) y \(\sum_{r}\delta_r=1\). Entonces, \(\delta_r^i=1\) representa que la \(i\)-ésima observación pertenece a la \(r\)-ésima clase.
\(q_r\) denota las proporciones correspondientes a los ponderadores de la mezcla, esto es \(P(\delta_r=1)=q_r\) y \[\sum_r q_r= 1\] Los valores \(q_r\) también se conocen como las probabilidades iniciales de pertenencia a la clase latente pues representan las probabilidades no condicionales de que un individuo pertenezca a cada clase antes de observar las respuestas \(y_{j,k}^i\).
En el caso del ejemplo (exámen de habilidad) observamos \(p = 4\) variables binarias \(Y_j\).
head(obs_items[, 3:6])## Source: local data frame [6 x 4]
##
## Y_1 Y_2 Y_3 Y_4
## 1 1 0 1 0
## 2 0 0 1 1
## 3 1 0 1 0
## 4 0 0 0 0
## 5 1 1 0 1
## 6 1 1 1 0Las probabilidades condicionales reales para cada variable \(Y_j\) (prueba) y cada clase \(k\) (gpo) son:
spread(pi, gpo, p)## Source: local data frame [4 x 4]
##
## prueba 1 2 3
## 1 1 0.05 0.45 0.80
## 2 2 0.10 0.50 0.90
## 3 3 0.20 0.80 0.75
## 4 4 0.15 0.90 0.45y las estimadas
round(rbind(mod_1$probs$Y_1[, 2], mod_1$probs$Y_2[, 2], mod_1$probs$Y_3[, 2],
mod_1$probs$Y_4[, 2]), 2)[, c(3, 2, 1)]## class 3: class 2: class 1:
## [1,] 0.06 0.44 0.77
## [2,] 0.14 0.47 0.90
## [3,] 0.26 0.80 0.72
## [4,] 0.09 0.94 0.42Finalmente, las probabilidades simuladas de pertenecer a cada clase son 0.25, 0.50 y 0.25, mientras que estimadas son:
mod_1$P[c(3, 2, 1)]## [1] 0.24 0.51 0.25Ahora, la probabilidad de que un individuo \(i\) perteneciente a la clase \(r\) produzca un conjunto de observaciones \(y^i\) es:
\[p(y^i|\delta_r=1) = \prod_{j=1}^p\prod_{k=1}^{K_j}(\pi_{jrk})^{y_{jk}^i}\]
y la función de densidad a través de todas las clases es
\[p(y^i) = \sum_{r=1}^R p(y^i|\delta_r=1) p(\delta_r^i=1)\] \[=\sum_{r=1}^R p(y^i|\delta_r^i=1)q_r\] \[ =\sum_{r=1}^R q_r \prod_{j=1}^p\prod_{k=1}^{K_j}(\pi_{jrk})^{y_{jk}^i}\]
Escribamos la log-verosimilitud respecto a \(q_r\) y \(\pi_{jrk}\):
\[\log\mathcal{L}=\sum_{i=1}^N \log \big\{\sum_{r=1}^R q_r \prod_{j=1}^p\prod_{k=1}^{K_j}(\pi_{jrk})^{y_{jk}^i}\big\} \]
Notamos que este problema es un buen candidato para EM ya que la clase de cada individuo es un dato faltante.
Iniciamos con valores arbitrarios \(\hat{q}_r\) y \(\hat{\pi}_{jrk}\).
En el paso de esperanza hace falta calcular la esperanza de las variables faltantes condicional a \(\hat{q}_r\) y \(\hat{\pi}_{jrk}\). Esto es muy similar al caso de mezclas gaussianas: \[\hat{\gamma}^i_r=E_{\hat{\theta}}(\delta^i_r|y^i)=P_{\hat{\theta}}(\delta^i_r=1|y^i),\] podemos usar la fórmula de Bayes para calcular la probabilidad de que cada individuo pertenezca a cada clase (condicional en las variables observadas):
\[\hat{\gamma}^i_r = \frac{p_{\hat{\theta}}(y^i|\delta_r^i=1)\hat{q}_r}{\sum_{l=1}^R {p_{\hat{\theta}}(y^i|\delta_l^i=1)\hat{q}_l}}\]
- En el paso de maximización actualizamos los valores de los parámetros maximizando la función de verosimilitud dadas \(\hat{\gamma}^i_r\):
\[\hat{q}_r = \frac{1}{N} \sum_{i=1}^N \gamma_r^i\] \[\hat{\pi}_{jr} = \frac{\sum_{i=1}^N y_{j}^i \cdot \gamma_r^i}{\sum_{l=1}^N \gamma_r^i}\]
Los pasos 2 y 3 se repiten hasta alcanzar convergencia.
Vale la pena notar que el modelo no determina el número de clases; sin embargo, a diferencia de otras técnicas para crear clusters, en los modelos de clase latente tenemos medidas de ajusre que nos pueden euxiliar en determinar el número de grupos. Una manera de proceder es comenzar con el modelo de un solo grupo e ir incrementando el número de clases latentes (y por tanto la complejidad del modelo) usando como criterio el AIC o BIC.
mod_1c <- poLCA(f, (obs_items[, 3:6] + 1), nclass = 1, verbose = FALSE)
mod_1c$aic## [1] 8098mod_1c$bic## [1] 8119mod_2c <- poLCA(f, (obs_items[, 3:6] + 1), nclass = 2, verbose = FALSE)
mod_2c$aic## [1] 7759mod_2c$bic## [1] 7806mod_3c <- poLCA(f, (obs_items[, 3:6] + 1), nclass = 3, verbose = FALSE)
mod_3c$aic## [1] 7693mod_3c$bic## [1] 7767mod_4c <- poLCA(f, (obs_items[, 3:6] + 1), nclass = 4, verbose = FALSE)
mod_4c$aic## [1] 7701mod_4c$bic## [1] 7802Es importante tener en cuenta que el número de parámetros estimados crece rápidamente con los valores de \(R\), \(J\) y \(K_j\), dados estos valores el número de parámetros es \(R\sum_j (K_j -1) + (R-1)\). Esto también limita la complejidad del modelo que se puede ajustar.
Observaciones
El análisis de clase latente es un tipo de modelo de mezclas finitas. Las distribuciones que componen las mezcla son tablas de contingencia de la misma dimensión que la tabla observada. Siguiendo el supuesto de independencia condicional, la frecuencia de cada celda en cada componente de la mezcla es simplemente el producto de las frecuencias marginales correspondientes a la clase dada.
Observaciones con conjuntos de respuestas similares tenderán a agruparse en las mismas clases latentes.
Recursos
poLCA. Artículo del paquete que utilizamos.
Blaydes y Linzer (2008). Modelos de clase latente con regresión, estos modelos generalizan la idea básica de modelos de clase latente permitiendo la inclusión de covariables para predecir la pertencia a una clase.
Otros ejemplos de clase latente en ciencias políticas: Hill y Kriesi (2001), Breen (200), Linzer (2011).
La teoría de respuesta al reactivo abarca los modelos de clase latente donde se observa un conjunto de variables categóricas \((Y_1, ...Y_p)\) que son condicionalmente independientes dado un conjunto de variables latentes continuas \((X_1,...,X_q)\) (con \(q<<p\)). \ Las variables no observables pueden ser inteligencia, habilidad matemática o verbal, actitud política o preferencias de consumidor, por tanto las aplicaciones se extienden a pruebas de educación, psicológicas, de sociología y marketing.\ Un ejemplo clásico de estos modelos es en la medición de inteligencia, analizaremos en particular el modelo de {} donde suponemos que la inteligencia de un individuo es una variable aleatoria \(X\) que se mide en una escala continua. Se realiza una prueba de inteligencia utilizando una batería de \(p\) preguntas, donde un individuo recibe un puntaje \(Y_j = 1\) si obtiene la respuesta correcta y \(Y_j = 0\) en otro caso (de manera similar al ejemplo anterior). La batería de preguntas se utilizará para estimar la inteligencia de cada individuo usando: \[E(X|Y_1=y_1,...,Y_p = y_p)\] como la estimación de inteligencia para un individuo con resultados \(y=(y_1,...,y_p)\). Suponemos que, \[\pi_j(x)=p(Y_j=1|X=x)=\frac{e^{x - \beta_j}}{1+ e^{x - \beta_j}}\] esto es, \[logit\pi_j(x) = x - \beta_i \] donde \(x\) es la inteligencia y \(\beta_j\) es la dificultad del reactivo \(j\), para la estimación suponemos además que \(X \sim N(0,1)\) y se procede utilizando el algoritmo EM. Analizaremos los resultados del exámen de admisión a la escuela de derecho (LSAT) de 1000 individuos en 5 preguntas. <