###--------------------------------------------------### ###--------------------------------------------------### ### Algumas funções para analise utilizando MCMC ### topicos: ### 1) Modelos com JAGS via R (Regressão linear, Logistica) ### 2) dclone juntamente com JAGS ### 3) MCMCsamp em modelos mistos ### 4) MCMCglmm ### 5) MCMCpack ###--------------------------------------------------### ###--------------------------------------------------### rm(list=ls()) ### JAGS require(runjags) ###--------------------------------------------------### #### Exemplo JAGS (função run-jags) ## Exemplo de regressão linear simples # Simulação de dados X <- 1:50 set.seed(1) Y <- rnorm(length(X), 2*X + 10, 6) plot(X,Y) ###--------------------------------------------------### ### Inferencia por minimos quadrados m0 <- lm(Y~X) summary(m0) ###--------------------------------------------------### ### Modelo y = a*x+b ### Ajustando com o JAGS ### Priori para a ~ dnorm(0,0.00001) ### Priori para b ~ dnorm(0,0.001) ### Priori para sigma ~ exp(1) ### Escrevendo o modelo model <- "model { for(i in 1 : N){ Y[i] ~ dnorm(y.est[i],sigma); y.est[i] <- (a * X[i]) + b; } a ~ dnorm(0,0.00001); b ~ dnorm(0,0.00001); sigma ~ dexp(1); }" ### Data e valor inicial data <- dump.format(list(X=X, Y=Y, N=length(X))) inits0 <- dump.format(list(a=1, b=10, sigma=1)) inits1 <- dump.format(list(a=0, b=0, sigma=1)) # Run the model m1 <- run.jags(model=model, monitor=c("a","b","sigma"), data=data, inits=c(inits0,inits1),n.chains=2, plots = TRUE) ### informações do objeto names(m1) ### Verificando a cadeia plot(m1$mcmc[[1]]) ### Intervalos de credibilidade m1$HPD ### Sumario m1$summary ###--------------------------------------------------### ###--------------------------------------------------### ### Regressão logistica (Livro Introdução Analise Bayesiana) ## pi (n/y) é a proporção de embriões com asas ## ti é a variável “tempo desde a deposição dos ovos” ## função de ligação logit t <- c(5,6,8,8,10,11,16,18) n <- c(34,33,33,35,30,27,33,39) y <- c(6,4,23,18,28,27,33,39) dados <- cbind(y,n-y) m0 <- glm(dados~t,family=binomial(link=logit)) summary(m0) datalist <- dump.format(list(metamorfose=y,total=n,tempo=t)) params <- c("beta0","beta1") inicial <- dump.format(list(beta0=coef(m0)[1], beta1=coef(m0)[2])) ###--------------------------------------------------### modmcmc <- "model{ for(i in 1:length(metamorfose)){ metamorfose[i] ~ dbin(p[i],total[i]) logit(p[i]) <- beta0+beta1*tempo[i] } beta0 ~ dnorm(0,0.001) beta1 ~ dnorm(0,0.001) }" ### rodando o modelo modfit <- run.jags(model=modmcmc,monitor=params,data=datalist, inits=inicial,n.chains=1,burnin=10000,thin=10,sample=19000,check.conv=TRUE) ### Verificando a cadeia plot(modfit$mcmc[[1]]) acf(modfit$mcmc[[1]]) ### Intervalos de credibilidade modfit$HPD ### Sumario modfit$summary ###--------------------------------------------------### ###--------------------------------------------------### ### Modelos com repetição efeito fixo x <- gl(5,4) y <- rnorm(20,10,5) plot(x,y) ###--------------------------------------------------### X <- model.matrix(~x) datalist <- dump.format(list(y=y,X=X,nb=ncol(X))) params <- c("beta1","tau") inicial <- dump.format(list(beta1=rep(0,nlevels(x)),tau=1)) ###--------------------------------------------------### modmcmc <- "model{ for(i in 1:length(y)){ y[i] ~ dnorm(yb[i],tau) yb[i] <- inprod(X[i,], beta1) } tau ~ dunif(0,1000); for(j in 1:nb){ beta1[j] ~ dnorm(0,0.001) } }" ### rodando o modelo modfit <- run.jags(model=modmcmc,monitor=params,data=datalist, inits=inicial,check.conv=TRUE) ### Verificando a cadeia plot(modfit$mcmc[[1]]) ### Intervalos de credibilidade modfit$HPD ### Sumario modfit$summary #### Minimos quadrados summary(lm(y~x)) ###--------------------------------------------------### ###--------------------------------------------------### ### JAGS com Dclone require(dclone) ## simple regression example from the JAGS manual jfun <- function() { for (i in 1:N) { Y[i] ~ dnorm(mu[i], tau) mu[i] <- alpha + beta * (x[i] - x.bar) } x.bar <- mean(x[]) alpha ~ dnorm(10.0, 1.0E-4) beta ~ dnorm(10.0, 1.0E-4) sigma <- 1.0/sqrt(tau) tau ~ dgamma(1, 1) } ## data generation set.seed(1234) N <- 100 alpha <- 1 beta <- -1 sigma <- 0.5 x <- runif(N) linpred <- model.matrix(~x) %*% c(alpha, beta) Y <- rnorm(N, mean = linpred, sd = sigma) plot(x,Y) ## list of data for the model jdata <- list(N = N, Y = Y, x = x) ## what to monitor jpara <- c("alpha", "beta", "sigma") ## fit the model with JAGS regmod <- jags.fit(jdata, jpara, jfun, n.chains = 3) ## model summary summary(regmod) ## data cloning dcdata <- dclone(jdata, 5, multiply = "N") dcmod <- jags.fit(dcdata, jpara, jfun, n.chains = 3) summary(dcmod) #### opções de computação paralela ?jags.parfit ?dc.fit m1 <- dc.fit(dcdata, jpara, jfun,n.clones=c(1,5,10),multiply = "N", n.chains = 3,flavour = "jags") plot(m1) summary(m1) m2 <- dctable(m1) plot(m2) plot(m2,type="log") ###--------------------------------------------------### ###--------------------------------------------------### ### Modelo misto - Amostrando a posteriorir via MCMC require(lme4) ### Resposta normal m0 <- lmer(Reaction ~ Days + (1|Subject) + (0+Days|Subject), sleepstudy) summary(m0) sampm0 <- mcmcsamp(m0, n = 1000) HPDinterval(sampm0) xyplot(sampm0)#ST1 sigma I / sigma e -- ST2 sigma day / sigma e qqmath(sampm0) densityplot(sampm0) ###--------------------------------------------------### ###--------------------------------------------------### browseURL("http://cran-r.c3sl.ufpr.br/web/packages/MCMCglmm/vignettes/CourseNotes.pdf") browseURL("http://cran-r.c3sl.ufpr.br/web/packages/MCMCglmm/vignettes/Overview.pdf") browseURL("http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.160.5098&rep=rep1&type=pdf") require(MCMCglmm) ?MCMCglmm data(PlodiaPO) ### Modelo normal model1<-MCMCglmm(PO~1, random=~FSfamily, data=PlodiaPO, verbose=TRUE) summary(model1) plot.MCMCglmm(model1) ### MOdelo binomial data(PlodiaR) model2 <- MCMCglmm(cbind(Pupated, Infected) ~ 1,random=~FSfamily, family = "multinomial2", data = PlodiaR, verbose = FALSE) plot.MCMCglmm(model2) ###--------------------------------------------------### ###--------------------------------------------------### require(MCMCpack) ### Modelo Poisson counts <- c(18,17,15,20,10,20,25,13,12) outcome <- gl(3,1,9) treatment <- gl(3,3) ?MCMCpoisson posterior <- MCMCpoisson(counts ~ outcome + treatment) plot(posterior) summary(posterior) ###--------------------------------------------------### ###--------------------------------------------------###