RStudio

En las entradas precedentes del blog  he usado, hasta ahora, R-Commander como interfaz que simplifica nuestra relación con R, especialmente al principio. Y sigo pensando que es un buen programa para principiantes, y continuaré explicando las opciones que ofrece.  Tengo además una deuda de gratitud con ese programa, porque mis primeros pasos con R habrían sido mucho más difíciles si no hubiera descubierto R-Commander (y el paquete R-UCA de la Universidad de Cádiz).

Sin embargo, en los últimos tiempos me he dejado convencer por otra interfaz gráfica de R, que se llama R-Studio (el enlace es la página principal del proyecto). Para presentarlo, aquí tenemos una captura de pantalla de RStudio, en una sesión de trabajo con el código de la entrada sobre ANOVA para diseños factoriales (con tres factores en este caso).

RStudio001

RStudio ofrece muchas herramientas para el usuario avezado de R, convirtiéndose en un entorno de trabajo muy productivo. Esa es la mayor diferencia con R-Commander, que tal vez está más orientado a proporcionar al usuario una experiencia  similar a la de los programas comerciales clásicos, donde los análisis están disponibles mediante un sistema de menús y cuadros de diálogo, como este:

Rcmdr-UnMenu

En RStudio no encontraremos nada de eso. No hay un menú que nos diga cómo hacer un contraste de hipótesis, un intervalo de confianza. Para eso está R. En ese sentido, RStudio da por sentado que nosotros sabemos pedirle a R que haga esos análisis. Y por eso, para un principiante, puede resultar algo más intimidatorio.

Para el usuario con algo más de experiencia, como decíamos, la situación se invierte, y RStudio se muestra como una de las mejores opciones para el trabajo diario con R. Entre otras, quiero destacar estas características:

  • la posibilidad de trabajar con varios ficheros de comandos R simultáneamente, agrupados por pestañas (en la ventana superior izquierda, en la anterior captura de pantalla).
  • el editor de comandos es, en si mismo, una de las grandes virtudes de RStudio. Es realmente un editor pensado por y para programar en R. No echaremos en falta un editor externo, al que en general teníamos que recurrir con R-Commander.
  • autocompletado de comandos. Si no recordamos la sintaxis exacta de un comando de R, basta con escribir las primeras letras, pulsar la tecla Tab, y RStudio nos ofrece información sobre ese comando. Pero también sobre las variables, ficheros y demás objetos que hayamos creado.
  • gestión de los gráficos de R mucho más eficaz que, por ejemplo, en R-Commander.
  • facilidad de gestión para, por ejemplo, instalar paquetes, consultar ayuda, etcétera.
  • si se dispone de varias versiones de R instaladas en la máquina, es fácil seleccionar la que utilizará en cada sesión.

Pero además, RStudio ofrece la posibilidad de instalación en un servidor (bajo Linux), para posteriormente acceder al programa a través de una interfaz web, desde cualquier máquina, aunque no esté instalado R. Esa interfaz web es esencialmente idéntica a la que hemos mostrado más arriba, como puede verse en esta figura:

RStudio002-InterfazWeb

Esto hace posible (aunque desde luego, no muy cómodo) acceder a RStudio desde un teléfono móvil. Y desde luego es especialmente interesante para usuarios de tablets, (tanto Ipads como Android).

Para finalizar esta primera visita, en la página del proyecto RStudio se puede ver otro producto, Shiny, del que hablaremos en una futura entrada y que permite convertir los análisis que hacemos con R en páginas web interactivas, para que cualquier pueda experimentar con esos análisis, y modificar los parámetros utilizando deslizadores, cajas de diálogo, etc.  Una muestra, en este enlace.

Shiny  es una forma de avanzar en la idea de mostrar a los demás nuestro trabajo, cada vez de forma más trasparente. Una idea que está en la base del movimiento hacia lo que se ha llamado Investigación reproducible (Reproducible Research, página en la Wikipedia). En esencia, la idea es que los trabajos de investigación, a la hora de publicar sus resultados, debería acompañar la publicación tradicional con el conjunto completo de datos experimentales y los ficheros de código del software que se haya utilizado para el análisis de esos datos.

R tiene una herramienta muy poderosa para contribuir a esta idea, llamada SWeave. Y en esta entrada no puedo dejar de mencionar que RStudio (a diferencia de R-Commander) ofrece mucha ayuda para el trabajo con SWeave. Pero ese un tema que merece ser tratado con mucho más detenimiento. En una próxima entrada volveremos sobre la investigación reproducible (y sus conexiones con la política y la economía, nada menos), sobre SWeave, y sobre la noción de Programación Literaria, que hizo posible todas estas ideas. Y (spoiler alert para profesores) sobre la forma de conectar todo esto con Moodle y usarlo en nuestras clases.

Gracias por la atención.

Anuncios

Actualizar R en Windows

Y una entrada mínima para complementar la información de la anterior, esta vez sobre la actualización en Windows. Las instrucciones, en inglés, están aquí.

Una advertencia: si, por ejemplo, has instalado R y R-Commander usando el paquete R-UCA de la universidad de Cadiz (como expliqué en esta entrada), entonces es posible que te hayas acostumbrado a comunicarte con R a través de R-Commander, y entonces apenás usarás la ventana RGui. En tal caso, y para evitar problemas, no ejecutes la instalación desde R-Commander.  Es mejor que cierres la ventana de R-Commander y ejecutes la instalación desde esa ventana RGui.

Contrastes de hipótesis y primeros pasos en la programación: estructuras If-else.

Contraste de hipótesis sobre igualdad de varianzas

Aunque esta entrada se presenta como continuación de la 14ª sesión con Rcmdr, en la que hablábamos de contrastes de hipótesis, en realidad es más bien una excusa para empezar a hablar en el blog de las estructuras de  programación de R. Si ya tienes alguna experiencia escribiendo programas (en cualquier lenguaje de programación), la lectura de esta entrada te resultará muy, muy fácil. Pero si nunca lo has hecho, entonces esta entrada es para ti, he tratado de hacerla lo más fácil de digerir posible.

En aquella entrada hicimos algún contraste de hipótesis usando, por un lado, funciones de R específicas para ese tarea (como t.test), y por otro lado, haciendo “a mano” el contraste. Es decir, calculando los valores necesarios para el contraste paso a paso. Y al final de aquella entrada dejamos sin cubrir el caso del contraste de igualdad de varianzas en dos poblaciones normales. En parte para no hacerla demasiado extensa, y en parte porque quería reservar ese caso para esta entrada. Vamos con ello, entonces. En uno de los últimos exámenes que hemos hecho, pedíamos a nuestros alumnos que hicieran este ejercicio:

Dos marcas de analgésicos a firman en su publicidad que sus productos proporcionan un
alivio muy rapido del dolor de cabeza. Para decidir cual de los dos es más rápido, hemos
seleccionado un grupo de 22 pacientes y los hemos dividido aleatoriamente en dos grupos. Cada grupo ha recibido tratamiento con uno de los dos analgésicos. Los tiempos, en minutos, que los pacientes de cada grupo tardaron en sentir aliviado su dolor fueron estos:
Grupo 1={ 9.9, 8.0, 9.5, 5.9, 12.2, 13.5, 9.6, 11.5, 12.5, 9.5, 10.3, 11.9 }
Grupo 2={ 8.2,17.3, 10.1, 10.2, 9.1, 10.5, 9.7, 9.0, 15.2, 11.6 }
Utiliza estos datos para decidir cuál de los dos medicamentos parece producir efecto más
rápidamente. Ten en cuenta que no se sabe nada sobre la posible igualdad de varianzas en las poblaciones involucradas. Utiliza un nivel de confi anza del 95% en todos los contrastes de hipótesis que realices, y calcula además el p-valor de esos contrastes.

Es un ejercicio de contraste de medias en dos poblaciones (que se asumen normales), pero, como indica la frase subrayada, y dado que se trata de muestras pequeñas, es conveniente llevar a cabo previamente un contraste de igualdad de varianzas, con hipótesis nula:

H_0=\{\sigma_1^2=\sigma_2^2\}

Que es precisamente el caso que habíamos dejado pendiente.

El estadístico que vamos a utilizar en este caso es el cociente de las cuasivarianzas muestrales:

\dfrac{s_1^2}{s_2^2}

Utilizamos R para calcular el valor del estadístico en las muestras de ambas poblaciones:

> muestra1=c(9.9, 8.0, 9.5, 5.9, 12.2, 13.5, 9.6, 11.5, 12.5, 9.5, 10.3, 11.9)
> muestra2=c(8.2,17.3, 10.1, 10.2, 9.1, 10.5, 9.7, 9.0, 15.2, 11.6)

> (s1=sd(muestra1))
[1] 2.115509
> (s2=sd(muestra2))
[1] 2.914504

> (n1= length(muestra1))
[1] 12
> (n2= length(muestra2))
[1] 10

> (Estadistico=s1^2/s2^2)
[1] 0.5268664

Si la hipótesis nula fuera cierta, y las varianzas de las dos poblaciones fuesen iguales,  la distribución muestral de este estadístico sería una F de Fisher.

Así que para calcular el p-valor del contraste basta con hacer

> (pValor=2*(pf(Estadistico,df1=n1-1,df2=n2-1)))
[1] 0.3139952

