Pesca y estrategias de Muestreo y Análisis. Sesión 1
Master en Gestión Sostenible de Recursos Pesqueros
Author
Angelo Santana
Published
October 6, 2022
Presentación
A continuación se presenta la actividad realizada durante las sesiones que componen la parte dedicada a la estadística en el curso Pesca y estrategias de Muestreo y Análisis del Master en Gestión Sostenible de Recursos Pesqueros.
Este archivo se ha elaborado utilizando Quarto un sistema de publicación científica y técnica desarrollado por Rstudio y basado en Pandoc que es un conversor universal de documentos. Quarto permite combinar texto, código (en este caso R), resultados del código y gráficos en un único documento formateado para su lectura como pdf, en la web o como documento .docx para su edición en Word.
A lo largo de este documento se usa código R que por defecto no se muestra, aunque puede verse pinchando donde se indica “Mostrar código R”.
Mostrar código R
## Librerías necesarias para el análisis de los datos utilizados en esta presentación## -----------------------------------------------------------------------------------library(tidyverse)library(rstatix)library(qqplotr)library(flextable)library(DescTools)# Algunas funciones de utilidad# ------------------------------## Media y desviación típica:meanSd <-function(x) sprintf("%.2f (%.2f)", mean(x,na.rm=TRUE),sd(x,na.rm=TRUE))## Mediana y cuartiles:medianIQR <-function(x){ q <-quantile(x,probs=c(0.5,0.25,0.75),na.rm=TRUE)sprintf("%.2f (%.2f; %.2f)",q[1],q[2],q[3])}## Formateo del p-valor para que muestre "<0.0001" cuando p es menor que dicho valor:formatPval=function(p) ifelse(p<0.0001,"<0.0001",sprintf("%.4f",p))
Material and methods
Durante las sesiones de este curso veremos un resumen de algunos de los métodos estadísticos más usuales en muchas publicaciones científicas. Toda publicación científica tiene una sección de “Material y métodos” donde se describe qué experimentos u observaciones se han hecho, y de qué manera se han realizado, de tal forma que puedan ser replicadas si se desea. Además de dicha descripción, esta sección debe incluir también una subsección en la que se especifique el tipo de análisis estadístico realizado.
Statistical Analysis
El contenido de esta sección en general será similar al siguiente:
“Continuous variables were described as mean (sd) for normal variables and as median and quartiles for non normal ones. Comparison between two groups were made using the Student’s t-test or the Wilcoxon-Mann-Whitney test depending on whether the variable was normal or not in both groups. When more than two groups were compared, the Fisher’s F-test or the Kruskal Wallis test was used depending on the normality or non-normality of the variables. When significant differences were detected, Tukey post hoc test (for normal variables) or Conover test (for non-normal variables) was performed to identify between which groups differences occurred and to evaluate their magnitude. Normality was tested using the Shapiro-Wilk test. A significance level of 5% was used for all contrasts. General statistical analysis and figures were performed using R software version 4.2.1 (R Core Team, 2022).”
Debe siempre citarse el software utilizado. En el caso de utilizar R, deben citarse también los paquetes específicos que se hayan empleado. Para ello R dispone de la función citation(), que produce la salida siguiente:
Para citar R:
Mostrar código R
citation()
To cite R in publications use:
R Core Team (2022). R: A language and environment for statistical
computing. R Foundation for Statistical Computing, Vienna, Austria.
URL https://www.R-project.org/.
A BibTeX entry for LaTeX users is
@Manual{,
title = {R: A Language and Environment for Statistical Computing},
author = {{R Core Team}},
organization = {R Foundation for Statistical Computing},
address = {Vienna, Austria},
year = {2022},
url = {https://www.R-project.org/},
}
We have invested a lot of time and effort in creating R, please cite it
when using it for data analysis. See also 'citation("pkgname")' for
citing R packages.
Para citar paquetes de R (en este ejemplo tidyverse):
Mostrar código R
citation("tidyverse")
To cite package 'tidyverse' in publications use:
Wickham et al., (2019). Welcome to the tidyverse. Journal of Open
Source Software, 4(43), 1686, https://doi.org/10.21105/joss.01686
A BibTeX entry for LaTeX users is
@Article{,
title = {Welcome to the {tidyverse}},
author = {Hadley Wickham and Mara Averick and Jennifer Bryan and Winston Chang and Lucy D'Agostino McGowan and Romain François and Garrett Grolemund and Alex Hayes and Lionel Henry and Jim Hester and Max Kuhn and Thomas Lin Pedersen and Evan Miller and Stephan Milton Bache and Kirill Müller and Jeroen Ooms and David Robinson and Dana Paige Seidel and Vitalie Spinu and Kohske Takahashi and Davis Vaughan and Claus Wilke and Kara Woo and Hiroaki Yutani},
year = {2019},
journal = {Journal of Open Source Software},
volume = {4},
number = {43},
pages = {1686},
doi = {10.21105/joss.01686},
}
Procedimientos estadísticos de uso común.
Comparación de medias en dos grupos
El archivo calamares.csv contiene datos de una muestra de 200 calamares en los que se han medido las variables longitud del manto (manto.mm), peso (peso.gr), diámetro del estatolito (statol.mm) y sexo.
Objetivo del estudio: ¿Hay diferencias significativas entre machos y hembras en los valores medios de las tres variables biométricas consideradas (longitud, peso y diámetro del estatolito)?
Descripción de datos.
Comenzamos con una estadística descriptiva de las variables para tener una idea de su comportamiento.
En primer lugar leemos los datos en R mediante la siguiente sintaxis:
Mostrar código R
# Lectura de datoscalamares <-read.csv2("calamares.csv")
Evaluamos si las variables son normales para cada sexo:
Podemos inspeccionar gráficamente la distribución de la variable manto.mm:
La distribución no parece normal para ninguno de los dos sexos. Aplicamos el test de Shapiro-Wilk para realizar el contraste de hipótesis de normalidad y tomar nuestra decisión basándonos en el p-valor:
El test de Shapiro-Wilk contrasta las siguientes hipótesis:
\[\begin{cases}
H_{0}: & \textrm{La variable es normal}\\
H_{1}: & \textrm{La variable NO es normal}
\end{cases}\]
Como vemos, en ambos sexos el p-valor es menor que el nivel de significación de referencia (0.05); por tanto no puede admitirse que la variable longitud del manto sea normal en ninguno de los dos sexos. En consecuencia describiremos esta variable mediante su mediana y cuartiles:
Como vemos, solo para el estatolito puede admitirse la normalidad en ambos sexos. Por esta razón esta variable se describe en términos de media y desviación típica:
Para comparar las medias de estas variables entre machos y hembras deberemos usar el test de Wilcoxon para la longitud del manto y el peso (por ser no normales) y el t-test para el diámetro del estatolito (que sí es normal):
Los p-valores nos indican que no existen diferencias significativas en la longitud del manto entre machos y hembras, y sí que existen en el peso y el diámetro del estatolito. En un artículo este resultado se redactaría de forma similar a la siguiente:
Significant differences have been detected between male and female squids in statolite diameter (Student’s t-test p<0.0001) and weight (Wilcoxon test p<0.0001). No differences have been detected in mantle length (Wilcoxon test p=0.526).
También es frecuente que los resultados se recojan en el artículo en una tabla como la siguiente, que reúne los estadísticos descriptivos y los p-valores de las comparaciones realizadas:
Para este ejercicio utilizaremos los datos del archivo sargos.csv correspondientes a una muestra de sargos capturados en las distintas islas del archipiélago canario. Mostramos a continuación los primeros datos del archivo:
Queremos decidir si la longitud media de los sargos es la misma en todas las islas o si hay diferencias significativas entre ellas. Un diagrama de cajas (boxplot) nos permite visualizar la situación:
El p-valor del test confirma que efectivamente podemos aceptar que la variable long es normal en todas las islas. Describimos la longitud en términos de media y desviación típica:
Dado que son siete grupos a comparar y la variable es normal, podemos utilizar el F-test (test F de Fisher); es conveniente convertir en factor la variable que define los grupos, ya que muchos de los procedimientos de R relacionados con este tratamiento de datos (aov, tukeyHSD, levene_test) requieren que dicha variable sea de tipo factor:
Mostrar código R
sargos <- sargos %>%mutate(isla=factor(isla))
Para la validez del F-test debemos comprobar además la homoscedasticidad (igualdad de varianzas en todos los grupos). Esto lo comprobamos mediante el test de Levene:
Mostrar código R
sargos %>%levene_test(long~isla)
# A tibble: 1 × 4
df1 df2 statistic p
<int> <int> <dbl> <dbl>
1 6 193 1.47 0.192
Puesto que el p-valor es mayor que 0.05 podemos admitir la homoscedasticidad de los datos. Como además ya hemos comprobado la normalidad (estos son los dos requisitos del F-test), podemos proceder a su realización:
Df Sum Sq Mean Sq F value Pr(>F)
factor(isla) 6 58.9 9.81 0.642 0.697
Residuals 193 2950.1 15.29
El p-valor 0.697 es mayor que 0.05 por lo que no hay evidencia de que la longitud media de los sargos difiera entre islas.
Si hubiese fallado la normalidad o la homoscedasticidad deberíamos haber aplicado el test de Kruskal-Wallis. No es necesario en este caso, pero mostramos a continuación como se realizaría. Simplemente:
Mostrar código R
kruskal.test(long~factor(isla),data=sargos)
Kruskal-Wallis rank sum test
data: long by factor(isla)
Kruskal-Wallis chi-squared = 5.0138, df = 6, p-value = 0.542
El p-valor (0.542) indica también que el peso de los sargos no muestra diferencias significativas entre islas.
Comparaciones múltiples
En el ejemplo anterior no hemos detectado diferencias entre islas por lo que no es necesario ningún análisis adicional. A continuación utilizaremos los datos del archivo Crabs.txt correspondientes a una muestra de cangrejos en los que se han medido diversas variables. Mostramos los primeros datos del archivo:
En particular nos interesa determinar si el peso medio de estos cangrejos depende de su color. Como en el caso anterior, comenzamos el análisis realizando un diagrama de cajas que nos describa la situación:
Supongamos por el momento que la distribución de la variable peso fuera normal y homoscedástica (con igual varianza) en los cuatro grupos. En tal caso podríamos realizar el test F de Fisher:
Df Sum Sq Mean Sq F value Pr(>F)
color 3 3766792 1255597 3.966 0.00917 **
Residuals 169 53502001 316580
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Vemos que en este caso hay diferencias significativas entre grupos (pues p<0.05), pero el test F no nos indica qué grupo es distinto de los demás. Para ello necesitamos llevar a cabo un test de Tukey:
Como vemos, la única diferencia significativa es la existente entre los cangrejos del color 4 y los del color 2. Concretamente la tabla nos indica que en esta muestra los cangrejos del color 4 pesaron por término medio 363.95 gramos menos que los del color 2; y que con un 95% de confianza lo que podemos afirmar a partir de los datos es que, en la población los cangrejos del color 4 pesan entre 709.78 y 18.52 gramos menos que los del color 2.
Pero en realidad no sabemos si el análisis anterior es correcto, pues no hemos comprobado normalidad ni homoscedasticidad. Si lo hacemos obtenemos:
Mostrar código R
shapiro.test(residuals(lm(weight~color,Crabs)))
Shapiro-Wilk normality test
data: residuals(lm(weight ~ color, Crabs))
W = 0.96842, p-value = 0.0005711
Por tanto el peso de estos cangrejos no es normal en los distintos grupos.
Mostrar código R
Crabs %>%levene_test(weight~color)
# A tibble: 1 × 4
df1 df2 statistic p
<int> <int> <dbl> <dbl>
1 3 169 1.87 0.137
El peso sí que es homoscedástico. En la práctica, el análisis de la varianza es más sensible a heteroscedasticidad que a la falta de normalidad; como en este caso se cumple la homoscedasticidad y los datos no presentan fuertes asimetrías (la no normalidad no es muy acusada), el resultado del análisis de la varianza es posiblemente válido (y muchas revistas nos lo aceptarían si pretendiésemos publicar el resultado).
No obstante, si la falta de normalidad fuese muy acusada, el test F anterior no sería correcto. En su lugar debemos aplicar el test de Kruskal-Wallis:
Mostrar código R
kruskal.test(weight~color,data=Crabs)
Kruskal-Wallis rank sum test
data: weight by color
Kruskal-Wallis chi-squared = 12.563, df = 3, p-value = 0.005685
El p-valor (0.0057) confirma que existen diferencias significativas en el peso entre los cangrejos de distinto color. El test de Tukey que hemos utilizado antes para localizar entre qué colores hay diferencias solo es adecuado cuando la variable bajo consideración es normal. Dado que no es el caso, dicho test no puede emplearse, y en su lugar podemos utilizar el test de Conover:
Este test nos confirma la existencia de diferencias significativas entre los colores 2 y 4 (p=0.0314) y además detecta diferencias significativas entre el 1 y el 4 (p=0.0457).
Dada la falta de normalidad de los valores de la variable peso, la describiríamos en términos de mediana y cuartiles:
El problema de la comparación de grupos mediante el test de Kruskal-Wallis y el posterior test de Conover es que los valores obtenidos no son directamente interpretables. Sí que tenemos un p-valor que nos indica si la diferencia entre grupos es o no significativa, pero no tenemos un indicador de la magnitud de dicha diferencia. Así, el test de Conover en el ejemplo anterior nos dice que la diferencia de peso entre cangrejos del color 2 y del color 4 es significativa (p=0.0314); pero la magnitud de dicha diferencia se evalúa en términos de diferencia media entre rangos, con el valor -32.54498; esto solo nos indica que el peso en los cangrejos de color 4 es inferior al peso de los cangrejos del color 2, pero no nos dice en cuanto, ya que la escala de rangos no tiene fácil interpretación. La tabla anterior, que nos muestra las medianas y rangos intercuartílicos en cada grupo nos permite ver que la mediana en el grupo 4 (2125) es menor que en el grupo 2 (2500), es decir, la diferencia entre ambas medianas es de 375 gramos en estas muestras; pero no sabemos cuál es el margen de error en el que nos movemos si pretendemos trasladar este resultado a la población.
No obstante, si convertimos el peso de los calamares a escala logarítmica, podemos comprobar que en dicha escala es admisible la normalidad y homoscedasticidad del peso en los distintos grupos:
# A tibble: 1 × 4
df1 df2 statistic p
<int> <int> <dbl> <dbl>
1 3 169 1.59 0.194
Dado que se cumplen las hipótesis de normalidad y homoscedasticidad podemos proceder con el análisis de la varianza. El resultado, en este caso, nos indica que hay diferencia significativas entre los pesos medios (medidos en escala logarítmica), pues p=0.0068:
Mostrar código R
summary(model)
Df Sum Sq Mean Sq F value Pr(>F)
color 3 0.670 0.22320 4.199 0.00678 **
Residuals 169 8.984 0.05316
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
El test de Tukey es también válido en este caso, y nos proporciona el siguiente resultado:
Como vemos, se confirma la diferencia significativa entre los cangrejos del color 4 y los del color 2. La diferencia entre los pesos medios en escala logarítmica entre ambas clases es -0.148. Pero nuevamente este es un valor difícil de interpretar, ya que sus unidades no son directamente convertibles en gramos. No obstante, el paquete emmeans de R permite obtener una comparación más entendible. En principio, obtenemos el mismo resultado para el test de Tukey (salvo que realiza las comparaciones en otro orden):
contrast estimate SE df lower.CL upper.CL t.ratio p.value
color1 - color2 0.0528 0.0706 169 -0.13049 0.236 0.747 0.8777
color1 - color3 0.1525 0.0751 169 -0.04238 0.347 2.030 0.1810
color1 - color4 0.2010 0.0827 169 -0.01367 0.416 2.430 0.0754
color2 - color3 0.0997 0.0420 169 -0.00944 0.209 2.370 0.0868
color2 - color4 0.1482 0.0546 169 0.00668 0.290 2.717 0.0362
color3 - color4 0.0486 0.0602 169 -0.10764 0.205 0.807 0.8512
Results are given on the log (not the response) scale.
Confidence level used: 0.95
Conf-level adjustment: tukey method for comparing a family of 4 estimates
P value adjustment: tukey method for comparing a family of 4 estimates
Pero ahora podemos especificar que queremos el resultado en la misma escala de la variable respuesta, y en tal caso obtenemos:
Mostrar código R
pairs(em.model,infer=TRUE, type="response")
contrast ratio SE df lower.CL upper.CL null t.ratio p.value
color1 / color2 1.05 0.0745 169 0.878 1.27 1 0.747 0.8777
color1 / color3 1.16 0.0875 169 0.959 1.42 1 2.030 0.1810
color1 / color4 1.22 0.1012 169 0.986 1.52 1 2.430 0.0754
color2 / color3 1.10 0.0465 169 0.991 1.23 1 2.370 0.0868
color2 / color4 1.16 0.0633 169 1.007 1.34 1 2.717 0.0362
color3 / color4 1.05 0.0632 169 0.898 1.23 1 0.807 0.8512
Confidence level used: 0.95
Conf-level adjustment: tukey method for comparing a family of 4 estimates
Intervals are back-transformed from the log scale
P value adjustment: tukey method for comparing a family of 4 estimates
Tests are performed on the log scale
Como vemos, en este caso, en lugar de calcular la diferencia media entre los pesos de los cangrejos de colores 2 y 4, el procedimiento calcula el cociente entre dichas medias, que resulta ser 1.16, con un intervalo de confianza al 95% dado por [1.007, 1.34]; esto significa que el peso medio de los cangrejos de color 2 es entre un 0.7% y un 34% mayor que el peso medio de los de color 4, siendo dicha diferencia significativa (p=0.0362). No tenemos la diferencia en gramos, pero sí la tenemos evaluada en porcentaje lo que, en cualquier caso, es mejor y desde luego más interpretable, que tenerla evaluada simplemente en “rango medio” como ocurría con el test de Conover.
