#########################################################################
#### ####
#### Rotina Modelos Nao-Lineares ####
#### ####
####-----------------------------------------------------------------####
#### ####
#### + Logistico ####
#### + Gompertz ####
#### ####
#########################################################################
#########################################################################
#### Definicao Modelo ####
#########################################################################
#### Logistico
Logistico <- function(x,b0,b1,b2){b0/(1+exp(b1-b2*x))}
#### Gompertz
Gompertz <- function(x,b0,b1,b2){b0*exp(-exp(b1-b2*x))}
#########################################################################
#### Load Packages ####
#########################################################################
library(MASS)
library(lmtest)
library(car)
library(nortest)
library(manipulate) ## Permite criar gráficos "móveis"
library(readxl)
#########################################################################
#### Load Dados ####
#########################################################################
dados <- read.table("C:/Resultados R pos doc/MFF.txt",h=T, dec=".")
str(dados)
head(dados)
summary(dados)
###Determinação das estimativas iniciais, para cada modelo.
#########################################################################
#### Estimarivas Iniciais ####24
#########################################################################
#### Logistico
b0=200 #Assintota com chute inicial "bem errado"
b1=3.5 #Parâmetro de localização com chute inicial "bem errado"
b2=0.005 #Taxa de crescimento com chute inicial "bem errado"
plot(y~x,data=dados)
curve(Logistico(x,b0,b1,b2),add=T,col=2)
Chute_logist <- list() ## Lista em branco onde serão armazenados os chutes
### Programação específica do gráfico "movel"
manipulate({
  plot(y~x,data=dados)
  curve(Logistico(x,b0=b0,b1=b1,b2=b2),add=TRUE)
  Chute_logist <<-list(b0=b0,b1=b1,b2=b2)},
  b0=slider(0,350,initial=0),
  b1=slider(0,6.1,initial=0),
  b2=slider(0.001,0.9,initial=0.001)
)
Chute_logist ## onde estão armazenados os chutes
#### Gompertz
b0=200 #Assintota com chute inicial "bem errado"
b1=3.5 #Parâmetro de localização com chute inicial "bem errado"
b2=0.005 #Taxa de crescimento com chute inicial "bem errado"
plot(y~x,data=dados)
curve(Gompertz(x,b0,b1,b2),add=T,col=2)
Chute_gomp <- list() ## Lista em branco onde serão armazenados os chutes
### Programação específica do gráfico "movel"
manipulate({
  plot(y~x,data=dados)
  curve(Gompertz(x,b0=b0,b1=b1,b2=b2),add=TRUE)
  Chute_gomp <<-list(b0=b0,b1=b1,b2=b2)},
  b0=slider(0,350,initial=0),
  b1=slider(0,6.1,initial=0),
  b2=slider(0.001,0.9,initial=0.001))
Chute_gomp ## onde estão armazenados os chutes

#########################################################################
####                        Ajuste Modelo                            ####
#########################################################################
#### Logistico
Modelo_Logist <- nls(y~b0/(1+exp(b1-b2*x)), data=dados,start=Chute_logist)
ResumoModel_Logist <- summary(Modelo_Logist)
b0_logist <- summary(Modelo_Logist)$parameters[1,1]
b1_logist <- summary(Modelo_Logist)$parameters[2,1]
b2_logist <- summary(Modelo_Logist)$parameters[3,1]
IC_logist <- confint.default(Modelo_Logist)
#### Gompertz
Modelo_Gompertz <- nls(y~b0*exp(-exp(b1-b2*x)), data=dados,start=Chute_gomp)
ResumoModel_Gompertz <- summary(Modelo_Gompertz)
b0_gomp <- summary(Modelo_Gompertz)$parameters[1,1]
b1_gomp <- summary(Modelo_Gompertz)$parameters[2,1]
b2_gomp <- summary(Modelo_Gompertz)$parameters[3,1]
IC_gomp <- confint.default(Modelo_Gompertz, level = .95)