Las preguntas estándar que se hace cualquier aprendiz de Estadística (y que tienen una contestación más extensa en los apuntes de teoría) son:  ¿por qué el 2? y ¿por qué usamos pf en lugar de 1-pf? Multiplicamos por dos porque se trata de un contraste bilateral, y debemos tener en cuenta la probabilidad de las dos colas de la distribución. Y usamos pf porque estamos calculando la probabilidad correspondiente a la cola izquierda de la distribución. Y ese es el punto clave (que hace de este el ejemplo ideal para esta entrada, como veremos) y que muchas veces causa dificultades a los principiantes. ¿Cómo sabemos que hay que usar la cola izquierda y no la cola derecha? Una gráfica de la distribución de Fisher nos ayudará en la discusión:2012_13_EstadisticaBiologia-ExamenEnero-Ejercicio01-ContrasteSi la hipótesis nula fuera cierta, los valores del estadístico cercanos a uno serían los más probables (porque el estadístico vale uno cuando las cuasivarianzas son iguales). Los valores que ponen en aprietos a la hipótesis nula son los que se alejan de uno, hacia una u otra de las colas de la distribución. Y como nosotros hemos obtenido un estadístico aproximadamente igual a 0.527, que es menor que uno, debemos usar la cola izquierda. Para tratar de dejarlo aún más claro, si intercambiamos los nombres de las muestras (como si las hubiéramos ordenado al revés), tendríamos:

> muestra2=c(9.9, 8.0, 9.5, 5.9, 12.2, 13.5, 9.6, 11.5, 12.5, 9.5, 10.3, 11.9)
> muestra1=c(8.2,17.3, 10.1, 10.2, 9.1, 10.5, 9.7, 9.0, 15.2, 11.6)

> (s1=sd(muestra1))
[1] 2.914504
> (s2=sd(muestra2))
[1] 2.115509

> (n1= length(muestra1))
[1] 10
> (n2= length(muestra2))
[1] 12

> (Estadistico=s1^2/s2^2)
[1] 1.898014

Y naturalmente, \dfrac{1}{1.898014}=0.5268665. Y si usamos  esta segunda ordenación de las muestras, puesto que ahora el valor del estadístico es mayor que uno, tenemos que usar la cola derecha de la distribución, como ilustra esta figura (no es la misma F de antes; los grados de libertad también han intercambiado sus posiciones):2012_13_EstadisticaBiologia-ExamenEnero-Ejercicio01-Contraste-aY ahora usamos 1-pf  para calcular el p-valor. En cualquier caso el resultado es, naturalmente,el mismo:

> (pValor=2*(1-pf(Estadistico,df1=n1-1,df2=n2-1)))
[1] 0.3139952

Y ahora que hemos discutido esto, vamos a lo que más nos interesa en esta entrada. Hemos visto que, dependiendo del valor del estadístico (es decir, dependiendo de cómo se ordenen las muestras), tendremos que decidir si usar la cola derecha o la cola izquierda para obtener el p-valor.  ¿Podemos automatizar este proceso de decisión? En el fondo, de eso trata la programación: de automatizar las decisiones. Por eso en el próximo apartado vamos a dar nuestro primer paso en el mundo de la programación en R.

Tomando decisiones. Estructuras If-else.

Si el objetivo (confesable) de todas las sesiones precedentes era acompañar al lector que está usando R en un curso de introducción a la Estadística, a partir de esta entrada vamos a  adentrarnos en un terreno que excede ese objetivo (tan modesto, por otra parte).  Si R ha llegado a la posición que ahora ocupa, y si su radio de acción no deja de crecer, es sobre todo porque no se limita a ser un programa de Estadística más al uso, con cada vez más opciones en los menús. Ante todo, R es un lenguaje de programación, con un fuerte sesgo -desde luego- hacia el análisis de datos. Para sacarle todo el partido a R hay que aprender algo de programación. Y otro de los objetivos de este blog es allanar ese camino a quien quiera hacerlo.

Sin más preámbulos, ¿cómo se hace esto? La decisión más básica consiste siempre en elegir entre dos caminos posibles (en general, dos opciones; pero la imagen del camino ayuda a pensar).  Y el esquema de una decisión básica siempre es este:

Si se cumple cierta condición, entonces vamos por el camino A. Y si no se cumple, vamos por el camino B.

Las frases que hemos marcado en azul son las que cambiarán en cada caso concreto, porque en cada ejemplo la condición será distinta, y distintos serán también los caminos A y B. Pero la forma de la frase es la misma siempre. Para traducir una frase de esta forma  al lenguaje de R disponemos de un estructura especial que, esquemáticamente, dejando todavía sin traducir las partes en azul de la frase, funciona así:

if(se sumple cierta condicion){
  vamos por el camino A
}else{
  vamos por el camino B
}

Veamos un ejemplo de decisión más concreto. Vamos a lanzar un dado (usando la función sample de R).  Si el resultado es menor o igual que  tres, entonces tú ganas. Y si eso no se cumple (es decir, si es mayor que tres),  yo gano. He coloreado en azul las partes correspondientes para resaltar la estructura de la frase. En R esto se convierte en:

(dado=sample(1:6,1))

if(dado<=3){
decision='tu ganas'
}else{
decision='yo gano'
}

decision

La parte interesante es el grupo central de instrucciones; el fragmento de código empieza con el lanzamiento del dado (y se muestra el resultado para ayudar a entender lo que sucede a continuación). Y el código termina mostrando el resultado, la decisión que se ha tomado.  Si es la primera vez que ves una de estas estructuras if-else, te conviene ejecutar el código varias veces para hacerte una idea de como funciona.

Volvamos al contraste sobre igualdad de varianzas. ¿Cuál es la estructura de la decisión que teníamos que tomar? Es decir, ¿cuál es la condición, y cuáles son los dos caminos? Vamos a intentar escribirlo de forma que la estructura sea evidente:

Si el estadístico es menor que uno, entonces calculamos el p-valor usando la cola izquierda de la F. Y si eso no se cumple (si el estadístico es >=1), entonces calculamos el p-valor usando la cola derecha.

Traducido a R el código que continúa con el ejemplo del examen es así (repetimos la sentencia en la que se calcula el Estadísitico para que se vea cómo encaja  esto con lo anterior):

(Estadistico=s1^2/s2^2)

if(Estadistico<1){
  pValor=2*pf(Estadistico,df1=n1-1,df2=n2-1)
}else{
  pValor=2*(pf(Estadistico,df1=n1-1,df2=n2-1))
)

pValor

Al ejecutar este código se obtiene el resultado correcto ¡sea cual sea el orden en que se coloquen las muestras! Esto nos permite escribir un fichero de comandos R que calcula el p-valor correcto, sin tener que preocuparnos del orden en que se colocan las muestras. Es muy razonable decir que es un programa en lenguaje R. Muy sencillo, pero un programa al fin y al cabo.

En próximas entradas vamos a aprender más sobre estos grupos if-else, y vamos a presentar más estructuras de programación,  que nos van a permitir crear programas en R cada vez más sofisticados, que por un lado nos permitirán sacar el mayor partido posible del sistema, y por otro lado nos harán capaces de entender los programas que otros han escrito y tal vez adaptarlos a nuestras necesidades concretas.

En la próxima entrada de esta serie, y antes de avanzar más, vamos a asegurarnos de que el terreno que pisamos es firme. Para ello, nos entretener en la forma de escribir la condición que usamos en una estructura de tipo if-else.  ¿Qué significa exactamente una condición? ¿Y como se escriben en R? Pronto lo veremos.

Muchas gracias por la atención.

15ª sesión con Rcmdr: Factores y ANOVA unifactorial en R.

(Modificado el 20/01/2013 para corregir un error en la primera aparición de tapply)

Factores

En la octava sesión (se recomienda una relectura rápida de esa entrada antes de leer esta) hablamos de los datos de tipo character en R, que podían usarse, por ejemplo, para representar valores de variables cualitativas. Una situación frecuente en Estadística es esta: hemos medido una serie de valores de una variable cuantitativa X (es decir, los valores de X son números), pero esos valores aparecen agrupados de manera natural. Por ejemplo, cuando estamos midiendo la respuesta X (un número) de una serie de pacientes frente a varios tratamientos, es evidente que lo natural es agrupar los resultados según el tratamiento empleado. Una manera de hacer esto, en R, es usar un data.frame (ver la octava sesión) que contenga los valores de X, junto con los valores de una variable que indique el tipo de tratamiento empleado. Vamos a llamar

T_1, T_2, \ldots, T_k

a los distintos tratamientos. La terminología estándar en Estadística es decir que el tratamiento T es un factor, y que sus valores son los niveles del factor. Podríamos entonces pensar en recoger los resultados en un fichero de datos como este, llamado experimento-01.csv. y trabajar así, usando variables de tipo character para los factores. Pero eso tendría un impacto negativo en el rendimiento de R. Y como los factores son omnipresentes en Estadística, R ha optado por incluir un tipo especial de vector para representar los factores y sus niveles. Y, de hecho, cuando leemos un fichero como el anterior usando read.table (de nuevo, remitimos al lector a la octava sesión). el comportamiento por defecto de R es convertir las variables cualitativas en factores.

Vamos a ver esto en funcionamiento. Guardaremos el contenido del fichero en un data.frame llamado experimento, con dos campos:

  1. experimento$Respuesta, de tipo numérico que almacena los valores de la variable X (la respuesta individual de cada uno de los pacientes).
  2. experimento$Tratamiento, de tipo factor que almacena los valores de la variable T (el tipo de tratamiento empleado).

Usaremos la función read.table para leer (y mostrar) el fichero, y luego la función class para ver el tipo de objeto que corresponde a cada uno de los campos del data.frame:

(experimento=read.table(file="experimento-01.csv",header =T, sep=" "))
class(experimento$Respuesta)
class(experimento$Tratamiento)

La ventana de resultados de Rcmdr, aparte de mostrar el contenido del fichero, indica que el campo experimento$Tratamiento es, en efecto, un factor.

