PostData

Curso de Introducci´ on a la Estad´ıstica

Tutorial 07: Contraste de hip´ otesis. Factores en R. Atenci´ on: Este documento pdf lleva adjuntos algunos de los ficheros de datos necesarios. Y est´a pensado para trabajar con ´el directamente en tu ordenador. Al usarlo en la pantalla, si es necesario, puedes aumentar alguna de las figuras para ver los detalles. Antes de imprimirlo, piensa si es necesario. Los ´ arboles y nosotros te lo agradeceremos. Fecha: 2 de diciembre de 2013. Si este fichero tiene m´as de un a˜ no, puede resultar obsoleto. Busca si existe una versi´ on m´ as reciente.

´Indice 1. Contraste de hip´ otesis (una poblaci´ on).

1

1.1. Contrastes para µ y σ en pob. normales, usando R . . . . . . . . . . . . . . . . . .

1

1.2. La funci´ on t.test de R (y sus parientes) . . . . . . . . . . . . . . . . . . . . . . .

5

1.3. Contrastes de hip´ otesis para µ y σ (una poblaci´on normal) usando otros programas

8

2. Intervalos de confianza y contrastes de hip´ otesis para el par´ ametro p, proporci´ on binomial 10 2.1. Con R . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

10

2.2. Con otros programas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

11

3. Variables cualitativas (factores) en R

12

3.1. Factores . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

13

3.2. Usando los factores para explorar los datos . . . . . . . . . . . . . . . . . . . . . .

15

3.3. La funci´ on cut . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

17

1.

Contraste de hip´ otesis (una poblaci´ on).

Vamos a empezar este tutorial aprendiendo a utilizar R (y, en menor medida, otros programas) para llevar a cabo contrastes de hip´otesis como los que se discuten en el Cap´ıtulo 7 del curso. Nuestro objetivo, en todos los casos, es calcular el p-valor del contraste, y establecer los l´ımites de la regi´ on de rechazo de la hip´ otesis nula H0 .

1.1.

Contrastes para µ y σ en pob. normales, usando R

Empecemos por los contrastes para la media. Y vamos a suponer que el tama˜ no de la muestra es grande. La terminolog´ıa y notaci´ on que usaremos est´a en las Secciones 7.2 y del curso. El esquema del contraste es este: 1. Fijamos µ0 , y establecemos la hip´otesis nula y la alternativa. La forma de las hip´otesis depende de que estemos en un contraste bilateral o unilateral; y en este segundo caso, depende de cu´ al sea el lado.

1

2. Con los datos de la muestra, calculamos el estad´ıstico adecuado. Este es el paso clave. Puede ser u ´til consultar las tablas del Ap´endice B del curso en este paso. 3. Usando pnorm calculamos el p-valor, y usando qnorm calculamos los l´ımites de la regi´on de rechazo (aqu´ı interviene el nivel de significaci´on del contraste). Vamos a escribir un programa en R para resolver los ejercicios, t´ıpicos de libro de texto, que son habituales cuando se estudian los contrastes de hip´otesis. La parte no mec´anica de este tipo de ejercicios, la que no podemos programar en R, es aquella en la que analizamos el problema y decidimos el tipo de contraste que vamos a hacer. Casi todo lo dem´as es programable. Las decisiones que hay que tomar durante el proceso que hemos esbozado se pueden implementar a trav´es de estructuras condicionales de tipo if-else, como las que hemos visto en la Secci´on 9 del Tutorial05. Antes de seguir adelante, recomendamos encarecidamente al lector que repase, en el Cap´ıtulo 7 del curso los Ejemplos 7.2.1, p´ ag. 190, (y sus continuaciones, en la Secci´on 7.2), y muy especialmente el Ejemplo 7.4.1 (p´ ag. 202). Y para hacerlo, le proponemos un par de ejercicios. Ejercicios: 1. Calcula con R los p-valores y regiones de rechazo (al 95 %) que aparecen en esos ejemplos. Solo es necesario utilizar pnorm y qnorm, que ya conocemos. Te recuerdo los datos. En el Ejemplo 7.2.1, la hip´otesis nula es H0 = {µ ≤ µ0 }, con µ0 = 2.5. Los datos de la muestra son n = 100,

¯ = 2.65, X

s = 0.5

En el Ejemplo 7.4.1 todos los datos son iguales, salvo la media muestral, que ahora es ¯ = 2.35 X 2. Calcula tambi´en el p-valor y regi´on de rechazo (95 %) de un contraste (bilateral) para la media en una poblaci´ on normal, con hip´otesis nula: H0 = {µ = µ0 }, siendo µ0 = 3.1. La muestra que se usa para el contraste tiene estos datos: n = 250,

¯ = 2.6, X

s = 2.7

¡No sigas, si no has hecho este ejercicio!

2

