La velocidad de la luz (1)

Simon Newcomb

Los viajes de la luz son rápidos, pero no instantáneos. La luz tarda alrededor de un segundo en llegar desde la luna y sobre 10 billones de años desde el objeto más distante que se ha observado en nuestro expansivo universo. Como la radio y el radar también viajan a la velocidad de la luz, un valor ajustado de esta velocidad es muy importante para la comunicación con los astronautas y los satélites en órbita. Un valor ajustado de la velocidad de la luz también es muy importante para los diseñadores de computadoras, ya que las señales eléctricas viajan a esta velocidad.

La primera medida razonablemente ajustada de la velocidad de la luz se hizo hace más de 100 años gracias a los experimentos de A. A. Michelson y Simon Newcomb. La siguiente tabla contiene las 66 medidas hechas por Newcomb entre julio y septiembre de 1882.

28 22 36 26 28 28
26 24 32 30 27 24
33 21 36 32 31 25
24 25 28 36 27 32
34 30 25 26 26 25
-44 23 21 30 33 29
27 29 28 22 26 27
16 31 29 36 32 28
40 19 37 23 32 29
-2 24 25 27 24 16
29 20 28 27 39 23
Tabla. Medidas de Newcomb del lapso de la luz

Un conjunto de datos como estos no tiene sentido sin la información sobre su contexto. Debemos contestar algunas preguntas iniciales sobre cualquier conjunto de datos. Primero, "¿Qué variable se está midiendo?" Newcomb midió cuanto tiempo tardaba la luz en ir y volver desde su laboratorio en el río Potomac hasta un espejo en la base del monumento a Washington, una distancia de aproximadamente 7400 metros. Tal como nosotros calculamos la velocidad de un coche como el tiempo necesario para recorrer un kilómetro, Newcomb calculó la velocidad de la luz a partir del tiempo de ese viaje.

Contestar la pregunta ¿Qué variable se está midiendo? requiere una descripción del instrumento utilizado para hacer la medida. Entonces podremos juzgar si la variable medida es apropiada para nuestro propósito. Este juicio frecuentemente pide un conocimiento de experto en el campo particular de estudio. Por ejemplo, Newcomb inventó un nuevo y complicado aparato para medir el lapso de la luz. Nosotros como estadísticos aceptamos el juicio de los físicos en el sentido de que este instrumento es apropiado para la tarea encomendada y más preciso que los instrumentos anteriores.

El estudio de Newcomb de la velocidad de la luz mide una variable claramente definida y fácilmente comprensible. Preguntas sobre las medidas son mucho más difíciles de contestar en los campos de las Ciencias Sociales y Económicas que en las Físicas. Nosotros nos podemos poner de acuerdo fácilmente en el tipo de medida apropiado para calcular la altura de una persona, pero ¿como medimos la inteligencia?

Los usuarios de datos deberían ser conscientes de que considerar los números con su valor nominal, sin pensar en la variable medida y el proceso utilizado para medirla, puede producir serias malinterpretaciones.

Las dos preguntas que faltan sobre cualquier conjunto de datos deberían contestarse con mayor senzillez: ¿Cuales son las unidades de medida? y ¿Como están registrados los datos?

El paquete MASS contiene la base de datos newcomb con los datos de la tabla anterior, pero desordenados. Es mejor obtener los datos originales en la siguiente dirección:

http://people.reed.edu/~jones/141/Newcomb.html

La primera medida de Newcomb del lapso de la luz fué 0,000024828 segundos. De manera que su unidad de medida fueron los segundos. Pero los valores de la tabla inicial no se parecen a 0,000024828. Estos números son incómodos de escribir y difíciles de tratar aritméticamente. En consecuencia nosotros hemos movido la coma decimal nueve posiciones a la derecha, esto es 24828, y hemos registrado únicamente la desviación respecto a 24800. Así pues, 28 es el resumen de 0,000024828 y -2 significa 0,000024798. Este procedimiento se conoce como codificación de los datos. Se debe codificar cuando los datos originales contienen muchas cifras de las cuales sólo algunas varían de observación en observación. Los datos codificados son más fáciles de leer. Además, si utilizamos una calculadora o una computadora, reducir el número de cifras mejora los cálculos aritméticos y su precisión.