¿Cómo se trabaja con factores en R? La contestación no puede ser completa, porque los factores son un ingrediente clave del lenguaje de R, e intervienen en muchísimas construcciones del sistema. Pero, en cualquier caso, podemos empezar por lo más sencillo.

¿Qué aspecto tiene un vector de tipo factor? En el caso del data.frame experimento que acabamos de crear, al mostrar el campo experimento$Tratamiento se obtiene esta salida:

> experimento$Tratamiento
[1] T1 T1 T4 T1 T2 T5 T5 T2 T5 T2 T3 T1 T5 T3 T3 T4 T1 T4 T1 T5 T4 T3 T4 T4 T1 T4 T1 T2 T3 T5 T5 T4
[33] T2 T4 T4 T5 T1 T3 T3 T1 T3 T2 T2
Levels: T1 T2 T3 T4 T5

Como se ve, R muestra el contenido del vector, junto con sus niveles. Si sólo queremos mostrar los niveles del factor, podemos usar la función levels, cuya salida es esta:

> levels(experimento$Tratamiento)
[1] "T1" "T2" "T3" "T4" "T5"

¿Qué significan esas comillas? Para entenderlo, es preciso conocer la diferencia que hay entre el vector experimento$Tratamiento y un vector de tipo character “con el mismo contenido“. Podemos apreciar la diferencia pidiendo a R que convierta el factor en character mediante la función as.character:

 > as.character(experimento$Tratamiento)
[1] "T1" "T1" "T4" "T1" "T2" "T5" "T5" "T2" "T5" "T2" "T3" "T1" "T5" "T3" "T3" "T4" "T1" "T4" "T1"
[20] "T5" "T4" "T3" "T4" "T4" "T1" "T4" "T1" "T2" "T3" "T5" "T5" "T4" "T2" "T4" "T4" "T5" "T1" "T3"
[39] "T3" "T1" "T3" "T2" "T2"

Las comillas que aparecen en este caso indican precisamente que se trata de un vector de tipo character. Aunque los factores representan variables cualitativas, R los gestiona internamente mediante códigos numéricos. Al mostrarlos, a petición nuestra, el uso de las comillas permite distinguir entre character y factor.
Por cierto, cuando R carga un data.frame mediante read.table, los campos alfanuméricos (que contienen cadenas de texto, como el campo tratamientos de nuestro ejemplo) se convierten, por defecto en factores. Pero puede que, a veces, queramos preservar esos datos como variables de tipo character. Si es eso lo que queremos, basta con usar la opción stringsAsFactors=FALSE en la función read.table.

Los factores permiten modular el comportamiento de algunas funciones de R, de manera que la información que se obtiene tiene en cuenta esa organización en niveles de los valores que estamos analizando. Por ejemplo, en la quinta sesión vimos como usar la función summary para obtener un resumen descriptivo de un vector de datos (media y precentiles, básicamente). Ahora, que tenemos un data.frame con una clasificación por tratamientos, nos gustaría seguramente obtener esa misma información, pero para cada uno de los distintos tratamientos por separado. Lo primero que queremos hacer notar es que no sirve aplicar summary directamente al data.frame:

> summary(experimento)
Respuesta    Tratamiento
Min.   :2.600   T1:10
1st Qu.:4.540   T2: 7
Median :5.040   T3: 8
Mean   :4.791   T4:10
3rd Qu.:5.250   T5: 8
Max.   :6.060

Como se ve, se obtiene un resumen por separado para cada uno de los campos del data.frame: con Respuesta por un lado, y Tratamiento por otro, sin reconocer el vínculo entre ambas. Naturalmente, podemos ir separando cada uno de los grupos de tratamiento, y aplicarles la función summary a los vectores resultantes. Sería algo como esto (mostramos, por ejemplo el segundo grupo de tratamiento):ever

> (tratamiento1=experimento[experimento$Tratamiento=="T3",1])
[1] 5.41 4.87 5.68 5.34 5.71 5.24 5.06 4.52
> summary(tratamiento1)
Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
4.520   5.012   5.290   5.229   5.478   5.710

Pero como se puede apreciar, trabajar así es engorroso. De hecho, la razón para incluir aquí esto es para poder mencionar el símbolo == que aparece ahí. En una próxima entrada, examinaremos con mucho más detalle ese símbolo, y todo lo que tiene que ver con los valores y operadores lógicos, como primer paso hacia la programación en R. Pero, para obtener lo que queremos, hay una manera mucho más natural de proceder, usando la función tapply. Veamos como funciona:

> attach(experimento)
> tapply(Respuesta,Tratamiento, summary)
$T1
Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
4.700   4.955   5.080   5.079   5.228   5.370
$T2
Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
2.600   3.110   3.180   3.214   3.335   3.830
$T3
Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
4.520   5.012   5.290   5.229   5.478   5.710
$T4
Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
4.400   4.605   4.975   5.058   5.400   6.060
$T5
Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
4.560   4.935   5.015   5.039   5.182   5.390

El comando attach(experimento) le dice a R que vamos a usar el dataframe experimento, y nos permite usar las variables Respuesta y Tratamiento  directamente por sus nombres, en lugar de tener que andar escribiendo experimento$Respuesta y experimento$Tratamiento todo el rato.

El resultado, como se ve, es una tabla (de ahí la t inicial de tapply) que contiene, ordenados por filas, los resultados de aplicar la función summary a los valores de Respuesta, agrupadas según los distintos niveles del factor Tratamiento.

Otro ejemplo de las ventajas del lenguaje de factores se obtiene al pensar en la función boxplot. Si aplicamos esa función directamente al data.frame, ejecutando

>
boxplot(experimento)

se obtiene un gráfico que no es, desde luego, el más apropiado, porque de nuevo se tratan respuestas y tratamientos por separado:

BoxPlotDataFrameCompleto

Para obtener lo que queremos basta con escribir:

>
boxplot(Respuesta~Tratamiento)

y el resultado es un gráfico mucho adecuado para comparar los resultados de los tratamientos:BoxPlotDataFramePorFactorLa sintaxis Respuesta~Tratamiento que hemos usado en esta llamada es la que usa R para decir algo así como “la forma en que los valores de Respuesta dependen de los tratamientos“. En general, en R, el símbolo ~ se utiliza para indicar una dependencia o relación entre variables. En este caso, lo que estamos estudiando es si hay alguna relación entre la variable cualitativa (factor) Tratamiento, y la vaiable cuantitativa Respuesta. En particular, nos interesa saber si la respuesta media al tratamiento es distinta según el grupo de tratamiento que se considere. Se trata de un contraste de hipótesis, en la situación prototípica del método ANOVA unifactorial (¡porque sólo hay un factor, claro!) en el que, si llamamos

\mu_1,\mu_2,\ldots,\mu_k

a las respuestas medias de la población tratada con cada uno de los grupos de tratamiento, la hipótesis nula del contraste,sería

H_0=\{\mu_1=\mu_2=\cdots=\mu_k\}

Vamos a ver con más detalle algunas de las herramientas que R pone a nuestra disposición para llevar a cabo este tipo de contrastes ANOVA unifactoriales.

ANOVA unifactorial en R: hipótesis del modelo y tabla ANOVA

La teoría correspondiente se puede consultar en las notas de teoría, que están en esta entrada. Para empezar el trabajo de análisis con R, usaremos la función aov (de analysis of variance, claro) para calcular el modelo ANOVA y poder acceder a varias representaciones gráficas (por ejemplo, de los correspondientes residuos). La sintaxis sería algo como

aov(Respuesta~Tratamiento,datos)

siendo datos el dataframe que recibe como argumento, creado como hemos visto en el apartado anterior.  La ventana de resultados de Rcmdr muestra una respuesta bastante concisa, que no incluye lo que solemos llamar la tabla ANOVA, ni el p-valor del contraste de igualdad de medias:

> (datos.aov=aov(Respuesta~Tratamiento,datos))
Call:
aov(formula = Respuesta ~ Tratamiento, data = datos)
Terms:
Tratamiento Residuals
Sum of Squares    20.968566  5.693996
Deg. of Freedom           4        38

Residual standard error: 0.3870943
Estimated effects may be unbalanced

Pero si se almacenan el resultado de la función aov en una variable, como hemos hecho nosotros con datos.aov, entonces se puede llamar a la función plot así:

plot(datos.aov)

para obtener un análisis gráfico de los residuos del modelo, que proporcionan información adicional, útil para saber si se cumplen los requisitos de aplicación del ANOVA. En particular,  se obtiene una representación de los residuos frente a los valores predichos por el modelo (que son las medias), en la que debemos estar atentos a la aparición de posibles patrones en forma de cuña. En nuestro ejemplo, ese gráfico es ANOVA-ResiduosVsPredichosy no se aprecia en el ningún patrón que indique un fallo evidente en las hipótesis del modelo (lo que se aprecia, de hecho, es que la media de uno de los grupos -niveles- del factor aparece claramente separada de las cuatro restantes; el resto del análisis confirmará esto).

También se obtiene el llamado QQ-plot (de quantile vs quantile), en el que se comparan las distribuciones de los residuos concretos de nuestros datos, con los que se obtendría si la muestra  se ajustara exactamente a una distribución normal. Si las hipótesis del modelo se verifican, en este gráfico esperamos ver aparecer una nube de puntos situados muy aproximadamente a lo largo de una recta. En el ejemplo que nos ocupa, el gráfico QQ-plot es: ANOVA-QQplot-Residuosy como puede verse,  no es exactamente una recta, lo que sugiere que podemos tener algún problema con la hipótesis de normalidad en este ejemplo. ¿Qué hacemos, entonces? Bueno, por un lado, debemos tener presente que las muestras son pequeñas, así que, al igual que sucedía con los diagramas de cajas, el QQ-plot no es demasiado fiable como test de normalidad. En particular no sirve para garantizar que las poblaciones no son normales. Por otro lado, el contraste ANOVA es relativamente robusto frente a las desviaciones pequeñas de la normalidad, siempre y cuando las demás hipótesis del modelo se mantengan, y en particular cuando las varianzas de las poblaciones son homogéneas y todas las muestras son del mismo tamaño.  En este ejemplo, esa segunda condición no se cumple, lo cual nos obliga a ser especialmente cautos. Vamos a empezar por dedicar nuestra atención a la homogeneidad de las varianzas.

