6.1 Identificación de la estructura

Nuestro siguiente paso es describir la heurística de búsqueda para minimizar el score, que en lo que sigue suponemos que es el AIC.

Hay dos decisiones de diseño para decidir el algoritmo de aprendizaje de estructura:

Técnicas de busqueda.
* Hill-climbing
* Recocido simulado
* Algoritmos genéticos

Operadores de búsqueda Locales
* Agregar arista
* Eliminar arista
* Cambiar dirección de arista

Globales (ej. cambiar un nodo de lugar, más costoso)

Aunque hay varios algoritmos, uno utilizado comunmente es el de hill climbing.

6.1.1 Hill-climbing, escalada simple o ascenso de colina

  1. Iniciamos con una gráfica dada:
  • Gráfica vacía
  • Gráfica aleatoria
  • Conocimiento experto
  1. En cada iteración:
  • Consideramos el score para cada operador de búsqueda local (agregar, eliminar o cambiar la dirección de una arista)
  • Aplicamos el cambio que resulta en un mayor incremento en el score. Si tenemos empates elijo una operación al azar.
  1. Parar cuando ningún cambio resulte en mejoras del score.

Ejemplo. Eliminación de aristas Consideremos datos simulados de una red en forma de diamante:

set.seed(28)
n <- 600 # número de observaciones
a <- (rbinom(n, 1, 0.3)) # nodo raíz
b <- (rbinom(n, 1, a * 0.1 + (1 - a) * 0.8))
c <- (rbinom(n, 1, a * 0.2 + (1 - a) * 0.9))
d <- (rbinom(n, 1, b * c * 0.9 + (1 - b * c) * 0.1))
dat <- data.frame(a = factor(a), b = factor(b), c = factor(c), d = factor(d))
pander(head(dat), caption = "Muestra de la tabla de datos recien creada")
Muestra de la tabla de datos recien creada
a b c d
0 1 1 1
0 1 1 1
0 1 1 1
1 0 0 0
0 1 1 1
1 0 0 0

6.1.2 complejidad y score AIC

Una solución al problema de selección de modelos es usar una función de score por una que penalice la verosimilitud según el número de parámetros en el modelo. Una medida popular de este es el AIC (Aikaike information criterion). El AIC se define como

El score AIC de un modelo se define como \[AIC({\mathcal G}, \theta_{\mathcal G}) = -\frac{2}{N}loglik + \frac{2d}{N}=Dev + \frac{2d}{N},\] donde \(d\) es el número de parámetros en \(\theta_{\mathcal G}.\)


Nótemos que bajo este criterio, agregar una arista no necesariamente representa una mejora, pues aunque \(loglik\) no aumente o incluso disminuya, \(d\) definitivamente aumenta (añadimos variables). Es decir, el AIC es una combinación de una medida de ajuste del modelo con una penalización por complejidad, donde medimos complejidad por el número de parámetros del modelo.

¿Qué modelo debemos elegir de acuerdo al AIC?

No hay garantía de escoger el modelo óptimo usando el AIC (según la experiencia tiende a escoger modelos quizá demasiado complejos), pero es una guía útil para controlar complejidad.

Otra alternativa útil es el BIC, que también es un tipo de verosimilitud penalizada:

En la práctica se utilizan AIC y BIC. El AIC tiende a producir modelos más complejos con algún sobreajuste, mientras que el BIC tiende a producir modelos más simples y a veces con falta de ajuste. No hay acuerdo en qué medida es mejor en general, pero el balance se puede considerar como sigue: cuando es importante predecir o hacer inferencia, y no nos preocupa tanto obtener algunas aristas espurias, preferimos el AIC. Cuando buscamos la parte más robusta de la estructura de variación de las variables, aún cuando perdamos algunas dependencias débiles, puede ser mejor usar el BIC.

Ahora, comenzamos el proceso con una gráfica vacía:

# Esta función produce una especie de "devianza basada en AIC", útil para comparar modelo.
aic_dev <- function(fit, data){
  -2 * AIC(fit, data = data) / nrow(data)
}
grafica_0 <- empty.graph(c('a','b','c','d'))
fit_0 <- bn.fit(grafica_0, dat)
logLik(fit_0, data = dat)
## [1] -1543.696
AIC(fit_0, data = dat) # cuatro parámetros
## [1] -1547.696
aic_dev(fit_0, data = dat)
## [1] 5.158988

Consideramos agregar \(a\to d\), la nueva arista que mejora el AIC, y escogemos este cambio. Notemos que esta arista no existe en el modelo que genera los datos,