Ahora que hemos entendido lo que significan los datos de la tabla y su procedencia, podemos empezar a mirarlos más de cerca y estudiarlos más a fondo. En cualquier caso, y aunque parezca que los datos ya se entienden, conviene siempre contestar las tres preguntas preliminares.

(continuará...)

Bibliografía

D.S. Moore & G.P. McCabe, Introduction to the Practice of Statistics, W.H. Freeman & Company.

Instalar R y Bioconductor en Ubuntu

Estos días he tenido que instalar Ubuntu en algunos ordenadores portátiles del departamento. Como es habitual los portátiles vienen con Windows 7 instalado, pero a la mayoría de compañeros y compañeras nos gusta trabajar con Linux y su distribución más popular es Ubuntu o Kubuntu (versión de Ubuntu con el escritorio KDE). La verdad es que ya quedan muy lejos los días en los que era muy difícil instalar una distribución de Linux. Ahora con Ubuntu es realmente muy sencillo.
Después de instalar el sistema operativo y sus aplicaciones por defecto, llega la hora de instalar las aplicaciones para trabajar en serio y, entre ellas, nuestro amado R. La cosa es tan sencilla como ir al Centro de software de Ubuntu del menú Aplicaciones y buscar en el apartado de Ciencia e ingeniería -> Matemáticas el famoso programa R Commander o, como alternativa RKWard (mejor en Kubuntu). Se instalan con un sólo click. Fabuloso.
Sin embargo, si deseamos tener siempre la última versión, es mejor añadir el repositorio oficial del R-project a nuestro conjunto de fuentes de programas. Para ello seguiremos las instrucciones del siguiente enlace

http://cran.es.r-project.org/bin/linux/ubuntu/

Si no hay ninguna dificultad, así tendremos la última versión, tanto para 32 como para 64 bits, y cuando salga la siguiente el propio sistema de Ubuntu nos avisará para actualizarnos.
El único detalle que puede sorprender al principio es la utilización de un APT seguro con las claves de Vincent Goulet. Basta con ejecutar las dos instrucciones que se muestran en una Terminal y ya está.

Una vez instalados todos los paquetes de R, ya que siempre se necesita alguno de ellos, queremos instalar Bioconductor.
Bioconductor es un conjunto de paquetes de R, más de 400, para el análisis de datos de genómica que ha tenido un brillante desarrollo en los últimos años. De hecho, el proyecto de Bioconductor es uno de los más activos y ha impulsado notablemente el desarrollo del propio R.
Instalar el Bioconductor básico es muy sencillo. En una sesión de R (mejor si la hemos iniciado como administrador, con "sudo R"), debemos introducir las instrucciones

> source("http://bioconductor.org/biocLite.R")
> biocLite()

Pero, tras esperar un buen rato que se bajen, descompriman y compilen un motón de paquetes, nuestra frustración será grande cuando veamos que algunos paquetes no se han instalado con un mensaje final del tipo

installation of package 'XML' had non-zero exit status

El problema es que la instalación por defecto de Ubuntu no permite compilar algunos paquetes de Bioconductor. La solución es instalar previamente todo el software necesario.
En el blog de James Reid está muy bien explicado. Los paquetes de Ubuntu necesarios se instalan así:

sudo apt-get install build-essential curl graphviz-dev libcurl4-gnutls-dev libboost-dev libgd2-xpm-dev libglu1-mesa-dev libgtk2.0-dev libmysqlclient-dev libxml2-dev mesa-common-dev sun-java6-jdk unixodbc-dev

De paso, dejo estas instrucciones como recordatorio para nuestras propias instalaciones.

El paquete HistData

Gracias a dos alumnos del postgrado de Bioestadística de la UOC he descubierto el paquete HistData que contiene, entre otros, los datos de Florence Nightingale y que utilicé en el artículo anterior.
El paquete HistData recoge algunos de los más famosos conjuntos de datos de la historia de la Estadística, como los datos de Sir Francis Galton que sirvieron para entrever la normal bivariante y los conceptos de correlación y regresión. Pero además, el paquete contiene la reproducción con R de gráficos famosos asociados a esos datos. De esta forma los docentes podemos hacer memoria histórica estadística, que siempre va bien.
Así, debo rectificar parcialmente mi afirmación en el sentido de que no hay una función de R que reproduzca el coxcomb o rosa de Nightingale. Si bien es estrictamente cierta, el código que acompaña los datos de Nightingale permite generar el gráfico que veis al principio de este artículo y que compara las frecuencias de soldados muertos por diversas causas, antes y después de aplicar las mejoras sanitarias.
Este código utiliza el lenguaje gráfico del paquete ggplot2 que es una auténtica maravilla. Aprender este lenguaje no es trivial, pero si se domina se pueden hacer gráficos que combinan los mejores aspectos de los gráficos base de R y los lattice. Ya tengo otro buen tema de estudio para este verano.

