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
stanModel <- rstan::stan_model("script_mod20-21.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.001753 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 17.53 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: 345.389 seconds (Warm-up)
## Chain 1: 427.534 seconds (Sampling)
## Chain 1: 772.923 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.48658 -0.52534 0.28168 0.22009 -1.36840 -0.41532
## beta6 theta[,1,1]
## 0.02744 0.38893
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: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.053468 0.23215 3.283e-03 3.283e-03
## beta1 0.996602 0.13663 1.932e-03 1.932e-03
## beta2 -0.715563 0.19676 2.783e-03 2.612e-03
## beta3 -0.001526 0.00105 1.485e-05 1.485e-05
## beta4 -0.025432 0.02926 4.138e-04 4.138e-04
## beta5 0.126811 0.11414 1.614e-03 1.614e-03
## beta6 0.154766 0.01659 2.346e-04 2.346e-04
## theta[,1,1] 0.597663 0.02880 4.074e-04 4.074e-04
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## beta0 -2.509508 -2.207111 -2.056641 -1.898019 -1.5887939
## beta1 0.732816 0.903740 0.994998 1.087726 1.2667077
## beta2 -1.115439 -0.847945 -0.716334 -0.580783 -0.3272180
## beta3 -0.003566 -0.002229 -0.001517 -0.000824 0.0005614
## beta4 -0.084236 -0.044861 -0.025481 -0.006323 0.0327793
## beta5 -0.101093 0.050172 0.127015 0.202315 0.3507811
## beta6 0.122536 0.143583 0.155361 0.165942 0.1871160
## theta[,1,1] 0.542775 0.577481 0.597370 0.616544 0.6559154
##
##
## $hpd
## lower upper
## beta0 -2.517009274 -1.5990365610
## beta1 0.725946069 1.2573081791
## beta2 -1.097567235 -0.3122274182
## beta3 -0.003513325 0.0006020531
## beta4 -0.084621788 0.0317487819
## beta5 -0.096886363 0.3525467249
## beta6 0.121728804 0.1861774284
## theta[,1,1] 0.541195726 0.6537816823
## attr(,"Probability")
## [1] 0.95
##
## $effective_size
## beta0 beta1 beta2 beta3 beta4 beta5
## 5000.000 5000.000 5674.654 5000.000 5000.000 5000.000
## 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.48658 -0.52534 0.28168 0.22009 -1.36840 -0.41532
## beta6 theta[,1,1]
## 0.02744 0.38893