Modelos Gam aplicados a estudios de pesquerías

Los modelos GAM (Generalized Additive Models) suelen utilizarse en análisis de CPUE o abundancia en pesquerías para modelar relaciones no lineales entre la respuesta y las variables explicativas mediante funciones suaves (splines).

Un modelo GAM es una extensión del modelo lineal generalizado:

\[ g(\mathbb{E}[Y_i]) = \beta_0 + f_1(x_{i1}) + f_2(x_{i2}) + \dots + f_p(x_{ip}) \]

donde:

Cuando se una una función de enlace logarítmica, se modela la respuesta media en escala logarítmica, muy útil para datos positivos como la CPUE.

Ejemplo

Queremos modelar la CPUE (captura por unidad de esfuerzo) de una especie pelágica en función de:

  • Temperatura superficial del mar (SST)
  • Profundidad del muestreo

En este ejemplo no incluimos efectos aleatorios (como el año o el barco). Un efecto aleatorio es una variable de la que no nos interesa evaluar el efecto particular de cada valor, sino tener en cuenta su contribución a la variabilidad general observada en los datos. Por ejemplo, si el muestreo se realiza de forma tal que los registros de capturas son realizados por una muestra de barcos del total que realiza actividades pesqueras, normalmente no nos interesa estudiar las posibles diferencias entre barcos, sino descontar la variabilidad que se pueda haber introducido por utilizar barcos distintos.

Simulación y ajuste del modelo en R

Mostrar código R 
# Instalación de paquetes si es necesario
#install.packages(c("mgcv", "marginaleffects", "ggplot2", "dplyr"))

library(mgcv)
library(marginaleffects)
library(ggplot2)
library(dplyr)

set.seed(42)

# Simular datos
n <- 400
SST <- runif(n, 18, 28)        # temperatura superficial del mar (°C)
prof <- runif(n, 10, 200)      # profundidad (m)

# Efectos de la temperatura y profundidad que vamos a simular:
f1 <- function(x) sin(x/3) + 0.2*(x-23)^2/25   # efecto SST
f2 <- function(x) -0.004*x + 0.002*(x-100)^2/100 # efecto profundidad
  • Visualizamos las funciones suaves “verdaderas”:
Mostrar código R 
x <- seq(18, 28, length.out = 100)
plot(x, f1(x), type='l', main="Efecto no lineal de SST", ylab="f1(SST)", xlab="SST (°C)")

Mostrar código R 
x <- seq(10, 200, length.out = 100)
plot(x, f2(x), type='l', main="Efecto no lineal de profundidad", ylab="f2(prof)", xlab="Profundidad (m)")

  • Simulamos los datos:
Mostrar código R 
# Modelo verdadero (en log escala)
eta <- 2 + f1(SST) + f2(prof)
CPUE <- rlnorm(n, eta, 0.3)    # respuesta positiva y asimétrica

datos <- data.frame(CPUE, SST, prof)
  • Ajuste del GAM
Mostrar código R 
modelo_gam <- gam(
  log(CPUE) ~ s(SST, k=10) + s(prof, k=10),
  data = datos
)

summary(modelo_gam)

Family: gaussian 
Link function: identity 

Formula:
log(CPUE) ~ s(SST, k = 10) + s(prof, k = 10)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.26389    0.01485   152.4   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Approximate significance of smooth terms:
          edf Ref.df     F p-value    
s(SST)  3.507  4.345 97.64  <2e-16 ***
s(prof) 2.082  2.595 67.72  <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.624   Deviance explained =   63%
GCV = 0.089745  Scale est. = 0.088267  n = 400
  • Predicciones con marginaleffects

Creamos predicciones suavizadas manteniendo las demás variables constantes en su media. De esta forma podemos interpretar el modelo ajustado:

Mostrar código R 
# Predicciones para SST (manteniendo prof en su media)
ef_sst <- predictions(
  modelo_gam,
  variables = list(SST = seq(18, 28, length.out = 100)),
  newdata = datagrid(prof = mean(datos$prof))
)

