Modelos de Markov Gaussianos

Danna L. Cruz
Profesora Asistente de Carrera
Universidad del Rosario
Vicerrectoría de investigación y consultoría
Escuela de medicina

Efecto aleatorio

Resultado sin efecto aleatorio

Mapa de Colombia

Definición de vecindad

Grafo de Colombia

Campos aleatorios gaussianos de Markov

Sea un grafo \(\mathcal{G} = (\mathcal{V}, \mathcal{E})\) con vértices \(n\) donde cada vértice representa uno de los componentes del vector \(\boldsymbol{\theta} = (\theta_1, \theta_2, \dots, \theta_n)\) y las aristas conectan nodos que tienen algún tipo de asociación. Un GMRF supone que \(\boldsymbol{\theta} = (\theta_1, \theta_2, \dots, \theta_n)^{t} \sim N (\boldsymbol {\mu}, \boldsymbol{\boldsymbol{\Sigma}})\) y que las aristas del grafo conectan nodos \(i\) y \(j\) si y solo si \(\theta_i \not\perp \theta_j | \boldsymbol {\theta}_{- ij}\), es decir, si \(\theta_i\) es independiente de \(\theta_{j}\), dados los componentes de \(\boldsymbol{\theta}\) excepto \(\theta_i\) y \(\theta_j\).

Si \(\boldsymbol{\theta} \sim N(\boldsymbol{\mu}, \boldsymbol{ Q})\), entonces para \(i\neq j\),

\[ \theta_i \perp \theta_j | \boldsymbol{\theta}_{-ij} \Leftrightarrow Q_{ij}=0. \]

Campos aleatorios gaussianos de Markov

Un vector aleatorio \(\boldsymbol{\theta}=(\theta_1, \theta_2, \dots, \theta_n)^{t} \in \mathbb{R}^{n}\) es llamado GMRF correspondiente a un grafo \(\mathcal{G}=(\mathcal{V}, \mathcal{E})\) con media \(\boldsymbol{\mu}\) y matriz de precisión \(\boldsymbol{ Q} >0\), si y solamente si la función de densidad de probabilidad de \(\boldsymbol{\theta}\) tiene la siguiente forma:

\[ \pi(\boldsymbol{\theta})=(2 \pi)^{-n/2}|\boldsymbol{ Q}|^{1/2}\exp\left(-\frac{1}{2} (\boldsymbol{\theta}-\boldsymbol{\mu})^{t} \boldsymbol{ Q}(\boldsymbol{\theta}-\boldsymbol{\mu}) \right), \]

donde la matriz \(\boldsymbol{ Q}\) cumple la siguiente condición: \[ Q_{ij}\neq 0 \Leftrightarrow \{i, j \}\in \mathcal{E}, \forall i \neq j. \]

Propiedades

Sea un grafo \(\mathcal{G}=(\mathcal{V}, \mathcal{E})\) que representa \(\boldsymbol{\theta}\), con media \(\boldsymbol{\mu}\) y matriz de precisión \(\boldsymbol{ Q}\) simétrica y definitiva positiva. Luego, la distribución de cada componente \(\theta_i\) de \(\boldsymbol{\theta}\), dado el vector \(\boldsymbol{\theta}_{-i}\) formado por todos los componentes de \(\boldsymbol{\theta}\), excepto \(\theta_{i}\) es una distribución normal tal que:

\[ E(\theta_i | \theta_{-i})=\mu_i-\frac{1}{Q_{ii}} \sum_{j: j\sim i} Q_{ij}(\theta_j - \mu_j) \]

\[ \text{Prec} (\theta_i | \theta_{-i}) = Q_{ii}, \] \[ Corr(\theta_i, \theta_j | \boldsymbol{\theta}_{-ij}) = -\frac{Q_{ij}}{\sqrt{Q_{ii}Q_{jj}}}, \,\,\,\,\, i \neq j, \]

donde \(i\sim j\) denota que el nodo \(j\) es vecino de nodo \(i\).

