R

Area:	Estadística
Web del proyecto:	r-project

R es un software de análisis estadístico de código abierto basado en el lenguaje S, desarrollado a principio de los años 90. El programa S-Plus utiliza el mismo lenguaje S, pero se trata de un software propietario y con un precio prohibitivo para la mayoría de los mortales.

Instalación y ejecución

Sistemas M$

Los binarios para la instalación en sistemas M$ se pueden descargar desde aquí. Para entrar al programa hay que ejecutar el acceso directo a Rgui.exe situado en el menú inicio.

Sistemas GNU/Linux

En Ubuntu la instalación se realiza mediante la aplicación Synaptic, buscando la librería r-base y aceptando la instalación de todas las librerías propuestas por el programa para resolver las dependencias. El programa se ejecuta desde una ventana Terminal tecleando simplemente

>R

Nota: ">" es el simbolo del sistema. asi que no hay que teclearlo

Una interfaz gráfica (similar al Rgui disponible para sistemas M$) es r-gnome, se instala desde Synaptic; para lanzar R con la interfaz grafica r-gnome debemos escribir en una ventana Terminal:

>R --gui=GNOME

Lenguaje encarado a objetos

Los lenguajes estadísticos clásicos (SAS, Stata, SPSS) parten de unos datos a los que se aplican distintos procedimientos para obtener unos resultados.

R no actua de la misma forma. En R existen los denominados objetos, sobre los que se aplican funciones para transformarlos en otros objetos, sobre quienes, a su vez, pueden aplicarse nuevas funciones para obtener nuevos objetos y asi sucesivamente.

Una función es cualquier código de R al que se le pasan uno o más argumentos para transformar uno o más objetos.

Por lo tanto en R no existen propiamente datos y resultados, sino distintos tipos de objetos y una colección de operadores y funciones para modificarlos.

Los nombres de los objetos pueden contener mayúsculas, minúsculos, números y puntos. La única restricción es que deben empezar por una letra.

Si a un objeto no se le asigna ningún nombre, simplemente será impreso en pantalla y se perderá. Es muy interesante, pues, guardar con un nombre todos los objetos creados (u obtenidos como resultado de la aplicación de funciones sobre otros objetos) para poder usarlos así en cálculos posteriores, lo cual constituye la principal ventaja del lenguaje S sobre otros lenguajes estadísticos

Convenciones básicas de la sintaxis

• Las mayúsculas y minúsculas representan distintos caracteres.
• Los espacios no representan nada y R los ignora.
• Los missings se representan por las siglas NA si los datos son numéricos, o por un espacio blanco si se trata de caracteres.
• El valor null indica que un objeto no contiene datos.
• El signo # provoca que la linea no se evalue, y por lo tanto permite realizar comentarios.
• Los decimales deben indicarse con el signo . .
• Los parentesis ( ) sirven para especificar los argumentos de las funciones.
• Los diferentes argumentos de una función deben separarse con , .
• Cuando un argumento de una función es una fórmula, se especifica con vardep~verindep1+varindep2... .
• Los corchetes [ ] sirven para especificar una posición dentro de un objeto.
• Los claudators {} se utilizan para abrir y cerrar procesos condicionales (if, for...), utilizados especialmente en la construcción de nuevas funciones.
• Todas las funciones son anidables, de tal forma que un argumento de una función puede ser a su vez una función.
• El signo $ se usa para referirse a una parte concreta y diferenciada de un objeto. Las distintas partes de un objeto pueden conocerse con la función names(nombreobjeto).
• Para asignar un nombre a un objeto se utiliza la combinación <-. Esta instrucción es fundamental.

El sistema de librerías

Una librería o paquete es un conjunto de funciones sobre una temática común. El sistema R-base trae las librerías diseñadas por el núcleo central de desarrolladores (core). Sin embargo, el modelo descentralizado de producción de R conlleva que muchas funciones interesantes no se instalen con el paquete básico sino que estén contenidas en paquetes o librerías externas que hay que instalar por separado. Aunque son descarganles desde muchos sitios de internet, se recomienda su descarga desde el repositorio Comprehensive R Archive Network (CRAN) porque se garantiza que su funcionamiento ha sido revisado y aprobado por el núcleo de desarrolladores.

Existen en R cuatro tipos de librerías:

1. las que se instalan con el sistema R-base y se cargan por defecto al inicio (no requieren pasos adicionales para ser utilizadas)
2. las que se instalan con el sistema R-base y no se cargan por defecto al inicio (requieren ser cargadas para ser utilizadas)
3. las que no se instalan con el sistema R-base y hay que bajar del CRAN (requieren instalación desde CRAN y carga)
4. las que no se instalan con el sistema R-base y no están disponible en el CRAN. (Disponibles en ficheros .zip que se encuentran habitualmente en páginas personales de internet)

Nota: La nomenclatura utilizada para los tipos de librerias es propia y no estándar.