Un primer paso puede ser un diagrama de varianzas frente a medias. En este ejemplo se obtiene: ANOVA-VarianzasVsMediasLas varianzas de los grupos parecen ser bastante distintas, lo que podría preocuparnos, pero en este caso lo que más nos interesa es destacar que no hay relación aparente entre el tamaño de la medias y el tamaño de las varianzas. Esta información se une a la del diagrama de residuos frente a valores predichos. En este punto, y para obtener una representación alternativa de la dispersión, podemos volver a los boxplots paralelos por niveles del factor, que no contradicen a los otros gráficos. Como se ve, usamos toda la información disponible para ver si se detecta algún problema evidente, aunque en casos como este, con muestras tan pequeñas, es difícil ser concluyentes.

Antes de seguir adelante, una observación: puede pensarse en utilizar los diagramas de cajas para comprobar si la hipótesis de normalidad  de las poblaciones del modelo se verifican. Pero debe tenerse en cuenta que en el caso de muestras pequeñas, no hay tal cosa como un diagrama de cajas representativo de una distribución normal. El siguiente fragmento de código R ilustra lo que decimos. Se calcula una muestra aleatoria (la media y la desviación típica también se eligen de forma aleatoria) de tamaño 10 de una distribución normal, y se dibuja su diagrama de cajas. Si se ejecuta este código, se verá que la forma de ese diagrama de cajas varía tanto de unas muestras a otras, que realmente no tiene sentido tratar de usar la forma de ese diagrama como comprobación de normalidad.

(media=runif(1,min=-10,max=10))
#(desvTip=sample(seq(from=0, to=abs(media)/2,length.out=100),1))
(desvTip=runif(1,min=0,max=abs(media)/2))
k=10
(muestra=rnorm(k,mean=media,sd=desvTip))
boxplot(muestra)

Y si quieres verlo de otra manera, aquí tienes un applet Java (hecho con GeoGebra, del que ya hablaremos en otra entrada) que ilustra el uso de los boxplots para comprobar la hipótesis de normalidad de una población a partir de sus muestras. Un poco más abajo volveremos sobre el tema del chequeo de la hipótesis de normalidad. El diagrama de boxplots paralelos tampoco parece ser de mucha utilidad para comprobar la homogeneidad de las varianzas. Es fácil usar R para comprobar que aproximadamente el 20% de las muestras de tamaño 10 de la normal estándar tienen desviaciones típicas que se alejan de 1 (la desviación típica de la población) en más de 0.3 unidades.

Volviendo a las hipótesis del modelo ANOVA, ¿cuál es la conclusión? La conclusión, como ya hemos dicho, es que es difícil ser concluyentes con muestras tan pequeñas. ¿Cómo vamos a proceder, entonces? Vamos a usar la tabla ANOVA para calcular el p-valor del contraste. Si ese p-valor es concluyente (en uno u otro sentido), confiaremos en la robustez del método frente a pequeñas desviaciones de las hipótesis. Naturalmente, si el p-valor se sitúa en el umbral de la región de rechazo, deberíamos ser mucho más cuidadosos, mirar con escepticismo las conclusiones y plantearnos la posibilidad de obtener muestras de mayor tamaño.

Para obtener la tabla ANOVA tenemos dos opciones. La primera es usar summary junto con la función aov que ya hemos visto. Recuérdese que hicimos

(datos.aov=aov(Respuesta~Tratamiento,datos))

La tabla ANOVA obtenida mediante summary(datos.aov) es:

> summary(datos.aov)
Df Sum Sq Mean Sq F value   Pr(>F)
Tratamiento  4 20.969   5.242   34.98 2.91e-12 ***
Residuals   38  5.694   0.150
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

La otra opción es usar la sintaxis lm, de linear model, que en R se usa de forma genérica para representar la relación entre variables. De hecho, esta estructura lm juega un papel tan importante en R que, en breve, le dedicaremos una entrada. Para obtener la tabla ANOVA de esta manera hacemos:

anova(lm(Respuesta~Tratamiento,datos))

y el resultado es

> anova(lm(Respuesta~Tratamiento,datos))
Analysis of Variance Table
Response: Respuesta
Df Sum Sq Mean Sq F value    Pr(>F)
Tratamiento  4 20.969  5.2421  34.984 2.906e-12 ***
Residuals   38  5.694  0.1498
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Como se ve, de una u otra manera, resultados son esencialmente iguales en este ejemplo. Y, en particular, el p-valor  (2.906e-12) del  contraste es extremadamente pequeño. Así que podemos rechazar la hipótesis nula, concluir que las medias son distintas, y a la vez estar relativamente seguros de que esta conclusión no depende de pequeñas desviaciones de las hipótesis del  modelo.

En cualquier caso, el método ANOVA no agota, desde luego, el repertorio de recursos para realizar un contraste de igualdad de medios. Existen otros métodos, en especial métodos no paramétricos, que no dependen de las hipótesis de normalidad y homogeneidad de varianzas en la misma medida en que lo hace ANOVA. En R podemos usar uno de estos métodos para confirmar nuestras conclusiones, usando la función oneway.test. En concreto la forma de hacerlo es esta (se muestra la ventana de resultados):

> oneway.test(Respuesta~Tratamiento,datos)
One-way analysis of means (not assuming equal variances)
data:  Respuesta and Tratamiento
F = 36.7064, num df = 4.000, denom df = 17.373, p-value = 2.946e-08

El p-valor (2.946e-08) de este método alternativo confirma la conclusión que habíamos obtenido utilizando ANOVA.  ¿Por qué, entonces, no usamos siempre este método en lugar de ANOVA? La razón principal para hacerlo es que, cuando sus hipótesis se cumplen, el método ANOVA es, en general, más potente que este otro método.

Resultado significativo. ¿Y ahora qué?

Hemos rechazado la hipótesis nula, y por lo tanto, creemos que las medias de los distintos grupos son diferentes entre sí. El siguiente paso natural es, en muchos casos, tratar de comparar por parejas las medias. En nuestro caso eso significa comparar T1 con T2, T1 con T3, etcétera, hasta T4 con T5; un total de diez comparaciones, ya que

\binom{5}{2}=10.

Pero no siempre queremos hacer esto; a veces sólo nos interesa comparar una media concreta con las restantes, en lugar de comparar todas las parejas posibles.  En R se puede hacer esto usando summary combinado con lm de esta forma:

summary(lm(Respuesta~Tratamiento,datos))

y el resultado es:

> summary(lm(Respuesta~Tratamiento,datos))
Call:
lm(formula = Respuesta ~ Tratamiento, data = datos)
Residuals:
Min      1Q  Median      3Q     Max
-0.7087 -0.1590 -0.0180  0.1710  1.0020
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)        5.07900    0.12241  41.492  < 2e-16 ***
Tratamiento[T.T2] -1.86471    0.19076  -9.775  6.4e-12 ***
Tratamiento[T.T3]  0.14975    0.18361   0.816    0.420
Tratamiento[T.T4] -0.02100    0.17311  -0.121    0.904
Tratamiento[T.T5] -0.04025    0.18361  -0.219    0.828
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

¿Cómo se interpreta esto? La parte más importante para nosotros es la tabla que aparece  bajo el epígrafe Coefficients. En esa tabla aparecen los contrastes de igualdad de medias de T1 (el primer nivel del factor; en el orden en que nosotros los hemos ordenado, claro) con el resto de los niveles. El valor que aparece al lado de (Intercept) es simplemente la media del primer grupo (el de T1),  mientras que en las restantes filas de la columna Estimate aparecen las diferencias entre las medias de T1 y las de los restantes niveles. Por ejemplo, la diferencia entre las medias de los tratamientos T1 y T4 es -0.02100. Y el p-valor para el contraste de diferencia de medias entre T1 y T4 está en la última columna, y vale 0.904. Los resultados de esta tabla indican que la media de T1 es significativamente distinta de T2 (tal vez sea buena idea volver a mirar los boxplots paralelos), mientras que no hay diferencias significativas entre la media de T1 y la de los tratamientos T3, T4 y T5.

Naturalmente, esta forma de trabajar privilegia a T1 frente a los demás tratamientos, y no responde a preguntas como ¿hay diferencia significativa entre T3 y T5? Para obtener la información de los diez contrastes por parejas posibles utilizamos la función pairwise.t.test de R, de esta manera (se muestra la ventana de resultados):

> pairwise.t.test(datos$Respuesta,datos$Tratamiento, p.adj="bonferroni")
Pairwise comparisons using t tests with pooled SD
data:  datos$Respuesta and datos$Tratamiento
T1      T2      T3 T4
T2 6.4e-11 -       -  -
T3 1       2.9e-11 -  -
T4 1       8.7e-11 1  -
T5 1       4.3e-10 1  1
P value adjustment method: bonferroni

