Estudio Mortalidad Patos Merganneta armata

Introducción

El pato Merganneta armata o pato contra corriente es una especie que vive en los Andes por todo su extensión en Sudamérica, desde los 4500 metros de altura hasta el nivel del mar y como su nombre común lo advierte gusta de ríos con corriente fuertes. Esta especie presenta dimorfismo sexual, los machos tienen manchas blancas y las hembras manchas anaranjadas, esto se puede apreciar en la figura 1, el pato de la derecha corresponde a un especimen juvenil. Estos patos se encuentran actualmente en la lista de animales de menor preocupación dentro de “IUCN Red List of Threatened Species”, a pesar de esto, sus números han ido en declinación con el pasar de los años y cuando la mortalidad alcance 30% serán movidos a la categoría de vulnerables (BirdLife International, 2016).

Figura 1: Patos Merganneta armata. Fuente: Duck Duck Go images

Metodología

En primer lugar, para comenzar el estudio se realiza una matriz de confusión, donde se muestran las hipotesis y predicciones para este estudio y los casos en los que se puede cometer error.

Tabla 1: Matriz de confusión para presente estudio. En la primera fila se muestran las hipótesis y en la primera columna las predicciones.

.	No hay diferencias en las tasas de mortalidad entre ríos	Al menos una de las tasas de mortalidad es diferente
No hay diferencias en las tasas de mortalidad entre ríos	No hay error	Error tipo II (ß)
Al menos una de las tasas de mortalidad es diferente	Error tipo I (a)	No hay error

Para el presente estudio la hipotesis nula correspondería a que no existe diferencias en mortalidad entre ríos, ecuación 1. Por otra parte, la hipótesis alternativa es que al menos una de las mortalidades es diferente, ecuación 2.

Ecuación 1: \(Ho: μ1 = μ2 = μ3\)

Ecuación 2: \(Ha: μ1 ≠ μ2 = μ3\) ó \(Ha: μ1 = μ2 ≠ μ3\)

Para este estudio beta representa la probabilidad de no encontrar diferencias de mortalidad de patos entre ríos cuando en realidad existe una diferencia (aceptar Ho cuando es falsa), o sea, los patos están muriendo por culpa de la contaminación y el estudio determina que no mueren. Por otra parte, alfa representa la probabilidad de encontrar una diferencia de mortalidad de patos entre ríos cuando en realidad no la hay (rechazar Ho cuando es verdadera), o sea, los patos no mueren por la contaminación y el estudio determina que en realidad si mueren.
En este estudio se considera más importante que se encuentre si los patos mueren de verdad, por esto, el error tipo II o beta debería ser el menor posible y es que mantendrá constante, ya que no encontrar si los patos están muriendo podría conllevar, en el peor de los casos, a la extinción de la especie en este río.

Para realizar los análisis es necesario establecer el directorio de trabajo con la función setwd. En mi caso es de la siguiente manera.

setwd("G:/Mi unidad/5to/Bioes 2/Practico")

Luego se cargan las librerías necesarias para trabajar, tidyverse para hacer la visualización y ordenar los datos fácilmente y pwr2 para poder realizar los análisis de poder necesarios.

library(tidyverse)
library(pwr2)

Se leen los archivos necesarios, en este caso Patos.csv corresponde al conjunto completo de los datos a usar después de determinar n y MuestraPatos.csv a un conjunto menor de observaciones que se usará para las pruebas de poder a priori.

patos <- readr :: read_csv("Tarea/Patos.csv")
muestra_patos <- readr :: read_csv("Tarea/MuestraPatos.csv")

Resultados

Prueba de poder a priori

Para poder realizar un prueba de poder se necesita establecer cinco parámetros.
1. Los grupos pre establecidos son rio Contaminado, No Contaminado 1 y No Contaminado 2, por lo que k es tres (k = 3).
2. El alfa aceptado mundialmente para estudios científicos biológicos es 0.05 (alpha = 0.05).
3. Para calcular el sigma, se utilizan los datos obtenidos en el primer estudio (MuestraPatos) de la siguiente manera:
3.1. Se realiza un ANOVA a los datos.
3.2. Del resultado se obtienen los residuales.
3.3. Se calcula la desviación estándar de los resultados y se almacena para uso posterior.

#lo anterior en un solo comando
DE <- sd(broom::augment(aov(mortalidad ~ rio, data = muestra_patos))$.resid)

## [1] "Valor de sigma a usar es" "7.70880277430675"

4. Para determinar delta, se observa como se comportan los datos de muestreo previo.

Tabla 2: Medias de mortalidad inicial en patos en los ríos muestreados

