########################################################################################## ## Laplace integration ################################################################### ########################################################################################## laplace <- function(funcao, otimizador,n.dim, ...){ integral <- -sqrt(.Machine$double.xmax) inicial <- rep(0,n.dim) temp <- try(optim(inicial,funcao,...,method=otimizador,hessian=TRUE,control=list(fnscale=-1))) if(class(temp) != "try-error"){ integral <- exp(temp$value)* (exp((n.dim/2)*log(2*pi) - 0.5*determinant(-temp$hessian)$modulus))} return(integral) } ########################################################################################## ## Gauss Hermite n-dimensional ########################################################### ########################################################################################## gauss.hermite <- function(funcao, n.dim, n.pontos, ... ){ normaliza <- function(x){exp(-t(as.numeric(x))%*%as.numeric(x))} pontos <- gauss.quad(n.pontos,kind="hermite") nodes <- matrix(rep(pontos$nodes,n.dim),ncol=n.dim) pesos <- matrix(rep(pontos$weights,n.dim),ncol=n.dim) lista.nodes <- list() lista.pesos <- list() for(i in 1:ncol(nodes)){ lista.nodes[[i]] <- nodes[,i] lista.pesos[[i]] <- pesos[,i]} nodes = do.call(expand.grid,lista.nodes) pesos = do.call(expand.grid,lista.pesos) pesos.grid = apply(pesos,1,prod) norma = apply(nodes,1,normaliza) integral <- sum(pesos.grid*(funcao(nodes,...)/norma)) return(integral) } ########################################################################################## ## Quasi Monte Carlo n-dimensional ####################################################### ########################################################################################## quase.mc <- function(funcao, n.dim, n.pontos, tipo, ...){ if( tipo == "Halton"){ pontos <- rnorm.halton(n.pontos,n.dim)} if( tipo == "Sobol"){ pontos <- rnorm.sobol(n.pontos,n.dim)} norma <- apply(pontos,1,dmvnorm) integral <- mean(funcao(pontos,...)/norma) return(integral) } ########################################################################################## ## Adaptative Gauss-Hermite ############################################################## ########################################################################################## adaptative.gauss.hermite <- function(funcao, n.dim, n.pontos, otimizador, ... ){ normaliza <- function(x){exp(-t(as.numeric(x))%*%as.numeric(x))} pontos <- gauss.quad(n.pontos,kind="hermite") integral <- -sqrt(.Machine$double.xmax) inicial <- rep(0,n.dim) temp <- try(optim(inicial,funcao,...,method=otimizador,hessian=TRUE,control=list(fnscale=-1))) #temp <- try(optim(inicial,funcao,beta.fixo=beta.fixo, prec.pars=prec.pars, X = X, Z=Z, Y = dados.id[[1]]$y, # method=otimizador,hessian=TRUE,control=list(fnscale=-1))) if(class(temp) != "try-error"){ z.chapeu <- temp$par sd.chapeu <- sqrt(diag(solve(-temp$hessian))) mat.nodes <- matrix(NA, ncol=n.dim,nrow=n.pontos) mat.pesos <- matrix(NA,ncol=n.dim,nrow=n.pontos) for(i in 1:length(z.chapeu)){ mat.nodes[,i] <- z.chapeu[i] + sd.chapeu[i]*pontos$nodes mat.pesos[,i] <- sd.chapeu[i]* (exp(-mat.nodes[,i]^2)/exp(-pontos$nodes^2))*pontos$weights } lista.nodes <- list() lista.pesos <- list() for(i in 1:ncol(mat.nodes)){ lista.nodes[[i]] <- mat.nodes[,i] lista.pesos[[i]] <- mat.pesos[,i]} nodes = do.call(expand.grid,lista.nodes) pesos = do.call(expand.grid,lista.pesos) pesos.grid = apply(pesos,1,prod) norma = apply(nodes,1,normaliza) integral <- sum(pesos.grid*(exp(funcao(nodes,...))/norma)) } return(integral) } ########################################################################################## ## Adaptative Quase Monte Carlo ########################################################## ########################################################################################## adaptative.monte.carlo <- function(funcao, n.pontos, n.dim, tipo, otimizador, ... ){ if(tipo == "Halton") {pontos <- rnorm.halton(n.pontos, n.dim)} if(tipo == "Sobol") {pontos <- rnorm.sobol(n.pontos, n.dim)} integral <- -sqrt(.Machine$double.xmax) inicial <- rep(0,n.dim) temp <- try(optim(inicial,funcao,...,method=otimizador,hessian=TRUE,control=list(fnscale=-1))) #temp <- try(optim(inicial,funcao,beta.fixo=beta.fixo, prec.pars=prec.pars, X.fixo = X.mod, X.ale = X.ale , # Y = dados.id[[1]]$y, log=TRUE,method=otimizador,hessian=TRUE,control=list(fnscale=-1))) if(class(temp) != "try-error"){ z.chapeu <- temp$par H <- solve(-temp$hessian) sd.chapeu <- sqrt(diag(H)) mat.nodes <- matrix(NA, ncol=n.dim,nrow=n.pontos) for(i in 1:length(z.chapeu)){ mat.nodes[,i] <- z.chapeu[i] + sd.chapeu[i]*pontos[,i] } norma <- dmvnorm(mat.nodes,mean=z.chapeu,sigma=H,log=TRUE) integral = mean(exp(funcao(mat.nodes,...) - norma)) } return(integral) } ########################################################################################## ## Generic function to evaluate the marginal likelihood ################################## ########################################################################################## inv.logit <- function(x){exp(x)/(1+exp(x))} verossimilhanca <- function(modelo, formu.X, formu.Z, beta.fixo, prec.pars, integral, pontos, otimizador, n.dim, dados){ dados.id <- split(dados, dados$ID) ll <- c() for(i in 1:length(dados.id)){ X <- model.matrix(as.formula(formu.X),data=dados.id[[i]]) Z <- model.matrix(as.formula(formu.Z),data=dados.id[[i]]) if(integral == "LAPLACE"){ ll[i] <- laplace(modelo,otimizador=otimizador,n.dim=n.dim, X=X, Z=Z , Y=dados.id[[i]]$y, beta.fixo=beta.fixo, prec.pars=prec.pars,log=TRUE)} if(integral == "GH"){ ll[i] <- gauss.hermite(modelo,n.pontos=pontos, n.dim=n.dim, X=X, Z=Z, Y=dados.id[[i]]$y, beta.fixo=beta.fixo, prec.pars=prec.pars,log=FALSE)} if(integral == "QMH"){ ll[i] <- quase.mc(modelo,n.pontos=pontos, n.dim=n.dim, tipo= "Halton", X=X, Z=Z, Y=dados.id[[i]]$y, beta.fixo=beta.fixo, prec.pars=prec.pars,log=FALSE)} if(integral == "QMS"){ ll[i] <- quase.mc(modelo,n.pontos=pontos, n.dim=n.dim, tipo= "Sobol", X=X, Z=Z, Y=dados.id[[i]]$y, beta.fixo=beta.fixo, prec.pars=prec.pars,log=FALSE)} if(integral == "AGH"){ ll[i] <- adaptative.gauss.hermite(modelo,n.pontos=pontos, n.dim=n.dim, otimizador=otimizador, X=X, Z=Z, Y=dados.id[[i]]$y, beta.fixo=beta.fixo, prec.pars=prec.pars,log=TRUE) } if(integral == "AQMH"){ ll[i] <- adaptative.monte.carlo(modelo,n.pontos=pontos, n.dim=n.dim, otimizador=otimizador, tipo="Halton", X=X, Z=Z, Y=dados.id[[i]]$y, beta.fixo=beta.fixo, prec.pars=prec.pars,log=TRUE)} if(integral == "AQMS"){ ll[i] <- adaptative.monte.carlo(modelo,n.pontos=pontos, n.dim=n.dim, otimizador=otimizador, tipo="Sobol", X=X, Z=Z, Y=dados.id[[i]]$y, beta.fixo=beta.fixo, prec.pars=prec.pars,log=TRUE)} } return(sum(log(ll))) } ########################################################################################## ## Beta mixed models - Random intercept JAGS implementation ############################## ########################################################################################## beta.estruturado <- 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(0.1,0.01) tau ~ dgamma(0.1,0.01) } ########################################################################################## ## Beta mixed models - Random intercept and slope JAGS implementation #################### ########################################################################################## beta.estruturado.slope <- 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 + v[j])*mat.renda[j,i] } b[j] ~ dnorm(0,tau.U) v[j] ~ dnorm(0,tau.V) } 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.U ~ dgamma(0.1,0.01) tau.V ~ dgamma(0.1,0.01) } ########################################################################################## ## Beta mixed models - IQA example ####################################################### ########################################################################################## beta.estruturado.IQA <- 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]) <- X.int[j,i]*beta0 + X.reser[j,i]*beta1 + X.jusa[j,i]*beta2 + X.trim2[j,i]*beta3 + X.trim3[j,i]*beta4 + X.trim4[j,i]*beta5 + Z.int[j,i]*b0[j] + Z.trim1[j,i]*b1[j] + Z.trim2[j,i]*b2[j] + Z.trim3[j,i]*b3[j] + Z.trim4[j,i]*b4[j] } b0[j] ~ dnorm(0,tau.U) b1[j] ~ dnorm(0,tau.UV) b2[j] ~ dnorm(0,tau.UV) b3[j] ~ dnorm(0,tau.UV) b4[j] ~ dnorm(0,tau.UV) } beta0 ~ dnorm(0,0.001) beta1 ~ dnorm(0,0.001) beta2 ~ dnorm(0,0.001) beta3 ~ dnorm(0,0.001) beta4 ~ dnorm(0,0.001) beta5 ~ dnorm(0,0.001) phi ~ dgamma(1,0.01) tau.U ~ dgamma(1,0.01) tau.UV ~ dgamma(1,0.01) } ########################################################################################## ## Beta mixed models - IQA example model 2 ############################################### ########################################################################################## beta.estruturado.IQA2 <- 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]) <- X.int[j,i]*beta0 + X.reser[j,i]*beta1 + X.jusa[j,i]*beta2 + X.trim2[j,i]*beta3 + X.trim3[j,i]*beta4 + X.trim4[j,i]*beta5 + Z.trim1[j,i]*b1[j] + Z.trim2[j,i]*b2[j] + Z.trim3[j,i]*b3[j] + Z.trim4[j,i]*b4[j] } b1[j] ~ dnorm(0,tau.UV) b2[j] ~ dnorm(0,tau.UV) b3[j] ~ dnorm(0,tau.UV) b4[j] ~ dnorm(0,tau.UV) } beta0 ~ dnorm(0,0.001) beta1 ~ dnorm(0,0.001) beta2 ~ dnorm(0,0.001) beta3 ~ dnorm(0,0.001) beta4 ~ dnorm(0,0.001) beta5 ~ dnorm(0,0.001) phi ~ dgamma(1,0.01) tau.UV ~ dgamma(1,0.01) }