Campos aleatorios gaussianos de Markov

Sea un grafo \(\mathcal{G}=(\mathcal{V}, \mathcal{E})\) con \(n\) vértices \(\mathcal{V}\) que representa al vector \(\boldsymbol{\theta}=(\theta_1, \theta_2, \dots, \theta_n)\) y \(\mathcal{E}\) el conjunto de aristas que conecta dos vértices \(\theta_i\) e \(\theta_j\). Así,

\[ f(\theta_i | \boldsymbol{\theta}_{-i})=f(X_i | \boldsymbol{\theta}_{j \sim i}), \]

con \(\boldsymbol{\theta}_{j \sim i}\) el vector formado por todos los componentes de \(\boldsymbol{\theta}\) quienes son vecinos de \(i\).

Un ejemplo de Campos de Markov Gaussianos: el modelo CAR

Considere:

• Una región espacial que está dividida en subregiones \(n\) indexadas por enteros \(1, 2, \dots, n\).

• Dotada de un sistema de vecindad \(\{V_i: i: 1, \dots, n\}\),

\[ V_{i} =\{j: \text{subregiones i y j que comparten borde}\}, \]

para \(i \in \{1, 2, \dots, n\}\).

Un ejemplo de Campos de Markov Gaussianos: el modelo CAR

Entonces las distribuciones condicionales completas para efectos aleatorios dados por \(\boldsymbol{\theta}=\theta_i | \theta_{-i}\), \(i = 1, \dots, n\) son:

\[ \theta_i | \theta_{-i} \sim N(\mu_i + \rho \overline{\theta_i}, \frac{\sigma ^2} {d _i} ), \] donde: \(\overline{\theta_i}=\sum_{\mathcal{E} } (d_i^{ \mathcal{G}})^{-1} (\theta_j)\) y \(\mathcal{E }=\{(i,j) \in E(\mathcal{G}): j \sim i \}\) es el conjunto de aristas que pertenecen al grafo \(\mathcal{G}\).

Un ejemplo de Campos de Markov Gaussianos: el modelo CAR

Sea la matriz de adyacencia:

\[ A = \begin{cases} 0, \,\,\,\,\,\text{si}\,\,\,\,\,i =j\\ 1, \,\,\,\,\,\text{si}\,\,\,\,\,i \sim j\\ 0, \,\,\,\,\,\text{si}\,\,\,\,\,i \not\sim j\\ \end{cases} \]

tome \(M = diag\{d_1 , d_2 , \dots, d_n \}\).

\[ \boldsymbol{\theta} \sim N(\boldsymbol{0}, \sigma^2 \boldsymbol{\Sigma} ), \]

donde \(\boldsymbol{\Sigma}= \sigma ^{2} (M -\rho A )^{-1}\). Para que la matriz de covarianza sea definida positiva, es necesario que \(\rho <\frac{1}{\lambda_1}\) donde \(\lambda_1\) es el valor propio más pequeño de la matriz \(M ^{-1/2} A M^{-1/2}\).

Cargar datos a R

library(spdep)
library(rgdal)
head(data.frame(shp))

##   OBJECTID SCACODIGO SCATIPO            SCANOMBRE
## 0        1    002605       0              BOLONIA
## 1        2    005103       0           SAN MIGUEL
## 2        3    005615       0     LA SOLEDAD NORTE
## 3        4    205318       2 SAN BERNARDINO XVIII
## 4        5    008104       0     LA PERSEVERANCIA
## 5        6    008109       0           SAN MARTIN
##                                 GLOBALID Shape_Leng Shape_Area
## 0 {2A535970-1BC4-4BDE-BE28-45C6A0EF1F6C}  2822.0357   136874.8
## 1 {9AC3041D-B72E-4C53-A969-00AA1F8F0B50}  2091.4844   247110.2
## 2 {19209EE3-4C37-4214-8A63-D9CB4FBECE38}  2175.3285   283351.2
## 3 {51A8B7DD-DAC9-4377-9865-53D24E70BB9F}   941.3076    55498.3
## 4 {B7BE793E-AAB1-44D6-A3E4-FCD2F749AB96}  1407.0933   104365.6
## 5 {B31AA042-875A-47D2-904C-066C6EFF6C42}  1625.5442   138656.9

