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

Aplicar una función a una matriz o array

La función apply

Para aplicar una función a una matriz o array, R dispone de la función apply que se utiliza con tres parámetros: el objeto array (una matriz es un array de dos dimensiones), la dimensión sobre la que actuaremos y la función que se aplicará. Como con la función sapply, es posible añadir al final todos los argumentos que precise la función que aplicaremos. En el caso de las matrices, si el segundo argumento es un 1, significa que la función se aplicará a las filas, mientras que si es un 2, se aplicará a las columnas.
En muchos casos se puede utilizar una instrucción apply en vez de un bucle. Además, el resultado es un vector o una matriz con nombres (si los había) y dicen que a menudo es más eficiente que un bucle.
Si queremos aplicar un procedimiento, primero definiremos la función a aplicar. Por ejemplo, si queremos saber el número de datos en cada una de las variables (columnas) que forman una base de datos (matriz), primero definiremos la función

n.datos <- function(x) n <- sum(!is.na(x))

y la aplicaremos sobre la base de datos

apply(base.de.datos, 2, n.datos)

La función apply se puede aplicar sobre elementos que no son arrays, pero debemos saber que entonces trata de convertirlos mediante un as.matrix o un as.array, dependiendo de las dimensiones.
Supongamos ahora que queremos estandarizar cada una de las variables de una base de datos, pero utilizando la mediana como estadístico de centralidad y la MAD como medida de dispersión.
El truco consiste en utilizar la función scale pero con los parámetros adecuados:

datos.estand <- scale(datos, center=apply(datos, 2, median), scale=apply(datos, 2, mad))

Sin embargo, si lo que queremos es calcular la suma o la media, es mejor utilizar funciones como rowSums, colSums, rowMeans o colMeans que son más eficientes. Además tienen un argumento na.rm= que permite prescindir de los valores faltantes.

La función sweep

Otra situación se presenta cuando queremos procesar una matriz por filas o por columnas, pero cada fila o columna de forma distinta, dependiendo de los valores de un vector dado. En estos casos utilizaremos la función sweep que tiene básicamente 4 argumentos. Como en la función apply, los dos primeros son la matriz y la dimensión sobre la que aplicaremos la función. El tercero es el vector de valores para procesar cada fila o columna de forma distinta y finalmente, el cuarto argumento es la función a aplicar.
Como funciones a aplicar podemos utilizar los operadores binarios como la suma "+", la resta "-" (por defecto), el producto "*" y el cociente "/", siempre con comillas. Por ejemplo, para dividir cada columna de unos datos por su máximo haremos

maximos <- apply(datos, 2, max)

sweep(datos, 2, maximos, "/")


Ahora bien, para otro tipo de funciones deberemos asegurarnos que la función trabaje correctamente con sweep. No todas lo hacen.
Si la función a aplicar calcula correctamente lo que queremos para cada vector de la matriz, pero no lo hace con sweep, podemos pensar en la función mapply.
La función mapply es la versión multivariante de sapply.
mapply aplica la función (primer argumento) sobre los primeros elementos de cada argumento, los segundos elementos, los terceros y así sucesivamente. El resultado se trata de simplificar, normalmente es una lista o vector.

mapply(rep, 1:4, 4:1)

Si la función necesita más argumentos, se los podemos pasar así:

mapply(rep, 1:4, MoreArgs=list(x=5))

Con un data.frame datos podemos hacer

mapply("/", datos, maximos)

y obtendremos un resultado similar al de sweep.

Aplicar una función a un vector o una lista

Ciertamente R está preparado para operar todos los elementos de un vector de forma casi "natural". La mayoría de funciones son vectoriales. No ocurre lo mismo con una lista. Aunque muchas funciones devuelven una lista, que permite mucha flexibilidad, no se puede aplicar una función cualquiera directamente a una lista. Para ello disponemos de las funciones lapply y sapply.
Las dos funciones admiten, básicamente, dos argumentos: el primero es el objeto vector o lista sobre el que se aplica la función expresada en el segundo argumento. La diferencia entre ambas es que mientras lapply devuelve una (l)ista, sapply siempre tratará de (s)implificar el resultado en un vector, si es posible.

