martes, 21 de agosto de 2012

Clase 4: Contrastes de hipótesis en R-project

Inferencia: Contraste de hipótesis


Continuando con el breve curso de R que he presentado (Clase 1, Clase 2, Clase 3), en esta entrada trataremos la Estadística Inductiva o Inferencia Estadística es un conjunto de métodos que se fundamentan en la Teoría de la Probabilidad y que tienen por finalidad generalizar los resultados, obtenidos mediante una muestra, a toda una población. Es decir, se trata de procedimientos para aceptar o rechazar una hipótesis (una afirmación de la realidad) acerca de un parámetro o característica poblacional a partir de una muestra de la población.
Los contrastes pueden ser paramédicos o no paramédicos según asuman una determinada distribución de los datos o no. Los contrastes presentan una hipótesis nula H0 (la que queremos contrastar) y una hipótesis alternativa H1 (lo que ocurre si no es cierta la hipótesis nula). Se dice que un contraste es significativo cuando se rechaza H0; y es no-significativo cuando se acepta H0 (o mejor dicho, cuando no rechazamos H0).

La elección de la prueba de contraste que debemos aplicar depende de nuestros datos y las preguntas que nos hacemos sobre ellos. Aquí les dejo un esquema muy ilustrativo sobre ello:


Inferencial


Contrastes paramétricos

Asumen la distribución Normal de los datos.
Existen contrastes para cualquiera de los parámetros poblacionales: media, varianza, proporción; para 1 o para más muestras.


Contrastes no paramétricos

No asumen Normalidad.


Aplicación en R

Contraste para el promedio de una muestra

Ejemplo:

data(sleep)
sleep
 
#contraste de una muestra: promedio o media
t.test(sleep$extra)
t.test(sleep$extra,mu=1.4, conf.level = 0.95) #two.sided
 t.test(sleep$extra,mu=1,alternative="greater")
 t.test(sleep$extra,mu=1,alternative="less")
 
shapiro.test(sleep$extra) #para probar si existe normalidad en los datos, es decir, si los datos se distribuyen según una Normal.
wilcox.test(sleep$extra,alternative="less",mu=1)  #sin asumir normalidad
 
 

Contraste para la proporción de una muestra
 
#contraste para una muestra: proporción
heads <- rbinom(1, size=100, prob = .5)
prop.test(heads, 100)
 
#contraste para dos muestras: igualdad de proporción
smokers  <- c( 83, 90, 129, 70 )
patients <- c( 86, 93, 136, 82 )
prop.test(smokers, patients)



Contraste de igualdad de medias


Ejemplo
#contraste de dos muestras: igualdad de medias
t.test(extra ~ group, data = sleep, var.equal=TRUE) #t.test para muestras independientes
#t.test para muestras relacionadas
wilcox.test(extra ~ group, data = sleep, alternative = "less")   #sin asumir normalized



Contraste para la igualdad de varianzas

Ejemplo
#contraste de dos muestras: cociente de varianzas
var.test(extra ~ group, data = sleep, alternative = "two.sided", conf.level = 0.95)
levene.test(acero$pr.galv1, acero$averias) #sin asumir normalidad
bartlett.test(consumo ~ temperatura, data = acero) #para más de 2 muestras




Contraste de medias para más de dos muestras:  




> Anova1 <- aov(consumo ~ temperatura, data = acero) > summary(Anova1)


#comparaciones múltiples
> comparacion <- glht(Anova1, linfct = mcp(temperatura = "Tukey")) > summary(comparacion)


#si no hay homocedasticidad
kruskal.test(consumo ~ linea, data = acero)


Aplicación en R mediante el paquete ctest


###########################################################################
#                 Analisis de datos con R                                 #
#                           2012                                          #
#                        Msc. Rosana Ferrero                              #
#             http://statisticalecology.blogspot.com/                     #
###########################################################################
# 4. Contrastes de hip?tesis                                       #
###########################################################################
 
