rm(list=ls())
##
## Regressão Não Linear
## (slides da aula)
##
## Primeiro Exemplo:
## Crescimento de Planta: Modelo logístico
##

##
## Comandos preliminares
##
plantas <- data.frame(x = 1:7,
                      y = c(5, 13, 16, 23, 33, 38, 40)
                      )
plot(y ~ x, data = plantas, xlab = "Idade da planta (semanas)",
     ylab = "Altura (cm)", pch = 20, cex = 1.4, ylim = c(0, 45))

## modelo de regressão linear simples (linear)
m0 <- lm(y ~ x, data = plantas)
summary(m0)

## modelo logístico (não linear)
n0 <- nls(y ~ SSlogis(x, A, C, S), data = plantas)
summary(n0)
confint(n0)

## predições dos modelos
pred <- data.frame(x = seq(0, 10, length.out = 101))
pred$fit_m0 <- predict(m0, newdata = pred, interval = "confidence")
pred$fit_n0 <- data.frame(fit = predict(n0, newdata = pred))

## cálculos para bandas de confiança do modelo não linear
## (poderia usar uma função "pronta", mas vamos fazer na mão)
partials <- deriv(~A/(1 + exp(-(x - C)/S)),
                  namevec = c("A", "C", "S"),
                  function.arg = function(x, A, C, S) {
                      NULL
                  })
X <- do.call(what = partials,
             args = c(list(x = pred$x),
                      as.list(coef(n0))))
X <- attr(X, "gradient")

qtl <- qt(0.975, df = df.residual(n0))
se <- sqrt(diag(X %*% vcov(n0) %*% t(X)))
pred$fit_n0 <- transform(pred$fit_n0,
                  lwr = fit - qtl * se,
                  upr = fit + qtl * se)

rm(partials,X,qtl,se)

##
## Visualizações nos slides
##
par(mar=c(3.5, 3.5, 0,0), mgp=c(1.8, 0.8,0))
plot(y ~ x, data = plantas, xlab = "Idade da planta (semanas)",
     ylab = "Altura (cm)", pch = 20, cex = 1.4, xlim = range(pred$x), ylim = c(0, 50))
lines(fit_m0[,"fit"] ~ x, data = pred, col = "olivedrab4", lwd = 2)
lines(fit_n0[,"fit"] ~ x, data = pred, col = "red", lwd = 2)

par(mfrow = c(1, 3), mar=c(3.5, 3.5, 0,0), mgp=c(1.8, 0.8,0))
#layout(1)
plot(y ~ x, data = plantas, type = "n",
     xlab = "Idade da planta (semanas)",
     ylab = "Altura (cm)", pch = 20, cex = 1.4, xlim = range(pred$x),
     ylim = c(0, 50))
with(pred, matlines(x, y = fit_m0,
                    lty = c(1, 2, 2), lwd = c(2, 1, 1), col = "olivedrab4"))
with(pred, matlines(x, y = fit_n0,
                    lty = c(1, 2, 2), lwd = c(2, 1, 1), col = "red"))
points(y ~ x, data = plantas, pch = 19)
##
## Ajustes e bandas
##
plot(y ~ x, data = plantas, type = "n",
     xlab = "Idade da planta (semanas)",
     ylab = "Altura (cm)", pch = 20, cex = 1.4, xlim = range(pred$x),
     ylim = c(0, 50))
##names(pred)
with(pred,
     polygon(x = c(x, rev(x)),
             y = c(fit_m0[,"upr"], rev(fit_m0[,"lwr"])),
             col = "gray90", border = "olivedrab4"))
points(y ~ x, data = plantas, pch = 19)
lines(fit_m0[,"fit"] ~ x, data = pred, col = "olivedrab4")
plot(y ~ x, data = plantas, type = "n",
     xlab = "Idade da planta (semanas)",
     ylab = "Altura (cm)", pch = 20, cex = 1.4, xlim = range(pred$x),
     ylim = c(0, 50))
