Duración del amplexo

Introducción.

La medición de la duración del amplexo fue uno de los aspectos que consideré para el proyecto de tesis que desarrollé. Para medir el amplexo utilicé el método descrito en el artículo de crustáceos que recientemente publicaron. Para conseguirlo, tomé y modifiqué un algoritmo del libro “The art of R programming” de Norman Matloff el cual usé con la información de presencia/ausencia de amplexos que recolecté a lo largo de mis experimentos. Los resultados indican que la duración del amplexo aumenta cuando la proporción de sexos se acerca a 0.5 y disminuye cuando la proporción de sexos está sesgada.

Programa para la duración del amplexo

Para medir la duración del amplexo utilicé un programa en R. De la sección 2.5.1 “Extended Example: Finding Runs of Consecutive Ones” del libro “The art of R programming”, tomé y modifiqué el algoritmo para medir la duración del amplexo como lo describen en el artículo de crustáceos reciéntemente publicado. El código del programa es el siguiente:

dura <- function(x) {#La función recibe un vector x formado por una secuencia
#de 1 y 0 (presencia o ausencia de amplexos). Esta función admite vectores de
#cualquier longitud, por lo que permite medir la duración del amplexo en una
#secuencia de cualquier cantidad de días
  n <- length(x)
  corridas <- numeric(1)
  dias <- numeric(1)
  dur <- numeric(1)
  k <- 1
  if (all(x[1:n]==0)) {
    dur <- 0
  } else {
    
    for (j in 1:n) {
      for (i in 1:(n-k+1)) {
        if (all(x[i:(i+k-1)]==1)) 
        {dias <- dias+k
        corridas <- corridas+1
        } else {}
      }
      k <- k+1
    }
    dur <- dias/corridas
  }  
  return(dur)
}

Una vez escrita la función a la que llamé dura(), hice pruebas para verificar que los cálculos fuesen correctos. Para fines ilustrativos de este reporte, comparé la estimación de las siguientes secuencias:

secuencia	corridas tamaño n=1	corridas tamaño n=2	corridas tamaño n=3	corridas tamaño n=4	corridas tamaño n=5	corridas tamaño n=6	corridas tamaño n=7	corridas tamaño n=8
10111	4	2	1	0	0	NA	NA	NA
1110111	6	4	2	0	0	0	0	NA
10011101	5	2	1	0	0	0	0	0

Haciendo los cálculos manuales para las secuencias de la tabla anterior, la estimación de la duración del amplexo es la siguiente para las tres secuencias: \[1° secuencia: (4corridas*1día + 2corridas*2días + 1corrida*3días)/(4corridas+2corridas+1corrida) = 1.571429días\] \[2° secuencia: (6corridas*1día + 4corridas*2días + 2corridas*3días)/(6corridas+4corridas+2corridas) = 1.666667días\] \[3° secuencia: (5corridas*1día + 2corridas*2días + 1corrida*3días)/(5corridas+2corridas+1corrida) = 1.5días\]

Ahora, para mostrar que el algoritmo funciona correctamente, simplemente tengo que llamar a la función dura() para las tres secuencias que estoy usando como ejemplo:

paste("La duración del amplexo estimada para la primer secuencia es de: ", dura(c(1,0,1,1,1)), "días")

## [1] "La duración del amplexo estimada para la primer secuencia es de:  1.57142857142857 días"

paste("La duración del amplexo estimada para la segunda secuencia es de: ", dura(c(1,1,1,0,1,1,1)), "días")

## [1] "La duración del amplexo estimada para la segunda secuencia es de:  1.66666666666667 días"

paste("La duración del amplexo estimada para la tercer secuencia es de: ", dura(c(1,0,0,1,1,1,0,1)), "días")

## [1] "La duración del amplexo estimada para la tercer secuencia es de:  1.5 días"

Comparando los resultados calculados a mano con los obtenidos mediante la función, se puede observar que el programa funciona correctamente. Aún así, para tener mayor seguridad, usé el programa para calcular la duración del amplexo de todas las permutaciones posibles de secuencias binarias de 8 días, como lo hicieron en el artículo de crustáceos reciéntemente publicado y obtuve el mismo histograma del modelo nulo.

library(gtools)
library(png)
library(grid)
library(ggplot2)
y <- 8
subsets <- permutations(n= 2, r= y, v= c(1,0), repeats.allowed = TRUE)
corridas <- function(x) {
  n <- length(x)
  runs <- vector(mode = "numeric", length = y)
  k <- 1
  for (i in 1:n){
  for (i in 1:(n-k+1)) {
    if (all(x[i:(i+k-1)]==1)) runs[k] <- runs[k] + 1
  }
  k <- k+1
  }
  return(runs)
}


totalruns <- function(a) {
  total <- vector(mode = "numeric", length = y)
  i <- 1
  for (i in 1:nrow(a)) {
    total <- total + corridas(c(a[i, 1:y]))
  }
  return(total)
}
cada_corrida <- totalruns(subsets)
sum(cada_corrida)

## [1] 1793