#Existe en R una librer?a de funciones para realizar contrastes estad?sticos, llamada ctest.
#Para cargarla, escribiremos:
 library(ctest)
#Dicha librer?a contiene bastantes funciones, pero s?lo nos vamos a fijar en algunas que implementan
#contrastes muy com?nmente utilizados. Las funciones que vamos a ver son las siguientes:
 
    # binom.test  Test exacto sobre el par?metro de una binomial
    # cor.test  Test de asociaci?n entre muestras apareadas
    # wilcox.test  Test de suma de rangos de Wilcoxon para una y dos muestras
    # prop.test  Test de igualdad de proporciones
    # chisq.test  Test de la chi-cuadrado para datos de conteo
    # fisher.test  Test exacto de Fisher para datos de conteo
    # ks.test  Test de Kolmogorov-Smirnov para ajuste de datos a distribuciones dadas
    # shapiro.test  Test de Shapiro para comprobar ajuste de datos a una distribuci?n normal
    # oneway.test  Test para comprobar la igualdad de medias entre varios grupos de datos
    # var.test  Test para comprobar la igualdad de varianzas entre dos grupos de datos
 
#Vamos a poner ejemplos de cada contraste, para ver qu? tipo de problemas resuelven y tener una base para su aplicaci?n.
 
#### binom.test ################################################################
#Lleva a cabo un contraste exacto (no aproximado) sobre el valor de la probabilidad de ?xito en un experimento de Bernoulli.
 
  #Por ejemplo, tiramos una moneda al aire 600 veces, y nos salen 300 caras. Y queremos contrastar la hip?tesis de que la moneda est? equilibrada, esto es, la probabilidad de que p=0.5. Para ello diremos cu?ntas veces sali? cara y cu?ntas veces sali? cruz:
   binom.test(c(300,300),p=0.5)
    #        Exact binomial test
    # data:  300 and 600
    # number of successes = 300, number of trials = 600, p-value = 1
    # alternative hypothesis: true probability of success is not equal to 0.5
    # 95 percent confidence interval:
    #  0.4592437 0.5407563
    # sample estimates:
    # probability of success
    #                    0.5
   ## Como el p-valor ha salido 1 (es lo m?s alto posible), no rechazaremos la hip?tesis nula. Esto es, daremos por bueno que la moneda estaba equilibrada. Observamos que tambi?n obtenemos un intervalo de confianza al 95% para el verdadero valor de la probabilidad de obtener cara. Dicho intervalo es (0.459, 0.540) en nuestro ejemplo.
 
   #Otro ejemplo: tiramos la moneda al aire y obtenemos 30 caras y 100 cruces:
    binom.test(c(30,100),p=0.5)
    #        Exact binomial test
    #  data:  30 and 130
    #  number of successes = 30, number of trials = 130, p-value = 5.421e-10
    #  alternative hypothesis: true probability of success is not equal to 0.5
    #  95 percent confidence interval:
    #   0.1614375 0.3127614
    #  sample estimates:
    #  probability of success
    #               0.2307692
   ##En este caso, el p-valor es pr?cticamente cero, lo que quiere decir que debemos rechazar la hip?tesis de que la moneda est? equilibrada. Como el intervalo de confianza resulta ser (0.16, 0.31), concluiremos que la probabilidad de obtener cara para la moneda utilizada est? probablemente dentro de dicho intervalo.
 