En la tabla que aparece hay, en efecto, diez p-valores (se evitan las duplicidades mostrando sólo los valores de la parte de la tabla bajo la diagonal). Con esos p-valores podemos establecer qué parejas de medias son relamente distintas. En este ejemplo la conclusión (que ya avanzaban los boxplots) es que la media de T2 es significativamente distinta de todas las demás, pero que no hay diferencias significativas entre T1, T3, T4 y T5 (los p-valores correspondientes son – aproximadamente – 1). Sabemos (ver las notas de teoría), que es preciso ser cuidadosos a la hora de calcular esos p-valores, porque la simple repetición de diez contrastes aumenta mucho la probabilidad de cometer errores de tipo I. Así que R aplica alguno de los métodos habituales (hay varias estrategias, de nuevo remitimos a las notas de teoría y la bibliografía que allí aparece) para ajustar esos p-valores de manera que el error de tipo I se mantenga bajo control. Aquí hemos optado por pedirle a R que utilice la más clásica y conocida de esas estrategias, la de  Bonferroni.

Estos contrastes por parejas, incluso cuando se usan correcciones como la de Bonferroni, dependen de la hipótesis de normalidad de las poblaciones correspondientes a los factores. Para reducir la dependencia de nuestras conclusiones de esa hipótesis,  podemos utilizar algún método no paramétrico, de forma similar a lo que hicimos con la tabla ANOVA.  En este caso, el método más conocido para realizar el contraste es el de Tukey. Sin entrar en detalles, la forma de aplicarlo en R es esta (se muestra la ventana de resultados):

> summary(glht(datos.aov, linfct = mcp(Tratamiento = "Tukey")))
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Tukey Contrasts
Fit: aov(formula = Respuesta ~ Tratamiento, data = datos)
Linear Hypotheses:
Estimate Std. Error t value Pr(>|t|)
T2 - T1 == 0 -1.86471    0.19076  -9.775   <1e-05 ***
T3 - T1 == 0  0.14975    0.18361   0.816    0.924
T4 - T1 == 0 -0.02100    0.17311  -0.121    1.000
T5 - T1 == 0 -0.04025    0.18361  -0.219    0.999
T3 - T2 == 0  2.01446    0.20034  10.055   <1e-05 ***
T4 - T2 == 0  1.84371    0.19076   9.665   <1e-05 ***
T5 - T2 == 0  1.82446    0.20034   9.107   <1e-05 ***
T4 - T3 == 0 -0.17075    0.18361  -0.930    0.883
T5 - T3 == 0 -0.19000    0.19355  -0.982    0.861
T5 - T4 == 0 -0.01925    0.18361  -0.105    1.000
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)

Que es, en esencia, de nuevo una tabla con los diez p-valores, que viene a confirmar lo que habíamos obtenido anteriormente (observar los astericos a la derecha de la última columna, y su significado bajo la tabla; hay cuatro contrastes significativos, todos con T2 y uno de los otros grupos).

Representaciones gráficas. Bargraphs y otras opciones.

Después de un análisis como el que hemos hecho, los gráficos son una de las formas más comunes de presentar los resultados.  Una de las representaciones gráficas más populares, pese a sus inconvenientes son los gráficos de columnas adaptados para incluir líneas de error, como este para el ejemplo que nos ocupa: ANOVA-BargraphParaMedias

Pero hay que tener en cuenta que los gráficos de columnas se concibieron para representar frecuencias, no medias. Y en estos gráficos de columnas adaptados aparece un nivel de referencia (la base de las columnas, que en este caso R ha situado de forma automática en 3), nivel que a menudo no tiene ningún significado real en el contexto del problema que estamos estudiando.  Como puede verse, las barras, cuya altura se corresponde con la media de cada uno de los grupos, se acompañan de una línea o barra de error adicional. Un problema añadido de estos diagramas es que, en muchos casos, no está claro si esas líneas de error representan intervalos de confianza (y cuál es en ese caso el nivel de confianza), o errores estándar, o alguna otra medida de dispersión.  Hay que subrayar que debe dejarse siempre claro, en el texto que acompañe a la figura, cuál es la medida de dispersión que se ha utilizado. Una alternativa es eliminar las columnas, y representar siempre en paralelo intervalos de confianza para la media, como en este gráficoque corresponde a nuestro ejemplo (la idea procede del capítulo 7 de  Introductory Statistics with R, 2nd Ed, de P.Daalgard):ANOVA-GraficoIntervalosConfianzaSe muestra la nube de puntos correspondiente a cada muestra, y barras que corresponden con intervalos de confianza para cada una de esas medias. En este caso la única información necesaria para la correcta interpretación del gráfico es el nivel de confianza de esos intervalos.

Formato del fichero de datos y conclusión

Durante todo el trabajo anterior hemos usado un fichero de datos, llamado experimento-01.csv,.en el que los datos de la respuesta por grupos se presentan en dos columnas. La primera columna contiene el valor numérico de la respuesta, y la segunda el nivel (grupo) del factor Tratamiento.  En otras ocasiones esa misma información se puede presentar de otra manera, con cada uno de los grupos/niveles ocupando su propia columna, como en este fichero experimento-01a.csv, que contiene los mismos datos del ejemplo que hemos usado. El capítulo correspondiente de las notas de teoría contiene un fichero de código (y esa será la versión que en el futuro trataré de mantener actualizada) para el análisis ANOVA unifactorial, que permite realizar todos los pasos del análisis que hemos descrito en esta entrada. Este fichero trata de reconocer automáticamente el formato de entrada del fichero de datos (en las dos situaciones más comunes, las que hemos descrito),  para realizar el análisis correcto. Es conveniente comprobar, en cualquier caso, que los datos se han reconocido correctamente y que no aparece ningún mensaje de error en la lectura. Ese reconocimiento automático utiliza estructuras de programación de R. La existencia de esas estructuras de programación es una de las razones que hacen de R una herramienta tan versatil y potente, así que volveremos sobre ellas en próximas entradas.

Gracias por la atención.

14ª sesión con Rcmdr: Contraste de hipótesis y muestras precocinadas.

Muestras precocinadas, al gusto del cliente

No voy a usar esta entrada para cantar las alabanzas de R como plataforma de análisis de datos. Eso, a estas alturas, es un clamor que ha llegado hasta lugares tan inesperados como las páginas del New York Times. Pero sí quiero hablar del uso de R en la enseñanza de la Estadística, porque eso ayudará a entender la primera parte de esta entrada:. R no se diseñó pensando en la enseñanza, y eso se le nota (“pero mucho”, me parece oír decir a algunos de mis alumnos…).

Por ejemplo, cuando los estudiantes aprenden el concepto de contraste de hipótesis, los ejercicios típicos de los libros de texto se parecen a este:

Se ha recibido un envío de latas de conserva, de las que se afirma que el peso medio son 1000gramos. Al examinar una muestra aleatoria de 5 latas se obtuvo un peso medio muestral de 995 gramos, con una cuasivarianza muestral s²= 19.6. Al nivel de confianza 95%, ¿se puede aceptar que el peso medio son 1000 gramos? Obtener también el p-valor de este contraste.

Naturalmente, al tratarse de un contraste sobre la media, con una muestra de tamaño pequeño, tendremos que usar la distribución t de Student. Y, como no podía ser de otra manera, R incluye una función, llamada t.test, diseñada para realizar este tipo de contrastes. Pero si se intenta utilizar la función t.test para este ejercicio, descubrimos que no es demasiado fácil. La ayuda de R para esta función dice que la función se usa escribiendo

t.test(x, ...)

donde x es un vector de datos numéricos, que contiene los datos de la muestra, claro. Pero no hay nada previsto para introducir directamente los estadísticos muestrales (la media \bar X y la cuasivarianza s^2, y el tamaño n de la muestra), y pedirle a R que haga el contraste usando esos valores sin disponer de la muestra.

Hay que suponer que los programadores de R piensan seguramente que el punto de partida siempre son los datos. Y que nadie hace (en el mundo real, fuera de las aulas y los libros de texto) cálculos estadísticos sin partir de los datos.  Es una suposición muy razonable para el mundo real, pero choca con la forma habitual de enseñar la Estadística.

Otros programas ofrecen la posibilidad de hacer directamente este tipo de ejercicios. ¿Es eso bueno? Podríamos enzarzarnos en una discusión pedagógica, pero lo cierto es que no me interesa demasiado. Una de las principales razones (aparte, naturalmente, del hecho de que es gratuito) por las que usamos R en las clases es porque creemos que es una inversión rentable para nuestros alumnos, porque van a poder seguir usándolo en el futuro, en aplicaciones profesionales. Y, a su vez, una de las razones que explican su éxito en entornos profesionales es su enorme flexibilidad. Si hay algo que echas en falta en R, siempre puedes programarlo tu mismo.  Y eso es lo que ha hecho alguien por nosotros en este caso. Para poder hacer los problemas de los libros de texto directamente, vamos a utilizar una función de R que prepara muestras “al gusto del consumidor”. Es decir, nosotros le decimos a R cuáles son  los valores de

\bar X, s^2 \mbox{ y } n

y R nos “cocinará” una muestra

x_1,x_2,\ldots,x_n

con exactamente esos parámetros, para que podamos dárselos como alimento a la función t.test. La función que hace esto por nosotros es mvrnorm. Por ejemplo, estos comandos de R fabrican una muestra de una población aproximadamente normal, cuya longitud, media y cuasivarianza se han elegido de antemano. Y a continuación comprobamos que esos valores son los que queríamos.

# Fijamos los parámetros deseados de la muestra,
n=10
mediaMuestral=8.2
cuasiVarianzaMuestral=1.3
# formamos la muestra,
(muestra=mvrnorm(n, mediaMuestral,cuasiVarianzaMuestral, empirical = TRUE))
# y comprobamos que responde a nuestros deseos.
mean(muestra)
var(muestra)

La parte interesante del resultado es la comprobación, claro:

> mean(muestra)
[1] 8.2
> var(muestra)
[,1]
[1,]  1.3

Con esto estamos preparados para ver como funciona la función t.test de R, y el resto de funciones que realizan contrastes de hipótesis.

