Rjags_Laboratorio
require(readxl)
require(tidyverse)
require(rjags)
require(coda)
require(MCMCpack)
require(bayesplot)Em construção: problema com convergencia
- algumas tentativas feitas - sem sucesso
- Transformação box-cox para tx.latrc;
- Modelo glm com dist=gamma & link=exp (implementa-se na sintaxe do Rjags)
- Modelo lognormal para tx.latrc.
Planilha de dados
- https://rstudio-pubs-static.s3.amazonaws.com/34971_bf47fd9d008246a78d2511810f34922e.html
- https://stackoverflow.com/questions/53262766/jags-multiple-linear-regression-with-yi-gamma-bayesian
- https://people.stat.sc.edu/hansont/stat740/jags_user_manual.pdf
dados=read_excel("dados_renato.xlsx",sheet=1)
names(dados)=c("Ano","UF","Cod.UF","tx.latrc","tx.pres","gini.ibge","perc.jov.1524","perc.hom","pbf","densidade.urbana1","densidade.urbana2","taxa.casamentos","Taxa.desligamentos","raz.2020","raz.1040")
l1<-1.16673
l2 <-1.931932
dados = dados %>%
mutate(ytrans = ((tx.latrc + l2)^l1-1)/l1,
taxa.casamentos=taxa.casamentos*100000, #por 100 mil habitantes!
Taxa.desligamentos=Taxa.desligamentos*100) #taxa por 100 habitantes!
head(dados)## # A tibble: 6 x 16
## Ano UF Cod.UF tx.latrc tx.pres gini.ibge perc.jov.1524 perc.hom pbf
## <dbl> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2005 Rond~ 1 1.30 298. 0.544 19.7 49.6 62.7
## 2 2006 Rond~ 1 1.60 272. 0.537 20.2 50.7 84.6
## 3 2007 Rond~ 1 1.51 354. 0.472 19.7 49.7 110.
## 4 2008 Rond~ 1 1.14 400. 0.478 19.1 50.2 119.
## 5 2009 Rond~ 1 0.997 465. 0.49 19.6 49.9 139.
## 6 2010 Rond~ 1 1.67 476. NA NA NA 145.
## # ... with 7 more variables: densidade.urbana1 <dbl>, densidade.urbana2 <dbl>,
## # taxa.casamentos <dbl>, Taxa.desligamentos <dbl>, raz.2020 <dbl>,
## # raz.1040 <dbl>, ytrans <dbl>cat("model {
for( i in 1:n)
{
eta[i] <- b[1]*tx.pres[i]+b[2]*gini.ibge[i]+b[3]*perc.jov.1524[i]+b[4]*perc.hom[i]+b[5]*pbf[i]+b[6]*densidade.urbana1[i]+b[7]*Taxa.desligamentos[i]+b[8]*Ano[i]+b[9]*taxa.casamentos[i]+u[UFnum[i],1]*Ano[i]+u[UFnum[i],2]
tx.latrc[i] ~ dnorm( eta[i], tau)
}
#MVN Prior for random intercepts and slopes
for (j in 1:nUF)
{
u[j, 1:2] ~ dmnorm(meanu[], prec.Sigma[,])
}
meanu[1]~dnorm(0, .001)
meanu[2]~dnorm(0, .001)
prec.Sigma[1:2, 1:2] ~ dwish(Omega[,], 2)
Sigma[1:2, 1:2]<-inverse(prec.Sigma[,])
rho12<-Sigma[1,2]/ sqrt(Sigma[1,1]* Sigma[2,2])
#Valores iniciais para a matrix de covariancias
for (j in 1:2){ for (k in 1:2){ Omega[j,k] <-equals(j,k)*.1 } }
#Priori multivariada normal para os coeficientes
b[1:9]~dmnorm(meanb[], prec.beta[,])
for(j in 1:9){
meanb[j]~dnorm(0, .0001)
}
#A Wishart é generalização da Gamma para o caso multivariado
prec.beta[1:9,1:9]~dwish(Obeta[,], 9)
for(j in 1:9){ for(k in 1:9){ Obeta[j,k] <- equals(j,k)*.1 } }
#Priori Gamma para a precisão da normal (dos dados)
tau ~ dgamma(0.001,0.001)
}",file="modelo1_lab.txt")
dados1<-dados[complete.cases(dados),]
nUF<-table(dados1$Cod.UF)
head(nUF) #Number of UF##
## 1 2 3 4 5 6
## 8 8 8 7 8 7dados1$UFnum<-rep(1:length(unique(dados1$Cod.UF)), nUF)
dados2<-list(tx.latrc=dados1$tx.latrc,Ano=dados1$Ano-2000, tx.pres=dados1$tx.pres,gini.ibge=dados1$gini.ibge,perc.jov.1524=dados1$perc.jov.1524,perc.hom=dados1$perc.hom,pbf=dados1$pbf,densidade.urbana1=dados1$densidade.urbana1,Taxa.desligamentos=dados1$Taxa.desligamentos,taxa.casamentos=dados1$taxa.casamentos,n=length(dados1$tx.latrc), UFnum=dados1$UFnum,nUF=length(unique(dados1$Cod.UF)))
lapply(dados2 , summary)## $tx.latrc
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.06595 0.67274 0.92725 1.06438 1.30328 4.78987
##
## $Ano
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 5.000 7.000 9.000 9.631 12.000 14.000
##
## $tx.pres
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 51.71 156.34 224.67 247.80 328.76 634.53
##
## $gini.ibge
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.4290 0.4800 0.5040 0.5074 0.5280 0.6150
##
## $perc.jov.1524
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 9.711 16.821 18.300 18.079 19.400 23.500
##
## $perc.hom
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 46.80 48.40 48.90 49.02 49.70 52.40
##
## $pbf
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 22.04 72.46 135.08 156.68 223.49 423.79
##
## $densidade.urbana1
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 34.29 45.03 63.84 73.63 99.63 160.29
##
## $Taxa.desligamentos
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.196 3.666 5.443 6.899 10.103 17.911
##
## $taxa.casamentos
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 210.8 396.2 485.9 485.8 559.1 801.1
##
## $n
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 225 225 225 225 225 225
##
## $UFnum
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.00 8.00 15.00 14.41 21.00 27.00
##
## $nUF
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 27 27 27 27 27 27#Set some initial values
init.rng1<-list(".RNG.seed" = 1234, ".RNG.name" = "base::Mersenne-Twister", tau=1)
init.rng2<-list(".RNG.seed" = 5678, ".RNG.name" = "base::Mersenne-Twister", tau=5)
#initialize the model
mod1<-jags.model(file = "modelo1_lab.txt", data=dados2, n.chains=2,inits =list(init.rng1, init.rng2) )## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 225
## Unobserved stochastic nodes: 42
## Total graph size: 4586
##
## Initializing model#burn in
update(mod1, 5000)
#collect samples of the parameters
samps<-coda.samples(mod1, variable.names=c( "u","b", "tau", "Sigma", "rho12"), n.iter=100000, n.thin=5)
#effective sample size for each parameter
effectiveSize(samps)## Sigma[1,1] Sigma[2,1] Sigma[1,2] Sigma[2,2] b[1] b[2]
## 60060.03933 50398.78096 50398.78096 31407.89343 11625.82401 410.64111
## b[3] b[4] b[5] b[6] b[7] b[8]
## 23673.32431 99.65523 12821.49782 4432.76172 7793.60050 34.01867
## b[9] rho12 tau u[1,1] u[2,1] u[3,1]
## 17738.80176 72793.52566 51507.01575 34.19321 35.26782 33.37459
## u[4,1] u[5,1] u[6,1] u[7,1] u[8,1] u[9,1]
## 33.64053 32.52244 31.18955 34.01279 34.03364 35.39632
## u[10,1] u[11,1] u[12,1] u[13,1] u[14,1] u[15,1]
## 32.37113 34.12017 33.20623 31.95757 33.54143 36.80153
## u[16,1] u[17,1] u[18,1] u[19,1] u[20,1] u[21,1]
## 34.40408 33.73431 34.47381 32.51550 36.28141 33.80605
## u[22,1] u[23,1] u[24,1] u[25,1] u[26,1] u[27,1]
## 33.77891 33.41121 32.21172 33.93159 32.49731 32.33780
## u[1,2] u[2,2] u[3,2] u[4,2] u[5,2] u[6,2]
## 52.53121 52.26985 50.01875 55.01452 48.00119 50.45733
## u[7,2] u[8,2] u[9,2] u[10,2] u[11,2] u[12,2]
## 52.41599 48.13428 48.10095 54.25454 55.05785 52.74327
## u[13,2] u[14,2] u[15,2] u[16,2] u[17,2] u[18,2]
## 58.31251 55.69996 57.53209 51.09583 49.07004 49.57876
## u[19,2] u[20,2] u[21,2] u[22,2] u[23,2] u[24,2]
## 46.68085 47.65631 46.87304 50.47525 46.37216 49.07358
## u[25,2] u[26,2] u[27,2]
## 50.57598 48.35832 47.78949#Numerical summary of each parameter:
summary(samps, quantiles = c(.025, .05, .95, .975))##
## Iterations = 5001:105000
## Thinning interval = 1
## Number of chains = 2
## Sample size per chain = 1e+05
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## Sigma[1,1] 0.0131751 0.0044173 9.877e-06 1.804e-05
## Sigma[2,1] -0.0990381 0.0407061 9.102e-05 1.814e-04
## Sigma[1,2] -0.0990381 0.0407061 9.102e-05 1.814e-04
## Sigma[2,2] 1.3429230 0.4747631 1.062e-03 2.679e-03
## b[1] 0.0005835 0.0008397 1.878e-06 7.794e-06
## b[2] -2.9586538 1.8140545 4.056e-03 8.886e-02
## b[3] -0.0135856 0.0277427 6.203e-05 1.805e-04
## b[4] 0.0935164 0.0660053 1.476e-04 6.601e-03
## b[5] -0.0012602 0.0020310 4.541e-06 1.798e-05
## b[6] 0.0018992 0.0048291 1.080e-05 7.253e-05
## b[7] -0.0530504 0.0482953 1.080e-04 5.479e-04
## b[8] -0.2682581 0.7176049 1.605e-03 1.162e-01
## b[9] 0.0003209 0.0006678 1.493e-06 5.013e-06
## rho12 -0.7335677 0.1014565 2.269e-04 3.761e-04
## tau 6.3785553 0.7039709 1.574e-03 3.112e-03
## u[1,1] 0.1767952 0.7132692 1.595e-03 1.145e-01
## u[2,1] 0.4098966 0.7166177 1.602e-03 1.165e-01
## u[3,1] 0.3492482 0.7159102 1.601e-03 1.194e-01
## u[4,1] 0.3390226 0.7214067 1.613e-03 1.182e-01
## u[5,1] 0.2275859 0.7140559 1.597e-03 1.211e-01
## u[6,1] 0.3144307 0.7154721 1.600e-03 1.214e-01
## u[7,1] 0.3076919 0.7208290 1.612e-03 1.202e-01
## u[8,1] 0.3316518 0.7192204 1.608e-03 1.161e-01
## u[9,1] 0.3668030 0.7183958 1.606e-03 1.160e-01
## u[10,1] 0.2532712 0.7171159 1.604e-03 1.218e-01
## u[11,1] 0.2938173 0.7204021 1.611e-03 1.183e-01
## u[12,1] 0.2567256 0.7206316 1.611e-03 1.219e-01
## u[13,1] 0.0452871 0.7191123 1.608e-03 1.214e-01
## u[14,1] 0.4262527 0.7203427 1.611e-03 1.185e-01
## u[15,1] 0.4142366 0.7175641 1.605e-03 1.149e-01
## u[16,1] 0.3506178 0.7210824 1.612e-03 1.171e-01
## u[17,1] 0.2574542 0.7164915 1.602e-03 1.186e-01
## u[18,1] 0.3593535 0.7175202 1.604e-03 1.146e-01
## u[19,1] 0.2400627 0.7162367 1.602e-03 1.203e-01
## u[20,1] 0.3014717 0.7161569 1.601e-03 1.139e-01
## u[21,1] 0.2751692 0.7162111 1.601e-03 1.174e-01
## u[22,1] 0.3303303 0.7153043 1.599e-03 1.151e-01
## u[23,1] 0.2894581 0.7166352 1.602e-03 1.174e-01
## u[24,1] 0.3467643 0.7165647 1.602e-03 1.188e-01
## u[25,1] 0.2955389 0.7160196 1.601e-03 1.181e-01
## u[26,1] 0.3750559 0.7155612 1.600e-03 1.179e-01
## u[27,1] 0.1953586 0.7139816 1.597e-03 1.151e-01
## u[1,2] -1.0612331 3.5124052 7.854e-03 4.824e-01
## u[2,2] -3.4312486 3.5633245 7.968e-03 4.907e-01
## u[3,2] -2.7415641 3.5432958 7.923e-03 4.996e-01
## u[4,2] -3.2139587 3.6892142 8.249e-03 4.954e-01
## u[5,2] -0.2355212 3.5114290 7.852e-03 5.053e-01
## u[6,2] -1.7490985 3.5127352 7.855e-03 4.930e-01
## u[7,2] -2.2526716 3.6431645 8.146e-03 5.024e-01
## u[8,2] -2.2479455 3.5675986 7.977e-03 5.122e-01
## u[9,2] -3.0344206 3.5504230 7.939e-03 5.106e-01
## u[10,2] -1.4578350 3.5490599 7.936e-03 4.814e-01
## u[11,2] -2.1279031 3.6734031 8.214e-03 4.932e-01
## u[12,2] -1.8161046 3.5935473 8.035e-03 4.936e-01
## u[13,2] 0.8985677 3.5876285 8.022e-03 4.679e-01
## u[14,2] -2.7458182 3.6228836 8.101e-03 4.834e-01
## u[15,2] -2.9376063 3.5562892 7.952e-03 4.675e-01
## u[16,2] -2.4683810 3.6335063 8.125e-03 5.066e-01
## u[17,2] -2.0715560 3.5237555 7.879e-03 5.019e-01
## u[18,2] -2.8554186 3.5580922 7.956e-03 5.033e-01
## u[19,2] -1.1482265 3.4069035 7.618e-03 4.968e-01
## u[20,2] -2.2862453 3.4971791 7.820e-03 5.047e-01
## u[21,2] -2.0379131 3.5199254 7.871e-03 5.123e-01
## u[22,2] -2.4367671 3.5253595 7.883e-03 4.941e-01
## u[23,2] -1.7732653 3.4968612 7.819e-03 5.116e-01
## u[24,2] -2.6133254 3.5493172 7.937e-03 5.044e-01
## u[25,2] -1.4549802 3.6146560 8.083e-03 5.065e-01
## u[26,2] -2.4081407 3.4996987 7.826e-03 5.011e-01
## u[27,2] 0.1839307 3.4157130 7.638e-03 4.925e-01
##
## 2. Quantiles for each variable:
##
## 2.5% 5% 95% 97.5%
## Sigma[1,1] 6.902e-03 0.0075609 0.021411 0.023930
## Sigma[2,1] -1.975e-01 -0.1746043 -0.046473 -0.039924
## Sigma[1,2] -1.975e-01 -0.1746043 -0.046473 -0.039924
## Sigma[2,2] 6.646e-01 0.7359750 2.224191 2.495970
## b[1] -1.066e-03 -0.0007974 0.001962 0.002222
## b[2] -6.453e+00 -5.9187199 0.077959 0.698624
## b[3] -6.813e-02 -0.0592878 0.031919 0.040828
## b[4] -3.778e-02 -0.0159789 0.198787 0.220002
## b[5] -5.269e-03 -0.0046113 0.002064 0.002698
## b[6] -7.829e-03 -0.0061351 0.009670 0.011126
## b[7] -1.477e-01 -0.1324413 0.026477 0.041735
## b[8] -1.327e+00 -1.2149939 1.158751 1.377456
## b[9] -9.863e-04 -0.0007776 0.001420 0.001631
## rho12 -8.847e-01 -0.8685809 -0.544319 -0.492599
## tau 5.076e+00 5.2661417 7.576784 7.834061
## u[1,1] -1.457e+00 -1.2319746 1.121948 1.230695
## u[2,1] -1.228e+00 -1.0161688 1.360617 1.471061
## u[3,1] -1.285e+00 -1.0755821 1.295616 1.407346
## u[4,1] -1.319e+00 -1.1005661 1.291663 1.402093
## u[5,1] -1.406e+00 -1.1920288 1.173190 1.283529
## u[6,1] -1.323e+00 -1.1089594 1.259732 1.372227
## u[7,1] -1.347e+00 -1.1324977 1.259967 1.369962
## u[8,1] -1.314e+00 -1.1029122 1.281773 1.394006
## u[9,1] -1.282e+00 -1.0654829 1.319727 1.427677
## u[10,1] -1.394e+00 -1.1732302 1.203084 1.312171
## u[11,1] -1.364e+00 -1.1446101 1.247166 1.356616
## u[12,1] -1.401e+00 -1.1819715 1.210163 1.319882
## u[13,1] -1.609e+00 -1.3864355 0.998420 1.106120
## u[14,1] -1.231e+00 -1.0140825 1.380866 1.489753
## u[15,1] -1.229e+00 -1.0147199 1.362467 1.475799
## u[16,1] -1.307e+00 -1.0916502 1.302121 1.413118
## u[17,1] -1.388e+00 -1.1656478 1.204790 1.314215
## u[18,1] -1.290e+00 -1.0661373 1.308226 1.418139
## u[19,1] -1.404e+00 -1.1820103 1.187517 1.296480
## u[20,1] -1.342e+00 -1.1169519 1.247341 1.357017
## u[21,1] -1.368e+00 -1.1456876 1.222569 1.331126
## u[22,1] -1.311e+00 -1.0848399 1.275810 1.384522
## u[23,1] -1.358e+00 -1.1324119 1.234824 1.345680
## u[24,1] -1.297e+00 -1.0785440 1.292563 1.403353
## u[25,1] -1.348e+00 -1.1250967 1.240837 1.351757
## u[26,1] -1.270e+00 -1.0460499 1.321235 1.431077
## u[27,1] -1.437e+00 -1.2173766 1.139019 1.250273
## u[1,2] -7.782e+00 -6.6284091 4.989643 6.211596
## u[2,2] -1.023e+01 -9.0317417 2.732629 3.937803
## u[3,2] -9.549e+00 -8.3404721 3.358385 4.530828
## u[4,2] -1.026e+01 -9.0454472 3.135060 4.402625
## u[5,2] -6.995e+00 -5.7868915 5.824094 7.030410
## u[6,2] -8.500e+00 -7.2918839 4.294099 5.484548
## u[7,2] -9.215e+00 -8.0251958 4.030086 5.263622
## u[8,2] -9.110e+00 -7.8925312 3.900965 5.103167
## u[9,2] -9.851e+00 -8.6378403 3.101585 4.349864
## u[10,2] -8.307e+00 -7.0601407 4.658784 5.898439
## u[11,2] -9.221e+00 -7.9214946 4.196840 5.491607
## u[12,2] -8.752e+00 -7.4717419 4.370511 5.648082
## u[13,2] -6.044e+00 -4.7530386 7.072301 8.352060
## u[14,2] -9.719e+00 -8.4689011 3.487261 4.771750
## u[15,2] -9.796e+00 -8.5542036 3.175363 4.424958
## u[16,2] -9.464e+00 -8.1967958 3.779300 5.046218
## u[17,2] -8.831e+00 -7.6272414 4.004332 5.222245
## u[18,2] -9.674e+00 -8.4684517 3.288630 4.496681
## u[19,2] -7.691e+00 -6.5202668 4.741497 5.907003
## u[20,2] -8.989e+00 -7.7829201 3.743758 4.934541
## u[21,2] -8.770e+00 -7.6057873 4.039051 5.256644
## u[22,2] -9.173e+00 -8.0174155 3.638764 4.830611
## u[23,2] -8.485e+00 -7.2931352 4.250986 5.445516
## u[24,2] -9.400e+00 -8.2047126 3.512052 4.744964
## u[25,2] -8.379e+00 -7.1621115 4.777746 6.004306
## u[26,2] -9.092e+00 -7.9302572 3.640119 4.841649
## u[27,2] -6.340e+00 -5.2002297 6.097174 7.306433Segundo a literatura, o valor de R não pode ultrapassar 1.2.
#GBR diagnostic
gelman.diag(samps, multivariate = F)## Potential scale reduction factors:
##
## Point est. Upper C.I.
## Sigma[1,1] 1.00 1.00
## Sigma[2,1] 1.00 1.00
## Sigma[1,2] 1.00 1.00
## Sigma[2,2] 1.00 1.00
## b[1] 1.00 1.00
## b[2] 1.03 1.15
## b[3] 1.01 1.03
## b[4] 1.17 1.58
## b[5] 1.00 1.01
## b[6] 1.01 1.03
## b[7] 1.00 1.02
## b[8] 2.34 5.99
## b[9] 1.00 1.01
## rho12 1.00 1.00
## tau 1.00 1.00
## u[1,1] 2.36 6.04
## u[2,1] 2.35 5.98
## u[3,1] 2.36 6.05
## u[4,1] 2.31 5.91
## u[5,1] 2.36 6.06
## u[6,1] 2.35 6.04
## u[7,1] 2.31 5.91
## u[8,1] 2.33 5.96
## u[9,1] 2.32 5.94
## u[10,1] 2.33 5.97
## u[11,1] 2.31 5.90
## u[12,1] 2.31 5.89
## u[13,1] 2.32 5.92
## u[14,1] 2.31 5.89
## u[15,1] 2.34 6.01
## u[16,1] 2.31 5.92
## u[17,1] 2.34 6.00
## u[18,1] 2.33 5.96
## u[19,1] 2.34 6.02
## u[20,1] 2.34 6.01
## u[21,1] 2.34 5.99
## u[22,1] 2.35 6.02
## u[23,1] 2.34 5.99
## u[24,1] 2.34 6.01
## u[25,1] 2.34 6.02
## u[26,1] 2.34 6.01
## u[27,1] 2.36 6.06
## u[1,2] 1.25 1.80
## u[2,2] 1.24 1.78
## u[3,2] 1.23 1.76
## u[4,2] 1.24 1.77
## u[5,2] 1.24 1.78
## u[6,2] 1.24 1.78
## u[7,2] 1.24 1.78
## u[8,2] 1.24 1.78
## u[9,2] 1.25 1.81
## u[10,2] 1.25 1.81
## u[11,2] 1.25 1.80
## u[12,2] 1.25 1.81
## u[13,2] 1.25 1.80
## u[14,2] 1.24 1.79
## u[15,2] 1.24 1.78
## u[16,2] 1.24 1.79
## u[17,2] 1.24 1.79
## u[18,2] 1.24 1.79
## u[19,2] 1.24 1.79
## u[20,2] 1.24 1.79
## u[21,2] 1.24 1.79
## u[22,2] 1.24 1.78
## u[23,2] 1.24 1.79
## u[24,2] 1.24 1.79
## u[25,2] 1.24 1.78
## u[26,2] 1.24 1.79
## u[27,2] 1.25 1.81gelman.plot(samps)
mcmc_trace(samps)
mcmc_dens(samps)
plot(samps)
#DIC
dic.samples(mod1, n.iter = 1000,type = "pD")## Mean deviance: 222.3
## penalty 56.93
## Penalized deviance: 279.2