Diversidad de especies versus grupo de nivel de zinc

Medley y Clements (1998) muestrearon varias estaciones en seis arroyos que se sabía que estaban contaminados por metales pesados en la región de las Montañas Rocosas de Colorado, Estados Unidos. Registraron la concentración de zinc y la riqueza de especies y diversidad de especies de la comunidad de diatomeas y proporción de células de diatomeas que fueron la especie de sucesión temprana, Achanthes minutissima. Los datos se encuentran en la hoja “Datos” del archivo EJEMPLO_DCA_est.xlsx

El primer análisis compara la diversidad media de especies de diatomeas (variable de respuesta) en los cuatro grupos de niveles de zinc (variable predictiva categórica : BACK, HIGH, LOW y MED), nivel de zinc tratado como factor fijo.

La H0 fue que no había diferencia en la diversidad media de especies de diatomeas entre grupos de niveles de zinc.

Tabla ANOVA

Para probar las hipótesis dada utilizaremoas ANOVA para un diseño completo al azar (DCA), en el cual descompone la variabilidad total de los datos en sus dos componentes: la variabilidad debida a tratamientos y la que corresponde al error aleatorio.

library(readxl)
medley <- read_excel("EJEMPLO_DCA_est.xlsx", 
    sheet = "Datos") # leer los datos
## New names:
## * `` -> ...5
medley
# obtencion de la tabla ANOVA
fit <- aov(DIVERSTY~ZINC, data= medley)
summary(fit)
##             Df Sum Sq Mean Sq F value Pr(>F)  
## ZINC         3  2.567  0.8555   3.939 0.0176 *
## Residuals   30  6.516  0.2172                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

La validez de los resultados obtenidos queda supeditado a que los supuestos del modelo se cumplan, a continuación realizamos la verificación de los supuestos

Verificación de supuestos

Supuesto de normalidad de residuos

Para verificar la normalidad de los residuos utilizaremos la prueba de Shapiro-Wilks cuyo script es el siguiente:

shapiro.test(residuals(fit))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(fit)
## W = 0.96883, p-value = 0.43

Con el resultado anterior con una confianza del 95 % no podemos rechazar la idea de que los residuos provengan de una distribución normal \((W=0.9688, p-value > 0.05)\)

Homocedasticidad de Varianzas

Para verificar el supuesto de homocedasticidad de las varianzas utilizaremos la prueba de Bartlett script es el siguiente:

bartlett.test (DIVERSTY~ZINC, data= medley)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  DIVERSTY by ZINC
## Bartlett's K-squared = 0.25294, df = 3, p-value = 0.9686

Con una confianza del 95 % no podemos rechazar la idea de que las varianzas sean homogeneas \((chi^2= 4.0967, df = 4, p-value>0.05)\)

Independencia

plot(fit$residuals)

El grafico no muestra evidencia de falta de independencia.

Dado que se han verificado los supuestos del modelo ahora podriamos concluir que con una confianza del 95% rechamos la hipótesis nula \((F= 3.939,df(3,30),p-value < 0.05)\), es decir, la diversidad media de especies de diatomeas en por lo menos un grupo de niveles de zinc es estadisticamente diferente a los otros grupo.

Para averiguar cuál de los grupos es el que difiere procedemos a realizar un test de rangos multiples.

Test de rangos multiples

Test LSD

library(agricolae)
 LSD.test(fit, "ZINC", console = TRUE, group= TRUE, alpha = 0.05)
## 
## Study: fit ~ "ZINC"
## 
## LSD t Test for DIVERSTY 
## 
## Mean Square Error:  0.2172137 
## 
## ZINC,  means and individual ( 95 %) CI
## 
##      DIVERSTY       std r       LCL      UCL  Min  Max
## BACK 1.797500 0.4852613 8 1.4609789 2.134021 0.76 2.27
## HIGH 1.277778 0.4268717 9 0.9605026 1.595053 0.63 1.90
## LOW  2.032500 0.4449960 8 1.6959789 2.369021 1.40 2.83
## MED  1.717778 0.5030104 9 1.4005026 2.035053 0.80 2.19
## 
## Alpha: 0.05 ; DF Error: 30
## Critical Value of t: 2.042272 
## 
## Groups according to probability of means differences and alpha level( 0.05 )
## 
## Treatments with the same letter are not significantly different.
## 
##      DIVERSTY groups
## LOW  2.032500      a
## BACK 1.797500      a
## MED  1.717778     ab
## HIGH 1.277778      b
 LSD.test(fit, "ZINC", console = TRUE, group= FALSE, alpha = 0.05)
