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.

Anuncios

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

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

  2. Pingback: Tablas de contingencia y contraste de independencia (test chi cuadrado) | PostData

Responder

Introduce tus datos o haz clic en un icono para iniciar sesión:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión /  Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión /  Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión /  Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión /  Cambiar )

Conectando a %s