rm(list=ls())
par.ori <- par(no.readonly = TRUE)
##
## Introdução a regressão não linear
##

##
## Parte I : uma análise "básica" de regressão não linear
##

## Exemplo 1: Y = \beta_0 + \beta_1 * 2^{-x/\beta_2} + \epsilon

##
## Definindo a função
##
fun1 <- function(x, beta0, beta1, beta2) {
  beta0 + beta1 * 2^(-x / beta2)
}
fun1(6, 100, 80, 25)
fun1(1:10, 100, 80, 25)

##
## Visualizando a função e a interpretação dos parâmetros
##
curve(fun1(x, 80, 40, 30), from = 0, to = 150, ylim = c(70, 120))
abline(h = 80, lty=3)
text(0, 80, expression(beta[0]), pos = 4)
##
arrows(0, 80, 0, 80+40, lty = 3, code = 3, length = 0.1, col = "red")
text(0,100, expression(beta[1]), pos = 4, col = "red")
##
segments(0, 100, 30, fun1(30, 80, 40, 30), lty = 3, col = "blue")
arrows(30, fun1(30, 80, 40, 30), 30, 70, lty = 3, code = 3, length = 0.1, col = "blue")
text(30, 70, expression(beta[2]), pos = 4, col = "blue")

##
## Simulando um conjunto de dados
##
x <- sort(sample(1:150, 20))
y <- fun1(x, 80, 40, 30) + rnorm(20, sd = 2)
points(x, y, pch = 20, col = "black")

##
## Ajustando o modelo não linear aos dados simulados
## (várias formas de obter o ajuste
##
plot(x, y, pch = 20, ylim = c(70, 120))

## valores iniciais por inspeção visual do gráfico
ini <- c(beta0 = 70, beta1 = 50, beta2 = 40)

ajuste <- nls(y ~ beta0 + beta1*2^(-x/beta2),
              start = ini, trace = TRUE)

## Examinando e avaliando o ajuste 
ajuste
coef(ajuste)
summary(ajuste)
c(lLik = logLik(ajuste), AIC = AIC(ajuste), BIC = BIC(ajuste))

##
vcov(ajuste)

## IC para parâmetros
confint.default(ajuste)
cbind(
    coef(ajuste) - 1.96 * sqrt(diag(vcov(ajuste))),
    coef(ajuste) + 1.96 * sqrt(diag(vcov(ajuste)))
)

confint(ajuste)
confint.default(ajuste)

plot(profile(ajuste, which = "beta2"))
par(mfrow=c(1,3), mar=c(3,3,0.5,0.5), mgp = c(1.8, 0.8, 0))
plot(profile(ajuste))
par(par.ori)

## Predição:
## i) nos pontos observados
predict(ajuste)
plot(x, y, pch = 20, ylim = c(70, 120))
lines(x, predict(ajuste), type = "b", col = "red", pch = 19)

## ii) em um grid de valores de x
pred <- data.frame(x =  seq(0, 150, length.out = 151))
pred$y <- predict(ajuste, newdata = pred)
pred$y

plot(x, y, pch = 20, ylim = c(70, 120))
with(pred, lines(x, y, col = "blue", type = "b", pch=20))
plot(x, y, pch = 20, ylim = c(70, 120))
with(pred, lines(x, y, col = "blue", type = "l"))

## Bandas de confiança
pred$yic <- investr::predFit(ajuste, newdata = pred, interval = "confidence")
head(pred$yic)
with(pred, matplot(x, yic, type = "l", col = "blue", lty = c(1,2,2), add = TRUE))
## Bandas de predição
pred$yip <- investr::predFit(ajuste, newdata = pred, interval = "prediction")
head(pred$yip)
with(pred, matplot(x, yip, type = "l", col = "blue", lty = c(1,3,3), add = TRUE))

##
## Fim da parte I
##

##
## Parte II
##
## Algumas opções adicionais para o ajuste do modelo
##
fits <- list()

## Ajuste "básico"
fits$fit1 <- nls(y ~ beta0 + beta1*2^(-x/beta2),
            start = ini, trace = TRUE)
fits$fit1

## Ajuste usando a função definida
fits$fit2 <- nls(y ~ fun1(x, beta0, beta1, beta2),
                 start = ini, trace = TRUE)
fits$fit2

## Ajuste usando a função definida com gradiente
## derivadas:
## d/d beta0 = 1
## d/d beta1 = 2^(-x/beta2)
## d/d beta2 = (log(2) * x beta1 * 2^(-x/beta2)) / beta2^2

fun1g <- function(x, beta0, beta1, beta2) {
    termo <- 2^(-x / beta2)
    GR <- cbind(1, 2^(-x / beta2), log(2) * x * beta1 * termo / (beta2^2))
    dimnames(GR) <- list(NULL, c("beta0", "beta1", "beta2"))
    res <- beta0 + beta1 * termo
    attr(res, "gradient") <- GR
    return(res)
}

fits$fit3 <- nls(y ~ fun1g(x, beta0, beta1, beta2),
                 start = ini, trace = TRUE)
fits$fit3

## Ajuste usando a função definida com gradiente fornecido por deriv() 
fun1gd <- function(x, beta0, beta1, beta2) {
    termo <- 2^(-x / beta2)
    GR <- deriv(y ~ beta0 + beta1 * 2^(-x/beta2),
                namevec = c("beta0", "beta1", "beta2"),
                function.arg = c("x", "beta0", "beta1", "beta2")
                )
    res <- beta0 + beta1 * termo
    attr(res, "gradient") <- GR
    return(res)
}
fits$fit4 <- nls(y ~ fun1g(x, beta0, beta1, beta2),
                 start = ini, trace = TRUE)
fits$fit4


## Definindo função SelfStart

initfun1 <- function(mCall, data, LHS) {
  xy <- sortedXyData(mCall[["x"]], LHS, data)
  x <- xy[,"x"]
  y <- xy[,"y"]
  ini <- c(beta0 = 0.9 * min(y),
           beta1 = 1.1*(max(y) - min(y)),
           beta2 = (max(x) - min(x))/2)
  return(ini)
}
SSfun1 <- selfStart(model = fun1,
                    initial = initfun1,
                    parameters = c("beta0", "beta1", "beta2"))
fits$fit5 <- nls(y ~ SSfun1(x, beta0, beta1, beta2), trace = TRUE)
fits$fit5


## Ajuste aproveitando a linearidade em beta0 e beta1
fits$fit6 <- nls(y ~ cbind(1, 2^(-x/beta2)),
                 start = list(beta2 = 40), trace = TRUE, algorithm = "plinear")
fits$fit6   ## notar a ordem dos parâmetros