corridas_total <- function(z) {
  total <- NULL
  i <- 1
  for (i in 1:length(z)) {
    total <- c(total, rep(i,z[i]))
  }
  return(total)
}


simulacion <- function(n) {
  btrp <- vector(length = length(n))
  for (i in 1:n) {
  btrp[i] <- mean(sample(corridas_total(cada_corrida), replace = TRUE))
  }
  return(btrp)
}

dts <- round(simulacion(10000), digits = 2)

data <- as.data.frame(dts, rep("duración",10000))
a <- ggplot(data, aes(x = dts)) +
geom_bar() 
a + labs(  title = "Modelo de simulación",
  x = "Duración del amplexo en días",
       y = "Frecuencia",
       caption = "n = 10000")

Con esto, asumí que el programa que hice funciona correctamente.

Duración del amplexo estimado en los experimentos.

En mis experimentos, 7 observadores ajenos a la investigación (1 diferente para cada uno de los 7 bloques) registraron el número de amplexos observados en cada tratamiento tres veces al día, durante 5 días. Con esta información obtuve secuencias de números de amplexos observados en 5 días para cada observación las cuales transformé en secuencias binarias. Por ejemplo, la secuencia 32201 representa 3 amplexos observados el primer día, 2 el segundo día y así sucesivamente hasta 1 amplexo el quinto día. Para poder usar la función dura() definida anteriormente, transformé las secuencias obtenidas en secuencias binarias. El ejemplo 32201 entonces lo transformé a 11101 y usé la función dura() para estimar la duración del amplexo.

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(readxl)
library(gridExtra)

## Warning: package 'gridExtra' was built under R version 4.3.3

## 
## Attaching package: 'gridExtra'

## The following object is masked from 'package:dplyr':
## 
##     combine

duración <- read_excel("datos.xlsx", sheet = "Observaciones")

duración <- duración |>
  mutate(Amplx_obs = case_when(Amplx_obs >= 1 ~ 1,
                               Amplx_obs == 0 ~ 0))

#Lo siguiente que hice fue aplicar la función de la duración para cada bloque y tratamiento
#y las separé por observación.

dur_x_obs <- vector()

for (i in 1:3) { #Bucle para separar por observación
  for (j in 1:7) { #Bucle para separar por bloque
      for (k in 1:11) { #Bucle para separar por tratamiento
        ind <- duración$Obs == i & duración$Bloq == j & duración$Prop_inicial == k
        dur_x_obs <- c(dur_x_obs, dura(duración$Amplx_obs[ind]))
      }
  }
}

#Debido a que la información de la duración estimada es de un solo dato para el tratamiento de cada bloque, añadí está información al dataframe que contiene la base de datos en la hoja "Exito_amplx"
Amplexo <- read_excel("datos.xlsx", sheet = "Exito_amplx")
Amplexo <- Amplexo |> mutate(Dur_obs_1 = dur_x_obs[1:77], Dur_obs_2 = dur_x_obs[78:154],
                             Dur_obs_3 = dur_x_obs[155:231])

#Por último, grafiqué las duraciones con respecto a cada observación.

ggduración_1 <- ggplot(data = Amplexo, mapping = aes(Prop_final, Dur_obs_1)) +
  geom_point(alpha = 0.3) +
  geom_smooth(colour = "red") +
  labs(title = "Observación 1",
       x = "Proporción de sexos",
       y = "Duración del amplexo",
       caption = "n = 7")

ggduración_2 <- ggplot(data = Amplexo, mapping = aes(Prop_final, Dur_obs_2)) +
  geom_point(alpha = 0.3) +
  geom_smooth(colour = "red") +
  labs(title = "Observación 2",
       x = "Proporción de sexos",
       y = "Duración del amplexo",
       caption = "n = 7")

ggduración_3 <- ggplot(data = Amplexo, mapping = aes(Prop_final, Dur_obs_3)) +
  geom_point(alpha = 0.3) +
  geom_smooth(colour = "red") +
  labs(title = "Observación 3",
       x = "Proporción de sexos",
       y = "Duración del amplexo",
       caption = "n = 7")

grid.arrange(ggduración_1, ggduración_2, ggduración_3, ncol=2)

## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

Además de las estimaciones de la duración del amplexo para cada observación de forma individual, también hice la medición de forma global. En este caso, reuní la información de las tres observaciones para cada tratamiento por día y lo resumí en una sola secuencia.

#Para trabajar con los datos, resumí la información. Hice un nuevo dataframe
#en donde usé el promedio de crustáceos observados por día,
#el promedio de la temperatura por día y la presencia/ausecia de amplexos por día (independientemente de cuántos y en qué momento fueron encontrados.

#Primero defino los vectores y contadores que usaré para los bucles de escritura.
bloque <- vector()
día <- vector()
prop <- vector()
crust <- vector()
temp <- vector()
amplexos <- vector()