> mi.lista <- list(a=1:10, b=letters[1:3], cc=c(TRUE,FALSE,TRUE,FALSE))

> length(mi.lista)
[1] 3

> lapply(mi.lista, length)
$a
[1] 10

$b
[1] 3

$cc
[1] 4

> sapply(mi.lista, length)
a b cc
10 3 4

Observamos que mientras la función length aplicada a la lista nos devuelve el número de elementos de la lista (3), aplicada mediante un lapply o sapply nos da la longitud de sus elementos.

Cuando el objeto del primer argumento no es un vector o una lista, se forzará automáticamente que sea una lista mediante la función as.list. De este modo podemos aplicar estas funciones a un data.frame que no es estrictamente una lista, pero su conversión es sencilla.

> mi.df <- data.frame(a=1:10, b=letters[1:10], ca=runif(10), cb=rnorm(10))
> class(mi.df)
[1] "data.frame"

> sapply(mi.df,class)
a b ca cb
"integer" "factor" "numeric" "numeric"

> mi.df.num <- mi.df[ ,sapply(mi.df,class)=="numeric"]

> sapply(mi.df.num, mean)
ca cb
0.59009349 0.02423536

Cuando la función a aplicar (segundo argumento) necesita fijar sus propios argumentos, éstos se pueden incluir en las funciones lapply o sapply.

> sapply(mi.df.num, mean, trim = 0.05)

Por último, como estas funciones sirven para repetir el mismo cálculo sobre los elementos de un vector o lista, podemos pensar en utilizarlas siempre que podamos en lugar de un bucle (loop).
Además, tenemos una función replicate que permite simulaciones como la siguiente:

hist(replicate(100,mean(rexp(10))))

Fórmulas matemáticas en un gráfico

Algunas de las funciones que permiten introducir texto en un gráfico son text, mtext, axis, title,...
Si en el argumento que explicita el texto de alguna de estas funciones se escribe una expresión, R la interpretará como una expresión matemática y le dará formato al estilo TeX. Dichas expresiones son el resultado de una función expression().
En una expresión de R ciertos nombres se interpretan como símbolos matemáticos, por ejemplo, alpha será la letra griega, sum será el símbolo de sumatorio, etc.
Para ver las diversas posibilidades, se puede consultar la ayuda help(plotmath) o la demostración demo(plotmath).

El gráfico inicial se consigue con el siguiente código:

curve(dnorm(x),-3,3,axes=F)
box()
axis(2)
axis(1,at=0,labels=c(expression(mu)))
title("Densidad normal")
text(2, 0.35, expression(paste(f(x)==frac(1, sigma*sqrt(2*pi)), " ",
                          plain(e)^{frac(-(x-mu)^2, 2*sigma^2)})),
   cex = 1.2)

En ciertas situaciones, como por ejemplo para definir una función, se necesita combinar texto con valores y variables. En ese caso no es posible utilizar la función expression() ya que las variables se tratarían literalmente como texto. Para solucionar este caso se utiliza la función substitute().
Ejemplo:

mifunc <- function(media) {           
          text(2, 3, substitute(paste("El valor de ", bar(x), " es ", media)))
          } 

Este tipo de fórmulas matemáticas se puede reproducir en cualquiera de los dispositivos gráficos de pantalla como X11, Windows y Quartz y, gracias a la información de las fuentes Adabe Tipo 1 estándar que tiene R, también en PostScript y PDF.

Clases S3 y S4

El sistema de clases de los objetos en R proporciona alguno de los mecanismos de la programación orientada a objetos como el despacho del método (method dispatch) y la herencia.
El método es la implementación de un algoritmo que representa una operación o función que un objeto realiza. El conjunto de los métodos de un objeto determinan el comportamiento del objeto. En R, el despacho del método consiste en examinar la clase de los argumentos de una función para decidir (despachar) la versión adecuada para los objetos de esa clase. No todas las funciones de R tienen despacho del método. Las que sí lo tienen se llaman funciones genéricas.
La herencia permite a los programadores crear nuevas clases, similares a otras ya existentes. Únicamente deberán proporcionar métodos adecuados para las nuevas clases o mantener los heredados. Un objeto de R que hereda las propiedades de un objeto ya definido, tiene como atributo de clase un vector que contiene la clase de ese objeto (en primer lugar), junto con las clases del objeto del que hereda.

