########################################################################################## ## Bayesian analysis for a class of Beta mixed models #################################### ########################################################################################## ## Carregando os pacotes extras require(betareg) require(INLA) ## Carregando a base de dados do IQVT ## Lendo os dados dados <- read.table("dadosBeta.txt",header=TRUE)[,c(1,5,7,8)] levels(dados$ESTADO) <- c("AM","CE","DF","MT","MS","PB","PR","RO","RR") dados$RENDA <- log(dados$RENDA) - mean(log(dados$RENDA)) names(dados) <- c("y","RENDA","PORTE","ID") ## Analysis descriptive mod.graf <- betareg(y ~ RENDA,data=dados) RENDA <- data.frame(RENDA = seq(-0.5,1.1,l=100)) preditos = predict(mod.graf,newdata=RENDA) pch.id <- unique(as.numeric(dados$ID)) mat <- matrix(c(1,1,2,3),1,4) layout(mat) pch.id <- unique(as.numeric(dados$ID)) par(mar=c(2.7, 2.5, 0.1, 0.5),mgp = c(1.7, 0.7, 0)) plot(y ~ RENDA,pch=pch.id,data=dados,xlim=c(-0.5,1.1),ylab="IQVT",xlab="log(Income)") lines(preditos ~ as.numeric(as.matrix(RENDA))) legend("bottomright",legend=levels(dados$ID),pch=pch.id, cex=0.8) levels(dados$PORTE) <- c("Large","Medium","Small") boxplot(y ~ PORTE,varwidth=TRUE,ylab="",data=dados,xlab="Size") boxplot(y ~ ID,varwidth=TRUE,data=dados,xlab="States",ylab="") ########################################################################################## ## Fitting the models with R-INLA package ################################################ ########################################################################################## ## Model 1 - y ~ Intecerpt formu1 <- y ~ 1 mod1 <- inla(formu1, family="beta", control.family=(list(hyper=list(theta = list(prior = "loggamma", param = c(1,0.01))))), control.compute = list(dic = TRUE, cpo = TRUE), data=dados) ## Model 2 - y ~ size formu2 <- y ~ PORTE mod2 <- inla(formu2, family = "beta", control.family=(list(hyper=list(theta = list(prior = "loggamma", param = c(1,0.01))))), control.compute = list(dic = TRUE, cpo = TRUE), data=dados) ## Model 3 - y ~ size + income formu3 <- y ~ PORTE + RENDA mod3 <- inla(formu3, family = "beta", control.family=(list(hyper=list(theta = list(prior = "loggamma", param = c(1,0.01))))), control.compute = list(dic = TRUE, cpo = TRUE),data=dados) ## Model 4 - y ~ size + income + b_i dados$id <- NA dados$id2 <- NA dados[which(dados$ID == "AM"),]$id <- 1 dados[which(dados$ID == "CE"),]$id <- 2 dados[which(dados$ID == "DF"),]$id <- 3 dados[which(dados$ID == "MT"),]$id <- 4 dados[which(dados$ID == "MS"),]$id <- 5 dados[which(dados$ID == "PB"),]$id <- 6 dados[which(dados$ID == "PR"),]$id <- 7 dados[which(dados$ID == "RO"),]$id <- 8 dados[which(dados$ID == "RR"),]$id <- 9 dados[which(dados$ID == "AM"),]$id2 <- 10 dados[which(dados$ID == "CE"),]$id2 <- 11 dados[which(dados$ID == "DF"),]$id2 <- 12 dados[which(dados$ID == "MT"),]$id2 <- 13 dados[which(dados$ID == "MS"),]$id2 <- 14 dados[which(dados$ID == "PB"),]$id2 <- 15 dados[which(dados$ID == "PR"),]$id2 <- 16 dados[which(dados$ID == "RO"),]$id2 <- 17 dados[which(dados$ID == "RR"),]$id2 <- 18 formu4 <- y ~ PORTE + RENDA + f(id, model="iid", param=c(0.5, 0.001487)) mod4 <- inla(formu4, family="beta", control.family=(list(hyper=list(theta = list(prior = "loggamma", param = c(1,0.01))))), control.compute = list(dic = TRUE, cpo = TRUE), data=dados) ## Model 5 - y ~ size + income + b_{i1} + b_{i2}*RENDA N = 18 formu5 <- y ~ PORTE + RENDA + f(id, model="iid2d", param=c(5, 0.001487, 0.005, 0), n = N, diagonal=0) + f(id2, RENDA , copy="id") mod5 <- inla(formu5, family="beta", control.family=(list(hyper=list(theta = list(prior = "loggamma", param = c(1,0.01))))), control.compute = list(dic = TRUE, cpo = TRUE), data=dados) summary(mod5) ## Comparando os modelos ## Verossimilhança marginal data.frame("M1" = mod1$mlik[1], "M2" = mod2$mlik[1], "M3" = mod3$mlik[1], "M4" = mod4$mlik[1], "M5" = mod5$mlik[1]) ## DIC data.frame("M1" = mod1$dic$dic, "M2" = mod2$dic$dic, "M3" = mod3$dic$dic, "M4" = mod4$dic$dic, "M5" = mod5$dic$dic) ## Cpo data.frame("M1" = -mean(log(mod1$cpo$cpo)), "M2" = -mean(log(mod2$cpo$cpo)), "M3" = -mean(log(mod3$cpo$cpo)), "M4" = -mean(log(mod4$cpo$cpo)), "M5" = -mean(log(mod5$cpo$cpo))) ## PIT par(mfrow=c(2,2)) hist(mod2$cpo$pit,breaks=seq(0,1,l=50), main="", xlab="Model 2") hist(mod3$cpo$pit,breaks=seq(0,1,l=50),main="", xlab="Model 3") hist(mod4$cpo$pit,breaks=seq(0,1,l=50),main="", xlab="Model 4") hist(mod5$cpo$pit,breaks=seq(0,1,l=50),main="", xlab="Model 5") ########################################################################################## ## Ajustando os modelos no JAGS ########################################################## ########################################################################################## ## Ajustando via JAGS require(dclone) ## Carregando os modelos e os dados source("modelJAGS.R") bayes.mod4 <- jags.fit(dat.jags,c("beta1","beta2","beta3","beta4","phi", "tau"),modelo4,n.adapt = 10000, n.update = 5000, thin = 100, n.iter= 500000, n.chains=3) summary(bayes.mod4) load("BetaBayes1502.RData") ########################################################################################## ## Comparando as posteriori ############################################################## ########################################################################################## par(mfrow=c(2,3),mar=c(2.7, 2.5, 0.1, 0.5),mgp = c(1.7, 0.7, 0)) ## Beta 0 plot(mod4$marginals.fixed$"(Intercept)"[,1],mod4$marginals.fixed$"(Intercept)"[,2], ylab="Density", xlab= expression(beta[0]),type="l") hist(c(bayes.mod4[[1]][,1],bayes.mod4[[2]][,1],bayes.mod4[[3]][,1]),add=TRUE, prob=TRUE, breaks=seq(0.1,0.8,l=80)) abline(v=c(0.2918,0.40, 0.4978), col="red", lty=2,lwd=2) ## Beta 1 plot(mod4$marginals.fixed$"PORTEMedium"[,1], mod4$marginals.fixed$"PORTEMedium"[,2], ylab="Density", xlab= expression(beta[1]),type="l") hist(c(bayes.mod4[[1]][,2],bayes.mod4[[2]][,2],bayes.mod4[[3]][,2]), prob=TRUE, add=TRUE, breaks = seq(-0.3,0.20,l=80)) abline(v=c(-0.1275,-0.07,-0.0172), col="red", lty=2,lwd=2) ## Beta 2 plot(mod4$marginals.fixed$"PORTESmall"[,1],mod4$marginals.fixed$"PORTESmall"[,2], ylab="Density", xlab= expression(beta[2]),ylim=c(0,14),type="l") hist(c(bayes.mod4[[1]][,3],bayes.mod4[[2]][,3],bayes.mod4[[3]][,3]), prob=TRUE, add=TRUE, breaks = seq(-0.3,0.1,l=80)) abline(v=c(-0.1910,-0.13,-0.0741), col="red", lty=2,lwd=2) ## Beta 3 plot(mod4$marginals.fixed$"RENDA"[,1],mod4$marginals.fixed$"RENDA"[,2], ylab="Density", xlab= expression(beta[3]),type="l") hist(c(bayes.mod4[[1]][,4],bayes.mod4[[2]][,4],bayes.mod4[[3]][,4]),add=TRUE, prob=TRUE, breaks=seq(0.2,0.8,l=80)) abline(v=c(0.3931,0.47,0.5480), col="red",lty=2,lwd=2) ## Phi plot(mod4$marginals.hyperpar[[1]][,1],mod4$marginals.hyperpar[[1]][,2], ylab="Density", xlim=c(60,130), xlab= expression(phi),type="l") hist(c(bayes.mod4[[1]][,5],bayes.mod4[[2]][,5],bayes.mod4[[3]][,5]),prob=TRUE,add=TRUE,breaks=seq(50,150,l=80)) abline(v=c(81.08,94.19,108.6460), col="red",lty=2,lwd=2) ## Tau plot(mod4$marginals.hyperpar[[2]][,1],mod4$marginals.hyperpar[[2]][,2], ylab="Density", xlim=c(0,300),xlab= expression(tau),type="l") hist(c(bayes.mod4[[1]][,6],bayes.mod4[[2]][,6],bayes.mod4[[3]][,6]),add=TRUE,prob=TRUE,breaks=seq(0,500,l=80)) abline(v=c(19.7383,62.35,156.4794), col="red", lty=2,lwd=2) ########################################################################################## ## Fim do código IQVT Bayes ############################################################## ########################################################################################## ########################################################################################## ## Comparando as coberturas do IC de perfil com o MCMC e INLA ############################ ########################################################################################## source("auxfunctions.R") ## data frame with profile interval ic.perfil <- data.frame("Param" = c("b0","b1", "b2","b3","phi", "tau"), "Lower" = c(0.2918,-0.1275,-0.1910,0.3931,81.08,19.7283),"Upper" = c(0.4978,-0.0172,-0.0741,0.5480,108.6460,156.4794)) ## MCMC mc.b0 <- cobertura.mcmc(obj.mcmc = bayes.mod4, cod.param = 1, ic.perfil = ic.perfil) mc.b1 <- cobertura.mcmc(obj.mcmc = bayes.mod4, cod.param = 2, ic.perfil = ic.perfil) mc.b2 <- cobertura.mcmc(obj.mcmc = bayes.mod4, cod.param = 3, ic.perfil = ic.perfil) mc.b3 <- cobertura.mcmc(obj.mcmc = bayes.mod4, cod.param = 4, ic.perfil = ic.perfil) mc.phi <- cobertura.mcmc(obj.mcmc = bayes.mod4, cod.param = 5, ic.perfil = ic.perfil) mc.tau <- cobertura.mcmc(obj.mcmc = bayes.mod4, cod.param = 6, ic.perfil = ic.perfil) ## INLA inla.b0 <- cobertura.inla(marginal.inla = mod4$marginals.fixed[[1]], cod.param = 1, ic.perfil) inla.b1 <- cobertura.inla(marginal.inla = mod4$marginals.fixed[[2]], cod.param = 2, ic.perfil) inla.b2 <- cobertura.inla(marginal.inla = mod4$marginals.fixed[[3]], cod.param = 3, ic.perfil) inla.b3 <- cobertura.inla(marginal.inla = mod4$marginals.fixed[[4]], cod.param = 4, ic.perfil) inla.phi <- cobertura.inla(marginal.inla = mod4$marginals.hyperpar[[1]], cod.param = 5, ic.perfil) inla.tau <- cobertura.inla(marginal.inla = mod4$marginals.hyperpar[[2]], cod.param = 6, ic.perfil) ## Samples MCMC in credibility interval from INLA c(as.numeric(mod4$summary.fixed[,3]),as.numeric(mod4$summary.hyperpar[,3])) ic.inla <- data.frame("Param" = c("b0","b1", "b2","b3","phi", "tau"), "Lower" = as.numeric(c(mod4$summary.fixed[,3],mod4$summary.hyperpar[,3])), "Upper" = as.numeric(c(mod4$summary.fixed[,5],mod4$summary.hyperpar[,5]))) mc.inla.b0 <- cobertura.mcmc(obj.mcmc = bayes.mod4, cod.param = 1, ic.perfil = ic.inla) mc.inla.b1 <- cobertura.mcmc(obj.mcmc = bayes.mod4, cod.param = 2, ic.perfil = ic.inla) mc.inla.b2 <- cobertura.mcmc(obj.mcmc = bayes.mod4, cod.param = 3, ic.perfil = ic.inla) mc.inla.b3 <- cobertura.mcmc(obj.mcmc = bayes.mod4, cod.param = 4, ic.perfil = ic.inla) mc.inla.phi <- cobertura.mcmc(obj.mcmc = bayes.mod4, cod.param = 5, ic.perfil = ic.inla) mc.inla.tau <- cobertura.mcmc(obj.mcmc = bayes.mod4, cod.param = 6, ic.perfil = ic.inla) ## Comparing the coverage rate data.frame("MCMC-Profile" = c(mc.b0,mc.b1,mc.b2,mc.b3,mc.phi,mc.tau),"INLA-Profile" = c(inla.b0,inla.b1,inla.b2,inla.b3,inla.phi,inla.tau), "INLA-MCMC" = c(mc.inla.b0,mc.inla.b1,mc.inla.b2,mc.inla.b3,mc.inla.phi,mc.inla.tau)) ########################################################################################## ## Analysis sensibility for Beta mixed models ############################################ ########################################################################################## ## Loading extra functions source("BetaEspacial/sensibilityFunctions.R") ########################################################################################## ## I used the priori G(a = 1, b = 0.01) for phi and G(a=0.5,b = 0.001487) for tau ######## ## Choosing the priori's distributions for the sensibility study ######################### ## The priori default that I propose is G(a1=1, b = 0.01), ############################### ## I will choosing priori with divergence between 0.1 and 0.6 ############################ ########################################################################################## hellinger.gamma(a1=1,b1=0.01,a2=1,b2=0.0135) hellinger.gamma(a1=1,b1=0.01,a2=1,b2=0.0178) hellinger.gamma(a1=1,b1=0.01,a2=1,b2=0.0242) hellinger.gamma(a1=1,b1=0.01,a2=1,b2=0.0338) hellinger.gamma(a1=1,b1=0.01,a2=1,b2=0.05) hellinger.gamma(a1=1,b1=0.01,a2=1,b2=0.0765) ## First I will fix the priori default for the tau and vary the priori for the phi a <- c(1,1+0.01,1-0.01) b <- c(c(0.01,0.0135,0.0178,0.0242,0.0338,0.05,0.0765),c(0.01,0.0135,0.0178,0.0242,0.0338,0.05,0.0765)+0.0001, c(0.01,0.0135,0.0178,0.0242,0.0338,0.05,0.0765)-0.0001) gridi.phi <- expand.grid(a,b) gridi.phi ## Fitting the model modelo.phi <- list() formu4 <- y ~ PORTE + RENDA + f(id, model="iid", param=c(0.5, 0.001487)) for(i in 1:dim(gridi.phi)[1]){ modelo.phi[i][[1]] <- inla(formu4, family="beta", control.family=(list(hyper=list(theta = list(prior = "loggamma", param = c(gridi.phi[i,1],gridi.phi[i,2]))))), control.compute = list(dic = TRUE, cpo = TRUE), data=dados) } ## Evaluating the divergence between the posteriori distributions hl.post <- matrix(NA,ncol=6,nrow=63) for(i in 1:63){ b0 = hellinger.inla(approx.1 = mod4$marginals.fixed[[1]], approx.2 = modelo.phi[i][[1]]$marginals.fixed[[1]], lower=-Inf, upper=Inf) b1 = hellinger.inla(approx.1 = mod4$marginals.fixed[[2]], approx.2 = modelo.phi[i][[1]]$marginals.fixed[[2]], lower=-Inf, upper=Inf) b2 = hellinger.inla(approx.1 = mod4$marginals.fixed[[3]], approx.2 = modelo.phi[i][[1]]$marginals.fixed[[3]], lower=-Inf, upper=Inf) b3 = hellinger.inla(approx.1 = mod4$marginals.fixed[[4]], approx.2 = modelo.phi[i][[1]]$marginals.fixed[[4]], lower=-Inf, upper=Inf) phi = hellinger.inla(approx.1 = mod4$marginals.hyperpar[[1]], approx.2 = modelo.phi[i][[1]]$marginals.hyperpar[[1]], lower=0, upper=5000) tau = hellinger.inla(approx.1 = mod4$marginals.hyperpar[[2]], approx.2 = modelo.phi[i][[1]]$marginals.hyperpar[[2]], lower=0, upper=5000) hl.post[i,] <- c(b0,b1,b2,b3,phi,tau) } ## The measure the sensibility is defined as the ration between the Hellinger divergence between the posteriori and DH between the priori hl.priori <- c() for(i in 1:63){ hl.priori[i] <- hellinger.gamma(a1=1,b1=0.01,a2=gridi.phi$Var1[i],b2=gridi.phi$Var2[i])} resul.phi <- data.frame(gridi.phi, "Hl Post" = hl.post[,5], "Hl priori" = hl.priori, "S" = hl.post[,5]/hl.priori) ## Table comparative between the divergence tab.phi <- resul.phi[c(4,7,10,13,16,19),] save(tab.phi,file = "tab.phi.RData") ########################################################################################## ### Graphical comparing priori with posteriori ########################################### ########################################################################################## y.phi <- c(modelo.phi[[1]]$marginals.hyperpar[[1]][,2],modelo.phi[[4]]$marginals.hyperpar[[1]][,2],modelo.phi[[7]]$marginals.hyperpar[[1]][,2],modelo.phi[[10]]$marginals.hyperpar[[1]][,2], modelo.phi[[13]]$marginals.hyperpar[[1]][,2], modelo.phi[[16]]$marginals.hyperpar[[1]][,2], modelo.phi[[19]]$marginals.hyperpar[[1]][,2]) x.phi <- c(modelo.phi[[1]]$marginals.hyperpar[[1]][,1],modelo.phi[[4]]$marginals.hyperpar[[1]][,1],modelo.phi[[7]]$marginals.hyperpar[[1]][,1],modelo.phi[[10]]$marginals.hyperpar[[1]][,1], modelo.phi[[13]]$marginals.hyperpar[[1]][,1], modelo.phi[[16]]$marginals.hyperpar[[1]][,1], modelo.phi[[19]]$marginals.hyperpar[[1]][,1]) y.tau <- c(modelo.phi[[1]]$marginals.hyperpar[[2]][,2],modelo.phi[[4]]$marginals.hyperpar[[2]][,2],modelo.phi[[7]]$marginals.hyperpar[[2]][,2],modelo.phi[[10]]$marginals.hyperpar[[2]][,2], modelo.phi[[13]]$marginals.hyperpar[[2]][,2], modelo.phi[[16]]$marginals.hyperpar[[2]][,2], modelo.phi[[19]]$marginals.hyperpar[[2]][,2]) x.tau <- c(modelo.phi[[1]]$marginals.hyperpar[[2]][,1],modelo.phi[[4]]$marginals.hyperpar[[2]][,1],modelo.phi[[7]]$marginals.hyperpar[[2]][,1],modelo.phi[[10]]$marginals.hyperpar[[2]][,1], modelo.phi[[13]]$marginals.hyperpar[[2]][,1], modelo.phi[[16]]$marginals.hyperpar[[2]][,1], modelo.phi[[19]]$marginals.hyperpar[[2]][,1]) HL.phi <- rep(c(0,0.1,0.2,0.3,0.4,0.5,0.6),each=61) y.priori.phi <- c(curve(dgamma(x, shape=1, scale=1/0.01),0,150)$y,curve(dgamma(x, shape=1, scale=1/0.0135),0,150)$y,curve(dgamma(x, shape=1, scale=1/0.0178),0,150)$y, curve(dgamma(x, shape=1, scale=1/0.0242),0,150)$y,curve(dgamma(x, shape=1, scale=1/0.0338),0,150)$y,curve(dgamma(x, shape=1, scale=1/0.05),0,150)$y, curve(dgamma(x, shape=1, scale=1/0.0765),0,150)$y) x.priori.phi <- c(curve(dgamma(x, shape=1, scale=1/0.01),0,150)$x,curve(dgamma(x, shape=1, scale=1/0.0135),0,150)$x,curve(dgamma(x, shape=1, scale=1/0.0178),0,150)$x, curve(dgamma(x, shape=1, scale=1/0.0242),0,150)$x,curve(dgamma(x, shape=1, scale=1/0.0338),0,150)$x,curve(dgamma(x, shape=1, scale=1/0.05),0,150)$x, curve(dgamma(x, shape=1, scale=1/0.0765),0,150)$x) HL.priori.phi <- rep(c(0,0.1,0.2,0.3,0.4,0.5,0.6),each=101) Dist <- c(rep("Posteriori",length(y.phi)),rep("Priori",length(y.priori.phi))) data.tau <- data.frame("y.tau" = y.tau, "x.tau" = x.tau, "HL.phi" = HL.phi) data.graf <- data.frame("y.phi" = c(y.phi, y.priori.phi), "x.phi" = c(x.phi,x.priori.phi), "HL.phi" = c(HL.phi, HL.priori.phi), "Dist" = Dist) my.col <- c(1,1,1,1,"gray50","gray50","gray50") my.lty <- c(1,1,2,3,1,2,3) my.lwd <- c(2,1,1,1,1,1,1) ps <- list(superpose.line=list(col=my.col, lty=my.lty, lwd = my.lwd)) p.phi <- xyplot(y.phi ~ x.phi | Dist, groups= HL.phi, data=data.graf, type="l", scales=list(x="free"),auto.key=list(columns=4,points=FALSE, lines=TRUE), ylab="Density", xlab=expression(phi), par.settings=ps, strip = strip.custom(bg="gray90")) p.tau <- xyplot(y.tau ~ x.tau,groups=HL.phi, type="l", scales=list(x="free"), ylab="", xlab=expression(tau), par.settings=ps, strip = strip.custom(bg="gray90")) ########################################################################################## ## Sensibility for tau ################################################################### ########################################################################################## ## For the parameter tau, precision of the random effect the priori default is G(0.5,0.001487), ## the same form I will choosing priori with divergence the 0.1 until 0.6 hellinger.gamma(a1=0.5,b1=0.001487, a2=0.5,b2=0.00225) hellinger.gamma(a1=0.5,b1=0.001487, a2=0.5,b2=0.0035) hellinger.gamma(a1=0.5,b1=0.001487, a2=0.5,b2=0.0055) hellinger.gamma(a1=0.5,b1=0.001487, a2=0.5,b2=0.0088) hellinger.gamma(a1=0.5,b1=0.001487, a2=0.5,b2=0.016) hellinger.gamma(a1=0.5,b1=0.001487, a2=0.5,b2=0.033) ## I will fix the priori default for the tau and vary the priori for the tau a <- c(0.5,0.5+0.01,0.5-0.01) b <- c(c(0.001487,0.00225,0.0035,0.0055,0.0088,0.016,0.033),c(0.001487,0.00225,0.0035,0.0055,0.0088,0.016,0.033)+0.00001,c(0.001487,0.00225,0.0035,0.0055,0.0088,0.016,0.033)-0.00001) gridi.tau <- expand.grid(a,b) ## Fitting the model modelo.tau <- list() for(i in 1:63){ formu4 <- y ~ PORTE + RENDA + f(id, model="iid", param=c(gridi.tau[i,1], gridi.tau[i,2])) modelo.tau[i][[1]] <- inla(formu4, family="beta", control.family=(list(hyper=list(theta = list(prior = "loggamma", param = c(1,0.01))))), control.compute = list(dic = TRUE, cpo = TRUE), data=dados) } ## Evaluating the divergence between the posteriori distributions hl.post.tau <- matrix(NA,ncol=6,nrow=63) for(i in 1:63){ b0 = hellinger.inla(approx.1 = mod4$marginals.fixed[[1]], approx.2 = modelo.tau[i][[1]]$marginals.fixed[[1]], lower=-Inf, upper=Inf) b1 = hellinger.inla(approx.1 = mod4$marginals.fixed[[2]], approx.2 = modelo.tau[i][[1]]$marginals.fixed[[2]], lower=-Inf, upper=Inf) b2 = hellinger.inla(approx.1 = mod4$marginals.fixed[[3]], approx.2 = modelo.tau[i][[1]]$marginals.fixed[[3]], lower=-Inf, upper=Inf) b3 = hellinger.inla(approx.1 = mod4$marginals.fixed[[4]], approx.2 = modelo.tau[i][[1]]$marginals.fixed[[4]], lower=-Inf, upper=Inf) phi = hellinger.inla(approx.1 = mod4$marginals.hyperpar[[1]], approx.2 = modelo.tau[i][[1]]$marginals.hyperpar[[1]], lower=0, upper=5000) tau = hellinger.inla(approx.1 = mod4$marginals.hyperpar[[2]], approx.2 = modelo.tau[i][[1]]$marginals.hyperpar[[2]], lower=0, upper=5000) hl.post.tau[i,] <- c(b0,b1,b2,b3,phi,tau) } hl.post.tau ## The measure the sensibility is defined as the ration between the Hellinger divergence between the posteriori and DH between the priori hl.priori.tau <- c() for(i in 1:63){ hl.priori.tau[i] <- hellinger.gamma(a1=0.5,b1=0.001487,a2=gridi.tau$Var1[i],b2=gridi.tau$Var2[i])} resul.tau <- data.frame(gridi.tau, "Hl Post" = hl.post.tau[,6], "Hl priori" = hl.priori.tau, "S" = hl.post.tau[,6]/hl.priori.tau) xyplot(S ~ Var2, grouped = resul.tau$Var1, type="p", data=resul.tau) ## Table comparative between priori and posteriori resul.tau[c(4,7,10,13,16,19),] save(modelo.phi,file = "modelo.phi.RData") save(modelo.tau,file="modelo.tau.RData") ########################################################################################## ### Graphical comparing priori with posteriori ########################################### ########################################################################################## y.tau1 <- c(modelo.tau[[1]]$marginals.hyperpar[[2]][,2],modelo.tau[[4]]$marginals.hyperpar[[2]][,2],modelo.tau[[7]]$marginals.hyperpar[[2]][,2],modelo.tau[[10]]$marginals.hyperpar[[2]][,2], modelo.tau[[13]]$marginals.hyperpar[[2]][,2], modelo.tau[[16]]$marginals.hyperpar[[2]][,2], modelo.tau[[19]]$marginals.hyperpar[[2]][,2]) x.tau1 <- c(modelo.tau[[1]]$marginals.hyperpar[[2]][,1],modelo.tau[[4]]$marginals.hyperpar[[2]][,1],modelo.tau[[7]]$marginals.hyperpar[[2]][,1],modelo.tau[[10]]$marginals.hyperpar[[2]][,1], modelo.tau[[13]]$marginals.hyperpar[[2]][,1], modelo.tau[[16]]$marginals.hyperpar[[2]][,1], modelo.tau[[19]]$marginals.hyperpar[[2]][,1]) y.phi1 <- c(modelo.tau[[1]]$marginals.hyperpar[[1]][,2],modelo.tau[[4]]$marginals.hyperpar[[1]][,2],modelo.tau[[7]]$marginals.hyperpar[[1]][,2],modelo.tau[[10]]$marginals.hyperpar[[1]][,2], modelo.tau[[13]]$marginals.hyperpar[[1]][,2], modelo.tau[[16]]$marginals.hyperpar[[1]][,2], modelo.tau[[19]]$marginals.hyperpar[[1]][,2]) x.phi1 <- c(modelo.tau[[1]]$marginals.hyperpar[[1]][,1],modelo.tau[[4]]$marginals.hyperpar[[1]][,1],modelo.tau[[7]]$marginals.hyperpar[[1]][,1],modelo.tau[[10]]$marginals.hyperpar[[1]][,1], modelo.tau[[13]]$marginals.hyperpar[[1]][,1], modelo.tau[[16]]$marginals.hyperpar[[1]][,1], modelo.tau[[19]]$marginals.hyperpar[[1]][,1]) HL.tau <- rep(c(0,0.1,0.2,0.3,0.4,0.5,0.6),each=61) y.priori.tau <- c(curve(dgamma(x, shape=1, scale=1/0.001487),0,200)$y,curve(dgamma(x, shape=1, scale=1/0.00225),0,200)$y,curve(dgamma(x, shape=1, scale=1/0.0035),0,200)$y, curve(dgamma(x, shape=1, scale=1/0.0055),0,200)$y,curve(dgamma(x, shape=1, scale=1/0.0088),0,200)$y,curve(dgamma(x, shape=1, scale=1/0.016),0,200)$y, curve(dgamma(x, shape=1, scale=1/0.033),0,200)$y) x.priori.tau <- c(curve(dgamma(x, shape=1, scale=1/0.001487),0,200)$x,curve(dgamma(x, shape=1, scale=1/0.00225),0,200)$x,curve(dgamma(x, shape=1, scale=1/0.0035),0,200)$x, curve(dgamma(x, shape=1, scale=1/0.0055),0,200)$x,curve(dgamma(x, shape=1, scale=1/0.0088),0,200)$x,curve(dgamma(x, shape=1, scale=1/0.016),0,200)$x, curve(dgamma(x, shape=1, scale=1/0.033),0,200)$x) HL.priori.tau <- rep(c(0,0.1,0.2,0.3,0.4,0.5,0.6),each=101) Dist <- c(rep("Posteriori",length(y)),rep("Priori",length(y.priori))) data.graf <- data.frame("y" = c(y.tau1, y.priori.tau), "x" = c(x.tau1,x.priori.tau), "HL" = c(HL.tau, HL.priori.tau), "Dist" = Dist) data.tau <- data.frame("y.phi" = y.phi1, "x.phi" = x.phi1, "HL.tau" = HL.tau) my.col <- c(1,1,1,1,"gray50","gray50","gray50") my.lty <- c(1,1,2,3,1,2,3) my.lwd <- c(2,1,1,1,1,1,1) ps <- list(superpose.line=list(col=my.col, lty=my.lty, lwd = my.lwd)) p1.tau <- xyplot(y ~ x | Dist, groups= HL, data=data.graf,xlim=c(0,200), type="l", scales=list(x="free"), #auto.key=list(columns=4,points=FALSE, lines=TRUE), ylab="Density", xlab=expression(tau), par.settings=ps, strip = strip.custom(bg="gray90")) p1.phi <- xyplot(y.phi ~ x.phi, groups= HL.tau, data=data.tau, type="l", scales=list(x="free"), ylab="Density", xlab=expression(phi), par.settings=ps, strip = strip.custom(bg="gray90")) ## plot(p.phi,position=c(0,0,0.63,0.5),more=TRUE) plot(p.tau,position=c(0.6,0,1,0.45),more=TRUE) plot(p1.tau,position=c(0,0.5,0.63,1),more=TRUE) plot(p1.phi,position=c(0.6,0.5,1,0.95),more=FALSE) ########################################################################################## ## Repetindo as análises com as priori's que mais causaram impacto ####################### ########################################################################################## resul.phi[c(4,7,10,13,16,19),] tab.tau <- resul.tau[c(4,7,10,13,16,19),] save(tab.tau,file="tab.tau.RData") ## Comparando as estimativas pontuais tt <- data.frame("Default" = c(mod4$summary.fixed[,1],mod4$summary.hyperpar[,1]), "Std1" = c(mod4$summary.fixed[,2],mod4$summary.hyperpar[,2]), "phi" = c(modelo.phi[[19]]$summary.fixed[,1],modelo.phi[[19]]$summary.hyperpar[,1]), "Std2"= c(modelo.phi[[19]]$summary.fixed[,2],modelo.phi[[19]]$summary.hyperpar[,2]), "tau" = c(modelo.tau[[19]]$summary.fixed[,1],modelo.tau[[19]]$summary.hyperpar[,1]), "Std3" = c(modelo.tau[[19]]$summary.fixed[,2],modelo.tau[[19]]$summary.hyperpar[,2])) rownames(tt) <- c("b0","b1", "b2","b3","phi", "tau") tt save(tt,file="tt.RData")