Modelado en R

B0305 – Lab Eco Gral
https://ucr-ecologia.github.io/B0305-lab-eco-gral/

Prof. G. Avalos y T. Nakov

Escuela de Biologia, Universidad de Costa Rica

01 Aug 2025

Que vamos a ver:

• Introducción al Modelado
• Modelo 1: Regresión Lineal
• Modelo 2: Regresión Log-Lineal
• Modelo 3: Crecimiento Logístico
• Evaluación de Modelos
• Interpretación Biológica
• Conclusiones

Introducción al Modelado

¿Qué es un modelo?

Definición

Un modelo es una representación simplificada de un proceso o sistema real que nos ayuda a:

• Entender patrones en los datos

• Predecir valores futuros

• Probar hipótesis sobre mecanismos biológicos

Ejemplo: Crecimiento de levaduras

Imaginemos que cultivamos levaduras en el laboratorio y medimos su densidad poblacional a lo largo del tiempo…

¿Qué tipo de patrón esperaríamos ver? 🤔

Patrón típico del crecimiento microbiano

Fase lag (retraso)

Primeras horas: Las células se adaptan al medio, casi no hay crecimiento

Fase exponencial (logarítmica)

Crecimiento rápido: Recursos abundantes, las células se dividen a máxima velocidad

• Duplicación constante: 1 → 2 → 4 → 8 → 16…

• ¡Patrón exponencial! \(N(t) = N_0 \cdot e^{rt}\)

Fase estacionaria

Recursos limitados: Nutrientes se agotan, se acumulan toxinas

• Capacidad de carga: Densidad máxima que el medio puede sostener

• Crecimiento se desacelera hasta detenerse

Resultado: Curva sigmoidea (forma de S)

¿Qué modelo matemático describe mejor este patrón? 🤔

La ecuación de crecimiento exponencial

Durante la fase exponencial, el crecimiento se describe con:
\[N(t) = N_0 \cdot e^{rt}\]

Donde:
- \(N(t)\) = densidad poblacional en el tiempo \(t\)
- \(N_0\) = densidad poblacional inicial (en \(t = 0\))
- \(r\) = tasa intrínseca de crecimiento (por unidad de tiempo)
- \(e\) = base del logaritmo natural (≈ 2.718)
- \(t\) = tiempo

Interpretación biológica:
- Si \(r > 0\): población crece exponencialmente
- Si \(r = 0\): población constante
- Si \(r < 0\): población decrece exponencialmente

Problema: Que pasa cuando se agotan los recursos?

Aggregando capacidad de carga

Problema del modelo exponencial: ¡Crecimiento infinito es imposible!

En la realidad, los recursos son limitados:
- Espacio físico finito
- Nutrientes se agotan
- Acumulación de desechos tóxicos
- Competencia entre individuos

Solución: Introducir la capacidad de carga (\(K\))

\[N(t) = \frac{K}{1 + \left(\frac{K-N_0}{N_0}\right)e^{-rt}}\]

Parámetros del modelo logístico:
- \(K\) = capacidad de carga (densidad máxima sostenible)
- \(r\) = tasa intrínseca de crecimiento (cuando N << K)
- \(N_0\) = densidad inicial

Comportamiento del modelo:
- Cuando \(N \ll K\): crecimiento ≈ exponencial
- Cuando \(N \approx K/2\): máxima tasa de crecimiento
- Cuando \(N \to K\): crecimiento se desacelera hacia cero

Datos de crecimiento de levaduras

Simulemos datos realistas de crecimiento de Saccharomyces cerevisiae:

# Crear datos de crecimiento de levaduras
set.seed(123)
tiempo <- seq(0, 24, by = 2)  # horas
# Crecimiento logístico con ruido
K <- 1000  # capacidad de carga
r <- 0.3   # tasa de crecimiento
N0 <- 10   # población inicial

densidad_real <- K / (1 + ((K - N0) / N0) * exp(-r * tiempo))
# Añadir ruido realista
densidad <- densidad_real * exp(rnorm(length(tiempo), 0, 0.1))

# Crear data frame
levaduras <- data.frame(tiempo = tiempo, densidad = densidad)
head(levaduras, 8)

  tiempo   densidad
1      0   9.454942
2      2  17.661370
3      4  37.921560
4      6  57.995942
5      8 101.493323
6     10 200.220597
7     12 282.633370
8     14 354.657568

Visualizar los datos

library(ggplot2)

# Configurar tema global para todos los gráficos
theme_set(theme_bw())

p1 <- ggplot(levaduras, aes(x = tiempo, y = densidad)) +
  geom_point(size = 3, color = "steelblue") +
  labs(x = "Tiempo (horas)", 
       y = "Densidad (células/mL)",
       title = "Crecimiento de Saccharomyces cerevisiae")