Rio	Media_Mortalidad
Contaminado	45.85630
No Contaminado 1	38.23249
No Contaminado 2	40.78101

En la tabla se observa una diferencia entre contaminado y no contaminado entre 5 y 7 unidades, por esto se determina un delta de 6 para la primera exploración de poder.

5. El n inicial a probar es 42 que es la respuesta a la pregunta de la vida, el universo y todo (Adams, 1979).
   Finalmente queda por establecer un poder mínimo que se quiere obtener, se establece un valor de 0,95.

#prueba 1 con 42 por grupo
pwr.1way(k = 3, n = 42, alpha = 0.05, delta = 6, sigma = DE)

## 
##      Balanced one-way analysis of variance power calculation 
## 
##               k = 3
##               n = 42
##           delta = 6
##           sigma = 7.708803
##     effect.size = 0.3177523
##       sig.level = 0.05
##           power = 0.8940889
## 
## NOTE: n is number in each group, total sample = 126 power = 0.894088947677986

Se obtiene un poder de 0,89 que esta aún bajo el límite que nos impusimos. Se aumenta el n a 52 para poder llegar al poder establecido.

#prueba 2 con 52 por grupo
pwr.1way(k = 3, n = 52, alpha = 0.05, delta = 6, sigma = DE)

## 
##      Balanced one-way analysis of variance power calculation 
## 
##               k = 3
##               n = 52
##           delta = 6
##           sigma = 7.708803
##     effect.size = 0.3177523
##       sig.level = 0.05
##           power = 0.9500056
## 
## NOTE: n is number in each group, total sample = 156 power = 0.950005632087449

Con esto se determina que el n necesario es de 52 unidades experimentales por tratamiento para tener un poder de 0,95, detectar diferencias de 6 unidades entre medias y con un alfa de 0,05.

ANOVA

Con esta información se procede a realizar un muestreo aleatorizado del archivo completo disponible. Con el muestreo al azar se procede a observar los datos , determinar si cumplen los supuestos del ANOVA y realizar el ANOVA si estos se cumplen.

El primer supuesto es la independiencia de datos, no hay modo de verificar esto ahora, pero se asume que los datos fueron tomados de zonas independientes entre sí.
El segundo supuesto a verificar es la aleatorización para lograr esto se realiza el muestreo de 52 unidades experimentales de forma aleatoria usando la función sample_n de R. Para que este muestreo sea reproducible, se usa la función set.seed.

set.seed(13)
patos1 <- patos %>% group_by(rio) %>% sample_n(52)

El tercer supuesto a verificar es la homocedasticidad de los datos. Para esto se usan dos métodos, el primero es el visual donde se espera observar que las cajas y bigotes del boxplot sean parecidos. El segundo es el test de Bartlett, el cual provee un valor de p que si es menor a 0,05 las varianzas son distintas.

bartlett.test(data = patos1, mortalidad ~ rio)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  mortalidad by rio
## Bartlett's K-squared = 2.1628, df = 2, p-value = 0.3391

Visualmente las varianzas se parecen, en conjunto con el valor de p demuestran que la varianza entre grupos es la misma, esto significa que nuestros datos son homocedásticos.

Finalmente se verifica si los datos son normales, para esto se verifica la normalidad de los residuales. Primero se consiguen los residuales de los datos.

Residuales <- aov(mortalidad ~ rio, data = patos1)$residuals

Luego se puede verificar visualmente usando la función hist que genera un histograma y con qqnorm que genera un gráfico de línea de tendencia que si es linear representa normalidad.

hist(Residuales, col = 'green')

Figura 4: Histograma de residuales de los datos muestreados.

qqnorm(Residuales, col = 'red')

Figura 5: grafico de probabilidades de residuales de los datos muestreados.

Visualmente los residuales parecen normales, la última manera y más precisa es realizar un test de normalidad. En este caso se usa el test Shapiro-Wilk, en este test si el valor de p es mayor a 0,05 los datos son normales.

shapiro.test(Residuales)

## 
##  Shapiro-Wilk normality test
## 
## data:  Residuales
## W = 0.99317, p-value = 0.6722

Este resultado demuestra que nuestros datos son normales. Con toda esta información de que los supestos se cumplen se procede a realizar el ANOVA; cabe notar que ANOVA es un test muy robusto, por lo que si algunos de los últimos dos supuestos no se cumplen aún es posible realizar el test.

summary(aov(mortalidad ~ rio, data = patos1))