# Predicciones para profundidad (manteniendo SST en su media)
ef_prof <- predictions(
  modelo_gam,
  variables = list(prof = seq(10, 200, length.out = 100)),
  newdata = datagrid(SST = mean(datos$SST))
)

(a) Efecto de la temperatura superficial (SST)

Mostrar código R 
ggplot(ef_sst, aes(x = SST, y = estimate)) +
  geom_point(data=datos, aes(x = SST, y = log(CPUE)),
             color = "gray60", alpha = 0.85) +
  geom_line(color = "darkblue", linewidth = 1.2) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high),
              fill = "skyblue", alpha = 0.3) +
  labs(
    title = "Efecto de la temperatura superficial (SST) sobre log(CPUE)",
    y = "log(CPUE) predicho",
    x = "Temperatura superficial (°C)"
  ) +
  theme_minimal(base_size = 13)

(b) Efecto de la profundidad

Mostrar código R 
ggplot(ef_prof, aes(x = prof, y = estimate)) +
  geom_point(data=datos, aes(x = prof, y = log(CPUE)),
             color = "gray60", alpha = 0.85) +
  geom_line(color = "darkgreen", linewidth = 1.2) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high),
              fill = "palegreen3", alpha = 0.3) +
  labs(
    title = "Efecto de la profundidad sobre log(CPUE)",
    y = "log(CPUE) predicho",
    x = "Profundidad (m)"
  ) +
  theme_minimal(base_size = 13)

Pendientes locales con slopes()

Podemos explorar cómo cambia la tasa de variación local de la CPUE respecto a cada covariable.

Mostrar código R 
# Pendiente local de SST
slope_sst <- slopes(
  modelo_gam,
  variables = "SST",
  newdata = datagrid(SST = seq(18, 28, length.out = 100),
                     prof = mean(datos$prof))
)

ggplot(slope_sst, aes(x = SST, y = estimate)) +
  geom_line(color = "firebrick", linewidth = 1.2) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high),
              fill = "salmon", alpha = 0.3) +
  labs(
    title = "Pendiente local: efecto de la SST sobre log(CPUE)",
    y = "d(log(CPUE)) / d(SST)",
    x = "Temperatura superficial (°C)"
  ) +
  theme_minimal(base_size = 13)

Interpretación en contexto pesquero

  • El modelo muestra cómo la CPUE varía suavemente con la temperatura y la profundidad, sin imponer una forma funcional específica (como lineal o cuadrática).

  • El pico en la curva de SST puede interpretarse como la temperatura óptima donde la especie es más abundante.

  • El descenso con la profundidad puede reflejar una preferencia por aguas someras.

  • El análisis de pendientes locales ayuda a detectar zonas de cambio rápido (por ejemplo, umbrales térmicos).

  • Estos resultados pueden informar estrategias de manejo pesquero, como definir áreas protegidas o temporadas de pesca basadas en condiciones ambientales favorables para la especie objetivo.

Diferencias regionales

En la práctica, la relación entre la CPUE y variables como la temperatura superficial (SST) o la profundidad puede variar según la región o zona de pesca, debido a diferencias ambientales o ecológicas.

Vamos a ver cómo tratar ese caso con GAM y marginaleffects, paso a paso.

Para ello ajustamos un modelo del tipo:

\[ \log(\text{CPUE}_i) = \beta_0 + f_1(\text{SST}_i, \text{región}_i) + f_2(\text{profundidad}_i, \text{región}_i) + \varepsilon_i \]

Es decir:

  • los efectos de SST y profundidad pueden cambiar según la región.

  • las funciones (f_1) y (f_2) incluyen suavizados distintos por región.

Esto se conoce en mgcv como el uso de “by-smooths”.

Simulación de datos

Creamos tres regiones (región A, B, C) con diferentes respuestas a la SST o la profundidad.

Mostrar código R 
library(mgcv)
library(marginaleffects)
library(ggplot2)
library(dplyr)

