1. Cargar los datos
set.seed(2023)
t<-rep(seq(0,77,7),each=3)
W<- c(0,0,0,5.69,5.74,5.68,6.52,6.58,6.51,6.96,7.03,6.95,7.26,7.32,7.25,7.47,7.54,7.46,7.64,7.71,7.63,7.78,7.84,7.77,7.89,7.96,7.88,7.99,8.05,7.98,8.07,8.14,8.06,8.14,8.21,8.14)
W<-W+rnorm(36,0.01,0.03)
plot(W~t,pch=16)
abline(v=t,col="grey",lty=2)

medias = tapply(W,t,mean)
medias
##           0           7          14          21          28          35 
## -0.01941795  5.71604500  6.54355313  6.98446092  7.29289164  7.51746316 
##          42          49          56          63          70          77 
##  7.68057307  7.81835025  7.89078219  8.05055745  8.10190797  8.16210642
plot(W~t,pch=16)
points(medias~unique(t),pch=16,col="red",cex=1)
abline(v=t,col="grey",lty=2)

2. Medias en cada punto

medias_diff = diff(medias)
  1. Tasa de crecmiento
TAC = medias_diff/7
  1. Se ajusta el modelo Brody growth model Minimos cuadrados no lineales
model_2=nls(W~a*(1-exp(-k*t)),start = list(a=10,k=0.3))
summary(model_2)
## 
## Formula: W ~ a * (1 - exp(-k * t))
## 
## Parameters:
##   Estimate Std. Error t value Pr(>|t|)    
## a 7.766891   0.067412  115.22   <2e-16 ***
## k 0.159746   0.009696   16.48   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3317 on 34 degrees of freedom
## 
## Number of iterations to convergence: 7 
## Achieved convergence tolerance: 8.676e-06

\[ W = \alpha * ( 1 - \exp(-kt))\]

plot(W~t,pch=16)
curve(coef(model_2)[1] * (1 - exp(-coef(model_2)[2] * x)),add=T,col="red")

predicciones_2 = fitted.values(model_2)
Metrics::rmse(actual = W,predicted = predicciones_2)
## [1] 0.322335
  1. La derivada del modelo
f2= expression(a *(1 - exp(-k*t)))
D(f2,"t")
## a * (exp(-k * t) * k)
curve((coef(model_2)[1]*coef(model_2)[2])* exp(-coef(model_2)[2]*x),from = 0,to = 77,ylab = "D")
abline(v=c(14,21))

TAC_i2 = function(x){
  y=(coef(model_2)[1]*coef(model_2)[2])* exp(-coef(model_2)[2]*x)
  return(y)
}

x=seq(14,21,length.out=1000)

y= TAC_i2(x)

hist(y);abline(v=mean(y),col="red",lwd=2);abline(v=quantile(y,0.25),col="blue",lwd=2)

# Finding distance matrix
y_mat2 <- dist(y, method = 'euclidean')

Hierar_cl2 <- hclust(y_mat2, method = "ward.D2")

 
# Cutting tree by no. of clusters
fit <- cutree(Hierar_cl2, k = 6)
df_c=data.frame(y,cluster=fit)
boxplot(df_c$y~df_c$cluster)
points(x=1:6,y=tapply(df_c$y,df_c$cluster,mean),pch=16,col="red")
abline(h=mean(df_c$y),lwd=2,col="blue",lty=2)

plot(x=x,y=y,col=as.factor(df_c$cluster))

pesos=table(df_c$cluster)/1000
medias=tapply(df_c$y,df_c$cluster,mean)
mp = t(pesos)%*%medias
boxplot(df_c$y~df_c$cluster)
points(x=1:6,y=tapply(df_c$y,df_c$cluster,mean),pch=16,col="red")
abline(h=mean(df_c$y),lwd=2,col="blue",lty=2)
abline(h=mp,lwd=2,col="orange",lty=2)

## Comparando entre los días 7 a 14

curve((coef(model_2)[1]*coef(model_2)[2])* exp(-coef(model_2)[2]*x),from = 0,to = 77,ylab = "D")
abline(v=c(8,14))

#############################################################33

Comparación entre los 42 a 56

curve((coef(model_2)[1]*coef(model_2)[2])* exp(-coef(model_2)[2]*x),from = 0,to = 77,ylab = "D")
abline(v=c(42,56))

TAC_i2 = function(x){
  y=(coef(model_2)[1]*coef(model_2)[2])* exp(-coef(model_2)[2]*x)
  return(y)
}

x=seq(42,56,length.out=1000)

y= TAC_i2(x)

hist(y);abline(v=mean(y),col="red",lwd=2);abline(v=quantile(y,0.25),col="blue",lwd=2)

# Finding distance matrix
y_mat2 <- dist(y, method = 'euclidean')

Hierar_cl2 <- hclust(y_mat2, method = "ward.D2")

 
# Cutting tree by no. of clusters
fit <- cutree(Hierar_cl2, k = 6)
df_c=data.frame(y,cluster=fit)
boxplot(df_c$y~df_c$cluster)
points(x=1:6,y=tapply(df_c$y,df_c$cluster,mean),pch=16,col="red")
abline(h=mean(df_c$y),lwd=2,col="blue",lty=2)

Se evidencia que el cambio es más pequeño en sitios donde la pendite es casi cero

plot(x=x,y=y,col=as.factor(df_c$cluster))

