TAREA AJUSTE MODELO NO LINEAL

1. Modelo de crecimiento Schnute

En la Ecuación 1 se muestra la fórmula del modelo de crecimiento Schnute. Este modelo es comunmente utilizado para medir el crecimiento de peces.

\[ y(t) = [r_0 - \beta\ e^{(kt)}]^m \ \ \ \ (1) \]

1.1. Derivadas parciales

En la ecuación 2 se presenta la derivada del modelo de Schnute respecto a la variable tiempo (t). En las ecuaciones 3, 4, 5, y 6 se presentan las derivadas parciales para \(r_0\), \(\beta\), \(k\), y \(m\), respectivamente. Se observa que la ecuación es no lineal con respecto a todos los parámetros, ya que una vez derivada parcialmente, ésta sigue dependiendo de los mismos.

1.1.1. Respecto a la variable (t)

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

1.1.2. Respecto a parámetros

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

2. Uso de la librería ‘manipulate’ - Modelo asignado: Blumberg

Esta herramineta no funciona en markdown. Código bloqueado. En la ecuación 7 se presenta la fórmula del modelo de crecimiento Blumberg.

\[ y(t) = \frac{\alpha * (t + t_0)^m}{w_0 + (t + t_0)^m} \ \ \ \ (7) \]

#Blumbg <- function (alfa,w,x,x0,m) (alfa*(x+x0)^m)/(w+(x+x0)^m)
#curve(Blumbg(x, alpha = 1, w = 1, x0 = 1, m = 1), from = 0, to = 50, ylab = "y(x)")
#y = Blumbg(0:50, alpha = 1, w = 1, x0 = 1, m = 1)
#data = data.frame(x = 0:50, y)
#library(manipulate)
#manipulate(plot(data,xlim=c(x.min,x.max),ylim=c(0,200),
#type=tipo,axes=axes,ann=label,col=color),
#x.min=slider(0,10),
#x.max=slider(10,50),
#alpha=slider(1,15),
#w = slider(0,50),
#xo = 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","pink","black")
#)

3. Ejemplo comportamiento del modelo no lineal Schunte

library(growthmodels)
growth <- schnute(0:10, 10, 5, 0.5, 0.5)
plot(growth)

4. Generación de datos y comportamiento del modelo Schnute

Inicialmente se usó los datos dados en clase para el Modelo Glumberg. Se hizo la gráfica pero el modelo mostraba un comportamiento inicial muy pequeño y en los últimos tiempos el aumento en el crecimiento era muy elevado y no se evidenciaba el comportamiento de los datos. Por lo tanto, se deció generar unos datos con una escala menor. De acuerdo a lo anterior, se generaron 20 datos de tiempo. El modelo Schunete es usado para medir el crecimiento de peces; la variable W representa el crecimiento de peces y la variable t los tiempos de medicion del crecimiento. Al generar los datos de W se usaron los datos del ejemplo anterior, asumiendo los parámetros k y m con valor 1. La ecuación quedó de la siguiente manera:

\[ y(t) = r_0 - \beta\ e^{t} \ \ \ \ (8) \]

set.seed(1964)
t=rep(c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20),each=4)
mediast=tapply(as.numeric(t),as.factor(t),mean);mediast
##  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 
##  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20
W <-(runif(4,15,20)+(runif(4,1,6)*(exp(t))+rnorm(80,0.1,0.3)));length(W)
## [1] 80
mediasW=tapply(W,as.factor(t),mean);mediasW
##            1            2            3            4            5            6 
## 2.871132e+01 4.744611e+01 9.875071e+01 2.389975e+02 6.185814e+02 1.651993e+03 
##            7            8            9           10           11           12 
## 4.460010e+03 1.209342e+04 3.284328e+04 8.924690e+04 2.425678e+05 6.593375e+05 
##           13           14           15           16           17           18 
## 1.792235e+06 4.871770e+06 1.324281e+07 3.599767e+07 9.785178e+07 2.659887e+08 
##           19           20 
## 7.230322e+08 1.965405e+09
plot(t,W,pch=19,cex=0.6,ylab="Peso(g)",xlab="tiempo(días)")
points(mediast,mediasW,col="darkblue",pch=19)
points(18,mediasW[18],col="darkgreen",cex=1.2,pch=19)
points(19,mediasW[18],col="darkgreen",cex=1.2,pch=19)
points(19,mediasW[19],col="darkgreen",cex=1.2,pch=19)

xx = c(18, 19, 19)
yy = c(mediasW[18], mediasW[18], mediasW[19])
polygon(xx, yy, col = "lightblue")
text(18.7, (mediasW[13] + mediasW[19])/2, 'A', cex = 0.75)
text(17, mediasW[18], bquote(paste('(t'['18']*',','W'['18'],')')), cex = 0.75)
text(20, mediasW[18],bquote(paste('(t'['19']*',','W'['18'],')')), cex = 0.75)
text(18, mediasW[19], bquote(paste('(t'['19']*',','W'['19'],')')), cex = 0.75)
lines(mediast, mediasW, lty = 1, col = "darkred", lwd = 1); grid(15, 15)

5. Ajuste del modelo Schnute a datos observados

5.1. Cálculo de tasa de crecimiento media

Se calculó la tasa de crecimiento media para los intervalos de tiempo.

