Carregando pacotes
## Warning: package 'mlbench' was built under R version 4.3.3## Warning: package 'bayestestR' was built under R version 4.3.3## Warning: package 'bayesplot' was built under R version 4.3.3## This is bayesplot version 1.11.1## - Online documentation and vignettes at mc-stan.org/bayesplot## - bayesplot theme set to bayesplot::theme_default()## * Does _not_ affect other ggplot2 plots## * See ?bayesplot_theme_set for details on theme setting## Warning: package 'insight' was built under R version 4.3.3## Warning: package 'broom' was built under R version 4.3.3## Warning: package 'rstanarm' was built under R version 4.3.3## Carregando pacotes exigidos: Rcpp## This is rstanarm version 2.32.1## - See https://mc-stan.org/rstanarm/articles/priors for changes to default priors!## - Default priors may change, so it's safest to specify priors, even if equivalent to the defaults.## - For execution on a local, multicore CPU with excess RAM we recommend calling## options(mc.cores = parallel::detectCores())## Warning: package 'ggplot2' was built under R version 4.3.3## 'data.frame': 506 obs. of 14 variables:
## $ crim : num 0.00632 0.02731 0.02729 0.03237 0.06905 ...
## $ zn : num 18 0 0 0 0 0 12.5 12.5 12.5 12.5 ...
## $ indus : num 2.31 7.07 7.07 2.18 2.18 2.18 7.87 7.87 7.87 7.87 ...
## $ chas : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
## $ nox : num 0.538 0.469 0.469 0.458 0.458 0.458 0.524 0.524 0.524 0.524 ...
## $ rm : num 6.58 6.42 7.18 7 7.15 ...
## $ age : num 65.2 78.9 61.1 45.8 54.2 58.7 66.6 96.1 100 85.9 ...
## $ dis : num 4.09 4.97 4.97 6.06 6.06 ...
## $ rad : num 1 2 2 3 3 3 5 5 5 5 ...
## $ tax : num 296 242 242 222 222 222 311 311 311 311 ...
## $ ptratio: num 15.3 17.8 17.8 18.7 18.7 18.7 15.2 15.2 15.2 15.2 ...
## $ b : num 397 397 393 395 397 ...
## $ lstat : num 4.98 9.14 4.03 2.94 5.33 ...
## $ medv : num 24 21.6 34.7 33.4 36.2 28.7 22.9 27.1 16.5 18.9 ...Selecionar novas variáveis para análise
## medv age dis nox rm
## 1 24.0 65.2 4.0900 0.538 6.575
## 2 21.6 78.9 4.9671 0.469 6.421
## 3 34.7 61.1 4.9671 0.469 7.185
## 4 33.4 45.8 6.0622 0.458 6.998
## 5 36.2 54.2 6.0622 0.458 7.147
## 6 28.7 58.7 6.0622 0.458 6.430## 'data.frame': 506 obs. of 5 variables:
## $ medv: num 24 21.6 34.7 33.4 36.2 28.7 22.9 27.1 16.5 18.9 ...
## $ age : num 65.2 78.9 61.1 45.8 54.2 58.7 66.6 96.1 100 85.9 ...
## $ dis : num 4.09 4.97 4.97 6.06 6.06 ...
## $ nox : num 0.538 0.469 0.469 0.458 0.458 0.458 0.524 0.524 0.524 0.524 ...
## $ rm : num 6.58 6.42 7.18 7 7.15 ...modelo com novas variáveis
##
## Call:
## lm(formula = medv ~ age + dis + nox + rm, data = bost)
##
## Residuals:
## Min 1Q Median 3Q Max
## -18.045 -2.961 -0.558 2.030 39.178
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6.61135 4.15733 -1.590 0.112
## age -0.06933 0.01568 -4.422 1.20e-05 ***
## dis -1.08527 0.22376 -4.850 1.65e-06 ***
## nox -22.10858 4.02962 -5.487 6.52e-08 ***
## rm 8.00052 0.40705 19.655 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.112 on 501 degrees of freedom
## Multiple R-squared: 0.5619, Adjusted R-squared: 0.5584
## F-statistic: 160.7 on 4 and 501 DF, p-value: < 2.2e-16Regressão Bayesiana Ajustamos um modelo bayesiano para comparação com o modelo de regressão tradicional. Usamos stan_glm()para especificar um modelo bayesiano de regressão com variável dependente medve as mesmas variáveis preditoras.
modelo_bayes <- stan_glm(medv ~ age + dis + nox + rm,
data = bost,
prior = normal(0, 1),
prior_intercept = normal(0, 5),
seed = 111) ##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 0.000127 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.27 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.163 seconds (Warm-up)
## Chain 1: 0.115 seconds (Sampling)
## Chain 1: 0.278 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 1.9e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.19 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.152 seconds (Warm-up)
## Chain 2: 0.117 seconds (Sampling)
## Chain 2: 0.269 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 1.9e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.19 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.141 seconds (Warm-up)
## Chain 3: 0.145 seconds (Sampling)
## Chain 3: 0.286 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 1.7e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.17 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.156 seconds (Warm-up)
## Chain 4: 0.176 seconds (Sampling)
## Chain 4: 0.332 seconds (Total)
## Chain 4:## stan_glm
## family: gaussian [identity]
## formula: medv ~ age + dis + nox + rm
## observations: 506
## predictors: 5
## ------
## Median MAD_SD
## (Intercept) -13.260 3.099
## age -0.102 0.016
## dis -0.475 0.198
## nox -1.331 0.977
## rm 7.200 0.385
##
## Auxiliary parameter(s):
## Median MAD_SD
## sigma 6.321 0.209
##
## ------
## * For help interpreting the printed output see ?print.stanreg
## * For info on the priors used see ?prior_summary.stanregmcmc_dens(modelo_bayes, pars = c("age", "dis", "nox", "rm")) +
geom_vline(xintercept = 0, color = "blue", linetype = "dashed")
Extraímos as configurações pós-ajuste ( a posteriori) do modelo Bayesiano.
Calculamos as médias, medianas e MAP (estimativas de máxima verossimilhança a posteriori) das parâmetros, valores centrais para descrever as parâmetros a posteriori.
mean_vals <- purrr::map_dbl(post, mean)
median_vals <- purrr::map_dbl(post, median)
map_estimates <- purrr::map_dbl(post, ~ as.numeric(map_estimate(.x)))Visualizamos a densidade do parâmetro age, com linhas verticais decrescentes a média (amarelo), a mediana (vermelho) e o MAP (verde), facilitando a interpretação da distribuição a posteriori.
mcmc_dens(modelo_bayes, pars = c("age")) +
geom_vline(xintercept = median_vals["age"], color = "red") +
geom_vline(xintercept = mean_vals["age"], color = "yellow") +
geom_vline(xintercept = map_estimates["age"], color = "green")
O gráfico mostra a distribuição do parâmetro ageno modelo Bayesiano, com
três linhas verticais: a mediana (vermelha), a média (amarela) e o valor
mais provável (MAP, em verde). A linha azul tracejada no zero serve de
referência para ver o efeito de ageincluir zero. Como a distribuição
está toda à esquerda de zero, isso indica que agehá um efeito negativo
significativo sobre medvnenhum modelo.
Intervalo de Credibilidade e ROPE Calculamos o intervalo de alta densidade (HDI) e o intervalo de substituição eti (Equal-tailed Interval), usado para identificar o intervalo de valores mais prováveis dos parâmetros.
Calculamos o ROPE (Região de Equivalência Prática) para avaliar se as parâmetros estão próximas de zero, o que ajuda a determinar a significância prática dos coeficientes a posteriori.
rope_age <- rope(post$age)
rope_nox <- rope(post$nox)
rope_intercept <- rope(post$`(Intercept)`)
rope_range_total <- rope_range(modelo_bayes)Exibimos os valores das variações de alteração calculada (CrI) e do HDI para uma visão geral da incerteza em torno dos parâmetros estimados.
## [1] "Intervalo de Credibilidade (CrI) e HDI:"## Highest Density Interval
##
## Parameter | 95% HDI
## -----------------------------
## (Intercept) | [-19.45, -7.49]
## age | [ -0.13, -0.07]
## dis | [ -0.87, -0.09]
## nox | [ -3.10, 0.74]
## rm | [ 6.46, 7.97]## Equal-Tailed Interval
##
## Parameter | 95% ETI
## -----------------------------
## (Intercept) | [-19.17, -7.16]
## age | [ -0.13, -0.07]
## dis | [ -0.87, -0.09]
## nox | [ -3.25, 0.62]
## rm | [ 6.44, 7.95]esses resultados mostram as diferenças de variação para as parâmetros de um modelo Bayesiano ajustado, que representam uma incerteza em torno dos coeficientes estimados. Dois tipos de intervalos são mostrados: o Intervalo de Maior Densidade (HDI) e o Intervalo Equal-Tailed (ETI) , ambos no nível de 95%.
O IDH indica o intervalo mais provável onde os valores dos parâmetros estão concentrados, enquanto o IDH considera 2,5% de cada cauda da distribuição a posteriori, sendo simétrico em relação às caudas e menos sensível à forma da distribuição.
Os intervalos dos intervalos agee disnão rmincluem zero, indicando que esses preditores têm efeito significativo sobre medv. Já o parâmetro noxpossui um intervalo que cruza zero, traz incerteza sobre seu efeito real em medv.
Concluimos
As áreas como medicina adotaram a análise bayesiana para permitir interpretações probabilísticas diretas e atualizáveis com novas evidências. Contudo, esta abordagem traz desafios, como a subjetividade na definição das distribuições a priori, que afeta os resultados, e a dificuldade de convergência dos algoritmos MCMC em problemas complexos, especialmente na ausência de a priori conjugada.