library(metodosMultivariados2017) library(datasets) library(knitr)
data("HairEyeColor") tbl <- apply(HairEyeColor, c(1,2), sum) dimnames(tbl) <- list(paste(dimnames(tbl)[[1]], "Hair"), paste(dimnames(tbl)[[2]], "Eyes")) kable(tbl)
X <- contingency_dummies(tbl)[[1]] Y <- contingency_dummies(tbl)[[2]]
kable(head(X))
kable(tail(Y))
Recuperar la tabla de contingencia
kable(t(X) %*% Y)
Total de observaciones en X
kable(t(X) %*% X)
Total de observaciones en Y
kable(t(Y) %*% Y)
Matriz de frecuencias relativas
kable(t(X) %*% Y / sum(as.numeric(t(X) %*% Y)))
Otra forma de construirla matricialmente
ones_row <- rep(1, ncol(X)) ones_col <- rep(1, ncol(Y)) N <- as.numeric(ones_row %*% t(X) %*% Y %*% ones_col) P <- t(X) %*% Y / N kable(P)
Vamos a construir la matriz de perfiles por fila. Consiste del conteo en cada celda dividiendo cada fila por la suma de la fila. Recordemos que las suma de las filas se obtienen de la diagonal de $X^tX$.
Ademas reconocemos la estructura conocida de $$(X^tX)^{-1}X^tY$$. La entrada $(r,c)$ representa $P(Y=c|X=r)$.
# (X'X)^(-1)X'Y Pr <- solve(t(X) %*% X, t(X) %*% Y) kable(Pr)
Visualizamos
library(ggtern)
plot_dat <- data.frame(Pr, check.names = FALSE) ggtern(plot_dat, aes(x = `Brown Eyes`, y = `Green Eyes`, z = `Hazel Eyes`, colour = `Blue Eyes`)) + geom_point() + geom_text(aes(label = row.names(plot_dat)))
Para poder medir distancias adecuadamente entre cada perfil de fila, necesitamos medir la diferencia en distribuciĆ³n pesada inversamente proporcional al peso o masa de cada columna. Esto se conoce como distancia chi-cuadrada.
chi_dist <- function(row1, row2, tbl) { col_masses <- as.numeric(apply(tbl, 2, sum)) / sum(as.numeric(tbl)) Pr <- solve(t(X) %*% X, t(X) %*% Y) dif <- as.numeric(Pr[row1, ]) - as.numeric(Pr[row2, ]) sqrt(sum((1 / col_masses)*dif^2)) }
col_mass <- apply(tbl, 2, sum) / sum(as.numeric(tbl)) kable(col_mass)
Por ejemplo, distancia entre brown hair y black hair es
chi_dist(1, 2, tbl)
y la diferencia entre brown hair y blond hair es
chi_dist(1, 4, tbl)
Matriz de distancia entre perfiles de fila
hair_types <- dimnames(tbl)[[1]] n_types <- length(hair_types) row_dist_mat <- matrix(0, ncol = n_types, nrow = n_types, dimnames = list(hair_types, hair_types)) for (i in 1:(n_types - 1)) { for (j in i:n_types) { row_dist_mat[i, j] <- chi_dist(i, j, tbl) row_dist_mat[j, i] <- row_dist_mat[i, j] } } kable(row_dist_mat)
Tenemos que hacer una descomposicion de valores singulares generalizada como sigue.
Primero vamos a voltear a ver la matriz de diferencias respecto al "tipo de ojo esperado del total de la poblacion del mundo".
col_mass <- apply(tbl, 2, sum) / sum(as.numeric(tbl)) kable(col_mass)
Vamos a calcular las diferencias respecto a perfil del "mundo"
# restar a cada columna su masa (el porcentaje obsevado de ese color de ojos) Pr_tilde <- Pr - ones_col %*% t(col_mass) kable(Pr_tilde)
Notemos
kable(ones_col %*% t(col_mass))
apply(tbl, 2, sum) / sum(as.numeric(tbl))
n <- sum(as.numeric(tbl)) P <- tbl/n row_mass <- as.numeric(apply(P, 1, sum)) col_mass <- as.numeric(apply(P, 2, sum)) Dr <- diag(row_mass) Dc <- diag(col_mass) RowProfile <- solve(Dr, P) RowProfileC <- RowProfile for (i in 1:nrow(RowProfile)) RowProfileC[i, ] <- RowProfileC[i, ] -col_mass covR <- RowProfileC %*% solve(Dc) %*% t(RowProfileC) res <- eigen(covR) row_factor1 <- res$vectors[ ,1] * sqrt(res$values[1]) row_factor2 <- res$vectors[ ,2] * sqrt(res$values[2]) row_names <- dimnames(tbl)[[1]]
data <- data.frame(row_factor1, row_factor2, row_names) ggplot(data, aes(x = row_factor1, y= row_factor2, colour = row_names)) + geom_point(size = 2, alpha = 0.6) + geom_text(aes(x = row_factor1, y= row_factor2, label = row_names))
# varianza explicada var_exp <- cumsum(res$values) / sum(res$values) scales::percent(var_exp)
ColProfile <- P %*% solve(Dc) ColProfileC <- ColProfile for (j in 1:ncol(ColProfileC)) ColProfileC[ ,j] <- ColProfileC[ ,j] - row_mass covC <- t(ColProfileC) %*% solve(Dr) %*% ColProfileC res <- eigen(covC) col_factor1 <- res$vectors[ ,1] * sqrt(res$values[1]) col_factor2 <- res$vectors[ ,2] * sqrt(res$values[2]) col_names <- dimnames(tbl)[[2]]
data <- data.frame(col_factor1, col_factor2, col_names) ggplot(data, aes(x = col_factor1, y= col_factor2, colour = col_names)) + geom_point(size = 2, alpha = 0.6) + geom_text(aes(x = col_factor1, y= col_factor2, label = col_names))
data_row <- data.frame(factor1 = row_factor1, factor2 = row_factor2, variable = row_names, type = "hair") data_col <- data.frame(factor1 = - col_factor1, factor2 = col_factor2, variable = col_names, type = "eyes") data_all <- rbind(data_row, data_col) ggplot(data_all, aes(x = factor1, y= factor2, colour = variable, shape = type)) + geom_point(size = 2, alpha = 0.6) + geom_text(aes(x = factor1, y= factor2, label = variable))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.