data <- read_csv("C:/Users/thais/Desktop/Monografia/dados2021.csv")
## Rows: 853 Columns: 12
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): mun
## dbl (11): cod_municipio, latitude, longitude, H21, cad21, edu21, gpc21, rend...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
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.001064 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 10.64 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: 138.493 seconds (Warm-up)
## Chain 1: 339.709 seconds (Sampling)
## Chain 1: 478.202 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)
##
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5
##
## beta0 beta1 beta2 beta3 beta4 beta5 beta6 theta[1]
## -2.0194 0.1255 -0.2291 -1.1688 0.5188 0.6574 1.6992 -1.4593
posterior <- as.array(aux)
mcmc_trace(posterior, pars = c("beta0", "beta1", "beta2", "beta3",
"beta4", "beta5", "beta6", "theta[1]"))