Contrastes de hipótesis en R (t.test y sus parientes)

Vamos a ver como usar la función t.test de R para realizar un contraste de hipótesis sobre la media. Usaremos este ejemplo (que apreciará cualquiera que haya estado parado en un paso a nivel canadiense, esperando a que termine de pasar el mercancías de turno):

Una compañía ferroviaria canadiense afirma que sus trenes de mercancías no bloquean los pasos a nivel durante más de 8 minutos, en promedio. Una muestra aleatoria de 10 tiempos de bloqueo dio como resultado estos valores (en minutos):

10.1, 9.5, 6.5, 8.0, 8.8, 12,  7.2, 10.5, 8.2, 9.3

Empezamos por observar que en este caso tenemos todos los valores de la muestra. Llamando μ al tiempo medio de bloqueo, queremos usar estos valores para contrastar la hipótesis nula:

H_0=\{ \mu \leq \mu_0=8 \}

Y naturalmente la hipótesis alternativa es:

H_a=\{ \mu > \mu_0=8 \}

El estadístico del contraste es

\dfrac{\bar X-\mu_0}{\dfrac{s}{\sqrt{n}}}

Hagamos primero los cálculos del contraste “a mano”, sin recurrir a t.test. Vamos a fijar un nivel de significación del 5%, es decir, 0.05. Puesto que se trata de una muestra pequeña (n=10), usaremos la distribución t de Student para el cálculo del p-valor.

datos=c(10.1, 9.5, 6.5, 8.0, 8.8, 12,  7.2, 10.5, 8.2, 9.3)
mu0=8
(xBarra=mean(datos))
(s=sd(datos))
(n=length(datos))
EstadisticoContraste = ( xBarra - mu0 ) / ( s / sqrt(n) )
(pValor=1-pt(EstadisticoContraste,df=n-1))

El resultado es este:

> datos=c(10.1, 9.5, 6.5, 8.0, 8.8, 12,  7.2, 10.5, 8.2, 9.3)
> mu0=8
> (xBarra=mean(datos))
[1] 9.01
> (s=sd(datos))
[1] 1.631938
> (n=length(datos))
[1] 10
> (EstadisticoContraste = ( xBarra - mu0 ) / ( s / sqrt(n) ))
[1] 1.957121
> pValor=1-pt(EstadisticoContraste,df=n-1)
> (pValor=1-pt(EstadisticoContraste,df=n-1))
[1] 0.04101152

De manera que, con un nivel de significación 0.05 (mayor que el p-valor), tenemos evidencia empírica para rechazar la hipótesis nula y concluir que los trenes bloquean el paso a nivel más tiempo del que dice la empresa. Veamos ahora como hacer el contraste usando t.test:

datos=c(10.1, 9.5, 6.5, 8.0, 8.8, 12,  7.2, 10.5, 8.2, 9.3)
mu0=8
t.test(datos,mu=mu0, alternative = c("greater"), conf.level = 0.95)

Como se ve, aparte del vector de datos, le hemos indicado a R el valor de \mu_0, y mediante las opciones  alternative = c(“greater”)conf.level = 0.95 hemos seleccionado un contraste de cola derecha (greater, o mayor que) y el nivel de significación deseado (conf.level, R usa la misma terminología que para los intervalos de confianza) .

El resultado es este:

> t.test(datos,mu=mu0, alternative = c("greater"), conf.level = 0.95)
One Sample t-test
data:  datos
t = 1.9571, df = 9, p-value = 0.04101
alternative hypothesis: true mean is greater than 8
95 percent confidence interval:
8.063996      Inf
sample estimates:
mean of x
9.01

Y como se ve, la respuesta contiene tanto el valor del estadístico de contraste ( t = 1.9571), como el p-valor = 0.04101. Además, para que la interpretación del resultado sea más fácil, y para que podamos comprobar que estamos haciendo lo que deseamos, R describe la hipótesis alternativa del contraste. Como subproducto se obtiene lo que R llama un intervalo de confianza, pero que en el caso de un contraste unilateral como este, no tiene la forma que seguramente espera un estudiante de un curso de introducción a la Estadística. Para obtener la forma del intervalo de confianza que aparece en los libros introductorios, basta con ejecutar t.test para hacer un contraste bilateral (que es la opción por defecto):

> t.test(datos,mu=mu0, alternative = c("two.sided"), conf.level = 0.95)
One Sample t-test
data:  datos
t = 1.9571, df = 9, p-value = 0.08202
alternative hypothesis: true mean is not equal to 8
95 percent confidence interval:
7.842582 10.177418
sample estimates:
mean of x
9.01

El intervalo de confianza que obtenemos seguramente se parece más a lo que se esperaba.

Si has estudiado los  contrastes de hipótesis en el orden habitual, seguramente te estarás preguntado ¿y el contraste para la media con la Z, la normal estándar?  Lo cierto es que esos contrastes Z son casi exclusivamente “ejemplos de libro de texto”, que no se usan en las aplicaciones reales. Y no se incluyen en R por defecto (¿hemos dicho ya que R no se diseñó pensando en la enseñanza?). Pero eso no significa que no estén disponibles. Basta con cargar una librería, cuyo revelador nombre es TeachingDemos, y listos: ya tenemos disponible una función z-test, con la que podemos hacer cosas como esta: suponemos que los datos del ejemplo anterior proceden de una población normal con desviación típica σ conocida, e igual a 1.6, y realizamos el mismo contraste para la media que hicimos antes:

datos=c(10.1, 9.5, 6.5, 8.0, 8.8, 12,  7.2, 10.5, 8.2, 9.3)
z.test(datos, mu =mu0, sd=1.6, alternative = c("greater"), conf.level = 0.95)

El resultado es:

> datos=c(10.1, 9.5, 6.5, 8.0, 8.8, 12,  7.2, 10.5, 8.2, 9.3)
> library("TeachingDemos")
> z.test(datos, mu =mu0, sd=1.6, alternative = c("greater"), conf.level = 0.95)
One Sample z-test
data:  datos
z = 1.9962, n = 10.000, Std. Dev. = 1.600, Std. Dev. of the sample mean = 0.506, p-value = 0.02296
alternative hypothesis: true mean is greater than 8
95 percent confidence interval:
8.177763      Inf
sample estimates:
mean of datos
9.01

Es un ejercicio saludable repetir estas cuentas “a mano”, usando pnorm, para comprobar que somos capaces de obtener el p-valor de este contraste de esa forma, sin recurrir a z.test.

Volviendo a t.test, esa misma función puede usarse para contrastes de medias entre dos poblaciones, simplemente introduciendo sendos vectores con las muestras correspondientes. Por ejemplo, dado este ejercicio:

Una empresa quiere introducir un nuevo método de trabajo, y quiere saber si ese método influye sobre el tiempo que se emplea en llevar a cabo una determinada tarea. Se han seleccionado aleatoriamente dos grupos, cada uno de cinco empleados, a los que se les ha pedido que lleven a cabo la tarea por uno de los dos métodos. La tabla de tiempos empleados es esta:

\begin{tabular}{|l|c|c|c|c|c|} \hline\mbox{M\'etodo antiguo}&20&18&18&19&17\\ \hline \mbox{M\'etodo nuevo}&18&15&17&16&17\\ \hline \end{tabular}

¿Se puede afirmar, con un nivel de significación del 5%, que el  tiempo que se emplea en esa tarea es el mismo para ambos métodos?

El resultado de llevar a cabo este contraste en R, utilizando la función t.test, es el siguiente:

> muestra1=c(20,18,18,19,17)
> muestra2=c(18,15,17,16,17)
> t.test(muestra1,muestra2,alternative=c("two.sided"),conf.level=0.95)
  Welch Two Sample t-test
data:  muestra1 and muestra2
t = 2.4962, df = 8, p-value = 0.03716
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.1371168 3.4628832
sample estimates:
mean of x mean of y
18.4      16.6

Algunos comentarios sobre este ejemplo: en este caso el contraste es bilateral, y por eso la opción que usamos es alternative=c(“two.sided”). Si fuera unilateral, usaríamos, como en el ejemplo previo, o bien alternative=c(“greater”), o bien alternative=c(“less”), y en ambos casos las muestras se usan en el orden en que se escriben. Es decir, que si usáramos muestra1,muestra2,alternative=c(“greater”), estaríamos utilizando como hipótesis alternativa que la media de la población 1 es mayor (greater) que la media de la población 2.
Además, el lector atento habrá observado que no le hemos indicado a R si estamos suponiendo que las varianzas de las poblaciones se suponen iguales o distintas (son muestras de tamaño muy pequeño, así que es necesario hacer esa distinción). La opción var.equal=FALSE (o bien, TRUE, claro) nos permite elegir explícitamente el comportamiento de R. Por defecto R utiliza, lógicamente, el valor FALSE. En particular, se usa una corrección de tipo Welch para los grados de libertad. Y digo una porque esta mañana he estado echando un vistazo en mi biblioteca, y he encontrado al menos tres expresiones distintas (todas etiquetadas como “de Welch”) para el cálculo de los grados de libertad en esta situación. R usa una de ellas, y eso significa ¡ATENCIÓN! que el intervalo de confianza que calcula R puede no coincidir exactamente con el que se obtiene usando las fórmulas de tu libro o de tus apuntes de clase. SI quieres asegurarte esa coincidencia, calcula el intervalo a mano (usando qt, y la fórmula de Welch que corresponda).

Los contrastes de proporciones también pueden hacerse en R, utilizando en este caso la función prop.test. Por ejemplo, si tenemos dos muestras independientes, de tamaños n1=100 y n2=250, en las que el número de éxitos es  respectivamente, e1=33 y e2=80, la forma de usar prop.test para calcular un contraste de igualdad de proporciones se ilustra en este ejemplo:

> prop.test(c(33,80),c(100,250),alternative = c("two.sided"),conf.level = 0.95, correct = FALSE)
2-sample test for equality of proportions without continuity
correction
data:  c(33, 80) out of c(100, 250)
X-squared = 0.0327, df = 1, p-value = 0.8566
alternative hypothesis: two.sided
95 percent confidence interval:
-0.09879831  0.11879831
sample estimates:
prop 1 prop 2
0.33   0.32

Quizá lo mas peculiar es la forma en la que se introducen los datos muestrales, de manera que el primer vector contiene los tamaños de ambas muestras, y el segundo los recuentos de éxitos. El resto de las opciones son fáciles de comprender, con la excepción de correct = FALSE.  Con esta opción desactivamos una corrección de continuidad que aplica R, para que los resultados sean (menos precisos pero) más parecidos a los que se aprende a calcular en un curso de introducción a la Estadística, usando la distribución normal estándar Z.

Obsérvese que en este caso de proporciones, la sintaxis de R hace innecesario el “cocinado” de muestras que hemos explicado en la primera parte de la entrada.

También podemos llevar a cabo contrastes sobre la varianza (y, por consiguiente, la desviación típica). Para una sola población, disponemos , en la librería TeachingDemos (de la que ya hemos hablado), de la función sigma.test. En este ejemplo vemos como usarla para un contraste unilateral sobre la desviación típica de una población (se muestra la salida en la ventana de resultados):

> muestra=c(9.6,11.4,8.5,10.8,10.5,10.7,10.1,9.8,8.9,11.6,11.4,12.9)
> library("TeachingDemos")
> sigma.test(muestra,sigma=2,alternative=c("less"),conf.level=0.95)
One sample Chi-squared test for variance
data:  muestra
X-squared = 4.1842, df = 11, p-value = 0.03578
alternative hypothesis: true variance is less than 4
95 percent confidence interval:
0.000000 3.658437
sample estimates:
var of muestra
1.521515

Para los aprendices de la Estadística es bueno leer la ayuda de la función sigma.test para entender porque está función aparece “relegada” a la librería TeachingDemos.

Si tenemos que hacer un contraste de igualdad de varianzas entre dos poblaciones independientes, podemos hacerlo empleando la función var.test (no es necesario cargar librerías adicionales para usarla. ni siquiera TeachingDemos). Para no extendernos más, remitimos al lector  a la ayuda de esa función, y ponemos aquí punto final a esta entrada.

Gracias por la atención.

13ª sesión con Rcmdr: un momento de respiro, y algo de ayuda.

Esta semana, por diversas razones,  ha sido complicada, y además la mayoría de mis alumnos se examinan la semana que viene. Así que seguramente todos, ellos y yo, agradeceremos una entrada del blog más relajada que las anteriores.  Además, estoy escribiendo una entrada más larga, con los primeros pasos en R como lenguaje de programación, y todavía no está acabada. Así que, mientras ese trabajo llega a puerto, haremos una excursión breve a algunos temas sueltos que han ido quedando por el camino. El denominador común es mejorar nuestra experiencia como usuarios de R.

Lo primero es hacer justicia con quienes nos abrieron el camino

Aunque este blog aspira a proporcionar una lectura autosuficiente sobre R, lo cierto es que hay muchas fuentes de información que yo usé, en su momento,  para aprender a usar R y a las que sigo acudiendo en busca de lo mucho, muchísimo, que siempre queda por delante. Y os animo a que las exploréis. Hay gente haciendo cosas increíbles, que nos hacen soñar con llegar a ser algún día la décima parte de buenos en lo que hacemos.

La primera de esas fuentes no hay que ir a buscarla muy lejos. En cualquier instalación de R hay disponible un documento, al que todos nos referimos como R-Intro.  En Windows, por ejemplo, para localizarlo, vamos a la ventana RGui (donde aparecen los gráficos, no la de Rcmdr), y en el menú Ayuda seleccionamos Manuales (en pdf) y luego An Introduction to R, como en esta figura.

Este documento, en inglés (en este enlace hay una traducción al castellano de una versión algo más antigua),  es relativamente fácil de leer y constituye un curso de choque a algunas de las posibilidades del programa. Es el primer paso para empezar a convertirse en un usuario avanzado de R, aunque -pese a su título- a veces da la sensación de que está más pensado como una guía de referencia que como una ayuda al aprendizaje.

Aparte de este, y otros documentos parecidos, hay varios (muchos) blogs disponibles con un espíritu similar a este.  Sin embargo, no conozco ninguno, en español, que se plantee como objetivo enseñar a usar R a un completo novato, como acompañamiento para un primer curso en Estadística. Además, incluso en los blogs en inglés, muchas veces se presentan muchos detalles técnicos de R antes de empezar con las aplicaciones estadísticas. Y por esas razones, entre otras, existe este blog.

Volviendo a los blogs en inglés, hay dos que he consultado bastantes veces, y que tienen nombres casi idénticos :

  • R tutorial. Contiene a su vez varios temas, de los cuales los más relacionados con este blog son por un lado Elementary Satistics with R y por otro R Introduction. Contiene numerosos ejemplos de problemas estadísticos resueltos con R.
  • R tutorials. Probablemente el blog más parecido a este en sus intenciones. Pero no usa Rcmdr, sino la interfaz más espartana de R (y yo no me atrevo a infligir semejante castigo a mis ya suficientemente sufridos lectores). En cualquier caso, contiene un montón de información útil y bastante bien organizada.

Pero el universo de los blogs sobre R es mucho más amplio. De hecho, existe una página (R bloggers) que se dedica a recopilar información sobre blogs que hablan de R (y tienen más de 400 registrados, mayoritariamente en inglés, pero también en otros idiomas…). Así que puedes buscar el que más te guste.

Ayuda directa sobre comandos desde el propio R

Por otra parte, el propio R, como cualquier programa mínimamente decente, incluye documentación interna sobre sus propios comandos. Por ejemplo, ya conocemos el comando sample,  que hemos usado en la novena sesión para fabricar muestras. Para acceder a la ayuda interna sobre el comando, tecleamos en R su nombre precedido de un paréntesis, así:

?sample

Al ejecutar este comando se abrirá nuestro navegador de internet (pero no hace falta estar conectados, la ayuda reside en nuestro ordenador), y nos mostrará  la documentación del comando sample (en inglés):De forma equivalente, puedes usar el comando help(sample) para obtener el mismo documento.

Toda la documentación de R está disponible también en internet, con copias en múltiples servidores. Así que muchas veces, a pesar de todos los sistemas de ayuda que hemos mencionado,  lo más rápido es usar un buscador de internet, teclear algo como “Binomial distribution with R“, y hojear los resultados.

Asignaciones con flechas

Es algo que he querido mencionar en el  blog desde hace tiempo, pero que nunca parecía encontrar su sitio. Hasta hoy, claro. SI vais a consultar el código r escrito por otras personas, descubriréis que, donde yo pondría:

datos = c(2,5,12,3,4)

mucha otra gente escribe:

datos <- c(2,5,12,3,4)

Lo que estamos haciendo, en cualquier caso, es una asignación, en este caso guardando un vector de números en la variable datos. Y la notación tradicional en R para las asignaciones utiliza esa especie de flecha <-  en lugar del símbolo igual. Es una buena idea, en principio, porque el símbolo = ya se utiliza en matemáticas para indicar que dos cosas son iguales, y no con este sentido de “guarda esto en la variable…” Pero en la práctica hay tres razones que me llevan a usar = en lugar de la flecha.

  • la primera, es que el símbolo = se utiliza en muchos lenguajes de programación para las asignaciones (no en todos, pero sí en muchos de los más usados).
  • la segunda es que el símbolo elegido, la flecha <- , no me parece demasiado afortunado. En otros lenguajes se usan alternativas como := (dos puntos seguidos de un igual), que me parecen más logradas.
  • y finalmente, ¡da igual! Puedes cambiar las flechas por iguales o viceversa, y R no protestará (al menos, ninguna versión moderna de R).

Gracias por la atención.

Duodécima sesión con Rcmdr: funciones de densidad para variables aleatorias continuas.

Definición de funciones en R

En las dos últimas sesiones hemos visto sendos ejemplos de la forma en la que R gestiona una distribución discreta (la binomial, en la 10ª sesión), y una distribución continua (la normal, en la 11ª sesión). Se trata en ambos casos, de distribuciones muy importantes, con nombre propio, por así decirlo, y a las que R conoce bien. Por eso existen funciones específicas como dbinom, pnorm, etc., para trabajar con estas distribuciones en concreto.

Por otra parte, en la novena sesión,  hemos visto como utilizar R para tratar con variables aleatorias discretas genéricas (que no tienen nombre propio, como la binomial), definidas mediante una tabla de densidad de probabilidad.  Y la pregunta es evidente: ¿qué ocurre con las variables aleatorias genéricas de tipo continuo? ¿Cómo las podemos manejar desde R?

Como sabemos, una variable aleatoria continua X se caracteriza mediante una función de densidad, f(x), que permite asignar un valor de probabilidad a cada intervalo de valores de la variable X. Para un intervalo de la forma (a,b), la probabilidad de que X tome valores en ese intervalo se calcula mediante:

\int_a^b f(x)dx

Para poder realizar este tipo de cálculos en R, vamos a aprender dos cosas:

  1. En primer lugar, aprenderemos a definir funciones en R.
  2. A continuación, veremos como calcular sus integrales en un intervalo de valores.

Empecemos con la definición de una función en R. Queremos subrayar desde el principio que, aunque aquí las vamos a usar como funciones de densidad, las funciones se pueden usar para muchas otras cosas en R (veremos ejemplos de esto en próximas entradas).  La definición de funciones (de una o más variables) es una de las características que hacen de R un programa tan flexible y potente.

Supongamos que la función de densidad de nuestra variable aleatoria X es

f(x)=6x(1-x)

en el intervalo (0,1) (más abajo veremos que es realmente una función de densidad). Para definir esta función en R usamos este comando:

f = function(x){ 6*x*(1-x) }

Es decir, f es el nombre que nosotros hemos elegido para la función. A continuación viene el comando function de R, seguido del nombre de la variable (o las variables, puede haber varias separadas por comas). Y finalmente, lo que llamaremos el cuerpo de la función, que es la parte en la que explicamos como calcular el valor de la función a partir de la variable o variables. Aunque no es imprescindible en todos los casos, es bueno acostumbrarse a escribir el cuerpo de la función entre llaves, porque eso hace más fácil la lectura de los ficheros de código que escribimos.

Una vez definida la función, podemos usarla como las otras funciones de R. Por ejemplo, para saber el valor de f(x) en x=1/2  hacemos

f(1/2)

y en este ejemplo se obtiene 1.5 como respuesta. Insistimos, estas funciones que definimos son como las funciones del propio R. Podemos por ejemplo, aplicar la función f(x) a un vector:

datos=c(1,3,-1,2,5)
f(datos)

y se obtiene como respuesta:

> datos=c(1,3,-1,2,5)
> f(datos)
[1]    0  -36  -12  -12 -120

Para dibujar la gráfica de la función de densidad (o de cualquier función de una variable, en realidad) podemos usar el comando curve que vimos en la 11ª sesión. En concreto, para la f(x) que estamos usando de ejemplo, en el intervalo (0,1), hacemos:

curve(f,from=0,to=1,col="red",ylim=range(0, 1.5),lwd=6)

(ver la anterior sesión para la explicación de las opciones) y se obtiene:

Integrales. Cálculo de la media y varianza de una variable aleatoria continua.

Ahora que ya sabemos definir la función de densidad, vamos a pasar a la segunda parte del plan. Aprenderemos a calcular (siempre aproximadamente) la integral de la función en un intervalo. Y como aplicación de esto, vamos a calcular la media y desviación típica de la variable aleatoria definida por esa función de densidad.

La integración se lleva a cabo mediante el comando integrate de R Debemos indicar el nombre de la función, y los extremos del intervalo de integración, mediante las opciones lower y upper. Por ejemplo, la integral de nuestra función de ejemplo f(x) en el intervalo (1/3,2/3) se obtiene mediante:

integrate( f, lower=1/3, upper=2/3)

La respuesta que se obtiene es:

> integrate( f, lower=1/3, upper=2/3 )
0.4814815 with absolute error < 5.3e-15

Como se ve, R no se limita a calcular el valor aproximado de la integral sino que además nos da información sobre la calidad de esa aproximación. Eso está muy bien, pero si queremos utilizar el resultado en otra operación, puede ser un engorro. afortunadamente el remedio es muy  fácil. Sólo hay que añadir una pequeña modificación al final de la llamada al comando integrate:

integrate( f, lower=1/3, upper=2/3 )$value

y ahora se obtiene como salida un número (el valor aproximado de la integral),

> integrate( f, lower=1/3, upper=2/3 )$value
[1] 0.4814815

que podemos usar en otras operaciones. Como segundo ejemplo, vamos a comprobar algo que hemos dejado pendiente. Veamos que la función f(x) que estamos usando es, realmente, una función de densidad en el intervalo (0,1), comprobando que su integral en ese intervalo es 1. Ejecutamos este comando:

integrate(f, lower=0, upper=1)$value

y la respuesta es, en efecto, 1.

¿Y para calcular la media \mu de la variable X cuya densidad es f(x)? Recordemos que, si la densidad de X en (a,b), es una función f(x), entonces la media de X es:

\mu_X=\int_a^b x f(x) dx

Para aplicar esto a la función f(x), como primer paso creamos una función auxiliar, que llamaremos  aux_f (elegimos ese nombre para recordar su papel; podría ser cualquier otro), cuyo valor es el de f(x) multiplicado por x.

aux_f = function(x){ x * f(x) }

Es interesante observar que la definición de esta función incluye una llamada a la propia función f(x), en lugar de copiar directamente su definición. De esa manera, esta función auxiliar sigue siendo válida incluso aunque la función f(x) original cambie,

Y ahora, para calcular la media \mu de f basta con calcular la integral de esta función auxiliar entre 0 y 1.

mu=integrate(aux_f, lower=0, upper=1)$value
mu

Se obtiene \mu=2 como valor de la media.

El siguiente paso lógico es obtener la varianza \sigma (y desviación típica) de X.  Ahora, recordando que la fórmula para la varianza es

\sigma^2_X=\int_a^b (x-\mu)^2 f(x) dx

ya debería estar claro que el primer paso es definir una nueva función auxiliar:

aux2_f = function(x){ (x-mu)^2 * f(x) }

e integrarla entre 0 y 1.

varianza=integrate(aux2_f, lower=0, upper=1)$value
varianza

El valor que se obtiene de la varianza es 0.05=\frac{1}{20}. La desviación típica es simplemente la raíz cuadrada de este número (calculada, por ejemplo, con la función sqrt de R).

Función de distribución (acumulada)

La función de distribución (acumulada) de una variable aleatoria continua X, definida en el intervalo (a,b) por la función de densidad f(x) es

F(x)=P(X\leq x)=\int_a^x f(t)dt

Para aclararlo un poco: si, por ejemplo, (a,b)=(1,5), y queremos calcular F(3), entonces tenemos que integrar la densidad f(x) desde 1 (el principio del intervalo) hasta 3 (el valor en el que calculamos F).

Esta función de distribución F se puede definir en R  de forma sencilla, recurriendo al comando integrate. En lugar de F, que ya se usa para el valor lógico FALSE en R, vamos a llamarla Df en R. Eso sí, primero debemos almacenar en una variable el comienzo del intervalo donde está definida X:

a=0 #extremo inferior del intervalo
Df = function(x){ integrate(f, lower=a, upper=x)$value }

Obsérvese el uso de $value, para obtener el valor de la integral. Ahora podemos pedirle a R que calcule cualquier valor de la función de distribución (acumulada) Df, Por ejemplo, evidentemente Df1)=1 en este ejemplo; pero también se obtiene Df(1/3)=0.2592593, lo cual significa que:

