[R-br] Como escrever o modelode Michaelis-Menten na linguagem BUGS - R2WinBUGS

[Augusto Ribas]

Opa.
Alguem usa o winbugs com o R?
Eu estava tentando ajustar um Modelo de michelis-menten e não consegui.
Se alguem puder dar uma olhadinha e ver onde estou errado.
Segue um exemplo do que estou fazendo:

library(R2WinBUGS)
#################################################################
#Relação Riqueza~estudos
#################################################################


#Modelo do que eu acho que acontece
#riqueza~(c*estudos)/(d+estudos)

#Fazendo um levantamento do numero de especies de parasitas que
parasitam algumas especies de sapos, eu supuz que qd uma especie tinha
poucos estudos
#tava sempre bem subestimado a riqueza de parasitas associadas aquela especie.
#dai como a unica informação que fico é com o número de estudos,
pensei que deveria ter uma relação onde qt mais estudos, mais especies
eu devo encontrar
#até um limite, o limite é o parametro c nesse caso, e d representa
como os estudos vão contribuindo pra eu chegar ao valor limite de
especies.
#existem alguns problemas acredito, como especies podem ser mais
faceis ou mais dificies de estudar, e o limite pra cada especie
depende de mais coisas e tipos de estudos.
#mas ficando so nesse caso, eu fiz esse exemplo

set.seed(13)
#gerando dados exemplo
estudos<- rpois(48,5)
#parametros originais
#então o limite de especies de parasitos seria 10
c<-10
#e eu preciso de mais ou menos 5 estudos pra ver legal a riqueza de especies,
d<-2.5
#resposta
riqueza<-round(((c*estudos)/(d+estudos))+rnorm(length(estudos),0,1),0)

plot(riqueza~estudos,xlim=c(0,15),ylim=c(0,15),xlab="Riqueza de
espécies",ylab="Número de estudos")
curve(c*x/(d+x),0,15,add=T,col="black")

#ajustando com modelo não linear
modelo.estudos<-nls(riqueza~((a*estudos)/(b+estudos)),start=list(a=8,b=1.5))
summary(modelo.estudos)

points(predict(modelo.estudos)[order(predict(modelo.estudos))]~estudos[order(estudos)],type="l",col="red",lty=2)

#eu li que é possivel lineariza essa função e ajusta com glm da seguinte forma:
#ajustando com glm
modelo.menten<- glm(riqueza~c(1/estudos),family=gaussian(link = "inverse"))
summary(modelo.menten)
a1<-1/modelo.menten$coef[1]
b1<-a1*modelo.menten$coef[2]
a1
b1
curve(a1*x/(b1+x),0,15,add=T,col="blue",lty=2)

### Ajuste com winbugs
# Definir modelo
sink("michaelis.menten.test.txt")
cat("
model {

# Priors
alpha ~ dnorm(0,0.001)
beta ~ dnorm(0,0.001)
sigma ~ dunif(0, 100)

# Likelihood
 for (i in 1:n) {
    riqueza[i] ~ dnorm(mu[i], tau)
    mu[i] <- alpha + beta * estudos[i]
}
# Derived quantities
 tau <- 1/ (sigma * sigma)
}
",fill=TRUE)
sink()

# Juntar os dados - Note que eu não entendi como fazer a função de
ligação dentro do modelo winbugs
#dai eu registrei os dados como 1/riqueza e 1/estudos aqui, acho que
isso que o pepino deve ta nisso
win.data <- list(riqueza=c(1/riqueza),estudos=c(1/estudos),n=length(estudos))

# Valores iniciais
inits <- function(){ list(alpha=rlnorm(1), beta=rlnorm(1),sigma = rlnorm(1))}

# Paramstros a estimar
params <- c("alpha","beta","tau")

# MCMC settings
nc <- 3
ni <- 3000
nb <- 1000
nt <- 2

# Começa o winbugs
out <- bugs(data=win.data, inits=inits, parameters.to.save=params,
model.file="michaelis.menten.test.txt", n.thin=nt, n.chains=nc, n.burnin=nb,
n.iter=ni, debug = TRUE)


###Olhar o modelo
print(out, dig = 3)

#retornoar os parametros como no glm acima
a2<-1/out$summary[1,1]
b2<-a*out$summary[2,1]
a2
b2

curve(a2*x/(b2+x),0,15,add=T,col="green",lty=2)

legend("topright",legend=c("Original","Ajustado com nls","Ajustado com
glm","Ajustado com WinBUGS"),
lty=c(1,2,2,2),col=c("black","red","blue","green"))


#Minhas duvidas são:

#01. como montar o modelo na linguagem bugs pra fazer esse tipo de analise?
#procurando pela internet eu não achei nenhum exemplo pra esse caso especifico.

#02. como eu colocaria dentro do modelo esse função de ligação inversa.
#eu tentei usar a função pow(mu[i],-1) mas não rolou,

#03. e o que estou fazendo no modelo winBUGS tem muito problemas?
entrar com os dados daquele jeito, ao invez de especificar no modelo?
#não consigo imaginar a implicação de não definir no modelo a função de ligação.

Bom se alguem tiver esperiencia com winbugs e puder dar uma olhadinha
no que estou fazendo errado
Ou alguma dica se o que estou fazendo tem algum problema critico

-- 
Grato
Augusto C. A. Ribas

Site Pessoal: http://augustoribas.heliohost.org
Lattes: http://lattes.cnpq.br/7355685961127056