Obtendo distribuições posteriores conjuntas via INLA

1 Introdução

Este documento tem a finalidade de exemplificar a obtenção de distribuições posteriores conjuntas via “Integrated Nested Laplace Approximation” (INLA), utilizando o pacote de mesmo nome no R.

Como exemplo de motivação é utilizada uma regressão logística, onde o objetivo é estimar, além dos parâmetros \(\beta_0\) e \(\beta_1\), uma outra medida que é uma função linear destes dois parâmetros.

Do ponto de vista frequentista, a estimativa pontual de uma função linear de parâmetros é obtida de maneira direta. No entanto, para se obter uma estimativa do erro-padrão dessa estimativa, é necessário recorrer à teoria assintótica para se obter a distribuição dessa função.

No caso bayesiano, a distribuição conjunta dos parâmetros que formam essa função pode ser utilizada de imediato para se obter tanto estimativas pontuais quanto medidas de variação.

A obtenção da distribuição conjunta dos parâmetros de um modelo pode ser obtida facilmente quando se realiza a inferância bayesiana via métodos de simulação como MCMC. No entanto, quando a inferência é realizada através de métodos de aproximação, como é o caso do INLA, a obtenção de posteriores conjuntas não é direta.

Portanto, aqui serão mostrados três métodos de obtenção de estimativas pontuais e de variação para os parâmetros de uma regressão logística e para uma função linear destes parâmetros: método da máxima verossimilhança, método bayesiano via MCMC e método bayesiano via INLA.

2 Dados

Os dados utilizados aqui como exemplo são do livro de Kinas and Andrade (2010), capítulo 8, e são referentes ao comprimento e maturidade de fêmeas do peixe-galo (Zenopsis conchifera).

Classe <- paste0("[", seq(10, 50, 5), ",", seq(15, 55, 5), ")")
Total <- c(97, 89, 73, 69, 62, 55, 59, 60, 40)
Maturas <- c(4, 2, 9, 52, 62, 53, 59, 60, 40)
dados <- data.frame(Classe, Total, Maturas)

Table 2.1: Número total e de fêmeas maturas por classe de comprimento (cm)

Classe	Total	Maturas
[10,15)	97	4
[15,20)	89	2
[20,25)	73	9
[25,30)	69	52
[30,35)	62	62
[35,40)	55	53
[40,45)	59	59
[45,50)	60	60
[50,55)	40	40

Com estes dados podemos calcular então a proporção de fêmeas maturas (Prop) e definir o ponto médio de cada classe (PM) como o valor de comprimento. Além disso, vamos centralizar o ponto médio (PMc) para facilidade de interpretação. A tabela de dados completa fica dessa forma:

Prop <- Maturas/Total
PM <- seq(10, 50, 5) + 2.5
PMc <- PM - mean(PM)
dados <- cbind(dados, Prop, PM, PMc)

Table 2.2: Tabela de dados completa

Classe	Total	Maturas	Prop	PM	PMc
[10,15)	97	4	0.0412371	12.5	-20
[15,20)	89	2	0.0224719	17.5	-15
[20,25)	73	9	0.1232877	22.5	-10
[25,30)	69	52	0.7536232	27.5	-5
[30,35)	62	62	1.0000000	32.5	0
[35,40)	55	53	0.9636364	37.5	5
[40,45)	59	59	1.0000000	42.5	10
[45,50)	60	60	1.0000000	47.5	15
[50,55)	40	40	1.0000000	52.5	20

Através do gráfico ao lado, podemos verificar que é razoável a suposição de que a proporção de fêmeas maturas aumenta conforme o tamanho.

library(ggplot2)
ggplot(dados, aes(PM, Prop)) + geom_line() +
    xlab("Comprimento (cm)") +
    ylab("Proporção de fêmeas maturas")

Figure 2.1: Proporção de fêmeas maturas por comprimento.

Normalmente, existe o interesse em se determinar o “tamanho de primeira maturação” das fêmeas, que geralmente é o comprimento em que 50% das fêmeas estão maturas (LT-50). A inspeção visual do gráfico já da uma ideia de qual seria esse comprimento. No entanto, precisamos modelar esses dados e obter uma forma mais robusta de se obter essa estimativa.

3 Modelo

Seja \(p_i\) a probabilidade de uma fêmea estar matura na classe de comprimento \(i\). Temos ao total \(n_i\) fêmeas, das quais \(y_i\) são maturas. Dessa forma, é natural que a variável aleatória \(Y_i\) (número de fêmeas maturas) seja considerada como proveniente de uma distribuição binomial, ou seja \(Y_i \sim \text{Bin}(n_i, p_i)\).

