Introducción

El presente análisis se basa en un conjunto de datos de conejos europeos (Oryctolagus Cuniculus) recopilados por Dudzinski & Mykytowycz (1961) y posteriormente estudiados por Ratkowsky (1983) y Wei (1998). El objetivo de este estudio es modelar el peso de las lentes de los ojos de los conejos en función de su edad (medida en días).

La relación entre el peso de las lentes oculares y la edad de los conejos presenta un patrón no lineal, donde el peso aumenta rápidamente en las primeras etapas de vida, pero tiende a estabilizarse a medida que los conejos envejecen, alrededor de los 400 días. Esta estabilización es un fenómeno común en muchos procesos biológicos.

En este ejercicio, ajustaremos varios modelos, incluyendo un modelo no lineal inverso y diferentes modelos lineales generalizados (GLM), para identificar cuál de ellos captura mejor la relación entre la edad y el peso de las lentes oculares. Compararemos estos modelos utilizando el Criterio de Información de Akaike (AIC), que nos permite evaluar la calidad del ajuste de los modelos penalizando su complejidad.

Datos

Vamos a utilizar los datos que contienen la edad de los conejos y el peso de las lentes oculares.

# Datos de edad y peso
edad <- c(15, 15, 15, 18, 28, 29, 37, 37, 44, 50, 50, 60, 61, 64, 65, 65, 72, 75, 75, 82, 85, 91, 91, 97, 98, 125, 142, 142, 147, 150, 159, 165, 183, 192, 195, 218, 219, 224, 225, 227, 232, 237, 246, 258, 276, 285, 300, 301, 305, 312, 317, 338, 347, 354, 357, 375, 394, 513, 535, 554, 591, 648, 660, 705, 723, 756, 768, 860)
peso <- c(21.66, 22.75, 22.30, 31.25, 44.79, 40.55, 50.25, 46.88, 52.03, 63.47, 61.13, 81.00, 73.09, 79.09, 79.51, 65.31, 71.90, 86.10, 94.60, 92.50, 105.00, 101.70, 102.90, 110.00, 104.30, 134.90, 130.68, 140.58, 155.30, 144.50, 142.15, 139.81, 153.22, 145.72, 161.10, 174.18, 173.54, 178.86, 177.68, 173.73, 159.98, 187.07, 176.13, 183.40, 186.26, 189.66, 186.09, 186.70, 186.80, 195.10, 216.41, 203.23, 188.38, 189.70, 195.31, 202.63, 224.82, 203.30, 209.70, 233.90, 234.70, 244.30, 231.00, 242.40, 230.77, 242.57, 232.12, 246.70)

El gráfico muestra la relación entre el peso de las lentes oculares (en miligramos) y la edad de los conejos (en días). Se observa una clara tendencia no lineal, donde el peso de las lentes aumenta rápidamente durante los primeros 200 días de vida de los conejos, y posteriormente tiende a estabilizarse alrededor de los 250 mg a partir de los 400 días. Este patrón es consistente con muchos fenómenos biológicos, donde el crecimiento inicial es rápido y se desacelera conforme el organismo alcanza la madurez.

Modelo No Lineal Inverso

El modelo no lineal inverso utilizado en este análisis sigue la estructura propuesta por Wei (1998), este modelo no lineal busca capturar la relación asintótica observada en los datos: el peso de las lentes aumenta rápidamente en los primeros días de vida de los conejos y se estabiliza a medida que los conejos envejecen, alrededor de los 400 días. Esta estabilización del crecimiento es una característica común en procesos biológicos.

El modelo no lineal inverso utilizado es de la siguiente forma:

\[ \mu_i = \beta_1 - \frac{\beta_2}{x_i + \beta_3} \]

donde:

Este modelo asume que los pesos \(Y_i\) siguen una distribución normal inversa con media \(\mu_i\) y varianza \(\sigma^2\), es decir:

\[ Y_i \sim NI(\mu_i, \sigma^2) \]

