## mun cod_municipio latitude longitude
## Length:853 Min. :3100104 Min. :-22.85 Min. :-50.69
## Class :character 1st Qu.:3119500 1st Qu.:-21.09 1st Qu.:-45.46
## Mode :character Median :3137304 Median :-19.88 Median :-43.93
## Mean :3136982 Mean :-19.58 Mean :-44.10
## 3rd Qu.:3154903 3rd Qu.:-18.48 3rd Qu.:-42.62
## Max. :3172202 Max. :-14.27 Max. :-39.94
## H21 cad21 edu21 gpc21
## Min. : 0.000 Min. :0.1549 Min. :0.000 Min. : 0.000
## 1st Qu.: 0.000 1st Qu.:0.3843 1st Qu.:0.027 1st Qu.: 0.000
## Median : 1.000 Median :0.4832 Median :0.104 Median : 6.080
## Mean : 3.002 Mean :0.4981 Mean :0.118 Mean : 8.701
## 3rd Qu.: 2.000 3rd Qu.:0.6171 3rd Qu.:0.188 3rd Qu.: 11.280
## Max. :354.000 Max. :0.9689 Max. :0.536 Max. :153.920
## rend21 pmpc21 pop21 lpop
## Min. :1.345 Min. :0.0771 Min. : 771 Min. : 6.648
## 1st Qu.:1.761 1st Qu.:0.4094 1st Qu.: 4843 1st Qu.: 8.485
## Median :1.951 Median :0.5586 Median : 8288 Median : 9.023
## Mean :2.053 Mean :0.5496 Mean : 25102 Mean : 9.232
## 3rd Qu.:2.209 3rd Qu.:0.6782 3rd Qu.: 17850 3rd Qu.: 9.790
## Max. :5.419 Max. :0.9523 Max. :2530701 Max. :14.744
observed = data$H21
x1 = data$cad21
x2 = data$edu21
x3 = data$gpc21
x4 = data$rend21
x5 = data$pmpc21
x6 = data$lpop
population = data$pop21
n=length(observed)
r = sum(observed)/sum(population)
expected = population * r
# Especificações a priori
m0 <- 0
v0 <- 10
m1 <- 0
v1 <- 10
m2 <- 0
v2 <- 10
m3 <- 0
v3 <- 10
m4 <- 0
v4 <- 10
m5 <- 0
v5 <- 10
m6 <- 0
v6 <- 10
# Organizando os elementos para o Stan
data <- list(
n = n,
observed = observed,
expected = expected,
x1 = x1,
x2 = x2,
x3 = x3,
x4 = x4,
x5 = x5,
x6 = x6,
m0 = m0,
v0 = v0,
m1 = m1,
v1 = v1,
m2 = m2,
v2 = v2,
m3 = m3,
v3 = v3,
m4 = m4,
v4 = v4,
m5 = m5,
v5 = v5,
m6 = m6,
v6 = v6
)
# Nomes dos parâmetros a serem estimados
pars <- c("beta0", "beta1", "beta2", "beta3", "beta4", "beta5", "beta6", "theta")
# Chutes iniciais
initvals <- list()
initvals[[1]] <- list(
beta0 = 0,
beta1 = 0,
beta2 = 0,
beta3 = 0,
beta4 = 0,
beta5 = 0,
beta6 = 0
)
# Configurações para o MCMC
iter <- 15000 # Total de iterações (incluindo warm-up)
warmup <- 5000
chains <- 1
stanModel <- rstan::stan_model("script_mod1.stan")
output <- sampling(
stanModel,
data = data,
iter = iter,
warmup = warmup,
chains = chains,
pars = pars,
init = initvals,
verbose = FALSE
)
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 0.001636 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 16.36 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 15000 [ 0%] (Warmup)
## Chain 1: Iteration: 1500 / 15000 [ 10%] (Warmup)
## Chain 1: Iteration: 3000 / 15000 [ 20%] (Warmup)
## Chain 1: Iteration: 4500 / 15000 [ 30%] (Warmup)
## Chain 1: Iteration: 5001 / 15000 [ 33%] (Sampling)
## Chain 1: Iteration: 6500 / 15000 [ 43%] (Sampling)
## Chain 1: Iteration: 8000 / 15000 [ 53%] (Sampling)
## Chain 1: Iteration: 9500 / 15000 [ 63%] (Sampling)
## Chain 1: Iteration: 11000 / 15000 [ 73%] (Sampling)
## Chain 1: Iteration: 12500 / 15000 [ 83%] (Sampling)
## Chain 1: Iteration: 14000 / 15000 [ 93%] (Sampling)
## Chain 1: Iteration: 15000 / 15000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 268.335 seconds (Warm-up)
## Chain 1: 557.759 seconds (Sampling)
## Chain 1: 826.094 seconds (Total)
## Chain 1:
samp = extract(output)
sbeta0 = samp$beta0
sbeta1 = samp$beta1
sbeta2 = samp$beta2
sbeta3 = samp$beta3
sbeta4 = samp$beta4
sbeta5 = samp$beta5
sbeta6 = samp$beta6
stheta = samp$theta
aux = cbind(sbeta0, sbeta1, sbeta2, sbeta3, sbeta4, sbeta5, sbeta6, stheta[,1])
colnames(aux) = c("beta0", "beta1", "beta2", "beta3", "beta4", "beta5", "beta6", "theta[1]")
aux = as.mcmc(aux)
#geweke.diag(aux)
posterior <- as.array(aux)
mcmc_trace(posterior, pars = c("beta0", "beta1", "beta2", "beta3",
"beta4", "beta5", "beta6", "theta[1]"))
traceplot(aux) # plotes um a um
diagnostics <- list(
summary = summary(aux),
hpd = HPDinterval(aux),
effective_size = effectiveSize(aux),
# gelman = gelman.diag(aux) #duas cadeias
geweke = geweke.diag(aux)
)
diagnostics
## $summary
##
## Iterations = 1:10000
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 10000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## beta0 -0.996394 0.256726 2.567e-03 2.567e-03
## beta1 -0.004545 0.001702 1.702e-05 1.702e-05
## beta2 -0.781942 0.258026 2.580e-03 2.502e-03
## beta3 0.037228 3.116740 3.117e-02 3.066e-02
## beta4 -0.021658 0.041957 4.196e-04 4.196e-04
## beta5 0.256747 0.164292 1.643e-03 1.643e-03
## beta6 0.090625 0.020653 2.065e-04 1.997e-04
## theta[1] 0.766958 0.046821 4.682e-04 4.526e-04
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## beta0 -1.499447 -1.169600 -0.992760 -0.822931 -0.492423
## beta1 -0.007988 -0.005659 -0.004524 -0.003387 -0.001247
## beta2 -1.287019 -0.955866 -0.784165 -0.605884 -0.277785
## beta3 -6.116946 -2.046722 0.033626 2.120816 6.167277
## beta4 -0.102668 -0.050191 -0.021685 0.006279 0.060420
## beta5 -0.069102 0.145740 0.256135 0.367465 0.577908
## beta6 0.050068 0.076358 0.090578 0.104668 0.131049
## theta[1] 0.679848 0.733982 0.765098 0.798025 0.862091
##
##
## $hpd
## lower upper
## beta0 -1.476518896 -0.473159320
## beta1 -0.007827059 -0.001107242
## beta2 -1.288242607 -0.279497371
## beta3 -5.942045228 6.267224087
## beta4 -0.104850306 0.057938483
## beta5 -0.065905668 0.578883843
## beta6 0.051188314 0.131792458
## theta[1] 0.681755426 0.863518455
## attr(,"Probability")
## [1] 0.95
##
## $effective_size
## beta0 beta1 beta2 beta3 beta4 beta5 beta6 theta[1]
## 10000.00 10000.00 10637.94 10337.03 10000.00 10000.00 10695.86 10700.40
##
## $geweke
##
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5
##
## beta0 beta1 beta2 beta3 beta4 beta5 beta6 theta[1]
## -1.1869 0.8827 0.1263 1.5733 0.3613 2.3030 0.3930 1.1656