Com isso, podemos assumir um modelo logístico para relacionar a proporção de fêmeas maturas com o comprimento, através das seguinte relação:

\[ \log \left( \frac{p_i}{1 - p_i} \right) = \beta_0 + \beta_1 (x_i - \bar{x}) \]

onde a função \(\log(\frac{p_i}{1-p_i})\) é chamada de logit e é a função de ligação canônica entre a variável resposta (proporção de fêmeas maturas, Prop) e o preditor linear, \(\beta_0\) e \(\beta_1\) são os parâmetros a serem estimados, e \((x_i - \bar{x})\) é a variável resposta centrada (ponto médio das classes de comprimento, PMc).

Para encontrarmos o LT-50, basta lembrar que esse é o comprimento em que 50% das fêmeas estão maturas. Portanto, fazemos \(p_i = 0.5\) na equação acima e resolvemos para \(x_i\), que agora chamaremos de LT-50. Portanto, temos que o LT-50 pode ser estimado com:

\[\begin{align*} \log \left( \frac{0.5}{1 - 0.5} \right) &= \beta_0 + \beta_1 (\text{LT-50} - \bar{x}) \\ 0 &= \beta_0 + \beta_1 (\text{LT-50} - \bar{x}) \\ -\frac{\beta_0}{\beta_1} &= \text{LT-50} - \bar{x} \\ & \Rightarrow \text{LT-50} = -\frac{\beta_0}{\beta_1} + \bar{x} \end{align*}\]

4 Métodos de estimação

4.1 Método da máxima verossimilhança

As estimativas pontuais dos parâmetros do modelo pelo método da máxima verossimilhança serão obtidas através de um GLM, especificando a família de distribuição da variável resposta como binomial e função de ligação logit.

m.glm <- glm(Prop ~ PMc, data = dados, weights = Total,
             family = binomial(link = "logit"))
summary(m.glm)
# 
# Call:
# glm(formula = Prop ~ PMc, family = binomial(link = "logit"), 
#     data = dados, weights = Total)
# 
# Deviance Residuals: 
#     Min       1Q   Median       3Q      Max  
# -2.3136  -0.7782   0.1192   1.0183   3.1793  
# 
# Coefficients:
#             Estimate Std. Error z value Pr(>|z|)    
# (Intercept)  2.89039    0.33099   8.733   <2e-16 ***
# PMc          0.41010    0.03839  10.681   <2e-16 ***
# ---
# Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 
# (Dispersion parameter for binomial family taken to be 1)
# 
#     Null deviance: 625.993  on 8  degrees of freedom
# Residual deviance:  27.295  on 7  degrees of freedom
# AIC: 48.008
# 
# Number of Fisher Scoring iterations: 6

Através deste resultado, obtemos as estimativas pontuais e o erro-padrão para cada parâmetro. Com isso, já temos ideia da precisão dessas estimativas.

Para calcular o LT-50, basta então aplicar a equação já obtida acima.

## Estimativas pontuais
b0.glm <- coef(m.glm)[1]
b1.glm <- coef(m.glm)[2]

## Cálculo do LT50
(lt50.glm <- -(b0.glm/b1.glm) + mean(dados$PM))
# (Intercept) 
#    25.45201

Assim, obtemos uma estimativa do LT-50 como sendo 25.45. No entanto, não sabemos a precisão dessa estimativa pois para calcular o erro-padrão associado precisamos da distribuição assintótica do LT-50.

Mas antes de obtê-la, vamos analisar a distribuição assintótica de \(\beta_0\) e \(\beta_1\).

Mostrar que asintoticamente cada parâmetro tem distribuição normal com média \(\beta\) e variância x. E que a matriz de variância-covariância das estimativas é a matriz de informação de Fisher…

## Variâncias dos parâmetros
varcov <- vcov(m.glm)
var.b0 <- vcov(m.glm)[1,1]
var.b1 <- vcov(m.glm)[2,2]

## Distribuição assintótica de b0
x.b0 <- seq(-3, 3, length = 100) + b0.glm
dens.b0 <- dnorm(x.b0, mean = b0.glm, sd = sqrt(var.b0))
## Distribuição assintótica de b1
x.b1 <- seq(-2, 2, length = 100) + b1.glm
dens.b1 <- dnorm(x.b1, mean = b1.glm, sd = sqrt(var.b1))