#########################################################################
#### Pressupostos Modelo ####
#########################################################################
#### Logitico
# Normalidade
SW_logist <- shapiro.test(residuals(Modelo_Logist))
LL_logist <- lillie.test(residuals(Modelo_Logist))
KS_logist <- ks.test(residuals(Modelo_Logist), "pnorm", mean(residuals(Modelo_Logist)), sd(residuals(Modelo_Logist)))
AD_logist <- ad.test(residuals(Modelo_Logist))
# Homogeneidade de variacia
gradiente_mod_logist <- attr(Modelo_Logist$m$fitted(),"gradient") # obtem matriz gradiente, sobre a qual será feito o diagnóstico
m0_mod_logist <- lm(y~-1+gradiente_mod_logist, data=dados) # passando a m atriz gradiente para a lm(), importante remover intercepto (-1)
BP_logist <- bptest(m0_mod_logist) # teste de Breusch-Pagan (homogeneidade)
BP_logist
BAR_logist <- bartlett.test(residuals(Modelo_Logist)~x,data=dados)
BAR_logist
# Idependencia
DW_logist <- durbinWatsonTest(m0_mod_logist) # teste de DW (independência)
DW_logist
# Não linearidade
## Função deriv3 retorna a matriz hessiana necessária para calcular a não linearidade
## pelo método de Bates e Watts
prmt_mod_logist <- deriv3(~b0/(1+exp(b1-b2*x)),c("b0","b1","b2"),function
                          (x,b0,b1,b2)NULL)
nls_logist <- nls(y~prmt_mod_logist(x,b0,b1,b2),data=dados, start=Chute_logist)
rms.curv(nls_logist)
parm_R1_logist <- rms.curv(nls_logist)$pe
intrinseca_R1_logist <- rms.curv(nls_logist)$ic
#### Gompertz
# Normalidade
####### Analisando os resíduos
SW_gomp <- shapiro.test(residuals(Modelo_Gompertz))
LL_gomp <- lillie.test(residuals(Modelo_Gompertz))
KS_gomp <- ks.test(residuals(Modelo_Gompertz), "pnorm", mean(residuals(Modelo_Gompertz)), sd(residuals(Modelo_Gompertz)))
AD_gomp <- ad.test(residuals(Modelo_Gompertz))
# Homogeneidade
gradiente_mod_gomp <- attr(Modelo_Gompertz$m$fitted(),"gradient") # obtem matriz gradiente, sobre a qual será feito o diagnóstico
m0_mod_gomp <- lm(y~-1+gradiente_mod_gomp, data=dados) # passando a matriz gradiente para a lm(), importante remover intercepto (-1)
BP_gomp <- bptest(m0_mod_gomp) # teste de Breusch-Pagan (homogeneidade)
BP_gomp
BAR_gomp <- bartlett.test(residuals(Modelo_Gompertz)~x,data=dados)
BAR_gomp
# Independencia
DW_gomp <- durbinWatsonTest(m0_mod_gomp) # teste de DW (independência)
DW_gomp
# Não linearidade
## Função deriv3 retorna a matriz hessiana necessária para calcular a não
linearidade
## pelo método de Bates e Watts
prmt_mod_gomp <- deriv3(~b0*exp(-exp(b1-b2*x)),c("b0","b1","b2"),function
                        (x,b0,b1,b2)NULL)
nls_gomp <- nls(y~prmt_mod_gomp (x,b0,b1,b2),data=dados, start=Chute_gomp)
rms.curv(nls_gomp)
parm_R1_gomp <- rms.curv(nls_gomp)$pe
intrinseca_R1_gomp <- rms.curv(nls_gomp)$ic

#########################################################################
#### Pontos Criticos ####
#########################################################################
# Logistico
PIX_logist <- summary(Modelo_Logist)$parameters[2,1]/summary(Modelo_Logist)$parameters[3,1]
PIY_logist <- summary(Modelo_Logist)$parameters[1,1]/2

PAMX_logist <- (summary(Modelo_Logist)$parameters[2,1]-log(2+sqrt(3)))/summary(Modelo_Logist)$parameters[3,1]
PAMY_logist <- summary(Modelo_Logist)$parameters[1,1]/(3+sqrt(3))
PDMX_logist <- (summary(Modelo_Logist)$parameters[2,1]-log(2-sqrt(3)))/summary(Modelo_Logist)$parameters[3,1]
PDMY_logist <- summary(Modelo_Logist)$parameters[1,1]/(3-sqrt(3))
PDAX_logist <- (summary(Modelo_Logist)$parameters[2,1]-log(5-2*sqrt(6)))/summary(Modelo_Logist)$parameters[3,1]
PDAY_logist <- summary(Modelo_Logist)$parameters[1,1]/(2*(3-sqrt(6)))