w <- poly2nb(shp)

Mapa de Bogotá

par(mai=c(0,0,0,0))
plot(shp, col=2:7)
xy <- coordinates(shp)
points(xy, cex=0.2, pch=20, col='white')

Grafo de Bogotá

plot(w, xy, col='red', lwd=0.1)

Resultados

leaflet(map) %>%
  addTiles() %>%
  addPolygons(color = "white", fillColor = ~ pal(Freq),
    fillOpacity = 1, popup = state_popup) %>%
  addLegend(pal = pal, values = ~Freq, opacity = 1)

Problemas con el modelo CAR

• En problemas de alta dimensión, la inversión de \(\boldsymbol{\Sigma}= \sigma ^{2} (M -\rho A )^{-1}\) es un desafío.
• Es difícil interpretar el parámetro \(\rho\).
• Problemas con alta heterogeneidad.

Propuesta

\[ \boldsymbol{y | \boldsymbol{ \theta}, \tau_{\boldsymbol{y}}} \sim N( \boldsymbol{X\beta}+\boldsymbol{ \theta}, \frac{1}{\tau_{\boldsymbol{y}}} \boldsymbol{I}) \]

\[ \boldsymbol{\theta} | \boldsymbol{\Sigma} \sim N(\boldsymbol{0}, (\nu-n-1) \boldsymbol{\Sigma}) \] \[ \boldsymbol{\Sigma} \sim Inversa-Wishart_n(\nu, \tau_{\boldsymbol{\theta}}^{-1}\boldsymbol{C}(\boldsymbol{M}-\gamma \boldsymbol{A})^{-1}\boldsymbol{C}^t) \]

Recuperación de imagenes

• Suponga que un espacio bidimensional se divide en \(N\) píxeles, etiquetados por los números enteros \(1, 2, \dots, N\). Las ubicaciones de los píxeles formarán una red cuadrada regular.

• Cada píxel es una variable aleatoria \(\theta_i\), \(1 \leq i \leq N\) que puede tomar cualquier valor real, por tanto \(\theta_i \in \mathbb{R}\). Los valores de las variables que representan píxeles se denominan intensidades.

• Un ruido arbitrario se indicará con \(\boldsymbol{y} = \{y_1, y_2, \dots, y_ {N} \}\), así \(\boldsymbol{y} \in \mathbb{R}^{N}\).

• En general, no es posible observar \(\boldsymbol{\theta}\) directamente, sino que la imagen observada \(\boldsymbol{y}\) es una copia degradada de \(\boldsymbol{\theta}\).

Recuperación de imagenes

\[ \boldsymbol{y}=\beta+\boldsymbol{\theta}+\epsilon \]

donde \(\beta\) es una constante y para cada \(i \in \{1, \dots, N \}\), \(\epsilon\sim N (0,\sigma_y^{-1})\) y \(\epsilon_i\) e \(\epsilon_j\) son independientes cuando \(i\neq j\). Entonces, el vector \(\boldsymbol{y} = (y_1, \dots, y_N)\) tiene como función de verosimilitud:

\[ L(\boldsymbol{\theta}, \sigma, \sigma_y \rho, \boldsymbol{y}) =(2 \pi \sigma_y^2)^{-n/2} \exp\left\{-\frac{1}{2} (\boldsymbol{y}-\boldsymbol{\theta}-\beta)^{t}(\boldsymbol{y}-\boldsymbol{\theta}-\beta) \right\} \]

Recuperación de imagenes

Imagen de ejemplo: Caso homogéneo

Imagen de ejemplo: Caso heterogéneo

Gracias