Gráficos com as distribuiçãos assintoticas dos parâmetros do modelo logístico 4.1.

library(gridExtra)
p1 <- qplot(x.b0, dens.b0, geom = "line",
            xlab = expression(beta[0]), ylab = "Densidade")
p2 <- qplot(x.b1, dens.b1, geom = "line",
            xlab = expression(beta[1]), ylab = "Densidade")
grid.arrange(p1, p2, ncol = 2)

Figure 4.1: Distribuição assintótica dos betas.

Resultados assintóticos

Dentre as várias propriedades de um estimador, duas são fundamentais: a consistência e a eficiência. A consistência se refere ao fato de que um estimador deve deve convergir para o verdadeiro valor do parâmetro, a medida que o tamanho da amostra aumenta. Já a eficiência diz respeito à variância assintótica de um estimador. Quando dois (ou mais) estimadores não viesados são comparados, o melhor (mais eficiente) é aquele que possui a menor variância possível. Também podemos dizer que um estimador é eficiente se sua variância atinge um certo limite, ou a mínima variância possível.

Os estimadores de máxima verossimilhança (EMV) possuem estas duas propriedades: são consistentes e eficientes. Isso significa que, à medida que o tamanho da amostra aumenta, a estimativa pontual fica cada vez mais próxima do verdadeiro valor, e a variância dessa estimativa é a menor possível entre todos os estimadores.

Formalizando, dizemos que \(\tau(\hat\theta)\) é um estimador consistente para \(\tau(\theta)\) se

\[ \lim_{n \to \infty} P[|\tau(\hat\theta) - \tau(\theta)| \geq \epsilon] = 0 \]

para \(\epsilon > 0\). Note que \(\tau(\cdot)\) é uma função qualquer (diferenciável) de \(\theta\).

Como mencionado anteriormente, podemos definir uma “cota” ou “limite” inferior que a variância de um estimador deve atingir para ser considerado um estimador eficiente. Esse limite inferior é o Limite Inferior de Cramér-Rao (LICR), definido por (ver Casella and Berger 2011, 299–301):

