Manual rápido de métodos estadísticos

Nunca comiences con pruebas estadísticas, primero analiza los datos con una figura.

Primero: Imagina el patrón esperado de tus resultados.

Segundo: Diagrama la figura que debería cumplir con este patrón.

Tercero: Identificas el mejor modelo que explique las diferencias esperadas.

Cuarto: Evalúa los supuestos del modelo.

Quinto: Evalúa el modelo e interpreta los resultados.

Sexto: Integrar los resultados del modelo en la figura original.

Objetivos

  1. Chi-cuadrado y tablas de contingencias

  2. Prueba t de dos muestras

  3. ANOVA

  4. Regresión linear

Paquetes

library(dplyr)
library(ggfortify)
library(ggplot2)

Datos a emplear

1. Chi-cuadrado y trablas de contingencias

Esta prueba analiza tablas de contingencia de datos de cuentas. Evaluando la asociaión entre dos o más variables categóricas.

El ejercicio evalúa la asociación del color de aves “Lady birds” (factor 2 niveles: black y red) y su hábitat (factor 2 niveles: Rural e Industrial).

Biológicamente, tratamos de responder la pregunta si alguna característica del hábitat está asociado con las frecuencias de diferentes colores de las aves.

Buscamos evaluar si las frecuencias de colores de ladybirds difieren entre los tipos de hábitat.

Exploramos la base de datos:

str(lady)
## 'data.frame':    20 obs. of  4 variables:
##  $ Habitat     : Factor w/ 2 levels "Industrial","Rural": 2 2 2 2 2 2 2 2 2 2 ...
##  $ Site        : Factor w/ 10 levels "R1","R2","R3",..: 1 2 3 4 5 1 2 3 4 5 ...
##  $ morph_colour: Factor w/ 2 levels "black","red": 1 1 1 1 1 2 2 2 2 2 ...
##  $ number      : int  10 3 4 7 6 15 18 9 12 16 ...

Calculamos totales:

totals <- lady %>% group_by(Habitat, morph_colour) %>% summarise(total.number = sum(number))
totals
## # A tibble: 4 x 3
## # Groups:   Habitat [2]
##   Habitat    morph_colour total.number
##   <fct>      <fct>               <int>
## 1 Industrial black                 115
## 2 Industrial red                    85
## 3 Rural      black                  30
## 4 Rural      red                    70

Creamos un diagrama de barras que ejemplifique los conteos de color por zona, adicionalmente le añadimos colores que ayuden a ejemplificar mejor los datos.

ggplot(totals, aes(x = Habitat, y = total.number, fill = morph_colour)) + geom_bar(stat = "identity", 
    position = "dodge") + scale_fill_manual(values = c(black = "black", red = "red"))

La hipótesis nula es la NO ASOCIACIÓN entre la frecuencia de colores en cada hábitat. Recordar que las pruebas Chi-cuadrado comparar PROPORCIONES y no números absolutos.

Empleamos la función chisq.test(), es necesario tabular los datos en una matriz de doble entrada con la funcion xtabs().

lady.mat <- xtabs(number ~ Habitat + morph_colour, data = lady)
lady.mat
##             morph_colour
## Habitat      black red
##   Industrial   115  85
##   Rural         30  70

Prueba de chi-cuadrado:

lady.chi <- chisq.test(lady.mat)

El valor de \(p\) es muy pequeño, lo que indica una muy baja probabilidad de que lo observado suceda por casualidad. El resultado permite RECHAZAR la hipótesis nula y concluir que hay alguna asociación. La figura sugiere que los colores negros son más frecentes en los hábitats industriales.

Se expresa:

“Evaluamos la hipótesis que hay asociación del color de las aves y los hábitats industrial y rual, Los colores no están distribuidos (Chi2 = 19.1, df = 1, p<0.001), con las formas negras siendo más frecuentes en hábitats industriales y las formas rojas en hábitats rurales.”

Los valores esperados son calculados mediante:

names(lady.chi)
## [1] "statistic" "parameter" "p.value"   "method"    "data.name" "observed" 
## [7] "expected"  "residuals" "stdres"
lady.chi$expected
##             morph_colour
## Habitat         black       red
##   Industrial 96.66667 103.33333
##   Rural      48.33333  51.66667

