Práctica: Ecología reproductiva de organismos acuáticos

Objetivos:

  • Brindar un enfoque general de los métodos de evaluación en ecología reproductiva de organismos acuáticos con un enfoque práctico en pesquerías.
  • Modelar la curva de madurez sexual de dos poblaciones.

La estrategia reproductiva de una especie es la asociación de rasgos fenotípicos que se definen en base a relaciones complejas con el medio ambiente y la historia evolutiva del individuo. Esta estrategia reproductiva forma parte de la historia de vida de la especie la misma que define la forma en que los organismos optimizan sus recursos para maximizar su supervivencia y la de su descendencia.

El estudio de la ecología reproductiva de especie comprende el análisis de los que son llamados rasgos de historia de vida y como éstos interactúan, entre los más relevantes que definen la estrategia reprductiva están el tamaño del nacimiento, edad y tamaño de madurez, número, tamaño y sexo de la descendencia; el esfuerzo reproductivo, la tasa de supervivencia de la descendencia y los efectos maternales.

El cálculo de la edad y/o tamaño en la que los organismos alcanzan la madurez es un parámetro importante en el manejo pesquero ya que es un factor a tomar en consideración para establecer la talla mínima de captura, bajo el enfoque de asegurar que el recurso se reproduzca al menos una vez antes de formar parte de la captura de este modo “asegurando” la sostenibilidad.

Asimismo, los efectos maternales se dan cuando las características fenotípicas de la hembra influyen directamente sobre el fenotipo de la descendencia, la forma más típica de estos efectos maternales se da en relación a la fluctuaciones en el tamaño del huevo y el número de huevos producidos por el individuo. En términos generales, se espera que hembras más grandes produzcan huevos más grandes y un mayor número de estos, sin embargo, también es posible observar “intercambios” entre el tamaño del huevo y el número de crías en orden como una forma de adaptación a las condiciones ambientales.

Paquetes requeridos

Las funciones necesarias para la práctica se encuentran en los siguientes paquetes:

library(FSA)
library(magrittr)
library(dplyr)
library(lubridate)
library(car)
library(ggplot2)

Base de datos empleada

Las variables a analizar corresponden a longitud total (longitud; en cm), edad (años), y madurez (Inmaduro y Maduro) de las hembras del pez piedra (Sebastes ruberimus) registrados en la base de datos PracticaRepro.csv (Descargar)

df <- read.csv("PracticaRepro.csv")
str(df)
## 'data.frame':    158 obs. of  6 variables:
##  $ Fecha   : Factor w/ 71 levels "10/01/2003","10/02/2003",..: 64 7 39 27 49 4 38 55 35 40 ...
##  $ Longitud: int  31 32 32 32 32 33 33 34 34 34 ...
##  $ Edad    : int  10 6 11 11 13 9 10 8 10 11 ...
##  $ Madurez : Factor w/ 2 levels "Inmaduro","Maduro": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Estado  : int  1 1 1 2 2 1 1 1 1 1 ...
##  $ Sexo    : Factor w/ 2 levels "Hembra","Macho": 1 1 1 1 1 1 1 1 1 1 ...

Adicionalmente se ha registrado la fecha en la que los datos fueron registrados, la fecha de referencia es el año 2002. A continuación, se transformará la fecha en un factor de dos niveles: Pre-2002 y Post-2002 este último incluye aquellos individuos capturados en 2002.

headtail(df)
##          Fecha Longitud Edad  Madurez Estado   Sexo
## 1    9/02/2003       31   10 Inmaduro      1 Hembra
## 2   10/07/2002       32    6 Inmaduro      1 Hembra
## 3    7/18/2000       32   11 Inmaduro      1 Hembra
## 156  8/18/2002       67   50   Maduro      6  Macho
## 157 10/07/2002       68   88   Maduro      6  Macho
## 158  4/23/2001       70   66   Maduro      4  Macho
df %<>% mutate(Fecha = as.POSIXct(Fecha, format = "%m/%d/%Y"))
str(df)
## 'data.frame':    158 obs. of  6 variables:
##  $ Fecha   : POSIXct, format: "2003-09-02" "2002-10-07" ...
##  $ Longitud: int  31 32 32 32 32 33 33 34 34 34 ...
##  $ Edad    : int  10 6 11 11 13 9 10 8 10 11 ...
##  $ Madurez : Factor w/ 2 levels "Inmaduro","Maduro": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Estado  : int  1 1 1 2 2 1 1 1 1 1 ...
##  $ Sexo    : Factor w/ 2 levels "Hembra","Macho": 1 1 1 1 1 1 1 1 1 1 ...
headtail(df)
##          Fecha Longitud Edad  Madurez Estado   Sexo
## 1   2003-09-02       31   10 Inmaduro      1 Hembra
## 2   2002-10-07       32    6 Inmaduro      1 Hembra
## 3   2000-07-18       32   11 Inmaduro      1 Hembra
## 156 2002-08-18       67   50   Maduro      6  Macho
## 157 2002-10-07       68   88   Maduro      6  Macho
## 158 2001-04-23       70   66   Maduro      4  Macho

