Notas de aula de Estatística Bayesiana
Lia Hanna Martins Morita
2020
1 Aplicações com o pacote RJAGS
- JAGS é sigla para “Just Another Gibbs Sampler”, para análise de modelos hierárquicos Bayesianos usando simulação de Markov Chain Monte Carlo (MCMC);
- CODA é o pacote de “Output Analysis and Diagnostics for MCMC” para ser utilizado juntamente com o JAGS para criação de gráficos de amostras simuladas via Markov Chain Monte Carlo (MCMC), além de testes de diagnóstico de convergência para a distribuição de equilíbrio das cadeias de Markov;
- Instalação dos pacotes:
#install.packages("rjags")
#install.packages("coda")
require(rjags)
require(coda)
require(DT)
require(tidyverse)
require(car)2 Retomando os dados do Exemplo 3.1:
- A entrada dos dados é em formato de lista;
- A sintaxe do modelo é armazenada em um documento .txt;
- Prioris normais não informativas para \(\beta_0\) e \(\beta_1\) com média igual a zero e variância igual a 100000;
- Priori gamma não informativa para \(\tau\) com média igual a 1 e variância igual a 100.
- Atribui-se valores iniciais para os parâmetros desconhecidos \(\beta_0\), \(\beta_1\) e \(\tau\): que terão distribuição a priori;
##### Entrada dos dados em formato de lista
dados=list(n=10,x=c(3,3,4,5,6,6,7,8,8,9),
y=c(9,5,12,9,14,16,22,18,24,22))
##### Sintaxe do modelo
cat("model {
for (i in 1:n)
{
y[i] ~ dnorm(beta0 + x[i] *beta1,tau)
y_hat[i] <- beta0 + x[i] *beta1
}
beta0 ~ dnorm(0,0.00001)
beta1 ~ dnorm(0,0.00001)
tau ~ dgamma(0.01,0.01)
sigma2 = 1/tau
}",file="modelo1.txt")
##### Valores iniciais
inits = list(beta0=0,beta1=0,tau=1)
#### Carrega o modelo
jags.m <- jags.model( file = "modelo1.txt", data=dados, inits=inits, n.chains=2, n.adapt=500 )## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 10
## Unobserved stochastic nodes: 3
## Total graph size: 43
##
## Initializing model## Especifica as quantidades a serem monitorados (algumas vêm direto dos parâmetros que foram amostrados, e outras são funções dos parâmetros)
params <- c("y_hat","beta0","beta1","tau","sigma2")
## Implementa o modelo e armazena as amostras da distribuição a posteriori
samps <- coda.samples( jags.m, params, n.iter=10000 )
## Pede os sumários a posteriori
summary(samps)##
## Iterations = 1:10000
## Thinning interval = 1
## Number of chains = 2
## Sample size per chain = 10000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## beta0 -1.1035 3.18337 0.0225098 0.0960576
## beta1 2.7456 0.51063 0.0036107 0.0152924
## sigma2 10.4695 6.91125 0.0488699 0.0870838
## tau 0.1263 0.06282 0.0004442 0.0006819
## y_hat[1] 7.1334 1.80050 0.0127314 0.0495064
## y_hat[2] 7.1334 1.80050 0.0127314 0.0495064
## y_hat[3] 9.8790 1.40976 0.0099685 0.0339763
## y_hat[4] 12.6246 1.12008 0.0079201 0.0188991
## y_hat[5] 15.3702 1.02137 0.0072222 0.0066711
## y_hat[6] 15.3702 1.02137 0.0072222 0.0066711
## y_hat[7] 18.1158 1.16332 0.0082259 0.0130001
## y_hat[8] 20.8614 1.47815 0.0104521 0.0272085
## y_hat[9] 20.8614 1.47815 0.0104521 0.0272085
## y_hat[10] 23.6070 1.88097 0.0133005 0.0426371
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## beta0 -7.25657 -3.11291 -1.1792 0.8143 5.528
## beta1 1.68547 2.43701 2.7575 3.0678 3.732
## sigma2 3.63695 6.20507 8.6364 12.4663 28.545
## tau 0.03503 0.08022 0.1158 0.1612 0.275
## y_hat[1] 3.61299 6.00161 7.1066 8.2190 10.854
## y_hat[2] 3.61299 6.00161 7.1066 8.2190 10.854
## y_hat[3] 7.10604 9.00024 9.8624 10.7226 12.790
## y_hat[4] 10.41382 11.93508 12.6142 13.3017 14.898
## y_hat[5] 13.34338 14.74660 15.3754 15.9956 17.395
## y_hat[6] 13.34338 14.74660 15.3754 15.9956 17.395
## y_hat[7] 15.77227 17.40673 18.1250 18.8262 20.445
## y_hat[8] 17.90995 19.95466 20.8712 21.7768 23.817
## y_hat[9] 17.90995 19.95466 20.8712 21.7768 23.817
## y_hat[10] 19.81602 22.46487 23.6230 24.7827 27.325- Traço: é a série dos valores dos parâmetros ao longo das iterações do algoritmo;
- Densidade: é o gráfico da densidade dos parâmetros;
plot(samps)
# Análise de dados com regressão linear múltipla:
- Pendente em verificar a fonte e explicação dos dados!
- Verificar Transformação Box Cox
- Referências: Box, G. E. P. and Cox, D. R. (1964) An analysis of transformations. Journal of the Royal Statisistical Society, Series B. 26 211-46.
- <a href=“https://www.ime.usp.br/~abe/lista/pdfm9cJKUmFZp.pdf” < Material sobre Box Cox
tempos=c(195,255,195,255,225,225,195,255,225,225,225,225,230,230)
ganhos=c(1004,1636,852,1506,1275,1270,1269,903,1555,1260,1146,1276,1225,1321)
doses=c(4,4,4.6,4.6,4.2,4.1,4.6,4.3,4.3,4,4.7,4.3,4.72,4.3)
dados=data.frame(tempos,doses,ganhos)
datatable(dados)3 Ajuste do modelo linear múltiplo
- Os resíduos são aproximadamente normais, mas o gráfico de resíduos x preditos apresenta problemas;
- Sugestão: transformações nos dados.
- Transformação Box Cox
- Referências: Box, G. E. P. and Cox, D. R. (1964) An analysis of transformations. Journal of the Royal Statisistical Society, Series B. 26 211-46;
- <a href=“https://www.ime.usp.br/~abe/lista/pdfm9cJKUmFZp.pdf” < Material sobre Box Cox
- A conclusão do modelo é que os tempo & dose não são estatisticamente significativos para explicar os ganhos.
model1=lm(ganhos ~ tempos + doses,data=dados)
summary(model1)##
## Call:
## lm(formula = ganhos ~ tempos + doses, data = dados)
##
## Residuals:
## Min 1Q Median 3Q Max
## -496.38 -50.20 5.42 117.26 304.41
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 629.382 1246.225 0.505 0.624
## tempos 4.960 2.938 1.688 0.120
## doses -115.052 227.600 -0.506 0.623
##
## Residual standard error: 215.6 on 11 degrees of freedom
## Multiple R-squared: 0.2317, Adjusted R-squared: 0.09199
## F-statistic: 1.658 on 2 and 11 DF, p-value: 0.2347anova(model1)## Analysis of Variance Table
##
## Response: ganhos
## Df Sum Sq Mean Sq F value Pr(>F)
## tempos 1 142313 142313 3.0615 0.1080
## doses 1 11878 11878 0.2555 0.6232
## Residuals 11 511338 46485plot(model1)
shapiro.test(residuals(model1))##
## Shapiro-Wilk normality test
##
## data: residuals(model1)
## W = 0.92649, p-value = 0.2722dados = dados %>%
mutate(y_hat=fitted(model1),
res = residuals(model1))
datatable(dados)#Transformação Box Cox
a=boxCox(model1)
maximo=cbind(a$x,a$y)
lambda=maximo[which.max(maximo[,2])]
#1.59596
model2 = lm((ganhos^lambda-1)/lambda ~ tempos + doses,data=dados)
summary(model2)##
## Call:
## lm(formula = (ganhos^lambda - 1)/lambda ~ tempos + doses, data = dados)
##
## Residuals:
## Min 1Q Median 3Q Max
## -33733 -4142 -402 8089 22042
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 11162.9 85883.1 0.130 0.899
## tempos 355.6 202.5 1.756 0.107
## doses -8239.5 15685.0 -0.525 0.610
##
## Residual standard error: 14860 on 11 degrees of freedom
## Multiple R-squared: 0.246, Adjusted R-squared: 0.1089
## F-statistic: 1.794 on 2 and 11 DF, p-value: 0.2116anova(model2)## Analysis of Variance Table
##
## Response: (ganhos^lambda - 1)/lambda
## Df Sum Sq Mean Sq F value Pr(>F)
## tempos 1 731370389 731370389 3.3128 0.09603 .
## doses 1 60921234 60921234 0.2759 0.60979
## Residuals 11 2428463475 220769407
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(model2)
shapiro.test(residuals(model2))##
## Shapiro-Wilk normality test
##
## data: residuals(model2)
## W = 0.9344, p-value = 0.35134 O modelo Bayesiano para os dados das doses
##### Entrada dos dados em formato de lista
dados=list(n=14,tempos=c(195,255,195,255,225,225,195,255,225,225,225,225,230,230),
ganhos=c(1004,1636,852,1506,1275,1270,1269,903,1555,1260,1146,1276,1225,1321),
doses=c(4,4,4.6,4.6,4.2,4.1,4.6,4.3,4.3,4,4.7,4.3,4.72,4.3))
##### Sintaxe do modelo
cat("model {
for (i in 1:n)
{
ganhos[i] ~ dnorm(beta0 + tempos[i]*beta1 + doses[i]*beta2 ,tau)
y_hat[i] <- beta0 + tempos[i]*beta1 + doses[i]*beta2
}
beta0 ~ dnorm(0,0.00001)
beta1 ~ dnorm(0,0.00001)
beta2 ~ dnorm(0,0.00001)
tau ~ dgamma(0.01,0.01)
sigma2 = 1/tau
}",file="modelo2.txt")
##### Valores iniciais
inits = list(beta0=0,beta1=0,beta2=0,tau=1)
#### Carrega o modelo
jags.m <- jags.model( file = "modelo2.txt", data=dados, inits=inits, n.chains=2, n.adapt=1000 )## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 14
## Unobserved stochastic nodes: 4
## Total graph size: 75
##
## Initializing model## Especifica as quantidades a serem monitorados (algumas vêm direto dos parâmetros que foram amostrados, e outras são funções dos parâmetros)
params <- c("y_hat","beta0","beta1","beta2","tau","sigma2")
## Implementa o modelo e armazena as amostras da distribuição a posteriori
samps <- coda.samples( jags.m, params, n.iter=10000 )
## Pede os sumários a posteriori
summary(samps)##
## Iterations = 1:10000
## Thinning interval = 1
## Number of chains = 2
## Sample size per chain = 10000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## beta0 3.520e+01 2.996e+02 2.118e+00 1.444e+01
## beta1 5.539e+00 2.289e+00 1.619e-02 1.930e-01
## beta2 -8.365e+00 1.203e+02 8.508e-01 1.031e+01
## sigma2 5.167e+04 2.524e+04 1.784e+02 3.613e+02
## tau 2.307e-05 9.247e-06 6.538e-08 1.240e-07
## y_hat[1] 1.082e+03 6.982e+01 4.937e-01 2.406e+00
## y_hat[2] 1.414e+03 1.197e+02 8.462e-01 8.967e+00
## y_hat[3] 1.077e+03 1.148e+02 8.120e-01 8.799e+00
## y_hat[4] 1.409e+03 7.877e+01 5.570e-01 2.724e+00
## y_hat[5] 1.246e+03 6.129e+01 4.334e-01 8.333e-01
## y_hat[6] 1.247e+03 6.461e+01 4.569e-01 1.843e+00
## y_hat[7] 1.077e+03 1.148e+02 8.120e-01 8.799e+00
## y_hat[8] 1.412e+03 9.466e+01 6.694e-01 5.489e+00
## y_hat[9] 1.245e+03 6.023e+01 4.259e-01 4.259e-01
## y_hat[10] 1.248e+03 6.987e+01 4.941e-01 2.819e+00
## y_hat[11] 1.242e+03 7.752e+01 5.482e-01 3.472e+00
## y_hat[12] 1.245e+03 6.023e+01 4.259e-01 4.259e-01
## y_hat[13] 1.270e+03 7.433e+01 5.256e-01 2.861e+00
## y_hat[14] 1.273e+03 6.214e+01 4.394e-01 8.278e-01
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## beta0 -5.493e+02 -1.612e+02 3.610e+01 2.270e+02 6.296e+02
## beta1 1.259e+00 3.990e+00 5.520e+00 7.086e+00 1.009e+01
## beta2 -2.468e+02 -8.637e+01 -9.649e+00 7.202e+01 2.293e+02
## sigma2 2.256e+04 3.505e+04 4.587e+04 6.134e+04 1.151e+05
## tau 8.686e-06 1.630e-05 2.180e-05 2.853e-05 4.433e-05
## y_hat[1] 9.414e+02 1.036e+03 1.082e+03 1.127e+03 1.219e+03
## y_hat[2] 1.183e+03 1.335e+03 1.415e+03 1.493e+03 1.653e+03
## y_hat[3] 8.484e+02 1.001e+03 1.076e+03 1.152e+03 1.304e+03
## y_hat[4] 1.252e+03 1.358e+03 1.409e+03 1.461e+03 1.566e+03
## y_hat[5] 1.124e+03 1.207e+03 1.247e+03 1.286e+03 1.368e+03
## y_hat[6] 1.118e+03 1.206e+03 1.247e+03 1.289e+03 1.374e+03
## y_hat[7] 8.484e+02 1.001e+03 1.076e+03 1.152e+03 1.304e+03
## y_hat[8] 1.227e+03 1.350e+03 1.411e+03 1.474e+03 1.598e+03
## y_hat[9] 1.125e+03 1.207e+03 1.245e+03 1.284e+03 1.365e+03
## y_hat[10] 1.108e+03 1.203e+03 1.248e+03 1.293e+03 1.387e+03
## y_hat[11] 1.089e+03 1.192e+03 1.242e+03 1.293e+03 1.397e+03
## y_hat[12] 1.125e+03 1.207e+03 1.245e+03 1.284e+03 1.365e+03
## y_hat[13] 1.122e+03 1.221e+03 1.270e+03 1.318e+03 1.418e+03
## y_hat[14] 1.149e+03 1.233e+03 1.273e+03 1.313e+03 1.396e+03plot(samps)