for(j in 1:7) { #Primero un bucle que representa el número de bloques
  
 for (k in 1:5) { #Un segundo bloque anidado que representa el número de días

   for (l in 1:11) { #Un tercer bucle anidado que representa la cantidad de tratamientos
     prop <- c(prop, l)
     ind <- duración$Bloq == j & duración$Dia == k & duración$Prop_inicial == l
     crust <- c(crust, mean(duración$Crust_obs[ind]))
     temp <- c(temp, max(duración$Temp[ind]))
     amplexos <- c(amplexos, mean(duración$Amplx_obs[ind]))
     día <- c(día, k)
     bloque <- c(bloque, j)
   }
  l <- 1
 }
k <- 1
}


resumen <- data.frame(bloque, día, prop, crust, temp, amplexos) #Lo añadí todo dentro de un dataframe


#Después, bajé la resolución de la medición al mínimo. Cualquier medición
#diferente de cero representará un valor de TRUE o 1 y los valores de 0 representarán False

resumen <- mutate(resumen, amplexos = amplexos != 0)

#Usé la función dura() en bucle para escribir un vector con las duraciones estimadas
#para cada tanque por bloque

Duración <- vector()

for (i in 1:7) {
  for (j in 1:11){
    ind <- resumen$bloque == i & resumen$prop == j
    x <- resumen$amplexos[ind]
    Duración <- c(Duración, dura(x))
  }
}

#Finalmente, grafiqué la duración del amplexo en función de la proporción de sexos.

ggduración <- ggplot(data = Amplexo, mapping = aes(Prop_final, Duración)) +
  geom_point(alpha = 0.3) +
  geom_smooth() +
  labs(title = "Duración del amplexo",
       x = "Proporción de sexos",
       y = "Duración del amplexo",
       caption = "n = 7")
ggduración

## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

Discusión

Las estimaciones difieren en la escala, pero el patrón es el mismo. En la siguiente gráfica se puede apreciar que la estimación de la duración del amplexo para las tres observaciones indivuduales no es diferente y que el patrón es similar tanto para las observaciones individuales como para las observaciones globales (por día). Sin embargo, las curvas de las observaciones individuales están por debajo de la curva de las observaciones globales dando como resultado una estimación de menor duración para el amplexo.

ggplot(data = Amplexo, aes(x = Prop_inicial)) +
  geom_smooth(aes(y = Dur_obs_1), colour = "red") +
  geom_smooth(aes(y = Dur_obs_2), colour = "red") +
  geom_smooth(aes(y = Dur_obs_3), colour = "red") +
  geom_smooth(aes(y = Duración)) +
  labs(title = "Comparación de las duraciones estimadas",
       x = "Proporción de sexos",
       y = "Duración del amplexo en días",
       caption = "n = 7")

## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

Personalmente, considero que el método de estimación usando la observación global es más adecuado porque considera toda la información disponible de presencia/ausencia de amplexos por día.

El método de estimación de la duración del amplexo depende de la cantidad de días observados. Algo de lo que me pude percatar, es que el método de medición depende del numero de días de la secuencia binaria de presencia/ausencia de amplexos. Los siguientes histogramas demuestran que al aumentar el número de días de las secuencias binarias, la estimación de la duración también aumenta al usar todas las permutaciones posibles.

y <- 5
subsets <- permutations(n= 2, r= y, v= c(1,0), repeats.allowed = TRUE)
cada_corrida <- totalruns(subsets)
dts <- round(simulacion(10000), digits = 2)
data <- as.data.frame(dts, rep("duración",10000))
b <- ggplot(data, aes(x = dts)) +
geom_bar() +
  labs(title = "5 días",
       x = "Duración",
       y = "Frecuencia",
       caption = "n = 10000")

a <- a + labs(title = "8 días",
       x = "Duración",
       y = "Frecuencia",
       caption = "n = 10000")

y <- 11
subsets <- permutations(n= 2, r= y, v= c(1,0), repeats.allowed = TRUE)
cada_corrida <- totalruns(subsets)
dts <- round(simulacion(10000), digits = 2)
data <- as.data.frame(dts, rep("duración",10000))
c <- ggplot(data, aes(x = dts)) +
geom_bar() +
    labs(title = "11 días",
       x = "Duración",
       y = "Frecuencia",
       caption = "n = 10000")

y <- 14
subsets <- permutations(n= 2, r= y, v= c(1,0), repeats.allowed = TRUE)
cada_corrida <- totalruns(subsets)
dts <- round(simulacion(10000), digits = 2)
data <- as.data.frame(dts, rep("duración",10000))
d <- ggplot(data, aes(x = dts)) +
geom_bar() +
    labs(title = "14 días",
       x = "Duración",
       y = "Frecuencia",
       caption = "n = 10000")

grid.arrange(b, a, c, d, ncol = 2)

En mi caso, solo registré los amplexos durante 5 días (porque así lo acordamos), a diferencia de el artículo publicado, en donde registraron los amplexos durante 8 días. Esto es una explicación posible de que las duraciones que yo estimé en mis experimentos sean menores en comparación a las reportadas en el artículo.

Conclusión

El método de medición que usé para medir la duración del amplexo me permitió estimar las duraciones para cada proporción de sexos al menos de forma relativa. Esto fue suficiente para observar que la duración del amplexo aumenta cuando la proporción de sexos se acerca a 0.5.