El primer mecanismo de la programación orientada a objetos en R es el conjunto de clases S3 o del viejo estilo (old-style), donde el despacho del método se produce a través de las funciones genéricas del siguiente modo:

Supongamos que tenemos un objeto con el nombre de cdr de la clase Cuadrado, el cual es una subclase de Rectangulo que a su vez es una subclase de Forma. El mecanismo de despacho del método consiste en S3 y la función UseMethod. Utilizando S3 definimos un método Area para la clase Rectangulo como

  Area.Rectangulo <- function(objeto) { attr(objeto, "ladoA") * attr(objeto, "ladoB"); } 

dado que un objeto Rectangulo tiene dos atributos ladoA y ladoB (ahora se puede acceder directamente a los argumentos con el operador @, es decir, objeto@ladoA y objeto@ladoB, respectivamente). Entonces definimos la función genérica Area así:

  Area <- function(objeto, ...) UseMethod("Area"); 

Cuando esta función se aplica en el objeto con Area(cdr), UseMethod despachará el método basado en la clase del primer argumento, es decir, objeto. Como cdr es de la clase Cuadrado, primero buscará un método llamado Area.Cuadrado. Si no existe, como en este caso, probará con Area.Rectangulo y así sucesivamente. Si no existiera ningún método específico para cualquiera de las clases del objeto, entonces se recurrirá al método Area.default que siempre debemos tener.

Observemos que las funciones genéricas S3 se pueden reconocer por la utilización de UseMethod en su código. Esto es importante, ya que las páginas de ayuda para una combinación de método y objeto depende de su nombre completo del tipo "function.class". Por ejemplo, la página de ayuda de la función summary no explica absolutamente nada sobre su actuación cuando le pasamos un objeto de la clase factor. Será mejor buscar la ayuda de la función summary.factor, aunque para que la función actúe sólo hay que escribir summary(objeto).

Cuando se crea una clase, asociadas a ella, deberemos crear un conjunto de funciones para extraer datos o información sobre los objetos de esa clase. Dada la convención explicada sobre los nombres de las funciones en las clases S3, podemos utilizar la función apropos para hallar los métodos disponibles para una clase:

> apropos('.*\\.factor$')
 [1] "all.equal.factor"     "as.character.factor"  "as.data.frame.factor"
 [4] "as.Date.factor"       "as.factor"            "as.list.factor"  
 [7] "as.POSIXlt.factor"    "as.vector.factor"     "codes.factor"    
[10] "[<-.factor"           "[.factor"             "[[.factor"             
[13] "format.factor"        "is.factor"            "is.na<-.factor"        
[16] "length<-.factor"      "levels<-.factor"      "Math.factor"           
[19] "Ops.factor"           "print.factor"         "rep.factor"            
[22] "summary.factor"       "Summary.factor"       "xtfrm.factor"          

Así descubriremos que existen algunas funciones específicas para la clase factor.

Ahora bien, como el despachado del método que proporcionan las clases S3 está limitado al primer argumento de la función, y como las convenciones sobre los nombres que hemos explicado pueden crear alguna confusión, se decidió crear un nuevo sistema de clases S4 o de nuevo estilo (new-style). Éste es ahora el sistema preferido en el desarrollo de R.
En las clases S4, las funciones genéricas se reconocen por la llamada a una función standardGeneric en su definición.
Las funciones necesarias para trabajar con las clases S4 se hallan en el paquete methods. Por ejemplo, para saber si un objeto utiliza el nuevo estilo podemos hacer:

> temp <- c(20,15,15,20,20,30,25,25,30,15,25,30,30)  
> temp <- factor(temp)  
> isS4(temp)
[1] FALSE