La experiencia que proporciona ese ejercicio debe ayudar a entender algunas de las decisiones que hay que tomar en un contraste, y que son las que el programa que vamos a describir toma por nosotros. Nosotros (el usuario humano de R) tenemos los datos del problema, y hemos tomado ya la decisi´ on del tipo de contraste que vamos a realizar. ¿Cuales son, entones, esas decisiones que tiene que tomar nuestro programa? A partir del tipo de contraste, el programa tiene que decidir si utiliza la cola derecha o izquierda, en el caso de contrastes unilaterales. Si es un contraste bilateral, tiene que calcular el valor absoluto del Estad´ıstico antes de calcular con la cola derecha (aqu´ı es donde el Ejemplo 7.4.1 resulta relevante). El listado de ese programa, llamado Tut07-Contraste-Media-UsandoZ.R aparece en la Tabla 1 (p´ ag. 4). El fichero no funcionar´ a hasta que hayas introducido todos los datos. En particular, debes indicar con un c´ odigo num´erico (de 1 a 3) el tipo de contraste que estamos haciendo. Con esos datos se calcula el p-valor y el l´ımite de la regi´on de rechazo. Es muy conveniente leer una vez detenidamente el listado de este fichero y practicar su uso con ejemplos, de los que hayamos hecho las cuentas de forma independiente. Ficheros R para contrastes sobre los par´ ametros µ y σ de una poblaci´ on normal De la misma forma, usando los m´etodos del Cap´ıtulo 7 del curso se pueden escribir ficheros de comandos R para realizar los contrastes sobre µ y σ en poblaciones normales. Hemos incluido aqu´ı esos ficheros. Como hicimos en el caso de los intervalos de confianza, distinguimos entre el ¯ s) y el caso en el que disponemos caso en el que disponemos de los estimadores de la muestra (n, X, de todos los datos de la muestra (muestra “en bruto”). Contrastes para la media en poblaciones normales (o aprox. normales) • Muestra grande o el caso de σ conocida. ◦ Estad´ısticos de la muestra: Tut07-Contraste-Media-UsandoZ.R ◦ Datos en bruto: Tut07-Contraste-Media-UsandoZ-DatosEnBruto.R • Muestra peque˜ na ◦ Estad´ısticos de la muestra: Tut07-Contraste-Media-UsandoT.R ◦ Datos en bruto: Tut07-Contraste-Media-UsandoT-DatosEnBruto.R Contrastes para la varianza o desviaci´on t´ıpica en poblaciones normales (o aprox. normales) ◦ Estad´ısticos de la muestra: Tut07-Contraste-Varianza.R ◦ Datos en bruto: Tut07-Contraste-Varianza-DatosEnBruto.R Para practicar el uso de estos ficheros, aqu´ı tienes unos cuantos ejercicios. Ejercicios: 1. Comprueba los c´ alculos del Ejemplo 7.2.1, p´ag. 190 del curso. 2. Utiliza los datos del fichero Tut07-Contraste-Media-UsandoZ-datos.csv para contrastar la hip´ otesis nula H0 = {µ ≤ µ0 }, siendo µ0 = 7. Usa un nivel de significaci´on del 95 %. Primero puedes suponer que σ es conocida, y vale 1.5. Despu´es repite el contraste suponiendo que no conocemos σ. ¿Llegas a la misma conclusi´on? (Soluci´on en la p´ag. 6) 3. Comprueba los c´ alculos del Ejemplo 7.5.1, p´agina 204 del curso. Ten en cuenta que se hacen dos contrastes en ese ejemplo. 4. Utiliza los datos del fichero Tut07-Contraste-Media-UsandoT-datos.csv para contrastar la hip´ otesis nula H0 = {µ = µ0 }, siendo µ0 = 2.2. Usa un nivel de significaci´on del 95 %. 5. Comprueba los c´ alculos del Ejemplo 7.6.1, p´agina 206 del curso. 6. Usando los datos del fichero Tut07-Contraste-Varianza-datos.csv, contrasta (al 95 %) la hip´otesis nula H0 = σ ≥ 0.56. (Soluci´on al final de la Secci´on 1.2.1)

3

# ############################################################### # # Estadistica , Grados Biologia y Bio . Sanitaria UAH 2013 / 2014 # # Fichero de instrucciones R para calcular # un contraste de hipotesis para la media de una # poblacion normal N ( mu , sigma ) , a partir de una # muestra con n datos , su media muestral y valor de s o sigma . # # El fichero no funcionara si no introduces todos los datos . # # ############################################################### # ################################################################ # CASO : sigma conocida o desconocida , pero muestra grande n >30. # ################################################################ rm ( list = ls () ) # Numero de elementos en la muestra ( n = ) # SE SUPONE QUE LA MUESTRA ES GRANDE , salvo que se conozca sigma # Media muestral ( xbar = ) # Cuasidesviacion tipica muestral ( o sigma , si fuera conocida ) (s= ) # Valor a contrastar de la media ( aparece en la hipotesis nula ) ( mu0 = ) # ¿Que tipo de contraste estamos haciendo ? # Escribe 1 si la HIP . ALTERNATIVA es mu > mu0 , 2 si es mu < mu0 , 3 si es mu = mu0 TipoContraste = # Nivel de significacion ( nSig = ) # ############################################## # NO CAMBIES NADA DE AQU ´ I PARA ABAJO # ############################################## ( alfa =1 - nSig ) # Calculo del estadistico del contraste ( Estadistico =( xbar - mu0 ) / ( s / sqrt ( n ) ) ) # Funcion para el calculo del p - valor pValor = function ( EstadCon , tipoCon ) { if ( tipoCon ==1) { ( pV =1 - pnorm ( EstadCon ) ) } if ( tipoCon ==2) { ( pV = pnorm ( EstadCon ) ) } if ( tipoCon ==3) { pV =2 * (1 - pnorm ( abs ( EstadCon ) ) ) } return ( paste ( " El p - Valor es " ,pV , sep = " " , collapse = " " ) ) } # Funcion para el calculo del l´ ı mite de la regi´ o n de rechazo RegionRechazo = function ( alfa , tipoCon ) { if ( tipoCon ==1) { ( regionRech = paste ( " Valores del Estadistico mayores que " , qnorm (1 alfa ) ) ) } if ( tipoCon ==2) { ( regionRech = paste ( " Valores del Estadistico menores que " , qnorm ( alfa )) ) } if ( tipoCon ==3) { ( regionRech = paste ( " Valores del Estadistico mas alejados del origen que " , qnorm (1 - alfa / 2) ) ) } regionRech = paste ( " La region de rechazo la forman los " , regionRech , sep = " " , collapse = " " ) return ( regionRech ) } # Y ahora se aplican ambas funciones para mostrar los resultados pValor ( Estadistico , TipoContraste ) Estadistico RegionRechazo ( alfa , TipoContraste )

Tabla 1: Contraste de hip´ otesis para la media, pob. normal, muestras grandes (o σ conocida), con el fichero Tut07-Contraste-Media-UsandoZ.R

4

1.2.

La funci´ on t.test de R (y sus parientes)

Vamos a ver como usar la funci´ on t.test de R (que ya conocimos en la p´agina 7.1 del Tutorial06) para realizar un contraste de hip´ otesis sobre la media. Usaremos este ejemplo (que apreciar´a cualquiera que haya estado parado en un paso a nivel canadiense, esperando a que termine de pasar el mercanc´ıas de turno). Una compa˜ n´ıa ferroviaria canadiense afirma que sus trenes de mercanc´ıas no bloquean los pasos a nivel durante m´ as 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. Si llamamos µ al tiempo medio de bloqueo, queremos usar estos valores para contrastar la hip´otesis nula: H0 = {µ ≤ µ0 = 8} Y naturalmente la hip´ otesis alternativa es: Ha = {µ > µ0 = 8} Vamos a fijar un nivel de significaci´ on del 95 %, es decir, α = 0.05. Puesto que se trata de una muestra peque˜ na (n = 10), usaremos la distribuci´on t de Student para el c´alculo del p-valor. Hagamos primero los c´ alculos del contraste utilizando el fichero Tut07-Contraste-Media-UsandoT-DatosEnBruto.R, sin recurrir a t.test. Ejercicio: Haz esto. Comprueba si obtienes estos valores: n = 10, El estad´ıstico del contraste vale

¯ ≈ 9.01, ,X

s ≈ 1.631938.

¯ − µ0 X ≈ 1.957121. s √ n

Y el p-valor es: pV alor ≈ 0.04101152

De manera que, con un nivel de significaci´on 0.05 (mayor que el p-valor), tenemos evidencia emp´ırica para rechazar la hip´ otesis nula y concluir que los trenes bloquean el paso a nivel m´as tiempo del que dice la empresa. Veamos ahora como hacer este mismo contraste usando t.test (se muestra la salida): > 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="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 Como se ve, aparte del vector de datos, le hemos indicado a R el valor de µ0 y, mediante las opciones alternative = c("greater") y conf.level = 0.95, hemos seleccionado un contraste de cola derecha (greater) y el nivel de significaci´on deseado (conf.level, R usa aqu´ı la misma terminolog´ıa que para los intervalos de confianza). Si quieres hacer un contraste de otro tipo, 5

con la cola izquierda o bilateral, debes usar alternative = c("less") o bien alternative = c("two.sided"), respectivamente. La respuesta de R contiene tanto el valor del estad´ıstico de contraste en la forma t = 1.9571, como el p-valor, en p-value = 0.04101. Adem´as, para que la interpretaci´on del resultado sea m´ as f´ acil, y para que podamos comprobar que estamos haciendo lo que deseamos, R describe la hip´ otesis alternativa del contraste. Como subproducto se obtiene lo que R llama un intervalo de confianza. Ten en cuenta, en cualquier caso que nosotros no hemos visto en el curso este caso de los intervalos de confianza unilaterales. 1.2.1.

La librer´ıa TeachingDemos. Contrastes para µ y σ

Despu´es de conocer t.test seguramente te estar´as preguntado ¿y el contraste para la media con la Z, la normal est´ andar? 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˜ n´ o pensando en la ense˜ nanza?). Pero eso no significa que no est´en disponibles. Basta con cargar una librer´ıa, cuyo revelador nombre es TeachingDemos (tendr´as que instalarla previamente, claro, sin no se ha hecho previamente), y con eso ya tenemos disponible la funci´ on z-test, con la que podemos hacer esos contrastes. Vamos a usar esa funci´ on z-test para hacer el Ejercicio 2 de la p´agina 3 de este tutorial. Para la primera versi´ on, suponemos que σ es conocida y vale 1.5. Recuerda que debes seleccionar como directorio de trabajo aquel que contiene el fichero: Tut07-Contraste-Media-UsandoT-datos.csv Una vez hecho esto, vamos a mostrar el c´odigo que permite realizar el contraste (se muestra tambi´en la salida), y a continuaci´ on lo comentaremos: > muestra = scan(file="Tut07-Contraste-Media-UsandoZ-datos.csv",sep=" ",dec=".") Read 325 items > z.test(muestra, mu = 7, stdev = 1.5, alternative="greater", conf.level = 0.95) One Sample z-test data: muestra z = 2.457, n = 325.000, Std. Dev. = 1.500, Std. Dev. of the sample mean = 0.083, p-value = 0.007006 alternative hypothesis: true mean is greater than 7 95 percent confidence interval: 7.067571 Inf sample estimates: mean of muestra 7.204431 Puedes ver que el p-valor (que es aproximadamente 0.007) permite rechazar H0 . Como ves, la llamada a la funci´ on z.test incluye el argumento stdev=1.5, que representa el valor de σ, la desviaci´ on t´ıpica de la poblaci´ on, que suponemos conocida. El argumento mu=7 se usa para indicarle a R el valor que nosotros llamamos µ0 en los contrastes. Enseguida volveremos con la segunda parte de este ejercicio, en la que se supone que σ es desconocida. Pero antes, algunas preguntas sobre este primer c´alculo. Ejercicios: 1. En la salida de z.test para este ejemplo se incluye un intervalo de confianza unilateral, que es (7.204431, +∞). ¿Qu´e relaci´ on hay entre este intervalo y la regi´ on de rechazo para este tipo de contrastes, que aparece en la Ecuaci´ on 7.4 (p´ ag. 195) del curso? 2. ¿Qu´e ocurre si (sin tener en cuenta el tama˜ no de la muestra) haces este contraste usando la t de Student (con el fichero adecuado de la p´ag. 3)? ¿Qu´e p-valor obtienes? Volvamos a la segunda parte del Ejercicio 2 de la p´agina 3. Ahora ya no suponemos σ conocido 6