Mediante el paquete lubridate extraemos el año con year(), y la clasificamos en Pre-2002 y Post-2002 el cual es luego añadido a la base de datos como una variable era.

df %<>% mutate(year = year(Fecha), era = ifelse(year < 2002, "Pre-2002", "Post-2002"), 
    era = factor(era, levels = c("Pre-2002", "Post-2002")))
str(df)
## 'data.frame':    158 obs. of  8 variables:
##  $ Fecha   : POSIXct, format: "2003-09-02" "2002-10-07" ...
##  $ Longitud: int  31 32 32 32 32 33 33 34 34 34 ...
##  $ Edad    : int  10 6 11 11 13 9 10 8 10 11 ...
##  $ Madurez : Factor w/ 2 levels "Inmaduro","Maduro": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Estado  : int  1 1 1 2 2 1 1 1 1 1 ...
##  $ Sexo    : Factor w/ 2 levels "Hembra","Macho": 1 1 1 1 1 1 1 1 1 1 ...
##  $ year    : num  2003 2002 2000 2001 2000 ...
##  $ era     : Factor w/ 2 levels "Pre-2002","Post-2002": 2 2 1 1 1 2 1 2 1 1 ...
headtail(df)
##          Fecha Longitud Edad  Madurez Estado   Sexo year       era
## 1   2003-09-02       31   10 Inmaduro      1 Hembra 2003 Post-2002
## 2   2002-10-07       32    6 Inmaduro      1 Hembra 2002 Post-2002
## 3   2000-07-18       32   11 Inmaduro      1 Hembra 2000  Pre-2002
## 156 2002-08-18       67   50   Maduro      6  Macho 2002 Post-2002
## 157 2002-10-07       68   88   Maduro      6  Macho 2002 Post-2002
## 158 2001-04-23       70   66   Maduro      4  Macho 2001  Pre-2002

Removemos, los valores perdidos NA presentes en la base de datos, estos se han registrado sólo en la variable Madurez.

df %<>% filterD(!is.na(Madurez))

Análisis exploratorio de la madurez

Los datos de madurez son registrados como Maduro e Inmaduro, estos datos son comparados respecto de la Longitud o la Edad. Además es posible comparar el comportamiento de la madurez respecto de la época de muestreo (Era) o el Sexo.

Los datos de madurez deben ser resumidos en base a la proporción de individuos maduros en cada categoría de edad o longitud. Estas categorías de longitud son intervalos de tamaño constante, el tamaño del invervalo es variable y debe tener suficientes individuos para calcular una proporción confiable.

Los datos presentados son clasificados en categorías de 2 cm lo cual asegura que en cada intervalo hay aproximadamente 10 peces o más, estas categorías son calculadas con lencat.

df %<>% mutate(lcat2 = lencat(Longitud, w = 2))
headtail(df)
##          Fecha Longitud Edad  Madurez Estado   Sexo year       era lcat2
## 1   2003-09-02       31   10 Inmaduro      1 Hembra 2003 Post-2002    30
## 2   2002-10-07       32    6 Inmaduro      1 Hembra 2002 Post-2002    32
## 3   2000-07-18       32   11 Inmaduro      1 Hembra 2000  Pre-2002    32
## 146 2002-08-18       67   50   Maduro      6  Macho 2002 Post-2002    66
## 147 2002-10-07       68   88   Maduro      6  Macho 2002 Post-2002    68
## 148 2001-04-23       70   66   Maduro      4  Macho 2001  Pre-2002    70

Calculamos la frecuencia de individuos maduros en cada categoría de longitud con xtabs(). Las frecuencias son calculadas con prop.table, y estas proporciones son graficadas.