# Mostrar el gráfico
p1

¿Qué patrón observas? ¿Cómo modelarías este crecimiento? 🤔

Modelo 1: Regresión Lineal

Intentemos un modelo lineal

Hipótesis ingenua: La densidad aumenta linealmente con el tiempo

# Ajustar modelo lineal
modelo_lineal <- lm(densidad ~ tiempo, data = levaduras)
summary(modelo_lineal)

Call:
lm(formula = densidad ~ tiempo, data = levaduras)

Residuals:
     Min       1Q   Median       3Q      Max 
-122.785  -95.839    2.071   73.387  166.855 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -157.400     54.507  -2.888   0.0148 *  
tiempo        45.346      3.854  11.765 1.42e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 104 on 11 degrees of freedom
Multiple R-squared:  0.9264,    Adjusted R-squared:  0.9197 
F-statistic: 138.4 on 1 and 11 DF,  p-value: 1.425e-07

Visualizar el ajuste lineal

p2 <- ggplot(levaduras, aes(x = tiempo, y = densidad)) +
  geom_point(size = 3, color = "steelblue") +
  geom_smooth(method = "lm", se = FALSE, color = "red", size = 1.2) +
  labs(x = "Tiempo (horas)", 
       y = "Densidad (células/mL)",
       title = "Modelo Lineal: ¿Buen ajuste?")

p2

¿Qué problemas observas con este modelo? 🚨

Diagnósticos del modelo lineal

# Gráficos de diagnóstico
par(mfrow = c(1, 2))
plot(modelo_lineal, which = 1:2)

predicciones <- predict(modelo_lineal)
tiempo_negativo <- levaduras$tiempo[predicciones < 0]
cat(
  "Tiempos con predicciones negativas:", 
  tiempo_negativo
)

Tiempos con predicciones negativas: 0 2

Problemas evidentes:
- Residuos no aleatorios (patrón curvo)
- Heteroscedasticidad (varianza no constante)
- ¡El modelo predice densidades negativas!

Modelo 2: Regresión Log-Lineal

Transformación logarítmica

Nueva hipótesis: El crecimiento es exponencial (al menos inicialmente)

Si: \(N(t) = N_0 \cdot e^{rt}\)

Entonces: \(\log(N(t)) = \log(N_0) + rt\)

# Crear variable transformada
levaduras$log_densidad <- log(levaduras$densidad)

# Ajustar modelo log-lineal
modelo_log <- lm(log_densidad ~ tiempo, data = levaduras)
summary(modelo_log)

Call:
lm(formula = log_densidad ~ tiempo, data = levaduras)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.66649 -0.33371  0.05152  0.26777  0.49981 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.84004    0.20683   13.73 2.88e-08 ***
tiempo       0.19596    0.01462   13.40 3.72e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.3946 on 11 degrees of freedom
Multiple R-squared:  0.9423,    Adjusted R-squared:  0.937 
F-statistic: 179.5 on 1 and 11 DF,  p-value: 3.716e-08

Visualizar en escala logarítmica

p3 <- ggplot(levaduras, aes(x = tiempo, y = log_densidad)) +
  geom_point(size = 3, color = "darkgreen") +
  geom_smooth(method = "lm", se = FALSE, color = "red", size = 1.2) +
  labs(x = "Tiempo (horas)", 
       y = "log(Densidad)",
       title = "Modelo Log-Lineal")

p3

¿Mejor ajuste en los primeros puntos, pero qué pasa después? 🤔

Comparar modelos en escala original

# Predicciones del modelo log-lineal
levaduras$pred_log <- exp(predict(modelo_log))

p4 <- ggplot(levaduras, aes(x = tiempo, y = densidad)) +
  geom_point(size = 3, color = "steelblue") +
  geom_line(aes(y = pred_log), color = "darkgreen", size = 1.2) +
  labs(x = "Tiempo (horas)", 
       y = "Densidad (células/mL)",
       title = "Modelo Log-Lineal vs Datos Reales")

p4

Problema: ¡Predice crecimiento exponencial infinito! 📈

Modelo 3: Crecimiento Logístico

El modelo biológicamente realista

Hipótesis biológica: Las poblaciones tienen una capacidad de carga limitada

\[N(t) = \frac{K}{1 + \left(\frac{K-N_0}{N_0}\right)e^{-rt}}\]

Donde:
- \(K\) = capacidad de carga (densidad máxima)
- \(r\) = tasa intrínseca de crecimiento
- \(N_0\) = densidad inicial

