rm(list=ls())
##
## Outros dados para o modelo de Michaelis Menten 
## (do arquivo de trabalhos de Est Comp)

##
## Visualizando o Modelo
##
par(mar = c(3, 3, 1, 1), mgp = c(2, 0.5, 0))
funMM <- function(x, Vmax, Km) Vmax*x/(Km + x)
curve(funMM(x, Vmax = 10, Km = 1), from = 0, to = 12, xlab = "S", ylab = "v", col = 1)
curve(funMM(x, Vmax = 10, Km = 2), from = 0, to = 12, add = TRUE, col = 7)
curve(funMM(x, Vmax = 10, Km = 5), from = 0, to = 12, add = TRUE, col = 4)
abline(v=c(1,2,5), col = c(1,7,4), lty = 2)
abline(h = 5, lty = 3)
legend("bottomright", legend = c("Km = 1", "Km = 2", "Km = 5"), col = c(1, 7, 4), lty = 1)
par(mfrow=c(1,1))

## x: concentração de substrato (conc)
## y: velocidade da reação (rate)
funMM <- function(x, Vmax, Km) {(Vmax * x)/(Km + x)}
curve(funMM(x, Vmax = 1, Km = 2), from = 0, to = 20, ylim = c(0, 1.2))
abline(h = 1, lty=3)
text(0, 1, expression(V[max]), pos = 4)
abline(v = 2, h = 0.5, lty=3)
text(2, 0, expression(K[m]), pos = 3)

## Um arquivo de dados
dadosMM <- data.frame(conc = c(0.1, 0.2, 0.5, 1, 2, 5, 10, 20),
                      rate = c(0.05, 0.09, 0.18, 0.30, 0.45, 0.60, 0.75, 0.85))

## ajuste com nls()
iniMM <- with(dadosMM, list(Vmax = max(rate), Km = 0.25 * max(conc)))
fitMM <- nls(rate ~ Vmax * conc / (Km + conc), 
             data = dadosMM,  start = iniMM)
fitMM
summary(fitMM)

## Comando alternativo para ajuste usando a função selfStart SSmicmen()
nls(rate ~ SSmicmen(conc, Vmax, Km), data = dadosMM)

## visualizando modelo ajustado
plot(rate ~ conc, data = dadosMM, pch = 20,
     xlab = "Concentração de substrato", ylab = "Velocidade da reação")
curve(funMM(x, Vmax = coef(fitMM)[1], Km = coef(fitMM)[2]), 
      from = 0, to = 20, add = TRUE, col = "blue", lwd = 2)

c(lLik = logLik(fitMM), AIC = AIC(fitMM), BIC = BIC(fitMM))

par(mfrow=c(1,2))
plot(profile(fitMM))


## Para estudos de simulação é preciso simular do modelo
## e tentar garantir bons valores iniciais para os ajustes do modelo.

##
## Simulando dados para o modelo de Michaelis Menten
##
par(mar = c(3, 3, 1, 1), mgp = c(2, 0.5, 0))
set.seed(2024)
beta1 = 10; beta2 = 2
#x <- sort(round(runif(15, 0, 12), dig = 1))
x <-  1:15
y <- beta1*x/(beta2 + x) + rnorm(15, 0, 0.5)
df <- data.frame(x = x, y = y); rm(x,y)
with(df, plot(y ~ x, xlim = c(0, 12), ylim = c(0, 11)))
(fit <- nls(y ~ beta1*x/(beta2 + x), start = list(beta1 = 10, beta2 = 2), data = df))

funMM <- function(x, Vmax, Km) Vmax*x/(Km + x)
curve(funMM(x, Vmax = coef(fit)[1], Km = coef(fit)[2]), from = 0, to = 12,
      xlab = "S", ylab = "v", col = 4, add = TRUE, lwd = 2)

## Modelo linearizado (para fornecer ajustes iniciais para o modelo não linear)
## (pode ser feito com lm(),  ou com nls() usando a função selfStart SSmicmen())

fit.ini <- lm(1/y ~ I(1/x), data = df)
beta.ini <- coef(fit.ini)
Vmax.ini <- 1/beta.ini[1]
Km.ini <- beta.ini[2] * Vmax.ini
c(Vmax.ini, Km.ini)
(fit <- nls(y ~ beta1*x/(beta2 + x), start = list(beta1 = Vmax.ini, beta2 = Km.ini), data = df))

(fitSS <- nls(y ~ SSmicmen(x, Vm, K), data = df))

## e ainda corresponde a um GLM !
fitGLM <- glm(y ~ I(1/x), data = df, family = gaussian(link = "inverse"))
fitGLM
coef(fitGLM)
c(1/coef(fitGLM)[1], coef(fitGLM)[2]/coef(fitGLM)[1])


par(mar = c(3, 3, 1, 1), mgp = c(2, 0.5, 0))
with(df, plot(y ~ x, xlim = c(0, 12), ylim = c(0, 11)))
nd <- data.frame(x = seq(0,12, by = 0.25))
nd$pred <- investr::predFit(fit, newdata = nd, interval = "prediction")
with(nd, matlines(x, pred, lty = c(1,2,2), col = 4, lwd = 2))

##
## Outras opções para predição
##

nd <- data.frame(x = seq(0,12, by = 0.25))
#require(propagate)   [Monte Carlo (default) ou aproximação e segunda ordem]
Pred1 <- propagate::predictNLS(fit, newdata = nd, interval = "prediction")
matlines(nd$x, Pred1$summary[,5:6], lty = 1, col = 1)
Pred2 <- propagate::predictNLS(fit, newdata = nd, interval = "prediction", do.sim = FALSE)
names(Pred2)
matlines(nd$x, Pred2$summary[,5:6], lty = 2, col = 2)
## require(investr) [Delta]
Pred3 <- investr::predFit(fit, newdata = nd, interval = "prediction")
head(Pred3)
matlines(nd$x, Pred3[,2:3], lty = 3, col = 4, lwd = 2)
## require(nlraa) : predict_nls() [Monte Carlo] e predict2_nls() [Delta]
Pred4 <- predict_nls(fit, newdata = nd, interval = "prediction")
Pred5 <- predict2_nls(fit, newdata = nd, interval = "prediction")
names(Pred5)
matlines(nd$x, Pred5[,3:4], lty = 1, col = 2, lwd = 2)
## e tem mais glsnls::  nlsr:: minpack.lm:: drc::drm ...e tem mais