Unidade IV - Exemplo II
Com pacotes extras para fazer gráfico
require(rjags)
require(coda)
require(DT)
require(tidyverse)
require(car)
#install.packages("MCMCpack")
require(MCMCpack)
#install.packages("bayesplot")
require(bayesplot)
#install.packages("ggrepel")
require(grid)Exemplo 3:
Modelo de regressão não linear
Dugongues são animais da mesma ordem dos peixes-bois (ordem Sirenia). Dentre suas principais características, temos: narinas no topo da cabeça, lábio superior voltado para baixo e nadadeira dividida em duas partes (como a das baleias e golfinhos). Os Dugongues são encontrados na costa leste da Àfrica, Ìndia, Indonésia, Malásia e Austrália.
- referência na página 57 .
- Abaixo temos os tamanhos & idades de 27 dugongues. Copiar os dados em um txt no seu computador
- Carlin & Gelfand (1991) propõem o seguinte modelo para estes dados: \[ Y_i \sim N(\mu_i,\sigma^2)\\ \mu_i = \alpha - \beta \gamma^{x_i} \]
Leitura dos dados & representação gráfica
- gráfico de dispersão de y versus x
dugongues=read.table("https://raw.githubusercontent.com/liamorita/estatisticabayesiana/master/exemplo3_dugongues.txt",sep=";",head=TRUE)
datatable(dugongues)dugongues %>%
ggplot(aes(x=x,y=y)) +
geom_point()
Ajuste do modelo não linear pelo método de máxima verossimilhança, utilizando a função optim do R - minimiza a função log-verossimilhança e realiza a estimação numérica (não analítica). O algoritmo utilizado é o “BFGS”
O algoritmo Broyden – Fletcher – Goldfarb – Shanno (BFGS) é um método iterativo para resolver problemas de otimização não linear irrestrita.
O método BFGS pertence aos métodos quasi-Newton, uma classe de técnicas de otimização que busca um ponto estacionário de uma função (de preferência duas vezes continuamente diferenciável).
“https://stat.ethz.ch/R-manual/R-devel/library/stats/html/optim.html”> Saiba mais sobre optim no R ;
Houve NaNs produzidos que podem ser investigados (não foram exibidos aqui porém se os alunos “rodarem” irão enxergar;
geralmente ocorrem quando ao longo das iterações os valores dos parâmetros assumem valores nulos ou negativos, quando não são por natureza;
No entanto o método apresentou convergência (código igual a zero).
x=dugongues$x
y=dugongues$y
logL=function(par){
alpha=par[1]
beta=par[2]
gamma=par[3]
sigma=par[4]
media = alpha-beta*gamma^x
L = prod(dnorm(y,media,sigma))
(-1)*log(L)
}
nlm_model = optim(logL,par=c(1,1,0.5,5),method="BFGS")
nlm_model## $par
## [1] 2.67067170 0.97575073 0.87382246 0.09115797
##
## $value
## [1] -26.36358
##
## $counts
## function gradient
## 115 29
##
## $convergence
## [1] 0
##
## $message
## NULLEntrada dos dados e valores iniciais
dados=list("n"=nrow(dugongues),"x"=dugongues$x,"y"=dugongues$y)
iniciais=list(
list(alpha=1,beta=1,gamma=0.4,tau=0.5),
list(alpha=2,beta=2,gamma=0.8,tau=1)
)Sintaxe e inicialização do modelo
- o n.adapt é o período de adaptação, ou aquecimento, em inglês: burn-in. Estas amostras serão descartadas pois são muito dependentes dos valores iniciais.
cat("
model{
for (i in 1 : n) {
y[i] ~ dnorm(alpha-beta*pow(gamma,x[i]), tau)
y_pred[i] ~ dnorm(alpha-beta*pow(gamma,x[i]), tau)
}
alpha ~ dnorm(0,0.000001)
beta ~ dnorm(0,0.000001)
gamma ~ dunif(0,100000)
tau ~ dgamma(0.001,0.001)
sigma2 = 1/tau
y_new ~ dnorm(alpha-beta*pow(gamma,26), tau)
}
",file="modelo5.txt")
modelo=jags.model(file="modelo5.txt",data=dados,inits=iniciais,n.chains=2, n.adapt=1000)## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 27
## Unobserved stochastic nodes: 32
## Total graph size: 157
##
## Initializing modelEspecifica as quantidades a serem monitorados
algumas vêm direto dos parâmetros que foram amostrados, e outras são funções dos parâmetros
- Implementa o modelo e armazena as amostras da distribuição a posteriori
- Pode-se pedir o objeto samps - irá retornar todos os valores gerados ao longo das iterações, já considerando o salto (thin) para retirar a dependência, ou autocorrelação entre valores gerados.
params = c("alpha","beta","gamma","sigma2","y_pred","y_new")
samps = coda.samples( modelo, params, n.iter=20000,thin=5)
summary(samps)##
## Iterations = 1005:21000
## Thinning interval = 5
## Number of chains = 2
## Sample size per chain = 4000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## alpha 2.6613 0.070142 7.842e-04 1.907e-03
## beta 0.9832 0.076077 8.506e-04 1.033e-03
## gamma 0.8637 0.031040 3.470e-04 8.329e-04
## sigma2 0.0103 0.003236 3.618e-05 3.889e-05
## y_new 2.6307 0.114830 1.284e-03 1.669e-03
## y_pred[1] 1.8121 0.114127 1.276e-03 1.375e-03
## y_pred[2] 1.8704 0.112083 1.253e-03 1.253e-03
## y_pred[3] 1.8693 0.109338 1.222e-03 1.304e-03
## y_pred[4] 1.8710 0.110652 1.237e-03 1.237e-03
## y_pred[5] 1.9772 0.107078 1.197e-03 1.197e-03
## y_pred[6] 2.1093 0.105145 1.176e-03 1.264e-03
## y_pred[7] 2.1822 0.105758 1.182e-03 1.267e-03
## y_pred[8] 2.1838 0.106559 1.191e-03 1.241e-03
## y_pred[9] 2.2983 0.108235 1.210e-03 1.280e-03
## y_pred[10] 2.3435 0.105265 1.177e-03 1.225e-03
## y_pred[11] 2.3677 0.106802 1.194e-03 1.224e-03
## y_pred[12] 2.3868 0.104682 1.170e-03 1.247e-03
## y_pred[13] 2.4067 0.105873 1.184e-03 1.184e-03
## y_pred[14] 2.4065 0.105490 1.179e-03 1.179e-03
## y_pred[15] 2.4202 0.106682 1.193e-03 1.219e-03
## y_pred[16] 2.4774 0.103235 1.154e-03 1.165e-03
## y_pred[17] 2.4763 0.103397 1.156e-03 1.230e-03
## y_pred[18] 2.5023 0.104085 1.164e-03 1.206e-03
## y_pred[19] 2.5001 0.104711 1.171e-03 1.171e-03
## y_pred[20] 2.5294 0.104969 1.174e-03 1.187e-03
## y_pred[21] 2.5436 0.104689 1.170e-03 1.171e-03
## y_pred[22] 2.5460 0.105082 1.175e-03 1.188e-03
## y_pred[23] 2.5593 0.106283 1.188e-03 1.208e-03
## y_pred[24] 2.5667 0.104236 1.165e-03 1.243e-03
## y_pred[25] 2.6109 0.108815 1.217e-03 1.409e-03
## y_pred[26] 2.6395 0.113271 1.266e-03 1.710e-03
## y_pred[27] 2.6448 0.115402 1.290e-03 1.811e-03
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## alpha 2.53940 2.613345 2.656304 2.70208 2.81531
## beta 0.83417 0.933991 0.983091 1.03166 1.13606
## gamma 0.79263 0.846528 0.867086 0.88488 0.91423
## sigma2 0.00576 0.008039 0.009747 0.01183 0.01822
## y_new 2.39936 2.554749 2.632193 2.70582 2.85762
## y_pred[1] 1.58344 1.735360 1.812327 1.88738 2.03613
## y_pred[2] 1.65439 1.796430 1.869777 1.94328 2.09512
## y_pred[3] 1.65390 1.796076 1.869759 1.93959 2.08543
## y_pred[4] 1.65098 1.798606 1.871741 1.94469 2.09275
## y_pred[5] 1.77297 1.905153 1.977130 2.04549 2.19041
## y_pred[6] 1.90372 2.039097 2.109355 2.17823 2.31443
## y_pred[7] 1.97536 2.111871 2.181302 2.25063 2.39417
## y_pred[8] 1.97547 2.113528 2.182206 2.25260 2.39850
## y_pred[9] 2.08442 2.225636 2.298489 2.36888 2.51545
## y_pred[10] 2.13584 2.274250 2.342002 2.41239 2.54756
## y_pred[11] 2.15884 2.298064 2.365994 2.43699 2.58079
## y_pred[12] 2.18123 2.317433 2.387529 2.45446 2.59585
## y_pred[13] 2.20257 2.336441 2.407111 2.47537 2.61690
## y_pred[14] 2.20092 2.337026 2.406354 2.47714 2.61679
## y_pred[15] 2.20838 2.351378 2.420434 2.48949 2.62892
## y_pred[16] 2.27170 2.409261 2.478260 2.54611 2.67752
## y_pred[17] 2.26561 2.409075 2.477274 2.54433 2.67780
## y_pred[18] 2.30045 2.434011 2.501254 2.56911 2.70889
## y_pred[19] 2.29053 2.432591 2.501175 2.56801 2.70150
## y_pred[20] 2.32071 2.459909 2.530398 2.59763 2.73759
## y_pred[21] 2.33867 2.475504 2.542841 2.61234 2.75050
## y_pred[22] 2.34060 2.477554 2.544613 2.61637 2.75217
## y_pred[23] 2.34684 2.490368 2.559087 2.62915 2.76434
## y_pred[24] 2.35806 2.500366 2.567458 2.63403 2.77045
## y_pred[25] 2.39418 2.540086 2.611608 2.68271 2.82187
## y_pred[26] 2.41254 2.566873 2.639637 2.71322 2.86252
## y_pred[27] 2.41891 2.568296 2.643300 2.72177 2.872911 Outro modo de armazenar
utilizando a função as.mcmc.list do pacote coda
- desta forma o parâmetro “alpha” fica guardado na coluna samps2[,1]
- fica mais fácil para investigar as quantidades separadamente.
samps2=rbind(as.mcmc.list(samps)[[1]],as.mcmc.list(samps)[[2]])
summary(samps2)## alpha beta gamma sigma2
## Min. :2.436 Min. :0.6581 Min. :0.6596 Min. :0.004069
## 1st Qu.:2.613 1st Qu.:0.9340 1st Qu.:0.8465 1st Qu.:0.008039
## Median :2.656 Median :0.9831 Median :0.8671 Median :0.009747
## Mean :2.661 Mean :0.9832 Mean :0.8637 Mean :0.010297
## 3rd Qu.:2.702 3rd Qu.:1.0317 3rd Qu.:0.8849 3rd Qu.:0.011832
## Max. :3.063 Max. :1.3533 Max. :0.9476 Max. :0.036986
## y_new y_pred[1] y_pred[2] y_pred[3]
## Min. :2.167 Min. :1.296 Min. :1.366 Min. :1.461
## 1st Qu.:2.555 1st Qu.:1.735 1st Qu.:1.796 1st Qu.:1.796
## Median :2.632 Median :1.812 Median :1.870 Median :1.870
## Mean :2.631 Mean :1.812 Mean :1.870 Mean :1.869
## 3rd Qu.:2.706 3rd Qu.:1.887 3rd Qu.:1.943 3rd Qu.:1.940
## Max. :3.060 Max. :2.304 Max. :2.446 Max. :2.322
## y_pred[4] y_pred[5] y_pred[6] y_pred[7]
## Min. :1.409 Min. :1.483 Min. :1.633 Min. :1.770
## 1st Qu.:1.799 1st Qu.:1.905 1st Qu.:2.039 1st Qu.:2.112
## Median :1.872 Median :1.977 Median :2.109 Median :2.181
## Mean :1.871 Mean :1.977 Mean :2.109 Mean :2.182
## 3rd Qu.:1.945 3rd Qu.:2.045 3rd Qu.:2.178 3rd Qu.:2.251
## Max. :2.342 Max. :2.528 Max. :2.622 Max. :2.628
## y_pred[8] y_pred[9] y_pred[10] y_pred[11]
## Min. :1.718 Min. :1.854 Min. :1.875 Min. :1.938
## 1st Qu.:2.114 1st Qu.:2.226 1st Qu.:2.274 1st Qu.:2.298
## Median :2.182 Median :2.298 Median :2.342 Median :2.366
## Mean :2.184 Mean :2.298 Mean :2.343 Mean :2.368
## 3rd Qu.:2.253 3rd Qu.:2.369 3rd Qu.:2.412 3rd Qu.:2.437
## Max. :2.646 Max. :2.771 Max. :2.786 Max. :2.806
## y_pred[12] y_pred[13] y_pred[14] y_pred[15]
## Min. :2.004 Min. :1.967 Min. :1.964 Min. :1.935
## 1st Qu.:2.317 1st Qu.:2.336 1st Qu.:2.337 1st Qu.:2.351
## Median :2.388 Median :2.407 Median :2.406 Median :2.420
## Mean :2.387 Mean :2.407 Mean :2.407 Mean :2.420
## 3rd Qu.:2.454 3rd Qu.:2.475 3rd Qu.:2.477 3rd Qu.:2.489
## Max. :2.788 Max. :2.811 Max. :2.832 Max. :3.060
## y_pred[16] y_pred[17] y_pred[18] y_pred[19]
## Min. :2.016 Min. :2.064 Min. :1.924 Min. :1.986
## 1st Qu.:2.409 1st Qu.:2.409 1st Qu.:2.434 1st Qu.:2.433
## Median :2.478 Median :2.477 Median :2.501 Median :2.501
## Mean :2.477 Mean :2.476 Mean :2.502 Mean :2.500
## 3rd Qu.:2.546 3rd Qu.:2.544 3rd Qu.:2.569 3rd Qu.:2.568
## Max. :3.017 Max. :3.005 Max. :2.984 Max. :2.974
## y_pred[20] y_pred[21] y_pred[22] y_pred[23]
## Min. :1.994 Min. :2.052 Min. :2.039 Min. :2.004
## 1st Qu.:2.460 1st Qu.:2.476 1st Qu.:2.478 1st Qu.:2.490
## Median :2.530 Median :2.543 Median :2.545 Median :2.559
## Mean :2.529 Mean :2.544 Mean :2.546 Mean :2.559
## 3rd Qu.:2.598 3rd Qu.:2.612 3rd Qu.:2.616 3rd Qu.:2.629
## Max. :3.009 Max. :3.058 Max. :2.993 Max. :3.028
## y_pred[24] y_pred[25] y_pred[26] y_pred[27]
## Min. :2.174 Min. :2.152 Min. :2.195 Min. :2.237
## 1st Qu.:2.500 1st Qu.:2.540 1st Qu.:2.567 1st Qu.:2.568
## Median :2.567 Median :2.612 Median :2.640 Median :2.643
## Mean :2.567 Mean :2.611 Mean :2.640 Mean :2.645
## 3rd Qu.:2.634 3rd Qu.:2.683 3rd Qu.:2.713 3rd Qu.:2.722
## Max. :3.004 Max. :3.015 Max. :3.080 Max. :3.127summary(samps2[,1])## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 2.436 2.613 2.656 2.661 2.702 3.063Diagnóstico de Gelman & Rubin: através do pacote coda
- Relaciona as variâncias entre as cadeias e dentro das cadeias,
- O valor da estatística deve ser entre 1.0 e 1.2 para indicar convergência.
- Saiba mais sobre diagnósticos de Gelman & Rubin .
- Saiba mais sobre Métodos MCMC .
- Verificamos que houve convergência das cadeias
Referências na literatura:
[1] Gelman, A and Rubin, DB (1992) Inference from iterative simulation using multiple sequences, Statistical Science, 7, 457-511.
[2] Brooks, SP. and Gelman, A. (1998) General methods for monitoring convergence of iterative simulations. Journal of Computational and Graphical Statistics, 7, 434-455.
gelman.diag(samps, confidence = 0.95, transform=FALSE)## Potential scale reduction factors:
##
## Point est. Upper C.I.
## alpha 1 1.01
## beta 1 1.00
## gamma 1 1.01
## sigma2 1 1.00
## y_new 1 1.00
## y_pred[1] 1 1.00
## y_pred[2] 1 1.00
## y_pred[3] 1 1.00
## y_pred[4] 1 1.01
## y_pred[5] 1 1.00
## y_pred[6] 1 1.00
## y_pred[7] 1 1.00
## y_pred[8] 1 1.00
## y_pred[9] 1 1.00
## y_pred[10] 1 1.00
## y_pred[11] 1 1.00
## y_pred[12] 1 1.00
## y_pred[13] 1 1.01
## y_pred[14] 1 1.01
## y_pred[15] 1 1.00
## y_pred[16] 1 1.00
## y_pred[17] 1 1.02
## y_pred[18] 1 1.00
## y_pred[19] 1 1.00
## y_pred[20] 1 1.01
## y_pred[21] 1 1.00
## y_pred[22] 1 1.01
## y_pred[23] 1 1.00
## y_pred[24] 1 1.00
## y_pred[25] 1 1.00
## y_pred[26] 1 1.02
## y_pred[27] 1 1.01
##
## Multivariate psrf
##
## 1.01gelman.plot(samps)
Traços e densidades através do pacote bayesplot
- Problemas em mostrar os traços & densidade pelas funções mcmc_trace & mcmc_dens
mcmc_trace(samps)
mcmc_dens(samps)
Recorrendo ao plot para soltar os traços e densidades:
plot(samps)
Gráficos de autocorrelação do pacote coda
- os gráficos em preto são referentes a primeira cadeia
- os gráficos em vermelho são referentes a segunda cadeia
autocorr.plot(samps[[1]],col=1)
autocorr.plot(samps[[2]],col=2)
Gráfico de dispersão com observados & esperados:
sumarios=colMeans(samps2)
alpha=round(sumarios[1],4)
beta=round(sumarios[2],4)
gamma=round(sumarios[3],4)
preditos=as.numeric(sumarios[6:32])
temp = dugongues %>%
mutate(preditos=preditos)
datatable(temp)#para fazer anotação no gráfico
grob <- grobTree(textGrob(paste0("curva de regressão não linear: f(x)=",alpha,"-",beta,"*",gamma,"^","x"), x=0.2, y=0.1, hjust=0,
gp=gpar(col="red", fontsize=13, fontface="italic")))
temp %>%
ggplot(aes(x=x,y=preditos)) +
geom_line()+
geom_point(aes(x=x,y=y),color="blue")+
annotation_custom(grob)
Conclusões: Com a inferência Bayesiana e o pacote Rjags, é possível estabelecer as quantidades relacionadas ao tamanho dos dugongues em função da idade. Verificou-se convergência do algoritmo.
Pendências:
- Verificar se existe uma forma de estabelecer semente - igual set.seed(), testei nesta análise com rjags mas não funciona.
- Segundo este fórum precisa “setar” a semente de outra forma https://sourceforge.net/p/mcmc-jags/discussion/610037/thread/c6b2ee2b/
- Seria interessante agregar os sumários a posteriori na tabela dos dados de grama;
- Colocar uma tabela com as estatísticas do linear model (lm) junto às estimativas Bayesianas.