with(pred,
     polygon(x = c(x, rev(x)),
             y = c(fit_n0[,"upr"], rev(fit_n0[,"lwr"])),
             col = "gray90", border = "red"))
points(y ~ x, data = plantas, pch = 19)
lines(fit_n0[,"fit"] ~ x, data = pred, col = "red")


## Segundo Exemplo:
## Modelo de Michaelis Menten
##
## O modelo (hipérbole) tem origem em problemas de cinética enzimática sendo dado por:
## \[ \nu = \frac{V_{\max} \cdot [S]}{K_m + [S]} \]
## em que $\nu$ é a velocidade da reação (y),
## $[S]$ é a concentração do substrato (x),
## $V_{\max}$ é a velocidade máxima da reação e
## $K_m$ é a constante de Michaelis-Menten que quantifica a concentração de substrato necessária para atingir a metade do tempo máximo de reação.
## Os valores de
## $\nu$ e $[S]$ são observados
## $V_{\max}$ e $K_m$. parâmetros a estimar
##
rm(list=ls())

## Dados do exemplo
rquim <- data.frame(
        x = c(0.02, 0.06, 0.11, 0.22, 0.56, 1.10, 0.02, 0.06, 0.11, 0.22, 0.56, 1.10),
        y = c(47, 97, 123, 152, 191, 200, 76, 107, 139, 159, 201, 207)
)
par(mfrow=c(1,1), mar=c(3.5, 3, 0.3, 0.3), mgp=c(1.8,0.8,0))
plot(y ~ x, data=rquim,
     xlab="Concentração (ppm)", ylab=expression(Velocidade (c/min^2)),
     xaxp = c(0,1.1,11), yaxp = c(60,200,7), las=1)

## Ajuste do modelo de Michaelis-Menten
rquim.fit <- nls(y ~ b1*x/(x+b2), data = rquim,
                 start = list(b1=200, b2=0.2), trace = FALSE)
rquim.ndf <- data.frame(x = seq(0, 1.1, by=0.01))
rquim.ndf$y <-predict(rquim.fit, newdata=rquim.ndf, se.fit = TRUE)
with(rquim.ndf, lines(y ~ x, col=4))

rquim.fit <- nls(y ~ b1*x/(x+b2), data = rquim, start = list(b1=200, b2=0.2))
summary(rquim.fit)

## Modelo linearizado
b <- coef(rquim.fit)
b0nls <- 1/b[1]
b1nls <- b[2]/b[1]
rquim.lm <- lm(1/y ~ I(1/x), data=rquim)
bs <- coef(rquim.lm)
beta1.lm <- 1/bs[1]
beta2.lm <- bs[2]/bs[1]
##
ll.nls <- round(logLik(rquim.fit), 1)
ll.lm <- round(logLik(rquim.lm) + (-1-1)*sum(log(rquim$y)), 1)
##
par(mfrow=c(1,2), mar=c(3.2, 3.3, 3, 0.3), mgp=c(2,1,0))
plot(y ~ x, data=rquim,
     xlab="Concentração (ppm)", ylab=expression(Velocidade (c/min^2)),
     xaxp = c(0,1.1,11), yaxp = c(60,200,7), las=1, main="Escala Original")
with(rquim.ndf, lines(y ~ x, col=4))
MMfun <- function(x) beta1.lm*x/(x+beta2.lm)
curve(MMfun, from = 0, to = 1.1, col=2, add=TRUE)
text(0.8, 80, substitute(logLik == a, list(a = ll.nls)), col=4)
text(0.8, 70, substitute(logLik == a, list(a = ll.lm)), col=2)
##
plot(1/y ~ I(1/x), data=rquim, main="Escala transformada")
abline(bs, col=2)
abline(c(b0nls, b1nls), col=4)
legend("topleft", c("não linear","linearizado"), col=c(4,2), lty=1)