El ajuste de este modelo se realiza utilizando el método de mínimos cuadrados no lineales, buscando los parámetros \(\beta_1\), \(\beta_2\) y \(\beta_3\) que mejor describan la relación entre la edad y el peso de las lentes oculares de los conejos.

Además, la varianza del peso de las lentes está relacionada con su media a través de la siguiente ecuación:

\[ \text{Var}(Y_i) = \sigma^2 V(\mu_i) \quad \text{con} \quad V(\mu_i) = \mu_i^3 \]

Esta expresión refleja cómo la varianza crece junto con la media cúbica del peso, lo que es típico en datos biológicos donde se observa una mayor variabilidad en pesos más altos. Esto también implica que el modelo es capaz de capturar la heterogeneidad de la varianza de los datos, una característica importante en los datos observados.

En resumen, el modelo asume que los pesos \(\mu_i\) son valores esperados que siguen una relación no lineal inversa en función de la edad, mientras que la varianza crece cúbicamente con el peso. Esta estructura permite modelar adecuadamente la tendencia asintótica y la variabilidad observada en los datos.

Resultados del Ajuste del Modelo

El ajuste del modelo no lineal se ha llevado a cabo utilizando la función nlsLM() en R, que implementa el algoritmo de Levenberg-Marquardt. A continuación, se presenta un resumen de los parámetros estimados para este modelo:

# Ajuste del modelo no lineal inverso
modelo_nl <- nlsLM(peso ~ beta1 - beta2 / (edad + beta3), start = list(beta1 = 250, beta2 = 10000, beta3 = 10))

El resumen del modelo no lineal ajustado muestra los siguientes parámetros estimados:

  • \(\beta_1 = 290.87\): Peso asintótico de las lentes oculares.
  • \(\beta_2 = 45685.80\): Controla la velocidad de crecimiento del peso.
  • \(\beta_3 = 152.98\): Edad a partir de la cual el peso comienza a estabilizarse.

El modelo indica que el peso de las lentes oculares de los conejos aumenta rápidamente en los primeros días de vida y luego se estabiliza alrededor de los 290.87 mg a partir de los 153 días. Todos los parámetros son altamente significativos (\(p < 2e-16\)), lo que respalda la calidad del ajuste del modelo.

# Resumen del modelo ajustado
summary(modelo_nl)
## 
## Formula: peso ~ beta1 - beta2/(edad + beta3)
## 
## Parameters:
##        Estimate Std. Error t value Pr(>|t|)    
## beta1   290.868      5.431   53.55   <2e-16 ***
## beta2 45685.798   3553.435   12.86   <2e-16 ***
## beta3   152.984     11.179   13.69   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.1 on 65 degrees of freedom
## 
## Number of iterations to convergence: 6 
## Achieved convergence tolerance: 1.49e-08

Diagnóstico de Residuos

Generamos el gráfico de residuos y residuos estandarizados para diagnosticar el ajuste del modelo.

# Gráfico de residuos crudos
plot(fitted(modelo_nl), resid(modelo_nl), main = "Diagnóstico de Residuos", xlab = "Valores Ajustados", ylab = "Residuos")
abline(h = 0, col = "red")

El gráfico de residuos muestra que los residuos están distribuidos en torno a la línea de cero, lo cual es positivo, ya que indica que el modelo no presenta sesgo sistemático. Sin embargo, hay algunos aspectos críticos a tener en cuenta:

  1. Varianza no constante: La dispersión de los residuos parece aumentar para los valores ajustados más altos, lo que indica heterocedasticidad en los datos, algo de esperarse debido a la función de varianza. Esto significa que la variabilidad de los errores no es constante a lo largo del rango de valores ajustados.

  2. Outliers: Se observan algunos residuos extremos alejados del resto de los puntos.

# Residuos estandarizados
sigma <- summary(modelo_nl)$sigma
residuos_estandarizados <- resid(modelo_nl) / sigma
plot(fitted(modelo_nl), residuos_estandarizados, main = "Residuos Estandarizados", xlab = "Valores Ajustados", ylab = "Residuos Estandarizados")
abline(h = 0, col = "red")
abline(h = c(-2, 2), col = "blue", lty = 2)