freq <- xtabs(~lcat2 + Madurez, data = df)
freq
##      Madurez
## lcat2 Inmaduro Maduro
##    30        1      0
##    32        6      0
##    34        5      0
##    36        5      4
##    38        5      3
##    40        4      8
##    42        1     12
##    44        1     12
##    46        1     27
##    48        1     15
##    50        0      6
##    52        0     10
##    54        0      2
##    56        0      5
##    58        0      2
##    60        0      5
##    62        0      1
##    64        0      3
##    66        0      1
##    68        0      1
##    70        0      1
props <- prop.table(freq, margin = 1)
round(props, 3)  #Redondea los resultados a 3 decimales
##      Madurez
## lcat2 Inmaduro Maduro
##    30    1.000  0.000
##    32    1.000  0.000
##    34    1.000  0.000
##    36    0.556  0.444
##    38    0.625  0.375
##    40    0.333  0.667
##    42    0.077  0.923
##    44    0.077  0.923
##    46    0.036  0.964
##    48    0.062  0.938
##    50    0.000  1.000
##    52    0.000  1.000
##    54    0.000  1.000
##    56    0.000  1.000
##    58    0.000  1.000
##    60    0.000  1.000
##    62    0.000  1.000
##    64    0.000  1.000
##    66    0.000  1.000
##    68    0.000  1.000
##    70    0.000  1.000
props.df <- as.data.frame(props)
headtail(props.df)
##    lcat2  Madurez Freq
## 1     30 Inmaduro    1
## 2     32 Inmaduro    1
## 3     34 Inmaduro    1
## 40    66   Maduro    1
## 41    68   Maduro    1
## 42    70   Maduro    1
props.df %<>% filter(Madurez == "Maduro")
headtail(props.df)
##    lcat2 Madurez Freq
## 1     30  Maduro    0
## 2     32  Maduro    0
## 3     34  Maduro    0
## 19    66  Maduro    1
## 20    68  Maduro    1
## 21    70  Maduro    1
lcat1 <- as.numeric(levels(props.df$lcat2))
props.df$lcat1 <- lcat1
ggplot(props.df, aes(x = lcat2, y = Freq)) + geom_point() + labs(title = paste("Proporción de individuos aclanzando la madurez"), 
    x = "Longitud (cm)", y = "Frecuencia") + theme_bw()

Ejercicios:

  1. ¿Cómo interpreta la figura resultante?
  2. ¿Cómo cambia la figura cuando el intervalo de la categoría toma los valores de: 1 cm, 2.5 cm, 3 cm y 4 cm?

Modelando los datos

Ajuste del modelo de regresión logística

Usaremos un modelo linear generalizado mediante ajustando una curva logística, este modelo requiere de una variable de respuesta binomial (Maduro - Inmaduro), expresada como una variable cuantitativa. El objetivo es calcular la probabilidad de “éxito” del individuo de alcanzar la madurez (\(p\)) definida entre 0 y 1; mientras que la variable explicatoria es la longitud o edad. La relación entre \(p\) y la longitud generalmente no es linear.

En este sentido ajustamos una regresión logística con el comando glm(), la fórmula tine la forma: factor~cuantitativa, data = "matriz que contiene las variables", family = distribución a priori:

glm1 <- glm(Madurez ~ Longitud, data = df, family = binomial)

Luego extraemos los parámetros del modelo con coef()

coef(glm1)
## (Intercept)    Longitud 
## -16.9482593   0.4371786

Calculamos los intervalos de confianza mediante bootstrap con el comando bootCase del paquete car.

bcL <- bootCase(glm1, B = 100)
cbind(Ests = coef(glm1), confint(bcL))
##                    Ests     95% LCI     95% UCI
## (Intercept) -16.9482593 -32.5871502 -12.4382428
## Longitud      0.4371786   0.3262795   0.8287497

La pendiente estimada (0.4371) indica que con un incremento de 1 cm en la longitud del pez se incrementa por un 0.4371 la probabilidad de que el individuo haya alcanzado la madurez. Los intervalos de confianza al 95% superior (UCI) e inferior (LCI) indican entre qué valores se encuentran la pendiente y el intercepto del modelo.

Calculamos la probabilidad de que el individuo esté maduro

A continuación calculamos la probabilidad de que un individuo de Longitud o edad conocida haya alcanzado la madurez. Esta predicción se realiza mediante el comando predict(), este requiere el modelo glm1, un data frame con los valores e la variable explicatoria para la cual deseams realizar la predicción como segundo argumento y type = "response" fuerza el cálculo de la probabilidad de estar maduro.

predict(glm1, data.frame(Longitud = c(32, 42)), type = "response")
##         1         2 
## 0.0493342 0.8042766

Los resultados obtenidos indican que la probabilidad de que un individuo de 32 cm haya alcanzado la madurez es de 0.04 ó 4%, mientras que un individuo de 42 cm tendría una probabilidad de 0.80 ó 80%. A continuación calculamos los intervalos de fonfianza para estas predicciones. Para ello creamos una función que calcula los intervalos de confianza para cada una de las muestras del procedimiento bootstrap.