#### prop.test #################################################################
#Con esta funci?n contrastaremos si las proporciones en varios grupos son las mismas, o bien que dichas proporciones equivalen a unas determinadas.
 
  #Para empezar, tiraremos una moneda al aire 100 veces, y nos preguntaremos si la proporci?n de caras puede ser considerada el 50% :
       caras = rbinom(1, size=100, pr = .5)
       caras
      #[1] 51
       prop.test(caras,100)
      #  1-sample proportions test with continuity correction
      # data:  caras out of 100, null probability 0.5
      # X-squared = 0.01, df = 1, p-value = 0.9203
      # alternative hypothesis: true p is not equal to 0.5
      # 95 percent confidence interval:
      # 0.4086512 0.6105719
      # sample estimates:
      #   p
      # 0.51
     ## Como el p-valor es muy alto, 0.92, no rechazaremos la hip?tesis nula de que p=0.5. En este caso, el test funciona como en el caso de binom.test.
 
      #Supongamos ahora que queremos saber si dos muestras han salido de la misma poblaci?n.
      #Por ejemplo, queremos saber si las mujeres estudian un grupo de carreras diferentes en la misma proporcion. Para ello miramos la lista de matriculados en cuatro carreras y comprobamos si las proporciones de mujeres en las mismas coinciden:
        estudiantes = c(100,200,50,500)
        names(estudiantes) = c("quimica","derecho","matematicas","economia")
        estudiantes
       #    quimica     derecho matematicas    economia
       #        100         200          50         500
        estudiantes.mujeres = c(60,120,10,200)
        prop.test(estudiantes.mujeres,estudiantes)
       #      4-sample test for equality of proportions without continuity correction
       # data:  estudiantes.mujeres out of estudiantes
       # X-squared = 44.5373, df = 3, p-value = 1.160e-09
       # alternative hypothesis: two.sided
       # sample estimates:
       # prop 1 prop 2 prop 3 prop 4
       #   0.6    0.6    0.2    0.4
      ##   La conclusion es que no. Las mujeres, en este ejemplo, se matriculan de las carreras en distinta proporcion. Es aparente que se matriculan menos en Matematicas, en este ejemplo, pero el uso del contraste nos asegura con un alto grado de seguridad que esas diferencias observadas no son debidas al azar, sino que son estructurales.
 
 
#### chisq.test ################################################################
# Este contraste se utiliza para tablas de contingencia. El caso habitual ser? disponer de una tabla de contingencia que cruza dos variables, contando el n?mero de coincidencias seg?n las categor?as de ambas variables, y se trata de ver si las variables son o no son independientes.
 
   # Para poner un ejemplo, supongamos que tenemos una muestra de personas en las que hemos observado dos variables: el color de ojos, que puede ser "claro" y "oscuro", y el color de pelo, que puede ser "rubio" ? "moreno". Entonces hacemos una tabla de contingencia para cruzar dichas caracter?sticas, y aplicamos el test para ver si el "color de ojos" y el "color de pelo" son variables independientes:
   # Suponemos que los datos han sido recogidos previamente en una variable llamada 'datos':
     datos
    #        claro oscuro
    #  rubio     20      2
    #  moreno     3     18
     chisq.test(datos)
    #  Pearson's Chi-squared test with Yates' continuity correction
    #     data:  datos
    # X-squared = 22.3693, df = 1, p-value = 2.249e-06
    ##  Como el p-valor es pr?cticamente cero, concluiremos que hay que rechazar la hip?tesis nula, que es que hay independencia entre las variables observadas: el color de ojos y el color de pelo guardan correlaci?n en estos datos.
 
    #Repitamos el ejemplo, construyendo la tabla de contingencia para repasar c?mo se hace, para un conjunto de datos distinto y m?s peque?o que el anterior:
     color.ojos = c("claro","oscuro","oscuro","oscuro","claro","claro","claro", "oscuro","oscuro","claro","claro","claro","claro","oscuro","oscuro")
     color.pelo = c("rubio","moreno","moreno","rubio","moreno","rubio","rubio", "moreno","rubio","moreno","rubio","rubio","rubio","moreno","moreno")
     table(color.ojos,color.pelo)
    #          color.pelo
    # color.ojos moreno rubio
    #    claro       2     6
    #    oscuro      5     2
     chisq.test(table(color.ojos,color.pelo))
    #         Pearson's Chi-squared test with Yates' continuity correction
    #    data:  table(color.ojos, color.pelo)
    # X-squared = 1.637, df = 1, p-value = 0.2007
    #     Warning message:
    # Chi-squared approximation may be incorrect in: chisq.test(table(color.ojos, color.pelo))
   ##    En este caso, el p-valor es 0.2, que es muy alto: no rechazaremos la hip?tesis de que las variables color de ojos y color de pelo son independientes. El mensaje de aviso que da el test es porque el n?mero de casos utilizado es muy peque?o.
 