Gráfico del Ajuste del Modelo No Lineal

Graficamos el ajuste del modelo no lineal sobre los datos.

# Predicción usando el modelo ajustado
edad_nueva <- seq(min(edad), max(edad), length.out = 100)
predicciones <- predict(modelo_nl, newdata = data.frame(edad = edad_nueva))

# Gráfico de los datos observados y la curva ajustada
plot(edad, peso, main = "Peso vs Edad", xlab = "Edad (días)", ylab = "Peso (mg)", pch = 19)
lines(edad_nueva, predicciones, col = "#800080", lwd = 2)

GLM con Normal Inversa

Ajustamos un modelo lineal generalizado (GLM) con distribución normal inversa para ver si captura mejor la relación entre las variables.

# Ajuste del modelo GLM con distribución normal inversa
modelo_glm_inv <- glm(peso ~ edad, family = inverse.gaussian(link = "inverse"))

# Predicción y gráfico
predicciones_glm_inv <- predict(modelo_glm_inv, newdata = data.frame(edad = edad_nueva), type = "response")
plot(edad, peso, main = "Peso vs Edad (Modelo Normal Inversa)", xlab = "Edad (días)", ylab = "Peso (mg)", pch = 19)
lines(edad_nueva, predicciones_glm_inv, col = "#800080", lwd = 2)

Comparación de AIC entre Modelos

Compararemos los modelos utilizando el AIC para determinar cuál se ajusta mejor.

# Comparación de AIC
aic_nl <- AIC(modelo_nl)
aic_glm_inv <- AIC(modelo_glm_inv)

cat("AIC del Modelo No Lineal:", aic_nl, "\n")
## AIC del Modelo No Lineal: 482.3952
cat("AIC del Modelo Lineal con Normal Inversa:", aic_glm_inv, "\n")
## AIC del Modelo Lineal con Normal Inversa: 771.9755

Otros Modelos GLM

También ajustamos otros dos modelos lineales generalizados (GLM) con distribuciones Gamma y Log-Normal para comparar el ajuste.

Modelo GLM con Gamma

# Ajuste del modelo GLM con distribución Gamma y enlace logarítmico
modelo_glm_gamma <- glm(peso ~ edad, family = Gamma(link = "log"))
predicciones_glm_gamma <- predict(modelo_glm_gamma, newdata = data.frame(edad = edad_nueva), type = "response")

# Gráfico del ajuste
plot(edad, peso, main = "Peso vs Edad (Modelo Gamma con Enlace Log)", xlab = "Edad (días)", ylab = "Peso (mg)", pch = 19)
lines(edad_nueva, predicciones_glm_gamma, col = "#800080", lwd = 2)

Comparación de AIC para Modelo Gamma

aic_glm_gamma <- AIC(modelo_glm_gamma)
cat("AIC del Modelo GLM Gamma:", aic_glm_gamma, "\n")
## AIC del Modelo GLM Gamma: 725.2709

Modelo GLM con Log-Normal

# Ajuste del modelo GLM con distribución Log-Normal
modelo_glm_lognormal <- glm(log(peso) ~ edad, family = gaussian(link = "identity"))
predicciones_glm_lognormal <- exp(predict(modelo_glm_lognormal, newdata = data.frame(edad = edad_nueva)))

# Gráfico del ajuste
plot(edad, peso, main = "Peso vs Edad (Modelo Log-Normal)", xlab = "Edad (días)", ylab = "Peso (mg)", pch = 19)
lines(edad_nueva, predicciones_glm_lognormal, col = "#800080", lwd = 2)

Comparación de AIC para Modelo Log-Normal

aic_glm_lognormal <- AIC(modelo_glm_lognormal)
cat("AIC del Modelo Log-Normal:", aic_glm_lognormal, "\n")
## AIC del Modelo Log-Normal: 80.14779

Bibliografía