y por esa raz´ on tenemos que cambiar la forma en la que llamamos a z.test. La nueva versi´on es esta: > muestra = scan(file="Tut07-Contraste-Media-UsandoZ-datos.csv",sep=" ",dec=".") Read 325 items > z.test(muestra, mu = 7, sd = sd(muestra), alternative = "greater",conf.level = 0.95) One Sample z-test data: muestra z = 2.4206, n = 325.000, Std. Dev. = 1.523, Std. Dev. of the sample mean = 0.084, p-value = 0.007747 alternative hypothesis: true mean is greater than 7 95 percent confidence interval: 7.065517 Inf sample estimates: mean of muestra 7.204431 en la que, como puedes ver, hemos cambiado el argumento stdev = 1.5 por sd = sd(muestra), indic´ andole a la R que utilice la cuasidesviaci´on t´ıpica muestral en lugar de σ. El p-valor es ligeramente distinto, claro, pero no tanto como para suponer un cambio en nuestra decisi´on de rechazar H0 . Como hemos visto, la funci´ on z.test trabaja a partir de un vector de datos. Si lo que tenemos ¯ s, entonces debemos utilizar el son los estimadores (o descriptores) de una muestra, como n, X, m´etodo que vimos en la p´ agina 27 del Tutorial06, basado en la funci´on mvrnorm. Para terminar con esta visita a la librer´ıa TeachingDemos, vamos a presentar la funci´on sigma.test que, como su nombre sugiere, sirve para realizar un contraste de hip´otesis sobre la desviaci´on t´ıpica σ de una poblaci´ on normal. En el siguiente fragmento de c´ odigo R hemos usado esta librer´ıa para obtener el contraste que se ped´ıa en el Ejercicio 6 de la p´ agina 3. > muestra = scan(file="Tut07-Contraste-Varianza-datos.csv",sep=" ",dec=".") Read 250 items > sigma.test(muestra,sigma=0.56,alternative="less",conf.level=0.95) One sample Chi-squared test for variance data: muestra X-squared = 214.4613, df = 249, p-value = 0.05531 alternative hypothesis: true variance is less than 0.3136 95 percent confidence interval: 0.0000000 0.3150632 sample estimates: var of muestra 0.2701007 El modo de usar de la funci´ on sigma.test, como ves, es muy f´acil de entender. La u ´nica precauci´on que debemos tener es la de utilizar el argumento sigma= cuando la hip´otesis est´a formulada en t´erminos de la desviaci´ on t´ıpica, mientras que se usa sigmasq= cuando es la varianza σ02 la que aparece en las hip´ otesis del contraste. Por ejemplo, se obtiene exactamente el mismo resultado de antes si se usa esta otra versi´ on: sigma.test(muestra,sigmasq=0.56^2,alternative="less",conf.level=0.95) 1.2.2.