pesos=table(df_c$cluster)/1000
medias=tapply(df_c$y,df_c$cluster,mean)
mp = t(pesos)%*%medias
boxplot(df_c$y~df_c$cluster)
points(x=1:6,y=tapply(df_c$y,df_c$cluster,mean),pch=16,col="red")
abline(h=mean(df_c$y),lwd=2,col="blue",lty=2)
abline(h=mp,lwd=2,col="orange",lty=2)

b<-0:8

Datos de papa

  1. Cargar los datos, la curva de crecimiento es sigmoide.
library(readxl)
datos_entrega_1 <- read_excel("datos_entrega_1.xlsx", 
    sheet = "base")
population <- datos_entrega_1
plot(population$tiempo, population$peso)

  1. Encontrando los parametros para el modelo 1
pop.ss <- nls(peso ~ SSlogis(tiempo, phi1, phi2, phi3), data = population)
summary(pop.ss)
## 
## Formula: peso ~ SSlogis(tiempo, phi1, phi2, phi3)
## 
## Parameters:
##      Estimate Std. Error t value Pr(>|t|)    
## phi1  459.875     30.098  15.279 0.000107 ***
## phi2   79.288      2.959  26.796 1.15e-05 ***
## phi3   13.886      2.211   6.281 0.003282 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 21.06 on 4 degrees of freedom
## 
## Number of iterations to convergence: 1 
## Achieved convergence tolerance: 2.342e-06

2.2 Ensayando un segundo modelo

t <- c(35, 51, 66, 79, 92, 106, 12)
p <- c(3.275, 58.815, 146.475, 203.690 , 331.350, 421.200, 422.910)

# Ajustar un modelo de regresión linealizado del crecimiento logístico
modelo_4 <- nls(p ~ a * t / (b + t), start = list(a = 1, b = 1))

# Resumen del modelo
summary(modelo_4)
## 
## Formula: p ~ a * t/(b + t)
## 
## Parameters:
##   Estimate Std. Error t value Pr(>|t|)    
## a   10.907     27.250   0.400 0.705491    
## b  -75.548      9.312  -8.113 0.000462 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 302.9 on 5 degrees of freedom
## 
## Number of iterations to convergence: 20 
## Achieved convergence tolerance: 6.064e-06
  1. Graficando la curva que se ajusta al modelo
alpha <- coef(pop.ss)  
b2<- coef(modelo_4)
plot(peso ~ tiempo, data = population, main = "Crecimiento tuberculo clon 1", 
    xlab = "DDE", ylab = "Peso")  # Census data
  curve(alpha[1]/(1 + exp(-(x - alpha[2])/alpha[3])), add = T, col = "blue")

3. Predicciones del modelo

predict.pop.ss<- predict(pop.ss, data.frame(tiempo = c(79,92)))
predict.pop.ss
## [1] 227.5515 328.4038
## attr(,"gradient")
##           phi1     phi2       phi3
## [1,] 0.4948118 -8.27872  0.1718114
## [2,] 0.7141156 -6.76128 -6.1896642
predicciones_4 = fitted.values(modelo_4)
predicciones_5 = fitted.values(pop.ss)
  1. Comparación de los modelos
Metrics::rmse(actual = p,predicted = predicciones_4)
## [1] 255.9839
Metrics::rmse(actual = p,predicted = predicciones_5)
## [1] 15.92343

El modelo de mejor ajustes en el crecimiento del tubérculo es

\[ y = \frac{\phi1}{1 + \exp\left(-\frac{x - \phi2}{\phi3}\right)} \]

5.Hallando la derivada

f= expression(a/(1 + exp(-(x - b)/c)))
D(f,"x")
## a * (exp(-(x - b)/c) * (1/c))/(1 + exp(-(x - b)/c))^2
curve((alpha[1]*exp((alpha[2]-x)/alpha[3]))/
        (alpha[3]* exp((alpha[2]-x)+1)^2), from = 30,to = 130,ylab = "D")

abline(v=c(40,60))

TAC_i = function(x){
  y=((alpha[1]*exp((alpha[2]-x)/alpha[3]))/
        (alpha[3]* exp((alpha[2]-x)+1)^2))
  return(y)
}

x=seq(40,60,length.out=1000)

y= TAC_i(x)

hist(y);abline(v=mean(y),col="red",lwd=2);abline(v=quantile(y,0.25),col="blue",lwd=2)

# Finding distance matrix
y_mat <- dist(y, method = 'euclidean')

Hierar_cl <- hclust(y_mat, method = "ward.D2")


fit <- cutree(Hierar_cl, k = 3)
df_c=data.frame(y,cluster=fit)
boxplot(df_c$y~df_c$cluster)
points(x=1:3,y=tapply(df_c$y,df_c$cluster,median),pch=16,col="red")
abline(h=mean(df_c$y),lwd=2,col="blue",lty=2)

plot(x=x,y=y,col=as.factor(df_c$cluster))

pesos=table(df_c$cluster)/1000
medias=tapply(df_c$y,df_c$cluster,mean)
mp = t(pesos)%*%medias
boxplot(df_c$y~df_c$cluster)
points(x=1:3,y=tapply(df_c$y,df_c$cluster,mean),pch=16,col="red")
abline(h=mean(df_c$y),lwd=2,col="blue",lty=2)
abline(h=mp,lwd=2,col="orange",lty=2)