Análisis canónico de poblaciones

El análisis canónico de poblaciones es una técnica de análisis multivariante que tiene el objetivo de representar varios grupos de individuos de forma óptima mediante unos ejes canónicos ortogonales. Eso se consigue de manera que la dispersión entre esos grupos sea máxima con relación a la dispersión dentro de cada grupo. Además, en esa representación, la distancia euclídea entre dos individuos coincide con la distancia de Mahalanobis entre ellos en las variables originales.
Se trata pues de una técnica de representación de datos en dimensión reducida, normalmente los dos primeros ejes, que puede acompañar gráficamente un MANOVA de un factor (la población).

Vamos a poner un ejemplo muy conocido. Se trata de los datos sobre cráneos de varones egipcios de cinco épocas históricas que se pueden obtener en el siguiente enlace:

http://lib.stat.cmu.edu/DASL/Datafiles/EgyptianSkulls.html

También podemos bajarlos desde la página del libro de Everitt (2005) y así los podremos cargar directamente en R con la siguientes instrucciones:


skulls <- source("/(path)/chap5skulls.dat")$value

str(skulls)

attach(skulls)


Donde el path debe ser la dirección a la carpeta donde hemos dejado el archivo una vez descomprimido. Así tendremos la base de datos skulls con cinco variables. La primera variable es el factor EPOCH y las otras cuatro son las medidas biométricas del cráneo estudiadas.

En primer lugar podemos realizar un MANOVA para contrastar la diferencia de medias entre los niveles del factor o poblaciones. No entraremos aquí en la comprobación de las hipótesis de normalidad y de igualdad de las matrices de covarianzas.

skulls.manova <- manova(cbind(MB,BH,BL,NH) ~ EPOCH)

summary(skulls.manova, test="Wilks") # test="Pillai" or "Hotelling" or "Roy"

El test rechaza la igualdad de medias y, por lo tanto, justifica el análisis canónico de poblaciones.

El siguiente paso es obtener el paquete candisc para un Canonical discriminant analysis ya que los ejes canónicos también sirven ese tipo de análisis discriminante.

library(candisc)

La función candisc realiza el análisis canónico discriminante, pero necesita como objeto principal un modelo lineal:

skulls.mod <- lm(cbind(MB,BH,BL,NH) ~ EPOCH)

Anova(skulls.mod, test="Wilks") # Manova

skulls.can1 <- candisc(skulls.mod, term="EPOCH")

plot(skulls.can1, conf=0.90, type="n")

El gráfico que se obtiene es el que vemos al principio de este artículo. Como corresponde a un ejemplo de manual, los dos primeros ejes canónicos representan muy bien a los datos y las poblaciones se separan de forma cronológica. También se representan las variables en una especie de biplot que explica mejor las diferencias entre las poblaciones. Todo muy bonito, pero no es un auténtico gráfico del análisis canónico de poblaciones, donde los círculos deben ser regiones de confianza con los datos sin escalar. Son círculos y no elipses los que representan a las poblaciones, ya que los ejes son independientes.
Por suerte, entre los resultados del objeto skulls.can1 disponemos de los coeficientes sin escalar que proporcionan las variables canónicas. Podemos pues aprovecharlos para calcular los scores de todos los datos sin escalar. Observemos también la utilización de la función aggregate para calcular las medias de cada población por separado como expliqué en el artículo Aplicar una función según el grupo.



# raw scores

scores <- as.matrix(skulls[,-1]) %*% skulls.can1$coeffs.raw
plot(scores[,1],scores[,2],xlim=c(-1,6),ylim=c(20,25), xlab="1er. eje canónico",ylab="2o. eje canónico",pch=16)

medias<-aggregate(skulls[,-1],skulls["EPOCH"],mean)

scores.medias <- medias %*% skulls.can1$coeffs.raw
text(scores.medias[,1],scores.medias[,2],1:5,pch=15,col="red")




Ahora ya sólo queda añadir los círculos de confianza para un nivel de confianza especificado, por ejemplo el 90%. Para ello vamos a calcular los radios de cada círculo según el tamaño de la muestra de cada población y la función symbols. Observemos que el objeto r es un vector.

resumen <- table(EPOCH)

g <- length(resumen) # número de poblaciones

p <- dim(skulls[,-1])[2] # número de variables

n <- as.vector(resumen) # tamaño de las muestras en cada población

radios <- function(g,p,n,conf.level=0.95) {

N <- sum(n)

F <- qf(conf.level,p,N-g-p+1)

sqrt(F*(N-g)*p/((N-g-p+1)*n))

}

r <- radios(g,p,n,0.90)


plot.new()

plot.window(xlim=c(0.5,4),ylim=c(21,24.5))

axis(1)

axis(2)

box()

title(main="Análisis canónico", xlab="1er. eje", ylab="2o. eje")

text(scores.medias[,1],scores.medias[,2],labels=levels(EPOCH),col="red")

symbols(scores.medias[,1],scores.medias[,2],circles=r,inches=FALSE,add=T,lwd=2,fg="red")


Bibliografía

Cuadras, C.M. (2010). Nuevos métodos de Análisis Multivariante. CMC Editions, Barcelona.

Everitt, B. (2005). An R and S-Plus Companion to Multivariate Analysis, Springer, London.

Arthur Thomson & R. Randall-Maciver (1905). The Ancian Races of the Thebaid. Oxford University Press.