########################################################################################## ### Example 1 - Life quality index of industry workers ################################### ########################################################################################## rm(list=ls()) ## Loading extra function source("functions.R") ## Loading extra packages require(lattice) require(betareg) require(bbmle) require(mvtnorm) require(reshape) require(dclone) require(rjags) require(R2WinBUGS) ## Loading the data set dados <- read.table("datasetIQVT.txt",header=TRUE)[,c(1,5,7,8)] levels(dados$ESTADO) <- c("AM","CE","DF","MT","MS","PB","PR","RO","RR") ## Log-icome dados$RENDA <- log(dados$RENDA) - mean(log(dados$RENDA)) ## Curve of descriptive graphical mod.graf <- betareg(GLOBAL ~ RENDA,data=dados) RENDA <- data.frame(RENDA = seq(-0.5,1.1,l=100)) preditos = predict(mod.graf,newdata=RENDA) ## Centering income pch.id <- unique(as.numeric(dados$ESTADO)) mat <- matrix(c(1,2,1,3),2,2) layout(mat) par(mar=c(2.8, 2.8, 1.7, 1),mgp = c(1.8, 0.7, 0)) plot(GLOBAL ~ RENDA,pch=pch.id,col=pch.id,data=dados,xlim=c(-0.5,1.1),xlab="Log-Renda") lines(preditos ~ as.numeric(as.matrix(RENDA))) legend("bottomright",legend=levels(dados$ESTADO),col=pch.id,pch=pch.id) par(mar=c(2.8, 2.8, 1.7, 1),mgp = c(1.8, 0.7, 0)) boxplot(GLOBAL ~ PORTE,varwidth=TRUE,data=dados) par(mar=c(2.8, 2.8, 1.7, 1),mgp = c(1.8, 0.7, 0)) boxplot(GLOBAL ~ ESTADO,varwidth=TRUE,data=dados) ## Model beta with intercept and slope random vero.slope <- function(uv,beta.fixo, prec.pars, X, Z, Y,log=TRUE){ phi <- exp(prec.pars[1]) tau.U <- exp(prec.pars[2])^2 tau.V <- exp(prec.pars[3])^2 rho <- (exp(2*prec.pars[4]) - 1)/(exp(2*prec.pars[4])+1) cov.rho <- rho*(sqrt(tau.U)*sqrt(tau.V)) if(class(dim(uv)) == "NULL"){uv <- matrix(uv,1,2)} ll = apply(uv,1,function(uvi){ preditor <- X%*%beta.fixo + Z%*%as.numeric(uvi) mu <- inv.logit(preditor) sigma <- matrix(c(tau.U,cov.rho,cov.rho,tau.V),2,2) sum(dbeta(Y, mu*phi, (1-mu)*phi, log=TRUE)) + dmvnorm(uvi, c(0,0), sigma = sigma , log=TRUE)}) if(log == FALSE){ll <- exp(ll)} return(ll)} ########################################################################################## ## Fitting 6 models #### ################################################################# ## 1 - Without covariates and random effects ############################################# ## 2 - Covariate PORTE ################################################################## ## 3 - Covariate PORTE + RENDA ########################################################### ## 4 - Random intercept ############################################################### ## 5 - Random intercept and slope with correlation ####################################### ########################################################################################## dados$NADA1 <- 0 dados$NADA2 <- 0 names(dados) <- c("y","RENDA","PORTE","ID","NADA1","NADA2") ########################################################################################## ## Warning - Note that here all models are reparametrized, so you need take into account## ## to compare the standard error with paper !!!!!!!!!!!!! ############################### ########################################################################################## ########################################################################################## ## Model 1 - beta_0 ###################################################################### ########################################################################################## mod1 <- function(b1,phi,integral, pontos, otimizador, n.dim, dados){ ll = verossimilhanca(modelo = vero.slope, formu.X="~1",formu.Z="~NADA1 + NADA2 - 1", beta.fixo = b1, prec.pars=c(phi,-10,-10,0), integral=integral, pontos=pontos, otimizador=otimizador,n.dim=n.dim,dados=dados) print(round(c(b1,ll),2)) return(-ll)} ajuste.1 = mle2(mod1,start=list(b1=0,phi=log(70)), data=list(integral="LAPLACE",pontos=1,otimizador="BFGS", n.dim=2,dados=dados)) ########################################################################################## ## Model 2 - beta_0 + PORTE ############################################################## ########################################################################################## mod2 <- function(b1,b2,b3,phi,integral, pontos, otimizador, n.dim, dados){ ll = verossimilhanca(modelo = vero.slope, formu.X="~PORTE",formu.Z="~NADA1 + NADA2 - 1", beta.fixo = c(b1,b2,b3), prec.pars=c(phi,-10,-10,0), integral=integral, pontos=pontos, otimizador=otimizador,n.dim=n.dim,dados=dados) print(round(c(b1,b2,b3,ll),2)) return(-ll)} ajuste.2 = mle2(mod2,start=list(b1=0,b2=0,b3=0,phi=log(70)), data=list(integral="LAPLACE",pontos=1,otimizador="BFGS", n.dim=2,dados=dados)) ########################################################################################## ## Model 3 - PORTE + RAMO ############################################################## ########################################################################################## mod3 <- function(b1,b2,b3,b4,phi, integral, pontos, otimizador, n.dim, dados){ ll = verossimilhanca(modelo = vero.slope, formu.X="~ PORTE + RENDA",formu.Z="~NADA1 + NADA2 - 1", beta.fixo = c(b1,b2,b3,b4), prec.pars=c(phi,-10,-10,0), integral=integral, pontos=pontos, otimizador=otimizador,n.dim=n.dim,dados=dados) print(round(c(b1,b2,b3,b4,phi,ll),2)) return(-ll)} ajuste.3 = mle2(mod3,start=list(b1=0,b2=0,b3=0,b4=0, phi=log(70)), data=list(integral="LAPLACE",pontos=1,otimizador="BFGS", n.dim=2,dados=dados)) ########################################################################################## ## Model 4 - Random intercept ########################################################### ########################################################################################## mod4 <- function(b1,b2,b3,b4,phi,tau.U, integral, pontos, otimizador, n.dim, dados){ ll = verossimilhanca(modelo = vero.slope, formu.X="~ PORTE + RENDA",formu.Z="~NADA1", beta.fixo = c(b1,b2,b3,b4), prec.pars=c(phi,tau.U,-10,0), integral=integral, pontos=pontos, otimizador=otimizador,n.dim=n.dim,dados=dados) print(round(c(b1,b2,b3,b4,phi,tau.U,ll),2)) return(-ll)} ajuste.4 = mle2(mod4,start=list(b1=0,b2=0,b3=0,b4=0,phi=4.28,tau.U=0), method="BFGS", data=list(integral="LAPLACE",pontos=1,otimizador="BFGS", n.dim=2,dados=dados)) ########################################################################################## ## Model 5 - Intercept and slope com correlation########################################## ########################################################################################## mod5 <- function(b1,b2,b3,b4,phi,tau.U,tau.V,rho, integral, pontos, otimizador, n.dim, dados){ ll = verossimilhanca(modelo = vero.slope, formu.X="~PORTE + RENDA",formu.Z="~RENDA", beta.fixo = c(b1,b2,b3,b4), prec.pars=c(phi,tau.U,tau.V,rho), integral=integral, pontos=pontos, otimizador=otimizador,n.dim=n.dim,dados=dados) print(round(c(b1,b2,b3,b4,phi,tau.U,tau.V,rho,ll),2)) return(-ll)} ajuste.5 = mle2(mod5,start=list(b1=0,b2=0,b3=0,b4=0, phi=4.28,tau.U=0,tau.V=0,rho=0), data=list(integral="LAPLACE",pontos=1,otimizador="BFGS", n.dim=2,dados=dados)) ## Profile likelihood (Warning !! Time consuming) perf <- profile(ajuste.4) ########################################################################################## ## Comparing models ###################################################################### ########################################################################################## LogVero = cbind(logLik(ajuste.1),logLik(ajuste.2),logLik(ajuste.3),logLik(ajuste.4), logLik(ajuste.5)) Akaike = cbind(AIC(ajuste.1),AIC(ajuste.2),AIC(ajuste.3),AIC(ajuste.4), AIC(ajuste.5)) coef.MOD1 <- c(coef(ajuste.1)[1], NA, NA, NA, exp(coef(ajuste.1)[2]), NA,NA,NA) coef.MOD2 <- c(coef(ajuste.2)[1:3],NA, exp(coef(ajuste.2)[4]), NA,NA,NA) coef.MOD3 <- c(coef(ajuste.3)[1:4],exp(coef(ajuste.3)[5]), 1/exp(coef(ajuste.3)[6])^2,NA,NA) coef.MOD4 <- c(coef(ajuste.4)[1:4],exp(coef(ajuste.4)[5]), 1/exp(coef(ajuste.4)[6])^2,1/exp(coef(ajuste.4)[7])^2,NA) coef.MOD5 <- c(coef(ajuste.5)[1:4],exp(coef(ajuste.5)[5]), 1/exp(coef(ajuste.5)[6])^2,1/exp(coef(ajuste.5)[7])^2, (exp(2*coef(ajuste.5)[8])- 1)/(exp(2*coef(ajuste.5)[8]) + 1)) tab <- round(rbind(cbind(coef.MOD1,coef.MOD2,coef.MOD3,coef.MOD4, coef.MOD5),LogVero,Akaike),2) ########################################################################################## ## Fitting via dclone #################################################################### ########################################################################################## ## data-set to use inside jags 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))) } ########################################################################################## ## Random intercept model ################################################################ ########################################################################################## ## Fitting in JAGS 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) bayes <- jags.fit(dat.jags,c("beta1","beta2","beta3","beta4","phi","tau"),beta.estruturado,n.adapt = 2000, n.update = 500, thin = 5, n.iter= 5000, n.chains=3) plot(bayes) summary(bayes) ########################################################################################## ### Using clone data with dclone package ################################################# ########################################################################################## k <- c(1,5,10,20,30,40,50) clone.mod1 <- dc.fit(dat.jags,c("beta1","beta2","beta3","beta4","phi","tau"),beta.estruturado,n.clones=k,n.adapt=1000, n.update = 500, thin = 5, n.iter= 6500, multiply="n.bloco",unchanged = "n.rep") ## Summary of fitted model plot(clone.mod1) summary(clone.mod1) esti <- dctable(clone.mod1) plot(esti) plot(esti,type="log") ########################################################################################## ## Random intercept and slope ############################################################ ########################################################################################## k <- c(1,5,10,20,30,40,50) clone.data.slope <- dc.fit(dat.jags,c("beta1","beta2","beta3","beta4","phi","tau.U","tau.V"),beta.estruturado.slope,n.clones=k,n.adapt=1000, n.update = 500, thin = 5, n.iter= 5000, multiply="n.bloco",unchanged = "n.rep") ## Summary of fitted model plot(clone.data.slope) summary(clone.data.slope) esti <- dctable(clone.data.slope) plot(esti) plot(esti,type="log") ##########################################################################################