# Gompertz
PIX_gomp <- summary(Modelo_Gompertz)$parameters[2,1]/summary(Modelo_Gompertz)$parameters[3,1]
PIY_gomp <- summary(Modelo_Gompertz)$parameters[1,1]/2
PAMX_gomp <- (summary(Modelo_Gompertz)$parameters[2,1]-log(2+sqrt(3)))/summary(Modelo_Gompertz)$parameters[3,1]
PAMY_gomp <- summary(Modelo_Gompertz)$parameters[1,1]/(3+sqrt(3))
PDMX_gomp <- (summary(Modelo_Gompertz)$parameters[2,1]-log(2-sqrt(3)))/summary(Modelo_Gompertz)$parameters[3,1]
PDMY_gomp <- summary(Modelo_Gompertz)$parameters[1,1]/(3-sqrt(3))
PDAX_gomp <- (summary(Modelo_Gompertz)$parameters[2,1]-log(5-2*sqrt(6)))/summary(Modelo_Gompertz)$parameters[3,1]
PDAY_gomp <- summary(Modelo_Gompertz)$parameters[1,1]/(2*(3-sqrt(6)))

#########################################################################
#### Qualidade Modelo ####
#########################################################################
# Logistico
SQR_logist <- sum((dados$y-fitted(Modelo_Logist))^2) ; SQR_logist
SQT_logist <- sum(dados$y*dados$y) - (sum(dados$y)^2)/length(dados$y) ;
SQT_logist
R2_logist <- 1-SQR_logist/SQT_logist ; R2_logist
R2adj_logist <- 1 - (((1-R2_logist)*(dim(dados)[1]-1)/(((summary(Modelo_Logist))$df)[2]))); R2adj_logist
# df total=9 df residuo=6
QMe_logist <- SQR_logist/(((summary(Modelo_Logist))$df)[2]); QMe_logist
DPR_logist <- sqrt(QMe_logist) ; DPR_logist
# Gompertz
SQR_gomp <- sum((dados$y-fitted(Modelo_Gompertz))^2) ; SQR_gomp
SQT_gomp <- sum(dados$y*dados$y) - (sum(dados$y)^2)/length(dados$y) ; SQT_gomp
R2_gomp <- 1-SQR_gomp/SQT_gomp ; R2_gomp
R2adj_gomp <- 1 - (((1-R2_gomp)*(dim(dados)[1]-1)/(((summary(Modelo_Gompertz))$df)[2]))); R2adj_gomp
# df total=9 df residuo=6
QMe_gomp <- SQR_gomp/(((summary(Modelo_Gompertz))$df)[2]); QMe_gomp
DPR_gomp <- sqrt(QMe_gomp) ; DPR_gomp
#########################################################################
#### Arquivo Unico ####
#########################################################################
Resultado <- matrix(NA, nrow = nrow(dados)+1, ncol = 22)
Resultado[1,c(1:6)] <- c("X", "Y", "Predito_log", "Residuos_Log", "Predito_Gomp", "Residuos_Gomp")
Resultado[c(2:(nrow(dados)+1)),1] <- as.vector(dados$x)
Resultado[c(2:(nrow(dados)+1)),2] <- as.vector(dados$y)
Resultado[c(2:(nrow(dados)+1)),3] <- as.vector(predict(Modelo_Logist))
Resultado[c(2:(nrow(dados)+1)),4] <- as.vector(ResumoModel_Logist$residuals)
Resultado[c(2:(nrow(dados)+1)),5] <- as.vector(predict(Modelo_Gompertz))
Resultado[c(2:(nrow(dados)+1)),6] <- as.vector(ResumoModel_Gompertz$residuals)
Resultado[1, 9] <- "Modelo Logistico"
Resultado[4:6, 9:12] <- ResumoModel_Logist$coefficients
Resultado[4:6, 8] <- rownames(ResumoModel_Logist$coefficients)
Resultado[3, 9:12] <- colnames(ResumoModel_Logist$coefficients)
Resultado[10:12, 9:10] <- IC_logist
Resultado[10:12, 8] <- rownames(IC_logist)
Resultado[9, 9:10] <- colnames(IC_logist)
Resultado[15, 9] <- "Pressupostos"
Resultado[17, 9] <- "Normalidade"
Resultado[18, 9:11] <- c("Teste", "Estatistica", "Valor-p")
Resultado[19, 9:11] <- c(SW_logist$method, SW_logist$statistic, SW_logist$p.value)
Resultado[20, 9:11] <- c(LL_logist$method, LL_logist$statistic, LL_logist$p.value)
Resultado[21, 9:11] <- c(KS_logist$method, KS_logist$statistic, KS_logist$p.value)
Resultado[22, 9:11] <- c(AD_logist$method, AD_logist$statistic, AD_logist$p.value)
Resultado[25, 9] <- "Homogeneidade"
Resultado[26, 9:11] <- c(BP_logist$method, BP_logist$statistic, BP_logist$p.value)
Resultado[27, 9:11] <- c(BAR_logist$method, BAR_logist$statistic, BAR_logist$p.value)
Resultado[29, 9] <- "Independencia"
Resultado[30, 9:11] <- c(BP_logist$method, BP_logist$statistic, BP_logist
                         $p.value)