Para saber los métodos asociados a una clase S4 como los objetos mle, resultado de la función de estimación de máxima verosimilitud, hacemos:

> library(methods)
> showMethods(class='mle')

Descubriremos que se puede calcular la matriz de varianzas-covarianzas de un objeto mle con la función vcov.

Aunque no hay una función genérica print para las clases S4, la función show permite ver el contenido de un objeto.
Por otra parte, los elementos que componen un objeto S4 se guardan en los llamados slots. Para ver los tipos de slots en un objeto, podemos utilizar la función showClass.

> library(stats4)
> showClass("mle")
Class “mle” [package "stats4"]

Slots:
                                                                        
Name:       call      coef  fullcoef      vcov       min   details minuslogl
Class:  language   numeric   numeric    matrix   numeric      list  function
            
Name:     method
Class: character

Para acceder directamente a un slot de un objeto, utilizaremos el operador @ del mismo modo que el operador $ para una lista. La función slot también obtiene el mismo resultado.

> slot(objeto.lme,"vcov")

Primeras jornadas de R

Cuando me disponía a redactar un resumen de las Primeras jornadas de R celebradas en la Universidad de Murcia los pasados días 26 y 27 de noviembre, he descubierto con agrado que se me han adelantado.
En el blog Análisis y decisión podéis leer un amplio resumen de Carlos Gil Bellosta.
En fin, un ejemplo de trabajo colaborativo, él lo escribe y yo lo cito. ¡Ejem! Ya sé que no es eso.
Simplemente debo añadir que corroboro todas las buenas impresiones que nos llevamos las personas que asistimos del Departamento de Estadística de la Universidad de Barcelona: Miquel Calvo, Esteban Vegas y yo mismo. Además de conocernos y hablar de R y del software libre en general, también han quedado algunos compromisos de colaboraciones que poco a poco seguro que iremos concretando.
Un abrazo a todos y todas los que estuvisteis allí.

Factores numéricos

Población africana con SIDA

Frecuentemente es necesario convertir una variable numérica en factor ya que algunas funciones de R así lo esperan. Sin embargo, si disponemos únicamente de un factor no podremos calcular algunos estadísticos u otras operaciones numéricas, aunque los niveles (o sus etiquetas) sean aparentemente numéricos.

> temp <- c(20,15,15,20,20,30,25,25,30,15,25,30,30)
> temp <- factor(temp)
> temp
 [1] 20 15 15 20 20 30 25 25 30 15 25 30 30
Levels: 15 20 25 30
> mean(temp)
[1] NA
Warning message:
In mean.default(temp) : argument is not numeric or logical: returning NA

Si no disponemos de los datos numéricos originales y deseamos convertir el factor a datos numéricos podemos probar así:

> temp.n <- as.numeric(temp)
> temp.n
 [1] 2 1 1 2 2 4 3 3 4 1 3 4 4

Pero el resultado son los valores enteros en los que se codifica internamente el factor.
Mejor si primero convertimos el factor a caracteres con las etiquetas de los niveles y luego esos mismos a valores numéricos:

> temp.n <- as.numeric(as.character(temp))
> temp.n
 [1] 20 15 15 20 20 30 25 25 30 15 25 30 30

Por otra parte, para crear un factor a partir de una variable continua se utiliza la función cut.
Los siguientes datos corresponden a la tasa de mortalidad del SIDA por cada mil habitantes en los países africanos en el año 2007:

> tasa <-  c(7.21,1.15,10.49,2.86,2.37,3.79,2.49,4.88,0.81,4.70,2.10,1.97,0.49,0.65,
            0.89,8.96,1.30,5.84,2.53,1.29,0.65,8.76,0.62,1.38,0.80,1.70,0.47,2.46) 

> tasa.f <- cut(tasa,breaks=seq(0,12,2))   

> table(tasa.f)
tasa.f
(0,2]   (2,4]   (4,6]   (6,8]  (8,10] (10,12]
   14       7       3       1       2       1

> class(tasa.f)
[1] "factor"