La librer´ıa asbio

En el Tutorial06 tambi´en aprendimos a usar la librer´ıa asbio para obtener intervalos de confianza. Esa librer´ıa incluye, adem´ as, funciones para algunos de los contrastes de hip´otesis que estamos viendo. Concretamente, se incluyen las dos funciones: 7

one.sample.z, para contrastes sobre µ usando Z. one.sample.t, para contrastes sobre µ usando la t de Student. Vamos a dejar que el lector explore estas dos funciones por si mismo. Una ventaja de estas funciones ¯ de asbio es que son capaces de funcionar tanto con datos en bruto, como con los valores de n, X y s. Ejercicios: 1. Lee la descripci´ on de estas dos funciones en la ayuda de la librer´ıa asbio. ´ 2. Usalas para volver a hacer los ejercicios 1 a 4 de la p´ag. 27.

1.3.

Contrastes de hip´ otesis para µ y σ (una poblaci´ on normal) usando otros programas

Calc: desaconsejamos su uso Empecemos por lo m´ as f´ acil. Calc no incluye funciones para realizar los contrastes que hemos visto, no en el sentido de que sean m´ınimamente comparables con lo que podemos hacer en R. Existen dos funciones, PRUEBA.T y PRUEBA.Z, pero s´olo las mencionamos para recomendar al lector que no gaste demasiado tiempo tratando de aprender a usarlas: son muy limitadas. Las u ´ltimas versiones de otras hojas de c´ alculo m´ as sofisticadas, como Microsoft Office 2013, incluyen algunas funciones adicionales para estos contrastes. Pero estas tareas se realizan con mucha m´as facilidad usando software estad´ıstico como R. Wolfram Alpha. Por contra, Wolfram Alpha es capaz de realizar muchos de estos contrastes, con una sintaxis bastante sencilla. Prueba a utilizar el comando: z-test for population mean y llegar´ as a un cuadro de di´ alogo en el que puedes introducir los valores concretos de la muestra, como puedes ver en la Figura 1, en la que hemos usado los valores del Ejemplo 7.2.1 del curso. Por supuesto, tambi´en existe una interfaz similar para el contraste basado en la t de Student, que encontrar´ as usando: t-test for population mean No he encontrado informaci´ on sobre una implementaci´on equivalente para el contraste de hip´otesis sobre σ, usando χ2 . Naturalmente, es posible usar Wolfram Alpha para calcular la probabilidad de una cola de la distribuci´ on χ2 , como en este ejemplo donde el comando: P[X>23] for X chi-squared with 20 dof permite calcular la probabilidad P (χ220 > 23), como se ve en la Figura 2. A partir de aqu´ı, el c´ alculo de p-valores es f´ acil, aunque m´as laborioso, claro.

