datos <- read.csv2("AbarcoData.csv")
datos$D <- datos$D/10Aplicación de simulación Monte Carlo: crecimiento de abarco
Objetivo
Evaluar el efecto del tamaño de muestra en la estimación de los parámetros de la ecuación no lineal tipo exponencial negativa y en su incertidumbre, para la modelación del crecimiento diamétrico de árboles de abarco (Cariniana pyriformis).
Datos
Los datos consisten de 199 registros de edad y diámetro para árboles de abarco, obtenidas con técnicas dendrocronológicas.
par(lwd=2, family="sans", mfrow=c(1,1))
with(datos, plot(year, D, bty="L", xlab="Edad (años)", ylab="Diámetro (cm)",
cex.lab=1.2, las=1, cex.axis=1, col="#458D1B", pch=21,
cex=1.5, bg="white", lwd=2))
Se pueden presentar árboles de la edad entre 1 y 181 años, y diámetro entre 0.319 y 191.58 cm.
| Variable | Unidades | Mín | Máx | Media |
|---|---|---|---|---|
| Diámetro | cm | 0.319 | 191.58 | 49.23 |
| Edad | años | 1 | 181 | 65.98 |
Modelo
Debido a que los datos de diámetro \(D\) en función de edad \(t\) presentan heterogeneidad de varianzas, se considera la relación de \(log(D)\)~\(t\). El modelo no lineal a ajustar es el exponencial negativo (Panik, 2014), cuya expresión está dada por:
\[log(D) = \alpha(1-e^{-t/k} ) - \omega + \epsilon_{i},\]
donde \(\epsilon_{i}\)~\(N(0,\sigma^2)\).
Los parámetros estimados por mínimos cuadrados no lineales (nls en R) resultan en \(\alpha=4.74\), \(k=31\), y \(\omega=0.2\). Estos valores representan la asíntota, tasa intrínseca de crecimiento e intercepto, respectivamente.
modelo <- nls(log(D) ~ alpha*(1-exp(-year/k)) - w,
data = datos,
start = list(alpha=5, k=42, w=2))
summary(modelo)
Formula: log(D) ~ alpha * (1 - exp(-year/k)) - w
Parameters:
Estimate Std. Error t value Pr(>|t|)
alpha 4.7412 0.1132 41.897 <2e-16 ***
k 31.0702 1.7564 17.689 <2e-16 ***
w 0.2003 0.1195 1.676 0.0953 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.3983 on 196 degrees of freedom
Number of iterations to convergence: 6
Achieved convergence tolerance: 4.261e-06De la tabla resumen anterior podemos ver que los parámetros estimado son \(\hat{\alpha}=4.7412\), \(\hat{k}=31.0702\), \(\hat{w}=0.2003\) y \(\hat{\sigma}=0.3983\).
Los tres primeros parámetros del modelo los vamos a almacenar en el objeto betas para luego crear la curva correspondiente al modelo ajustado.
betas <- coef(modelo)
betas alpha k w
4.7411694 31.0702215 0.2003205 Para agregar el modelo ajustado al diagrama de dispersión original hacemos lo siguiente:
par(lwd=2, family="sans", mfrow=c(1,1), mar=c(5, 4, 1.5, 2))
with(datos, plot(year, log(D), bty="L", xlab="Edad (años)", ylab="log(D)",
cex.lab=1.2, las=1, cex.axis=1, col="#458D1B", pch=21,
cex=1.5, bg="white", lwd=2))
x <- seq(0,180,1)
y <- betas[1]*(1-exp(-x/betas[2])) - betas[3]
lines(x, y, col=1, lwd=2)
De la figura anterior vemos que el modelo ajustado logra explicar los datos.
Simulación
Con base en el vector de parámetros estimado \(\hat\Theta\) vamos a generar nuevas muestras aleatorias de tamaño \(n=10, 20, ..., 1000\) tomadas del modelo
\[log(D) = 4.7(1-e^{-t/31} ) - 0.2 + \epsilon_{i},\]
con \(\epsilon_{i}\)~\(N(0,0.4^2)\).
El proceso de simulación consiste de 1000 repeticiones por cada tamaño de muestra \(n_i\). Para cada repetición de la muestra simulada, se ajusta el modelo exponencial negativo y se extraen los parámetros \(\hat\alpha\) y \(\hat{k}\). Se calcula el promedio de los parámetros estimados para cada \(n_i\), y se evalúa su incertidumbre.
La simulación se realiza con la librería parSim de R:
library(parSim)
parSim(
### SIMULATION CONDITIONS
n = seq(from=10, to=1000, by=10),
reps = 1000, # Repetitions
write = TRUE, # Writing to a file
name = "results_parSim", # Name of the file
nCores = 1, # Number of cores to use
expression = {
# True parameter values
alpha = 4.7
k = 31
w = 0.2
sigma = 0.4
# Generate data:
edad <- runif(n=n, min=1, max=181)
ei <- rnorm(n=n, mean=0, sd=sigma)
log_D <- alpha*(1-exp(-edad/k)) - w + ei
dat <- data.frame(log_D=log_D, edad=edad)
# Run analysis:
mod <- nls(log_D ~ alpha*(1-exp(-edad/k)) - w,
data = dat,
start = list(alpha=4, k=40, w=2))
# Store coefficients:
alpha_hat <- coef(mod)[1]
k_hat <- coef(mod)[2]
w_hat <- coef(mod)[3]
sigma_hat <- summary(mod)$sigma
# Results list:
Results <- list(
alpha_hat = alpha_hat,
k_hat = k_hat,
w_hat = w_hat,
sigma_hat = sigma_hat
)
# Return:
Results
}
)Resultados
La base de datos resultante de la simulación consiste de 100000 registros, correspondientes a 1000 repeticiones de 100 tamaños de muestra distintos. El método nls falló 169 veces, de las cuales 132 correspondieron a \(n=10\), y 18 a \(n=20\). Posiblemente, los valores de inicio no fueron adecuados para algunas muestras simuladas, especialmente cuando el tamaño era pequeño, por lo tanto, el algoritmo no alcanzaba convergencia.
datos <- read.table("results_parSim.txt", header = TRUE)
require(dplyr)
dat <- datos %>% group_by(n) %>%
summarise(alpha = mean(alpha_hat, na.rm=TRUE),
alpha_rmse = sqrt(sum((alpha_hat[!is.na(alpha_hat)] - mean(alpha_hat, na.rm=TRUE))^2)/length(alpha_hat[!is.na(alpha_hat)])),
k = mean(k_hat, na.rm=TRUE),
k_rmse = sqrt(sum((k_hat[!is.na(k_hat)] - mean(k_hat, na.rm=TRUE))^2)/length(k_hat[!is.na(k_hat)])),
w = mean(w_hat, na.rm=TRUE),
sigma = mean(sigma_hat, na.rm=TRUE))library(kableExtra)
kable(dat[, -c(6,7)], format="html", digits=3, align="c",
col.names=c("n", "alpha", "RMSE (alpha)", "k", "RMSE (k)")) %>% kable_styling() %>%
scroll_box(height = "450px")| n | alpha | RMSE (alpha) | k | RMSE (k) |
|---|---|---|---|---|
| 10 | 5.160 | 4.129 | 37.677 | 31.563 |
| 20 | 4.883 | 1.186 | 32.390 | 11.053 |
| 30 | 4.743 | 0.541 | 31.960 | 6.537 |
| 40 | 4.725 | 0.369 | 31.348 | 4.620 |
| 50 | 4.717 | 0.383 | 31.361 | 4.349 |
| 60 | 4.730 | 0.314 | 31.228 | 3.774 |
| 70 | 4.714 | 0.272 | 31.438 | 3.364 |
| 80 | 4.721 | 0.258 | 31.165 | 3.075 |
| 90 | 4.711 | 0.238 | 31.177 | 2.877 |
| 100 | 4.712 | 0.224 | 31.135 | 2.727 |
| 110 | 4.718 | 0.218 | 31.023 | 2.632 |
| 120 | 4.706 | 0.189 | 31.153 | 2.443 |
| 130 | 4.713 | 0.191 | 31.083 | 2.331 |
| 140 | 4.703 | 0.182 | 31.173 | 2.291 |
| 150 | 4.705 | 0.178 | 31.062 | 2.183 |
| 160 | 4.694 | 0.161 | 31.157 | 2.084 |
| 170 | 4.697 | 0.164 | 31.107 | 2.023 |
| 180 | 4.706 | 0.155 | 31.144 | 2.101 |
| 190 | 4.701 | 0.156 | 31.137 | 1.982 |
| 200 | 4.702 | 0.144 | 31.185 | 1.892 |
| 210 | 4.706 | 0.145 | 31.106 | 1.809 |
| 220 | 4.696 | 0.141 | 31.055 | 1.767 |
| 230 | 4.700 | 0.142 | 31.052 | 1.763 |
| 240 | 4.698 | 0.134 | 31.001 | 1.767 |
| 250 | 4.712 | 0.128 | 30.995 | 1.627 |
| 260 | 4.698 | 0.136 | 31.135 | 1.634 |
| 270 | 4.702 | 0.129 | 31.163 | 1.556 |
| 280 | 4.703 | 0.127 | 31.066 | 1.608 |
| 290 | 4.701 | 0.123 | 31.016 | 1.585 |
| 300 | 4.699 | 0.124 | 31.114 | 1.617 |
| 310 | 4.704 | 0.121 | 31.045 | 1.483 |
| 320 | 4.701 | 0.116 | 31.055 | 1.445 |
| 330 | 4.701 | 0.117 | 31.034 | 1.415 |
| 340 | 4.708 | 0.109 | 31.083 | 1.410 |
| 350 | 4.710 | 0.112 | 31.050 | 1.401 |
| 360 | 4.700 | 0.113 | 31.023 | 1.365 |
| 370 | 4.699 | 0.105 | 31.027 | 1.320 |
| 380 | 4.700 | 0.104 | 31.057 | 1.324 |
| 390 | 4.697 | 0.109 | 31.058 | 1.402 |
| 400 | 4.704 | 0.104 | 31.023 | 1.302 |
| 410 | 4.697 | 0.102 | 31.039 | 1.308 |
| 420 | 4.702 | 0.103 | 31.028 | 1.307 |
| 430 | 4.699 | 0.100 | 31.089 | 1.252 |
| 440 | 4.703 | 0.100 | 31.006 | 1.217 |
| 450 | 4.703 | 0.100 | 31.076 | 1.235 |
| 460 | 4.703 | 0.098 | 31.024 | 1.226 |
| 470 | 4.702 | 0.099 | 30.972 | 1.256 |
| 480 | 4.704 | 0.094 | 31.006 | 1.200 |
| 490 | 4.698 | 0.092 | 31.008 | 1.147 |
| 500 | 4.704 | 0.091 | 31.049 | 1.198 |
| 510 | 4.706 | 0.093 | 30.974 | 1.183 |
| 520 | 4.703 | 0.092 | 30.961 | 1.162 |
| 530 | 4.706 | 0.092 | 31.016 | 1.192 |
| 540 | 4.701 | 0.090 | 31.015 | 1.123 |
| 550 | 4.702 | 0.090 | 31.061 | 1.180 |
| 560 | 4.703 | 0.085 | 31.001 | 1.104 |
| 570 | 4.703 | 0.084 | 31.024 | 1.057 |
| 580 | 4.704 | 0.083 | 31.025 | 1.109 |
| 590 | 4.701 | 0.084 | 31.080 | 1.110 |
| 600 | 4.700 | 0.082 | 31.016 | 1.026 |
| 610 | 4.702 | 0.083 | 30.998 | 1.076 |
| 620 | 4.699 | 0.086 | 31.096 | 1.059 |
| 630 | 4.701 | 0.080 | 30.998 | 1.043 |
| 640 | 4.706 | 0.082 | 30.990 | 0.973 |
| 650 | 4.703 | 0.081 | 31.032 | 1.004 |
| 660 | 4.700 | 0.081 | 30.971 | 0.997 |
| 670 | 4.697 | 0.079 | 31.120 | 1.044 |
| 680 | 4.703 | 0.080 | 30.982 | 0.991 |
| 690 | 4.699 | 0.077 | 31.021 | 0.984 |
| 700 | 4.696 | 0.079 | 31.041 | 0.992 |
| 710 | 4.700 | 0.078 | 31.021 | 0.979 |
| 720 | 4.702 | 0.077 | 31.017 | 1.002 |
| 730 | 4.701 | 0.075 | 31.040 | 0.935 |
| 740 | 4.698 | 0.082 | 31.087 | 0.992 |
| 750 | 4.704 | 0.078 | 31.008 | 1.004 |
| 760 | 4.705 | 0.076 | 31.012 | 0.951 |
| 770 | 4.699 | 0.075 | 31.025 | 0.940 |
| 780 | 4.705 | 0.075 | 30.986 | 0.937 |
| 790 | 4.705 | 0.074 | 31.002 | 0.912 |
| 800 | 4.700 | 0.072 | 31.034 | 0.929 |
| 810 | 4.701 | 0.070 | 31.046 | 0.909 |
| 820 | 4.703 | 0.072 | 30.991 | 0.942 |
| 830 | 4.702 | 0.074 | 31.033 | 0.922 |
| 840 | 4.700 | 0.071 | 31.011 | 0.869 |
| 850 | 4.700 | 0.074 | 31.059 | 0.906 |
| 860 | 4.701 | 0.071 | 30.999 | 0.877 |
| 870 | 4.704 | 0.073 | 30.980 | 0.886 |
| 880 | 4.701 | 0.069 | 31.011 | 0.877 |
| 890 | 4.702 | 0.071 | 31.007 | 0.907 |
| 900 | 4.697 | 0.067 | 31.015 | 0.857 |
| 910 | 4.700 | 0.070 | 31.066 | 0.855 |
| 920 | 4.704 | 0.068 | 30.946 | 0.857 |
| 930 | 4.700 | 0.071 | 31.042 | 0.861 |
| 940 | 4.699 | 0.068 | 31.005 | 0.864 |
| 950 | 4.701 | 0.066 | 31.024 | 0.839 |
| 960 | 4.700 | 0.066 | 30.981 | 0.837 |
| 970 | 4.702 | 0.070 | 30.999 | 0.866 |
| 980 | 4.700 | 0.067 | 30.988 | 0.843 |
| 990 | 4.699 | 0.067 | 31.023 | 0.860 |
| 1000 | 4.702 | 0.065 | 30.996 | 0.804 |
La media de los parámetros estimados se aleja de los valores verdaderos cuando se tienen muestras pequeñas, sin embargo, a medida que \(n\) crece, \(\hat\alpha\) y \(\hat{k}\) se parecen más a \(\alpha\) y \(k\), respectivamente. Las estimaciones son inestables para \(n\) pequeños, ya que la incertidumbre (RMSE) es mayor. Por ejemplo, para \(\alpha\) con \(n=10\), el RMSE es de 6.93, cuyo valor disminuye a 0.320 con \(n=60\). A \(n\) mayores, el RMSE disminuye gradualmente.
library(latex2exp)
# True parameter values
alpha <- 4.7
k <- 31
w <- 0.2
sigma <- 0.4
# Gráfico con incertidumbre para alpha
par(lwd=2, family="sans", mfrow=c(1,1), mar=c(4, 4, 0, 2))
with(dat, plot(x=n, y=alpha, type='l', las=1,
ylim=c(4,6), bty="L",
xlab="n", ylab=expression(alpha)))
# Valor de referencia
abline(h=alpha, col='dodgerblue2', lty='dashed')
# Intervalo
lines(dat$n, dat$alpha+dat$alpha_rmse, col=2)
lines(dat$n, dat$alpha-dat$alpha_rmse, col=2)
legend(x="topright",
legend=c(expression(hat(alpha)),
TeX(r'($\alpha$ verdadero)'),
TeX(r'(Intervalo con MSE para $\hat{\alpha}$)')),
col=c(1,'dodgerblue2', 2), lty=c(1,2,1), bty="n")
# Gráfico con incertidumbre para k
par(lwd=2, family="sans", mfrow=c(1,1), mar=c(4, 4, 0, 2))
with(dat, plot(x=n, y=k, type='l', las=1,
ylim=c(25,40), bty="L", xlab="n"))
# Valor de referencia
abline(h=k, col='dodgerblue2', lty='dashed')
# Intervalo
lines(dat$n, dat$k+dat$k_rmse, col=2)
lines(dat$n, dat$k-dat$k_rmse, col=2)
legend(x="topright",
legend=c(expression(hat(k)),
TeX(r'($k$ verdadero)'),
TeX(r'(Intervalo con MSE para $\hat{k}$)')),
col=c(1,'dodgerblue2', 2), lty=c(1,2,1), bty="n")
Conclusiones
El tamaño de la muestra tiene influencia en la estimación de los parámetros de la ecuación tipo exponencial negativa para crecimiento diamétrico. La incertidumbre disminuye con tamaño de muestra mayor.
La evaluación de la incertidumbre en la modelación del crecimiento diamétrico tiene implicaciones prácticas en el manejo forestal de la especie abarco. Por ejemplo, \(\alpha\) es un parámetro asociado a la talla máxima que pueden alcanzar los árboles, pero si se modela con pocas observaciones la incertidumbre es alta, lo que implica que se pueden obtener amplios valores de \(\alpha\), donde algunos llegan a ser probablemente erróneos. Asimismo, a partir de la talla máxima se estima la edad promedio de la especie, pero si su incertidumbre es alta, la edad estimada tendría poca precisión.