Comparar el modelo con el de Gompertz.

library(ggplot2)
set.seed(2023)
x<-rnorm(n = 1000,mean = 3,sd = 0.3)
hist(x);abline(v=mean(x),lwd=2,col="red")

y<-rbeta(n = 1000,shape1 =0.5 ,shape2 =0.5)
hist(y);abline(v=mean(y),lwd=2,col="red")

z<-rexp(n = 1000,rate = 1/3)
hist(z);abline(v=mean(z),lwd=2,col="red")

library(growthmodels)
growth <- gompertz(0:100, 10, 2, 0.5)
time <- 0:100
df_growth = data.frame(time,growth)
plot(df_growth$growth~df_growth$time)
text(x=80,y=4,paste("rs:",round(cor(df_growth,method="spearman"),3)))

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)
plot(medias~unique(t),pch=16)

plot(W~t,pch=16)
points(medias~unique(t),pch=16,col="red",cex=1)
abline(v=t,col="grey",lty=2)

medias_diff = diff(medias)
TAC = medias_diff/7

Minimos cuadrados no lineales

# Ajuste del modelo original
model_1 = nls(W~(a*t)/(b+t),start = list(a=1,b=1))
summary(model_1)
## 
## Formula: W ~ (a * t)/(b + t)
## 
## Parameters:
##   Estimate Std. Error t value Pr(>|t|)    
## a  8.41076    0.04266  197.14   <2e-16 ***
## b  3.69787    0.15152   24.41   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1256 on 34 degrees of freedom
## 
## Number of iterations to convergence: 5 
## Achieved convergence tolerance: 6.805e-06
plot(W~t,pch=16)
curve((coef(model_1)[1]*x)/(coef(model_1)[2]+x),add=T,col="blue")

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.704e-06
#Comparando con Gompertz
model_gompertz = nls(W ~ a * exp(-b * exp(-c * t)), start = list(a = 10, b = 0.1, c = 0.1))
summary(model_gompertz)
## 
## Formula: W ~ a * exp(-b * exp(-c * t))
## 
## Parameters:
##   Estimate Std. Error t value Pr(>|t|)    
## a  7.65368    0.08341  91.757  < 2e-16 ***
## b  3.36671    0.93555   3.599  0.00103 ** 
## c  0.32799    0.04359   7.524 1.19e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4431 on 33 degrees of freedom
## 
## Number of iterations to convergence: 15 
## Achieved convergence tolerance: 9.191e-06
plot(W ~ t, pch = 16)
curve(coef(model_gompertz)[1] * exp(-coef(model_gompertz)[2] * exp(-coef(model_gompertz)[3] * x)), add = TRUE, col = "blue")

plot(W~t,pch=16)
curve((coef(model_1)[1]*x)/(coef(model_1)[2]+x),add=T,col="blue")
curve(coef(model_2)[1] * (1 - exp(-coef(model_2)[2] * x)),add=T,col="red")
curve(coef(model_gompertz)[1] * exp(-coef(model_gompertz)[2] * exp(-coef(model_gompertz)[3] * x)), add = TRUE, col = "green")

predicciones_1 = fitted.values(model_1)
predicciones_2 = fitted.values(model_2)
predicciones_3= fitted.values(model_gompertz)
Metrics::rmse(actual = W,predicted = predicciones_1)
## [1] 0.1221091
Metrics::rmse(actual = W,predicted = predicciones_2)
## [1] 0.322335
Metrics::rmse(actual = W,predicted = predicciones_3)
## [1] 0.4242815

Viendo los valores obtenidos, se puede determinar que el modelo que más se ajusta en este caso es el 1. A pesar de utilizar el modelo Gompertz, el cual se utiliza comúnmente para describir el crecimiento de poblaciones biológicas y se ha aplicado a diversos campos como la epidemiología, el crecimiento de tumores y la economía. En este caso, según el valor RMSE, este no se adapta tanto como los demas, teniendo en cuenta que presenta el valor más alto

f= expression((a * t)/(b + t))
D(f,"t")
## a/(b + t) - (a * t)/(b + t)^2
curve((coef(model_1)[1]*coef(model_1)[2])/(coef(model_1)[2]+x^2),from = 0,to = 77,ylab = "D")
abline(v=c(7,14))

TAC_i = function(x){
  y=(coef(model_1)[1]*coef(model_1)[2])/(coef(model_1)[2]+x^2)
  return(y)
}

x=seq(0,14,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")

 
# Cutting tree by no. of clusters
fit <- cutree(Hierar_cl, 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)