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"

Valores perdidos

En R, para representar un valor perdido (missing value) o que falta se utiliza el valor NA (not available) sin comillas. Se puede asignar este valor a una variable de la forma

> x <- c(3,2,NA,6)

pero no se puede consultar así

> x[3] == NA
[1] NA

para ello tenemos la función is.na():

> is.na(x[3])
[1] TRUE

Por otra parte, al realizar un determinado cálculo, como dividir por cero o la raíz de un negativo, el resultado puede ser Inf o NaN (not a number). La función is.nan() permite dirimir si se trata de uno de estos valores.

Cuando se importan datos de algún tipo de archivo, los valores perdidos pueden ser un problema si no se tratan de forma adecuada. Una posibilidad es utilizar el argumento na.strings= de la función read.table, al que le podemos pasar un vector con todos los valores del tipo carácter que R debe considerar perdidos y dar valor NA. En todo caso vale la pena revisar el resultado de la importación. Si los datos provienen de una hoja de cálculo, siempre podemos substituir los valores perdidos por NA, antes de la exportación del archivo al formato CSV.

Trabajar con valores perdidos

Algunas funciones de R disponen de argumentos para trabajar con datos que contienen valores perdidos. Por ejemplo, la mayoría de funciones estadísticas del tipo mean, var, sum, min, max, etc.
tienen un argumento na.rm= que se puede pasar a TRUE para que suprima los valores perdidos del cálculo. Para otras funciones que no tengan este parámetro podemos eliminar dichos valores creando un nuevo vector:


> x[!is.na(x)]
[1] 3 2 6

También algunas funciones de modelización estadística, como lm, glm, etc., disponen de un argumento na.action= al que le podemos pasar la acción a realizar para los valores perdidos. La más común es na.action=na.omit con lo que prescindirá de dichos valores en su procedimiento.
Otra opción es seleccionar del vector, matriz o data.frame(s) los individuos o filas que no tienen valores perdidos mediante la función complete.cases que devuelve un vector lógico:

> complete.cases(x)
[1] TRUE TRUE FALSE TRUE

Finalmente, señalar que en la conversión de un vector a factor, los valores perdidos no se consideran.

> y <- factor(x)
> y
[1] 3 2 6
Levels: 2 3 6

Si deseamos formar un nivel con los valores perdidos haremos

> y <- factor(x, exclude=NULL)
> y
[1] 3 2 NA 6
Levels: 2 3 6 NA

De R a ODF: el paquete odfWeave

En el artículo LaTeX y R ya comenté las ventajas de combinar texto y fórmulas junto con instrucciones de R para obtener un informe dinámico que contendrá el texto, las fórmulas y los resultados numéricos y gráficos generados por R. Lamentablemente esto supone saber LaTeX (y R), y utilizar LaTeX no es sencillo. Sin embargo, la potencia de la función Sweave en R se ha extendido a otros formatos más amigables. Primero fué el paquete R2HTML para el formato html y ahora el paquete odfWeave para el formato ODF u OpenDocument.
El paquete odfWeave proporciona las funciones y el entorno para escribir informes automáticos en un formato libre y gratuito (como R), que se puede modificar a posteriori de forma muy sencilla y que permite su exportación a otros formatos como html o PDF.
La idea consiste en escribir un archivo ODF con un programa editor como OpenOffice con comentarios o explicaciones y código R para generar resultados y gráficos. La función odfWeave procesa ese archivo y nos entrega el documento ODF con el mismo texto y los resultados numéricos o gráficos insertados. Ese documento se puede modificar añadiendo otros elementos: vínculos, fotografías, encabezados y piés de página, etc. y finalmente imprimir o exportar.
Por ahora sólo funciona con documentos de texto, pero más adelante se puede ampliar a presentaciones y hojas de cálculo.

Veamos un ejemplo.

El documento de nombre ejemplo.odt contiene el siguiente texto:

---

Ejemplo

En este documento vamos a probar el paquete odfWeave.

Primero cargamos los datos de forma oculta para el documento final.

<<>>=

# Los datos son

coches <- mtcars[ ,1:6]

@

Los datos a estudiar contienen \Sexpr{dim(coches)[2]} variables medidas sobre \Sexpr{dim(coches)[1]} coches.

Podemos insertar una tabla resumen de los datos:

<<>>=

medias <- apply(coches,2,mean)

odfTable(medias, horizontal = TRUE)

@

También podemos añadir algún bonito gráfico:

<<>>=

pairs(coches)

@

Y así seguiría nuestro trabajo.

---

Si aplicamos el código:

library(odfWeave)
odfWeave(ejemplo.odt,resultado.odt)

obtenemos el documento resultado.odt

---

Ejemplo

En este documento vamos a probar el paquete odfWeave.

Primero cargamos los datos de forma oculta para el documento final.

Los datos a estudiar contienen 6 variables medidas sobre 32 coches.

Podemos insertar una tabla resumen de los datos:

mpg	cyl	disp	hp	drat	wt
20.091	6.188	230.722	146.688	3.597	3.217

También podemos añadir algún bonito gráfico:

Y así seguiría nuestro trabajo.

---

El formato de las imágenes y su tamaño se controlan con las funciones getImageDefs y setImageDefs.
Del mismo modo, el formato del texto o de las tablas se controla con las funciones de definición de estilo getStyleDefs y setStyleDefs, y las de asignación getStyle y setStyle para determinar algún elemento.

La clase odfTable se utiliza para convertir vectores, matrices y data frames a tablas ODF. La función odfCat sirve para escribir texto directamente en el formato ODF. También hay funciones para crear listas y para insertar imágenes externas.

Una vez tenemos el documento resultado, podemos utilizar OpenOffice para convertir manualmente este documento a otros formatos. Para ello se usa el menú Guardar como... o Exportar... del menú Archivo. En algunos sistemas operativos también disponemos de programas conversores desde la consola.

Para utilizar el paquete odfWeave es necesario que nuestro ordenador disponga de un compresor/descompresor de archivos ZIP. También deberemos tener en cuenta que los gráficos se graban en formato PNG, por si nuestro sistema no tiene el dispositivo adecuado.

Test de Levene para la igualdad de varianzas

El test de Levene (1960) se usa para contrastar si k muestras tienen la misma varianza, es decir, la homogeneidad de varianzas. Otros contrastes, como por ejemplo el análisis de la varianza, suponen que las varianzas son iguales para todos los grupos. De ahí la importancia de verificar con el test de Levene esa hipótesis.

Este test es una alternativa al test de Bartlett. El test de Levene es menos sensible a la falta de normalidad que el de Bartlett. Sin embargo, si estamos seguros de que los datos provienen de una distribución normal, entonces el test de Bartlett es el mejor.

El test de Levene se resuelve con un ANOVA de los valores absolutos de las desviaciones de los valores muestrales respecto a un estadístico de centralidad (media, mediana o media truncada) para cada grupo. La elección del estadístico de centralidad de los grupos detemina la robustez y la potencia del test. Por robustez se entiende la habilidad del test para no detectar erróneamente varianzas distintas, cuando la distribución no es normal y las varianzas son realmente iguales. La potencia significa la habilidad del test para señalar varianzas distintas, cuando efectivamente lo son.

El artículo original de Levene proponía la media como estadístico de centralidad. Brown y Forsythe (1974) extendieron este test al utilizar la mediana e incluso la media truncada al 10%. Sus estudios de Monte Carlo mostraron que la utilización de la media truncada mejoraba el test cuando los datos seguían una distribución de Cauchy (colas grandes) y la mediana conseguía mejorarlo cuando los datos seguían una (distribución asimétrica). Con la media se consigue el mejor test para distribuciones simétricas y con colas moderadas.

Así pues, aunque la elección óptima depende de la distribución de los datos, la definición del test basada en la mediana es la recomendación general ya que proporciona una buena robustez para la mayoría de distribuciones no normales y, al mismo tiempo, una aceptable potencia. Si conocemos la distribución de los datos, podemos optar por alguna otra de las opciones.

El paquete car Companion to Applied Regression de J. Fox trae una función levene.test con la mediana como medida de centralidad de los grupos. El mismo test está incorporado al paquete Rcmdr.

Aplicar el test de Levene a unos datos con la opción de la media truncada como medida de centralidad es un ejercicio sencillo:

datos <- c(rcauchy(50,0,10),rcauchy(50,0,20),rcauchy(50,0,30))

grupo <- gl(3,50)

boxplot(datos~grupo)

func <- function(x) mean(x,trim=0.1)
medias_trunc <- tapply(datos,grupo,func)
desv <- abs(datos-medias_trunc[grupo])
summary(aov(desv~grupo))

Referencias

Brown, M. B. and Forsythe, A. B. (1974), Journal of the American Statistical Association, 69, 364-367.

Fox, J. (2002), An R and S-PLUS Companion to Applied Regression, Sage.

Levene, H. (1960). In Contributions to Probability and Statistics: Essays in Honor of Harold Hotelling, I. Olkin et al. eds., Stanford University Press, pp. 278-292.