\[ Var[\tau(\hat\theta)] \geq \frac{[\tau'(\theta)]^2}{I_E(\theta)} \]

onde \(\tau'(\theta) = \frac{\partial \tau(\theta)}{\partial\theta}\), e

\[ I_E(\theta) = E \left[ \left(\frac{\partial l(\theta)}{\partial \theta} \right)^2 \right] = -E \left[ \frac{\partial^2 l(\theta)}{\partial \theta^2} \right] \]

é a Informação de Fisher esperada. Assim, se a variância de um estimador for igual ao LICR, então dizemos que ele é eficiente (e todos os EMVs são).

De acordo com o LICR, existe uma variância assintótica ótima. Com isso, podemos definir que a distribuição assintótica de \(\tau(\hat\theta)\) é

\[ \sqrt{n}[\tau(\hat\theta) - \tau(\theta)] \quad \rightarrow \quad N(0, V[\tau(\theta)]) \]

onde \(V[\tau(\theta)]\) é o LICR. Esse resultado nos diz então que o EMV \(\tau(\hat\theta)\) é um estimador consistente e assintoticamente eficiente de \(\tau(\theta)\). Também podemos escrever o mesmo resultado da seguinte forma:

\[ \tau(\hat\theta) \sim N(\tau(\theta), V[\tau(\theta)]) \]

Portanto, se um EMV é assintóticamente eficiente, sua variância assintótica é a variância do Método Delta. O Método Delta é uma generalização do Teorema do Limite Central, que diz o seguinte: se \(\hat\theta\) satisfaz \(\sqrt{n}[\hat\theta - \theta] \rightarrow N(0, \sigma^2)\), então dada uma função \(\tau(\cdot)\) diferenciável em \(\theta\), a distribuição de \(\tau(\hat\theta)\) fica (ver Casella and Berger 2011, 217)

\[ \sqrt{n}[\tau(\hat\theta) - \tau(\theta)] \quad \rightarrow \quad N(0, \sigma^2 [\tau'(\theta)]^2) \]

Portanto, podemos utilizar o LICR como uma aproximação para a variância verdadeira do EMV. Mas, note que no LICR o denominador é a informação de Fisher esperada. No entanto, nem sempre conseguimos calcular a esperança para chegar nessa informação. Sendo assim, a partir do Método Delta e da variância assintótica do EMV, podemos aproximar a variância de \(\tau(\hat\theta)\) por:

\[\begin{align*} V[\tau(\hat\theta)|\theta] &= \frac{[\tau'(\theta)]^2}{I_E(\theta)} \\ &\approx \frac{[\tau'(\theta)]^2}{I_O(\theta)} \end{align*}\]

onde

\[ I_O(\theta) = - \frac{\partial^2 l(\theta)}{\partial \theta^2} \]

é a informação de Fisher observada. Segundo Casella and Berger (2011), o uso do número de informações observado é superior ao número de informações esperado, que é o utilizado no LICR. (Até aqui é um resumo de Casella and Berger (2011), pp. 421-423). Além disso, essa aproximação tem uma consequência computacional. No processo de otimização de uma função qualquer, é possível obter a matriz Hessiana, que contém as segundas derivadas parciais (numéricas) da função em relação a cada parâmetro, que é por si própria a matriz de informação de Fisher observada. Portanto, para obter as variâncias dos estimadores, basta inverter a matriz Hessiana e selecionar os elementos de sua diagonal.

Resumindo, temos então que a distribuição assintótica de \(\tau(\hat\theta)\) é

\[ \tau(\hat\theta) \sim N \left( \tau(\theta), \frac{[\tau'(\theta)]^2}{I_O(\theta)} \right) \]

Note também que essa notação até aqui define todos os componentes de forma geral, ou seja para \(\tau(\theta)\), sendo \(\tau(\cdot)\) qualquer função diferenciável. Um caso particular (e mais simples) é quando \(\tau(\theta) = \theta\), ou seja, estamos interessados no próprio valor do parâmetro e não em uma função dele. Nessas condições, temos que \(\tau'(\theta) = 1\) e a distribuição do estimador \(\tau(\hat\theta) = \hat\theta\) fica

\[ \hat\theta \sim N \left( \theta, \frac{1}{I_O(\theta)} \right) \]

Ou seja, a variância assintótica de \(\hat\theta\) é simplesmente o inverso da informação de Fisher observada (ou o inverso da matriz Hessiana).

Outro detalhe importante é que o uso de \(I_E(\theta)\) ou \(I_O(\theta)\) geralmente leva ao mesmo resultado, pois se a segunda derivada da função em cada parâmetro não depender de \(Y\), então \(I_E(\theta) = I_O(\theta)\) (precisa de uma fonte aqui).

Voltando ao exemplo que estávamos utilizando, queremos agora chegar nos mesmos resultados obtidos anteriormente pela função lm() (que estima os parâmetros analiticamente pelo método dos mínimos quadrados), e calculados à mão. Ou seja, queremos as estimativas pontuais dos betas e da variância através da maximização da função de verossimilhança (ou minimização do log da verossimilhança), e as estimativas das variâncias destes estimadores obtidas através dos resultados assintóticos.

Para obter a distribuição assintotica do LT-50, precisamos:

1. Derivadas parciais da função \(-\frac{\beta_0}{\beta_1}\) em relação a \(\beta_0\) e \(\beta_1\)
2. Com isso calculamos a variância dessa função
3. Obtem a distribuição normal assintotica

## Distribuição assintótica do LT50 ------------------------------------

## Derivadas parciais de -(b0/b1)
## Em relação a b0
D(expression(-b0/b1), "b0")
# -(1/b1)
b0.deriv <- - (1/b1.glm)
## Em relação a b1
D(expression(-b0/b1), "b1")
# b0/b1^2
b1.deriv <- b0.glm/b1.glm^2

## Vetor de derivadas
u.theta <- cbind(b0.deriv, b1.deriv)

## Variância de -(b0/b1)
lt50.var <- as.numeric(u.theta %*% varcov %*% t(u.theta))

## Distribuição assintótica de LT50
x.lt50 <- seq(-3, 3, length = 100) + lt50.glm
dens.lt50 <- dnorm(x.lt50, mean = lt50.glm, sd = sqrt(lt50.var))

qplot(x.lt50, dens.lt50, geom = "line",
      xlab = "LT-50", ylab = "Densidade")

Figure 4.2: Distribuição assintótica do LT-50.

4.2 Inf. bayesiana via JAGS

##----------------------------------------------------------------------
## JAGS
library(runjags)

## Dados
datalist <- dump.format(list(y = dados$Maturas,
                             n = dados$Total,
                             x = dados$PMc))
params <- c("beta0", "beta1")
inicial <- dump.format(list(beta0 = b0.glm,
                            beta1 = b1.glm))

## Modelo
galomcmc <- "model{
for(i in 1:length(y)){
  logit(p[i]) <- beta0 + beta1 * x[i]
  y[i] ~ dbin(p[i], n[i])
}
beta0 ~ dnorm(0, 0.001)
beta1 ~ dnorm(0, 0.001)
}"

## Ajuste
m.jags <- run.jags(model = galomcmc, monitor = params,
                   data = datalist, inits = inicial,
                   n.chains = 3, burnin = 50000, thin = 3,
                   sample = 9000)
# Warning: Convergence cannot be assessed with only 1 chain

## Resultados
names(m.jags)
head(m.jags$mcmc)
b0.jags <- m.jags$mcmc[[1]][, 1]
b1.jags <- m.jags$mcmc[[1]][, 2]
summary(b0.jags)
summary(b1.jags)
par(mfrow = c(1,3))
plot(density(b0.jags))
plot(density(b1.jags))

## LT50
lt50.jags <- - (b0.jags/b1.jags) + mean(dados$PM)
summary(lt50.jags)
plot(density(lt50.jags))

par(mfrow = c(1,1))

4.3 Inf. bayesiana via INLA

##----------------------------------------------------------------------
## INLA
library(INLA)

## Modelo e ajuste
f <- Maturas ~ PMc
m.inla <- inla(f, family = "binomial", data = dados,
               Ntrials = dados$Total,
               control.compute = list(config = TRUE))

## Resultados
summary(m.inla)
# 
# Call:
#    c("inla(formula = f, family = \"binomial\", data = dados, Ntrials = 
#    dados$Total, ", " control.compute = list(config = TRUE))") 
# Time used:
#     Pre = 0.174, Running = 0.0398, Post = 0.0134, Total = 0.227 
# Fixed effects:
#              mean    sd 0.025quant 0.5quant 0.975quant  mode kld
# (Intercept) 2.904 0.331      2.292     2.89      3.591 2.863   0
# PMc         0.412 0.038      0.342     0.41      0.492 0.407   0
# 
# Expected number of effective parameters(stdev): 2.00(0.00)
# Number of equivalent replicates : 4.50 
# 
# Marginal log-Likelihood:  -29.57
names(m.inla)
#  [1] "names.fixed"                "summary.fixed"             
#  [3] "marginals.fixed"            "summary.lincomb"           
#  [5] "marginals.lincomb"          "size.lincomb"              
#  [7] "summary.lincomb.derived"    "marginals.lincomb.derived" 
#  [9] "size.lincomb.derived"       "mlik"                      
# [11] "cpo"                        "po"                        
# [13] "waic"                       "model.random"              
# [15] "summary.random"             "marginals.random"          
# [17] "size.random"                "summary.linear.predictor"  
# [19] "marginals.linear.predictor" "summary.fitted.values"     
# [21] "marginals.fitted.values"    "size.linear.predictor"     
# [23] "offset.linear.predictor"    "model.spde2.blc"           
# [25] "summary.spde2.blc"          "marginals.spde2.blc"       
# [27] "size.spde2.blc"             "model.spde3.blc"           
# [29] "summary.spde3.blc"          "marginals.spde3.blc"       
# [31] "size.spde3.blc"             "logfile"                   
# [33] "misc"                       "dic"                       
# [35] "mode"                       "neffp"                     
# [37] "joint.hyper"                "nhyper"                    
# [39] "version"                    "Q"                         
# [41] "graph"                      "ok"                        
# [43] "cpu.used"                   "all.hyper"                 
# [45] ".args"                      "call"                      
# [47] "model.matrix"
m.inla$marginals.fixed
# $`(Intercept)`
#                x            y
#  [1,] -0.4065977 2.067872e-29
#  [2,]  0.2553928 3.566060e-23
#  [3,]  0.9173833 6.148059e-12
#  [4,]  1.2483785 3.621000e-08
#  [5,]  1.4299908 2.003886e-06
#  [6,]  1.5454359 1.912332e-05
#  [7,]  1.5793737 3.558511e-05
#  [8,]  1.6696932 1.689339e-04
#  [9,]  1.8107451 1.456897e-03
# [10,]  1.9103690 5.462541e-03
# [11,]  1.9240775 6.468995e-03
# [12,]  1.9788145 1.232883e-02
# [13,]  2.1187039 5.168665e-02
# [14,]  2.1893420 9.518130e-02
# [15,]  2.2917255 2.028942e-01
# [16,]  2.3832994 3.532586e-01
# [17,]  2.4431758 4.782641e-01
# [18,]  2.4899999 5.872601e-01
# [19,]  2.5293444 6.834251e-01
# [20,]  2.5638603 7.688033e-01
# [21,]  2.5951611 8.451095e-01
# [22,]  2.6234687 9.116943e-01
# [23,]  2.6498830 9.705556e-01
# [24,]  2.6752624 1.023154e+00
# [25,]  2.6989708 1.068011e+00
# [26,]  2.7220372 1.107062e+00
# [27,]  2.7443201 1.139987e+00
# [28,]  2.7659749 1.167079e+00
# [29,]  2.7872801 1.188722e+00
# [30,]  2.8081404 1.204891e+00
# [31,]  2.8288535 1.215897e+00
# [32,]  2.8493679 1.221777e+00
# [33,]  2.8575545 1.222728e+00
# [34,]  2.8657412 1.222886e+00
# [35,]  2.8698346 1.222670e+00
# [36,]  2.8739309 1.222256e+00
# [37,]  2.8821484 1.220838e+00
# [38,]  2.8903660 1.218640e+00
# [39,]  2.8985836 1.215671e+00
# [40,]  2.9068011 1.211943e+00
# [41,]  2.9109099 1.209797e+00
# [42,]  2.9150493 1.207448e+00
# [43,]  2.9234246 1.202128e+00
# [44,]  2.9317999 1.196060e+00
# [45,]  2.9527382 1.177735e+00
# [46,]  2.9743623 1.154337e+00
# [47,]  2.9961195 1.126579e+00
# [48,]  3.0188428 1.093523e+00
# [49,]  3.0420306 1.056030e+00
# [50,]  3.0662025 1.013502e+00
# [51,]  3.0916069 9.657194e-01
# [52,]  3.1179214 9.136745e-01
# [53,]  3.1463701 8.554123e-01
# [54,]  3.1768919 7.916730e-01
# [55,]  3.2098673 7.226080e-01
# [56,]  3.2462631 6.475653e-01
# [57,]  3.2874977 5.657144e-01
# [58,]  3.3356165 4.763869e-01
# [59,]  3.3935437 3.798928e-01
# [60,]  3.4708479 2.720967e-01
# [61,]  3.5914171 1.508366e-01
# [62,]  3.7348971 6.720516e-02
# [63,]  3.8333249 3.614026e-02
# [64,]  3.8963404 2.362570e-02
# [65,]  4.0394602 8.298370e-03
# [66,]  4.1206606 4.360343e-03
# [67,]  4.2273356 1.772286e-03
# [68,]  4.2971664 9.504515e-04
# [69,]  4.5213016 1.074025e-04
# [70,]  4.5583308 7.295783e-05
# [71,]  4.7238784 1.181302e-05
# [72,]  4.8893261 1.648277e-06
# [73,]  4.9116810 1.248717e-06
# [74,]  5.5513165 1.387946e-10
# [75,]  6.2133070 3.859054e-16
# 
# $PMc
#                x            y
#  [1,] 0.02799081 4.508005e-31
#  [2,] 0.10478039 1.371643e-21
#  [3,] 0.18156997 1.195211e-11
#  [4,] 0.21996476 1.354145e-07
#  [5,] 0.24472410 1.876555e-05
#  [6,] 0.25709987 1.616732e-04
#  [7,] 0.25835955 1.989446e-04
#  [8,] 0.27124508 1.468177e-03
#  [9,] 0.28741668 1.316384e-02
# [10,] 0.29675434 3.987262e-02
# [11,] 0.30017823 5.817339e-02
# [12,] 0.30623929 1.094114e-01
# [13,] 0.32198368 4.552088e-01
# [14,] 0.32991624 8.336138e-01
# [15,] 0.34165896 1.789148e+00
# [16,] 0.35208464 3.109770e+00
# [17,] 0.35885564 4.194201e+00
# [18,] 0.36420477 5.145240e+00
# [19,] 0.36869732 5.980952e+00
# [20,] 0.37264230 6.721061e+00
# [21,] 0.37622523 7.381228e+00
# [22,] 0.37946896 7.955829e+00
# [23,] 0.38249758 8.462161e+00
# [24,] 0.38541199 8.913400e+00
# [25,] 0.38813152 9.296014e+00
# [26,] 0.39078574 9.628326e+00
# [27,] 0.39334510 9.906147e+00
# [28,] 0.39584158 1.013364e+01
# [29,] 0.39829288 1.031288e+01
# [30,] 0.40070148 1.044479e+01
# [31,] 0.40308956 1.053138e+01
# [32,] 0.40546112 1.057353e+01
# [33,] 0.40640715 1.057819e+01
# [34,] 0.40735319 1.057597e+01
# [35,] 0.40782621 1.057230e+01
# [36,] 0.40829923 1.056693e+01
# [37,] 0.40925031 1.055104e+01
# [38,] 0.41020224 1.052837e+01
# [39,] 0.41115416 1.049905e+01
# [40,] 0.41210608 1.046318e+01
# [41,] 0.41258204 1.044282e+01
# [42,] 0.41305801 1.042087e+01
# [43,] 0.41402956 1.037119e+01
# [44,] 0.41500218 1.031504e+01
# [45,] 0.41743372 1.014761e+01
# [46,] 0.41994219 9.936825e+00
# [47,] 0.42246397 9.689395e+00
# [48,] 0.42511539 9.394965e+00
# [49,] 0.42780911 9.064544e+00
# [50,] 0.43063546 8.689379e+00
# [51,] 0.43359463 8.271565e+00
# [52,] 0.43666305 7.818177e+00
# [53,] 0.44000001 7.309900e+00
# [54,] 0.44357012 6.757851e+00
# [55,] 0.44743379 6.161271e+00
# [56,] 0.45170571 5.514738e+00
# [57,] 0.45655226 4.811655e+00
# [58,] 0.46221110 4.047207e+00
# [59,] 0.46901897 3.225353e+00
# [60,] 0.47812522 2.308608e+00
# [61,] 0.49233337 1.280185e+00
# [62,] 0.50933285 5.688595e-01
# [63,] 0.52090559 3.068368e-01
# [64,] 0.52712307 2.154966e-01
# [65,] 0.54542613 6.972892e-02
# [66,] 0.55505314 3.653913e-02
# [67,] 0.56551786 1.736870e-02
# [68,] 0.57571608 8.073238e-03
# [69,] 0.60239278 8.972235e-04
# [70,] 0.60391265 7.850324e-04
# [71,] 0.62636097 9.821979e-05
# [72,] 0.64230744 1.989878e-05
# [73,] 0.64831364 1.062734e-05
# [74,] 0.71909701 2.250506e-09
# [75,] 0.79588659 8.619190e-15
b0.inla <- m.inla$marginals.fixed[[1]]
b1.inla <- m.inla$marginals.fixed[[2]]
par(mfrow = c(1,3))
plot(b0.inla[, 1], b0.inla[, 2], type = "l")
plot(b1.inla[, 1], b1.inla[, 2], type = "l")

## LT50
## Naive
lt50.inla.naive <- - (b0.inla[, 1]/b1.inla[, 1]) + mean(dados$PM)
## Na verdade, essa conta aqui não faz sentido nenhum, porque esses são
## apenas pontos escolhidos para fazer a curva (de uma posterior
## analítica, não simulada)
summary(lt50.inla.naive)
#    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#   24.69   25.31   25.45   25.87   25.63   47.03
plot(density(lt50.inla.naive))

par(mfrow = c(1,1))

5 Comparativos

##----------------------------------------------------------------------
## Comparativos
cbind(c(b0.glm, b1.glm, lt50.glm),
      c(mean(b0.jags), mean(b1.jags), mean(lt50.jags)),
      c(mean(b0.inla[, 1]), mean(b1.inla[, 1]), mean(lt50.inla.naive)))
#                   [,1]       [,2]       [,3]
# (Intercept)  2.8903923  2.9455219  2.9331178
# PMc          0.4101019  0.4177817  0.4157979
# (Intercept) 25.4520140 25.4524047 25.8662408

par(mfrow = c(2, 3))
plot(density(b0.jags))
plot(density(b1.jags))
plot(density(lt50.jags))
plot(b0.inla[, 1], b0.inla[, 2], type = "l")
plot(b1.inla[, 1], b1.inla[, 2], type = "l")
plot(density(lt50.inla.naive))

par(mfrow = c(1, 1))

## LT50
## AQUI ESTA O TRUQUE!
## posterior sample
## help(inla.posterior.sample)
m.inla.post <- inla.posterior.sample(n = 9000, m.inla, seed = 123)
## Ele retorna as predições para cada observação e as duas últimas
## linhas são os parâmetros
m.inla.post[[1]]$latent
#                  sample1
# Predictor:1   -6.1442145
# Predictor:2   -3.8674545
# Predictor:3   -1.5833894
# Predictor:4    0.7000430
# Predictor:5    2.9764060
# Predictor:6    5.2609576
# Predictor:7    7.5381843
# Predictor:8    9.8210664
# Predictor:9   12.0972381
# (Intercept):1  2.9895611
# PMc:1          0.4580167
betas.post <- t(sapply(m.inla.post, function(x) x$latent[10:11]))
b0.inla.post <- betas.post[, 1]
b1.inla.post <- betas.post[, 2]
summary(b0.inla.post)
#    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#   1.673   2.677   2.900   2.899   3.120   4.182
summary(b1.inla.post)
#    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#  0.2578  0.3859  0.4123  0.4118  0.4379  0.5572

lt50.inla.post <- - (b0.inla.post/b1.inla.post) + mean(dados$PM)
summary(lt50.inla.post)
#    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#   24.01   25.18   25.47   25.46   25.74   26.97

## Ver
#m.inla$misc$configs

##----------------------------------------------------------------------
## Comparativos
tab <- cbind(c(b0.glm, b1.glm, lt50.glm),
             c(mean(b0.jags), mean(b1.jags), mean(lt50.jags)),
             c(mean(b0.inla[, 1]), mean(b1.inla[, 1]), mean(lt50.inla.naive)),
             c(mean(b0.inla.post), mean(b1.inla.post), mean(lt50.inla.post)))
dimnames(tab) <- list(c("b0", "b1", "lt50"),
                      c("GLM", "JAGS", "INLA naive", "INLA sample"))
tab
#             GLM       JAGS INLA naive INLA sample
# b0    2.8903923  2.9455219  2.9331178   2.8990799
# b1    0.4101019  0.4177817  0.4157979   0.4118125
# lt50 25.4520140 25.4524047 25.8662408  25.4631084

tab.sum <- cbind(summary(as.numeric(b0.jags)),
                 summary(b0.inla.post),
                 summary(as.numeric(b1.jags)),
                 summary(b1.inla.post),
                 summary(as.numeric(lt50.jags)),
                 summary(lt50.inla.post))
colnames(tab.sum) <- c("b0 JAGS", "b0 INLA", "b1 JAGS", "b1 INLA",
                       "LT50 JAGS", "LT50 INLA")
tab.sum
#          b0 JAGS  b0 INLA   b1 JAGS   b1 INLA LT50 JAGS LT50 INLA
# Min.    1.619899 1.673084 0.2915252 0.2578092  24.00771  24.01088
# 1st Qu. 2.715195 2.676753 0.3906631 0.3858670  25.17208  25.18086
# Median  2.933627 2.899629 0.4163521 0.4122953  25.44913  25.46798
# Mean    2.945522 2.899080 0.4177817 0.4118125  25.45240  25.46311
# 3rd Qu. 3.165424 3.119573 0.4428232 0.4378938  25.72824  25.74015
# Max.    4.375082 4.182383 0.5848015 0.5571502  27.20434  26.96539

## Posteriores de b0, b1 e lt50 (JAGS, INLA naive e INLA sample)
par(mfrow = c(3, 3))
plot(density(b0.jags))
plot(density(b1.jags))
plot(density(lt50.jags))
plot(b0.inla[, 1], b0.inla[, 2], type = "l")
plot(b1.inla[, 1], b1.inla[, 2], type = "l")
plot(density(lt50.inla.naive))
plot(density(b0.inla.post), type = "l")
plot(density(b1.inla.post), type = "l")
plot(density(lt50.inla.post))

par(mfrow = c(1, 1))

## Posterior conjunta (JAGS e INLA)
par(mfrow = c(1,2))
plot(as.numeric(b0.jags), as.numeric(b1.jags),
     xlim = c(1.5, 5), ylim = c(0.2, 0.65))
plot(b0.inla.post, b1.inla.post,
     xlim = c(1.5, 5), ylim = c(0.2, 0.65))

par(mfrow = c(1,1))

Referências

Casella, George, and Roger L. Berger. 2011. Inferência Estatística. São Paulo: Cengage Learning.

Kinas, Paul Gerhard, and Humber Agrelli Andrade. 2010. Introdução à Análise Bayesiana (Com R). Porto Alegre: maisQnada.