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 una población en dos momentos diferentes.

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:

  1. Instalar paquetes:
Paquetes <- c("FSA", 
              "tidyverse", 
              "magrittr", 
              "lubridate", 
              "car")

install.packages(Paquetes)
  1. Extraer paquetes:
library("FSA")
library("tidyverse")
library("lubridate")
library("magrittr")
library("car")

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)

  • Paso 1: Fijar el directorio de trabajo, carpeta de su ordenador que contiene los datos descargados.
setwd("D:/Documentos/Clase_Repro")
  • Paso 2: Extraer los datos
df <- read.csv("PracticaRepro.csv")

Alternativamente usar el siguiente código (No recomendado, si la intención es trabajar con datos propios):

df <- read.csv("https://raw.githubusercontent.com/Fertapia/EcologiaReproductiva/master/PracticaRepro.csv")
str(df)
## 'data.frame':    158 obs. of  6 variables:
##  $ Fecha   : chr  "9/02/2003" "10/07/2002" "7/18/2000" "6/11/2001" ...
##  $ 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 : chr  "Inmaduro" "Inmaduro" "Inmaduro" "Inmaduro" ...
##  $ Estado  : int  1 1 1 2 2 1 1 1 1 1 ...
##  $ Sexo    : chr  "Hembra" "Hembra" "Hembra" "Hembra" ...

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 <- df %>% 
  mutate(Fecha = as.POSIXct(Fecha, format = "%m/%d/%Y"),
         Madurez = as.factor(Madurez),
         Sexo = as.factor(Sexo)) 
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 <- 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 <- 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 <- df %>% 
  mutate(lcat2=lencat(Longitud, w=2)) # el argumento w indica el ancho del intervalo en este caso es de 2 cm

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 <- 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

Graficamos la proporciones calculadas en intervalos de longitud de 2 cm:

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 longitud toma los valores de: 1 cm, 3 cm y 4 cm?
  3. ¿Qué pasaría si en lugar de Longitud consideramos la Edad?, Pista: el argumento w de la función lencat ahora tiene que manejarse en años.

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)

Nota: La variable categórica en este caso Madurez debe ser NECESARIAMENTE un factor en R, de no ser así transformar la variable a factor con la función as.factor.

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))

La pendiente estimada (0.4371) y el término independiente son básicos para el cálculo de la probabilidad en que un individuo al azar de la población estudiada 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 o 4%, mientras que un individuo de 42 cm tendría una probabilidad de 0.80 o 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.004722876 0.158841476

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

Ejercicios

  1. ¿Cuál es el significado de la función bootCase qué pasa si toma un valor de 1000?. Pista: Revisar la documentación de la función en https://www.rdocumentation.org/packages/car/versions/2.1-6/topics/bootCase para obtener más información de esta.

  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 asumida 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.

La función binomial_smooth permite ajustar la curva de forma logística en ggplot.

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 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.15767 40.02525
bL90 <- apply(bcL, 1, lrPerc, p=0.9)
(L90ci <- quantile(bL90, c(0.025, 0.975)))
##     2.5%    97.5% 
## 41.58129 45.47565

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. ¿Es distinta de la curva basada en la longitud?.
  3. Si tuviera que elegir entre la Talla de primera madurez y la Edad de primera madurez para proponer medidas de manejo pesquero para la conservación de este recurso ¿cuál escogería?. Justifique su respuesta (Máximo 200 palabras).

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 (variable era), evalúa la influencia de estas dos variables sobre los datos y si calcula una curva de madurez para cada periodo.

Nota: La variable categórica que separa los grupos de análisis debe ser un factor NECESARIAMENTE. De no ser así, transformar la variable con el argumento as.factor().

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

La variante está 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,]      -36.90385 0.9282711    27.427459           -0.67464480
## [2,]      -25.87431 0.6592949     7.079337           -0.17428860
## [3,]      -38.95628 0.9785044    24.749070           -0.61464428
##  [998,]   -27.58395 0.7099645    17.081991           -0.43268765
##  [999,]   -34.16044 0.8796750    -4.531490            0.09050686
## [1000,]   -45.16210 1.1326608    26.623216           -0.64797184
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.556

Estimamos los intervalos de confianza para ambos periodos

( ci.L50.pre <-  quantile(L50.pre,c(0.025,0.975)) )
##     2.5%    97.5% 
## 37.08374 40.93457
( ci.L50.post <- quantile(L50.post,c(0.025,0.975)) )
##     2.5%    97.5% 
## 35.34354 40.26091

Resumen de los resultados

(L50.pre.Fin <- c(L50.pre.m, ci.L50.pre))
##              2.5%    97.5% 
## 38.95635 37.08374 40.93457
(L50.post.Fin <- c(L50.post.m, ci.L50.post))
##              2.5%    97.5% 
## 38.02492 35.34354 40.26091

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 los resultados observados durante el ejercicio (entre eras Pre-2002 y Post-2002), si la longitud de madurez se incrementara fuera diferente entre estos periodos. ¿Cuál sería su significado evolutivo?. Justifique su respuesta (Máximo 200 palabras).

Fernando Tapia - Laboratorio de rerpoducción y Biología del Desarrollo UNMSM . 18 de Noviembre del 2018. Última actualización: 2020-10-03