2. Prueba t de dos muestras

Comparamos dos grupos de valores numéricos. Es apropiado cuando los tamaños de muestra en cada grupo son pequeños.

El modelo asume que los datos en cada grupo están normalmente distribuidos y que sus varianzas son iguales.

Caso:

Se analizan los niveles de ozono en jardines distribuidos en la ciudad. Estuvieron al oeste del centro o al este, el ozono en exceso puede dañar plantas de lechuga (> 8pphm).

La hipótesis nula es que no hay diferencia en los niveles medios de ozono entre grupos.

str(ozone)
## 'data.frame':    20 obs. of  3 variables:
##  $ Ozone          : num  61.7 64 72.4 56.8 52.4 44.8 70.4 67.6 68.8 53.7 ...
##  $ Garden.location: Factor w/ 2 levels "East","West": 2 2 2 2 2 2 2 2 2 2 ...
##  $ Garden.ID      : Factor w/ 20 levels "G1","G10","G11",..: 1 12 14 15 16 17 18 19 20 2 ...

Paso 1: El gráfico (Histogramas)

Consideramos un método que permite ver la tendencia central de los datos y su variabilidad, como sólo hay dos grupos se opta por un histograma.

ggplot(ozone, aes(x = Ozone)) + geom_histogram(binwidth = 10) + facet_wrap(~Garden.location, 
    ncol = 1) + theme_bw()

En base a la imspección visual es razonable pensar que el supuesto de normalidad e igualdad de varianza son cumplidos. Evaluamos los estadísticos:

ozone_sum <- ozone %>% group_by(Garden.location) %>% summarise(prom.ozone = mean(Ozone), 
    sd.ozone = sd(Ozone))

Paso 2: Análisis de dos muestras - Prueba-t

Ejecutamos la prueba-t:

t.test(Ozone ~ Garden.location, data = ozone)
## 
##  Welch Two Sample t-test
## 
## data:  Ozone by Garden.location
## t = 4.2363, df = 17.656, p-value = 0.0005159
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##   8.094171 24.065829
## sample estimates:
## mean in group East mean in group West 
##              77.34              61.26

Interpretación:

Primera línea: Welch Two Sample t-test Es la prueba efectuada. Welch quiere decir este método permite que los supuestos de homogeneidad de varianzas.

Segunda línea: Los datos empleados.

Tercera línea: Estadístico t, grados de libertad df y valor de p p-value

Cuarta línea: Explica la hipótesis alternativa

Quita línea: Muestra la diferencia entre las medias, si el intervalo incluye el cero entonces podemos interpretarlo como que no hay diferencias entre las medias de los grupos.

Ya que se usó la prueba de Welch, verificamos la homogeneidad de varianzas

var.test(Ozone ~ Garden.location, data = ozone)
## 
##  F test to compare two variances
## 
## data:  Ozone by Garden.location
## F = 0.75503, num df = 9, denom df = 9, p-value = 0.6823
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.1875386 3.0397437
## sample estimates:
## ratio of variances 
##          0.7550293

3. ANOVA de una vía

Caso:

Se desea saber si la tasa de crecimiento de Daphnia es alterada por la infección de parásitos Pregunta 1. También podemos comparar si los efectos de los parásitos difieren de un control Pregunta 2. Usaremos las variables growth.rate y parasite

str(daphnia)
## 'data.frame':    40 obs. of  3 variables:
##  $ parasite   : Factor w/ 4 levels "control","Metschnikowia bicuspidata",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ rep        : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ growth.rate: num  1.07 1.27 1.32 1.08 1.2 ...

Diagrama

Creamos un diagrama de cajas el cual es ideal para cuando se comparan variables categóricas como predictores y variables numéricas como variable de respuesta.

ggplot(daphnia, aes(x = parasite, y = growth.rate)) + geom_boxplot() + theme_bw() + 
    coord_flip()

El diagrama sugiere diferencias en algunos de los tratamientos de parásitos, con los que podemos suponer que al menos unos de los tratamientos tienen efectos (Pregunta 1) y ordenar visualmente cuáles tratamientos tienen mayor tasa de crecimiento (Pregunta 2).

Construyendo el modelo

