Modelos gaussianos no dirigidos
Distribución normal multivariada
Recordamos primero la distribución normal bivariada. Sean \((X_1,X_2)\) variables aleatorias continuas. Decimos que estas variables tienen distribución normal bivariada con media \((\mu_1,\mu_2)\) y matriz de varianza y covarianza
\[\Sigma=\left( \begin{array}{cc} \sigma_1^2 & \rho\sigma_1\sigma_2\\ \rho\sigma_1\sigma_2 & \sigma_2^2 \end{array} \right) \]
La función de densidad de \((X_1,X_2)\) está dada por \[f(x_1,x_2)\propto \exp \left\{- \frac{1}{2}(x_1-\mu_1, x_2-\mu_2)\Sigma^{-1} (x_1-\mu_1, x_2-\mu_2)^t \right\}.\]
Recordamos las siguientes propiedades:
- La marginal de \(X_i\) es normal con media \(\mu_i\) y varianza \(\sigma_1^2\).
- La condicional de \(X_2\) dado \(X_1=x_1\) es normal con media dada por la recta de regresión: \[\mu_2+\rho\frac{\sigma_2}{\sigma_1}(x_1-\mu_1)\] y varianza \((1-\rho^2)\sigma_2^2\).
- Si \(\rho=0\), entonces \(X_1\) y \(X_2\) son independentes.
- La dependencia de \(X_2\) de \(X_1\) sólo se da a través de la media condicional (la varianza de la condicional es fija).
Ahora consideramos \((X_1,\ldots, X_k)\) normal multivariada con media \(\mu\) y matriz de varianza y covarianza \(\Sigma\). \[f(x_1,\ldots,x_k)\propto \exp \left\{ -\frac{1}{2}(x-\mu)^t\Sigma^{-1} (x-\mu) \right\}.\] Si escribimos \(D=\Sigma^{-1}\) (a \(D\) se le llama a veces la matriz de precisión, que es simétrica positiva definida), entonces podemos también escribir la densidad conjunta como (agrupando términos y redefiniendo la constante de proporcionalidad):
\[f(x_1,\ldots,x_k)\propto \exp \left\{ -\frac{1}{2}\sum_{i,j}D_{ij}x_ix_j -\sum_{i} \gamma_i x_i \right\}.\]
Notamos ahora que esta conjunta se expresa naturalmente como un producto de factores donde cada factor depende de dos variables, notamos también que los factores que aparecen dependen de qué coeficientes son iguales a cero de la matriz de precisión \(D\). Según la teoría que hemos visto, si el coeficiente de un factor es cero esto implica que las dos variables correspondientes a dicho factor deben ser condicionalmente independientes dada el resto.
\small{ Estandarizamos las variables y examinamos las correlaciones:
library(bnlearn)## Warning: package 'bnlearn' was built under R version 3.1.2library(gRim)## Warning: package 'gRim' was built under R version 3.1.2## Loading required package: gRbase## Warning: package 'gRbase' was built under R version 3.1.2##
## Attaching package: 'gRbase'
##
## The following objects are masked from 'package:bnlearn':
##
## children, parents
##
## Loading required package: gRain## Warning: package 'gRain' was built under R version 3.1.2data(marks)
head(marks)## MECH VECT ALG ANL STAT
## 1 77 82 67 67 81
## 2 63 78 80 70 81
## 3 75 73 71 66 81
## 4 55 72 63 70 68
## 5 63 63 65 70 63
## 6 53 61 72 64 73marks_s <- data.frame(scale(marks, center=TRUE, scale=TRUE))
round(var(marks_s), 2)## MECH VECT ALG ANL STAT
## MECH 1.00 0.55 0.55 0.41 0.39
## VECT 0.55 1.00 0.61 0.49 0.44
## ALG 0.55 0.61 1.00 0.71 0.66
## ANL 0.41 0.49 0.71 1.00 0.61
## STAT 0.39 0.44 0.66 0.61 1.00Cuya matriz de precisión es
D <- solve(var(marks_s))
round(D, 2)## MECH VECT ALG ANL STAT
## MECH 1.60 -0.56 -0.51 0.00 -0.04
## VECT -0.56 1.80 -0.66 -0.15 -0.04
## ALG -0.51 -0.66 3.04 -1.11 -0.86
## ANL 0.00 -0.15 -1.11 2.18 -0.52
## STAT -0.04 -0.04 -0.86 -0.52 1.92Nótese que existen varios coeficientes cercanos a cero.
Más de la matriz de precisión
Los elementos de la diagonal de \(D=\Sigma^{-1}\) pueden interpretarse de la siguiente manera:
- Los elementos de la diagonal están relacionados con la \(R^2\) de la regresión de la variable correspondiente con respecto a todas las demás: \[D_{ii}= 1/(1-R_i^2)\]
fit_alg <- lm(ALG ~ MECH + VECT + ANL + STAT, data = marks_s)
# summary(fit_alg)
R2 <- summary(fit_alg)$r.squared
D_alg <- 1 / (1 - R2)
fit_anl <- lm(ANL ~ MECH + VECT + ALG + STAT, data = marks_s)
R2 <- summary(fit_anl)$r.squared
D_anl <- 1 / (1 - R2)
c(D_alg, D_anl)## [1] 3.0 2.2- Los elementos fuera de la diagonal están relacionados con la correlación parcial entre dos variables. La correlación parcial \(\rho_{ij|resto}\) entre \(X_i\) y \(X_j\) es la correlación entre los residuales de la regresión de \(X_i\) contra el resto de las variables y los residuales de \(X_j\) contra el resto de las variables. Se puede interpretar como la correlación que existe entre \(X_i\) y \(X_j\) cuando mantenemos el resto de las variables fijas. Cuando este coeficiente es mucho más chico que la correlación usual entre \(X_i\) y \(X_j\) quiere decir que una buena parte de la correlación de estas variables puede explicarse debido a variación entre el resto de las variables. Específicamente, \[\rho_{ij|resto}-\frac{D_{ij}}{\sqrt{D_{ii}D_{jj}}}.\]
fit_alg_2 <- lm(ALG ~ MECH + VECT + STAT, data = marks_s)
fit_anl_2 <- lm(ANL ~ MECH + VECT + STAT, data = marks_s)
cor(residuals(fit_alg_2), residuals(fit_anl_2))## [1] 0.43- D[3, 4] / sqrt(D[3, 3] * D[4, 4])## [1] 0.43Notemos que no es lo mismo que la correlación usual entre STAT y MECH.
cor(marks_s)[3, 4]## [1] 0.71Podemos examinar todas las correlaciones parciales haciendo:
D <- solve(cor(marks_s))
parciales_D <- -(t(D * (1/sqrt(diag(D)))) * (1 / sqrt(diag(D))))
round(parciales_D, 2)## MECH VECT ALG ANL STAT
## MECH -1.00 0.33 0.23 0.00 0.02
## VECT 0.33 -1.00 0.28 0.08 0.02
## ALG 0.23 0.28 -1.00 0.43 0.36
## ANL 0.00 0.08 0.43 -1.00 0.25
## STAT 0.02 0.02 0.36 0.25 -1.00
En cada una de las imágenes de abajo indica si la estructura de la gráfica del lado derecho corresponde a la matriz de correlación parcial (cuadro indica valor distinto a cero y vacío indica cero).
a. Arriba izquierda: Verdaderob. Arriba derecha: Verdadero
c. Abajo izquierda: Verdadero
d. Abajo derecha: Verdadero
Independencia condicional para la normal multivariada
De la discusión anterior, vemos que cuando \((X_1,\ldots, X_k)\) es normal multivariada, entonces \(X_i\) es independiente de \(X_j\) dado el resto de las variables si y sólo si la correlación parcial es cero, o de otra forma, cuando \(D_{ij}=0\).
Podemos ahora apelar a nuestros resultados anteriores (propiedades de Markov por pares, global y distribución de Gibbs) para demostrar el siguiente teorema (aunque también se puede resolver mediante cálculo):
round(cov(marks))## MECH VECT ALG ANL STAT
## MECH 306 127 102 106 117
## VECT 127 173 85 95 99
## ALG 102 85 113 112 122
## ANL 106 95 112 220 156
## STAT 117 99 122 156 298Si estandarizamos las variables para que la variabilidad sea comparable, obtenemos la matriz
cor(marks)## MECH VECT ALG ANL STAT
## MECH 1.00 0.55 0.55 0.41 0.39
## VECT 0.55 1.00 0.61 0.49 0.44
## ALG 0.55 0.61 1.00 0.71 0.66
## ANL 0.41 0.49 0.71 1.00 0.61
## STAT 0.39 0.44 0.66 0.61 1.00Cuya inversa es
inv_1 <- round(solve(cor(marks)), 2)
inv_1## MECH VECT ALG ANL STAT
## MECH 1.60 -0.56 -0.51 0.00 -0.04
## VECT -0.56 1.80 -0.66 -0.15 -0.04
## ALG -0.51 -0.66 3.04 -1.11 -0.86
## ANL 0.00 -0.15 -1.11 2.18 -0.52
## STAT -0.04 -0.04 -0.86 -0.52 1.92Donde vemos los coeficientes de esta inversa para los pares ANL-MECH, ANL-VECT, VEC-STAT, MECH-STAT son cercanos a cero. Esto implica que son variables cercanas a la independencia dado el resto de las variables. Así que un modelo apropiado para estos datos tiene todas las aristas, excepto la lista que mencionó arriba. Terminamos entonces con la gráfica
library(Rgraphviz)## Loading required package: graph## Warning: package 'graph' was built under R version 3.1.2##
## Attaching package: 'graph'
##
## The following objects are masked from 'package:bnlearn':
##
## degree, nodes, nodes<-data(marks)
ug_1 <- ug(~ ANL:STAT + ALG:STAT + ALG:ANL + MECH:ALG + VECT:ALG + VECT:MECH)
plot(ug_1)
Podemos graficar mostrando las correlaciones parciales:
nombres <- sapply(edgeList(ug_1), function(x){
nombre <- paste0('', paste(c(x[1], x[2]), collapse = "~"), '')
nombre
})
peso <- sapply(edgeList(ug_1), function(x){
peso <- as.character(round(parciales_D[x[1], x[2]], 2))
})
etiquetas <- peso
names(etiquetas) <- nombres
Rgraphviz::plot(ug_1, edgeAttrs = list(label = etiquetas))
Estimación de máxima verosimilitud con estructura conocida
* Si no hay una arista entre \(X_i\) y \(X_j\), entonces \(D_{ij}=D_{ji}=0\) donde \(D=\Sigma^{-1}\).
De esta forma, vemos que cuando tenemos una estructura no dirigida dada, nuestro trabajo es estimar la matriz de covarianza (o la precisión) con restricciones de ceros dadas por las independencias condicionales.
Podemos entonces maximizar la verosimilitud con las restricciones sobre la inversa de la matriz de varianza y covarianza. En general no es problema trivial.
Aprendizaje de estructura
Igual que en modelos log-lineales para datos categóricos, podemos ajustar modelos hacia atrás y hacia adelante (agregando/eliminando aristas) usando el BIC/AIC como criterio de selección, con distintas restricciones sobre la matriz de precisión.
Ejemplo En este ejemplo consideramos las estaturas de los dos padres y de dos de sus hijas seleccionados al azar (entre las familias que tienen al menos dos hijas).
library(plyr)##
## Attaching package: 'plyr'
##
## The following object is masked from 'package:graph':
##
## joinlibrary(gRim)
library(HistData)## Warning: package 'HistData' was built under R version 3.1.2data(GaltonFamilies)
dat <- GaltonFamilies
dat.2 <- ddply(dat, c('family'), function(df){
df$id_child <- sample(1:nrow(df), nrow(df))
df
})
dat.female <- subset(dat.2, gender=='female')
dat.female <- ddply(dat.female, c('family'), transform,
num.female=length(childHeight))
dat.female.2 <- subset(dat.female, num.female>1)
set.seed(28)
dat.3 <- ddply(dat.female.2, c('family'), function(df){
df$id_child <- sample(1:nrow(df), nrow(df))
subset(df, id_child<3)
})
library(reshape2)## Warning: package 'reshape2' was built under R version 3.1.2dat.4 <- dcast(dat.3, family+father+mother~ id_child, value.var='childHeight')
names(dat.4)[c(4,5)] <- c('h1','h2')
sigma <- var(dat.4[, c('father','mother', 'h1','h2')])
sigma## father mother h1 h2
## father 7.17 -0.19 2.3 2.0
## mother -0.19 4.66 1.2 1.1
## h1 2.29 1.23 5.5 1.5
## h2 2.04 1.11 1.5 4.9solve((sigma))## father mother h1 h2
## father 0.180 0.040 -0.066 -0.064
## mother 0.040 0.244 -0.056 -0.055
## h1 -0.066 -0.056 0.229 -0.029
## h2 -0.064 -0.055 -0.029 0.251mod.1 <- cmod(~.^1, data=dat.4[,c('father','mother', 'h1','h2')])
plot(mod_hijas <- forward(mod.1, type='unrestricted'))## . FORWARD: type=unrestricted search=all, criterion=aic(2.00), alpha=0.00
## . Initial model: is graphical=TRUE is decomposable=TRUE
## change.AIC -14.4621 Edge added: father h1
## change.AIC -12.5977 Edge added: father h2
## change.AIC -4.9947 Edge added: h1 mother
## change.AIC -7.9502 Edge added: h2 mother
## change.AIC -8.2410 Edge added: father mother
Podemos examinar la matriz de precisión ajustada:
D <- mod_hijas$fitinfo$K
D## father mother h1 h2
## father 0.186 0.044 -0.074 -0.073
## mother 0.044 0.250 -0.063 -0.062
## h1 -0.074 -0.063 0.228 0.000
## h2 -0.073 -0.062 0.000 0.250y para checar el ajuste comparamos la matriz de covarianzas empírica con la ajustada:
round(cor(dat.4[, -1]), 2)## father mother h1 h2
## father 1.00 -0.03 0.36 0.34
## mother -0.03 1.00 0.24 0.23
## h1 0.36 0.24 1.00 0.28
## h2 0.34 0.23 0.28 1.00round(cov2cor(solve(D)), 2)## father mother h1 h2
## father 1.00 -0.03 0.36 0.34
## mother -0.03 1.00 0.24 0.23
## h1 0.36 0.24 1.00 0.19
## h2 0.34 0.23 0.19 1.00Podemos comparar con
mod.hijas.2 <- update(mod_hijas, list(dedge = ~mother:father))
plot(mod.hijas.2)
mod.hijas.2## Model: A cModel with 4 variables
## graphical : TRUE decomposable : FALSE
## -2logL : 2059.38 mdim : 8 aic : 2075.38
## ideviance : 43.03 idf : 4 bic : 2097.41
## deviance : 6.77 df : 2mod.hijas.3 <- update(mod_hijas, list(dedge=~father:h1))
plot(mod.hijas.3)
mod.hijas.3## Model: A cModel with 4 variables
## graphical : TRUE decomposable : TRUE
## -2logL : 2072.68 mdim : 8 aic : 2088.68
## ideviance : 29.73 idf : 4 bic : 2110.71
## deviance : 20.07 df : 2Nótese que la estimación de la matriz de varianza y covarianza está cercana a la empírica:
D <- mod.hijas.2$fitinfo$K
round(cor(dat.4[,-1]),2)## father mother h1 h2
## father 1.00 -0.03 0.36 0.34
## mother -0.03 1.00 0.24 0.23
## h1 0.36 0.24 1.00 0.28
## h2 0.34 0.23 0.28 1.00round(cov2cor(solve(D)),2)## father mother h1 h2
## father 1.00 0.14 0.36 0.34
## mother 0.14 1.00 0.24 0.23
## h1 0.36 0.24 1.00 0.16
## h2 0.34 0.23 0.16 1.00Mientras que eliminar father-h1 muestra un desajuste más considerable:
D <- mod.hijas.3$fitinfo$K
round(cor(dat.4[,-1]),2)## father mother h1 h2
## father 1.00 -0.03 0.36 0.34
## mother -0.03 1.00 0.24 0.23
## h1 0.36 0.24 1.00 0.28
## h2 0.34 0.23 0.28 1.00round(cov2cor(solve(D)),2)## father mother h1 h2
## father 1.00 -0.03 -0.01 0.34
## mother -0.03 1.00 0.24 0.23
## h1 -0.01 0.24 1.00 0.06
## h2 0.34 0.23 0.06 1.00Ejemplo: datos bodyfat
Consideramos varias mediciones de dimensiones corporales de un conjunto de hombres, más una medición adicional de grasa corporal. En este ejemplo ilustramos
bodyfat <- read.table('bodyfat.txt', header=T)
head(bodyfat)## densidad grasacorp edad peso estatura cuello pecho abdomen cadera muslo
## 1 1.1 12.3 23 154 68 36 93 85 94 59
## 2 1.1 6.1 22 173 72 38 94 83 99 59
## 3 1.0 25.3 22 154 66 34 96 88 99 60
## 4 1.1 10.4 26 185 72 37 102 86 101 60
## 5 1.0 28.7 24 184 71 34 97 100 102 63
## 6 1.1 20.9 24 210 75 39 104 94 108 66
## rodilla tobillo biceps antebrazo muñeca
## 1 37 22 32 27 17
## 2 37 23 30 29 18
## 3 39 24 29 25 17
## 4 37 23 32 29 18
## 5 42 24 32 28 18
## 6 42 26 36 31 19cov(bodyfat)## densidad grasacorp edad peso estatura cuello pecho abdomen
## densidad 0.00036 -0.16 -0.067 -0.33 0.0068 -0.022 -0.11 -0.16
## grasacorp -0.15732 70.04 30.738 150.62 -2.7433 9.980 49.57 73.40
## edad -0.06659 30.74 158.811 -4.72 -7.9230 3.477 18.75 31.31
## peso -0.33227 150.62 -4.721 863.72 33.1856 59.348 221.55 281.41
## estatura 0.00682 -2.74 -7.923 33.19 13.4165 2.259 4.17 3.47
## cuello -0.02188 9.98 3.477 59.35 2.2591 5.909 16.08 19.77
## pecho -0.10952 49.57 18.746 221.55 4.1654 16.084 71.07 83.25
## abdomen -0.16396 73.40 31.310 281.41 3.4683 19.766 83.25 116.27
## cadera -0.08308 37.48 -4.544 198.10 4.4713 12.799 50.09 67.52
## muslo -0.05526 24.59 -13.238 134.03 2.8544 8.879 32.30 43.40
## rodilla -0.02272 10.27 0.532 60.47 2.5270 3.942 14.63 19.17
## tobillo -0.00854 3.77 -2.244 30.57 1.6436 1.969 6.90 8.28
## biceps -0.02801 12.47 -1.567 71.07 2.2998 5.370 18.54 22.32
## antebrazo -0.01352 6.11 -2.166 37.43 1.6923 3.063 9.88 10.97
## muñeca -0.00579 2.71 2.512 20.02 1.1013 1.690 5.20 6.24
## cadera muslo rodilla tobillo biceps antebrazo muñeca
## densidad -0.083 -0.055 -0.023 -0.0085 -0.028 -0.014 -0.0058
## grasacorp 37.483 24.587 10.267 3.7725 12.472 6.111 2.7078
## edad -4.544 -13.238 0.532 -2.2439 -1.567 -2.166 2.5122
## peso 198.099 134.032 60.473 30.5686 71.071 37.431 20.0230
## estatura 4.471 2.854 2.527 1.6436 2.300 1.692 1.1013
## cuello 12.799 8.879 3.942 1.9690 5.370 3.063 1.6904
## pecho 50.094 32.303 14.629 6.9013 18.540 9.883 5.1959
## abdomen 67.522 43.399 19.172 8.2832 22.316 10.967 6.2398
## cadera 51.324 33.715 14.228 6.7801 16.001 7.890 4.2142
## muslo 33.715 27.562 10.119 4.8032 12.078 6.013 2.7383
## rodilla 14.228 10.119 5.817 2.5001 4.946 2.709 1.4962
## tobillo 6.780 4.803 2.500 2.8727 2.483 1.435 0.8959
## biceps 16.001 12.078 4.946 2.4828 9.128 4.141 1.7830
## antebrazo 7.890 6.013 2.709 1.4352 4.141 4.083 1.1047
## muñeca 4.214 2.738 1.496 0.8959 1.783 1.105 0.8716precision <- round(100*solve(cor(bodyfat[,-1])))
round(100*cov2pcor(cov(bodyfat)))## densidad grasacorp edad peso estatura cuello pecho abdomen
## densidad 100 -96 5 8 0 1 7 -14
## grasacorp -96 100 8 4 -2 -3 6 4
## edad 5 8 100 -19 -13 10 3 30
## peso 8 4 -19 100 44 27 40 23
## estatura 0 -2 -13 44 100 -3 -21 -6
## cuello 1 -3 10 27 -3 100 3 10
## pecho 7 6 3 40 -21 3 100 39
## abdomen -14 4 30 23 -6 10 39 100
## cadera 6 3 -7 52 -25 -18 -19 26
## muslo -6 -3 -37 7 -20 12 -15 4
## rodilla -1 0 25 25 9 -15 -9 -6
## tobillo -10 -8 -10 19 -6 -9 -2 -12
## biceps -9 -7 7 19 -5 9 9 -15
## antebrazo -1 4 -18 3 -1 13 12 -13
## muñeca 6 0 35 11 9 26 -2 0
## cadera muslo rodilla tobillo biceps antebrazo muñeca
## densidad 6 -6 -1 -10 -9 -1 6
## grasacorp 3 -3 0 -8 -7 4 0
## edad -7 -37 25 -10 7 -18 35
## peso 52 7 25 19 19 3 11
## estatura -25 -20 9 -6 -5 -1 9
## cuello -18 12 -15 -9 9 13 26
## pecho -19 -15 -9 -2 9 12 -2
## abdomen 26 4 -6 -12 -15 -13 0
## cadera 100 33 4 -4 -7 -8 0
## muslo 33 100 27 -2 25 -1 -4
## rodilla 4 27 100 17 -8 9 7
## tobillo -4 -2 17 100 -5 -3 25
## biceps -7 25 -8 -5 100 28 6
## antebrazo -8 -1 9 -3 28 100 20
## muñeca 0 -4 7 25 6 20 100bodyfat.1 <- bodyfat[,c('grasacorp','abdomen','peso','estatura','biceps','antebrazo',
'rodilla','muñeca')]
summary(bodyfat.1)## grasacorp abdomen peso estatura biceps
## Min. : 0 Min. : 69 Min. :118 Min. :30 Min. :25
## 1st Qu.:12 1st Qu.: 85 1st Qu.:159 1st Qu.:68 1st Qu.:30
## Median :19 Median : 91 Median :176 Median :70 Median :32
## Mean :19 Mean : 93 Mean :179 Mean :70 Mean :32
## 3rd Qu.:25 3rd Qu.: 99 3rd Qu.:197 3rd Qu.:72 3rd Qu.:34
## Max. :48 Max. :148 Max. :363 Max. :78 Max. :45
## antebrazo rodilla muñeca
## Min. :21 Min. :33 Min. :15.8
## 1st Qu.:27 1st Qu.:37 1st Qu.:17.6
## Median :29 Median :38 Median :18.3
## Mean :29 Mean :39 Mean :18.2
## 3rd Qu.:30 3rd Qu.:40 3rd Qu.:18.8
## Max. :35 Max. :49 Max. :21.4bodyfat.1 <- subset(bodyfat.1, estatura > 30 & peso <300)
pairs(bodyfat.1)
Usamos el criterio bic ajustando modelos hacia adelante obtenemos:
mod.1 <- cmod(~.^1, data=bodyfat.1)
plot(mod.2 <- forward(mod.1, k=log(nrow(bodyfat.1)), type='unrestricted'))## . FORWARD: type=unrestricted search=all, criterion=aic(5.52), alpha=0.00
## . Initial model: is graphical=TRUE is decomposable=TRUE
## change.AIC -354.8385 Edge added: abdomen peso
## change.AIC -304.1353 Edge added: peso rodilla
## change.AIC -278.1341 Edge added: abdomen grasacorp
## change.AIC -234.1223 Edge added: peso biceps
## change.AIC -180.9975 Edge added: peso muñeca
## change.AIC -163.5689 Edge added: biceps antebrazo
## change.AIC -70.7975 Edge added: peso estatura
## change.AIC -211.9549 Edge added: estatura grasacorp
## change.AIC -83.0803 Edge added: abdomen estatura
## change.AIC -46.2559 Edge added: muñeca antebrazo
## change.AIC -26.0306 Edge added: muñeca grasacorp
## change.AIC -11.4932 Edge added: peso antebrazo
## change.AIC -9.9454 Edge added: estatura rodilla
## change.AIC -8.7992 Edge added: biceps estatura
## change.AIC -26.5494 Edge added: biceps abdomen
## change.AIC -4.1848 Edge added: peso grasacorp
## change.AIC -4.3748 Edge added: abdomen antebrazo
## change.AIC -2.7400 Edge added: estatura antebrazo
## change.AIC -0.4849 Edge added: abdomen muñeca
## change.AIC -1.0522 Edge added: rodilla muñeca
Podemos checar el ajuste del modelo verificando que la matriz ajustada de correlaciones es similar a la observada:
D <- mod.2$fitinfo$K
round(cor(bodyfat.1[,]),2)## grasacorp abdomen peso estatura biceps antebrazo rodilla muñeca
## grasacorp 1.00 0.82 0.62 -0.03 0.48 0.36 0.49 0.34
## abdomen 0.82 1.00 0.87 0.19 0.66 0.53 0.71 0.60
## peso 0.62 0.87 1.00 0.51 0.79 0.68 0.84 0.73
## estatura -0.03 0.19 0.51 1.00 0.32 0.32 0.51 0.40
## biceps 0.48 0.66 0.79 0.32 1.00 0.70 0.65 0.61
## antebrazo 0.36 0.53 0.68 0.32 0.70 1.00 0.58 0.60
## rodilla 0.49 0.71 0.84 0.51 0.65 0.58 1.00 0.66
## muñeca 0.34 0.60 0.73 0.40 0.61 0.60 0.66 1.00round(cov2cor(solve(D)),2)## grasacorp abdomen peso estatura biceps antebrazo rodilla muñeca
## grasacorp 1.00 0.82 0.62 -0.03 0.45 0.33 0.47 0.34
## abdomen 0.82 1.00 0.87 0.19 0.66 0.53 0.71 0.60
## peso 0.62 0.87 1.00 0.51 0.79 0.68 0.84 0.73
## estatura -0.03 0.19 0.51 1.00 0.32 0.32 0.51 0.40
## biceps 0.45 0.66 0.79 0.32 1.00 0.70 0.66 0.60
## antebrazo 0.33 0.53 0.68 0.32 0.70 1.00 0.58 0.60
## rodilla 0.47 0.71 0.84 0.51 0.66 0.58 1.00 0.66
## muñeca 0.34 0.60 0.73 0.40 0.60 0.60 0.66 1.00Y vemos que en en efecto estas dos matrices son muy similares.
En muchos casos, estos modelos presentan algunas correlaciones débiles que, aunque nos pueden interesar, pueden hacer difícil entender los rasgos importantes de nuestros datos. Vale la pena examinar ajustes penalizando por complejidad más fuertemente, y viendo qué tanto cambia la calidad del ajuste. Por ejemplo, si incrementamos la penalización (el bic da una penalización de 5):
plot(mod.3 <- forward(mod.1, k=20, type='unrestricted'))## . FORWARD: type=unrestricted search=all, criterion=aic(20.00), alpha=0.00
## . Initial model: is graphical=TRUE is decomposable=TRUE
## change.AIC -340.3599 Edge added: abdomen peso
## change.AIC -289.6568 Edge added: peso rodilla
## change.AIC -263.6555 Edge added: abdomen grasacorp
## change.AIC -219.6437 Edge added: peso biceps
## change.AIC -166.5189 Edge added: peso muñeca
## change.AIC -149.0904 Edge added: biceps antebrazo
## change.AIC -56.3189 Edge added: peso estatura
## change.AIC -197.4764 Edge added: estatura grasacorp
## change.AIC -68.6018 Edge added: abdomen estatura
## change.AIC -31.7774 Edge added: muñeca antebrazo
## change.AIC -11.5520 Edge added: muñeca grasacorp
D <- mod.2$fitinfo$K
parciales.D <- -(t(D*(1/sqrt(diag(D))))*(1/sqrt(diag(D))))
graf.1 <- ugList(mod.3$glist)
nombres <- sapply(edgeList(graf.1), function(x){
nombre <- paste0('',paste(c(x[1],x[2]), collapse="~"),'')
nombre
})
peso <- sapply(edgeList(graf.1), function(x){
peso <- as.character(round(parciales.D[x[1],x[2]],2))
})
etiquetas <- peso
names(etiquetas) <- nombres
plot(graf.1, edgeAttrs=list(label=etiquetas))
Este modelo establece, por ejemplo, que dado peso, biceps y estatura son independientes:
library(Hmisc)
library(ggplot2)
bodyfat.2 <- bodyfat.1
bodyfat.2$grupo.peso <- cut2(bodyfat.2$peso, g=5)
ggplot(bodyfat.2, aes(x=biceps, y=estatura)) +
geom_point() +
facet_wrap(~grupo.peso, nrow = 1) +
geom_smooth()
Sin condicionar observamos,
ggplot(bodyfat.2, aes(x=biceps, y=estatura)) +
geom_point() + geom_smooth()
Comparando las matrices de correlaciones.
D <- mod.3$fitinfo$K
round(cor(bodyfat.1[,]),2)## grasacorp abdomen peso estatura biceps antebrazo rodilla muñeca
## grasacorp 1.00 0.82 0.62 -0.03 0.48 0.36 0.49 0.34
## abdomen 0.82 1.00 0.87 0.19 0.66 0.53 0.71 0.60
## peso 0.62 0.87 1.00 0.51 0.79 0.68 0.84 0.73
## estatura -0.03 0.19 0.51 1.00 0.32 0.32 0.51 0.40
## biceps 0.48 0.66 0.79 0.32 1.00 0.70 0.65 0.61
## antebrazo 0.36 0.53 0.68 0.32 0.70 1.00 0.58 0.60
## rodilla 0.49 0.71 0.84 0.51 0.65 0.58 1.00 0.66
## muñeca 0.34 0.60 0.73 0.40 0.61 0.60 0.66 1.00round(cov2cor(solve(D)),2)## grasacorp abdomen peso estatura biceps antebrazo rodilla muñeca
## grasacorp 1.00 0.82 0.64 -0.03 0.49 0.35 0.54 0.34
## abdomen 0.82 1.00 0.87 0.19 0.68 0.52 0.74 0.58
## peso 0.64 0.87 1.00 0.51 0.79 0.61 0.84 0.73
## estatura -0.03 0.19 0.51 1.00 0.41 0.34 0.43 0.45
## biceps 0.49 0.68 0.79 0.41 1.00 0.70 0.66 0.62
## antebrazo 0.35 0.52 0.61 0.34 0.70 1.00 0.52 0.60
## rodilla 0.54 0.74 0.84 0.43 0.66 0.52 1.00 0.61
## muñeca 0.34 0.58 0.73 0.45 0.62 0.60 0.61 1.00Donde vemos que este modelo más simple recupera razonablemente bien la estructura de covarianza.
Utiliza los datos carcass para crear un modelo gráfico Gaussiano. ¿Qué relaciones de independencia condicional lees en la gráfica?
Supuesto de normalidad multivariada
El supuesto de normalidad multivariada es en general difícil de verificar. Generalmente, se considera que la pregunta interesante es si los datos son suficientemente cercanos a normalidad para aplicar un procedimiento dado. Existen diagnósticos, pero aquí consideramos que estos modelos gaussianos, típicamente, pueden utilizarse más como procedimientos exploratorios.