##rm(list=ls())
##
## Modelos não Lineares
## Exemplo:
## Ajustando Isoterma de Langmuir não linear no R
##

## Modelo:  E[Y] = Q*b*X/(1+b*X)
##
dt <- data.frame(x = c(10, 30, 50, 70, 100, 125),
                 y = c(155, 250, 270, 330, 320, 323)
                 )
plot(y ~ x, data = dt)
plot(y ~ x, data = dt, xlim = c(0, 150), ylim = c(100, 350))


## não converge com "má escolha" de valores iniciais
try(fit.bad <- nls(y ~ Q*b*x/(1+b*x),  data = dt,
               start = list(Q = 1, b = 0.5), algorithm = "port"))

## mas converge com
fit <- nls(y ~ Q*b*x/(1+b*x),  data = dt,
                 start = list(Q = 300, b = 1), algorithm = "port")
fit
coef(fit)

## Visualizando o ajuste
ff <- function(x, Q, b) Q*b*x/(1+b*x)
curve(ff(x, Q = coef(fit)["Q"], b = coef(fit)["b"]), add = TRUE, col = "blue", lwd = 2)

##
## Algumas opções para lidar com a especificação dos valores iniciais 
##

## a) Usando valores iniciais razoáveis (considerando a interpretação dos parâmetros)

## usando mean(dt$y) como valor inicial para Q funciona aqui
nls(y ~ Q*b*x/(1+b*x), data = dt, start = list(Q = mean(dt$y), b = 0.5))

## usando max(dt$y) como valor inicial para Q tb funciona aqui
nls(y ~ Q*b*x/(1+b*x), data = dt, start = list(Q = max(dt$y), b = 0.5))

## Se igualarmos Y a Q*b de forma que cancele e x a mean(dt$x) 
## podemos usar 1 - 1/mean(dt$x) como valor inicial para b
nls(y ~ Q*b*x/(1+b*x), data = dt, start = with(dt, list(Q = max(y), b = 1 - 1/mean(x))))


## b) fazendo um rpanel para selecionar valores iniciais visualmente

library(rpanel)

rp.langmuir <- function(data) {
    ultimos <- c(Q = NA, b = NA) ## vai guardar os últimos valores usados no gráfico
    plotlangmuir <- function(panel) {
        with(panel,
             plot(1, type = "n", xlim = c(0, 150), ylim = c(100, 350),
                  xlab = "Concentração (mg/L)", ylab = "Adsorção (mg/g)",
                  main = paste("Isoterma de Langmuir: Q =", round(Q, 2), ", b =", round(b, 2))
                  ))
        with(data, points(x, y, pch = 20))
        with(panel, curve(Q*b*x/(1+b*x), add = TRUE, col = "blue", lwd = 2))
        ultimos <<- c(Q = panel$Q, b = panel$b)
        panel
    }
    panel <- rp.control(title = "Ajuste Isoterma de Langmuir")
    rp.slider(panel, variable = Q, from = 1, to = 500, resolution = 1,
              initval = 300, showvalue = TRUE, action = plotlangmuir)
    rp.slider(panel, variable = b, from = 0.01, to = 1,
              resolution = 0.01, initval = 0.1, showvalue = TRUE, action = plotlangmuir)
    rp.do(panel, action = plotlangmuir)
    rp.button(panel, title = "Concluir", quitbutton = TRUE, action = I)
    rp.block(panel) ## espera a janela ser fechada (botão "Concluir" ou X)
    ultimos
}
ini <- rp.langmuir(dt)
ini

##
## c) Ajustando uma aproximação linear para obter valores iniciais
##    (note diferença no termo de erro: aqui usamos 1/y e não y)

## aproximação linearizando
## (1/y) = (1/Q) + (1/Qb)(1/x)

## E[y^*] = \beta_0 + \beta_1 x^*
## y^* = 1/y, x^* = 1/x,
## \beta_0 = 1/Q, \beta_1 = 1/(Qb)
## Q = 1/\beta_0, b = \beta_0/\beta_1
plot(I(1/y) ~ I(1/x), data = dt)
fitLin <- lm(I(1/y) ~ I(1/x), data = dt)
abline(fitLin, col = "blue", lwd = 2)
QLin <- 1/coef(fitLin)[1]               ## Q
bLin <- coef(fitLin)[1]/coef(fitLin)[2] ## b
nls(y ~ Q*b*x/(1+b*x), data = dt, start = list(Q = QLin, b = bLin))