8

Figura 1: Wolfram Alpha para un contraste de hip´otesis sobre la media.

Figura 2: Wolfram Alpha para c´alculos de probabilidad con χ2 .

9

2.

Intervalos de confianza y contrastes de hip´ otesis para el par´ ametro p, proporci´ on binomial

Esta secci´ on describe como obtener intervalos de confianza y realizar contrastes de hip´otesis como los que se discuten en la Secci´ on 8.1.1 del curso. Nuestro objetivo es, por lo tanto, calcular intervalos de confianza y hacer contrastes, para una variable aleatoria cuya distribuci´ on en la poblaci´on es de tipo Bernouilli, con par´ametro p. El intervalo de confianza, basado en la aproximaci´on normal a la binomial aparece en la Ecuaci´on 8.3 (p´ ag. 214) del curso: r r pˆ · qˆ pˆ · qˆ ≤ p ≤ pˆ + zα/2 . pˆ − zα/2 n n Para utilizar esta f´ ormula es muy importante que se cumplan, a la vez, estas condiciones: n > 30,

n · pˆ > 5,

n · qˆ > 5,

que, en esencia, dicen que n es bastante grande, y p y q no demasiado peque˜ nos. Vamos a ver como calcular estos intervalos y contrastes con varios de los programas que conocemos.

2.1.

Con R

Con R es muy f´ acil calcular los intervalos de confianza, usando un fichero como los que conocemos del anterior tutorial, para trabajar con los datos resumidos mediante n y pˆ, El fichero es este: Tut07-IntConf-Proporcion-UsandoZ-Estadisticos.R

En este caso no vamos a incluir un fichero para el caso en que trabajamos directamente con los datos de la muestra “en bruto” (raw data, en ingl´es). En R tambi´en se dispone de la funci´on prop.test para obtener estos intervalos de confianza para proporciones y para realizar contrastes de hip´otesis, con una sintaxis muy parecida a la que hemos visto para t.test. Vamos a aprender a usarla con un ejemplo. En una p´agina web1 ha aparecido recientemente una comparativa entre bater´ıas externas de respaldo para tel´efonos m´ oviles y aparatos similares. Adem´ as se realiz´o una encuesta entre los lectores para conocer cu´ al era la marca preferida de los lectores de esa web. De un total de 2499 opiniones, la marca ganadora obtuvo 1004 votos. Vamos a calcular un intervalo de confianza al 95 % para el porcentaje de lectores que prefieren esa marca. Para hacer esto usamos este comando en R (se muestra la salida) > prop.test(x=1004,n=2499,conf.level=0.95,correct=FALSE) 1-sample proportions test without continuity correction data: 1004 out of 2499, null probability 0.5 X-squared = 96.471, df = 1, p-value < 2.2e-16 alternative hypothesis: true p is not equal to 0.5 95 percent confidence interval: 0.3827042 0.4211188 sample estimates: p 0.4017607 Como puedes ver, la funci´ on est´ a dise˜ nada para contrates de hip´otesis, aunque produce el intervalo de confianza como subproducto. Vamos a explicar brevemente los argumentos que hemos usado para esta funci´ on. El argumento x se refiere al n´ umero de ´exitos que aparecen en la muestra, mientras que n es el tama˜ no de la muestra. El argumento conf.level es autoexplicativo, as´ı que s´ olo nos queda por comentar porque hemos incluido correct=FALSE. Si consultas la ayuda de prop.test, ver´ as que, con esa opci´ on, estamos desactivando la correcci´ on de continuidad de Yates. 1 Este

es el enlace.

10

Se trata de un ajuste similar al que hemos discutido en la Secci´on 5.5.2 (p´ag. 131) del curso. Si se usa esta correcci´ on se obtienen intervalos de confianza m´as precisos. Nosotros no lo hemos hecho para simplificar, y por eso, para obtener los resultados habituales de los problemas de libro de texto, debemos usar correct=FALSE. Ejercicios: 1. Usa el fichero Tut07-IntConf-Proporcion-UsandoZ-Estadisticos.R para comprobar las cuentas del Ejemplo 8.1.2 (p´ ag. 214 del curso), el de los araos embridados. 2. Haz lo mismo con prop.test. Y calcula un intervalo de confianza para los datos de araos embridados del a˜ no 2008 (138 embridados, 270 no embridados). 3. Comprueba las cuentas del contraste de hip´otesis que aparece en el Ejemplo 8.1.3 (p´ag. 215). 4. Seg´ un el Bar´ ometro del CIS de enero del 2013, en una muestra de 2452 personas residentes en Espa˜ na, 568 de ellas se declararon no creyentes. Usa estos datos para contrastar, al 95 %, la hip´ otesis de que la proporci´ on de no creyentes es superior al 20 %. 5. ¿Por qu´e no incluimos un fichero para trabajar con datos en bruto? ¿C´omo ser´ıa el fichero del Ejemplo 8.1.2?

2.2.

Con otros programas

Hemos incluido un fichero-plantilla de Calc para realizar estas mismas operaciones: Tut07-IntervaloConfianzaProporcionMuestraGrande.ods