data(Nightingale)
# For some graphs, it is more convenient to reshape death rates to long format
# keep only Date and death rates
require(reshape)
Night<- Nightingale[,c(1,8:10)]
melted <- melt(Night, "Date")
names(melted) <- c("Date", "Cause", "Deaths")
melted$Cause <- sub("\\.rate", "", melted$Cause)
melted$Regime <- ordered( rep(c(rep('Before', 12), rep('After', 12)), 3), levels=c('Before','After'))
Night <- melted

require(ggplot2)
cxc <- ggplot(Night, aes(x = factor(Date), y=Deaths, fill = Cause)) +
# do it as a stacked bar chart first
geom_bar(width = 1, position="identity", color="black") +
# set scale so area ~ Deaths
scale_y_sqrt() +
facet_grid(. ~ Regime, scales="free", labeller=label_both)
# A coxcomb plot = bar chart + polar coordinates
cxc + coord_polar(start=3*pi/2) +
opts(title="Causes of Mortality in the Army in the East") +
xlab("")

La idea fundamental es que un coxcomb es un diagrama de barras con áreas proporcionales a los datos y en coordenadas polares. Justamente su definición. ¡Fantástico!

Los datos de Florence Nightingale

Florence Nightingale (Florencia, 12 de mayo de 1820 - Londres, 13 de septiembre de 1910), fue una enfermera británica considerada una de les pioneras en la práctica de la enfermería moderna y creadora del primer modelo conceptual de enfermería. Pero además destacó desde muy joven en matemáticas, aplicando después la estadística a la epidemiología y explotando la estadística sanitaria. Fue la primera mujer admitida en la Royal Statistical Society británica, y miembro honorario de la American Statistical Association.

Sin embargo, esta mujer, estudiada y reconocida en enfermería, es poco conocida entre los estudiantes de estadística y de matemáticas. Con motivo del pasado Día Internacional de la Mujer (8 de marzo) compartí con Carmina Olivé un acto de CCOO en la Universitat Politècnica de Catalunya (conjunto con la UB) en la que ella explicó su importancia en la enfermería moderna y yo glosé sus dotes científico-matemáticas-estadísticas. Voy a explicar brevemente mis argumentos.

A través de su trabajo como enfermera en guerra de Crimea, Florence Nightingale fue pionera al establecer la importancia de la higiene en los hospitales, en particular los hospitales de campaña. Ella reunió muchos datos en relación al coste en vidas por la falta de limpieza y, con sus gráficos y razonamientos matemáticos, también fue una avanzada en estadística aplicada. Con todo ello consiguió convencer a las autoridades de la época y, en especial, a los entonces muy obtusos militares para mejorar las condiciones sanitarias de los hospitales militares y civiles.

Entre los resultados estadísticos de Nightingale el más famoso es este gráfico, conocido como coxcomb o rosa de Nightingale, en el que se comparan los datos de soldados muertos por diversas causas, antes y después de aplicar las mejoras sanitarias en los hospitales.

Una versión más moderna y dinámica de este gráfico se puede hallar en el siguiente enlace:

http://understandinguncertainty.org/node/213

y las matemáticas y los datos para generar dicho gráfico se pueden consultar aquí:

http://understandinguncertainty.org/node/214

Por desgracia y que yo sepa, R no dispone de este tipo de gráfico. Llegados a este punto y, aunque hay algunas técnicas univariantes descriptivas para ayudar a Florence Nightingale, vamos a utilizar R y un par de técnicas multivariantes para describir sus datos.

La enfermera Nightingale recogió el número de soldados muertos por diversas causas: Zymotic diseases (ZD), Wounds & injuries (WI) y All other causes (OC), con ellos calculó unos índices de mortalidad anual relativa por cada 1000. Son precisamente esos índices los que vamos a utilizar como variables para un análisis de componentes principales.

