########################################################################################## ## Modelos escritos em JAGS ############################################################## ########################################################################################## ## Montando as matrizes para passar pro JAGS ## Montando a matriz do modelo mat.Y <- matrix(NA,ncol=9,nrow=80) mat.intercepto <- matrix(NA,ncol=9,nrow=80) mat.renda <- matrix(NA,ncol=9,nrow=80) mat.Media <- matrix(NA,ncol=9,nrow=80) mat.Pequena <- matrix(NA,ncol=9,nrow=80) dados.id <- split(dados,dados$ID) for(i in 1:9){ temp <- dados.id[[i]]$y modelo <- model.matrix(~ dados.id[[i]]$PORTE + dados.id[[i]]$RENDA) completa = rep(NA,80-length(temp)) mat.Y[,i] <- c(temp,completa) mat.intercepto[,i] <- c(modelo[,1],rep(mean(modelo[,1]),length(completa))) mat.Media[,i] <- c(modelo[,2],rep(0,length(completa))) mat.Pequena[,i] <- c(modelo[,3],rep(0,length(completa))) mat.renda[,i] <- c(modelo[,4],rep(0,length(completa))) } dat.jags <- list(Y =t(mat.Y), mat.intercepto = t(mat.intercepto), mat.renda = t(mat.renda), mat.Media = t(mat.Media), mat.Pequena = t(mat.Pequena), n.bloco=9,n.rep=80) ## Modelo só com intercepto modelo1 <- function(){ for(j in 1:n.bloco){ for(i in 1:n.rep){ Y[j,i] ~ dbeta(mu[j,i]*phi,(1-mu[j,i])*phi) logit(mu[j,i]) <- mat.intercepto[j,i]*beta1 } } beta1 ~ dnorm(0,0.001) phi ~ dgamma(0.1,0.01) } ## Modelo com efeito de porte modelo2 <- function(){ for(j in 1:n.bloco){ for(i in 1:n.rep){ Y[j,i] ~ dbeta(mu[j,i]*phi,(1-mu[j,i])*phi) logit(mu[j,i]) <- mat.intercepto[j,i]*beta1 + mat.Media[j,i]*beta2 + mat.Pequena[j,i]*beta3 } } beta1 ~ dnorm(0,0.001) beta2 ~ dnorm(0,0.001) beta3 ~ dnorm(0,0.001) phi ~ dgamma(0.1,0.01) } ## Modelo com efeito de Porte e Renda modelo3 <- function(){ for(j in 1:n.bloco){ for(i in 1:n.rep){ Y[j,i] ~ dbeta(mu[j,i]*phi,(1-mu[j,i])*phi) logit(mu[j,i]) <- mat.intercepto[j,i]*beta1 + mat.Media[j,i]*beta2 + mat.Pequena[j,i]*beta3 + mat.renda[j,i]*beta4 } } beta1 ~ dnorm(0,0.001) beta2 ~ dnorm(0,0.001) beta3 ~ dnorm(0,0.001) beta4 ~ dnorm(0,0.001) phi ~ dgamma(0.1,0.01) } ## Modelo com intercepto aleatório modelo4 <- function(){ for(j in 1:n.bloco){ for(i in 1:n.rep){ Y[j,i] ~ dbeta(mu[j,i]*phi,(1-mu[j,i])*phi) logit(mu[j,i]) <- mat.intercepto[j,i]*beta1 + mat.Media[j,i]*beta2 + mat.Pequena[j,i]*beta3 + mat.renda[j,i]*beta4 + b[j] } b[j] ~ dnorm(0,tau) } beta1 ~ dnorm(0,0.001) beta2 ~ dnorm(0,0.001) beta3 ~ dnorm(0,0.001) beta4 ~ dnorm(0,0.001) phi ~ dgamma(1,0.01) tau ~ dgamma(0.5, 0.001487) } ## Modelo com intercepto aleatório + slope aleatório ## Modelo só com intercepto modelo1 <- function(){ for(j in 1:n.bloco){ for(i in 1:n.rep){ Y[j,i] ~ dbeta(mu[j,i]*phi,(1-mu[j,i])*phi) logit(mu[j,i]) <- mat.intercepto[j,i]*beta1 } } beta1 ~ dnorm(0,0.001) phi ~ dgamma(0.1,0.01) } ## Modelo com efeito de porte modelo2 <- function(){ for(j in 1:n.bloco){ for(i in 1:n.rep){ Y[j,i] ~ dbeta(mu[j,i]*phi,(1-mu[j,i])*phi) logit(mu[j,i]) <- mat.intercepto[j,i]*beta1 + mat.Media[j,i]*beta2 + mat.Pequena[j,i]*beta3 } } beta1 ~ dnorm(0,0.001) beta2 ~ dnorm(0,0.001) beta3 ~ dnorm(0,0.001) phi ~ dgamma(0.1,0.01) } ## Modelo com efeito de Porte e Renda modelo3 <- function(){ for(j in 1:n.bloco){ for(i in 1:n.rep){ Y[j,i] ~ dbeta(mu[j,i]*phi,(1-mu[j,i])*phi) logit(mu[j,i]) <- mat.intercepto[j,i]*beta1 + mat.Media[j,i]*beta2 + mat.Pequena[j,i]*beta3 + mat.renda[j,i]*beta4 } } beta1 ~ dnorm(0,0.001) beta2 ~ dnorm(0,0.001) beta3 ~ dnorm(0,0.001) beta4 ~ dnorm(0,0.001) phi ~ dgamma(0.1,0.01) } ## Modelo com intercepto aleatório modelo4 <- function(){ for(j in 1:n.bloco){ for(i in 1:n.rep){ Y[j,i] ~ dbeta(mu[j,i]*phi,(1-mu[j,i])*phi) logit(mu[j,i]) <- mat.intercepto[j,i]*beta1 + mat.Media[j,i]*beta2 + mat.Pequena[j,i]*beta3 + mat.renda[j,i]*beta4 + b[j] } b[j] ~ dnorm(0,tau) } beta1 ~ dnorm(0,0.001) beta2 ~ dnorm(0,0.001) beta3 ~ dnorm(0,0.001) beta4 ~ dnorm(0,0.001) phi ~ dgamma(1,0.1) tau ~ dgamma(1,0.00005) } ## Modelo com intercepto aleatório + slope aleatório modelo5 <- function(){ for(j in 1:n.bloco){ for(i in 1:n.rep){ Y[j,i] ~ dbeta(mu[j,i]*phi,(1-mu[j,i])*phi) logit(mu[j,i]) <- (mat.intercepto[j,i]*beta1 + b[j]) + mat.Media[j,i]*beta2 + mat.Pequena[j,i]*beta3 + (beta4 + d[j])*mat.renda[j,i] } b[j] ~ dnorm(0,tau) d[j] ~dnorm(0, tau2) } beta1 ~ dnorm(0,0.001) beta2 ~ dnorm(0,0.001) beta3 ~ dnorm(0,0.001) beta4 ~ dnorm(0,0.001) phi ~ dgamma(0.1,0.01) tau ~ dgamma(0.1,0.01) tau2 ~ dgamma(0.1,0.01) }