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