dades <- read.table(file="dades.csv",header=T,sep=";")
dades$treatment <- factor(c(rep(0,12),rep(1,12)))
levels(dades$treatment) <- c("No", "Yes")
dades$active <- dades$size - (dades$zymotic_diseases + dades$wounds.injuries + dades$all_other)
dades$zymotic_diseases_ar <- dades$zymotic_diseases*12000/dades$size
dades$wounds.injuries_ar <- dades$wounds.injuries*12000/dades$size
dades$all_other_ar <- dades$all_other*12000/dades$size
attach(dades)

library(FactoMineR)
dades.res <- dades[8:11]
rownames(dades.res) <- idx
res.pca <- PCA(dades.res, quali.sup = 1)
plot(res.pca,habillage=1)

Con este código se obtiene el siguiente gráfico que se explica con mucha facilidad:

En él observamos que, si bien el primer eje admite la típica explicación de tamaño, en este caso depende básicamente de la causa evitable ZD y de OC (menos importante). El otro eje se explica con la variable WI.

En el gráfico con los 24 meses observados sobre las dos primeras componentes principales, donde los 12 primeros son antes de aplicar sus nuevos métodos de cuidado en los hospitales militares, observamos que los meses con el tratamiento se sitúan claramente a la izquierda, indicando la mejora en las muertes por causas evitables.

Por otra parte, y como los datos iniciales son frecuencias, podemos realizar un análisis de correspondencias tomando los meses como perfiles.

tabla <- data.frame(zymotic_diseases, wounds.injuries, all_other, active)
library(ca)
my.ca <- ca(tabla)
plot(my.ca, map="rowprincipal", xlim=c(-0.2,0.8), ylim=c(-0.2,0.2))

El resultado es similar. La explicación de los ejes que proporcionan los vértices (causas de mortalidad) es la misma que con el PCA e indica que la mejora en el tratamiento de las enfermedades en los hospitales hace que dichos meses estén a la izquierda. El único detalle es que en el gráfico se ha tenido que hacer un zoom sobre los puntos, ya que todos ellos quedaban fuertemente agrupados junto al vértice de "activos" que es el más pesado (el triángulo rojo del gráfico).

Discriminador lineal de Fisher

Supongamos que queremos discriminar los datos entre dos poblaciones normales con diferente vector de medias y la misma matriz de covarianzas. Vamos a ver como hacerlo en el caso bivariante ya que así podemos representar el problema gráficamente.
En vez de realizar el problema con datos simulados o reales, lo vamos a tratar con dos distribuciones teóricas concretas.
Así pues, consideremos la distribución de una población normal bivariante con media mu1=(2.5,4) y matriz de covarianzas Sigma=(2,1,1,2) y otra población con la misma matriz de covarianzas, pero de media mu2=(6,3). El gráfico de arriba representa las densidades bivariantes de estas dos poblaciones. Un segmento de puntos une los dos puntos medios de las poblaciones.

La idea de Fisher consistió en hallar una combinación lineal de las variables originales (X1,X2) de la forma w1X1+w2X2 y tal que discrimine "el máximo posible" las dos poblaciones. Él mismo definió el criterio de máxima discriminación como maximizar la razón entre la suma de cuadrados entre grupos y la suma de cuadrados dentro de los grupos sobre la combinación lineal (discriminant scores). La solución viene dada por el vector w (o es proporcional al vector w)

w = Sigma^{-1} * (mu2 - mu1)

Es decir, entre todas las combinaciones lineales de la forma w1X1+w2X2, la que mejor discrimina es justamente esa. En el ejemplo propuesto, el vector solución w se ha dibujado en el gráfico como una flecha desde el punto mu1. Observemos que w tiene en cuenta la diferencia entre las medias pero también la relación de dependencia entre las variables.

Para visualizar este resultado he dibujado las densidades de las dos poblaciones si consideramos el vector w=(1,0), es decir, que la combinación lineal considerada sea simplemente X1.

Esta no es la solución óptima.

La solución exacta es w=(8/3,-11/6)

En el gráfico inicial, la recta que discrimina las dos poblaciones tiene un vector director ortogonal al vector w. Esa recta es la recta de puntos que equidistan de los dos puntos medios según la distancia de Mahalanobis.

El código para hacer todos los cálculos y los gráficos es el siguiente:


mu1 <- c(2.5,4)
mu2 <- c(6,3)
Sigma <- matrix(c(2,1,1,2),ncol=2)

dnormbv <- function(x1,x2,mu,Sigma) {
sigma1 <- sqrt(Sigma[1,1])
sigma2 <- sqrt(Sigma[2,2])
rho <- Sigma[1,2]/(sigma1*sigma2)
cte <- 1/(2*pi*sigma1*sigma2*sqrt(1-rho^2))
cte * exp((-1/2) * (1/(1-rho^2)) * (
(x1-mu[1])^2/sigma1^2 - 2*rho*(x1-mu[1])*(x2-mu[2])/(sigma1*sigma2) + (x2-mu[2])^2/sigma2^2)
)}

x1<-seq(0,6,0.1)
x2<-seq(1,7,0.1)
z <- matrix(rep(0,length(x1)*length(x2)),ncol=length(x2))
for (i in 1:length(x1)) {
for (j in 1:length(x2)) {
z[i,j] <- dnormbv(x1[i],x2[j],mu1,Sigma)
}
}
plot(outer(x1,x2),xlim=c(0,9),ylim=c(0,9),xlab="",ylab="",type="n")
contour(x1,x2,z,nlev=4,add=T)

x1<-seq(2,9,0.1)
x2<-seq(0,7,0.1)
z <- matrix(rep(0,length(x1)*length(x2)),ncol=length(x2))
for (i in 1:length(x1)) {
for (j in 1:length(x2)) {
z[i,j] <- dnormbv(x1[i],x2[j],mu2,Sigma)
}
}
contour(x1,x2,z,nlev=4,add=T)

segments(mu1[1],mu1[2],mu2[1],mu2[2],lty="dotted")

SigmaInv <- matrix(c(2/3,-1/3,-1/3,2/3),ncol=2)

w <- SigmaInv %*% (mu2-mu1)

arrows(mu1[1],mu1[2],mu1[1]+w[1],mu1[2]+w[2],lwd=2)

pm <- (mu1+mu2)/2
a <- w[1] * pm[1] / w[2] + pm[2]
b <- -w[1]/w[2]
abline(a,b,lty="dashed")


ee <- 3.5
x1 <- seq(mu1[1]-ee,mu1[1]+ee,by=0.1)
x2 <- dnorm(x1,mean=mu1[1],sd=sqrt(Sigma[1,1]))
plot.new()
plot.window(xlim=c(-1,9),ylim=c(0,0.5))
axis(1)
axis(2)
lines(x1,x2)
x1 <- seq(mu2[1]-ee,mu2[1]+ee,by=0.1)
x2 <- dnorm(x1,mean=mu2[1],sd=sqrt(Sigma[1,1]))
lines(x1,x2)



E1 <- t(w) %*% mu1
E2 <- t(w) %*% mu2
var.w <- as.numeric(t(w) %*% Sigma %*% w)

ee <- 7
x1 <- seq(E1-ee,E1+ee,by=0.1)
x2 <- dnorm(x1,mean=E1,sd=sqrt(var.w))
plot.new()
plot.window(xlim=c(-6,16),ylim=c(0,0.5))
axis(1)
axis(2)
lines(x1,x2)
x1 <- seq(E2-ee,E2+ee,by=0.1)
x2 <- dnorm(x1,mean=E2,sd=sqrt(var.w))
lines(x1,x2)



Con datos reales o simulados podemos utilizar la función lda() del paquete MASS.

Bibliografía

Lattin, J. et al, Analyzing Multivariate Data, Ed. Brooks/Cole, Belmont (2003).

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.

Aplicar una función según el grupo

Para aplicar una función a unas columnas de datos, separadamente según las clases o grupos que queramos, en R disponemos de dos funciones: tapply() y aggregate().
La función tapply() se aplica sobre un único vector de datos, cuyos valores se agrupan en función de una o unas variables o factores. Con los datos mtcars, por ejemplo, para saber el máximo valor de mpg según el número de cilindros hacemos

> data(mtcars)
> attach(mtcars)
> tapply(mpg, cyl, max)
4 6 8
33.9 21.4 19.2

En este caso el resultado es un vector, aunque también puede ser una lista si el resultado de la función aplicada no es un escalar:

> tapply(mpg,cyl,range)
$'4'
[1] 21.4 33.9

$'6'
[1] 17.8 21.4

$'8'
[1] 10.4 19.2