Con Wolfram Alpha tambi´en es muy f´acil calcular estos intervalos. En el Ejemplo 8.1.2 de los araos, era n = 456, y el n´ umero de ´exitos (araos embridados) era de 139. As´ı que pˆ =

139 ≈ 0.3048, 456

y qˆ =

317 ≈ 0.6952 456

Para calcular un intervalo de confianza para p a partir de estos datos, en Wolfram Alpha podmeos escribir: binomial confidence interval n=456, p-hat=139/456 O tambi´en binomial confidence interval n=456, p-hat=0.3048 o incluso, usando el n´ umero de ´exitos (successes): binomial confidence interval n=456, number of successes =139 Mostramos, en la siguiente figura, (parte de) el resultado del u ´ltimo de estos comandos:

11

Finalmente, si lo que deseamos es hacer un contraste sobre la proporci´on, podemos usar el comando proportion hypothesis test para llegar a un interfaz en el que introducir los valores necesarios para el contraste, que se muestra en esta figura:

3.

Variables cualitativas (factores) en R

En los u ´ltimos cap´ıtulos del curso, y en los tutoriales precedentes, nuestro trabajo se ha centrado en el an´ alisis de valores de variables cuantitativas; es decir, n´ umeros. Pero sabemos que a veces es necesario trabajar con variables cualitativas, que tambi´en hemos llamado factores (y sus valores se llaman niveles). Y de hecho esas variables van a tener un papel protagonista en varios de los cap´ıtulos de la cuarta parte del curso. Recuerda que una variable cualitativa se usa para establecer clasificaciones nominales en los datos. Por ejemplo, la clasificaci´ on en hombre o mujer entre los pacientes que siguen un cierto tratamiento. O la clasificaci´ on taxon´ omica por especies. Es cierto que siempre podr´ıamos codificar las variables cualitativas mediante n´ umeros. Pero no es menos cierto que lo m´ as c´ omodo (y prudente) es poder utilizar nombres como valores de las variables. Ya hemos visto (en la Secci´on 4 del Tutorial04) que R nos permite crear un cierto tipo de valores, concretamente los valores de tipo character, que son simplemente palabras o frases entrecomilladas (en general las llamamos cadenas de texto). Recomendamos una relectura r´apida de esa Secci´ on antes de seguir adelante. Por ejemplo, este vector podr´ıa representar el g´enero de los diez pacientes que se han sometido a un cierto tratamiento: pacientesPorGenero=c("mujer","mujer","hombre","mujer","hombre","hombre","hombre", "mujer","mujer","mujer") class(pacientesPorGenero) 12

Cuando ejecutes este c´ odigo,obtendr´ as como u ´ltimo resultado > class(pacientesPorGenero) [1] "character" Es decir, que los elementos del vector pacientesPorGenero son, para R, de tipo character. Ese es el tipo de dato que se utiliza en R para representar los valores de una variable cualitativa. Naturalmente, si intentas realizar operaciones num´ericas con ese vector, como estas, pacientesPorGenero^2 mean(pacientesPorGenero) R te obsequiar´ a con un surtido de insultos m´as o menos ofensivos en la consola de comandos. En el segundo caso (el intento fallido de calcular la media), el valor que devuelve R es interesante: se obtiene NA, que es el valor que R devuelve cuando un resultado num´erico es imposible de calcular, o no est´ a disponible (Not Available, de ah´ı el nombre). Queremos recordar tambi´en, que en el Tutorial04 vimos que los data frames de R son las estructuras de datos adecuadas para representar tablas en las que se mezclan datos cuantitativos y cualitativos (las matrices, sin embargo, almacenan datos que son todos del mismo tipo).

3.1.

Factores

Una situaci´ on 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´ umeros), pero esos valores aparecen agrupados de manera natural. Por ejemplo, cuando estamos midiendo la respuesta X (un n´ umero) de una serie de pacientes frente a varios tratamientos, es evidente que lo natural es agrupar los resultados seg´ un el tratamiento empleado. Una manera de hacer esto, en R, es usar un data.frame que contenga los valores de X, junto con los valores de una variable que indique el tipo de tratamiento empleado. Vamos a llamar T1 , T2 , . . . , Tk a los distintos tratamientos. La terminolog´ıa est´andar en Estad´ıstica (que ya vimos en la Secci´on 1.1.1, p´ ag. 5 del curso) consiste en decir que el tratamiento T es un factor, y que sus valores T1 , . . . , Tk son los niveles del factor. Podr´ıamos entonces pensar en recoger los resultados en un fichero de datos como el fichero adjunto, Tut07-Tratamiento-01.csv. Abre el fichero primero en un editor de texto (como el Bloc de Notas, en Windows), para hacerte una idea de su estructura. Podr´ıamos trabajar directamente con los nombres de los tratamientos, usando en R variables de tipo character para los factores. Pero eso tendr´ıa un impacto negativo en el rendimiento de R, porque manejar esas cadenas de caracteres consume muchos recursos de memoria y tiempo del procesador. Y, dado que los factores son omnipresentes en Estad´ıstica, los creadores de R han optado por incluir un tipo especial de datos, el tipo factor, para representar los factores y sus niveles. De hecho, si leemos un fichero como el anterior usando read.table (que tambi´en vimos en el Tutorial04), el comportamiento por defecto de R es convertir las variables cualitativas en factores. Vamos a ver esto en funcionamiento. Aseg´ urate de haber seleccionado el directorio de trabajo adecuado (el que contiene el fichero csv) y ejecuta estos comandos: experimento=read.table(file="Tut07-Tratamiento-01.csv",header =T, sep=" ") class(experimento$Respuesta) class(experimento$Tratamiento) Lo que hemos hecho aqu´ı es guardar el contenido del fichero en un data.frame llamado experimento, con dos campos, cuyos tipos hemos mostrado usando class: experimento$Respuesta, de tipo num´erico que almacena los valores de la variable X (la respuesta individual de cada uno de los pacientes). experimento$Tratamiento, de tipo factor que almacena los valores de la variable T (el tipo de tratamiento empleado). 13