##
## d) Definindo uma função selfStart para obter valores iniciais razoáveis automaticamente
##    (vamos usar a aproximação linear do item (c) para obter os valores iniciais)
##

## initial: recebe a fórmula (mCall), a resposta (LHS) e os dados, e
## retorna valores iniciais nomeados para os parâmetros do modelo
langmuirInit <- function(mCall, LHS, data, ...) {
    xy <- sortedXyData(mCall[["x"]], LHS, data)
    fitLin <- lm(I(1/y) ~ I(1/x), data = xy)
    b0 <- coef(fitLin)
    Q <- 1/b0[1]
    b <- b0[1]/b0[2]
    value <- c(Q, b)
    names(value) <- mCall[c("Q", "b")]
    value
}

SSlangmuir <- selfStart(
    function(x, Q, b) Q*b*x/(1+b*x),
    initial = langmuirInit,
    parameters = c("Q", "b")
)

## agora nls() obtém os valores iniciais sozinho, sem precisarmos informá-los
fit.ss <- nls(y ~ SSlangmuir(x, Q, b), data = dt)
summary(fit.ss)

##
## e) reparametrizando a função
##
## Q*b*x/(1+b*x) = Q*x/(x + 1/b) = Q*x/(x + K), com K = 1/b
##
## Esta é exatamente a forma do modelo de Michaelis-Menten (Q = Vmax,
## K = Km) -- ver observação junto à definição do modelo, no início
## do arquivo.
##
## Nessa parametrização K tem interpretação direta: é a concentração x
## na qual a adsorção atinge metade do máximo (y = Q/2). Além de mais
## fácil de interpretar, valores iniciais também ficam mais intuitivos
## (Q pelo valor máximo observado de y, K por um valor "típico" de x).

fit.K <- nls(y ~ Q*x/(x+K), data = dt, start = list(Q = max(dt$y), K = mean(dt$x)))
fit.K
coef(fit.K)

## voltando à parametrização original, se necessário
b.K <- 1/coef(fit.K)["K"]
b.K

## as duas parametrizações dão o mesmo ajuste (mesma soma de quadrados)
## e a mesma magnitude de correlação entre as estimativas dos parâmetros
## (o sinal se inverte porque K = 1/b)
cov2cor(vcov(fit))    ## Q, b
cov2cor(vcov(fit.K))  ## Q, K

##
## f) Usando a opção de "plinear"
##

## plinear: como o lado direito é linear em Q*b, basta escrever a fórmula
## como o termo que multiplica esse parâmetro linear (aqui, x/(1+b*x)) e
## deixar que o algoritmo "absorva" Q -- assim só precisamos de valor
## inicial para b. O ajuste retorna coeficientes "b" e ".lin", com
## Q = .lin/b.

fit.pl <- nls(y ~ x/(1+b*x), data = dt, start = list(b = 0.5), alg = "plinear")
(pars <- with(as.list(coef(fit.pl)), list(Q = .lin/b, b = b)))

## fit.pl já tem o ajuste.
## Entretanto, para obter uma chamada completa em termos de "Q" e "b" (ou invés de b and .lin)
## pode-se rodar nls() original com ajuste plinear como valor inicial
nls(y ~ Q*b*x/(1+b*x), data = dt, start = pars)


## g) Usar outro otimizador: por exemplo, optim()

rss <- function(par, data){
    ## data deve ser data-frame com variáveis x e y
    ## par[1] = Q, par[2] = b
    with(data, sum((y - par[2]*par[1]*x/(1+par[2]*x))^2))
}
optim(c(300, 1), rss, data = dt)[1:2]
optim(c(1, 0.5), rss, data = dt)[1:2]

##
## Comentário final
##
## Esta é a mesma função (hipérbole retangular) do modelo de
## Michaelis-Menten, E[V] = Vmax*S/(Km+S): basta escrever
## Q*b*X/(1+b*X) = Q*X/(X+K), com K = 1/b, para ver a correspondência
## Q <-> Vmax, X <-> S, K <-> Km.

## Vários domínios do conhecimento podem usar uma mesma função para modelar fenômenos diferentes, e cada domínio pode ter uma interpretação diferente para os parâmetros. Por isso, é importante conhecer a função de forma genérica e a interpretação.