Se presenta el ajuste de un modelo no lineal de la libreria Growthmodels. Para el desarrollo de este ejercicio se tomó el modelo Weibull, propuesto por Ernst Hjalmar Waloddi Weibull en 1951. Este modelo tiene muchas aplicaciones en el crecimiento de poblaciones, el crecimiento agrícola, para describir la supervivencia en casos de lesiones o enfermedades o en estudios de dinámica de poblaciones.

En este ejercicio se generaron datos uniformes para cada parametro de acuedo a lo dispuesto en Mahanta y Borah (2014).

Ecuación 1 del modelo de crecimiento Weibull.

\[ y(t) = \alpha\ - \beta\ e^{(-k*t^m)} \ \ \ \ (1)\]

1. Derivadas parciales

En la ecuación 2 se presenta la derivada del modelo de Weibull respecto a la variable tiempo. Mientras que en la ecuaciones 3, 4, 5, y 6 se presentan las derivadas parciales para \(\alpha\), \(\beta\), \(k\), y \(m\), respectivamente. En las ecuaciones 5 y 6 observamos nuevamente los parámetros \(k\), y \(m\) en la primer derivada, por lo tanto, podemos decir que el modelo es no lineal para los parámetros.

Respecto a variables

\[ \frac{\partial\ y}{\partial\ t} = \beta\ e^{-kt^m}k\ m\ t^{m-1}\ \ \ \ \ (2)\]

Respecto a parámetros

\[ \frac{\partial\ y}{\partial\ \alpha} = 1 \ \ \ \ \ (3) \] \[ \frac{\partial\ y}{\partial\ \beta} = e^{-kt^m}\ \ \ \ \ (4)\] \[ \frac{\partial\ y}{\partial\ k} = \beta\ t^m e^{-kt^m}\ \ \ \ \ (5)\] \[ \frac{\partial\ y}{\partial\ m} = -\beta\ t^m e^{-kt^m}\ln(t)\ \ \ \ \ (6)\]

2. Función manipulate

Se muestra el código bloqueado para la aplicación de la librería manipulate, ya que esta no funciona en markdown.

#weibull <- function(alpha, beta, k, x, m) (alpha - (beta*exp(-k*x^m)))
#curve(weibull(x, alpha = 1, beta = 1, k = 1, m = 1), from = 0, to = 50, ylab = "y(x)")
#y = weibull(0:50, alpha = 1, beta = 1, k = 1, m = 1)
#data = data.frame(x = 0:50, y)
#library(manipulate)
#manipulate(plot(data,xlim=c(x.min,x.max),ylim=c(0,alpha),
#                type=tipo,axes=axes,ann=label,col=color),
#           x.min=slider(0,10),
#           x.max=slider(10,50),
#           alpha=slider(1,15),
#           beta=slider(0,50),
#           k = slider(0,50),
#           m = slider(0,50),
#           tipo= picker("p","l","b","c","o","h","s","S","n"),
#           axes=checkbox(TRUE, "dibujar ejes"),
#           label=checkbox(FALSE, "dibujar ejes "),
#           color=picker("red","blue","green")
#)

3. Generación de datos

Ahora se presenta la generación de datos de crecimiento de plantas en 13 tiempos. La variable W representa el crecimiento de plantas y la variable t los tiempos de medicion del crecimiento. En un intento 1 para la generación de datos de crecimiento (W) se fijo valores constantes para \(\alpha\) y \(\beta\). Para los otros parámetros se generaron datos con la función runif, los valores mínimos y máximos fueron inventados, pero la tendencia de los datos no aparentaba seguir el modelo Weibull, Por lo tanto, se hizo un intento 2 generando datos de crecimiento para todas los parámetros del modelo, los valores mínimos y máximos fueron sacados de Mahanta y Borah (2014), En este último se observo una tendencia más similar a la descrita por el modelo de crecimiento de Weibull.

Intento 1 - Para parametros t=x, k y m