##              Df Sum Sq Mean Sq F value   Pr(>F)    
## rio           2   5505  2752.3   27.97 4.43e-11 ***
## Residuals   153  15054    98.4                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Los resultados muestran que las medias son significativamente distintas, al observar las medias de los datos en la tabla 3 se puede decir que la mortalidad de los patos en el rio contaminado es significativamente distinta y mayor que la de al menos uno de los ríos no contaminados.

Tabla 3: Medias de mortalidad en patos

Rio	Media_Mortalidad
Contaminado	51.86782
No Contaminado 1	38.74658
No Contaminado 2	39.86093

Prueba de poder a priori v2

La ONG puso una limitación, el presupuesto disponible para el proyecto es de $20.000.000. El costo por zona de río es de $500.000, esto significa $1.500.00 por cada n (tres ríos), esto límita a 13 unidades experimentales máximas. Usando el valor nuevo de n con los parámetro previamente establecidos.

#prueba 3 con 13 por grupo
pwr.1way(k = 3, n = 13, alpha = 0.05, delta = 6, sigma = DE)

## 
##      Balanced one-way analysis of variance power calculation 
## 
##               k = 3
##               n = 13
##           delta = 6
##           sigma = 7.708803
##     effect.size = 0.3177523
##       sig.level = 0.05
##           power = 0.3805215
## 
## NOTE: n is number in each group, total sample = 39 power = 0.380521543591419

Se obtiene un poder de 0,38, que es muy bajo para un estudio serio, además de haber previamente establecido que el poder del test debería ser el mayor posible. Con esta información se procede a hacer cambio en los parámetros.
Se opta por aumentar el valor de delta para poder conseguir un poder mayor con menor número de n, se usan datos de organizaciones que determinan la vulnerabilidad de especies (BirdLife International, 2016) y dos estudios de distribución de este pato en Argentina (Aragon y otros, 2011) y Chile (Pernollet y otros, 2013) . Con esto se establece que 12,5 es el nuevo valor de delta para el estudio.

#prueba 4 
pwr.1way(k = 3, n = 13, alpha = 0.05, delta = 12.5, sigma = DE)

## 
##      Balanced one-way analysis of variance power calculation 
## 
##               k = 3
##               n = 13
##           delta = 12.5
##           sigma = 7.708803
##     effect.size = 0.6619839
##       sig.level = 0.05
##           power = 0.9531252
## 
## NOTE: n is number in each group, total sample = 39 power = 0.953125219394817

Con los nuevos parámetros establecidos se alcanza y supera el poder esperado. Resumiendo nuestro estudio con un alfa de 0,05 y un poder de 0,96 es capaz de detectar diferencias de 13 unidades. Con esta información se procede a hacer el estudio nuevamente, pero en esta ocasión solo se obtienen 13 muestras por río.

ANOVA v2

Con los nuevos datos obtenidos se comprueba que cumplan con los mismos supuestos que antes para poder revisar el ANOVA. El primer supuesto se mantiene (independencia), para el segundo se realiza un nuevo muestreo al azar y reproducible de 13 unidades experimentales con la funciónes set.seed y sample_n.

set.seed(13)
patosf <- patos %>% group_by(rio) %>% sample_n(13)

El tercer supuesto a verificar es la homocedasticidad de los datos. Para esto se usan dos métodos, el primero es el visual donde se espera observar que las cajas y bigotes del boxplot sean parecidos. El segundo es el test de Bartlett, el cual provee un valor de p que si es menor a 0,05 las varianzas son distintas.

bartlett.test(data = patosf, mortalidad ~ rio)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  mortalidad by rio
## Bartlett's K-squared = 0.46663, df = 2, p-value = 0.7919

A pesar de que visualmente no lo pareciera, el valor de p de este demuestra que la varianza entre grupos es la misma, esto significa que nuestros datos son homocedásticos.

Finalmente se verifica si los datos son normales, para esto se verifica la normalidad de los residuales. Primero se consiguen los residuales de los datos.

Residualesf <- aov(mortalidad ~ rio, data = patosf)$residuals

Luego se puede verificar visualmente usando la función hist que genera un histograma y con qqnorm que genera un gráfico de línea de tendencia.

hist(Residualesf, col = 'green')

Figura 7: Histograma de residuales de los datos finales.

qqnorm(Residualesf, col = 'red')

Figura 8: grafico de probabilidades de residuales de los datos finales.

Visualmente los residuales parecen normales, la última manera y la más precisa es realizar un test de normalidad, en este caso se usa el test Shapiro-Wilk. En este test si el valor de p es mayor a 0,05 los datos son normales.

shapiro.test(Residualesf)