P(X<1/3)=0.2592593

Recordemos que en el caso de la distribución normal, la función de distribución acumulada  (con el mismo sentido que estamos usando aquí) se obtiene mediante pnorm. Y en general, para cualquier distribución “con nombre propio”, el prefijo p indica que calculamos la función de distribución. Lo que hemos aprendido a hacer es calcular esto para nuestras propias variables aleatorias continuas.

Muestreo de variables aleatorias continuas

Cuando, en la novena sesión,  vimos como utilizar R para tratar con variables aleatorias discretas genéricas, terminamos aquella sesión explicando como obtener muestras aleatorias de ese tipo de variables usando la función sample de R.

Desafortunadamente, sample no puede usarse en el caso continuo. Hay, no obstante una forma de obtener una muestra de una variable aleatoria continua dada mediante su función de densidad. El problema es que, para entender porque este método funciona, es necesario conocer más Estadística y más sobre el propio R de lo que hasta ahora hemos visto en el blog. Así que nos vamos a limitar a dar un código que funciona, y a dejar las explicaciones para más adelante.  Una advertencia, eso sí. El cálculo de muestras grandes, con funciones de densidad complicadas, puede llevar bastante tiempo (minutos, horas; en algunos caso puede resultar inviable).

El fichero R que cierra esta entrada contiene, resumidas, la mayor parte de las instrucciones que hemos comentado en esta entrada, incluido el muestreo.

Gracias por la atención.

##############################################################
#
# Estadística, Grado en Biología y Biologia Sanitaria
# Universidad de Alcala
#
# Fichero de instrucciones R para trabajar
# con variables aleatorias continuas genéricas.
##############################################################

#Definción de la función de densidad
f = function(x){ 6*x*(1-x) }

#Intervalo (a,b) de definición de la variable
a=0
b=1

#Comprobación de coherencia
integrate( f, lower=a, upper=b)$value
#El resultado debe ser 1, de lo contrario tenemos un problema.

#Cálculo de la media
aux_f = function(x){ x * f(x) }
mu=integrate(aux_f, lower=a, upper=b)$value
mu

#Cálculo de la varianza
aux2_f = function(x){ (x-mu)^2 * f(x) }
varianza=integrate(aux2_f, lower=a, upper=b)$value
varianza

# y desviación típica
desvTipica=sqrt(varianza)
desvTipica

#Definimos la función de distribución acumulada.
Df = function(x){ integrate(f, lower=a, upper=x)$value }

#Construimos una muestra aleatoria de la variable.

#Número de valores en la muestra
nValores = 1000;
puntosUniformes = runif(nValores);

#El código siguiente genera un vector muestra_f que contiene la muestra aleatoria.
sampling_f = function(x, u) {
return(Df(x) - u);
}

muestra_f = c();
for (i in 1:nValores) {
# find the root within (a,b)
r = uniroot(sampling_f, c(a,b), tol = 0.0001, u = puntosUniformes[i])$root;
muestra_f = c(muestra_f, r);
}
muestra_f

#Puede ser interesante terminar comparando el diagrama de frecuencias de la muestra con la función de densidad original.
h = hist(muestra_f, plot=F)
ylim = range(0, h$density, 0)
hist(muestra_f, freq=F, density=20,ylim=ylim)
curve(f,add=T,col="blue",lwd=4)