## Rows: 6,824
## Columns: 11
## $ ano <dbl> 2014, 2014, 2014, 2014, 2014, 2014, 2014, 2014, 2014, 20…
## $ mun <chr> "Abadia dos Dourados", "Abaeté", "Abre Campo", "Acaiaca"…
## $ cod_municipio <dbl> 3100104, 3100203, 3100302, 3100401, 3100500, 3100609, 31…
## $ hom <dbl> 1, 1, 1, 1, 3, 0, 0, 1, 12, 0, 10, 0, 1, 0, 2, 6, 1, 8, …
## $ cad <dbl> 0.4816, 0.5361, 0.6328, 0.7165, 0.6778, 0.7496, 0.4589, …
## $ gpc <dbl> 1.44, 2.81, 2.45, 7.09, 4.38, 0.78, 13.04, 10.64, 6.02, …
## $ edu <dbl> 0.2398, 0.1849, 0.1155, 0.0000, 0.1154, 0.0000, 0.1265, …
## $ rend <dbl> 1.37168, 1.23592, 1.13913, 1.09220, 1.10698, 1.04886, 1.…
## $ pop <dbl> 6913, 23224, 13551, 4003, 10087, 14662, 2042, 4286, 1903…
## $ pmpc <dbl> 0.95771, 0.56725, 0.60505, 0.78400, 0.93418, 1.51950, 0.…
## $ lpop <dbl> 8.841159, 10.052942, 9.514216, 8.294799, 9.219003, 9.593…
n_ano <- length(unique(data1$ano)) # 2014 a 2021 (8 anos)
n_municipio <- length(unique(data1$cod_municipio)) # 853 municÃpios
observed <- matrix(data1$hom, nrow=n_municipio, ncol=n_ano)
x1 <- matrix(data1$cad, nrow=n_municipio, ncol=n_ano)
x2 <- matrix(data1$edu, nrow=n_municipio, ncol=n_ano)
x3 <- matrix(data1$gpc, nrow=n_municipio, ncol=n_ano)
x4 <- matrix(data1$rend, nrow=n_municipio, ncol=n_ano)
x5 <- matrix(data1$pmpc, nrow=n_municipio, ncol=n_ano)
x6 <- matrix(data1$lpop, nrow=n_municipio, ncol=n_ano)
population <- matrix(data1$pop, nrow=n_municipio, ncol=n_ano)
expected <- matrix(NA, nrow=n_municipio, ncol=n_ano)
for(y in 1:n_ano) {
r <- sum(observed[,y])/sum(population[,y]) # Taxa para o ano y
expected[,y] <- population[,y] * r # Aplicar a todos os municÃpios no ano y
}
# 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_municipio = n_municipio,
n_ano = n_ano,
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 <- 10000 # Total de iterações (incluindo warm-up)
warmup <- 5000
chains <- 1
seed = 12345
stanModel <- rstan::stan_model("script_mod.stan")
output <- sampling(
stanModel,
data = data,
iter = iter,
warmup = warmup,
seed = seed,
chains = chains,
pars = pars,
init = initvals,
verbose = FALSE
)
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 0.006997 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 69.97 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 1: Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 1: Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 1: Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 1: Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 1: Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 1: Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 1: Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 1: Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 1: Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 1: Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 1: Iteration: 10000 / 10000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 1741.09 seconds (Warm-up)
## Chain 1: 1618.6 seconds (Sampling)
## Chain 1: 3359.69 seconds (Total)
## Chain 1:
samp = rstan::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,1])
colnames(aux) = c("beta0", "beta1", "beta2", "beta3", "beta4", "beta5", "beta6", "theta[,1,1]")
aux = as.mcmc(aux); geweke.diag(aux)
##
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5
##
## beta0 beta1 beta2 beta3 beta4 beta5
## 0.4458 0.7910 -1.9830 1.5213 1.5446 -0.7659
## beta6 theta[,1,1]
## -1.3022 -2.0956
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:5000
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 5000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## beta0 -2.416888 0.0772153 1.092e-03 1.092e-03
## beta1 0.920139 0.0522963 7.396e-04 7.396e-04
## beta2 -1.696812 0.0869023 1.229e-03 1.229e-03
## beta3 -0.002547 0.0003475 4.914e-06 4.914e-06
## beta4 0.008857 0.0137584 1.946e-04 1.852e-04
## beta5 -0.112982 0.0225983 3.196e-04 3.245e-04
## beta6 0.196607 0.0058292 8.244e-05 8.244e-05
## theta[,1,1] 0.476185 0.0090505 1.280e-04 1.280e-04
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## beta0 -2.564604 -2.469e+00 -2.417216 -2.363107 -2.265794
## beta1 0.816499 8.845e-01 0.919838 0.955217 1.021678
## beta2 -1.866652 -1.757e+00 -1.695522 -1.637772 -1.529026
## beta3 -0.003248 -2.778e-03 -0.002543 -0.002308 -0.001886
## beta4 -0.019325 -7.651e-05 0.009003 0.018034 0.035044
## beta5 -0.156822 -1.276e-01 -0.113412 -0.098376 -0.067308
## beta6 0.184947 1.927e-01 0.196701 0.200611 0.207848
## theta[,1,1] 0.458786 4.700e-01 0.476089 0.482183 0.494066
##
##
## $hpd
## lower upper
## beta0 -2.558192308 -2.259944374
## beta1 0.818801742 1.023794611
## beta2 -1.867872106 -1.532127123
## beta3 -0.003212337 -0.001852334
## beta4 -0.018286908 0.035816070
## beta5 -0.157450520 -0.068441916
## beta6 0.185067978 0.207926322
## theta[,1,1] 0.458669048 0.493636052
## attr(,"Probability")
## [1] 0.95
##
## $effective_size
## beta0 beta1 beta2 beta3 beta4 beta5
## 5000.000 5000.000 5000.000 5000.000 5520.197 4850.745
## beta6 theta[,1,1]
## 5000.000 5000.000
##
## $geweke
##
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5
##
## beta0 beta1 beta2 beta3 beta4 beta5
## 0.4458 0.7910 -1.9830 1.5213 1.5446 -0.7659
## beta6 theta[,1,1]
## -1.3022 -2.0956
## [1] -2.416888
## [1] 0.9201386
## [1] -1.696812
## [1] -0.00254695
## [1] 0.008856511
## [1] -0.1129817
## [1] 0.1966074
## [1] 0.7033233
saveRDS(rstan::extract(output), "output.rds")
posterior <- as.array(aux)
mcmc_trace(posterior, pars = c("beta0", "beta1", "beta2", "beta3",
"beta4", "beta5", "beta6", "theta[,1,1]"))