## 
## Study: fit ~ "ZINC"
## 
## LSD t Test for DIVERSTY 
## 
## Mean Square Error:  0.2172137 
## 
## ZINC,  means and individual ( 95 %) CI
## 
##      DIVERSTY       std r       LCL      UCL  Min  Max
## BACK 1.797500 0.4852613 8 1.4609789 2.134021 0.76 2.27
## HIGH 1.277778 0.4268717 9 0.9605026 1.595053 0.63 1.90
## LOW  2.032500 0.4449960 8 1.6959789 2.369021 1.40 2.83
## MED  1.717778 0.5030104 9 1.4005026 2.035053 0.80 2.19
## 
## Alpha: 0.05 ; DF Error: 30
## Critical Value of t: 2.042272 
## 
## Comparison between treatments means
## 
##              difference pvalue signif.         LCL          UCL
## BACK - HIGH  0.51972222 0.0289       *  0.05721821  0.982226233
## BACK - LOW  -0.23500000 0.3213         -0.71091270  0.240912701
## BACK - MED   0.07972222 0.7273         -0.38278179  0.542226233
## HIGH - LOW  -0.75472222 0.0023      ** -1.21722623 -0.292218212
## HIGH - MED  -0.44000000 0.0543       . -0.88869480  0.008694797
## LOW - MED    0.31472222 0.1748         -0.14778179  0.777226233
library(ggplot2)
ggplot(data = medley,mapping=aes(x = DIVERSTY,y = ZINC)) +
  geom_jitter(alpha=0.25,width=0.01) +
  stat_summary(fun.data= "mean_cl_normal",
               geom ="pointrange",
               size = 1.1,
               fatten = 3,
               pch = 21,
               fill = "red",
               colour = "black")+
  coord_flip()

Test Tuckey

library(agricolae)
HSD.test(fit, "ZINC", console = TRUE, group= TRUE, alpha = 0.05)
## 
## Study: fit ~ "ZINC"
## 
## HSD Test for DIVERSTY 
## 
## Mean Square Error:  0.2172137 
## 
## ZINC,  means
## 
##      DIVERSTY       std r  Min  Max
## BACK 1.797500 0.4852613 8 0.76 2.27
## HIGH 1.277778 0.4268717 9 0.63 1.90
## LOW  2.032500 0.4449960 8 1.40 2.83
## MED  1.717778 0.5030104 9 0.80 2.19
## 
## Alpha: 0.05 ; DF Error: 30 
## Critical Value of Studentized Range: 3.845401 
## 
## Groups according to probability of means differences and alpha level( 0.05 )
## 
## Treatments with the same letter are not significantly different.
## 
##      DIVERSTY groups
## LOW  2.032500      a
## BACK 1.797500     ab
## MED  1.717778     ab
## HIGH 1.277778      b
HSD.test(fit, "ZINC", console = TRUE, group= FALSE, alpha = 0.05)
## 
## Study: fit ~ "ZINC"
## 
## HSD Test for DIVERSTY 
## 
## Mean Square Error:  0.2172137 
## 
## ZINC,  means
## 
##      DIVERSTY       std r  Min  Max
## BACK 1.797500 0.4852613 8 0.76 2.27
## HIGH 1.277778 0.4268717 9 0.63 1.90
## LOW  2.032500 0.4449960 8 1.40 2.83
## MED  1.717778 0.5030104 9 0.80 2.19
## 
## Alpha: 0.05 ; DF Error: 30 
## Critical Value of Studentized Range: 3.845401 
## 
## Comparison between treatments means
## 
##              difference pvalue signif.         LCL        UCL
## BACK - HIGH  0.51972222 0.1219         -0.09606192  1.1355064
## BACK - LOW  -0.23500000 0.7457         -0.86863665  0.3986367
## BACK - MED   0.07972222 0.9847         -0.53606192  0.6955064
## HIGH - LOW  -0.75472222 0.0117       * -1.37050636 -0.1389381
## HIGH - MED  -0.44000000 0.2096         -1.03739837  0.1573984
## LOW - MED    0.31472222 0.5153         -0.30106192  0.9305064