predP <- function(cf, x) exp(cf[1] + cf[2] * x)/(1 + exp(cf[1] + cf[2] * x))
p32 <- apply(bcL, 1, predP, x = 32)
quantile(p32, c(0.025, 0.975))
##        2.5%       97.5% 
## 0.001738856 0.120809018

Entonces la probabilidad de que un individuo de 32 cm alcance la madurez está entre 0.001 y 0.138.

Ejercicios

  1. ¿Qué sucede con los intervalos de confianza si en lugar de 100 repeticiones empleamos, 10, 50, 500, y 1000 repeticiones?

  2. ¿Cuál es la probabilidad de que un individuo de 12 cm y otro de 65 cm hayan alcanzado la madurez? Calcular sus intervalos de confianza.

  3. ¿Porqué el tipo de distribución asumido es de tipo binomial?

Resumen gráfico

Mediante una figura ilustramos el ajuste de la regresión calculada en los pasos anteriores. Primero que nada recordar la variable Madurez corresponde a un factor por lo que debe ser transformada a una variable numérica de valores 0 (Inmaduro) y 1 (Maduro), notar que se le extrae el valor de 1 ya que al transformarlos a numéricos los valores que toman son de 1 y 2.

El paquete ggplot2 no contiene la opción de ajustar curvas binomiales por defecto, sin embargo es posible crear una función que le indique el tipo de modelo y distribución empleados.

binomial_smooth <- function(...) {
    geom_smooth(method = "glm", method.args = list(family = "binomial"), ...)
}
ggplot(df, aes(Longitud, as.numeric(Madurez) - 1)) + binomial_smooth() + geom_point(position = position_jitter(height = 0.03, 
    width = 0)) + xlab("Longitud (cm)") + ylab("Probabilidad de individuo maduro") + 
    theme_bw()

Longitud - o edad - de primera madurez

La longitud y/o la edad de primera madurez es una medida empleada como referencia en el manejo pesquero, está definida como la longitud o edad en la cual la probabilidad de ser maduro es del 50%. Otro valor usualmente calculado es al 90%. Estos cálculos son realizados mediante la función:

lrPerc <- function(cf, p) (log(p/(1 - p)) - cf[[1]])/cf[[2]]

La función toma los coeficientes de glm() como pirmer argumento y la probabilidad (\(p\)) como segundo argumento. Calculamos la longitud de madurez al 50% y al 90%.

(L50 <- lrPerc(coef(glm1), 0.5))
## [1] 38.76736
(L90 <- lrPerc(coef(glm1), 0.9))
## [1] 43.79328

Calculamos los intervalos de confianza para cada uno de los valores:

bL50 <- apply(bcL, 1, lrPerc, p = 0.5)
(L50ci <- quantile(bL50, c(0.025, 0.975)))
##     2.5%    97.5% 
## 37.41832 39.79788
bL90 <- apply(bcL, 1, lrPerc, p = 0.9)
(L90ci <- quantile(bL90, c(0.025, 0.975)))
##     2.5%    97.5% 
## 41.50798 45.24422

Entonces podemos decir que la longitud de madurez al 50% es de 38.76 cm y que esta se encuentra entre 37.67 y 40.02 cm con una confianza del 95%.

Graficamos nuevamente:

ggplot(df, aes(Longitud, as.numeric(Madurez) - 1)) + binomial_smooth() + geom_point(position = position_jitter(height = 0.03, 
    width = 0)) + geom_point(data = props.df, aes(x = lcat1, y = Freq), size = 2, 
    shape = 3, color = "darkgreen") + xlab("Longitud (cm)") + ylab("Probabilidad de individuo maduro") + 
    geom_segment(aes(x = L50, y = 0, xend = L50, yend = 0.5), color = "darkred") + 
    geom_segment(aes(x = min(df$Longitud), y = 0.5, xend = L50, yend = 0.5), 
        color = "darkred") + theme_bw()

Ejercicios

  1. ¿Porqué la talla de primera madurez es la longitud en la que \(p\) es 0.5?

  2. Calcular la edad de primera madurez

    2.1 ¿Cuál sería el tamaño de intervalo de edad adecuado?

    2.2 ¿Varía la curva de edad de primera madurez respecto de la curva basada en la longitud?

Comparando las tallas de madurez entre grupos

Ajustando el modelo

Según el tipo de estudio puede ser importante definir si existen diferencias entre los valores de longitud de madurez entre épocas, sexos, especies, etc. Para ello se puede determinar si los parámetros de la regresión logística difieren significativamente entre las situaciones. En el presente ejercicio compara dos periodos de manejo pesquero distitntos: Pre-2002 y Post-2002.

