Aula de Modelos Lineares I
Patrícia Lusié Velozo da Costa
2025-06-30
Modelo de Regressão Linear Múltipla - Ajustando MRLMs sob a abordagem bayesiana e selecionando o melhor modelo
1 Simulando um conjunto de dados
Simule um conjunto de dados considerando n=1000 (tamanho amostral), \(\beta_0=2\), \(\beta_1=5\), \(\beta_2=3\), \(\sigma^2=1\) e que \(X_{i1} \stackrel{iid}{\sim} N(0,1)\), \(X_{i2} \stackrel{iid}{\sim} t_{(5)}\) e o seguinte modelo \[Y_i \stackrel{indep}{\sim} N(\beta_0 + \beta_1 X_{i1} + \beta_2 X_{i2}, \sigma^2).\]
# Configurações iniciais
n = 1000 # Número de observações
# Simulando as covariáveis
x1 = rnorm(n, mean = 0, sd = sqrt(1)) # Covariável 1
x2 = rt(n, 5) # Covariável 2
x = cbind(1, x1=x1, x2=x2)
#definindo total de covariaveis
p = ncol(x)-1
# Definindo os coeficientes reais
beta_0 = 2 # Intercepto
beta_1 = 5 # Coeficiente para x1
beta_2 = 3 # Coeficiente para x2
beta = c(beta_0, beta_1, beta_2)
# Gerando a variável resposta com erro aleatório
sigma2 = 1
e = rnorm(n, mean = 0, sd = sqrt(sigma2)) # Erro aleatório
y = x %*% beta + e
# Criando um data frame com os dados
dados = data.frame(y, x1=x1, x2=x2)2 Ajustes usando rstanarm
Considere agora que os parâmetros são desconhecidos e independentes. Supondo que \(\beta_0 \sim N(0,1000)\), \(\beta_1 \sim N(0,1000)\), \(\beta_2 \sim N(0,1000)\) e \(\sigma^2 \sim Exp(1)\), ajuste os seguintes modelos usando o pacote rstanarm:
- \(Y_i \stackrel{indep}{\sim} N(\beta_0 + \beta_1 X_{i1} , \sigma^2).\)
- \(Y_i \stackrel{indep}{\sim} N(\beta_0 + \beta_2 X_{i2}, \sigma^2).\)
- \(Y_i \stackrel{indep}{\sim} N(\beta_0 + \beta_1 X_{i1} + \beta_2 X_{i2}, \sigma^2).\)
if(!require("rstanarm")) {install.packages("rstanarm")}
if(!require("parallel")) {install.packages("parallel")}
library(rstanarm)
library(parallel)
# Marcando o tempo inicial
tempo_inicial <- Sys.time()
# Fazendo os ajustes
ajusteBay_Stanarm_1 = stan_glm(y ~ x1 ,
data = dados,
family = gaussian(),
warmup = 500, chain = 2, iter = 2500,
prior = normal(0, 1000, autoscale = FALSE),
prior_intercept = normal(0, 1000, autoscale = FALSE),
prior_aux = exponential(1) ,
diagnostic_file = file.path(tempdir(), "diag1.csv"),
refresh = 0 # Suprime a saída de progresso
)
# Marcando o tempo final
tempo_final <- Sys.time()
# Calculando o tempo gasto
tempo_final - tempo_inicial## Time difference of 5.177693 secs# Marcando o tempo inicial
tempo_inicial <- Sys.time()
# Fazendo os ajustes
ajusteBay_Stanarm_2 = stan_glm(y ~ x2,
data = dados,
family = gaussian(),
warmup = 500, chain = 2, iter = 2500,
prior = normal(0, 1000, autoscale = FALSE),
prior_intercept = normal(0, 1000, autoscale = FALSE),
prior_aux = exponential(1) ,
diagnostic_file = file.path(tempdir(), "diag2.csv"),
refresh = 0 # Suprime a saída de progresso
)
# Marcando o tempo final
tempo_final <- Sys.time()
# Calculando o tempo gasto
tempo_final - tempo_inicial## Time difference of 1.742544 secs# Marcando o tempo inicial
tempo_inicial <- Sys.time()
# Fazendo os ajustes
ajusteBay_Stanarm_3 = stan_glm(y ~ x1 + x2,
data = dados,
family = gaussian(),
warmup = 500, chain = 2, iter = 2500,
prior = normal(0, 1000, autoscale = FALSE),
prior_intercept = normal(0, 1000, autoscale = FALSE),
prior_aux = exponential(1) ,
diagnostic_file = file.path(tempdir(), "diag3.csv"),
refresh = 0 # Suprime a saída de progresso
)
# Marcando o tempo final
tempo_final <- Sys.time()
# Calculando o tempo gasto
tempo_final - tempo_inicial## Time difference of 4.296973 secs2.1 Analisando os ajustes
# Resumo geral do modelo
summary(ajusteBay_Stanarm_1)##
## Model Info:
## function: stan_glm
## family: gaussian [identity]
## formula: y ~ x1
## algorithm: sampling
## sample: 4000 (posterior sample size)
## priors: see help('prior_summary')
## observations: 1000
## predictors: 2
##
## Estimates:
## mean sd 10% 50% 90%
## (Intercept) 2.0 0.3 1.8 2.0 2.2
## x1 5.0 0.1 4.9 5.0 5.2
## sigma 4.1 0.1 4.0 4.1 4.2
##
## Fit Diagnostics:
## mean sd 10% 50% 90%
## mean_PPD 2.0 0.3 1.8 2.0 2.3
##
## The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
##
## MCMC diagnostics
## mcse Rhat n_eff
## (Intercept) 0.0 1.0 176
## x1 0.0 1.0 3804
## sigma 0.0 1.0 239
## mean_PPD 0.0 1.0 212
## log-posterior 1.9 1.0 52
##
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).summary(ajusteBay_Stanarm_2)##
## Model Info:
## function: stan_glm
## family: gaussian [identity]
## formula: y ~ x2
## algorithm: sampling
## sample: 4000 (posterior sample size)
## priors: see help('prior_summary')
## observations: 1000
## predictors: 2
##
## Estimates:
## mean sd 10% 50% 90%
## (Intercept) 2.0 0.2 1.8 2.0 2.2
## x2 3.1 0.1 2.9 3.1 3.2
## sigma 5.0 0.1 4.9 5.0 5.2
##
## Fit Diagnostics:
## mean sd 10% 50% 90%
## mean_PPD 2.0 0.2 1.7 2.0 2.3
##
## The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
##
## MCMC diagnostics
## mcse Rhat n_eff
## (Intercept) 0.0 1.0 1295
## x2 0.0 1.0 4191
## sigma 0.0 1.0 1773
## mean_PPD 0.0 1.0 2112
## log-posterior 0.0 1.0 1030
##
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).summary(ajusteBay_Stanarm_3)##
## Model Info:
## function: stan_glm
## family: gaussian [identity]
## formula: y ~ x1 + x2
## algorithm: sampling
## sample: 4000 (posterior sample size)
## priors: see help('prior_summary')
## observations: 1000
## predictors: 3
##
## Estimates:
## mean sd 10% 50% 90%
## (Intercept) 2.0 0.0 2.0 2.0 2.0
## x1 5.0 0.0 4.9 5.0 5.0
## x2 3.0 0.0 3.0 3.0 3.0
## sigma 1.0 0.0 1.0 1.0 1.0
##
## Fit Diagnostics:
## mean sd 10% 50% 90%
## mean_PPD 2.0 0.0 1.9 2.0 2.1
##
## The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
##
## MCMC diagnostics
## mcse Rhat n_eff
## (Intercept) 0.0 1.0 1203
## x1 0.0 1.0 4087
## x2 0.0 1.0 3892
## sigma 0.0 1.0 976
## mean_PPD 0.0 1.0 1940
## log-posterior 0.0 1.0 1048
##
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).# Intervalos de credibilidade (posterior)
cbind(posterior_interval(ajusteBay_Stanarm_1, prob = 0.95),c(beta[1], beta[2], sigma2))## 2.5% 97.5%
## (Intercept) 1.713424 2.264000 2
## x1 4.785733 5.299947 5
## sigma 3.927059 4.316453 1cbind(posterior_interval(ajusteBay_Stanarm_2, prob = 0.95),c(beta[1], beta[3], sigma2))## 2.5% 97.5%
## (Intercept) 1.726394 2.360089 2
## x2 2.859945 3.323756 3
## sigma 4.838432 5.271319 1cbind(posterior_interval(ajusteBay_Stanarm_3, prob = 0.95),c(beta, sigma2))## 2.5% 97.5%
## (Intercept) 1.9407681 2.070654 2
## x1 4.8997797 5.025167 5
## x2 2.9672695 3.059690 3
## sigma 0.9701407 1.057230 1# Analise grafica
plot(ajusteBay_Stanarm_1, prob = 0.95)
plot(ajusteBay_Stanarm_2, prob = 0.95)
plot(ajusteBay_Stanarm_3, prob = 0.95)
library(bayesplot)
amostras_1 = as.matrix(ajusteBay_Stanarm_1)
amostras_2 = as.matrix(ajusteBay_Stanarm_2)
amostras_3 = as.matrix(ajusteBay_Stanarm_3)
# Traços para todos os parâmetros
mcmc_trace(amostras_1, pars = c("(Intercept)", "x1", "sigma"))
mcmc_trace(amostras_2, pars = c("(Intercept)", "x2", "sigma"))
mcmc_trace(amostras_3, pars = c("(Intercept)", "x1", "x2", "sigma"))
# Autocorrelação por parâmetro
mcmc_acf(amostras_1, pars = c("(Intercept)", "x1", "sigma"))
mcmc_acf(amostras_2, pars = c("(Intercept)", "x2", "sigma"))
mcmc_acf(amostras_3, pars = c("(Intercept)", "x1", "x2", "sigma"))
# Fazer thinning manual (ex: a cada 5 iterações)
amostras_1_thin = amostras_1[seq(1, nrow(amostras_1), by = 5), ]
mcmc_acf(amostras_1_thin, pars = c("(Intercept)", "x1", "sigma"))
amostras_2_thin = amostras_2[seq(1, nrow(amostras_2), by = 5), ]
mcmc_acf(amostras_2_thin, pars = c("(Intercept)", "x2", "sigma"))
amostras_3_thin = amostras_3[seq(1, nrow(amostras_3), by = 5), ]
mcmc_acf(amostras_3_thin, pars = c("(Intercept)", "x1", "x2", "sigma"))
amostras_1_thin = amostras_1
amostras_2_thin = amostras_2
amostras_3_thin = amostras_3
# Analise dos residuos
residuos_1 = residuals(ajusteBay_Stanarm_1)
valores_ajustados_1 = fitted(ajusteBay_Stanarm_1)
residuos_2 = residuals(ajusteBay_Stanarm_2)
valores_ajustados_2 = fitted(ajusteBay_Stanarm_2)
residuos_3 = residuals(ajusteBay_Stanarm_3)
valores_ajustados_3 = fitted(ajusteBay_Stanarm_3)
par(mfrow=c(3,3))
plot(valores_ajustados_1, residuos_1,
xlab = "Valores Ajustados", ylab = "Resíduos",
main = "Resíduos vs Ajustados", bty="n")
abline(h = 0, col = "red", lty = 2)
hist(residuos_1, breaks = 30, main = "Histograma dos Resíduos",
xlab = "Resíduos", col = "lightblue", freq=F, ylab="densidade")
qqnorm(residuos_1)
qqline(residuos_1, col = "red", lwd = 2)
plot(valores_ajustados_2, residuos_2,
xlab = "Valores Ajustados", ylab = "Resíduos",
main = "Resíduos vs Ajustados", bty="n")
abline(h = 0, col = "red", lty = 2)
hist(residuos_2, breaks = 30, main = "Histograma dos Resíduos",
xlab = "Resíduos", col = "lightblue", freq=F, ylab="densidade")
qqnorm(residuos_2)
qqline(residuos_2, col = "red", lwd = 2)
plot(valores_ajustados_3, residuos_3,
xlab = "Valores Ajustados", ylab = "Resíduos",
main = "Resíduos vs Ajustados", bty="n")
abline(h = 0, col = "red", lty = 2)
hist(residuos_3, breaks = 30, main = "Histograma dos Resíduos",
xlab = "Resíduos", col = "lightblue", freq=F, ylab="densidade")
qqnorm(residuos_3)
qqline(residuos_3, col = "red", lwd = 2)
3 Selecionando o melhor modelo
3.1 Validação Cruzada Leave-One-Out (LOO-CV)
A validação cruzada Leave-One-Out (LOO-CV) é uma técnica de avaliação de modelos estatísticos e preditivos em que, dado um conjunto de dados com \(n\) observações, o modelo é ajustado \(n\) vezes, cada vez deixando de fora uma observação distinta.
Seja \(\mathcal{D} = \{(x_i, y_i)\}_{i=1}^n\) o conjunto de dados. Em cada iteração \(i\):
A expressão algébrica geral do erro de LOO-CV é:
\[\text{LOO-CV} = \frac{1}{n} \sum_{i=1}^n L\left(y_i, \hat{f}^{(-i)}(x_i)\right)\] onde:
- \(\hat{f}^{(-i)}\) é a função ajustada sem a observação \(i\);
- \(L(\cdot,\cdot)\) é a função de perda (por exemplo, erro quadrático: \(L(y, \hat{y}) = (y - \hat{y})^2\));
- \(n\) é o número total de observações.
A principal vantagem do LOO-CV é que ele utiliza quase toda a base de dados para cada ajuste, produzindo estimativas menos enviesadas do erro de predição. No entanto, ele pode ser computacionalmente custoso para grandes conjuntos de dados, especialmente em modelos complexos.
Em métodos Bayesianos, técnicas como Pareto-smoothed importance sampling (PSIS) permitem uma aproximação eficiente do LOO-CV sem a necessidade de refazer o ajuste \(n\) vezes.
O critério central é o valor do (ELPD), que estima a soma da log-verossimilhança preditiva ponto a ponto: \[\widehat{\text{ELPD}}_{\text{LOO}} = \sum_{i=1}^n \log p(y_i \mid y_{-i})\] onde \(p(y_i \mid y_{-i})\) representa a densidade preditiva posterior para \(y_i\), condicionada a todos os outros dados exceto \(y_i\).
Sejam dois modelos \(M_1\) e \(M_2\), com ELPDs estimadas \(\widehat{\text{ELPD}}_1\) e \(\widehat{\text{ELPD}}_2\), e erro padrão da diferença \(\text{SE}(\Delta)\), a diferença é dada por: \[\Delta = \widehat{\text{ELPD}}_1 - \widehat{\text{ELPD}}_2.\] Se \(\Delta \gg \text{SE}(\Delta)\), há evidência favorável ao modelo com maior \(\widehat{\text{ELPD}}\).
- Diagnóstico de Pareto \(k\): Para que a aproximação LOO seja confiável, recomenda-se que os valores de \(k\) estejam abaixo de 0,7.
- Parcimônia: Em caso de modelos com ELPD semelhantes, deve-se preferir o modelo mais simples (com menos parâmetros).
- Visualização: Diferenças podem ser exploradas também com gráficos de densidade preditiva e resíduos preditivos.
library(loo)
loo_resultado_1 = loo(ajusteBay_Stanarm_1)
print(loo_resultado_1)##
## Computed from 4000 by 1000 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -2836.9 33.8
## p_loo 7.8 1.2
## looic 5673.8 67.6
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.0, 0.2]).
##
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.# Verificar se houve problemas na aproximação
pareto_k_ids_1 = which(loo_resultado_1$diagnostics$pareto_k > 0.7)
if (length(pareto_k_ids_1) == 0) {
print("A aproximação do LOO-CV foi considerada confiável (k < 0.7 para todas as observações).")
} else {
print(paste("Algumas observações apresentaram k > 0.7. Verifique os índices: ", pareto_k_ids, collapse = ", "))
}## [1] "A aproximação do LOO-CV foi considerada confiável (k < 0.7 para todas as observações)."loo_resultado_2 = loo(ajusteBay_Stanarm_2)
print(loo_resultado_2)##
## Computed from 4000 by 1000 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -3041.1 22.1
## p_loo 2.9 0.2
## looic 6082.2 44.3
## ------
## MCSE of elpd_loo is 0.0.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.3, 1.0]).
##
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.# Verificar se houve problemas na aproximação
pareto_k_ids_2 = which(loo_resultado_2$diagnostics$pareto_k > 0.7)
if (length(pareto_k_ids_2) == 0) {
print("A aproximação do LOO-CV foi considerada confiável (k < 0.7 para todas as observações).")
} else {
print(paste("Algumas observações apresentaram k > 0.7. Verifique os índices: ", pareto_k_ids, collapse = ", "))
}## [1] "A aproximação do LOO-CV foi considerada confiável (k < 0.7 para todas as observações)."loo_resultado_3 = loo(ajusteBay_Stanarm_3)
print(loo_resultado_3)##
## Computed from 4000 by 1000 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -1433.7 22.6
## p_loo 4.1 0.3
## looic 2867.4 45.2
## ------
## MCSE of elpd_loo is 0.0.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.2, 0.9]).
##
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.# Verificar se houve problemas na aproximação
pareto_k_ids_3 = which(loo_resultado_3$diagnostics$pareto_k > 0.7)
if (length(pareto_k_ids_3) == 0) {
print("A aproximação do LOO-CV foi considerada confiável (k < 0.7 para todas as observações).")
} else {
print(paste("Algumas observações apresentaram k > 0.7. Verifique os índices: ", pareto_k_ids, collapse = ", "))
}## [1] "A aproximação do LOO-CV foi considerada confiável (k < 0.7 para todas as observações)."loo_compare(loo_resultado_1, loo_resultado_2, loo_resultado_3)## elpd_diff se_diff
## ajusteBay_Stanarm_3 0.0 0.0
## ajusteBay_Stanarm_1 -1403.2 39.5
## ajusteBay_Stanarm_2 -1607.4 30.43.2 DIC – Deviance Information Criterion
O DIC () é um critério bayesiano para seleção de modelos que combina uma medida de ajuste com uma penalização pela complexidade do modelo, de forma semelhante ao AIC na abordagem frequentista.
Seja \(\theta\) o vetor de parâmetros do modelo, e seja \(p(y \mid \theta)\) a função de verossimilhança dos dados \(y\), a deviance é definida como: \[ D(\theta) = -2 \log p(y \mid \theta). \]
A esperança da deviance sob a distribuição a posteriori dos parâmetros é: \[ \overline{D} = \mathbb{E}_{\theta \mid y}[D(\theta)] = \frac{1}{S} \sum_{s=1}^S D(\theta^{(s)}) \] onde \(\theta^{(s)}\) são amostras da posteriori. A complexidade efetiva do modelo, \(p_D\), é estimada como: \[ p_D = \overline{D} - D(\hat{\theta}) \] onde \(\hat{\theta}\) é a média posterior de \(\theta\). O DIC é então dado por: \[ \mathrm{DIC} = D(\hat{\theta}) + 2p_D. \]
Modelos com menor valor de DIC são preferidos. Empiricamente, considera-se que diferenças maiores que 10 entre dois modelos indicam forte evidência a favor do modelo com menor DIC.
# Função para calcular DIC
calc_dic <- function(modelo) {
# Extrair as amostras da posterior
posterior <- as.matrix(modelo)
# Obter as predições médias (esperança posterior dos valores preditos)
X <- model.matrix(modelo)
beta_cols <- colnames(posterior)[-length(colnames(posterior))]#grep("^b", colnames(posterior), value = TRUE)
pred_mean <- X %*% colMeans(posterior[, beta_cols, drop = FALSE])
# Desvio padrão posterior do erro (sigma)
sigma_post_mean <- mean(posterior[, "sigma"])
# Deviance no ponto da média posterior
deviance_hat <- -2 * sum(dnorm(modelo$y, mean = pred_mean, sd = sigma_post_mean, log = TRUE))
# Deviance média posterior
loglik <- log_lik(modelo)
deviance_samples <- -2 * loglik
mean_deviance <- mean(rowSums(deviance_samples))
# Parâmetro de complexidade efetiva
pD <- mean_deviance - deviance_hat
# DIC final
dic <- deviance_hat + 2 * pD
return(dic)
}
dic1 <- calc_dic(ajusteBay_Stanarm_1)
dic2 <- calc_dic(ajusteBay_Stanarm_2)
dic3 <- calc_dic(ajusteBay_Stanarm_3)
dic1; dic2; dic3## [1] 5673.751## [1] 6082.149## [1] 2867.3683.3 Fator de Bayes
O Fator de Bayes (ou Bayes Factor) compara dois modelos \(M_1\) e \(M_2\) a partir da razão das evidências (ou marginais dos dados) sob cada modelo. A evidência de um modelo é dada pela integral da verossimilhança ponderada pela prior: \[ p(y \mid M_j) = \int_{\Theta_j} p(y \mid \theta_j, M_j) \, p(\theta_j \mid M_j) \, d\theta_j \]
O Fator de Bayes comparando \(M_1\) com \(M_2\) é definido como: \[ BF_{12} = \frac{p(y \mid M_1)}{p(y \mid M_2)}. \]
Sua interpretação é direta:
- \(BF_{12} > 1\): evidência a favor do modelo \(M_1\)
- \(BF_{12} < 1\): evidência a favor do modelo \(M_2\)
A escala de interpretação sugerida por Kass e Raftery (1995) é:
A principal limitação do Fator de Bayes está na sensibilidade às escolhas das distribuições a priori. Além disso, o cálculo da evidência marginal pode ser computacionalmente custoso, especialmente em modelos com alta dimensionalidade, sendo comum o uso de métodos como para sua estimativa.
if (!require("bridgesampling")) install.packages("bridgesampling")
library(bridgesampling)
bridge1 <- bridge_sampler(ajusteBay_Stanarm_1)## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4bridge2 <- bridge_sampler(ajusteBay_Stanarm_2)## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 6bridge3 <- bridge_sampler(ajusteBay_Stanarm_3)## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 6
## Iteration: 7# Comparando modelo 3 (x1 + x2) com os outros
bf_13 <- bf(bridge3, bridge1) # modelo3 vs modelo1
bf_23 <- bf(bridge3, bridge2) # modelo3 vs modelo2
print(bf_13)## Estimated Bayes factor in favor of bridge3 over bridge1: Infprint(bf_23)## Estimated Bayes factor in favor of bridge3 over bridge2: Inf