#### cor.test y wilcox.test ###############################################################
# Son dos tests para ver si hay asociaci?n entre muestras de datos. Por ejemplo, una serie de alumnos pueden ser sometidos a dos test de inteligencia distintos. Si hay 100 alumnos, tendr?amos dos vectores de longitud 100: el primero ser?a la puntuaci?n alcanzada por los alumnos en el primer test, y el segundo la puntuaci?n de dichos alumnos para el segundo test. Los contrastes cor.test y wilcox.test servir?an para ver si hay correlaci?n entre ambos tests.
   #  Como advertimos al principio del tema, se supone que el fundamento del contraste se conoce. De todos modos, digamos que el test de Wilcoxon es no param?trico, y tiene m?s en cuenta el orden correlativo de los datos que el valor de los datos mismos. Para apreciar las diferencias, pongamos unos ejemplos:
   #    Primero generamos una secuencia de datos, y a partir de la misma, otras dos, relacionadas de alguna manera, y veremos qu? dicen los tests al respecto:
      x = seq(-5,5)
      x
     #   [1] -5 -4 -3 -2 -1  0  1  2  3  4  5
     y = x^2
     y
     #  [1] 25 16  9  4  1  0  1  4  9 16 25
 
    # Vemos que entre x e y hay una relaci?n, aunque sea no lineal
      z = x^3
      z
     # [1] -125  -64  -27   -8   -1    0    1    8   27   64  125
 
    # Vemos si existe relaci?n entre x e y
       cor.test(x,y)
      #        Pearson's product-moment correlation
      # data:  x and y
      # t = 0, df = 9, p-value = 1
      # alternative hypothesis: true correlation is not equal to 0
      # sample estimates:
      # cor
      #   0
 
    # La relaci?n entre x e y no es detectada por el test cor.test: el p-valor ha salido 1.
       wilcox.test(x,y)
      #              Wilcoxon rank sum test with continuity correction
      #      data:  x and y
      # W = 17.5, p-value = 0.005106
      # alternative hypothesis: true mu is not equal to 0
      #       Warning message:
      # Cannot compute exact p-value with ties in: wilcox.test.default(x, y)
     ## La relaci?n entre x e y s? ha sido detectada por wilcox.test.
 
    # Vemos si existe relaci?n entre x y z
     cor.test(x,z)
      #            Pearson's product-moment correlation
      #    data:  x and z
      # t = 7.1257, df = 9, p-value = 5.51e-05
      # alternative hypothesis: true correlation is not equal to 0
      # sample estimates:
      #      cor
      # 0.921649
      ## En este caso, cor.test s? detecta la relaci?n, al mantenerse el orden de los datos.
 
 