Acceder a estos resultados significa utilizar las llamadas a elementos de una lista:

> rangos <- tapply(mpg, cyl, range)
> rangos[[1]]
[1] 21.4 33.9

> rangos$'4'
[1] 21.4 33.9

> rangos[['4']]
[1] 21.4 33.9

> rangos[["4"]]
[1] 21.4 33.9

Cuando se utiliza más de una variable de agrupación y el resultado de la función a aplicar no es un escalar, el valor que retorna tapply() es más difícil de gestionar.

> rangos2 <- tapply(mpg, mtcars[c("cyl","am")], range)
> rangos2
am
cyl 0 1
4 Numeric,2 Numeric,2
6 Numeric,2 Numeric,2
8 Numeric,2 Numeric,2

> rangos2["4","0"]
[[1]]
[1] 21.5 24.4

Al llamar la función tapply() sin el argumento función resulta un índice de agrupación, según los grupos determinados, que puede ser útil en algunos casos:

> idx <- tapply(mpg,mtcars[c("cyl","am")]) > idx
[1] 5 5 4 2 3 2 3 1 1 2 2 3 3 3 3 3 3 4 4 4 1 3 3 3 3 4 4 4 6 5 6 4

Por otra parte, si lo que se pretende es resumir una o más columnas de un data.frame o matriz mediante un estadístico o función escalar, entonces utilizaremos la función aggregate(), donde el segundo argumento debe ser una lista.

> aggregate(mtcars[c("mpg","hp","wt")], mtcars$cyl, mean)
Error in aggregate.data.frame(mtcars[c("mpg", "hp", "wt")], mtcars$cyl, :
'by' must be a list
Calls: aggregate -> aggregate.data.frame

> aggregate(mtcars[c("mpg","hp","wt")], mtcars["cyl"], mean)
cyl mpg hp wt
1 4 26.66364 82.63636 2.285727
2 6 19.74286 122.28571 3.117143
3 8 15.10000 209.21429 3.999214

> aggregate(mtcars[c("mpg","hp","wt")], mtcars[c("cyl","am")], mean)
cyl am mpg hp wt
1 4 0 22.90000 84.66667 2.935000
2 6 0 19.12500 115.25000 3.388750
3 8 0 15.05000 194.16667 4.104083
4 4 1 28.07500 81.87500 2.042250
5 6 1 20.56667 131.66667 2.755000
6 8 1 15.40000 299.50000 3.370000

Finalmente, cuando el problema es mas complejo y se trata de aplicar una función no escalar a más de un vector según una determinada agrupación, deberemos combinar dos procedimientos como split() para agrupar y sapply() o lapply() para aplicar.

> xx <- split(mtcars[c("mpg","hp","wt")], mtcars["cyl"])

> aggregate(xx$'4', list(cyl=rep(4, dim(xx$'4')[1])), min)

cyl mpg hp wt

1 4 21.4 52 1.513


> sapply(xx, cor)

4 6 8

[1,] 1.0000000 1.0000000 1.00000000

[2,] -0.5235034 -0.1270678 -0.28363567

[3,] -0.7131848 -0.6815498 -0.65035801

[4,] -0.5235034 -0.1270678 -0.28363567

[5,] 1.0000000 1.0000000 1.00000000

[6,] 0.1598761 -0.3062284 0.01761795

[7,] -0.7131848 -0.6815498 -0.65035801

[8,] 0.1598761 -0.3062284 0.01761795

[9,] 1.0000000 1.0000000 1.00000000

Otra opción es utilizar la función by() que generaliza la función tapply().

> by(mtcars[c("mpg","hp","wt")],mtcars["cyl"],cor)
cyl: 4
mpg hp wt
mpg 1.0000000 -0.5235034 -0.7131848
hp -0.5235034 1.0000000 0.1598761
wt -0.7131848 0.1598761 1.0000000
------------------------------------------------------------
cyl: 6
mpg hp wt
mpg 1.0000000 -0.1270678 -0.6815498
hp -0.1270678 1.0000000 -0.3062284
wt -0.6815498 -0.3062284 1.0000000
------------------------------------------------------------
cyl: 8
mpg hp wt
mpg 1.0000000 -0.28363567 -0.65035801
hp -0.2836357 1.00000000 0.01761795
wt -0.6503580 0.01761795 1.00000000