## 
##  Shapiro-Wilk normality test
## 
## data:  Residualesf
## W = 0.98747, p-value = 0.9353

Este resultado demuestra que nuestros datos son normales. Con toda esta información se procede a realizar el ANOVA.

summary(aov(mortalidad ~ rio, data = patosf))

##             Df Sum Sq Mean Sq F value Pr(>F)  
## rio          2    720   360.2   3.518 0.0402 *
## Residuals   36   3686   102.4                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Discusion

A pesar de tener una limitación de dinero para obtener muestras y tamaños de diferencia que se puede detectar para aumentar el poder del estudio. Con los resultados de este estudio, se puede determinar significativamente con un poder del 96% que los desechos están afectando la mortalidad de los patos Merganneta armata en los ríos estudiados. Como tanto alfa y beta se mantuvieron al 5%, existe una probabilidad de 1 en 20 por error que este estudio este arrojando un resultado errado.
Igual que en el caso anterior, usando la información mostrada en la tabla 4, se puede decir que la mortalidad de los patos en el rio contaminado es significativamente distinta y mayor que al menos una media de mortalidad de los ríos no contaminados. Para determinar a cuál de los ríos es significativamente distinto, es necesario realizar una prueba a-posteriori que afectará directamente alfa o crear un nuevo estudio considerando estas comparaciones desde antes o también llamado estudio de contrastes.

Tabla 4: Medias de mortalidad en patos con datos finales

Rio	Media_Mortalidad
Contaminado	49.34729
No Contaminado 1	42.85715
No Contaminado 2	38.92345

Conclusiones

Usando los datos provistos, se diseño un estudio bioestadístico poderoso que demostró que en el río contaminado la mortalidad de patos Merganneta armata ha aumentado al ser comparado con ríos sin contaminar. Se recomienda a la ONG RioSano a tomar medidas para remediar esta situación.

Bibliografía

Adams, D. (1979). Hitchhiker’s Guide To The Galaxy. London: PAN BOOKS.

Aragón, P.N.S., Rivera, L., Politi, N., 2011. VARIACIÓN DE LA ABUNDANCIA DEL PATO DE TORRENTE (MERGANETTA ARMATA) Y CARACTERISTICAS DEL HÁBITAT EN DOS RÍOS DE MONTAÑA DE LA PROVINCIA DE JUJUY, ARGENTINA. ORNITOLOGIA NEOTROPICAL 22, 589–599.

Baptiste Auguie (2016). gridExtra: Miscellaneous Functions for “Grid” Graphics. R package version 2.2.1. https://CRAN.R-project.org/package=gridExtra

BirdLife International. 2016. Merganetta armata. The IUCN Red List of Threatened Species 2016: e.T22680118A92844656. http://dx.doi.org/10.2305/IUCN.UK.2016-3.RLTS.T22680118A92844656.en. Downloaded on 04 April 2018.

David Robinson (2017). broom: Convert Statistical Analysis Objects into Tidy Data Frames. R package version 0.4.3. https://CRAN.R-project.org/package=broom

González, a L., & Fariña, J. M. (2013). Changes in the Abundance and Distribution of Black-necked Swans (Cygnus melanocoryphus) in the Carlos Andwanter Nature Sanctuary and Adjacent Wetlands, Valdivia, Chile. Waterbirds, 36(4), 507–514. https://doi.org/10.1675/063.036.0408

Hadley Wickham (2017). tidyverse: Easily Install and Load ‘Tidyverse’ Packages. R package version 1.1.1. https://CRAN.R-project.org/package=tidyverse

Hadley Wickham, Jim Hester and Romain Francois (2017). readr: Read Rectangular Text Data. R package version 1.1.1. https://CRAN.R-project.org/package=readr

Lambert, J. (2012). [Pato-Corta-Corriente.jpg].

Pengcheng Lu, Junhao Liu and Devin Koestler (2017). pwr2: Power and Sample Size Analysis for One-way and Two-way ANOVA Models. R package version 1.0. https://CRAN.R-project.org/package=pwr2

Pernollet, C.A., Pavez, E.F., Estades, C.F., 2013. Habitat Selection by Torrent Ducks (Merganetta armata armata) in Central Chile: Conservation Implications of Hydropower Production. Waterbirds 36, 287–299. doi:10.1675/063.036.0306

R Core Team (2017). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/.

RStudio Team (2016). RStudio: Integrated Development for R. RStudio, Inc., Boston, MA URL http://www.rstudio.com/.

Yihui Xie (2017). knitr: A General-Purpose Package for Dynamic Report Generation in R. R package version 1.18.