grafica_1 <- grafica_0
arcs(grafica_1) <- matrix(c('a', 'd'), ncol = 2, byrow = T)
fit_1 <- bn.fit(grafica_1, dat)
logLik(fit_1, data = dat)
## [1] -1458.222
aic_dev(fit_1, data = dat) 
## [1] 4.877407
graphviz.plot(grafica_1)
## Loading required namespace: Rgraphviz

Ahora agregamos \(a\to b\), que también mejora el AIC:

grafica_2 <- grafica_0
arcs(grafica_2) <- matrix(c('a','d','a','b'), ncol = 2, byrow = T)
fit_2 <- bn.fit(grafica_2, dat)
logLik(fit_2, data = dat)
## [1] -1288.1
aic_dev(fit_2, data = dat) 
## [1] 4.313667
graphviz.plot(grafica_2)

Igualmente, agregar \(a\to c\) merjoar el AIC:

grafica_3 <- grafica_0
arcs(grafica_3) <- matrix(c('a','d','a','b','a','c'), ncol = 2, byrow = T)
fit_3 <- bn.fit(grafica_3, dat)
logLik(fit_3, data = dat )
## [1] -1141.835
aic_dev(fit_3, data = dat) 
## [1] 3.829451
graphviz.plot(grafica_3)

Agregamos \(b\to d\) y \(c\to d\):

grafica_4 <- grafica_0
arcs(grafica_4) <- matrix(c('a','d','a','b','a','c','b','d'), ncol = 2, 
  byrow = T)
fit_4 <- bn.fit(grafica_4, dat)
logLik(fit_4, data = dat )
## [1] -1064.735
aic_dev(fit_4, data = dat) 
## [1] 3.579116
grafica_4 <- grafica_0
arcs(grafica_4) <- matrix(c('a','d','a','b','a','c','b','d','c','d'), ncol = 2, 
  byrow = T)
fit_4 <- bn.fit(grafica_4, dat)
logLik(fit_4, data = dat )
## [1] -1015.777
aic_dev(fit_4, data = dat) 
## [1] 3.429258
graphviz.plot(grafica_4)

Ahora nótese que podemos eliminar \(a\to d\), y mejoramos el AIC:

grafica_5 <- grafica_0
arcs(grafica_5) <- matrix(c('a','b','a','c','b','d','c','d'), ncol = 2, 
  byrow = T)
fit_5 <- bn.fit(grafica_5, dat)
logLik(fit_5, data = dat )
## [1] -1016.166
aic_dev(fit_5, data = dat) 
## [1] 3.417219
graphviz.plot(grafica_5)

Este última gráfica es el modelo original. La eliminación de arcos nos permitió recuperar el modelo original a pesar de nuestra decisión inicial temporalmente incorrecta de agregar \(a\to d\).

El algoritmo de ascenso de colina como está implementado en bn.learn resulta en:

graf_hc <- hc(dat, score='aic')
graphviz.plot(graf_hc)

Ejemplo: Cambios de dirección

Consideramos un ejemplo simple con un colisionador:

set.seed(28)
n <- 600
b <- (rbinom(n, 1, 0.4))
c <- (rbinom(n, 1, 0.7))
d <- (rbinom(n, 1, b*c*0.9+ (1-b*c)*0.1 ))
dat <- data.frame(factor(b),factor(c),factor(d))
names(dat) <- c('b','c','d')

Supongamos que comenzamos agregando la arista \(d\to c\) (sentido incorrecto).

grafica_0 <- empty.graph(c('b','c','d'))
arcs(grafica_0) <- matrix(c('d','c'), ncol=2, byrow=T)
graphviz.plot(grafica_0)

En el primer paso, agregamos \(b \to d\), que muestra una mejora grande:

graf_x <- hc(dat, start= grafica_0, score='aic', max.iter=1)
graphviz.plot(graf_x)

Pero en el siguiente paso nos damos cuenta que podemos mejorar considerablemente si construimos el modelo local de \(d\) a partir no sólo de \(b\) sino también de \(c\), y cambiamos dirección:

graf_x <- hc(dat, start= grafica_0, score='aic', max.iter=2)
graphviz.plot(graf_x)

Podemos examinar cada paso del algoritmo:

hc(dat, start = grafica_0, score='aic', debug=T)
## ----------------------------------------------------------------
## * starting from the following network:
## 
##   Random/Generated Bayesian network
## 
##   model:
##    [b][d][c|d] 
##   nodes:                                 3 
##   arcs:                                  1 
##     undirected arcs:                     0 
##     directed arcs:                       1 
##   average markov blanket size:           0.67 
##   average neighbourhood size:            0.67 
##   average branching factor:              0.33 
## 
##   generation algorithm:                  Empty 
## 
## * current score: -1105.878 
## * whitelisted arcs are:
## * blacklisted arcs are:
## * caching score delta for arc b -> c (51.958371).
## * caching score delta for arc b -> d (127.579833).
## * caching score delta for arc c -> b (-0.955426).
## * caching score delta for arc c -> d (25.648932).
## * caching score delta for arc d -> c (-25.648932).
## ----------------------------------------------------------------
## * trying to add one of 4 arcs.
##   > trying to add b -> c.
##     > delta between scores for nodes b c is 51.958371.
##     @ adding b -> c.
##   > trying to add b -> d.
##     > delta between scores for nodes b d is 127.579833.
##     @ adding b -> d.
##   > trying to add c -> b.
##     > delta between scores for nodes c b is -0.955426.
##   > trying to add d -> b.
##     > delta between scores for nodes d b is 127.579833.
## ----------------------------------------------------------------
## * trying to remove one of 1 arcs.
##   > trying to remove d -> c.
##     > delta between scores for nodes d c is -25.648932.
## ----------------------------------------------------------------
## * trying to reverse one of 1 arcs.
##   > trying to reverse d -> c.
##     > delta between scores for nodes d c is 0.000000.
## ----------------------------------------------------------------
## * best operation was: adding b -> d .
## * current network is :
## 
##   Bayesian network learned via Score-based methods
## 
##   model:
##    [b][d|b][c|d] 
##   nodes:                                 3 
##   arcs:                                  2 
##     undirected arcs:                     0 
##     directed arcs:                       2 
##   average markov blanket size:           1.33 
##   average neighbourhood size:            1.33 
##   average branching factor:              0.67 
## 
##   learning algorithm:                    Hill-Climbing 
##   score:                                 AIC (disc.) 
##   penalization coefficient:              1 
##   tests used in the learning procedure:  5 
##   optimized:                             TRUE 
## 
## * current score: -978.2983 
## * caching score delta for arc b -> d (-127.579833).
## * caching score delta for arc c -> d (78.562729).
## ----------------------------------------------------------------
## * trying to add one of 2 arcs.
##   > trying to add b -> c.
##     > delta between scores for nodes b c is 51.958371.
##     @ adding b -> c.
##   > trying to add c -> b.
##     > delta between scores for nodes c b is -0.955426.
## ----------------------------------------------------------------
## * trying to remove one of 2 arcs.
##   > trying to remove b -> d.
##     > delta between scores for nodes b d is -127.579833.
##   > trying to remove d -> c.
##     > delta between scores for nodes d c is -25.648932.
## ----------------------------------------------------------------
## * trying to reverse one of 2 arcs.
##   > trying to reverse b -> d.
##     > delta between scores for nodes b d is 0.000000.
##   > trying to reverse d -> c.
##     > delta between scores for nodes d c is 52.913797.
##     @ reversing d -> c.
## ----------------------------------------------------------------
## * best operation was: reversing d -> c .
## * current network is :
## 
##   Bayesian network learned via Score-based methods
## 
##   model:
##    [b][c][d|b:c] 
##   nodes:                                 3 
##   arcs:                                  2 
##     undirected arcs:                     0 
##     directed arcs:                       2 
##   average markov blanket size:           2.00 
##   average neighbourhood size:            1.33 
##   average branching factor:              0.67 
## 
##   learning algorithm:                    Hill-Climbing 
##   score:                                 AIC (disc.) 
##   penalization coefficient:              1 
##   tests used in the learning procedure:  7 
##   optimized:                             TRUE 
## 
## * current score: -925.3845 
## * caching score delta for arc b -> c (-0.955426).
## * caching score delta for arc b -> d (-180.493630).
## * caching score delta for arc c -> d (-78.562729).
## * caching score delta for arc d -> c (25.648932).
## ----------------------------------------------------------------
## * trying to add one of 2 arcs.
##   > trying to add b -> c.
##     > delta between scores for nodes b c is -0.955426.
##   > trying to add c -> b.
##     > delta between scores for nodes c b is -0.955426.
## ----------------------------------------------------------------
## * trying to remove one of 2 arcs.
##   > trying to remove b -> d.
##     > delta between scores for nodes b d is -180.493630.
##   > trying to remove c -> d.
##     > delta between scores for nodes c d is -78.562729.
## ----------------------------------------------------------------
## * trying to reverse one of 2 arcs.
##   > trying to reverse b -> d.
##     > delta between scores for nodes b d is -52.913797.
##   > trying to reverse c -> d.
##     > delta between scores for nodes c d is -52.913797.