delta.mediasW = diff(mediasW)
delta.mediast = diff(mediast)
vector.TAC = round(delta.mediasW/delta.mediast, 4); vector.TAC 
##            2            3            4            5            6            7 
## 1.873480e+01 5.130460e+01 1.402468e+02 3.795839e+02 1.033412e+03 2.808017e+03 
##            8            9           10           11           12           13 
## 7.633408e+03 2.074987e+04 5.640361e+04 1.533209e+05 4.167697e+05 1.132898e+06 
##           14           15           16           17           18           19 
## 3.079535e+06 8.371043e+06 2.275485e+07 6.185411e+07 1.681369e+08 4.570435e+08 
##           20 
## 1.242373e+09

Se obtuvo una tasa de crecimiento de \(4.570435*10^8\). La pregunta es: ¿Es representativo este valor para ese intervalo de tiempo?

5.2. Ajuste del modelo Schnute

La función nls se usa para estimar los mínimos cuadrados no lineales de los parámetros del modelo. Aunque el estimador Akaike (AIC) es relativo y requiere de más de un modelo para su interpretación, se usó para estimar la calidad relativa de los modelos estadísticos para un conjunto de datos dado.

m <- nls(W~ r+beta*(exp(t)), start = list(r = 10, beta = 5))
summary(m)
## 
## Formula: W ~ r + beta * (exp(t))
## 
## Parameters:
##       Estimate Std. Error t value Pr(>|t|)    
## r    1.748e+01  2.352e+07    0.00        1    
## beta 4.051e+00  2.119e-01   19.12   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 195300000 on 78 degrees of freedom
## 
## Number of iterations to convergence: 1 
## Achieved convergence tolerance: 1.031e-08
cf = c(coef(m))
cor(W, predict(m))
## [1] 0.9179846
AIC(m)
## [1] 3285.404
r2 = cor(W, predict(m))**2; r2*100
## [1] 84.26958
plot(t, W, pch=19, xlim=c(0,20), ylim=c(0,3e+09), axes = TRUE, xlab="t",
     lab=c(length(W), 9,12))
curve(cf[1]+cf[2]*exp(x), type = 'l', xlim= c(0, 90), ylim = c(0, 30), col = "darkblue", add = T)
text(14, 1500000000, expression(W(t) == paste(17.75 + 0.405*e^{(t)})))
grid(15, 15)

x =  seq(18, 19, 0.1)
fit = cf[2]*exp(x);head(fit) #derivada del modelo
## [1] 265988663 293962935 324879287 359047140 396808457 438541167
mean(fit)
## [1] 460795758
median(fit)
## [1] 438541167
summary(fit)
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 265988663 341963213 438541167 460795758 563802024 723032150
sd(fit)
## [1] 151247824
cv = sd(fit)*100/mean(fit); cv
## [1] 32.82318

Anteriormente, se obtuvo un valor de la tasa de crecimiento medio para el intervalo de 18 a 19 días de \(4.570435*10^8\) g/día. De acuerdo con el modelo Schunte, el valor medio de crecimiento para el intervalo de los 18 a 19 días es de \(4.60795758*10^8\), en un rango de \(2.65988663*10^8\) a \(7.23032150*10^8\) g/día**.

5.3. Representación de la primera derivada del modelo Schnute

La derivada de la ecuación 8 con respecto a t se muestra en la ecuación 9:

\[ y(t) = r_0 + \beta\ e^{t} \ \ \ \ (8) \] \[ \frac{\partial\ y}{\partial\ t} = \beta\ e^{t}\ \ \ \ \ (2)\]

library(extrafont)
## Registering fonts with R
curve(cf[2]*exp(x), type = 'l', col = 'darkblue', xlim = c(0, 20),
      ylim=c(0,2e+09), ylab = "W'(t)", xlab = 'Tiempo(días)', main = 'Tasa de crecimiento absoluto (TAC)')
text(15, 1.0e+09, expression(dW/dt==paste(4.05*e^{x})),
     col = "black")
grid(15, 15)

TAREA VARIANZAS DESIGUALES

  1. Los datos siguientes tienen varianzas iguales.
set.seed(2011)
gradosbrix<-replicate(1000,rnorm(48,22,0.8))
l1=list()
l2=list()
var<-gl(3,16,48,paste0("V",1:3))
for(i in 1:dim(gradosbrix)[2]){
  l1[[i]]=(aov(gradosbrix[,i]~var))$residuals
  l2[[i]]=shapiro.test(l1[[i]])$p.value
}
table(unlist(l2)>0.05) #Resultados de la prueba shapiro a los residuales
## 
## FALSE  TRUE 
##    35   965
#lapply(l1,shapiro.test)
  1. Se generaron varianzas desiguales
gradosbrix<-c(runif(16,20,24),runif(16,15,35),runif(16,20,21)) #Varianzas desiguales 
var<-gl(3,16,48,paste0("V",1:3))
boxplot(gradosbrix~var,main="Varianzas desiguales",col=rainbow(3))

set.seed(2011)
gradosbrix1<-replicate(1000,rnorm(16,20,0.8))
gradosbrix2<-replicate(1000,rnorm(16,22,1.6))
gradosbrix3<-replicate(1000,rnorm(16,24,2.4))
gradosbrixT<-rbind(gradosbrix1,gradosbrix2,gradosbrix3)
var<-gl(3,16,48,paste0("V",1:3))
l1=list()
l2=list()
for(i in 1:dim(gradosbrixT)[2]){
  l1[[i]]=(aov(gradosbrixT[,i]~var))$residuals
  l2[[i]]=shapiro.test(l1[[i]])$p.value
}
table(unlist(l2)>0.05) #Resultados de la prueba shapiro a los residuales
## 
## FALSE  TRUE 
##   253   747
#lapply(l1,shapiro.test)

Conclusión: La cantidad de falsos aumenta cuando las varianzas son desiguales.