set.seed(1964)
t=rep(c(0,7,14,21,28,35,42,49,56,63,70,77,84),each=4)
mediast=tapply(as.numeric(t),as.factor(t),mean);mediast
##  0  7 14 21 28 35 42 49 56 63 70 77 84 
##  0  7 14 21 28 35 42 49 56 63 70 77 84
W <- sort(1- 3*exp(-1*t^(runif(4,1,3)))+rnorm(52,1,4));length(W)
## [1] 52
mediasW=tapply(W,as.factor(t),mean);mediasW
##          0          7         14         21         28         35         42 
## -6.4166475 -3.5588335 -2.0852063 -0.7417551  0.6481116  1.8679434  2.5648219 
##         49         56         63         70         77         84 
##  3.6276395  4.5564358  4.8734982  5.7041642  6.8396964  9.4465955
plot(t,W,pch=19,cex=0.6,ylab="Altura(cm)",xlab="tiempo(días)")
points(mediast,mediasW,col="darkblue",pch=19)
points(7,mediasW[2],col="darkgreen",cex=1.2,pch=19)
points(14,mediasW[2],col="darkgreen",cex=1.2,pch=19)
points(14,mediasW[3],col="darkgreen",cex=1.2,pch=19)
xx = c(7,14,14)
yy = c(mediasW[2], mediasW[2], mediasW[3])
polygon(xx, yy, col = "lightblue")
text(12, (mediasW[2]+mediasW[3])/2, "A", cex = 0.9)
text(3, mediasW[2], bquote(paste('(t'['2']*',','W'['2'],')')), cex= 0.75)
text(19, mediasW[2], bquote(paste('(t'['3']*',','W'['2'],')')), cex= 0.75)
text(10, mediasW[3], bquote(paste('(t'['3']*',','W'['2'],')')), cex= 0.75)
lines(mediast, mediasW, lty = 1, col = 'darkred', lwd = 1); grid(10,10)

Intento 2 - Todos los parametros

En el intento 2 agregamos todos los parametros del modelo Weibull. Los valores generados con distribución uniforme provenienen del artículo de Mahanta y Borah (2014).

set.seed(1964)
t=rep(c(0,7,14,21,28,35,42,49,56,63,70,77,84),each=4)
mediast=tapply(as.numeric(t),as.factor(t),mean);mediast
##  0  7 14 21 28 35 42 49 56 63 70 77 84 
##  0  7 14 21 28 35 42 49 56 63 70 77 84
W=runif(4,21,28)-runif(4,14,22)*exp(-(runif(4,0.052,0.087))*t^runif(4,1.05,1.56))+rnorm(52,1,1);length(W)
## [1] 52
mediasW=tapply(W,as.factor(t),mean);mediasW
##         0         7        14        21        28        35        42        49 
##  6.403108 16.564918 22.978820 23.406832 25.369291 25.184458 25.330509 25.796734 
##        56        63        70        77        84 
## 25.705879 25.424430 25.312552 25.852404 25.916352
plot(t,W,pch=19,cex=0.6,ylab="Altura(cm)",xlab="tiempo(días)")
points(mediast,mediasW,col="darkblue",pch=19)
points(7,mediasW[2],col="darkgreen",cex=1.2,pch=19)
points(14,mediasW[2],col="darkgreen",cex=1.2,pch=19)
points(14,mediasW[3],col="darkgreen",cex=1.2,pch=19)
xx = c(7,14,14)
yy = c(mediasW[2], mediasW[2], mediasW[3])
polygon(xx, yy, col = "lightblue")
text(12, (mediasW[2]+mediasW[3])/2, "A", cex = 0.9)
text(3, mediasW[2], bquote(paste('(t'['2']*',','W'['2'],')')), cex= 0.75)
text(19, mediasW[2], bquote(paste('(t'['3']*',','W'['2'],')')), cex= 0.75)
text(10, mediasW[3], bquote(paste('(t'['3']*',','W'['2'],')')), cex= 0.75)
lines(mediast, mediasW, lty = 1, col = 'darkred', lwd = 1); grid(10,10)

4. Ajuste de un modelo no líneal a los datos observados.

Calculo de tasa de crecimiento media.

Antes de ajustar el modelo no líneal de Weibull a los datos generados en el intento 2. Se calculó el la tasa de crecimiento media, como se acostumbra a calcular en las ciencias agrícolas. Para el intervalo de 7 a 14 días se obtuvó una tasa de crecimiento de 0.916 cm/día. Sin embargo, la cuestión que surge es si este dato puntual representa la variación en la tasa de crecimiento de ese intervalo.