En R el modelo empleado es modelo de regresión linear. La variable explicatoria debe ser un factor o variable categórica.

model_grow <- lm(growth.rate ~ parasite, data = daphnia)

Revisamos los supuestos

Para revisar los supuestos empleamos la función autoplot

autoplot(model_grow, smooth.colour = NA) + theme_bw()

Todos los parámetros están bien, sin patrones y sin valore extremos, por otro lado el diagrama Q-Q puede sugerir que no hay normalidad, sin embargo el patrón permanece aún dentro de lo esperado. Para más detalles revisar la sección de supuestos de la regresión linear simple.

Modelo

Usamos la función anova()

anova(model_grow)
## Analysis of Variance Table
## 
## Response: growth.rate
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## parasite   3 3.1379 1.04597  32.325 2.571e-10 ***
## Residuals 36 1.1649 0.03236                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El resultado sugiere que hay evidencia de que el tratamiento de parásito ha producido un efecto.

La hipótesis nula es “todos los grupos de poseen la misma media”, los valores de F cuantifican que tan probable es que esto sea cierto. Un valor alto de F permite obtener un valor chico de \(p\) que a su vez permite rechazar la hipótesis nula de no diferencias.

Contraste de tratamientos

R emplea contrastes para resumir las diferencias entre tratamientos. Es necesario saber el método empleado para calcular los contrastes, y sus resultados pueden variar entre programas.

Para saber qué método estás empleado revisar el comando ?contr.treatment.

Para hacer contraste empleamos:

summary(model_grow)
## 
## Call:
## lm(formula = growth.rate ~ parasite, data = daphnia)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.41930 -0.09696  0.01408  0.12267  0.31790 
## 
## Coefficients:
##                                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                        1.21391    0.05688  21.340  < 2e-16 ***
## parasiteMetschnikowia bicuspidata -0.41275    0.08045  -5.131 1.01e-05 ***
## parasitePansporella perplexa      -0.13755    0.08045  -1.710   0.0959 .  
## parasitePasteuria ramosa          -0.73171    0.08045  -9.096 7.34e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1799 on 36 degrees of freedom
## Multiple R-squared:  0.7293, Adjusted R-squared:  0.7067 
## F-statistic: 32.33 on 3 and 36 DF,  p-value: 2.571e-10

Interpretación:

  1. (Intercept) este corresponde al grupo llamado control

  2. Los otros nombres corresponden a los niveles de tratamientos de parásitos.

  3. El contraste se hace en orden alfabético es decir el orden de presentación es \(control\) < \(M.bicuspidata\) < \(P.perplexa\) < \(P.ramosa\). y es orden presentado en el resultado de summary donde control corresponde a Intercept.

  4. Los contrastes reportan \(diferencias\) entre los niveles respecto del nivel de referencia, en nuestra caso el control, y los otros niveles. Tal que los números asociados a cada parásito son las diferencias entre las tasas de crecimiento de cada parásito y el control, estos es porqué todos son negativos.

  5. Los valores de p están asociados al grupo control, por lo tanto sirven para determinar si las diferencias entre los niveles y control son significativas.

4. Regresión linear simple

Caso:

Comparamos el crecimiento de la planta respecto de la humedad del suelo. La predicción subyacente es que a más humedad permitirá mayores tasas de crecimiento.

Las variables tienen dos características:

  1. Hay una clara relación entre las dos variables.

  2. La variable explicatoria es continua, variable numérica que no tiene categorías.

str(plant_gr)
## 'data.frame':    50 obs. of  2 variables:
##  $ soil.moisture.content: num  0.47 0.541 1.698 0.826 0.857 ...
##  $ plant.growth.rate    : num  21.3 27 39 30.2 37.1 ...

Gráfico exploratorio

Realizamos un diagrama de puntos:

ggplot(plant_gr, aes(x = soil.moisture.content, y = plant.growth.rate)) + geom_point() + 
    ylab("Plant Growth Rate (mm/week)") + theme_bw()

El gradiente (pendiente) es positivo, de hecho cuanto más humedad en el sueño mayor es la tasa de crecimiento.

Modelamos una regresión linear simple

Ajustamos el modelo, se interpreta como hipótesis que la tasa de crecimiento es una función de la humedad el suelo.