Resultado[32, 9] <- "Não-Linearidade"
Resultado[33, 10:11] <- c("Parametrica", "Intrinseca")
Resultado[34, 10:11] <- c(parm_R1_logist, intrinseca_R1_logist)
Resultado[1, 16] <- "Modelo Gompertz"
Resultado[4:6, 16:19] <- ResumoModel_Gompertz$coefficients
Resultado[4:6, 15] <- rownames(ResumoModel_Gompertz$coefficients)
Resultado[3, 16:19] <- colnames(ResumoModel_Gompertz$coefficients)
Resultado[10:12, 16:17] <- IC_gomp
Resultado[10:12, 15] <- rownames(IC_gomp)
Resultado[9, 16:17] <- colnames(IC_gomp)
Resultado[15, 16] <- "Pressupostos"
Resultado[17, 16] <- "Normalidade"
Resultado[18, 16:18] <- c("Teste", "Estatistica", "Valor-p")
Resultado[19, 16:18] <- c(SW_gomp$method, SW_gomp$statistic, SW_gomp$p.value)
Resultado[20, 16:18] <- c(LL_gomp$method, LL_gomp$statistic, LL_gomp$p.value)
Resultado[21, 16:18] <- c(KS_gomp$method, KS_gomp$statistic, KS_gomp$p.value)
Resultado[22, 16:18] <- c(AD_gomp$method, AD_gomp$statistic, AD_gomp$p.value)
Resultado[25, 16] <- "Homogeneidade"
Resultado[26, 16:18] <- c(BP_gomp$method, BP_gomp$statistic, BP_gomp$p.value)
Resultado[27, 16:18] <- c(BAR_gomp$method, BAR_gomp$statistic, BAR_gomp$p.value)
Resultado[29, 16] <- "Independencia"
Resultado[30, 16:18] <- c(BP_gomp$method, BP_gomp$statistic, BP_gomp$p.value)
Resultado[32, 16] <- "Não-Linearidade"
Resultado[33, 17:18] <- c("Parametrica", "Intrinseca")
Resultado[34, 17:18] <- c(parm_R1_gomp, intrinseca_R1_gomp)
Resultado[40, 9:11] <- c("", "Logistico", "Gompertz")
Resultado[41, 9:11] <- c("PIX", PIX_logist, PIX_gomp)
Resultado[42, 9:11] <- c("PIY", PIY_gomp, PIY_gomp)
Resultado[43, 9:11] <- c("PAMX", PAMX_logist, PAMX_gomp)
Resultado[44, 9:11] <- c("PAMY", PAMY_logist, PAMY_gomp)
Resultado[45, 9:11] <- c("PDMX", PDMX_logist, PDMX_gomp)
Resultado[46, 9:11] <- c("PDMY", PDMY_logist, PDMY_gomp)
Resultado[47, 9:11] <- c("PDAX", PDAX_logist, PDAX_gomp)
Resultado[48, 9:11] <- c("PDAY", PDAY_logist, PDAY_gomp)
Resultado[49, 9:11] <- c("SQR", SQR_logist, SQR_gomp)
Resultado[50, 9:11] <- c("SQT", SQT_logist, SQT_gomp)
Resultado[51, 9:11] <- c("QMe", QMe_logist, QMe_gomp)
Resultado[52, 9:11] <- c("DPR", DPR_logist, DPR_gomp)
Resultado[53, 9:11] <- c("R2", R2_logist, R2_gomp)
Resultado[54, 9:11] <- c("R2a", R2adj_logist, R2adj_gomp)

write.table(Resultado, file = "./TesteMFF.txt")
write.csv(Resultado, file = "./TesteMFF.csv")