#### fisher.test ###############################################################
# Este test tiene el mismo fin que el test  chisq.test, pero otro fundamento, y es ?til para muestras peque?as, caso para el cual chisq.test no es adecuado.
 
    # El ejemplo que trae la librer?a es muy adecuado. Cierta dama fina (inglesa, por supuesto), era aficionada a tomar te con leche (?c?mo no!), y presum?a de ser capaz de distinguir cuando su criada le echaba primero la leche y luego el t?, o bien echaba primero el t? y despu?s la leche. Tal habilidad era mirada con incredulidad por sus allegados. Para salir de dudas, Fisher le dio a probar a la dama una serie de tazas de t?, en las que unas veces se hab?a echado primero la leche y otras primero el t?. Fisher apunt? los aciertos de la dama y con la tabla de contingencia as? construida, veamos lo que pas?....
       TeaTasting =  matrix(c(3, 1, 1, 3), nr = 2, dimnames = list(Guess = c("Milk", "Tea"),  Truth = c("Milk", "Tea")))
       fisher.test(TeaTasting, alternative = "greater")
         # Fisher's Exact Test for Count Data
         #   data:  TeaTasting
         #   p-value = 0.2429
         #   alternative hypothesis: true odds ratio is greater than 1
         #   95 percent confidence interval:
         #    0.3135693       Inf
         #   sample estimates:
         #   odds ratio
         #     6.408309
        ## Dado el p-valor obtenido, el test no mostr? evidencia suficiente a favor de la habilidad de la dama. Sin embargo, vemos que de 4 veces que se le dio primero la leche, acert? 3, y lo mismo para el t?.
 
    # Otro ejemplo, con los datos de color de pelo y ojos que hemos utilizado m?s arriba:
       fisher.test(table(color.ojos,color.pelo))
      #        Fisher's Exact Test for Count Data
      # data:  table(color.ojos, color.pelo)
      # p-value = 0.1319
      # alternative hypothesis: true odds ratio is not equal to 1
      # 95 percent confidence interval:
      #  0.007811653 1.922966759
      # sample estimates:
      # odds ratio
      # 0.1563219
     ## En este caso, tampoco somos concluyentes, como no lo ?ramos al utilizar la funci?n chisq.test. A pesar de todo, para el tama?o de la muestra, que es peque?o, este test es m?s fiable.
 
#### ks.test ##################################################################
# Utilizaremos este contraste para comprobar si dos conjuntos de datos siguen la misma distribuci?n, o bien para ver si un determinado conjunto de datos se ajusta a una distribuci?n determinada. Las muestras no necesitan ser del mismo tama?o.
 
      # Por ejemplo, vamos a generar 50 datos de una distribuci?n normal, y 30 de una uniforme, a ver si el test es capaz de advertirnos que las distribuciones de origen no son las mismas:
       x = rnorm(50)
       y = runif(30)
       ks.test(x, y)
      #              Two-sample Kolmogorov-Smirnov test
      # data:  x and y
      # D = 0.52, p-value = 2.903e-05
      # alternative hypothesis: two.sided
      # Como el p-valor es pr?cticamente cero, rechazaremos la hip?tesis de que los datos vienen de la misma distribuci?n.
 
#### shapiro.test ##############################################################
# Es como ks.test, pero est? especializado en la distribuci?n normal.
 
      # Por ejemplo, con los mismos datos x generados en el ejemplo anterior, comprobemos si el test detecta que provienen de una distribuci?n normal:
       shapiro.test(x)
        #        Shapiro-Wilk normality test
        # data:  x
        # W = 0.9758, p-value = 0.3905
       ##  En efecto, el p-valor es grande (0.39), y no rechazaremos la hip?tesis nula de que los datos x provienen de una normal.
 