set.seed(101)

# Tamaño de muestra y regiones
n <- 600
region <- factor(rep(c("A", "B", "C"), each = n/3))

SST <- runif(n, 18, 28)
prof <- runif(n, 10, 200)

# Efectos diferentes por región
f_sst <- function(x, reg) {
  if (reg == "A") return(sin(x/3))
  if (reg == "B") return(0.8*sin(x/3 + 0.5))
  if (reg == "C") return(1.2*sin(x/3 - 0.3))
}

f_prof <- function(x, reg) {
  if (reg == "A") return(-0.004*x)
  if (reg == "B") return(-0.002*x + 0.001*(x-100)^2/100)
  if (reg == "C") return(-0.006*x + 0.002*(x-150)^2/150)
}

eta <- numeric(n)
for (i in 1:n) {
  eta[i] <- 2 + f_sst(SST[i], region[i]) + f_prof(prof[i], region[i])
}

CPUE <- rlnorm(n, eta, 0.3)
datos <- data.frame(CPUE, SST, prof, region)

Mostramos el efecto “real” por región:

Mostrar código R 
par(mfrow=c(1,2))
x <- seq(18, 28, length.out = 100)
plot(x, f_sst(x, "A"), type='l', col='blue', ylim=c(-1,1), ylab="f(SST)", xlab="SST ( °C )", main="Efecto de SST por región")
lines(x, f_sst(x, "B"), col='red')
lines(x, f_sst(x, "C"), col='green')
legend("bottomright", legend=c("Región A", "Región B", "Región C"), col=c("blue", "red", "green"), lty=1)

Mostrar código R 
x <- seq(10, 200, length.out = 100)
plot(x, f_prof(x, "A"), type='l', col='blue', ylim=c(-1,0.5), ylab="f(prof)", xlab="Profundidad (m)", main="Efecto de profundidad por región")
lines(x, f_prof(x, "B"), col='red')
lines(x, f_prof(x, "C"), col='green')
legend("topright", legend=c("Región A", "Región B", "Región C"), col=c("blue", "red", "green"), lty=1)

Ajuste del GAM con efectos por región

La sintaxis de mgcv usa el argumento by= dentro de s() para permitir efectos distintos por grupo.

Mostrar código R 
# Ajuste del GAM con efectos por región
modelo_reg <- gam(
  log(CPUE) ~ region + 
    s(SST, by = region, k = 10) + 
    s(prof, by = region, k = 10),
  data = datos
)

summary(modelo_reg)

Family: gaussian 
Link function: identity 

Formula:
log(CPUE) ~ region + s(SST, by = region, k = 10) + s(prof, by = region, 
    k = 10)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.14496    0.02074 103.407  < 2e-16 ***
regionB      0.16673    0.02945   5.662 2.37e-08 ***
regionC     -0.06215    0.02948  -2.108   0.0354 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Approximate significance of smooth terms:
                  edf Ref.df       F p-value    
s(SST):regionA  3.922  4.830  66.584 < 2e-16 ***
s(SST):regionB  7.652  8.536  28.438 < 2e-16 ***
s(SST):regionC  3.713  4.606 142.657 < 2e-16 ***
s(prof):regionA 1.000  1.000  91.149 < 2e-16 ***
s(prof):regionB 3.146  3.901   3.597 0.00616 ** 
s(prof):regionC 1.480  1.815 148.621 < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.754   Deviance explained = 76.3%
GCV = 0.088245  Scale est. = 0.084728  n = 600

¿Cómo saber si hay diferencias significativas entre regiones?

El summary(modelo_reg) muestra para cada región los términos suaves:

Mostrar código R 
summary(modelo_reg)

Family: gaussian 
Link function: identity 

Formula:
log(CPUE) ~ region + s(SST, by = region, k = 10) + s(prof, by = region, 
    k = 10)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.14496    0.02074 103.407  < 2e-16 ***