Recuerda que, para acceder a las distintas variables de un data.frame (piensa en ellas como las columnas de una tabla), utilizamos la notaci´on $, entre el nombre del data.frame y el de la variable. En R puedes ver (en una ventana nueva) el contenido del data.frame experimento usando el comando View(experimento) Esto es especialmente interesante cuando el data.frame es muy grande, y verlo en la consola no resulta pr´ actico.

¿C´ omo se usan los factores en R? La contestaci´on, en este tutorial, 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´as sencillo. ¿Qu´e aspecto tiene un vector de tipo factor? En el caso del data.frame experimento que acabamos de crear, al mostrar el vector experimento$Tratamiento se obtiene esto: > View(experimento) > 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 [24] T4 T1 T4 T1 T2 T3 T5 T5 T4 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´olo queremos mostrar los niveles del factor, podemos usar la funci´on levels, cuya salida es esta: > levels(experimento$Tratamiento) [1] "T1" "T2" "T3" "T4" "T5" ¿Qu´e 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´on as.character: > as.character(experimento$Tratamiento) [1] "T1" "T1" "T4" "T1" "T2" "T5" "T5" "T2" "T5" "T2" "T3" "T1" "T5" 14

[14] "T3" "T3" "T4" "T1" "T4" "T1" "T5" "T4" "T3" "T4" "T4" "T1" "T4" [27] "T1" "T2" "T3" "T5" "T5" "T4" "T2" "T4" "T4" "T5" "T1" "T3" "T3" [40] "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´ odigos num´ericos. Al mostrarlos, a petici´on nuestra, el uso de las comillas permite distinguir entre character (con comillas) y factor (sin comillas). Como hemos dicho, cuando R usa una funci´on como read.table para cargar un data.frame, los campos alfanum´ericos (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´on stringsAsFactors=FALSE en la funci´on read.table.

3.2.

Usando los factores para explorar los datos

Los factores permiten modular el comportamiento de algunas funciones de R, de manera que la informaci´ on que se obtiene tiene en cuenta esa organizaci´on en niveles de los valores que estamos analizando. Por ejemplo, en la Secci´on 5 del Tutorial03 vimos como usar la funci´on summary de R para obtener un resumen descriptivo de un vector de datos (media y percentiles, b´asicamente). Ahora, que tenemos un data.frame con una clasificaci´on por tratamientos, nos gustar´ıa seguramente obtener esa misma informaci´on, 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, puesto que ya hemos aprendido a seleccionar los elementos de un vector mediante condiciones (usando los corchetes [ ]), podemos ir separando cada uno de los grupos de tratamiento, y aplicarles la funci´ on summary a los vectores resultantes. Ser´ıa algo como esto (mostramos, por ejemplo el segundo grupo de tratamiento): > (tratamiento2=experimento[experimento$Tratamiento=="T2",1]) [1] 3.08 2.60 3.14 3.36 3.31 3.83 3.18 > summary(tratamiento2) Min. 1st Qu. Median Mean 3rd Qu. Max. 2.600 3.110 3.180 3.214 3.335 3.830 Esto funciona, pero como se puede apreciar, trabajar as´ı es engorroso. Para obtener lo que queremos, hay una manera mucho m´ as natural de proceder, usando la funci´on tapply. Veamos como funciona (primero hemos usado attach para poder referirnos a las variables Respuesta y Tratamiento sin tener que usar el $): > attach(experimento) > tapply(Respuesta,Tratamiento, summary) $T1 Min. 1st Qu. Median Mean 3rd Qu. 4.700 4.955 5.080 5.079 5.228 $T2 Min. 1st Qu.

Median

Mean 3rd Qu.

Max. 5.370

Max. 15

2.600

3.110

3.180

$T3 Min. 1st Qu. 4.520 5.012

3.214

3.335

3.830

Median 5.290

Mean 3rd Qu. 5.229 5.478

Max. 5.710

$T4 Min. 1st Qu. 4.400 4.605

Median 4.975

Mean 3rd Qu. 5.058 5.400

Max. 6.060

$T5 Min. 1st Qu. 4.560 4.935

Median 5.015

Mean 3rd Qu. 5.039 5.182

Max. 5.390

Ahora s´ı: como ves, R ha reconocido esa estructura en niveles del factor Tratamiento, y no est´a describiendo la variable Respuesta para cada uno de los grupos de tratamiento que conforman nuestros datos. El resultado es una tabla (de ah´ı la t inicial de tapply) que contiene, ordenados por filas, los resultados de aplicar la funci´ on summary a los valores de Respuesta, agrupadas seg´ un los distintos niveles del factor Tratamiento. Otro ejemplo de las ventajas del lenguaje de factores se obtiene al pensar en la funci´on boxplot. Si aplicamos esa funci´ on directamente al data.frame, ejecutando: boxplot(experimento) se obtiene un gr´ afico que es, de hecho un sinsentido:

F´ıjate en que se incluye un boxplot para la variable Tratamiento ¡qu´e es cualitativa! Evidentemente esto no es correcto. Para obtener lo que queremos basta con escribir: boxplot(Respuesta~Tratamiento) y el resultado es un gr´ afico mucho m´ as adecuado para comparar los resultados de los tratamientos:

16