Ejemplo simulado.

Comenzamos con una muestra relativamente chica, y utilizamos el BIC:

set.seed(280572)
n <- 300
a <- (rbinom(n, 1, 0.2))
b <- (rbinom(n, 1, a*0.1+(1-a)*0.8))
c <- (rbinom(n, 1, a*0.2+(1-a)*0.9))
d <- (rbinom(n, 1, b*c*0.9+ (1-b*c)*0.1 ))
e <- rbinom(n, 1, 0.4)
f <- rbinom(n, 1, e*0.3+(1-e)*0.6)
g <- rbinom(n, 1, f*0.2+(1-f)*0.8)
dat <- data.frame(factor(a),factor(b),factor(c),factor(d), factor(e), factor(f),
  factor(g))
names(dat) <- c('a','b','c','d','e','f','g')
grafica.1 <- hc(dat, score='bic')
graphviz.plot(grafica.1)

set.seed(280572)
n <- 300
a <- (rbinom(n, 1, 0.3))
b <- (rbinom(n, 1, a*0.1+(1-a)*0.8))
c <- (rbinom(n, 1, a*0.2+(1-a)*0.9))
d <- (rbinom(n, 1, b*c*0.9+ (1-b*c)*0.1 ))
e <- rbinom(n, 1, 0.4)
f <- rbinom(n, 1, e*0.3+(1-e)*0.6)
g <- rbinom(n, 1, f*0.2+(1-f)*0.8)
dat <- data.frame(factor(a),factor(b),factor(c),factor(d), factor(e), factor(f),
  factor(g))
names(dat) <- c('a','b','c','d','e','f','g')
grafica.1 <- hc(dat, score='aic')
graphviz.plot(grafica.1)

En este ejemplo, con el AIC obtenemos algunas aristas espurias, que en todo caso muestran relaciones aparentes débiles en los datos de entrenamiento. Nótese que AIC captura las relaciones importantes, y erra en cautela en cuanto a qué independencias están presentes en los datos.

6.1.3 Incorporando información acerca de la estructura

En algunos casos, tenemos información adicional de las posibles estructuras gráficas que son aceptables o deseables en los modelos que buscamos ajustar.

Esta información es muy valiosa cuando tenemos pocos datos o muchas variables (incluso puede ser crucial para obtener un modelo de buena calidad), y puede incorporarse en prohibiciones acerca de qué estructuras puede explorar el algoritmo.

Consideremos nuestro ejemplo anterior con considerablemente menos datos:

set.seed(28)
n <- 100
a <- (rbinom(n, 1, 0.2))
b <- (rbinom(n, 1, a*0.1+(1-a)*0.8))
c <- (rbinom(n, 1, a*0.2+(1-a)*0.9))
d <- (rbinom(n, 1, b*c*0.9+ (1-b*c)*0.1 ))
e <- rbinom(n, 1, 0.4)
f <- rbinom(n, 1, e*0.3+(1-e)*0.6)
g <- rbinom(n, 1, f*0.2+(1-f)*0.8)
dat <- data.frame(factor(a),factor(b),factor(c),factor(d), factor(e), factor(f),
  factor(g))
names(dat) <- c('a','b','c','d','e','f','g')
grafica.1 <- hc(dat, score='bic')
graphviz.plot(grafica.1)

Nótese que en este ejemplo BIC falla en identificar una dependencia, y afirma que hay una independencia condicional entre a y d dado c. AIC sin embargo captura la dependencia con un modelo demasiado complejo (tres flechas espurias):

grafica.1 <- hc(dat, score='aic')
graphviz.plot(grafica.1)

Sin embargo, si sabemos, por ejemplo, que no debe haber una flecha de c a f, y tiene que haber una de a a c, podemos mejorar nuestros modelos:

b.list <- data.frame(from=c('c','f'), to=c('f','c'))
w.list <- data.frame(from=c('a'), to=c('c'))
grafica.bic <- hc(dat, score='bic', blacklist=b.list, whitelist=w.list)
graphviz.plot(grafica.bic)

grafica.aic <- hc(dat, score='aic', blacklist=b.list, whitelist=w.list)
graphviz.plot(grafica.aic)

En este ejemplo estamos seguros de las aristas que forzamos. Muchas veces este no es el caso, y debemos tener cuidado:

  • Forzar la inclusión de una arista cuando esto no es necesario puede resultar en modelos demasiado complicados que incluyen estructuras espurias.

  • Exclusión de muchas aristas puede provocar también modelos que ajustan mal y no explican los datos.