Instalación de librerias adicionales en sistemas M$ (tipo 3)

Desde el programa (Rgui) existe un menú para instalar paquetes directamente desde el repositorio CRAN sin tenerlos que descargar previamente de forma manual. También se pueden instalar paquetes a partir de ficheros .zip a través del menú.

Instalación de librerias adicionales en sistemas GNU/Linux. (tipo 3)

Al igual que el sistema base, cualquier paquete de R puede instalarse desde Synaptic buscando la palabra CRAN.

Los paquetes se guardaran automaticamente en la ubicación:

/usr/lib/R/library/nombrepaquete

También se pueden instalar paquetes (tipo 4) manualmente descomprimiendo los ficheros .zip en dicha ubicación .

Cargado de librerías

Cuando se inicia R, sólo las funciones contenidas en determinadas librerías del núcleo básico están disponibles para ser utilizadas (tipo 1) . Incluso algunas librerías del núcleo básico no están cargadas al iniciar para ahorrar recursos (tipo 2). Para cargar una librería (tipos 2 y 3), hay que utilizar la función library(). Por ejemplo para cargar la librería 'survival' (para el análisis de supervivencia), que se instala con el paquete básico pero no se carga al inicio, habría que escribir:

library(survival)

Algunas librerías útiles

Las siguientes librerias *(tipo 2) y **(tipo 3) no se cargan directamente al iniciar R, y contienen funciones que pueden resultar interesantes:

• foreign*: Permite importar bases de datos en otros formatos (SPSS, SAS, Excel, Stata etc...).

• survival*: Funciones de análisis de supervivencia.

• car**: Recodificación de variables.

• xtable**: Permite exportar tablas a formato LaTeX.

• Rcmdr**: Presenta una interfaz gráfica más amigable para la edición de datos y ejecución de comandos desde menús. Puede ser útil para nuevos usuarios para facilitar la curva de aprendizaje, pero disminuye la gran flexibilidad del software.

El paquete Traba (tipo 4)

En el Trabachat se han desarrollado un conjunto de funciones agrupadas en el paquete Traba. El link anterior permite su descarga en formato .zip para poderlo instalar.

Se puede descargar también el código fuente para ver o modificar la sintaxis de las distintas funciones que contiene.

Editores de sintaxis

El programa R (a diferencia de SPSS, SAS o Stata) no cuenta con una ventana para escribir la sintaxis.

Cualquier editor de textos, como el Notepad de M$, sirve para almacenar los comandos, pero es recomendable trabajar con otros editores específicos que permitan, entre otros

• Colorear la sintaxis para una fácil detección de errores
• Enviar las instrucciones directamente para ser evaluadas en R (evitando asi tener que copiar/pegar cada instrucción desde el editor)

WinEdt

WinEdt es un software propietario, disponible sólo para plataformas M$, que permite la edición de distintos tipos de lenguajes, entre ellos Stata y R. Para poder utilizar el editor debe disponerse de licencia, y descargarse la librería de R RWinEdt, disponible en el CRAN.

WinEdt permite colorear la sintaxis de R, marcando las funciones correctas, y coloreando los paréntesis para evitar errores de apertura y cierre. Así mismo permite mandar una región de la sintaxis para su evaluación en R.

GNU/Emacs y XEmacs

Emacs es un anciano y robusto software para la edición de distintos lenguajes, con la posibilidad de ejecutar distintos programas en sus buffers lo cual permite evaluar directamente las expresiones dentro del propio programa. Por lo tanto Emacs es mucho más que un editor y coloreador de sintaxis, en algunos casos es poco menos que un estilo de vida.

El desarrollo de Emacs se ha dividido principalmente en dos proyectos: el programa original GNU Emacs (mas en consonacia con el movimiento software libre) y una ramificación del mismo XEmacs (mas en consonancia con el movimiento OpenSource).

Para la edición de lenguajes estadísticos en cualquiera de estos 2 programas debe usarse un modo de Emacs llamado Emacs Speaks Statistics (ESS).

Para sistemas M$ es aconsejable la instalación de XEmacs, cuyo instalador binario puede descargarse directamente desde la página web del proyecto, pues permite incorporar directamente cualquier modo, entre ellos Emacs Speaks Statistics (ESS).

En sistemas GNU/Linux se aconseja la instalación de GNU Emacs por su mayor facilidad de integración con R en estos sistemas. En Ubuntu y a través de Synaptic se puede instalar fácilmente tanto GNU Emacs como el modo Emacs Speaks Statistics (ESS), sin necesidad de mayores configuraciones.

Tipos de datos

Hay distintos tipos de datos en R, algunos de los cuales son:

• logical: sólo permite los valores T o F.
• numeric: sólo pueden contener números reales.
• character: permite introducir tanto número como cualquier caracter.

Para saber qué tipo de datos contiene un objeto:

mode(nombreobjeto)

Tipos de objetos

Los distintos tipos de datos pueden estar contenidos en los muchos tipos distintos de objetos que existen en R. Para conocer de qué tipo es un objeto:

class(nombreobjeto)

En R hay muchos tipos de objetos y con muchas propiedades distintas. Para una aproximación más detallada, consultar los distintos manuales relacionados al final de este artículo. A modo de introducción, algunos tipos de objetos son:

Vectores

Cadenas unidimensionales (es decir una sola columna o fila) de un tipo único de valores (númericos, caracteres, etc.).

Hay distintas funciones para crear vectores, la más simple de la cuales es c(), que simplemente crea un vector a partir de los valores que se introducen:

c(1,2,3,4,5) 

La instrucción anterior genera un vector que contiene 5 números (1 a 5), pero lo único que hará es imprimir en pantalla esta secuencia. Para guardar el vector como un objeto debemos escribir:

prueba<-c(1,2,3,4,5)

Ahora tenemos un objeto llamado prueba que contiene el vector de 5 números.

Una forma alternativa de obtener el mismo vector prueba es mediante una secuencia (sequence) de valores:

prueba<-seq(1:5)

Si por ejemplo ejecutamos:

prueba2<-prueba*2

obtendremos un objeto llamado prueba2 que corresponderá a un vector con los valores: 2, 4, 6, 8, 10.

La expresión:

prueba3<-prueba*prueba2

generará el vector prueba3 con los valores 2, 8, 18, 32, 50.

Para generar vectores con números repetidos se usa la función rep:

prueba4<-c(seq(1:5),rep(3,7))

que generará un vector con los valores 1, 2, 3, 4, 5, 3, 3, 3, 3, 3, 3, 3.

Para extraer un valor concreto de un vector, por ejemplo el que esté en tercera posición:

prueba3[3]

devolverá por pantalla el valor 18. Lógicamente este valor se podría haber asignado a un objeto.

Para saber el número de valores de un vector:

length(prueba3)

Matrices

Estructura bidimensional donde todos los datos deben ser del mismo tipo.

Para crear una matriz debemos usar las siguientes funciones:

• column bind

prueba4<-cbind(rep(72,5),c(2,7,9,3,5))

generará la siguiente matriz por la unión de las dos columnas (cbind) especificadas: rep(72,5) (repite 5 veces el valor 72) y c(2,7,9,3,5)

72	2
72	7
72	9
72	3
72	5

Si queremos unir filas en lugar de columnas:

• row bind

prueba5<-rbind(seq(1:5),c(2,7,9,3,5))

1	2	3	4	5
2	7	9	3	5

Para recuperar un valor concreto de la matriz anterior, debe especificarse, como es habitual, dicha posición entre corchetes. En primer lugar se marca la posición de la fila y en segunda lugar la posición de la columna.

prueba6<-prueba5[2,3]

Con la función anterior, asignamos el valor 9 al objeto prueba6.

La función:

prueba7<-prueba5[2,3:5]

asignará al objeto prueba7 un vector con los valores 9, 3, 5, correspondientes a las columnas 3, 4 y 5 de la segunda fila.

Para conocer las dimensiones de una matriz usaremos:

dim(prueba5)

que devolverá el valor 2 (filas) 5 (columas).

Para asignar nombre a las filas y a las columnas:

dimnames(prueba5)<-list(c("Fila1","Fila2"),c("Col1","Col2","Col3","Col4",Col5"))

Factores

Este tipo de objeto es fundamental para el análisis estadístico pues es la forma como se tratan las variables categóricas.

Para crear un factor:

factor1 <-factor(c(1,1,1,2,2,3,3,3,3,2,2,1,1,2,2,1,1),labels=c("Bajo","Medio","Alto"))

La función:

mode(factor1)

devolverá numeric porque los datos introducidos son numéricos, pero si usamos la función:

class(factor1)

devolverá factor indicando que el objeto es del tipo Factor.

Para saber los distintos niveles de un factor, que equivaldrían a sus etiquetas en otros sistemas:

levels(factor1)

Listas

Las listas son objetos genéricos que corresponden a colecciones de distintos objetos que pueden ser de tipos distintos.

Data Frames

Como una matriz pero pudiendo contener distintos tipos de datos. Es lo que clásicamente se ha llamado base de datos.

Si disponemos por ejemplo de dos vectores (que se pueden entender como variables) y queremos unirlos en un data frame:

base<-as.data.frame(cbind(vector1,vector2))

Otros tipos de objetos

Muchas funciones de R producen como resultado un objeto de un tipo particular y único de esa función, aunque en realidad se trata de listas que contienen como componentes los disitntos apartados de los resultados. PAra referirse a cada una de las partes de una lista se usa el caracter $.

Por ejemplo la función glm, que ajusta modelos de regresión, produce un objeto de tipo glm que contiene los distintos resultados.

Captura de datos externos

La forma más habitual de obtener datos para su análisis en R es su captura desde otros formatos, habitualmente M$ Excel, SPSS o M$ Access.

Cuando R captura una base de datos externos, en la mayoría de los casos detectará correctamente el tipo de datos que contiene, asignando como factores las variables categóricas y como vectores las variables cuantitativas.

Las mayoría de formatos externos se pueden capturar con funciones contenidas en la librería foreign (tipo 2), mientras que para la captura de bases de datos relacionales (como M$ Access) es necesaria la libreria RODBC (tipo 3).

Captura de SPSS

El siguiente ejemplo muestra cómo capturar una hoja de datos de SPSS:

library(foreign)
read.spss("C:/rutacompleta/nombrearchivo.sav")

Esta instrucción, sin embargo, sólo nos imprimirá los datos en pantalla. Para conseguir capturar la base de datos como un objecto de R:

library(foreign)
base1<-read.spss("C:/rutacompleta/nombrearchivo.sav")

Así tendremos un objeto llamado base1 del tipo data.frame. Para saber qué variables contiene:

names(base1)

Para referirse a una variable en concreto de la base de datos:

base1$nombrevar

Para evitar tener que indicar cada vez el nombre del objeto base de datos para referirse a una de sus variables, podemos usar la instrucción:

attach(base1)

Y entonces ya nos podremos referir a sus variables sin usar el prefijo de la base de datos:

nombrevar

Ya nos mostrará directamente el contenido de la variable nombrevar.

Captura de M$ Excel

R no permite capturar directamente los archivos .xls de M$ Excel. Para importar un archivo de M$ Excel debemos primero guardarlo con la extensión .dbf o .csv desde el propio M$ Excel, y luego ejecutar las funciones:

library(foreign)
base2<-read.dbf("C:/rutacompleta/nombrearchivo.dbf")
attach(base2)

o bien:

library(foreign)
base3<-read.csv("C:/rutacompleta/nombrearchivo.csv")
attach(base3)

Captura de M$ Access

El siguiente código muestra como capturar una tabla o una consulta desde una base de datos relacional creada con M$ Access:

library(RODBC)
canal<-odbcConnectAccess("c:/rutacompleta/nombrearchivo.mdb")
base4<-sqlQuery(canal,"select * from Hoja1")
attach(base4)

donde Hoja1 es el nombre de la tabla o consulta que se desea importar.

Manipulación de bases de datos

La función más simple para manipular una base de datos es la función subset, que creará un nuevo objeto que contendrá sólo parte de los datos del Data Frame original:

hombres<-subset(base1,sexo==1)

Sin embargo, hay muchas otras funciones para la manipulación de bases de datos. A continuación se muestran algunas de ellas mediante un ejemplo real que se puede reproducir paso a paso.

Se parte de dos tablas distintas en formato .dbf exportadas desde M$ Excel, y que se pueden descargar (AÚN NO!!!):

• La tabla bajas.dbf contiene un registro para cada una de las bajas laborales concedidas durante 2004 en una empresa X (a leer en catalán). Puede haber más de una baja por individuo. La base contiene la identificación del trabajador y las fechas de inicio y final de la baja.

• La tabla plantilla.dbf contiene los datos referentes a todos los trabajadores de la empresa X, con un sólo registro por individuo. Contiene la identificación del trabajador (permite el cruce con bajas.dbf), su fecha de nacimiento, el sexo y sus datos laborales (tipo de contrato, lugar de trabajo y categoría laboral).

En ambas bases de datos se han eliminado todas las variables que pudieran identificar a los individuos.

El primer procedimiento a realizar es la captura de las bases de datos mediante el procedimiento explicado anteriormente:

baixes<-read.dbf("c:/rutacompleta/baixes.dbf)
plantilla<-read.dbf("c:/rutacompleta/baixes.dbf)

Exploración de datos

Variables cuantitativas

R tiene distintas funciones para describir variables cuantitativas:

media<-mean(nombrevariable)
mediana<-median(nombrevariable)
desviacion<-sd(nombrevariable)
maximo<-max(nombrevariable)
minimo<-min(nombrevariable)

y muchas otras.

Para conocer el número de missings de una determinada variable podemos crear un subset de la variable que contenga sólo los valores missings (NA), y luego contar cuantos valores contiene ese subset:

length(subset(var,is.na(nombrevariable)==T))

La función Explore del paquete Traba devuelve todos estos índices y algunos más, de forma parecida al procedimiento Explore de SPSS.

Variables categóricas

Una variable categórica, o un Factor en el lenguaje de R, se explora mediante una distribución de frecuencias.

Para obtener el conteo de las distintas categorías de una variable usaremos:

tabla1<-table(nombrevariable)

Para obtener los porcentajes, usaremos:

prop.table(tabla1)

o para crearlo directamente desde la variable original:

prop.table(table(nombrevariable))

Para crear tablas de contingencia se usa el mismo sistema. Con:

tabla2<-table(nombrevariable1,nombrevariable2)

obtenemos los conteos correspondientes a cada una de las celdas.

Para saber los porcentajes:

prop.table(tabla2)

o directamente:

prop.table(table(nombredevariable1,nombredevariable2))

Esta función nos dará el porcentaje de cada celda respecto al total. Para obtener el porcentaje respecto al total de la cada fila:

prop.table(table(nombredevariable1,nombredevariable2),1)

Y el porcentaje respecto al total de cada columna:

prop.table(table(nombredevariable1,nombredevariable2),2)

Si asignamos las funciones anteriores a objetos, podremos capturar un valor concreto que nos pueda interesar:

prop1<-prop.table(table(nombredevariable1,nombredevariable2),2)
valor21<-prop1[2,1]

o directamente:

valor21<-prop.table(table(nombredevariable1,nombredevariable2),2)[2,1]

A continuación se detalla un ejemplo real de como obtener una tabla de contingencia a partir de dos factores:

var1<-factor(c(1,1,1,1,1,3,3,3,3,3,2,2,2,2,1,1,1,2,2,3,3,1),labels=c("Bajo","Medio","Alto"))
var2<-factor(c(1,2,2,2,1,1,1,2,2,2,1,1,1,1,1,2,2,2,1,1,2,1),labels=c("Hombre","Mujer"))
prop.table(table(var1,var2),2)

De esta forma obtendremos una tabla con el sexo en las columnas y el nivel socio-económico en las filas, y en cada celda habrá el porcentaje sobre el total de la columna.

Pruebas estadísticas clásicas

t de Student

La comparación de medias mediante la prueba de t de student se realiza con la función t.test, cuya sintaxis es la siguiente:

t.test(variablecuantitativa~variablebinaria)

Por ejemplo:

cuanti<-c(3,3,4,4,6,6,8,8,9,10,1,1,5,21,23,4,4,5,65,12)
catego<-c(2,1,1,1,2,2,1,1,1,1,1,1,2,2,2,1,1,2,1,2)
t.test(cuanti~catego)

Si en lugar de imprimir el resultado del test por pantalla, lo que queremos es guardarlo como un objeto:

test<-t.test(cuanti~catego)

Ejecutando:

names(test)

tendremos una relación de todas las partes del test. Si, por ejemplo, nos interesa solo el valor p:

p<-test$p.value

o directamente:

p<-t.test(cuanti~catego)$p.value

De esta forma, el objeto p contiene el valor p de la prueba, y lo podemos utilizar para realizar cálculos posteriores.

Ji-cuadrado

La prueba de Ji-cuadrado debe aplicarse sobre tablas de conteos creadas con los procedimientos explicados anteriormente. Por ejemplo, para obtener el valor del estadístico ji-cuadrado y guardarlo en el objeto ji:

var1<-factor(c(1,1,1,1,1,3,3,3,3,3,2,2,2,2,1,1,1,2,2,3,3,1),labels=c("Bajo","Medio","Alto"))
var2<-factor(c(1,2,2,2,1,1,1,2,2,2,1,1,1,1,1,2,2,2,1,1,2,1),labels=c("Hombre","Mujer"))
ji<-chisq.test(table(var1,var2))$statistic

Regresión

R posee una potente función, llamda glm, que permite ajustar cualquier modelo lineal generalizado. Los argumentos mínimos a especificar son la fórmula del modelo, la família a la que pertenece y la función de linkado. La siguiente tabla muestra las famílias y links más comúnmente utilizados en ciencias de la salud:

factores...

Los resultados de un modelo se almacenan en un objeto de tipo glm que.....

Gráficos

R es probablemente el sistema más potentes y flexible de creación de gráficos que existe actualmente. Existen tantas opciones que es imposible resumirlas en un espacio reducido.

La principal función para la creación de gráficos es plot, por lo que se aconseja consultar su ayuda cuando se quiera crear un gráfico:

?plot

y para consultar los múltiples parámetros que se pueden pasar a un gráfico:

?par

La función plot se denomina de alto nivel porque cualquier gráfico debe empezar con esta instrucción y sus argumentos indispensables. Esto abrirá una ventana de gráfico y empezará su creación.

Existen además muchas funciones gráficas denominadas de bajo nivel que van a modificar el gráfico presente en una ventana gráfica generada previamente por la función plot.

Además R permite guardar directamente un gráfico generado en un archivo de imagen para su uso posterior. Esto se puede hacer, por ejemplo, con la función jpeg.

Uno de los argumentos de la instrucción genérica plot es el tipo de gráfico que se desea crear. Sin embargo, para facilitar la creación de gráficos (aún a cuesta de perder flexibilidad) existen muchas funciones que crean gráficos concretos, como por ejemplo barplot, hist o boxplot:

boxplot(nombrevariable,color="orange")

La instrucción anterior creará un boxplot de color naranja de la variable nombrevariable

A continuación se muestran algunos ejemplos reales de gráficos creados en R, sin comentar el código::

• Ejemplo 1: Gráfico de barras con el sexo de los pacientes alérgicos al huevo vacunados con triple vírica:

jpeg("c:/documents and settings/acrida/escritorio/TV Ou/gràfics/sexe.jpg")
barplot(prop.table(table(sexe)),ylim=c(0,1),col="orange",width=.3,xlim=c(0,.7),names.arg=F,cex.axis=1.8)
dev.off()

• Ejemplo 2: Curva epidémica (histograma) de una toxinfección alimentaria:

jpeg("t:/rs_cdp/Epidemiologia/francesc/Brots/GEA SS Hospi/corba.jpg",width = 600, height = 600)
hist(base,breaks=seq(ISOdatetime(2004,11,6,0,0,0),ISOdatetime(2004,11,12,0,0,0),by="12      hours"),right=F,freq=T,col="orange",axes=F,xlab="Inici de simptomes",ylab="Casos",main="Corba epidemica")
       axis.POSIXct(base,side=1,at=seq(ISOdatetime(2004,11,6,0,0,0),ISOdatetime(2004,11,12,0,0,0),by="days"),format="%d-%b")
axis(side=2,at=seq(0,3,by=1))
polygon(c(ISOdatetime(2004,11,6,12,0,0),ISOdatetime(2004,11,6,12,0,0),ISOdatetime(2004,11,7,0,0,0),ISOdatetime(2004,11,7,0,0,0)),c(0,1,1,0),col="blue")
dev.off()

Programación de funciones

Otro punto fuerte de R es la posibilidad de crearse funciones propias, como es el caso de las contenidas en el paquete Traba. A continuación se explica el procedimiento de creación de una función mediante el ejemplo de las funciones intervals y montehall' de dicho paquete,

El código completo de ambas funciones puede encontrarse en el código fuente del paquete Traba.

El ejemplo de la función interval

Esta función va a calcular los intervalos de confianza para los parámetros de distintos modelos de regresión. Aunque existen funciones parecidas en otros paquetes, la intención al crear la función fue centralizar en una única expresión los cálculos para cualquiera de los modelos más usada en ciencas de la salud. Y aprender, claro :)

El primer paso consiste en crear un fichero con el nombre interval.R con cualquier editor (Emacs, WinEdt o cualquier editor de textos).

A continuación se irán mostrando las distintas partes del código con una breve descripción de los procedimientos que realiza, para que cada cual vaya descubriendo por si sólo la funcionalidad y sintaxis de las distinas partes:

1-Definimos el nombre de la función, los parámetros que deberá contener y sus valores por defecto en caso que los haya:

interval<-function(model,digits=2,nivell=0.95,alfa=F) {

2-Creamos algunas variables internas a partir de los parámetros recibidos

nivell2<--round(qnorm((1-nivell)/2),digits=2)
co<-summary(model)$coef
tot<-c()

3-Comprobamos que estamos recibiendo un modelo de regresión válido:

if(class(model)[1]!="lm" & class(model)[1]!="glm" & class(model)[1]!="coxph")
{
ex.warn <- options()$warn 
options(warn = 0)
warning("El objeto no es un modelo de regresion valido",call.=F)
options(warn = ex.warn)
}

4-Cálculos si hemos recibido un modelo de regresión lineal:

if(class(model)[1]=="lm")
{
df<-model$df
p<-(1-nivell)/2
punt<-round(co[,1],digits=digits)
inf<-round(confint(model,level=nivell)[,1],digits=digits)
sup<-round(confint(model,level=nivell)[,2],digits=digits)
tot<-cbind(punt,inf,sup)
dimnames(tot)<-list(row.names(summary(model)$coef),c("B","inf","sup"))
}

5-Cáculos si hemos recibido un modelo de riesgos proporcionales de Cox:

if(class(model)[1]=="coxph")
{
punt<-round(exp(co[,1]),digits=digits)
inf<-round(exp(co[,1]-nivell2*co[,3]),digits=digits)
sup<-round(exp(co[,1]+nivell2*co[,3]),digits=digits)
tot<-cbind(punt,inf,sup)
dimnames(tot)<-list(row.names(summary(model)$coef),c("HR","inf","sup"))
alfa<-T
}

6-Cálculos si hemos recibido un glm (logística, logbinomial, poisson):

if(class(model)[1]=="glm")
{
if(summary(model)$family[1]=="gaussian")
{
ex.warn <- options()$warn 
options(warn = 0)
warning("Utiliza la funcion lm en lugar de glm para los modelos de regresion lineal",call.=F)
options(warn = ex.warn)
}
else
{
punt<-round(exp(co[,1]),digits=digits)
inf<-round(exp(co[,1]-nivell2*co[,2]),digits=digits)
sup<-round(exp(co[,1]+nivell2*co[,2]),digits=digits)
tot<-cbind(punt,inf,sup)
if(summary(model)$family[1]=="binomial" & summary(model)$family[2]=="logit")
{
dimnames(tot)<-list(row.names(summary(model)$coef),c("OR","inf","sup"))
}
if(summary(model)$family[1]=="binomial" & summary(model)$family[2]=="log")
{
dimnames(tot)<-list(row.names(summary(model)$coef),c("RR","inf","sup"))
}
if(summary(model)$family[1]=="poisson")
{
dimnames(tot)<-list(row.names(summary(model)$coef),c("IRR","inf","sup"))
}
}
}

7-Preparamos e imprimimos los resultados:

if(!is.null(tot) & alfa==T)
{
return(tot)
}
if(!is.null(tot) & alfa==F)
{
tot<-tot[2:dim(tot)[1],]
if(length(tot)>3)
{
return(tot)
}
else
{
tot<-t(tot)
row.names(tot)<-row.names(summary(model)$coef)[2]
return(tot)
}
}
}

A continuación se presenta el código completo. Su edición con Emacs o WinEdt facilita su comprensión al colorear las funciones y los parentésis y claudators.

interval<-function(model,digits=2,nivell=0.95,alfa=F) {
nivell2<--round(qnorm((1-nivell)/2),digits=2)
co<-summary(model)$coef
tot<-c()
if(class(model)[1]!="lm" & class(model)[1]!="glm" & class(model)[1]!="coxph")
{
ex.warn <- options()$warn 
options(warn = 0)
warning("El objeto no es un modelo de regresion valido",call.=F)
options(warn = ex.warn)
}
if(class(model)[1]=="lm")
{
df<-model$df
p<-(1-nivell)/2
punt<-round(co[,1],digits=digits)
inf<-round(confint(model,level=nivell)[,1],digits=digits)
sup<-round(confint(model,level=nivell)[,2],digits=digits)
tot<-cbind(punt,inf,sup)
dimnames(tot)<-list(row.names(summary(model)$coef),c("B","inf","sup"))
}
if(class(model)[1]=="coxph")
{
punt<-round(exp(co[,1]),digits=digits)
inf<-round(exp(co[,1]-nivell2*co[,3]),digits=digits)
sup<-round(exp(co[,1]+nivell2*co[,3]),digits=digits)
tot<-cbind(punt,inf,sup)
dimnames(tot)<-list(row.names(summary(model)$coef),c("HR","inf","sup"))
alfa<-T
}
if(class(model)[1]=="glm")
{
if(summary(model)$family[1]=="gaussian")
{
ex.warn <- options()$warn 
options(warn = 0)
warning("Utiliza la funcion lm en lugar de glm para los modelos de regresion lineal",call.=F)
options(warn = ex.warn)
}
else
{
punt<-round(exp(co[,1]),digits=digits)
inf<-round(exp(co[,1]-nivell2*co[,2]),digits=digits)
sup<-round(exp(co[,1]+nivell2*co[,2]),digits=digits)
tot<-cbind(punt,inf,sup)
if(summary(model)$family[1]=="binomial" & summary(model)$family[2]=="logit")
{
dimnames(tot)<-list(row.names(summary(model)$coef),c("OR","inf","sup"))
}
if(summary(model)$family[1]=="binomial" & summary(model)$family[2]=="log")
{
dimnames(tot)<-list(row.names(summary(model)$coef),c("RR","inf","sup"))
}
if(summary(model)$family[1]=="poisson")
{
dimnames(tot)<-list(row.names(summary(model)$coef),c("IRR","inf","sup"))
}
}
}
if(!is.null(tot) & alfa==T)
{
return(tot)
}
if(!is.null(tot) & alfa==F)
{
tot<-tot[2:dim(tot)[1],]
if(length(tot)>3)
{
return(tot)
}
else
{
tot<-t(tot)
row.names(tot)<-row.names(summary(model)$coef)[2]
return(tot)
}
}
}

Una vez ejecutada la función, ya se puede utilizar:

interval(glm(vardep~varindep1+varindep2))

Para no tener que ejecutar cada vez la función, se recomienda incluirla en un paquete (como se hizo con el paquete Traba). El procedimiento para hacerlo se explica más abajo.

El ejemplo de la función montehall

Esta función reproduce la paradoja de Monte Hall. Se trata de un concursante que está en un concurso donde hay tres puertas. Detrás de una de ellas hay un coche y detrás de las otras dos una cabra (una en cada una de las otras dos puertas). Se supone que al concursante le hace más ilusión ganar el coche que la cabra. El presentador le hace escoger una puerta, pero antes de abrirla lo que hace es abrir una de las otras dos puertas (teniendo en cuenta que el presentador sabe dónde está el coche y abre siempre una puerta con una cabra). Luego le pregunta al concursante qué quiere hacer, si quedarse en la puerta que ha escogido o cambiar a la otra puerta que queda cerrada.

Intuitivamente la gente suele decir que es indiferente pues hay un 50% de probabilidad de acertar pues quedan dos puertas. Pero se demuestra que si se cambia de puerta, se gana en el 66,6% de los casos. El siguiente código lo demuestra, y además ilustra el uso de loops (función for) y funciones de generación de número aleatorios en R:

montehall<-function(repeticions,canvi)
{
cotxe<-trunc(runif(repeticions,1,3.99999999999999999999999999999999))
eleccio<-trunc(runif(repeticions,1,3.99999999999999999999999999999999))
res<-0
for(i in 1:repeticions)
{
if(cotxe[i]==eleccio[i])
{
if(cotxe[i]==1)
{
obrir<-trunc(runif(1,2,3.99999999999999999999999999999999))
}
if(cotxe[i]==2)
{
obrir<-trunc(runif(1,1,2.99999999999999999999999999999999))
if(obrir==2)
{
obrir=3
}
}
if(cotxe[i]==3)
{
obrir<-trunc(runif(1,1,2.99999999999999999999999999999999))
}
}
else
{
trampa<-seq(1:3)
obrir<-subset(trampa,trampa!=cotxe[i] & trampa!=eleccio[i])
}
if(canvi==T)
{
trampa2<-seq(1:3)
eleccio2<-subset(trampa2,trampa2!=obrir & trampa2!=eleccio[i])
}
else
{
eleccio2<-eleccio[i]
}
if(cotxe[i]==eleccio2)
{
res<-res+1
}
}
return(round(100*res/repeticions,digits=2))
}

Para realizar el experimento 1000 veces, cambiando siempre de puerta:

montehall (1000,T)

El resultado debería aproximarse a 66,6.

Documentación de funciones

Cada función debe ir acompañada de un archivo de ayuda que se presentará cuando el usuario teclee:

?nombrefunción

La forma más sencilla de hacerlo es crear un nuevo fichero llamado nombrefunción.Rd (por ejemplo intervals.Rd) y escribir directamente la ayuda. Se aconseja que los ficheros de ayuda contengan siempre los mismos apartados (nombre de la función, parámetros, resultados, etc.), como se puede ver en el código fuente del paquete Traba.

Construcción de paquetes

Una vez tenemos el código fuente tanto de las funciones como de las ayudas que deseamos incluir en el paquete, debemos efectuar una serie de procedimientos para empaquetarlo de forma que sea directamente instalable en R. Para hacerlo.......

Sistema de ayuda

Ayuda disponible en el propio sistema

Desde R podemos consultar la ayuda de todas las funciones contenidas en los paquetes que tengamos instalados. Hay tres formas de hacerlo:

• Si conocemos el nombre de la función que queremos utilizar pero queremos consultar sus argumentos u otros detalles:

?nombrefunción

• Si no conocemos el nombre de la función que queremos utilizar pero sabemos el nombre de la técnica:

help.search("búsquedaarealizar")

• Si queremos navegar por todas las funciones que están disponibles en nuestro sistema:

help.start()

Lo cual abrirá nuestro navegador y nos permitirá navegar por todos los paquetes y funciones de los que disponemos. Además desde la misma ventana podremos accedes a distintos manuales oficiales.

Ayuda externa

La comunidad R posee una potente lista de correo que permite la rápida solución de dudas y problemas que puedan surgir durante el uso del programa. Sin embargo, sólo de forma excepcional será necesario realizar una pregunta, pues seguro que alguien se la ha planteado antes que nosotros. Por tanto, una forma rápida de resolver una duda es buscando en Google la palabra clave r-help (identificativo de la mayoría de listas de correo de la comunidad) seguido de la duda que queremos resolver.

Por ejemplo, si no sabemos cómo realizar una regresión logística y aplicamos la búsqueda sugerida, encontramos más de 3.000 resultados.

Los usuarios de Firefox pueden configurarse un marcador Quick search con la palabra clave r y la descripción http://www.google.com/search?q=r-help%20%s, lo que permitirá realizar búsquedas en r-help escribiendo simplemente en la barra de direcciones r duda.

Este sistema puede ser muy útil cuando desconocemos qué técnica estadística aplicar, pero si ya sabemos por ejemplo que queremos usar una regresión logística pero no sabemos qué función utilizar per realizarla, lo más eficiente es buscar en R Help Center, que contiene las funciones y sus correspondientes ayudas de todos los paquetes disponibles en el CRAN.

Otros manuales

Manuales oficiales (en ingés)

La web de Victor Moreno con distintos manuales y clases en castellano.

Otra clase en castellano.