La sintaxis Respuesta~Tratamiento que hemos usado por primera vez 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´on entre variables (lo obtienes, en un teclado est´andar de PC, pulsando las teclas AltGr y 4 a la vez). En este caso, lo que estamos estudiando es si hay alguna relaci´on entre la variable cualitativa (factor) Tratamiento, y la variable cuantitativa Respuesta. En particular, una pregunta que a menudo nos interesa consiste en saber si la respuesta media al tratamiento es distinta seg´ un el grupo de tratamiento que se considere. Este tipo de problemas, en los que se investiga la relaci´on entre varias variables aleatorias (de distintos tipos), forman el contenido de la cuarta parte del curso, y por eso estamos empezando a preparar el terreno en este tutorial.

3.3.

La funci´ on cut

Otra de las primeras funciones de R que se aprenden al comenzar a trabajar con factores es la funci´ on cut, cortar en ingl´es. Y el nombre describe muy bien lo que hace esta funci´on que se usa para cortar en piezas los datos de una variable cuantitativa, con el fin de agruparlos en clases o intervalos, como discutimos en la Secci´on 1.1.3 del Cap´ıtulo 1 del curso. Por ejemplo, en el fichero adjunto Tut07-Edades.csv contiene las edades, en a˜ nos, de un grupo de personas. Imag´ınate que queremos agrupar esos datos, seg´ un la edad, en estas categor´ıas: Menores, con edades menores que 18 a˜ nos, Jovenes, con edades entre 18 y 30 a˜ nos, MedianaEdad, con edades entre 31 y 64 a˜ nos, Mayores, con edades de 65 a˜ nos en adelante. Vamos a ver como usar cut para hacer eso en R. Primero leemos el fichero, que debe estar en el 17

directorio de trabajo, y lo guardamos en un vector de datos, llamado edades, al que aplicaremos la funci´ on cut. edades=scan(file="Tut07-Edades.csv") Pero antes de cortar, tenemos que pensar bien donde vamos a dar los cortes. Los intervalos que R usa, en la funci´ on cut, son, por defecto, de la forma (a, b], incluyendo el l´ımite derecho pero no el izquierdo (esto se puede cambiar). Eso implica que so hacemos una elecci´on poco cuidadosa de los puntos de corte, como en: cortesEdad=c(0,18,30,65,100) terminaremos con un intervalo de edades como 18 < edad ≤ 30, que no se corresponde con lo que queremos, porque no incluye a las personas de 18 a˜ nos. Algo parecido sucede con el intervalo 30 < edad ≤ 65, que no es lo que queremos porque incluye a las personas de 65 a˜ nos. Una elecci´on algo m´ as meditada nos conduce a: cortesEdad=c(0,17,30,64,100) Esto est´ a casi bien, pero a´ un tenemos un problema. Para descubrirlo haz este: Ejercicio ¿Cu´ ales son los intervalos (a, b] determinados por este segundo vector cortesEdad?

¡No sigas, si no has hecho este ejercicio!

18

Como habr´ as podido ver, el primer intervalo es (0, 17]. Eso significa que si alguna de las edades es 0 (como, de hecho, sucede), no quedar´a incluida en ese intervalo. Para evitar ese problema, la funci´ on cut dispone de una opci´ on, include.lowest=TRUE, que permite cerrar a la izquierda el primero de los intervalos. Con eso estamos listos para usar la funci´on cut: factorEdad = cut(edades, breaks=cortesEdad, include.lowest=TRUE) El resultado, que hemos guardado en la variable factorEdad es un vector, de tipo factor (puedes usar class(factorEdad) para comprobarlo), que contiene, para cada elemento de edades el intervalo de edades (nivel del factor) al que pertenece. Para que lo veas con m´as claridad, vamos a ver, juntos, los primeros elementos de los vectores edades y factorEdad, usando head: > head(edades) [1] 62 89 13 15 36 17 > head(factorEdad) [1] (30,64] (64,100] [0,17] [0,17] (30,64] Levels: [0,17] (17,30] (30,64] (64,100]

[0,17]

Como ves, en cada posici´ on de factorEdad est´a el intervalo de edades (nivel del factor) al que pertenece el correspondiente elemento de edades. F´ıjate en que el intervalo del primer elemento es [0,17], cerrado en ambos extremos. Adem´ as, R nos recuerda al final el conjunto de nombres de los valores (niveles) que puede tomar el factor. Ese conjunto de valores de un factor se puede obtener usando la funci´on levels as´ı: > levels(factorEdad) [1] "[0,17]" "(17,30]"

"(30,64]"

"(64,100]"

Hemos dicho que estos son los nombres de los niveles, para recordar que, internamente, R utiliza c´ odigos num´ericos para trabajar con los niveles. Por eso es bueno usar una palabra como etiquetas para referirse a esos nombres (como en ingl´es, donde se usa labels). A menudo sucede que queremos cambiar estas etiquetas, para reemplazarlas con otras m´as u ´tiles para nosotros. La forma de hacerlo, en este ejemplo, es esta: levels(factorEdad)=c("menor","joven","medianaEdad","mayor") Ahora, si pruebas a ejecutar de nuevo head(factorEdad), ver´as el efecto que ha tenido ese cambio. Para lo que resta de trabajo, vamos a devolver al factor a sus etiquetas iniciales, que son m´ as concisas: levels(factorEdad)=c("[0,17]",

"(17,30]",

"(30,64]",

"(64,100]")

Podemos, con ayuda del factor, obtener muy f´acilmente una tabla de frecuencias de las edades (se muestra la salida): > table(factorEdad) factorEdad [0,17] (17,30] (30,64] (64,100] 63 79 230 128 y los res´ umenes (summary) o boxplots por niveles, etc., que ya hemos visto antes.

Muchas gracias por la atenci´ on.

19