########################################################################################## ## Example 2 - Water quality index ####################################################### ########################################################################################## ## Loading extra function source("functions.R") ## Loading extra packages require(betareg) require(mvtnorm) require(bbmle) require(statmod) require(fOptions) require(dclone) require(rjags) ## Lendo os dados dados <- read.table("datasetIQA.txt",header=TRUE) head(dados) ## 15 power plants table(dados$UHE) ## 3 points table(dados$LOCAL) ## 4 times table(dados$TRIM) ## Check classes sapply(dados,class) dados$TRIM <- as.factor(dados$TRIM) dados$LOCAL <- factor(dados$LOCAL, levels=c("MONT","RESER","JUSA")) names(dados) <- c("ID","y","LOCAL","TRIM") ########################################################################################## ## Descriptive analysis mat <- matrix(c(1,3,2,4),ncol=2) layout(mat) par(mar=c(2,2.5,1.5,1.5),mgp=c(1.8,0.8,0)) hist(dados$y,main="A",ylab="Frequency",xlab="IQA") par(mar=c(6.2,2.5,1.5,1.5),mgp=c(1.8,0.8,0)) boxplot(dados$y ~ dados$ID,xlab="",ylab="IQA",xaxt="n",main="B",notch=TRUE,outline=FALSE) labels <- levels(dados$ID) text(1:length(labels), par("usr")[3] - 0.05, srt = 90, adj = 1, labels = labels, xpd = TRUE, cex = 0.9) par(mar=c(2,2.5,1.5,1.5),mgp=c(1.8,0.8,0)) boxplot(dados$y ~ dados$LOCAL) par(mar=c(2,2.5,1.5,1.5),mgp=c(1.8,0.8,0)) boxplot(dados$y ~ dados$TRIM) ## Histogramas by level histogram(~ y | ID,data=dados) histogram(~ y | LOCAL,data=dados) histogram(~ y | TRIM,data=dados) ########################################################################################## ## Dispersion diagram #################################################################### ########################################################################################## require(lattice) xyplot(y ~ as.numeric(TRIM) | LOCAL, group=ID,data=dados,type="l") ########################################################################################## ## Writing the generic model in a suitable way ########################################### ########################################################################################## modelo <- function(uv,beta.fixo, prec.pars, X, Z, Y,log=TRUE){ phi <- exp(prec.pars[1]) tau.U <- exp(prec.pars[2]) tau.UT <- exp(prec.pars[3]) if(class(dim(uv)) == "NULL"){uv <- matrix(uv,1,length(uv))} ll = apply(uv,1,function(uvi){ preditor <- X%*%beta.fixo + Z%*%as.numeric(uvi) mu <- inv.logit(preditor) sum(dbeta(Y, mu*phi, (1-mu)*phi, log=TRUE)) + dnorm(uvi[1], 0, sd = tau.U , log=TRUE) + dnorm(uvi[2], 0,sd = tau.UT,log=TRUE) + dnorm(uvi[3],0, sd=tau.UT,log=TRUE) + dnorm(uvi[4],0, sd=tau.UT,log=TRUE) + dnorm(uvi[5],0, sd=tau.UT,log=TRUE) }) if(log == FALSE){ll <- exp(ll)} return(ll)} ########################################################################################## ## Rename collumns ans create instrumental collumns ###################################### ########################################################################################## dados$NADA1 <- 0 dados$NADA2 <- 0 dados$NADA3 <- 0 dados$NADA4 <- 0 dados$NADA5 <- 0 dados$UM <- 1 dados$TRIM <- as.factor(dados$TRIM) names(dados) <- c("ID","y","LOCAL","TRIM","NADA1","NADA2","NADA3","NADA4","NADA5","UM") ########################################################################################## ## Model 1 - Intercept ################################################################### ########################################################################################## mod1 <- function(b1,phi,integral, pontos, otimizador, n.dim, dados){ ll = verossimilhanca(modelo = modelo, formu.X="~1",formu.Z="~NADA1 + NADA2 + NADA3 + NADA4 + NADA5 - 1", beta.fixo = b1, prec.pars=c(phi,-100,-100), integral=integral, pontos=pontos, otimizador=otimizador,n.dim=n.dim,dados=dados) print(round(c(b1,phi,ll),2)) return(-ll)} ajuste.1 = mle2(mod1,start=list(b1=0,phi=log(100)), data=list(integral="LAPLACE",pontos=1,otimizador="BFGS", n.dim=5,dados=dados)) summary(ajuste.1) ########################################################################################## ## Model 2 - Intercept + LOCAL ########################################################## ########################################################################################## mod2 <- function(b1,b2,b3, phi,integral, pontos, otimizador, n.dim, dados){ ll = verossimilhanca(modelo = modelo, formu.X="~LOCAL",formu.Z="~NADA1 + NADA2 + NADA3 + NADA4 + NADA5 - 1", beta.fixo = c(b1,b2,b3), prec.pars=c(phi,-100,-100), integral=integral, pontos=pontos, otimizador=otimizador,n.dim=n.dim,dados=dados) print(round(c(b1,b2,b3,phi,ll),2)) return(-ll)} ajuste.2 = mle2(mod2,start=list(b1=0, b2=0,b3=0,phi=log(100)), data=list(integral="LAPLACE",pontos=1,otimizador="BFGS", n.dim=5,dados=dados)) summary(ajuste.2) ########################################################################################## ## Model 3 - Intercept + LOCAL + TRIM ################################################### ########################################################################################## mod3 <- function(b1,b2,b3,b4,b5,b6, phi,integral, pontos, otimizador, n.dim, dados){ ll = verossimilhanca(modelo = modelo, formu.X="~LOCAL + TRIM",formu.Z="~NADA1 + NADA2 + NADA3 + NADA4 + NADA5 - 1", beta.fixo = c(b1,b2,b3,b4,b5,b6), prec.pars=c(phi,-100,-100), integral=integral, pontos=pontos, otimizador=otimizador,n.dim=n.dim,dados=dados) print(round(c(b1,b2,b3,b4,b5,b6,phi,ll),2)) return(-ll)} ajuste.3 = mle2(mod3,start=list(b1=0,b2=0, b3=0,b4=0, b5=0, b6=0, phi=3), data=list(integral="LAPLACE",pontos=1,otimizador="BFGS", n.dim=5,dados=dados)) summary(ajuste.3) logLik(ajuste.3) ########################################################################################## ## Model 4 - Intercept + LOCAL + TRIM + UHE ############################################## ########################################################################################## mod4 <- function(b1, b2,b3,b4,b5,b6, phi, tau.U, integral, pontos, otimizador, n.dim, dados){ ll = verossimilhanca(modelo = modelo, formu.X="~LOCAL + TRIM",formu.Z="~NADA1 + NADA2 + NADA3 + NADA4", beta.fixo = c(b1,b2,b3,b4,b5,b6), prec.pars=c(phi,tau.U,-100), integral=integral, pontos=pontos, otimizador=otimizador,n.dim=n.dim,dados=dados) print(round(c(b1,b2,b3,b4,b5,b6,phi,ll),2)) return(-ll)} vetro <- c(1e-3,1e-3,1e-4,1e-3,1e-4,1e-4,1e-3,1e-3) ajuste.4 = mle2(mod4,start=list(b1=0, b2=0,b3=0,b4=0,b5=0,b6=0,phi=3,tau.U=0), method="BFGS",control=list(ndeps=vetro,lmm=5), data=list(integral="LAPLACE",pontos=1,otimizador="BFGS", n.dim=5,dados=dados)) summary(ajuste.4) logLik(ajuste.4) ########################################################################################## ## Model 5 - Intercept + LOCAL + TRIM + UHE:TRIM ######################################### ########################################################################################## mod5 <- function(b1, b2,b3,b4,b5,b6, phi, tau.UT, integral, pontos, otimizador, n.dim, dados){ ll = verossimilhanca(modelo = modelo, formu.X="~LOCAL + TRIM",formu.Z="~ UM + TRIM - 1", beta.fixo = c(b1,b2,b3,b4,b5,b6), prec.pars=c(phi,-100,tau.UT), integral=integral, pontos=pontos, otimizador=otimizador,n.dim=n.dim,dados=dados) print(round(c(b1,b2,b3,b4,b5,b6,phi,tau.UT,ll),2)) return(-ll)} vetro <- c(1e-3,1e-5,1e-5,1e-5,1e-5,1e-4,1e-3,1e-3) ajuste.5 = mle2(mod5,start=list(b1=0, b2=0,b3=0,b4=0,b5=0,b6=0,phi=3,tau.UT=-1), method="BFGS",control=list(ndeps=vetro,lmm=1), data=list(integral="LAPLACE",pontos=1,otimizador="BFGS", n.dim=5,dados=dados)) summary(ajuste.5) ########################################################################################## ## Model 6 - Int + LOCAL + TRIM + UHE + UHE:TRIM ######################################### ########################################################################################## mod6 <- function(b1, b2,b3,b4,b5,b6, phi, tau.U, tau.UT, integral, pontos, otimizador, n.dim, dados){ ll = verossimilhanca(modelo = modelo, formu.X="~LOCAL + TRIM",formu.Z="~ UM + TRIM -1", beta.fixo = c(b1,b2,b3,b4,b5,b6), prec.pars=c(phi,tau.U,tau.UT), integral=integral, pontos=pontos, otimizador=otimizador,n.dim=n.dim,dados=dados) print(round(c(b1,b2,b3,b4,b5,b6,phi,tau.U,tau.UT,ll),2)) return(-ll)} vetro <- c(1e-3,1e-4,1e-4,1e-5,1e-5,1e-5,1e-3,1e-3,1e-3) ajuste.6 = mle2(mod6,start=list(b1=0, b2=0,b3=0,b4=0,b5=0,b6=0, phi= 1,tau.U= 0,tau.UT= 0), method="BFGS",control=list(ndeps=vetro,lmm=5), data=list(integral="LAPLACE",pontos=1,otimizador="BFGS", n.dim=5,dados=dados)) summary(ajuste.6) logLik(ajuste.6) cbind(logLik(ajuste.1),logLik(ajuste.2),logLik(ajuste.3),logLik(ajuste.4),logLik(ajuste.5),logLik(ajuste.6)) ########################################################################################## ## Profile likelihood #################################################################### ########################################################################################## ## Time consuming perfil.phi <- profile(ajuste.5,which="phi") plot(perfil.phi) ## Time consuming perfil.tau.UT <- profile(ajuste.5,which="tau.UT") plot(perfil.tau.UT) ########################################################################################## ## Analysis with cloned data ############################################################# ########################################################################################## ## Data set in suitable way to use jags dados.id <- split(dados,dados$ID) dados.id[[4]] <- rbind(dados.id[[4]],c("Caxias",NA,"MONT",2,1)) dados.id[[4]]$UM <- as.numeric(dados.id[[4]]$UM) dados.id[[8]] dados.id[[8]] <- rbind(dados.id[[8]],c("Guaricana",NA,"JUSA",4,1)) dados.id[[9]]$TRIM[10:12] <- 4 dados.id[[7]]$UM <- as.numeric(dados.id[[7]]$UM) dados.id[[8]]$UM <- as.numeric(dados.id[[8]]$UM) ## Matrizes em branco X.int <- matrix(NA,ncol=16,nrow=12) X.reser <- matrix(NA,ncol=16,nrow=12) X.jusa <- matrix(NA,ncol=16,nrow=12) X.trim2 <- matrix(NA,ncol=16,nrow=12) X.trim3 <- matrix(NA,ncol=16,nrow=12) X.trim4 <- matrix(NA,ncol=16,nrow=12) #Z.int <- matrix(NA,ncol=16,nrow=12) Z.trim1 <- matrix(NA,ncol=16,nrow=12) Z.trim2 <- matrix(NA,ncol=16,nrow=12) Z.trim3 <- matrix(NA,ncol=16,nrow=12) Z.trim4 <- matrix(NA,ncol=16,nrow=12) Y <- matrix(NA,ncol=16,nrow=12) ## Model matrices for(i in 1:length(dados.id)){ Y[,i] <- dados.id[[i]]$y X <- model.matrix(~LOCAL + TRIM, data=dados.id[[i]]) Z <- model.matrix(~ UM + TRIM -1,data=dados.id[[i]]) X.int[,i] <- X[,1] X.reser[,i] <- X[,2] X.jusa[,i] <- X[,3] X.trim2[,i] <- X[,4] X.trim3[,i] <- X[,5] X.trim4[,i] <- X[,6] #Z.int[,i] <- Z[,1] Z.trim1[,i] <- Z[,2] Z.trim2[,i] <- Z[,3] Z.trim3[,i] <- Z[,4] Z.trim4[,i] <- Z[,5] print(i) } ########################################################################################## ## Data set suitable way to use jags like a list ######################################### ########################################################################################## dat.jags <- list(Y =t(Y), X.int = t(X.int), X.reser = t(X.reser),X.jusa = t(X.jusa), X.trim2 = t(X.trim2), X.trim3 = t(X.trim3), X.trim4 = t(X.trim4), Z.trim1 = t(Z.trim1), Z.trim2=t(Z.trim2), Z.trim3 = t(Z.trim3), Z.trim4 = t(Z.trim4),n.bloco=16,n.rep=12) ## Final model IQA.bayes <- jags.fit(dat.jags,c("beta0","beta1","beta2","beta3","beta4","beta5","phi","tau.UV"), beta.estruturado.IQA2,n.adapt = 1000, n.update = 500, thin = 5, n.iter= 5000, n.chains=3) plot(IQA.bayes) summary(IQA.bayes) #### Fitting via dclone k <- c(1,5,10,20,30,40,50) IQA.clone <- dc.fit(dat.jags,c("beta0","beta1","beta2","beta3","beta4","beta5","phi","tau.UV"), beta.estruturado.IQA2,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(IQA.clone) summary(IQA.clone) ## Estimability diagnostic esti <- dctable(IQA.clone) plot(esti) plot(esti,type="log")