set.seed(28)
n <- 600
b <- (rbinom(n, 1, 0.4))
c <- (rbinom(n, 1, 0.7))
d <- (rbinom(n, 1, b*c*0.9+ (1-b*c)*0.1 ))
dat.x <- data.frame(factor(b),factor(c),factor(d))
names(dat.x) <- c('b','c','d')

Supongamos que comenzamos agregando la arista \(d\to b\) (sentido incorrecto).

graphviz.plot(hc(dat.x, score='bic', whitelist=data.frame(from=c('d'), to=c('b'))))

Y no aprendimos nada, pues cualquier conjunta se factoriza de esta manera.

6.1.4 Sentido de las aristas

Los métodos de score a lo más que pueden aspirar es a capturar la clase de equivalencia Markoviana de la conjunta que nos interesa (es decir, gráficas que tienen las mismas independencias, y que cubren a exactamente las mismas conjuntas que se factorizan sobre ellas). Esto implica que hay cierta arbitrariedad en la selección de algunas flechas.

En la siguiente gráfica, por ejemplo, ¿qué pasa si cambiamos el sentido de la flecha entre e y f?

set.seed(28)
n <- 500
a <- (rbinom(n, 1, 0.2))
b <- (rbinom(n, 1, a*0.1+(1-a)*0.8))
c <- (rbinom(n, 1, a*0.2+(1-a)*0.9))
d <- (rbinom(n, 1, b*c*0.9+ (1-b*c)*0.1 ))
e <- rbinom(n, 1, 0.4)
f <- rbinom(n, 1, e*0.3+(1-e)*0.6)
g <- rbinom(n, 1, f*0.2+(1-f)*0.8)
dat <- data.frame(factor(a),factor(b),factor(c),factor(d), factor(e), factor(f),
  factor(g))
names(dat) <- c('a','b','c','d','e','f','g')
grafica.bic <- hc(dat, score='bic')
graphviz.plot(grafica.bic)

arcos <- grafica.bic$arcs
arcos
##      from to 
## [1,] "b"  "d"
## [2,] "a"  "b"
## [3,] "f"  "g"
## [4,] "a"  "c"
## [5,] "c"  "d"
## [6,] "e"  "f"
arcos[3,] <- c('g','f')
arcos[6,] <- c('f','e')
grafica.2 <- grafica.bic
arcs(grafica.2) <- arcos
graphviz.plot(grafica.2)

graphviz.plot(grafica.bic)

Vemos que no cambia la log-verosimilitud, ni ninguno de nuestros scores.

logLik(grafica.bic, data=dat)
## [1] -1648.022
logLik(grafica.2, data=dat)
## [1] -1648.022
BIC(grafica.bic, data=dat)
## [1] -1691.525
BIC(grafica.2, data=dat)
## [1] -1691.525
AIC(grafica.bic, data=dat)
## [1] -1662.022
AIC(grafica.2, data=dat)
## [1] -1662.022

Esto implica que la dirección de estas flechas no puede determinarse solamente usando los datos. Podemos seleccionar la dirección de estas flechas por otras consideraciones, como explicaciones causales, temporales, o de interpretación. Los modelos son equivalentes, pero tienen una parametrización destinta.

Mostrar que cambiar el sentido de una flecha que colisiona en \(d\) (que es un colisionador no protegido) no da scores equivalentes.

6.1.5 Variaciones de Hill-climbing

¿Cuál(es) de las siguientes opciones puede ser un problema para aprender la estructura de la red?
a. Máximos locales.
b. Pasos discretos en los scores cuando se perturba la estructura.
c. Eliminar un arco no se puede expresar como una operación atómica en la estructura.
d. Perturbaciones chicas en la estructura de la gráfica producen cambios muy chicos o nulos en el score (plateaux).

¿Por que consideramos el operador de cambiar dirección como candidato en cada iteración si es el resultado de elminar un arco y añadir un arco? Eliminar un arco en hill-climbing tiende a disminuir el score de tal manera que el paso inicial de eliminar el arco no se tomará.

Revertir la dirección es una manera de evitar máximos locales.

Algunas modificaciones de hill-climbing consisten en incluir estrategias:

  • Inicios aleatorios: Si estamos en una llanura, tomamos un número de pasos aleatorios para intentar encontrar una nueva pendiente y entonces comenzar a escalar nuevamente.

  • Tabu: Guardar una lista de los k pasos más recientes y la búsqueda no puede revertir estos pasos.