model_pgr <- lm(plant.growth.rate ~ soil.moisture.content, data = plant_gr)

Evaluamos los supuestos

Empleamos la función autoplot del paquete ggfortify, el cual produce lo cuatro plots para evaluar la normalidad basados en residuos y en su relación de suma de cuadrados y promedio de errores cuadrados.

autoplot(model_pgr, smooth.colour = NA) + theme_bw()

Interpretamos los gráficos:

  1. Arriba izquierda: Nos dice si los datos están apropiadamente ajustados a línea. Si los puntos toman forma de curva o valles las cosas no van bien. Este gráfico se relaciona a la igualdad de varianzas homo vs heterocedasticidad.

  2. Arriba a la derecha: La curva Normal Q-Q evalúa el supuesto de normalidad de los residuos, los puntos son los residuos y la línea punteada la distribución esperada bajo una distribución normal. Esta forma es más confiable que otros métodos incluso cuando la muestra es menor de 100 datos.

  3. Abajo a la izquierda: Evalúa el supuesto de igualdad de variaza. El eje y es un indicador estandarizado de la variación. El modelo linear asume que la varianza es constante a través de todos los valores predichos de la variable de respuesta. No debe haber patrón.

  4. Abajo a la derecha: Este evalúa la “influencia”, esta herramienta no solo detecta puntos influentes, algunos que mueven la tendencia más de lo que es esperado, y también detecta outliers.

Los supuestos de normalidad y homogeneidad se cumplen, en los gráficos de la izquierda no hay patrón, en el gráfico Q-Q se cumple con el ajuste a la línea, en el gráfico derecha abajo no hay un punto que destaque.

Modelo

Hipótesis nula: La humedad del suelo no tiene efectos en la tasa de crecimiento de las plantas.

Extraemos los resultados del modelo:

anova(model_pgr)
## Analysis of Variance Table
## 
## Response: plant.growth.rate
##                       Df  Sum Sq Mean Sq F value    Pr(>F)    
## soil.moisture.content  1 2521.15 2521.15  156.08 < 2.2e-16 ***
## Residuals             48  775.35   16.15                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova muestra los grados de libertad, suma de cudarados, y el valor de F, este último tiene un valor alto que indica que el error e la varianza es pequeño relativo a la varianaza atribuida a la variable explicatoria

summary(model_pgr)
## 
## Call:
## lm(formula = plant.growth.rate ~ soil.moisture.content, data = plant_gr)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.9089 -3.0747  0.2261  2.6567  8.9406 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             19.348      1.283   15.08   <2e-16 ***
## soil.moisture.content   12.750      1.021   12.49   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.019 on 48 degrees of freedom
## Multiple R-squared:  0.7648, Adjusted R-squared:  0.7599 
## F-statistic: 156.1 on 1 and 48 DF,  p-value: < 2.2e-16

El modelo indica los valores de los parámetros el intercepto y la pendiente.

Los resultados se interpretan como sigue:

“La humedad del suelo tiene un efecto positivo en el crecimiento de la planta. Por cada unidad de incremento de la humedad, la tasa de creciiento incrementa en 12.7 mm/week (pendiente = 12.7, t = 12.5, d.f. = 48, \(p\) < 0.001)”.

Complementamos la figura

Añadimos una línea de regresión:

ggplot(plant_gr, aes(x = soil.moisture.content, y = plant.growth.rate)) + geom_point() + 
    geom_smooth(method = "lm") + ylab("Plant Growth Rate (mm/week)") + theme_bw()

geom_smooth(method="lm") añade una línea ajustada, y el eror estándar del ajuste.

6. Notas

  1. Tener una ruta de análisis consistente.

  2. Siempre hacer un gráfico que responda a tu pregunta antes de cualquier prueba.

  3. Interpretar tanto como puedas, mirar pendientes, interceptos, contrastes, grados de libertad, etc desde el gráfico antes de la estadística.

  4. Revisar si los supuestos de los datos no se cumplen de alguna manera.

  5. Hacer un buen gráfico para comunicar los resultados.

  6. Crear conclusiones biológica y no estadísticamente centrados para describir los resultados.

  7. Usa R!.

Fernando Tapia Vilchez . 04 de Junio del 2019. Última actualización: 2019-06-05