glm2 <- glm(Madurez ~ Longitud * era, data = df, family = binomial)

La variante están en Longitud*era esto le indica al modelo que existe un factor que identifica a los grupos que nos interesa comparar era.

La significancia de las diferencias entre ambas épocas son calculadas mediante el comando Anova del paquete car

Anova(glm2)
## Analysis of Deviance Table (Type II tests)
## 
## Response: Madurez
##              LR Chisq Df Pr(>Chisq)    
## Longitud       68.599  1    < 2e-16 ***
## era             0.005  1    0.94541    
## Longitud:era    3.100  1    0.07831 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El Análisis de desviación se interpreta de abajo hacia arriba:

  • Longitud:era es el término de interacción, en este caso no es significativa.

  • era tampoco es significante sugiere que los interceptos no difieren entre ambas eras

  • Longitud presenta diferencias signficativas a nivel de longitudes.

Comparando Longitudes de madurez

Primeramente examinamos los parámetros del modelo de regresión

coef(glm2)
##           (Intercept)              Longitud          eraPost-2002 
##           -27.1314345             0.6956840            14.0137840 
## Longitud:eraPost-2002 
##            -0.3510082

Notar que Intercept y Longitud son los parámetros del grupo de referencia (Pre-2002) y los dos parámetros siguientes corresponden al intercepto y pendiente de Post-2002. El orden es alfabético y puede obtenerse con el siguiente comando.

levels(df$era)
## [1] "Pre-2002"  "Post-2002"

A fin de calcular los intervalos de confianza para los parámetros construimos muestras bootstrap a partir del glm().

bcL2 <- bootCase(glm2, B = 1000)
headtail(bcL2)
##         (Intercept)  Longitud eraPost-2002 Longitud:eraPost-2002
## [1,]      -31.87989 0.8379844     20.77550            -0.5299764
## [2,]      -29.99933 0.7854998     21.33652            -0.5416723
## [3,]      -39.14602 0.9797610     27.58623            -0.6658919
##  [998,]   -30.39729 0.8136865     16.76848            -0.4697490
##  [999,]   -28.71490 0.7212629     14.07330            -0.3510304
## [1000,]   -28.33026 0.7259271     16.87153            -0.4207113
L50.pre = apply(bcL2[, 1:2], 1, lrPerc, p = 0.5)
L50.pre.m <- mean(L50.pre)
L50.post = apply(bcL2[, 1:2] + bcL2[, 3:4], 1, lrPerc, p = 0.5)
L50.post.m <- mean(L50.post)

Si no hay diferencia entre las L50 entre eras, uno podría esperar que los promedios de estos grupos sean equivalentes o parecidos

L50.diff <- L50.pre - L50.post

Evaluamos la significancia de las diferencias entre L50 pre y post

(p.L50.diff <- 2 * min(c(mean(L50.diff > 0), mean(L50.diff < 0))))
## [1] 0.544

Estimamos los intervalos de confianza para ambos periodos

(ci.L50.pre <- quantile(L50.pre, c(0.025, 0.975)))
##     2.5%    97.5% 
## 37.39148 40.95884
(ci.L50.post <- quantile(L50.post, c(0.025, 0.975)))
##     2.5%    97.5% 
## 35.32764 40.36290

Resumen de los resultados

(L50.pre.Fin <- c(L50.pre.m, ci.L50.pre))
##              2.5%    97.5% 
## 39.00523 37.39148 40.95884
(L50.post.Fin <- c(L50.post.m, ci.L50.post))
##              2.5%    97.5% 
## 38.04047 35.32764 40.36290

Resumen gráfico

Finalmente graficamos las curvas de madurez para ambos periodos:

ggplot(df, aes(Longitud, as.numeric(Madurez) - 1, color = era)) + binomial_smooth() + 
    geom_point(position = position_jitter(height = 0.03, width = 0)) + xlab("Longitud (cm)") + 
    ylab("Probabilidad de individuo maduro") + theme_bw()

Ejercicios

  1. Calcular si existen diferencias entre la edad de madurez al 50% entre temporadas.

  2. Calcular si existen diferencias entre la longitud de madurez entre sexos.

  3. Según su criterio, a qué se deberían las diferencias observadas durante el ejercicio, si la longitud de madurez se incrementara ¿Cuál sería su significado evolutivo?

Fernando Tapia - Laboratorio de reproducción y Biología del Desarrollo UNMSM . 18 de Noviembre del 2018. Última actualización: 2019-05-23