regionB      0.16673    0.02945   5.662 2.37e-08 ***
regionC     -0.06215    0.02948  -2.108   0.0354 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Approximate significance of smooth terms:
                  edf Ref.df       F p-value    
s(SST):regionA  3.922  4.830  66.584 < 2e-16 ***
s(SST):regionB  7.652  8.536  28.438 < 2e-16 ***
s(SST):regionC  3.713  4.606 142.657 < 2e-16 ***
s(prof):regionA 1.000  1.000  91.149 < 2e-16 ***
s(prof):regionB 3.146  3.901   3.597 0.00616 ** 
s(prof):regionC 1.480  1.815 148.621 < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.754   Deviance explained = 76.3%
GCV = 0.088245  Scale est. = 0.084728  n = 600

Todos los p-valores son significativos, por lo que concluimos que los datos contienen evidencia suficiente para asegurar que los efectos de SST o profundidad cambian según la región.

Otra opción más formal es usar un test de razón de verosimilitud entre:

  • un modelo que asume mismo efecto para todas las regiones, y

  • otro con efectos diferentes (como el que acabamos de ajustar).

Mostrar código R 
modelo_noreg <- gam(
  log(CPUE) ~ region + s(SST, k=10) + s(prof, k=10),
  data = datos
)
# Test de razón de verosimilitud
anova(modelo_noreg, modelo_reg, test = "Chisq")
Analysis of Deviance Table

Model 1: log(CPUE) ~ region + s(SST, k = 10) + s(prof, k = 10)
Model 2: log(CPUE) ~ region + s(SST, by = region, k = 10) + s(prof, by = region, 
    k = 10)
  Resid. Df Resid. Dev     Df Deviance  Pr(>Chi)    
1    591.34     92.702                              
2    572.31     48.811 19.026   43.891 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El modelo con efectos por región (modelo_reg) mejora significativamente (p < 0.05), por lo que hay evidencia de diferencias entre regiones en el efecto de SST o profundidad.

Para visualizar estas diferencias, usamos marginaleffects.

Visualización con marginaleffects

Podemos visualizar las curvas predichas por región:

  • Efecto de la temperatura:
Mostrar código R 
# Predicciones de log(CPUE) frente a SST por región
ef_sst <- predictions(
  modelo_reg,
  variables = list(SST = seq(18, 28, length.out = 100)),
  newdata = datagrid(prof = mean(datos$prof), region = unique(datos$region))
)

ggplot(ef_sst, aes(x = SST, y = estimate, color = region, fill = region)) +
  geom_point(data=datos, aes(x = SST, y = log(CPUE), color = region),
             alpha = 0.6, size=1) +
  geom_line(linewidth = 1.2) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = 0.2, color = NA) +
  labs(
    title = "Efecto de la temperatura superficial (SST) sobre log(CPUE) por región",
    y = "log(CPUE) predicho",
    x = "Temperatura superficial (°C)"
  ) +
  theme_minimal(base_size = 13)

  • Efecto de la profundidad:
Mostrar código R 
ef_prof <- predictions(
  modelo_reg,
  variables = list(prof = seq(10, 200, length.out = 100)),
  newdata = datagrid(SST = mean(datos$SST), region = unique(datos$region))
)

ggplot(ef_prof, aes(x = prof, y = estimate, color = region, fill = region)) +
  geom_point(data=datos, aes(x = prof, y = log(CPUE), color = region),
             alpha = 0.6, size=1) +
  geom_line(linewidth = 1.2) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = 0.2, color = NA) +
  labs(
    title = "Efecto de la profundidad sobre log(CPUE) por región",
    y = "log(CPUE) predicho",
    x = "Profundidad (m)"
  ) +
  theme_minimal(base_size = 13)

####Comparación visual entre regiones

Para comparar directamente diferencias de predicciones entre regiones (por ejemplo, región B − región A) para cada valor de SST podemos usar el paquete gratia:

Mostrar código R 
library(gratia)
diff_sst <- difference_smooths(
  modelo_reg,
  select = "s(SST)",
  pairwise = TRUE
)