# Ajustar modelo logístico usando nls()
# Non-linear least squares 
modelo_logistico <- nls(
  densidad ~ K / (1 + ((K - N0) / N0) * exp(-r * tiempo)),
  data = levaduras,
  start = list(K = 800, r = 0.2, N0 = 15)  # valores iniciales
)

summary(modelo_logistico)

Formula: densidad ~ K/(1 + ((K - N0)/N0) * exp(-r * tiempo))

Parameters:
    Estimate Std. Error t value Pr(>|t|)    
K  1.107e+03  6.432e+01  17.213 9.26e-09 ***
r  2.794e-01  2.696e-02  10.362 1.15e-06 ***
N0 1.185e+01  4.272e+00   2.775   0.0196 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 34.83 on 10 degrees of freedom

Number of iterations to convergence: 8 
Achieved convergence tolerance: 3.168e-06

Visualizar el modelo logístico

levaduras$pred_logistico <- predict(modelo_logistico)

p5 <- ggplot(levaduras, aes(x = tiempo, y = densidad)) +
  geom_point(size = 3, color = "steelblue") +
  geom_line(aes(y = pred_logistico), color = "purple", size = 1.2) +
  labs(x = "Tiempo (horas)", 
       y = "Densidad (células/mL)",
       title = "Modelo Logístico: ¡Ajuste Biológicamente Realista!")

p5

Comparación de los tres modelos

# Preparar datos para comparación
levaduras$pred_lineal <- predict(modelo_lineal)
levaduras$pred_logistico <- predict(modelo_logistico)

p6 <- ggplot(levaduras, aes(x = tiempo, y = densidad)) +
  geom_point(size = 3, color = "steelblue") +
  geom_line(aes(y = pred_lineal, color = "Lineal"), 
            size = 1) +
  geom_line(aes(y = pred_log, color = "Log-lineal"), 
            size = 1) +
  geom_line(aes(y = pred_logistico, color = "Logístico"), 
            size = 1) +
  scale_color_manual(values = c("Lineal" = "red", 
                               "Log-lineal" = "darkgreen", 
                               "Logístico" = "purple")) +
  labs(x = "Tiempo (horas)", 
       y = "Densidad (células/mL)",
       title = "Comparación de Modelos",
       color = "Modelo")

p6

Evaluación de Modelos

Criterios de evaluación

¿Cómo decidimos cuál modelo es mejor?

1. R² (coeficiente de determinación)
2. Análisis de residuos
3. Interpretabilidad biológica
4. AIC/BIC (Otras opciones mas avanzadas…)

# Calcular R² para cada modelo
r2_lineal <- summary(modelo_lineal)$r.squared
r2_log <- summary(modelo_log)$r.squared

# Para modelo logístico, tenemos que calcular R2 manualmente
# ss_total = suma total de cuadrados (variación total en los datos)
ss_total <- sum((levaduras$densidad - mean(levaduras$densidad))^2)
# ss_residual = suma de cuadrados de residuos (variación no explicada)
ss_residual <- sum(residuals(modelo_logistico)^2)
# R2 = 1 - (variación no explicada / variación total)
r2_logistico <- 1 - (ss_residual / ss_total)

data.frame(
  Modelo = c("Lineal", "Log-lineal", "Logístico"),
  R_cuadrado = round(c(r2_lineal, r2_log, r2_logistico), 3)
)

      Modelo R_cuadrado
1     Lineal      0.926
2 Log-lineal      0.942
3  Logístico      0.992

Análisis de residuos

# Crear data frame con residuos
residuos_df <- data.frame(
  tiempo = rep(levaduras$tiempo, 3),
  residuos = c(residuals(modelo_lineal), 
               residuals(modelo_log), 
               residuals(modelo_logistico)),
  modelo = rep(c("Lineal", "Log-lineal", "Logístico"), 
               each = nrow(levaduras))
)

p7 <- ggplot(residuos_df) +
  aes(x = tiempo, y = residuos, color = modelo) +
  geom_point(size = 2) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  facet_wrap(~modelo, scales = "free_y") +
  labs(
    x = "Tiempo (horas)", 
    y = "Residuos",
    title = "Análisis de Residuos por Modelo") +
  theme(legend.position = "bottom")

p7

¿Qué modelo tiene los residuos más aleatorios? 🎯

⚠️ Cuidado con los residuos transformados

• Los residuos del modelo log-lineal están en escala logarítmica
• No son directamente comparables con residuos en escala original
• Para comparar correctamente, necesitamos calcular residuos en la misma escala

# Residuos en escala original para todos los modelos
residuos_originales <- data.frame(
  tiempo = levaduras$tiempo,
  lineal = levaduras$densidad - predict(modelo_lineal),
  log_lineal = levaduras$densidad - exp(predict(modelo_log)),  # ¡Transformar de vuelta!
  logistico = levaduras$densidad - predict(modelo_logistico)
)