#### oneway.test ###############################################################
#Este es el contraste de igualdad de medias que se estudia en los modelos ANOVA. Se supone que hay varios grupos de datos, que cada grupo sigue una distribuci?n normal (no necesariamente con la misma varianza), y se trata de ver si las medias son todas iguales o no. Este test aparece mucho en procesos de calidad, de dise?o de experimentos, en los que se fabrica un producto de varias maneras distintas, y se trata de ver si los procedimientos conducen al mismo resultado o no.
 
    # Para poner un ejemplo, consideraremos los datos de fabricaci?n de camisas de trabajo sacados del libro de Dise?o y Analisis de Experimentos de Montgomery, del grupo Editorial Iberoamericana. Se trata de fabricar camisas de trabajo que tengan la mayor resistencia posible a la rotura. Para ellos se fabrican 5 tipos de camisas, donde en cada grupo se da una proporci?n determinada de algod?n en la composici?n de la camisa. Las camisas del grupo 1 tienen un 15% de algod?n, las del grupo 2 un 20%, las del grupo 3 un 25%, las del grupo 4 un 30%, y las del grupo 5 un 35% de algod?n. Se fabrican 5 camisas de cada tipo, se rompen y se mide su resistencia. Se trata de determinar si la resistencia media de las camisas de cada grupo es la misma o no. Vamos all?:
       algodon = scan()      # introducimos los datos por teclado
          1: 7
          2: 7
          3: 15
          4: 11
          5: 9
          6: 12
          7: 17
          8: 12
          9: 18
          10: 18
          11: 14
          12: 18
          13: 18
          14: 19
          15: 19
          16: 19
          17: 25
          18: 22
          19: 19
          20: 23
          21: 7
          22: 10
          23: 11
          24: 15
          25: 11
          26:
          #Read 25 items
 
         porcentaje.algodon = c(15,15,15,15,15,20,20,20,20,20,25,25,25,25,25,30,30,30,30,30,35,35,35,35,35)
         datos.algodon = cbind(algodon,porcentaje.algodon)
         datos.algodon
        #      algodon porcentaje.algodon
        # [1,]       7                 15
        # [2,]       7                 15
        # [3,]      15                 15
        # [4,]      11                 15
        # [5,]       9                 15
        # [6,]      12                 20
        # [7,]      17                 20
        # [8,]      12                 20
        # [9,]      18                 20
        #[10,]      18                 20
        #[11,]      14                 25
        #[12,]      18                 25
        #[13,]      18                 25
        #[14,]      19                 25
        #[15,]      19                 25
        #[16,]      19                 30
        #[17,]      25                 30
        #[18,]      22                 30
        #[19,]      19                 30
        #[20,]      23                 30
        #[21,]       7                 35
        #[22,]      10                 35
        #[23,]      11                 35
        #[24,]      15                 35
        #[25,]      11                 35
 
 
        #Hasta aqu? tenemos los datos de resistencia de las 25 camisas, junto con la proporci?n de algod?n que hay en cada una.
        #Ahora aplicaremos el test de igualdad de medias. En la funci?n se especifica como par?metro  algodon~porcentaje.algodon  lo que significa que estudiaremos la variable algodon (su resistencia), dividiendola en grupos seg?n el porcentaje de algod?n presente.
 
       oneway.test(algodon~porcentaje.algodon,data=datos.algodon)
        #                 One-way analysis of means (not assuming equal variances)
        #        data:  algodon and porcentaje.algodon
        #F = 12.4507, num df = 4.000, denom df = 9.916, p-value = 0.0006987
       ##     Como vemos, el p-valor es pr?cticamente cero. Se rechaza, pues, la hip?tesis de que todos los grupos son similares. Esto es, variar el porcentaje de algod?n tiene consecuencias en la resistencia de la camisa. De hecho, podemos hacer un diagrama de caja (un boxplot) por grupos, para ver c?mo evoluciona la resistencia en funci?n del porcentaje de algod?n:
 
        boxplot(algodon~porcentaje.algodon)
 
 
#### var.test ##################################################################
# Este es un contraste para comprobar la igualdad de varianzas entre dos grupos de datos.
 
      # Por ejemplo, cuando tenemos dos m?quinas que fabrian lo mismo (por ejemplo, rodamientos), y queremos saber si lo hacen con la misma variabilidad. Para poner un ejemplo, generemos datos de dos distribuciones normales con distinta varianza y veamos si el test detecta esta situaci?n:
       x = rnorm(50, mean = 0, sd = 2)
       y = rnorm(30, mean = 1, sd = 1)
       var.test(x, y)
      #               F test to compare two variances
      # data:  x and y
      # F = 5.3488, num df = 49, denom df = 29, p-value = 7.546e-06
      # alternative hypothesis: true ratio of variances is not equal to 1
      # 95 percent confidence interval:
      #   2.687373 10.063368
      # sample estimates:
      # ratio of variances
      #           5.348825
     ## Como vemos, el p-valor sale pr?cticamente cero: rechazaremos la hip?tesis de que los conjuntos de datos x e y provienen de distribuciones con la misma varianza.
 

No hay comentarios:

Publicar un comentario

Libros para descargar (gratis) sobre Estadística Computacional