# Poner el directorio de trabajo con la función setwd()

# Carga y exploración de datos ####
delfines <- read.table("Delfines.txt", header=T) 

plot(delfines$anios, delfines$talla, pch=16, col="gray", ylim=c(0,300), xlim=c(0,45), ylab="Largo (cm)", xlab="Edad (años)", lwd=2)


# INFERENCIA MULTIMODELO SEGÚN FIGURA 2 DE LA AYUDA DIDÁCTICA

# (1) Planteamos múltiples modelos #####
# Ver Caja 2 en la ayuda didáctica

# vonBM = Modelo von Bertalanffy
# GptzM = Modelo Gompertz
# LogsM = Modelo Logístico


# (2) Estimamos los parámetros de cada modelo ####

# Los tres modelos son no lineales en sus parámetros y por lo tanto los ajustamos con la función nls
# Para saber más acerca de esta función podemos escribir en r 
?nls

# a1 - Modelo von Bertalanffy
vonBM <- nls(delfines$talla ~ a*(1-exp(-b*delfines$anios)), start=list(a=max(delfines$talla), b=0.5),  data=delfines)
summary(vonBM) 

# a2 - Modelo Gompertz 
GptzM <- nls(delfines$talla ~ a*exp(-exp(c-(b*delfines$anios))), start=list(a=max(delfines$talla), b=0.5, c=0),  data=delfines)
summary(GptzM)

# a3 - Modelo Logístico 
LogsM <- nls(delfines$talla ~ a/(1+exp(c-(b*delfines$anios))), 
          start=list(a=max(delfines$talla), b=0.5, c=0),  data=delfines)
summary(LogsM)


# (3) Verificamos los supuestos ####

resvonBM <- resid(vonBM) # residuos del modelo von Bertalanffy
predvonBM <- fitted(vonBM) # predichos del modelo von Bertalanffy

resGptzM <- resid(GptzM) # residuos del modelo Gompertz
predGptzM <- fitted(GptzM) # predichos del modelo Gompertz

resLogsM <- resid(LogsM) # residuos del modelo Logístico
predLogsM <- fitted(LogsM) # predichos del modelo Logístico


# 3.1 - Evaluamos la homocedasticidad, linealidad e independecia 

# Gráfico de residuos vs. predichos
op <- par(mfrow=c(1,3))
plot(predvonBM, resvonBM, ylab="Residuos von Bertalanfy", xlab="Predichos von Bertalanfy")
abline(a=0,b=0, lw=2)
plot(predGptzM, resGptzM, ylab="Residuos Gompertz", xlab="Predichos Gompertz")
abline(a=0,b=0, lw=2)
plot(predLogsM, resLogsM, ylab="Residuos Logistico", xlab="Predichos Logistico")
abline(a=0,b=0, lw=2)
par(op)

# Gráfico de residuos vs. edad 
op <- par(mfrow=c(1,3))
plot(delfines$anios, resvonBM,ylab="Residuos von Bertalanfy", xlab="Edad (años")
abline(a=0, b=0, lw=3)
plot(delfines$anios, resGptzM, ylab="Residuos Gompertz", xlab="Edad (años)")
abline(a=0, b=0, lw=3)
plot(delfines$anios, resLogsM, ylab="Residuos Logístico", xlab="Edad (años)")
abline(a=0, b=0, lw=3)
par(op)


# 3.2 - Evaluamos la normalidad 

# Histogramas y diagramas de caja y bigotes
op <- par(mfrow=c(3,2))
hist(resvonBM, col="grey", xlab="Residuos", ylab="Frecuencia", main="von Bertalanfy")
boxplot(resvonBM, bty="l", range=1.5, col="grey", horizontal=T, xlab="Residuos")
hist(resGptzM, col="grey", xlab="Residuos", ylab="Frecuencia", main="Gompertz")
boxplot(resGptzM, bty="l", range=1.5, col="grey", horizontal=T, xlab="Residuos")
hist(resLogsM, col="grey", xlab="Residuos", ylab="Frecuencia",main="Logístico")
boxplot(resLogsM, bty="l", range=1.5, col="grey", horizontal=T, xlab="Residuos")
par(op)

# Gráficos cuantil-cuantil
op <- par(las=1, mfrow=c(1,3))
qqnorm(resvonBM, main="von Bertalanfy", ylab="Cuantiles observados", xlab="Cuantiles teóricos") 
qqline(resvonBM)
qqnorm(resGptzM, main="Gompertz", ylab="Cuantiles observados", xlab="Cuantiles teóricos") 
qqline(resGptzM)
qqnorm(resLogsM, main="Logístico", ylab="Cuantiles observados", xlab="Cuantiles teóricos") 
qqline(resLogsM)
par(op)

# Asimetría y curtósis en  los residuos
library(moments)
skewness(resvonBM) # la asimetría de una distribución normal es cero
kurtosis(resvonBM) # la curtósis de una distribución normal es tres

skewness(resGptzM) 
kurtosis(resGptzM) 

skewness(resLogsM)
kurtosis(resLogsM) 

# Test de Kolmogorov-Smirnov para evaluar normalidad
ks.test(resvonBM,"pnorm" ,mean(resvonBM), sd(resvonBM))
ks.test(resvonBM,"pnorm", mean(resGptzM), sd(resGptzM))
ks.test(resvonBM,"pnorm", mean(resLogsM), sd(resLogsM))

# Modificacion de Lillefors al test de Kolmogorov-Smirnov
library(nortest) 
lillie.test(resvonBM)
lillie.test(resGptzM)
lillie.test(resLogsM)

# Test de Shapiro-Wilk para evaluar normalidad
shapiro.test(resvonBM)
shapiro.test(resGptzM)
shapiro.test(resLogsM)

# Graficamos los valores observados vs valores predichos
op <- par(mfrow=c(1,3))
plot(delfines$talla, predvonBM, ylab="Observados", xlab="Predichos", main="von Bertalanfy")
abline(a=0,b=1, lw=2)
plot(delfines$talla, predGptzM, ylab="Observados", xlab="Predichos", main="Gompertz")
abline(a=0,b=1, lw=2)
plot(delfines$talla, predLogsM, ylab="Observados", xlab="Predichos", main="Logístico")
abline(a=0,b=1, lw=2)
par(op)


# (4) Ordenamos los modelos según AIC ####
AIC(LogsM, GptzM, vonBM)

library("bbmle") 
AICtab(LogsM, GptzM, vonBM, weights = T, delta = TRUE, base=T, sort = TRUE, logLik = TRUE)


# (5) Presentamos los modelos ####
plot(delfines$anios, delfines$talla, pch=16, col="gray", ylim=c(0,300), xlim=c(0,45), ylab="Largo (cm)", xlab="Edad (años)", lwd=2)
curve(coef(vonBM)[1] * (1 - exp(-coef(vonBM)[2]  * x)), ylim=c(0,300), xlim=c(0,45), lwd=1, add=T) # von Bertalanffy 
curve(coef(GptzM)[1] * exp(-exp(coef(GptzM)[3] - (coef(GptzM)[2] * x))), ylim=c(0,300), lwd=1, xlim=c(0,45), lty=2, add=T) # Gompertz 
curve(coef(LogsM)[1] / (1 + exp(coef(LogsM)[3] - (coef(LogsM)[2] * x))), ylim=c(0,300), xlim=c(0,45), lwd=1, lty=3, add=T) # Logístico
legend("bottomright",legend=c("Modelo von Bertalanffy", "Modelo Gompertz", "Modelo Logístico"),  
       lty=c(1, 2, 3), bty = "n", x.intersp=1, y.intersp=1.15, cex = 0.8)