delta.mediasW = diff(mediasW)
delta.mediast = diff(mediast)
vector.TAC = round(delta.mediasW/delta.mediast, 4); vector.TAC 
##       7      14      21      28      35      42      49      56      63      70 
##  1.4517  0.9163  0.0611  0.2804 -0.0264  0.0209  0.0666 -0.0130 -0.0402 -0.0160 
##      77      84 
##  0.0771  0.0091

Ajuste modelo no lineal Weibull

Ahora se utilizó la función nls para la estimación de mínimos cuadrados no lineales de los parámetros del modelo no lineal. El criterio de Akaike mediante la función AIC es un estimador del error de predicción, usado para estimar la calidad relativa de los modelos estadísticos para un conjunto de datos dado. Sin embargo, este estimador es relativo y su interpretación depende de la comparación del ajuste de más de un modelo.

m <- nls(W~ a-b*exp(-c*t^m), start = list(a = 20, b = 10, c = 0.05, m = 1))
summary(m)
## 
## Formula: W ~ a - b * exp(-c * t^m)
## 
## Parameters:
##   Estimate Std. Error t value Pr(>|t|)    
## a 25.57519    0.59865  42.721  < 2e-16 ***
## b 19.20250    1.74590  10.999 1.02e-14 ***
## c  0.08568    0.08247   1.039  0.30405    
## m  1.13867    0.40463   2.814  0.00707 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.266 on 48 degrees of freedom
## 
## Number of iterations to convergence: 10 
## Achieved convergence tolerance: 9.285e-06
cf = c(coef(m))
cor(W, predict(m))
## [1] 0.8630646
AIC(m)
## [1] 276.5016
r2 = cor(W, predict(m))**2; r2*100
## [1] 74.48806
#Gráfico 
plot(t, W, pch=19, xlim=c(0,84), ylim=c(0,25), axes = TRUE, xlab="t",
     lab=c(length(W), 9,12))
curve((cf[1]-cf[2]*exp(-cf[3]*x^cf[4])), type = 'l', xlim= c(0, 90),
      ylim = c(0, 30), col = "darkblue", add = T)
text(42, 14, expression(W(t) == paste(25.75 - 12.02*e^{(-0.085*t^1.14)})))

Anteriormente, obtuvimos un valor de tasa de crecimiento medio para el intervalo de 7 a 14 días de 0.916 cm/día. De acuerdo con el modelo Weibull el valor medio para el intervalo de tiempo de 7 a 14 días es de 0.765 cm/día, pero el rango es de 0.4791 cm/día a 1.118 cm/día. Aquí evidenciamos que el rango de variación de la tasa de crecimiento es de 0.4791 cm/día a 1.118 cm/día, un rango bastante amplio que no es representado por la media de 0.916 cm/día.

x =  seq(7, 14, 0.1)
fit = c(cf[2]*exp(-cf[3]*x^cf[4])*cf[3]*cf[4]*x**(cf[4]-1));head(fit)
## [1] 1.118577 1.106536 1.094566 1.082671 1.070851 1.059107
mean(fit)
## [1] 0.7645294
median(fit)
## [1] 0.746307
summary(fit)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.4791  0.6003  0.7463  0.7645  0.9192  1.1186
sd(fit)
## [1] 0.1890405
cv = sd(fit)*100/mean(fit); cv
## [1] 24.72639

Representación gráfica primer derivada.

curve((cf[2]*cf[3]*cf[4]*x^(cf[4]-1)*exp(-cf[3]*x^cf[4])), type = 'l',
      col = 'darkblue', xlim = c(0, 84), ylab = "W'(t)", xlab = 'tiempo(días)',
      main = 'Tasa de crecimiento absoluto (TAC)')
text(50, 0.5, expression(dW/dt==paste(1.84*t^{0.14}*e^{-0.085*t^1.14})),
     col = "darkgreen")
grid(20, 20)

Bibliografía

Mahanta, D. J., Borah. (2014). Parameter Estimation of Weibull Growth Models in Forestry. International Journal of Mathematics Trends and Technology, 8(3), 157-163.