# Calcular suma de cuadrados de residuos en escala original
rss_original <- sapply(residuos_originales[,-1], function(x) sum(x^2))
print(round(rss_original, 2))

    lineal log_lineal  logistico 
 118961.09 1022167.84   12133.29 

Ahora sí podemos comparar correctamente ✅

Residuos en escala original - Comparación visual

Ahora podemos ver cuál modelo realmente ajusta mejor 👀

# Crear data frame para visualización
residuos_originales_plot <- data.frame(
  tiempo = rep(levaduras$tiempo, 3),
  residuos = c(residuos_originales$lineal, 
               residuos_originales$log_lineal, 
               residuos_originales$logistico),
  modelo = rep(c("Lineal", "Log-lineal", "Logístico"), 
               each = nrow(levaduras))
)

p8 <- ggplot(residuos_originales_plot) +
  aes(x = tiempo, y = residuos, color = modelo) +
  geom_point(size = 2) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  facet_wrap(~modelo) +
  labs(x = "Tiempo (horas)", 
       y = "Residuos (escala original)",
       title = "Residuos en Escala Original - Comparación Justa"
    ) +
  theme(legend.position = "bottom")

p8

Residuos en escala original - Comparación visual

Interpretación Biológica

Parámetros del modelo logístico

Tasa, Capacidad de carga, densidad \(t0\)

# Extraer parámetros estimados
params <- coef(modelo_logistico)
print(params)

          K           r          N0 
1107.075934    0.279384   11.853795 

Intervalos de confianza

# Calcular intervalos de confianza
conf_int <- confint(modelo_logistico)
print(conf_int)

          2.5%        97.5%
K  994.6076225 1283.0776710
r    0.2264351    0.3434993
N0   4.8573373   23.8982205

Interpretación biológica:

• K ≈ 1107 células/mL: Capacidad de carga del medio
  
• r ≈ 0.279 h⁻¹: Tasa intrínseca de crecimiento
  
• N₀ ≈ 11.9 células/mL: Densidad inicial
  

Parametros de simulacion:

• K = 1000 células/mL
• r = 0.3 h⁻¹
• N₀ = 10 células/mL

Tiempo de duplicación

# Calcular tiempo de duplicación
r_estimado <- params["r"]
tiempo_duplicacion <- log(2) / r_estimado

cat(
  "Tiempo de duplicación:", 
  round(tiempo_duplicacion, 2), "horas\n"
)

Tiempo de duplicación: 2.48 horas

# Tiempo para alcanzar 50% de la capacidad de carga
K_estimado <- params["K"]
N0_estimado <- params["N0"]
# Calcular la constanta logistica
A <- log((K_estimado - N0_estimado) / N0_estimado)
t_50_K <- A / r_estimado

cat(
  "Tiempo para alcanzar 50% de K:", 
  round(t_50_K, 2), "horas\n"
)

Tiempo para alcanzar 50% de K: 16.2 horas

Estas métricas son importantes para:
- Optimización de cultivos
- Predicción de rendimientos
- Comparación entre cepas

📚 Leer más:

El punto de inflexión de la curva logística ocurre exactamente en \(K/2\). La fórmula \(t_{50\%} = \frac{\ln((K-N_0)/N_0)}{r}\) se deriva resolviendo la ecuación logística para \(N(t) = K/2\).

Referencias recomendadas:
- Gotelli, N.J. (2008). A Primer of Ecology, 4th ed. Sinauer Associates. (Capítulo 2: Exponential and Logistic Population Growth)

Nota importante: La logística nunca alcanza exactamente \(K\), solo se aproxima asintóticamente.

Conclusiones

Lecciones aprendidas

1. Los modelos simples (lineal) pueden ser inadecuados para procesos biológicos complejos

2. Las transformaciones (log-lineal) pueden mejorar el ajuste pero no capturan toda la biología

3. Los modelos mecanísticos (logístico) incorporan conocimiento biológico y proporcionan parámetros interpretables

4. La evaluación de modelos requiere múltiples criterios, no solo R²

5. El contexto biológico es crucial para la selección del modelo apropiado

Próximos pasos

En sus proyectos de investigación, consideren:

• ¿Qué proceso biológico están estudiando?
• ¿Qué tipo de relación esperan teóricamente?
• ¿Cómo pueden validar sus modelos?
• ¿Los parámetros estimados tienen sentido biológico?

Bibliografia

Gotelli, N.J. (2008). A Primer of Ecology, 4th ed. Sinauer Associates. (Capítulo 2: Exponential and Logistic Population Growth)