Aula de Modelos Lineares I
Patrícia Lusié Velozo da Costa
2025-05-12
Modelo de Regressão Linear Simples - Inferência Bayesiana
Modelo
Considere o modelo de regressão linear simples: \[Y_i = \beta_0 + \beta_1 X_i + e_i, \quad e_i \stackrel{iid}{\sim} N(0,\sigma^2), \quad i=1,\ldots,n.\]
1ª tarefa - Simulando conjunto de dados
Simule um conjunto de dados do modelo de regressão linear simples considerando n=1000 (tamanho amostral), \(\beta_0=10\), \(\beta_1=2\), \(\sigma^2=0,2\) e que \(x_i \stackrel{iid}{\sim} N(0,1)\). Em seguida, faça o gráfico de dispersão de X e Y.
# Definindo parâmetros reais
n = 1000 # Tamanho da amostra
beta0 = 10 # Intercepto
beta1 = 2 # Coeficiente de regressão
sigma2 = 0.2
#Simulando os dados
x = rnorm(n, mean = 0, sd = 1) # Gerando variável independente X
e = rnorm(n, mean = 0, sd = sqrt(sigma2))
y = beta0 + beta1 * x + e # Gerando variável dependente Y com erro aleatório
# Gráfico de dispersão de X e Y
library(ggplot2)
ggplot(data = data.frame(x, y), aes(x = x, y = y)) +
geom_point(color = 'blue', alpha = 0.5) +
labs(title = "Gráfico de Dispersão entre X e Y",
x = "X",
y = "Y") +
theme_minimal()
plot(x,y, col="blue", bty="n", xlab="X", ylab="Y", main="Gráfico de Dispersão entre X e Y")
2ª tarefa - Ajuste 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)\) e \(\sigma^2 \sim Exp(1)\), ajuste o modelo usando o pacote rstanarm.
if(!require("rstanarm")) {install.packages("rstanarm")}
library(rstanarm)
library(parallel)
#criando um data frame com os dados
dados = data.frame(y,x)
# Marcando o tempo inicial
tempo_inicial <- Sys.time()
# Fazendo o ajuste
ajusteBay_Stanarm = stan_glm(y ~ x,
data = dados,
family = gaussian(),
chains = 2, iter = 1000,
prior = normal(0, 1000, autoscale = FALSE),
prior_intercept = normal(0, 1000, autoscale = FALSE),
prior_aux = exponential(1) ,
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.480852 secs# analisando o ajuste
summary(ajusteBay_Stanarm)##
## Model Info:
## function: stan_glm
## family: gaussian [identity]
## formula: y ~ x
## algorithm: sampling
## sample: 1000 (posterior sample size)
## priors: see help('prior_summary')
## observations: 1000
## predictors: 2
##
## Estimates:
## mean sd 10% 50% 90%
## (Intercept) 10.0 0.0 10.0 10.0 10.0
## x 2.0 0.0 2.0 2.0 2.0
## sigma 0.4 0.0 0.4 0.4 0.5
##
## Fit Diagnostics:
## mean sd 10% 50% 90%
## mean_PPD 10.1 0.0 10.0 10.1 10.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 236
## x 0.0 1.0 1220
## sigma 0.0 1.0 149
## mean_PPD 0.0 1.0 351
## log-posterior 0.1 1.0 145
##
## 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).prior_summary(ajusteBay_Stanarm)## Priors for model 'ajusteBay_Stanarm'
## ------
## Intercept (after predictors centered)
## ~ normal(location = 0, scale = 1000)
##
## Coefficients
## ~ normal(location = 0, scale = 1000)
##
## Auxiliary (sigma)
## ~ exponential(rate = 1)
## ------
## See help('prior_summary.stanreg') for more details# Extraindo amostras
amostras = as.array(ajusteBay_Stanarm)
# Acessando as cadeias
amostra_Stanarm_beta0_c1 = amostras[, 1, "(Intercept)"] # Primeira cadeia de beta0
amostra_Stanarm_beta0_c2 = amostras[, 2, "(Intercept)"] # Segunda cadeia de beta0
amostra_Stanarm_beta1_c1 = amostras[, 1, "x"] # Primeira cadeia de beta1
amostra_Stanarm_beta1_c2 = amostras[, 2, "x"] # Segunda cadeia de beta1
amostra_Stanarm_sigma2_c1 = amostras[, 1, "sigma"]^2 # Primeira cadeia de sigma2
amostra_Stanarm_sigma2_c2 = amostras[, 2, "sigma"]^2 # Segunda cadeia de sigma2
par(mfcol=c(3,3))
plot(amostra_Stanarm_beta0_c1, xlab="iteração", ylab=expression(beta[0]), type="l",bty="n", ylim=c(min(amostra_Stanarm_beta0_c1,amostra_Stanarm_beta0_c2),max(amostra_Stanarm_beta0_c1,amostra_Stanarm_beta0_c2)))
lines(amostra_Stanarm_beta0_c2, col=2)
#
plot(amostra_Stanarm_beta1_c1, xlab="iteração", ylab=expression(beta[1]), type="l",bty="n", ylim=c(min(amostra_Stanarm_beta1_c1,amostra_Stanarm_beta1_c2),max(amostra_Stanarm_beta1_c1,amostra_Stanarm_beta1_c2)))
lines(amostra_Stanarm_beta1_c2, col=2)
#
plot(amostra_Stanarm_sigma2_c1, xlab="iteração", ylab=expression(sigma^2), type="l",bty="n", ylim=c(min(amostra_Stanarm_sigma2_c1,amostra_Stanarm_sigma2_c2),max(amostra_Stanarm_sigma2_c1,amostra_Stanarm_sigma2_c2)))
lines(amostra_Stanarm_sigma2_c2, col=2)
hist(amostra_Stanarm_beta0_c1, ylab="densidade", xlab=expression(beta[0]), freq=F,main="")
abline(v=beta0, lwd=2, col="blue")
abline(v=mean(amostra_Stanarm_beta0_c1), lwd=2, col="red", lty=5)
abline(v=quantile(amostra_Stanarm_beta0_c1, 0.025), col="red", lty=3, lwd=2)
abline(v=quantile(amostra_Stanarm_beta0_c1, 0.975), col="red", lty=3, lwd=2)
#
hist(amostra_Stanarm_beta1_c1, ylab="densidade", xlab=expression(beta[1]), freq=F,main="")
abline(v=beta1, lwd=2, col="blue")
abline(v=mean(amostra_Stanarm_beta1_c1), lwd=2, col="red", lty=5)
abline(v=quantile(amostra_Stanarm_beta1_c1, 0.025), col="red", lty=3, lwd=2)
abline(v=quantile(amostra_Stanarm_beta1_c1, 0.975), col="red", lty=3, lwd=2)
#
hist(amostra_Stanarm_sigma2_c1, ylab="densidade", xlab=expression(sigma^2), freq=F,main="")
abline(v=sigma2, lwd=2, col="blue")
abline(v=mean(amostra_Stanarm_sigma2_c1), lwd=2, col="red", lty=5)
abline(v=quantile(amostra_Stanarm_sigma2_c1, 0.025), col="red", lty=3, lwd=2)
abline(v=quantile(amostra_Stanarm_sigma2_c1, 0.975), col="red", lty=3, lwd=2)
acf(amostra_Stanarm_beta0_c1, xlab="defasagem", main="")
acf(amostra_Stanarm_beta1_c1, xlab="defasagem", main="")
acf(amostra_Stanarm_sigma2_c1, xlab="defasagem", main="")
3ª tarefa - Ajuste usando brms
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)\) e \(\sigma^2 \sim Exp(1)\), ajuste o modelo usando o pacote brms.
# Instalar o pacote brms se necessário
if (!requireNamespace("brms", quietly = TRUE)) {
install.packages("brms")
}
library(brms)
#criando um data frame com os dados
dados = data.frame(y,x)
# Definindo total de cadeias
ncadeias = 2
# Marcando o tempo inicial
tempo_inicial <- Sys.time()
# Realizando o ajuste
ajusteBay_brms = brm(
formula = y ~ x, # Fórmula do modelo
data = dados,
family = gaussian(), # Família de distribuição para y
prior = c(set_prior("normal(0, 1000)", class = "b"), #para coeficientes
set_prior("normal(0, 1000)", class = "Intercept"), #para o intercepto
set_prior("exponential(1)", class = "sigma")), #Para o desvio padrao
chains = ncadeias, # Número de cadeias MCMC
iter = 1000, # Número de iterações
warmup = 100, # Número de iterações de aquecimento
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.101714 mins# Resumo do modelo ajustado
summary(ajusteBay_brms)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: y ~ x
## Data: dados (Number of observations: 1000)
## Draws: 2 chains, each with iter = 1000; warmup = 100; thin = 1;
## total post-warmup draws = 1800
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 9.99 0.01 9.96 10.02 1.00 2682 1333
## x 2.02 0.01 1.99 2.05 1.00 2464 1339
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.45 0.01 0.43 0.47 1.00 1051 1226
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).# Extraindo as amostras dos parâmetros
amostras = as.data.frame(as.matrix(ajusteBay_brms))
# Salvando as amostras em matrizes separadas
amostra_brms_sigma2 = matrix(amostras$sigma^2, ncol = ncadeias, byrow = FALSE)
amostra_brms_beta0 = matrix(amostras$`b_Intercept`, ncol = ncadeias, byrow = FALSE)
amostra_brms_beta1 = matrix(amostras$`b_x`, ncol = ncadeias, byrow = FALSE)
# Renomeando as colunas para indicar a cadeia
colnames(amostra_brms_sigma2) <- paste0("Cadeia_", 1:ncadeias)
colnames(amostra_brms_beta0) <- paste0("Cadeia_", 1:ncadeias)
colnames(amostra_brms_beta1) <- paste0("Cadeia_", 1:ncadeias)
par(mfcol=c(3,3))
plot(amostra_brms_beta0[,1], xlab="iteração", ylab=expression(beta[0]), type="l",bty="n")
lines(amostra_brms_beta0[,2], col=2)
plot(amostra_brms_beta1[,1], xlab="iteração", ylab=expression(beta[1]), type="l",bty="n")
lines(amostra_brms_beta1[,2], col=2)
plot(amostra_brms_sigma2[,1], xlab="iteração", ylab=expression(sigma^2), type="l",bty="n")
lines(amostra_brms_sigma2[,2], col=2)
#
hist(amostra_brms_beta0[,1], ylab="densidade", xlab=expression(beta[0]), freq=F,main="")
abline(v=beta0, lwd=2, col="blue")
abline(v=mean(amostra_brms_beta0[,1]), lwd=2, col="red", lty=5)
abline(v=quantile(amostra_brms_beta0[,1], 0.025), col="red", lty=3, lwd=2)
abline(v=quantile(amostra_brms_beta0[,1], 0.975), col="red", lty=3, lwd=2)
#
hist(amostra_brms_beta1[,1], ylab="densidade", xlab=expression(beta[1]), freq=F,main="")
abline(v=beta1, lwd=2, col="blue")
abline(v=mean(amostra_brms_beta1[,1]), lwd=2, col="red", lty=5)
abline(v=quantile(amostra_brms_beta1[,1], 0.025), col="red", lty=3, lwd=2)
abline(v=quantile(amostra_brms_beta1[,1], 0.975), col="red", lty=3, lwd=2)
#
hist(amostra_brms_sigma2[,1], ylab="densidade", xlab=expression(sigma^2), freq=F,main="")
abline(v=sigma2, lwd=2, col="blue")
abline(v=mean(amostra_brms_sigma2[,1]), lwd=2, col="red", lty=5)
abline(v=quantile(amostra_brms_sigma2[,1], 0.025), col="red", lty=3, lwd=2)
abline(v=quantile(amostra_brms_sigma2[,1], 0.975), col="red", lty=3, lwd=2)
acf(amostra_brms_beta0[,1], main="", xlab="defasagem")
acf(amostra_brms_beta1[,1], main="", xlab="defasagem")
acf(amostra_brms_sigma2[,1], main="", xlab="defasagem")
4ª tarefa - Ajuste usando rstan
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)\) e \(\sigma^2 \sim GI(2,1)\), ajuste o modelo usando o pacote rstan.
#crie um arquivo chamado ML1_Aula_06_p1_Stan_modelo.stan com os scripts abaixo
data {
int<lower=0> n; // Tamanho amostral
vector[n] x; // Variável preditora
vector[n] y; // Variável resposta
}
parameters {
real beta_0; // Intercepto
real beta_1; // Coeficiente
real<lower=0> sigma2; // Variância residual
}
model {
// Priori para os parâmetros
beta_0 ~ normal(0, sqrt(1000));
beta_1 ~ normal(0, sqrt(1000));
sigma2 ~ inv_gamma(2, 1); // Prior gama inversa para a variância
// Verossimilhança
y ~ normal(beta_0 + beta_1 * x, sqrt(sigma2));
}#Agora execute em um arquivo com extensão r os comandos abaixo
library(rstan)
# Definir os dados para o Stan
dados_stan <- list(n = n, x = x, y = y)
# Marcar o tempo inicial
tempo_inicial <- Sys.time()
# Compilar e ajustar o modelo
ajusteBayStan <- stan(
file = 'ML1_Aula_06_p1_Stan_modelo.stan', # Arquivo Stan
data = dados_stan, # Dados
chains = 2, # Número de cadeias MCMC
iter = 1000, # Número de iterações
warmup = 100, # Iterações de aquecimento
refresh = 0 # Suprime a saída de progresso
)
# Marcar o tempo final
tempo_final <- Sys.time()
tempo_final - tempo_inicial## Time difference of 1.059609 mins# Resumo dos parâmetros estimados
print(ajusteBayStan, pars = c("beta_0", "beta_1", "sigma2"))## Inference for Stan model: anon_model.
## 2 chains, each with iter=1000; warmup=100; thin=1;
## post-warmup draws per chain=900, total post-warmup draws=1800.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## beta_0 9.99 0 0.01 9.96 9.98 9.99 10.00 10.02 2571 1
## beta_1 2.02 0 0.01 1.99 2.01 2.02 2.03 2.05 2490 1
## sigma2 0.20 0 0.01 0.18 0.19 0.20 0.21 0.22 548 1
##
## Samples were drawn using NUTS(diag_e) at Mon May 12 12:12:02 2025.
## For each parameter, 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).# Extraindo os parametros estimados
amostra_posteriori = extract(ajusteBayStan, permuted = FALSE)
amostra_beta0 = amostra_posteriori[, , "beta_0"]
amostra_beta1 = amostra_posteriori[, , "beta_1"]
amostra_sigma2 = amostra_posteriori[, , "sigma2"]
par(mfcol=c(3,3))
plot(amostra_beta0[,1], xlab="iteração", ylab=expression(beta[0]), type="l",bty="n")
lines(amostra_beta0[,2], col=2)
plot(amostra_beta1[,1], xlab="iteração", ylab=expression(beta[1]), type="l",bty="n")
lines(amostra_beta1[,2], col=2)
plot(amostra_sigma2[,1], xlab="iteração", ylab=expression(sigma^2), type="l",bty="n")
lines(amostra_sigma2[,2], col=2)
hist(amostra_beta0[,1], ylab="densidade", xlab=expression(beta[0]), freq=F,main="")
abline(v=beta0, lwd=2, col="blue")
abline(v=mean(amostra_beta0[,1]), lwd=2, col="red", lty=5)
abline(v=quantile(amostra_beta0[,1], 0.025), col="red", lty=3, lwd=2)
abline(v=quantile(amostra_beta0[,1], 0.975), col="red", lty=3, lwd=2)
#
hist(amostra_beta1[,1], ylab="densidade", xlab=expression(beta[1]), freq=F,main="")
abline(v=beta1, lwd=2, col="blue")
abline(v=mean(amostra_beta1[,1]), lwd=2, col="red", lty=5)
abline(v=quantile(amostra_beta1[,1], 0.025), col="red", lty=3, lwd=2)
abline(v=quantile(amostra_beta1[,1], 0.975), col="red", lty=3, lwd=2)
#
hist(amostra_sigma2[,1], ylab="densidade", xlab=expression(sigma^2), freq=F,main="")
abline(v=sigma2, lwd=2, col="blue")
abline(v=mean(amostra_sigma2[,1]), lwd=2, col="red", lty=5)
abline(v=quantile(amostra_sigma2[,1], 0.025), col="red", lty=3, lwd=2)
abline(v=quantile(amostra_sigma2[,1], 0.975), col="red", lty=3, lwd=2)
acf(amostra_beta0[,1], main="", xlab="defasagem")
acf(amostra_beta1[,1], main="", xlab="defasagem")
acf(amostra_sigma2[,1], main="", xlab="defasagem")