# Visualizar
draw(diff_sst) +
  ggtitle("Diferencias entre suavizados de SST por región (con IC 95%)") +
  theme_minimal(base_size = 13)

Mostrar código R 
library(gratia)
diff_prof <- difference_smooths(
  modelo_reg,
  select = "s(prof)",
  pairwise = TRUE
)

# Visualizar
draw(diff_prof) +
  ggtitle("Diferencias entre suavizados de Prof por región (con IC 95%)") +
  theme_minimal(base_size = 13)

Modelo con efectos aleatorios.

Partimos del modelo anterior, añadiendo ahora un factor barco como efecto aleatorio:

Mostrar código R 
library(mgcv)
library(marginaleffects)
library(ggplot2)
library(dplyr)

set.seed(2025)

# Estructura del experimento
n <- 1200
region <- factor(rep(c("A", "B", "C"), each = n/3))
barco <- factor(rep(paste0("B", 1:20), each = n/20))

SST <- runif(n, 18, 28)
prof <- runif(n, 10, 200)

# Funciones suaves distintas por región
f_sst <- function(x, reg) {
  if (reg == "A") return(sin(x/3))
  if (reg == "B") return(0.8*sin(x/3 + 0.5))
  if (reg == "C") return(1.2*sin(x/3 - 0.3))
}

f_prof <- function(x, reg) {
  if (reg == "A") return(-0.004*x)
  if (reg == "B") return(-0.002*x + 0.001*(x-100)^2/100)
  if (reg == "C") return(-0.006*x + 0.002*(x-150)^2/150)
}

# Efectos aleatorios
b_barco <- rnorm(length(levels(barco)), 0, 0.4)

eta <- numeric(n)
for (i in 1:n) {
  eta[i] <- 2 +
    f_sst(SST[i], region[i]) +
    f_prof(prof[i], region[i]) +
    b_barco[barco[i]]
}

CPUE <- rlnorm(n, eta, 0.3)
datos <- data.frame(CPUE, SST, prof, region, barco)

El modelo completo se ajusta así:

Mostrar código R 
modelo_gamm2 <- gamm(
  log(CPUE) ~ region +
    s(SST, by = region, k = 10) +
    s(prof, by = region, k = 10),
  random = list(
      barco = ~1
  ),
  data = datos
)

summary(modelo_gamm2$gam)

Family: gaussian 
Link function: identity 

Formula:
log(CPUE) ~ region + s(SST, by = region, k = 10) + s(prof, by = region, 
    k = 10)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.24410    0.10574  21.223   <2e-16 ***
regionB      0.13146    0.07812   1.683   0.0927 .  
regionC     -0.13550    0.10541  -1.285   0.1989    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Approximate significance of smooth terms:
                  edf Ref.df      F p-value    
s(SST):regionA  5.180  5.180 105.77  <2e-16 ***
s(SST):regionB  4.963  4.963 101.19  <2e-16 ***
s(SST):regionC  6.194  6.194 198.25  <2e-16 ***
s(prof):regionA 1.894  1.894 109.30  <2e-16 ***
s(prof):regionB 2.219  2.219  16.13  <2e-16 ***
s(prof):regionC 2.175  2.175 271.74  <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.497   
  Scale est. = 0.092897  n = 1200
  • Evaluación de las diferencias regionales:
Mostrar código R 
modelo_comun <- gamm(
  log(CPUE) ~ region + s(SST, k = 10) + s(prof, k = 10),
  random = list( barco = ~1),
  data = datos
)

anova(modelo_comun$gam)

Family: gaussian 
Link function: identity 

Formula:
log(CPUE) ~ region + s(SST, k = 10) + s(prof, k = 10)

Parametric Terms:
       df     F p-value
region  2 5.011 0.00681

Approximate significance of smooth terms:
          edf Ref.df     F p-value
s(SST)  5.595  5.595 151.7  <2e-16
s(prof) 2.076  2.076 196.2  <2e-16