Embed Size (px)
Citation preview
Universidade Federal do Rio de Janeiro
Uma Comparacao entre Metodos de
Aproximacoes Determinısticas e
Estocastica para Inferencia Bayesiana em
Modelos Dinamicos Lineares
Generalizados
Teresa Villanueva Caballero
2013
Uma Comparacao entre Metodos deAproximacoes Determinısticas e
Estocastica para Inferencia Bayesiana emModelos Dinamicos Lineares
Generalizados
Teresa Villanueva Caballero
Dissertacao de Mestrado submetida ao Programa de Pos-Graduacao em Estatıstica do
Instituto de Matematica da Universidade Federal do Rio de Janeiro - UFRJ, como parte dos
requisitos necessarios a obtencao do grau de Mestre em Estatıstica.
Orientadora: Mariane Branco Alves
Rio de Janeiro
Novembro 2013
ii
iv
A minha famılia, em especial ao meus pais,
Juan e Flora.
v
Agradecimentos
A minha orientadora Mariane Branco Alves, obrigada pelo tempo que disponibilizou para
nossas reunioes, pela orientacao academica e paciencia ao longo deste trabalho.
A minha famılia, que sempre acreditou em mim. Aos meus pais, Juan e Flora, pelo amor,
carinho. Aos meus irmaos(as), pelo apoio incondicional, apesar da distancia.
Ao Alexei, pelo amor, compreensao e que sempre me deu forca nos momentos de desanimo,
estresse e cujo carinho e amor foi essencial ao dar um sentido na minha vida.
A Mariana Albi de Oliveira Souza, agradeco por me fornecer a programacao utilizada no
seu relatorio tecnico. Ao Thiago Guerrera Martins, pela ajuda dada, mesmo a distancia.
A todos meus amigos, e companheiros do DME que compartilharam comigo experiencias,
momentos de dificuldade e de alegria. Em especial, Mariana, Aniel, Cristian, Pamela, Renata,
Larissa, Kelly, Carlos e Arthur.
Agradeco a todos meus professores do programa de Pos-Graduacao do DME-UFRJ, pelo
valioso conhecimento transmitido, pelas maravilhosas aulas e toda a disponibilidade para
ajudar. Em especial, aos professores Nei Rocha e Alexandra Schmidt, pelo compartilhamento
de conhecimento das aulas didaticas que eles apresentam e fazem voce por mais vontade de
continuar estudando. Ao professor Dani Gamerman, pelo valioso conhecimento transmitido
no estagio docente e pela compressao, ajuda nos momentos de dificuldade.
Agradeco as professoras Alexandra Schmidt e Glaura da Conceicao Franco, por aceitarem
participar da banca e a professora Marina Silva Paez pela posicao de suplente na banca.
Por fim, agradeco a CAPES por ter financiado e possibilitaram o prosseguimento dos meus
estudos.
vi
Resumo
Nesta dissertacao, temos interesse em apresentar procedimentos de inferencia Bayesiana
na classe de modelos dinamicos lineares generalizados. Muitas vezes, as distribuicoes
de interesse nao sao possıveis de serem obtidas analiticamente, sendo necessario utilizar
metodos de aproximacao, tais como metodos determinısticos e estocasticos. Neste contexto,
apresentamos os metodos INLA (Integrated Nested Laplace Aproximation), Linear Bayes (LB)
e Monte Carlo via cadeia de Markov (MCMC). Particularmente, objetivamos comparar estes
metodos para um modelo dinamico Poisson com dados artificiais. Os tres metodos capturam
bem o comportamento da serie de dados artificiais, mas o metodo LB difere do modelo
ajustado por MCMC e INLA, pois nestes dois ultimos metodos, pressupoe variancias fixas
no tempo e enquanto o LB varia ao longo do tempo e especificadas por meio de fator de
desconto. Outra diferenca entre os metodos deve-se, a que o LB processa informacao em
tempo real, ja INLA e MCMC produzem inferencia condicional a toda informacao disponıvel.
Alem disso, estes ultimos produzem inferencia completa para os estados, diferentemente do
LB, em que tal inferencia resume-se a primeiro e segundo momentos, do vetor de estados.
Finalmente analisamos dois conjuntos de dados reais. O primeiro, trata do efeito de
poluentes atmosfericos sobre contagem de obitos de criancas menores de cinco anos por
doencas respiratorias, na cidade de Sao Paulo, usando os modelos Poisson e Poisson
inflacionado de zeros. O segundo conjunto de dados trata de efeito de volumes diarios de
chuva sobre nıveis de poluicao. Para sua estimacao usamos os modelos Gama e Bernoulli.
Palavras-Chaves: Metodos Integrated Nested Laplace Approximation, metodos de
Monte Carlo via cadeias de Markov, Linear Bayes, modelos de espaco de estados, inferencia
Bayesiana.
vii
Abstract
In this work, are interested in presenting procedures for Bayesian inference in the class
of generalized linear dynamic models. Often, distributions of interest are not available
analytically, approximated methods are needed, such as methods deterministic and stochastic.
We present methods Integrated Nested Laplace Aproximation (INLA), Linear Bayes (LB) and
Monte Carlo Markov Chain (MCMC). Particularly, we aimed to compare these methods for
a dynamic Poisson model with simulated data. The three methods capture the behavior the
series of artificial data, but the LB method differs from the adjusted model by MCMC and
INLA, since these last two methods, presupposes fixed variances in time and while the LB
varies over time and specified by the discount factor. Another difference between the methods
is due, LB processes information in real time, already INLA and MCMC processes conditional
inference all available information, Moreover, the latter produces full inference to the states,
unlike LB, in which such inference comes down to first and second moments of the state
vector.
Finally we analyze two real data sets. The first deals with the effect of air pollutants on
count of deaths of under five children with respiratory diseases in the city of Sao Paulo,using
the Poisson model and inflated Poisson models of zeros. The second set of data deals with
effect daily volumes of rain about levels of pollution. For its estimation we use the Gamma
and Bernoulli models.
Keywords: Method Integrated Nested Laplace Approximation, methods of Monte Carlo
Markov chain, Linear Bayes, state-space modeling, Bayesian inference.
viii
Sumario
1 Introducao 1
2 Modelos Dinamicos 5
2.1 Modelos Lineares Dinamicos . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.1.1 Inferencia Para Modelos Lineares Dinamicos . . . . . . . . . . . . . 6
2.2 Modelos Lineares Generalizados . . . . . . . . . . . . . . . . . . . . . . . . 7
2.3 Modelos Dinamicos Lineares Generalizados . . . . . . . . . . . . . . . . . . 8
2.4 Metodos de Aproximacao da Posteriori em Modelos Dinamicos Lineares
Generalizados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.4.1 Variacoes do Filtro de Kalman . . . . . . . . . . . . . . . . . . . . . 9
2.4.2 Aproximacao Linear Bayes . . . . . . . . . . . . . . . . . . . . . . . 10
3 Metodos de Aproximacao MCMC e INLA 21
3.1 Metodo de Aproximacao MCMC . . . . . . . . . . . . . . . . . . . . . . . . 22
3.1.1 Amostrador de Gibbs . . . . . . . . . . . . . . . . . . . . . . . . . 23
3.1.2 Algoritmo de Metropolis Hastings . . . . . . . . . . . . . . . . . . . 25
3.2 O Metodo de Aproximacao INLA . . . . . . . . . . . . . . . . . . . . . . . 27
3.2.1 Parametrizacao adequada do vetor parametrico e exploracao da grade 29
3.2.2 Aproximacao para p(θ|y) . . . . . . . . . . . . . . . . . . . . . . . 31
3.2.3 Aproximacao para p(xi|θ,y) . . . . . . . . . . . . . . . . . . . . . . 32
3.2.4 Algoritmo INLA . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
3.3 Modelo Poisson Dinamico com dados artificiais, exemplo . . . . . . . . . . . 34
3.3.1 Prioris para os parametros fixos . . . . . . . . . . . . . . . . . . . . 35
ix
3.3.2 Prioris para variaveis gaussianas latentes . . . . . . . . . . . . . . . 36
3.3.3 Aproximacao Gaussiana para a distribuicao Condicional Completa xt 38
3.3.4 Resultados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
4 Aplicacoes com Dados Reais 48
4.1 Efeito de Monoxido de Carbono sobre Obitos de Criancas em Sao Paulo . . . 48
4.1.1 Formulacao do Modelo Proposto . . . . . . . . . . . . . . . . . . . 51
4.1.2 Inferencia Bayesiana Utilizando INLA . . . . . . . . . . . . . . . . . 55
4.1.3 Escolha do melhor Modelo . . . . . . . . . . . . . . . . . . . . . . . 59
4.1.4 Resultados para o Modelo 1: Dinamica no Nıvel . . . . . . . . . . . 61
4.2 Efeito de Chuva sobre Nıveis de material Particulado no Rio de Janeiro . . . 68
4.2.1 Descricao dos Dados . . . . . . . . . . . . . . . . . . . . . . . . . . 68
4.2.2 Modelo Gama . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
4.2.3 Resultados do modelo . . . . . . . . . . . . . . . . . . . . . . . . . 71
4.2.4 Modelo Bernoulli . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
4.2.5 Resultados do modelo . . . . . . . . . . . . . . . . . . . . . . . . . 77
5 Conclusoes 85
A Codigos Usados para Dados Artificiais 87
A.1 Codigo usado pelo metodo Linear Bayes . . . . . . . . . . . . . . . . . . . 87
A.2 Codigo usado em WinBUGS . . . . . . . . . . . . . . . . . . . . . . . . . . 90
A.3 Codigo usando a Biblioteca INLA . . . . . . . . . . . . . . . . . . . . . . . 91
B Codigo Usados aos Dados de Contagem de Obitos 96
B.1 Codigo do modelo Poisson Tradicional . . . . . . . . . . . . . . . . . . . . 97
B.2 Codigo do modelo Poisson inflacao-zeros do tipo 0 . . . . . . . . . . . . . . 98
B.3 Codigo do modelo Poisson inflacao-zeros do tipo 1 . . . . . . . . . . . . . . 98
C Codigo Referentes aos Dados de Material Particulado 100
C.1 Codigo do modelo Gama . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
C.2 Codigo do modelo Bernoulli . . . . . . . . . . . . . . . . . . . . . . . . . . 102
x
Capıtulo 1
Introducao
Nos ultimos anos, varios estudos tem evidenciado associacao significativa entre a
exposicao a poluicao do ar e eventos adversos a saude, com foco em exposicao de curto
prazo. Numerosos estudos epidemiologicos tem encontrado associacoes positivas entre os
poluentes, tais como material particulado (PM10) e monoxido de carbono (CO) e mortalidade
ou morbidade, sendo muitas destas associacoes relacionadas com nıveis de poluicao que
ultrapassam limiares de seguranca, veja por exemplo, Vedal S e J. (2003), Dominici F e
J. (2002) e Alves et al. (2010). Recentemente, a Agencia Internacional de Pesquisas sobre
o Cancer (IARC), vinculada a Organizacao Mundial da Saude (OMS), classificou a poluicao
do ar exterior como uma causa de cancer. Estes estudos sao tipicamente baseados em dados
diarios de uma regiao especıfica e perıodo de tempo e a analise e efetuada utilizando metodos
de regressao de series temporais. Caso os dados de saude estejam disponıveis apenas como
contagens diarias, no contexto epidemiologico, o modelo Poisson linear generalizado e modelos
aditivos sao o metodo padrao de analise. Pode-se ter interesse, ainda, na modelagem do nıvel
diario de certo poluente atmosferico, podendo-se utilizar para tal fim, por exemplo, um modelo
Gama linear generalizado. Um outro interesse, no mesmo contexto, pode ser a explicacao de
uma resposta binaria, como por exemplo, a ultrapassagem de um limiar de seguranca, por
um certo poluente atmosferico.
Na literatura estatıstica, muitos modelos sao construıdos sob a suposicao de normalidade
da variavel resposta. Alternativas sao necessarias para o tratamento de dados que nao
satisfacam essa restricao. Como descrito por Alves (2006), Nelder e Wedderburn (1972)
1
propuseram uma classe, denominada modelos lineares generalizados (MLG), permitindo que
a distribuicao da variavel resposta pertenca a famılia exponencial. A funcao de reposta media
relaciona-se a um preditor linear por meio de uma transformacao monotona e diferenciavel,
conhecida como funcao de ligacao. Apesar da grande flexibilizacao permitida pelos MLGs - se
comparados aos tradicionais modelos lineares - estes ainda supoem independencia da variavel
resposta sobre diferentes unidades observacionais. Por outro lado, os modelos dinamicos
lineares (MDL) West e Harrison (1997), que sao um caso particular da classe de modelos de
espaco de estados, Franco et al. (2009) pressupoem normalidade da variavel resposta, mas
tratam formalmente a autocorrelacao tıpica de dados de serie temporal, ao permitir evolucao
aos parametros que controlam o preditor linear. A evolucao desses parametros e tipicamente
descrita por relacoes estocasticas markovianas. West et al. (1985) estendem tanto os MLGs
quanto os MDLs, ao combinar uma estrutura observacional nao necessariamente Gaussiana -
mais especificamente, pertencente a famılia exponencial, como nos MLGs - a uma estrutura
de evolucao dinamica para os parametros, como nos MDLs.
Do ponto de vista de realizacao de inferencia bayesiana, a classe de modelos dinamicos
lineares generalizados apresenta dificuldades, pois nao e possıvel a obtencao analıtica de
distribuicao a posteriori de todas as quantidades latentes de interesse. Ha varias alternativas
na literatura para aproximacao ou resumo da distribuicao a posteriori nesse contexto. Na
decada de 80, eram primordialmente utilizadas aproximacoes determinısticas, que se tornam
mais complexas a medida em que a dimensao do vetor parametrico aumenta. West et al.
(1985) propoem a metodologia linear bayes para realizacao de inferencia em MDLGs. Sua
abordagem baseia-se na adocao de uma distribuicao a priori conjugada para a resposta media
e a avaliacao incompleta das distribuicoes a priori e a posteriori do vetor de estados, apenas em
termos de momentos de primeira e segunda ordens, evitando assim esforco computacional para
integracao ou otimizacao. Alem do reduzido tempo computacional, permitindo realizacao de
inferencia em tempo real, outra vantagem do metodo e a obtencao de distribuicoes preditivas
com forma analıtica fechada, devido a propriedades de conjugacao da famılia exponencial
(Migon e Gamerman 1999, pp 62-70). A perda em relacao a metodos que exigem maior
esforco computacional, como MCMC, reside no fato de nao se obter a distribuicao a posteriori
para o campo latente de forma completa, mas apenas sua media e matriz de covariancia, o
2
que permite a obtencao de estimativas intervalares. No metodo linear Bayes, todo o ciclo
de inferencia e baseado no conhecimento de hiperparametros ou na sua especificacao por
meio de alguma racionalizacao, como por exemplo o uso de fatores de desconto (West e
Harrison 1997, pp 193-202), para especificacao de variancias/covariancias evolucionais. A
especificacao desses fatores pode nao ser trivial.
Fahrmeir (1992), por outro lado, apresenta uma generalizacao do filtro de Kalman
estendido em modelos dinamicos lineares generalizados multivariadas, para estimar os
parametros de estado atraves de modas a posteriori.
A partir da decada de 90, com avancos computacionais, metodos aproximados baseados
em simulacao - em particular os metodos de Monte Carlo via cadeias de Markov (MCMC),
detalhados em Gamerman e Lopes (2006) - dominaram o cenario de aproximacoes para
distribuicoes a posteriori, nos casos em que estas sao analiticamente intrataveis. Tais metodos
buscam, a partir de nucleos de transicao convenientes, a construcao iterativa de uma cadeia
de Markov homogenea, irredutıvel, ergodica, que tenha como distribuicao estacionaria a
posteriori de interesse. No caso dos modelos de espaco de estados, que pressupoem correlacao
temporal entre seus parametros, a convergencia de metodos MCMC para a distribuicao
estacionaria pode ser bastante lenta.
Devido ao elevado custo computacional dos metodos MCMC no contexto abordado, busca-
se alternativas a estes, de forma a tornar a realizacao de inferencia bayesiana aproximada mais
rapida e eficiente. Tem despertado grande interesse o trabalho de Rue et al. (2009), propondo
a realizacao de inferencia bayesiana por meio de aproximacoes determinısticas para modelos
de espaco de estados com campos latentes Gaussianos, ou seja, aqueles em que se supoe que a
evolucao estocastica dos parametros de estado e ditada por uma distribuicao Gaussiana (mas
a resposta, nao necessariamente). Os autores relatam a obtencao de estimativas acuradas
de hiperparametros e do campo latente, com tempos computacionais bastante reduzidos,
em comparacao a longas cadeias obtidas via MCMC. Resende (2011) propoe uma extensao
desse metodo para modelos de espaco de estados com campos latentes nao Gaussianos,
apresentando a base teorica da proposta, entretanto relata problemas computacionais que
impediram a exemplificacao do metodo ali proposto.
No presente trabalho, propomos uma comparacao do metodo Linear Bayes, MCMC e
3
aproximacao determinıstica, como sugerem Rue et al. (2009), a duas aplicacoes de interesse
no contexto de epidemiologia ambiental.
Este documento esta organizado da seguinte forma: no capıtulo 2, e apresentada
a estrutura dos modelos dinamicos lineares generalizados e discute-se as dificuldades
relacionadas a sua estimacao, sob abordagem bayesiana e os metodos de aproximacao da
distribuicao a posteriori em MDLG, como o Linear Bayes. Ja os metodos MCMC e INLA sao
descritos no capıtulo 3, isto devido a importancia neste trabalho, apresentando um exemplo
de aplicao a dados Poisson artificialmente gerados, com base em um preditor estruturado em
termos de um nıvel e uma covariavel (CO) com efeito dinamico. As estimativas obtidas
via INLA sao comparadas aquelas obtidas com os metodos LB e MCMC. No capıtulo
4, apresentam-se dois conjuntos de dados reais com a metodologia descrita. Na secao
4.1 apresenta-se um modelo de regressao dinamica Poisson, em que se busca descrever o
impacto de poluentes atmosfericos e variaveis climaticas sobre desfechos epidemiologicos, com
diferentes estruturas preditivas. Na secao 4.2 e apresentado um modelo Gama para quantificar
o efeito cumulativo de volumes diarios de chuva sobre o nıveis de material particulado e em
seguida o modelo com resposta Bernoulli, para analisar fatores associados a ultrapassagem de
um limiar de seguranca no nıvel de material particulado. O capıtulo 5 conclui este trabalho.
4
Capıtulo 2
Modelos Dinamicos
2.1 Modelos Lineares Dinamicos
Na literatura bayesiana, os Modelos Lineares Dinamicos (MLD) sao conhecidos tambem
como modelos de espaco de estados. Tais modelos foram introduzidos por Harrison e
Stevens (1976), estao bem documentados em West e Harrison (1997) e constituem uma
ampla classe parametrica, com parametros variando no tempo, em que tanto a variacao dos
parametros quanto a informacao a respeito de quantidades observaveis sao descritas de uma
forma probabilıstica. Os modelos lineares dinamicos possuem estrutura hierarquica e sao uma
metodologia flexıvel para tratar problemas em analises de series temporais, caracterizando-se
atraves das seguintes equacoes:
Yt = F′txt + vt, vt ∼ N(0,Vt) (2.1a)
xt = Gtxt−1 + ωt, ωt ∼ N(0,Wt), , (2.1b)
em que para t = 1, . . . , Y t = (y1, y2, . . . , yn) e o vetor de observacoes; xt e um vetor
p−dimensional denominado vetor de estados; F t e uma matriz de p×n de variaveis regressoras
ou variaveis explicativas, cujos elementos sao conhecidos; Gt e uma matriz quadrada de
ordem p que descreve a evolucao dos parametros de estado no tempo. As matrizes de
covariancia Vt e Wt, de ordem n e p, estao associadas ao erro observacional vt e ao erro
de evolucao dos estados ωt, respectivamente. Assume-se que os erros vt e ωt, sejam serial e
mutuamente independentes. O modelo completa-se com uma densidade a priori (x1|D0) ∼
5
N(a,R), em que D0 denota a informacao inicial disponıvel ao analista. O modelo descrito em
(2.1) e completamente especificado pela quadrupla {F t,Gt,Vt,Wt} e de uma distribucao
a priori assumida para os parametros de estados. A equacao (2.1a) e denominada equacao
da observacao e relaciona o vetor de observacoes e componentes estruturais (como nıvel,
tendencia, sazonalidade etc.), tendo a forma de uma regressao multivariada e a equacao
(2.1b) e denominada equacao de estados ou do sistema, responsavel pela evolucao do vetor
de coeficientes de regressao (ou parametros de estado) ao longo do tempo.
De acordo com Migon et al. (2005), modelos dinamicos podem ser vistos como uma
generalizacao de modelos de regressao, permitindo alteracoes nos valores de parametros ao
longo do tempo, por meio da introducao de uma equacao que rege a evolucao temporal da
coeffcientes regressao.
2.1.1 Inferencia Para Modelos Lineares Dinamicos
Seja Y o vetor de obsevacoes e x o vetor de parametros. De acordo com o paradigma
bayesiano, assume-se uma distribuicao a priori p(x)1 , a qual representa a incerteza inicial
acerca do vetor de parametros, antes de que Y seja observado, e a funcao de verossimilhanca
do modelo, p(Y |x). A especificacao de p(x) e p(Y |x) fornece um modelo probabilıstico,
p(Y ,x) = p(Y |x)p(x).
Tendo observado os dados Y que contem informacao acerca de x, pode-se usar Y para
atualizar a informacao acerca de x. Atraves do teorema de Bayes, encontra-se a distribucao
a posteriori de x, que contem toda informacao probabilıstica de interesse sobre x, dada por
p(x|Y ) =p(Y |x)p(x)∫p(Y |x)p(x)dx
(2.2)
Em modelos lineares dinamicos, a inferencia segue os passos usuais em inferencia bayesiana
e e realizada en forma sequencial, combinando duas operacoes principais: evolucao para
construir, a cada instante, a priori e atualizacao, para incorporar a nova observacao no tempo
t. Seja Dt = Dt−1 ∪ yt a informacao disponıvel no instante t. Entao, para cada tempo t, a
distribuicao a priori, preditiva a um passo e posteriori sao, respectivamente:
1De fato, p(x|D0), mas para simplicidade de notacao, omitimos o condicionamento no conjunto inicial de
informacao, D0.
6
p(xt|Dt−1) =
∫p(xt|xt−1)p(xt−1|Dt−1)dxt−1 (2.3)
p(yt|Dt−1) =
∫p(yt|xt)p(xt|Dt−1)dxt (2.4)
p(xt|Dt) ∝ p(yt|xt)p(xt|Dt−1), (2.5)
sendo a equacao (2.5) obtida via teorema de Bayes. Essa forma simplificada do teorema de
Bayes sera util em problemas que envolvam estimacao de parametros, ja que o denominador e
apenas uma constante normalizadora, as vezes facilmente obtida. Isso ocorre, em particular,
no caso em que (F ,G,V,W) sao todos conhecidos e assumido-se normalidade dos erros.
O algoritmo resultante, neste caso, e conhecido como filtro de Kalman (Anderson e Moore,
1979).
Em geral, o medelo descrito em (2.1) e completamente especificado pela quadrupla
{Ft, Gt, Vt,Wt} e de uma distribucao a priori assumida para os parametros de estados. Mas,
geralmente, Vt,Wt e em alguns casos elementos de Ft e Gt nao sao conhecidos, o que implica
que a inferencia nao pode ser feita de forma analıtica. Estas quantidades desconhecidas sao
chamadas de hiperparametros.
2.2 Modelos Lineares Generalizados
A classe dos modelos lineares tem por objetivo analisar a influencia de covariaveis em uma
determinada variavel resposta atraves de uma relacao linear nos parametros que governam os
impactos de tais regressoras. Uma suposicao usual, porem frequentemente inadequada, e a
de que as variaveis resposta a serem modeladas seguem distribuicao Normal.
Uma extensao dos modelos lineares permite modelar observacoes descritas por membros da
famılia exponencial. Esta classe de modelos e conhecida como Modelos Lineares Generalizados
(MLG), introduzida por Nelder e Wedderburn (1972). A ideia basica consiste em ampliar a
gama de opcoes para a distribuicao da variavel resposta, sendo a mesma pertencente a famılia
exponencial, por um conjunto de covariaveis independentes, as quais e aplicada uma estrutura
linear e dar flexibilidade para a relacao funcional entre a media da variavel resposta e o preditor
7
linear. A media passa a ser relacionada a um preditor linear apos passar por uma transformacao
monotona e diferenciavel, denominada funcao de ligacao g.
Considere-se Yt, para cada t = 1, . . . , T , a variavel resposta e F t o vetor de covariaveis
no instante t. A estrutura para o modelo linear generalizado univariado e dada por:
p(yt|ηt, φ) = exp[φ {ytηt − a(ηt)}]b(yt, φ), (2.6)
em que ηt e o parametro natural da distribucao de yt, satisfazendo
E[Yt|ηt, φ] = µt = a(ηt) (2.7a)
V [Yt|ηt, φ] = a(ηt)/φ (2.7b)
e φ e denominado parametro de escala. Um modelo linear generalizado e composto pela
estrutura observacional (2.6), combinada a um preditor linear λt, determinado por um vetor
(p× 1) de regressoras conhecidas Ft:
g(ηt) = λt = F′tx, (2.8)
sendo x um vetor latente de parametros a estimar, de ordem (p× 1), e g(.) uma funcao
de ligacao monotona e diferenciavel. O modelo completa-se com a hipotese de que os Y ′t s,
condicionalmente a ηt, t = 1, . . . , T e φ, sejam independentes e identicamente distribuıdos.
2.3 Modelos Dinamicos Lineares Generalizados
West et al. (1985) formalizaram uma extensao dos modelos lineares dinamicos (MLD) para
observacoes que pertencam a famılia exponencial, baseados no modelo linear generalizado de
Nelder e Wedderburn (1972), fazendo possıvel a utilizacao destes modelos para variados tipos
de problemas.
Os modelo dinamicos lineares generalizados (MDLG) contornam a restricao gaussiana
do modelo linear e atribuem tratamento formal a autocorrelacao serial, ao substutituir a
especificacao do preditor linear em termos de quantidades latentes x estaticas, como em
(2.8), pela dinamica:
g(ηt) = λt = F′txt, (2.9)
8
adicionando ainda a formulacao do modelo uma equacao de evolucao ou sistema, descrevendo
a forma de atualizacao do campo latente dinamico:
xt = Gtxt−1 + ωt, ωt ∼ N [0,Wt], (2.10)
em que Gt e uma matriz de transicao, suposta conhecida e de ordem (p× p) e Wt a matriz
de covariancias de ordem (p× p), associada aos erros de evolucao, ωt, dos estados ou campo
latente xt.
No caso de MLGD as integrais em (2.3),(2.4) e (2.5) nao podem ser obtidas
analiticamente, e assim a inferencia nao pode ser feita de forma exata. Muitas propostas
para resolver este problema tem sido apresentadas na literatura. Nas subsecoes seguintes
apresentam-se algumas delas.
2.4 Metodos de Aproximacao da Posteriori em Modelos
Dinamicos Lineares Generalizados
Modelos dinamicos introduzidos na secao 2.1.1 permitem a inferencia completa apenas
quando o F t, Gt e W t sao totalmente conhecidas e, ainda, sob suposicao de normalidade
dos erros. Em geral, quando estas quantidades ou outras quantidades sao desconhecidas
(hiperparametros) e a inferencia sobre eles devem basear-se na distribuicao a posteriori, essa
distribuicao nao tem solucao analıtica.
No que segue apresentamos uma revisao de alguns metodos adotados para aproximacoes
de distribuicoes a posteriori para os MLGD.
2.4.1 Variacoes do Filtro de Kalman
Fahrmeir (1992) apresenta uma generalizacao do filtro de Kalman estendido em modelos
dinamicos lineares generalizados multivariado, para estimar os parametros de estado atraves
da moda a posteriori.
O algoritmo e aplicado sequencialmente e proporciona uma aproximacao da moda a
posteriori. A utilizacao do estimador da moda a posteriori, e apenas para evitar a integracao
9
numerica. Para estimar o hiperparametro e proposto um procedimento baseado em um
algoritmo tipo-EM, Junger (2002).
Para estimar os parametros de espaco de estado, Singh e Roberts (1992) propuseram
uma aplicacao iterativa do filtro Kalman linear a modelos dinamicos lineares generalizados,
modificando a equacao observacional (2.6) por:
yt = F′txt + vt, vt ∼ N(0, Vt), (2.11)
em que yt sao observacoes modificadas, dadas por uma aproximacao linear das observacoes,
segundo:
yt = ηt + (yt − µt)g(µt) (2.12)
e com variancias associadas:
Vt = Vt(xt) =a
φt(ηt)[g
′(µt)]2 (2.13)
com g e a indicando a primeira e a segunda derivada das funcoes g e a, respectivamente.
Estas observacoes e variancias modificadas sao definidas a cada iteracao usando os valores de
xt estimados em iteracoes anteriores, pelo filtro de Kalman. Singh e Roberts (1992) estimam
W t = W utilizando uma abordagem baseada em momentos.
Fahrmeir (1997) tambem trabalhou na obtencao da moda a posteriori dos parametros
de estado para MLGD multivariado. Eles mostraram que o algoritmo proposto por Singh
e Roberts (1992) leva a moda a posteriori dos parametros de estado condicionado em um
valor fixo W . Eles tambem mostraram que a generalizacao do filtro de Kalman estendido de
Fahrmeir (1992) e um caso especial deste algoritmo com apenas uma iteracao e uma escolha
conveniente dos valores iniciais. Eles sugerem a utilizacao de um procedimento com base no
criterio de validacao cruzada generalizada para estimar hiperparametros.
2.4.2 Aproximacao Linear Bayes
West et al. (1985) propuseram uma aproximacao baseada em linear Bayes. Esta ideia foi
tambem descita por Migon e Harrison (1985) dentro do contexto de modelos nao-lineares
10
dinamicos normais e e um dos metodos aplicados neste trabalho, portanto passamos a
descreve-lo mais detalhadamente.
Suponha-se que o modelo de amostragem no tempo t tem a forma (2.6) e admita-se que
uma distribuicao a priori conjugada, denotada por (ηt|Dt−1) ∼ CP [rt, st] seja adotada para
o parametro natural ou canonico ηt
p(ηt|Dt−1) = C(rt, st)exp [rtηt − sta(ηt)] (2.14)
para algum par rt e st. A extensao dinamica natural de um modelo linear generalizado
pressupoe g(ηt) = λt = F′txt, mas tal especificacao imporia severas restricoes a priori de ηt.
Ao inves disso, West et al. (1985) utilizam a ligacao entre g(ηt) e λt apenas como um guia
para formar a priori para ηt, passando a denotar tal relacao guia por g(ηt) ≈ λt.
Ainda, suponha-se que as distribuicoes a priori e a posteriori do vetor de estados agora
nao sejam necessariamente normais, mas que, por analogia ao modelo Gaussiano, sejam
especificadas apenas pelos momentos de primeira e segunda ordens do vetor de estados xt,
dadas por
(xt−1|Dt−1) ∼ [mt−1,Ct−1] , (2.15)
(xt|Dt−1) ∼ [at, Rt] , (2.16)
em que:
at = Gtmt−1 and Rt = GtCt−1G′t +W t. (2.17)
Nesse ponto, os autores sugerem o uso de fatores de desconto para contornar o problema
de especificacao ou estimativa de W t. A ideia de fatores de desconto e especificar uma
quantidade que descreva a perda do valor de observacoes passadas para a inferencia a cada
instante. Mais especificamente, segundo West e Harrison (1997), observando-se que
V [xt−1|Dt−1] = Ct−1
e
V [xt|Dt−1] = GtCt−1G′t +W t,
denotando-se a primeira parcela no lado direito da equacao acima por P t, tem-se
Rt = P t +W t (2.18)
11
e W t, portanto, representa a inflacao na incerteza ao se passar do instante t− 1 ao instante
t, condicionalmente ao mesmo conjunto de informacao, Dt−1. Assim, tomando-se 0 < δ ≤ 1,
tal inflacao poderia ser representada por
Rt =P t
δ. (2.19)
Igualando-se (2.18) e (2.19), tem-se
W t =1− δδP t.
Portanto, condicional a Pt e arbitrando-se δ, W t fica completamente especificada. West e
Harrison (1997)[pp. 196-8] estendem essa ideia, permitindo especificacao de diferentes fatores
de desconto δj para cada bloco estrutural em um preditor, permitindo trajetorias mais suaves
(δj ≈ 1) ou mais volateis, como em West et al. (1985), que definem uma matriz diagonal
Bt, de dimensao p× p, cujos elementos sao 1√δj
, 0 < δj ≤ 1, j = 1, . . . , p.
Logo, reescrevendo (2.17), temos
at = Gtmt−1 and Rt = BtGtCt−1G′tBt, (2.20)
com a matriz de transicao Gt e matriz de descontos Bt conhecidos. Note-se que a
representacao (2.10) pode ser utilizada, mas evidentemente, ωt nao e necessariamente
normal. Alem disso,a distribuicao completa do vetor de estados nao e especificada; apenas
a media e matriz de covariancia sao assumidas. Finalmente, a distribuicao a priori para
g(ηt) = λt = F′txt e dada por
λt|Dt−1 ∼ [ft, qt], (2.21)
em que
ft = E [λt|Dt−1] = F ′tat ,
qt = V [λt|Dt−1] = F ′tRtF t ,
e
St = C [λt,xt|Dt−1] = RtF t .
Neste ponto, a priori para o parametro natural ηt esta apenas parcialmente especificada,
tendo a forma (2.14), sem qualquer restricao sobre os valores de rt e st. Estes valores sao
12
escolhidos com base na relacao g(ηt) ≈ λt, que fixa os dois primeiros momentos de g(ηt)
determinando rt e st. A relacao guia sugere os valores de ft e qt, para estes momentos e St,
para a covariancia entre g(ηt) e xt.
Com base nessa especificacao, e imediato que a distribuicao preditiva tem a forma:
p(yt|Dt−1, φ) =c(rt, st)
c(rt + φyt, st + φ)b(yt, φ), (2.22)
podendo ser diretamente obtida, e a distribucao a posteriori para (ηt|Dt) e a congujada
atualizada da forma ηt|Dt ∼ CP (rt + φyt, st + φ). Uma analise bayesiana completa requer
tambem a posteriori para (xt|Dt), mas esta nao esta disponıvel porque a priori para (xt|Dt−1)
e apenas parcialmente especificada e o modelo nao fornece verossimilhanca para xt. O modelo
desenvolvido ate agora, no entanto, nao requer a especificacao completa para prosseguir para
o tempo (t + 1), apenas a media e matriz de covariancia de (xt|Dt) sao necessarias e estas
satisfazem as identidades
mt = E [E [xt|ηt, Dt]] (2.23)
e
Ct = V [E [xt|ηt, Dt]] + E [V [xt|ηt, Dt]] (2.24)
Alem disso, assim como no caso normal, (xt|ηt, Dt) e condicionalmente independente de
It = {Yt,Ft}, e como Dt = {It, Dt−1}, tem-se que os momentos condicionais nas esperancas
internas em (2.23) e (2.24) sao (xt|ηt, Dt−1). Em geral, estes momentos serao desconhecidos,
funcoes nao lineares de ηt, sendo a unica informacao disponıvel aquela que diz respeito aos
momentos conjuntos de (g (ηt) ,x′t|Dt−1),
g (ηt)
xt|Dt−1
∼ ft
at
,
qt S′t
St Rt
, (2.25)
onde a matriz de covariancia completa e singular. Na base desta informacao por si, uma
abordagem alternativa e necessaria para que a informacao em It possa ser filtrada de volta
para xt.
13
Atualizacao do campo latente
O metodo linear Bayes pode ser aplicado no modelo anterior para fornecer feedback das
informacoes em It para xt. A densidade de p(xt|ηt, Dt−1) e a distribuicao desconhecida
preditiva de xt, dado ηt; a media e o preditor otimo, no sentido de minimizar o traco da
funcao de risco quadratica [At (d)] em relacao a d, em que:
At (d) = E[(xt − d) (xt − d)′ |ηt, Dt−1
]. (2.26)
A matriz de covariancia da distribuicao e o valor de At(d) na media.
Agora, uma vez que a media e desconhecida, um preditor alternativo e procurado atraves
da abordagem linear Bayes. Tendo em vista a relacao (2.9) e a construcao de ηt, e natural
que se adote uma funcao linear de g(ηt) como preditor de xt. Especificamente, suponha que
d deva ser escolhido de tal modo que d = d0 + d1g(ηt) para algum d0 e d1 e que, em vez
de (2.26), d minimize o risco global quadratico (ou soma de variancias) dado por
rt (d) = tracoE [At (d) |Dt−1 ] ,
onde a esperanca e com relacao a p(ηt|Dt−1).
Neste modelo, os momentos conjuntos (2.25) sao suficientes para determinar o preditor
requerido. Diretamente minimizando rt(d) em relacao a d0 e d1, obtem-se um unico mınimo
em d = at, onde
at = at + St (g (ηt)− ft) /qt (2.27)
o valor de E [At (d) |Dt−1 ] no mınimo e dado por
Rt = Rt − StS′
t/qt (2.28)
Os valores at e Rt, fornecem um preditor linear otimo de xt|ηt, Dt−1 e a medida do risco
associado e um problema nao linear. A alimentacao da informacao de It, pode agora ser
completada substituindo-se a media condicional e matriz de covariancia em (2.23) e (2.24)
por at e Rt para se obter o preditor esperado e risco, dados por:
14
mt =at + St (gt − ft) /qt
Ct =Rt − StS′t (1− pt/qt) /qt ,
onde gt = E [g (ηt) |Dt] e pt = V [g (ηt) |Dt] sao calculadas pela posteriori conjugada de
(ηt, Dt).
Modelo Poisson Dinamico
No capıtulo 3, o metodo proposto por West et al. (1985) sera aplicado a contagens
epidemiologicas. Suponha-se, em particular, que tais observacoes sigam uma distribuicao
Poisson com media λt. O modelo e definido pelas seguintes quatro componentes: equacao
de observacao, distribuicao a priori, funcao de ligacao e a evolucao de estados.
Considere os seguintes componentes essenciais do analise para o modelo dinamico Poisson:
• Modelo observacional
yt ∼ Poisson(λt)
p(yt|λt) = exp [ytlog(λt)− λt]1
yt!, (2.29)
em que φ = 1, ηt = logλt com ηt parametro natural e a(ηt) = λt = eηt , sendo a
media e variancia E[yt|ηt, φ] = µt = a′(ηt) = eηt = λt e V [yt|ηt, φ] = a′′(ηt)/φ = eηt
respectivamente.
• Priori para (ηt|Dt−1) ∼ CP [rt, st].
No caso Poisson, especificamos uma priori log-Gama para ηt, ou seja, uma priori Gama
para λt:
λt|Dt−1 ∼ CP [rt, st] = gama(rt, st)
O par (αt, βt) e deduzido usando propriedades da famılia exponencial.
• A funcao de ligacao e a equacao do sitema:
ηt = log(λt) = F′
txt
xt = Gtxt−1 + wt ∼ [0,Wt],
15
A fim de iniciar o procedimento de estimacao sequencial precisamos a informacao inicial
de x0.
• Informacao inicial:
(x0|D0) ∼ (m0,C0),
A natureza sequencial de modelos dinamicos e conseguida atraves da ciclagem de tres
passos: evolucao, equalizacao dos parametros e atualizacao, a partir de t = 1, · · · , T . As
distribuicoes sao apenas parcialmente especificadas em termos de seus momentos. Para um
determinado tempo t, os passos (1)-(3) sao descritos a seguir.
1. Evolucao:
• Prioris para o parametro de estado e do preditor linear:
xt|Dt−1 ∼ [at, Rt]
ηt|Dt−1 ∼ [ft, qt]
• Priori para λt: ja que o parametro λt > 0 e real positivo, uma escolha natural
para a priori e a famılia gama: (λt|Dt−1) ∼ CP [rt, st] = gama(rt, st), em que
rt, st > 0. Os seus dois primeiros momentos sao conhecidos e serao utilizados na
solucao de um sistema nao linear simples, a fim de obter os valores dos parametros
(rt, st) consistentes com (ft, qt), os momentos de (ηt|Dt−1). Os detalhes sobre a
solucao do sistema nao-linear sao descritos no proximo passo.
2. Equalizacao dos parametros:
Considerando-se que o preditor linear esta relacionado com a media da distribuicao
observacional por meio de uma funcao de ligacao, alguma aproximacao e necessaria
para determinar os hiperparametros rt e st da distribuicao a priori de λt.
Da priori de (ηt|Dt−1) e a transformacao ηt = log(λt) obtemos a priori (λt|Dt−1) como
uma distribuicao gama, isto e (λt|Dt−1) = gama(rt, st) com densidade
16
p(λt|Dt−1) =srtt
Γ(rt)λrt−1exp (−stλt) ,
O par (rt, st) e deduzido usando propriedades da famılia exponencial, isto e,
p(λt|Dt−1) = exp
[(rt − 1)log(λt)− stλt + log
(srtt
Γ(rt)
)]e ηt = log(λt), entao temos:
p(ηt|Dt−1) = exp
rt log(λt)︸ ︷︷ ︸ηt
−st λt︸︷︷︸exp(ηt)
+log
(srtt
Γ(rt)
)com
T = (T1(λt), T2(λt)) = (logλt,−λt)
b(rt, st) = −rtlog(st) + logΓ(rt).
Entao,
E[T1] =∂b
∂rt= −log(st) + ψ(rt)
E[T2] =∂b
∂st= −rt
st
V [T1] =∂2b
∂r2t
= ψ′(rt)
V [T2] =∂2b
∂s2t
=rts2t
Cov[T2] =∂2b
∂st∂rt= − 1
st,
com ψ(z) a funcao digamma, definida por ψ(z) = dlog(Γ(z))dz
e ψ′(z) = dψ(z)
dz, a funcao
trigamma (Abramowitz e Stegun (1964)). Da teoria associada a famılia exponencial,
temos
17
ft = E[ηt|Dt−1] = E[log(λt)|Dt−1]
= E[T1] = −log(st) + ψ(rt)
qt = V [ηt|Dt−1] = V [log(λt)|Dt−1]
= V [T1] = ψ′(rt)
Com base na avaliacao da media e variancia de log(λt) e uma aproximacao numerica
da funcao digamma dada por ψ(z) ≈ log(z) e ψ′(z) ≈ z−1 segundo Abramowitz e
Stegun (1964), temos
ft ≈ −log(st) + log(rt) = log
(rtst
)(2.30)
qt ≈1
rt(2.31)
Resolvendo as equacoes (2.33) e (2.31) , temos
rt =1
qt, st = exp
(−ftqt
)(2.32)
com isso temos a priori conjugada para λt e completamente especificada e (ηt|Dt−1) ∼[ft = log
(rtst
), qt = 1
rt
].
3. Distribuicao preditiva um passo a frente:
A distribuicao incondicional da distribuicao preditiva um passo a frente e obtida atraves
da integracao de λt:
18
p(yt|Dt−1) =
∫p(yt, λt|Dt−1)dλt =
∫p(yt|λt, Dt−1)p(λt|Dt−1)dλt
∝∫
1
yt!exp(−λt)λytt
srtt λrt−1t
Γ(rt)e−stλtdλt
∝ srttΓ(rt)yt!
∫λ
(yt+rt)−1t e−(st+1)λtdλt
∝ srttΓ(rt)yt!
× Γ(yt + rt)t(1 + st)yt+rt
, entao
p(yt|Dt−1) =Γ(yt + rt)t
Γ(yt + 1)Γ(rt)
(st
st + 1
)rt ( 1
st + 1
)yt,
que e uma distribuicao binomial negativa, denotada por yt|Dt−1 ∼ Bin neg(rt,
1st+1
).
A media e a variancia da distribuicao preditiva podem ser calculadas usando esperancas
condicionais, isto e,
E(yt|Dt−1) = E (E(yt|λt)|Dt−1) =rtst
V (yt|Dt−1) = E (V (yt|λt)|Dt−1) + V (E(yt|λt)|Dt−1) =rt(st + 1)
s2t
.
4. Atualizacao:
• Posteriori para λt: A distribuicao posterior de λt e obtida usando o teorema de
Bayes. Seja
p(λt|Dt) =p(yt|λt, Dt−1)p(λt|Dt−1)
p(yt|Dt−1)
∝ p(yt|λt, Dt−1)p(λt|Dt−1)
∝ 1
yt!exp(−λt)λytt
srtt λrt−1t
Γ(rt)e(−stλt)
∝ srttΓ(yt + 1)Γ(rt)
λyt+rt−1t exp(−(st + 1)λt)
que e a distribuicao gama, denotada por λt|Dt ∼ gama (yt + rt, st + 1). Da
teoria associada a famılia exponencial e por analogia, o calculo de gt e pt que sao
a media e a variancia a posteriori do preditor linear ηt, respectivamente, temos
19
gt = E[ηt|Dt] = E[log(λt)|Dt] = −log(st + 1) + ψ(yt + rt)
pt = V [ηt|Dt] = V [log(λt)|Dt] = ψ′(yt + rt),
que podem ser calculados recursivamente, pois:
ψ(z) = ψ(z + 1)− z−1
ψ′(z) = ψ
′(z + 1) + z−2.
Utilizando a aproximacao numerica da funcao digamma dada por ψ(z) ≈ log(z)+
(2z)−1 e ψ′(z) ≈ 1
z− 1
2z2segundo Abramowitz e Stegun (1964), temos
gt = −log(st + 1) + ψ(yt + rt)
≈ −log(st + 1) + log(yt + rt) +1
2(yt + rt)
= log
(yt + rtst + 1
)+
1
2(yt + rt)
pt = ψ′(yt + rt) ≈
1
(yt + rt)− 1
2(yt + rt)2
=2(yt + rt)− 1
2(yt + rt)2
• Atualizacao dos estados: A distribuicao conjunta de xt e ηt e parcialmente
especificada e obtida a partir dos resultados anteriores. O metodo de estimacao
linear bayesiana West e Harrison (1997) podem ser utilizados para obter xt|Dt ∼
[mt,Ct].
com mt = at + St(gt − ft)/qt e Ct = Rt − StS′t(1− pt/qt)/qt.
• Transicao de estado: at = Gtmt−1 e Rt = BtGtCt−1G′tBt.
Outros metodos de aproximacao de inferencia bayesiana sao os metodos de simulacao
estocastica, em particular os metodos de Monte Carlo via Cadeias de Markov (MCMC) e o
metodo determinıstico Integrated Nested Laplace Approximation (INLA). Estos metodos sao
centrais nesta dissertacao e sao desenvolvidos de forma mas detalhada no capıtulo 3.
20
Capıtulo 3
Metodos de Aproximacao MCMC e
INLA
Neste capıtulo, descrevemos os metodos de aproximacao MCMC e INLA para realizar
inferencia bayesiana completa em uma classe de modelos de espaco de estados. Em termos
gerais os metodos de Monte Carlo via Cadeias de Markov (MCMC), baseados em simulacao
estocastica que estao relacionados ao processo de obtencao de amostras da distribuicao a
posteriori para sumarizar informacao e que sao descritos de forma detalhada por Gamerman
e Lopes (2006). Por outro lado o metodo de aproximacao determinıstico, Integrated Nested
Laplace Approximation (INLA), proposto por Rue et al. (2009), combinando aproximacoes
Laplace e integracao numerica tornando este metodo eficiente (ver Rue e Martino (2007),
Rue et al. (2009), para um tratamento mais extenso).
O metodo de aproximacao INLA, calcula diretamente aproximacoes muito precisas
para as marginais a posteriori de interesse, nao passando por atualizacao recursiva, como
metodos baseados em variacoes do filtro de Kalman ou por procedimento iterativo, como
metodos MCMC. Em comum com esses ultimos, o INLA fornece a posteriori dos estados e
hiperparametros com respeito a toda a amostra observada, mas seu principal benefıcio e o
tempo computacional reduzido, em comparacao a metodos MCMC. Os metodos MCMC, em
contrapartida, aplicam-se a classes mais abrangentes, nao sujeitas as imposicoes descritas a
seguir sobre a forma dos modelos trataveis via INLA.
Na secao 3.1 apresentamos o metodo MCMC e na secao 3.2 apresentamos a metodologia
21
INLA, e para a ilustracao dos metodos, exibimos uma aplicacao com dados artificias, em
que a variavel resposta nao-gaussiana segue uma distribuicao Poisson. No dado artificial,
apresentamos comparacao entre os metodos INLA, MCMC e Linear Bayes (LB) comentando
restricoes do LB em relacao aos demais, e ganho de eficiencia computacional ao se usar o
metodo INLA, ao inves de MCMC.
3.1 Metodo de Aproximacao MCMC
Em modelos dinamicos, como vimos no capıtulo anterior, no caso que Ft,Gt,Wt sao
totalmente conhecidas, sob normalidade dos erros e se uma forma conjugada e imposta a
Vt = V , ∀t, entao tem-se inferencia bayesiana completa, de forma analıtica. E natural
assumir-se, entretanto, que Wt seja desconhecida. Uma alternativa e a especificacao de Wt
por meio de fatores de desconto, como descrito na secao anterior, mas pode-se ter interesse
na estimacao de Wt ou de quantidades desconhecidas em Ft e Gt. Ainda, a suposicao de
normalidade dos erros pode nao ser valida. Nesse caso, nao ha forma analıtica fechada para
distribuicoes a priori, preditiva e posteriori.
Em particular, no caso MDLG, devido a verossimilhanca construıda com base na famılia
exponencial, associada a prioris nao conjugadas, nao se obtem forma fechada para a densidade
a posteriori de diversos parametros, ao contrario do que ocorre nos modelos dinamicos normais
Alves (2006).
Em inferencia bayesiana, os problemas nao solucionados analiticamente podem ser
resolvidos usando metodos de simulacao que estao relacionados ao processo de obtencao
de amostras de distribuicoes a posteriori. Os metodos de Monte Carlo via Cadeias de
Markov (MCMC) sao metodos de simulacao estocastica, amplamente utilizados na inferencia
bayesiana nas duas ultimas decadas, quando se tem interesse em simular amostras de uma
determinada distribuicao a posteriori, a qual nao possui forma analıtica conhecida.
A ideia basica do metodo MCMC consiste em construir uma cadeia de Markov que, por
meio de escolhas adequadas de nucleos de transicao, tenha como distribuicao estacionaria a
distribuicao de interesse: no contexto bayesiano, a distribuicao a posteriori. Tais metodos
requerem ainda que a cadeia de Markov seja homogenea (as probabilidades de transicao de
22
um estado para outro sao invariantes), irredutıvel (cada estado pode ser atingido a partir de
qualquer outro em um numero finito de iteracoes) e aperiodica(nao haja estados absorventes),
cuja distribuicao estacionaria seja igual a distribuicao de interesse.
Uma vez que a convergencia da cadeia tenha sido atingida, as amostras estarao sendo
geradas da distribuicao estacionaria.
A principal vantagem desta abordagem e a possibilidade de se fazer a analise bayesiana
completa, o que significa tratamento formal da incerteza devida ao fato de que os
hiperparametro θ sao desconhecidos, sendo possıvel integrar θ a fim de apresentar inferencia
sobre (x1, · · · ,xT ). Alem disso, a estimacao pontual e a estimacao por intervalo de θ podem
ser feitas com base na distribuicao a posteriori.
Quando a distribuicao condicional completa de um parametro de interesse esta disponıvel
para amostragem, usualmente adota-se o amostrador de Gibbs, caso particular de algoritmo
MCMC, descrito a seguir. Em MDLGs, entretanto, nao se consegue amostrar a condicional
completa de xt. Existem algumas propostas de implementacao do amostrador de Gibbs para
casos particulares e o algoritmo Metropolis Hastings e indicado para as aplicacoes em geral,
tais algoritmos serao apresentados na subsecoes seguintes. Detalhes sobre metodos MCMC
podem ser vistos em Gamerman e Lopes (2006).
A difusao da aplicacao destes metodos foi iniciada com o trabalho de Gelfand e Smith
(1990), no qual foi feita uma comparacao entre o amostrador de Gibbs, proposto inicialmente
por Geman e Geman (1984), com outros esquemas de simulacao estocastica. Ate entao,
os trabalhos desenvolvidos eram baseados principalmente em aproximacoes numericas e
analıticas. O avanco computacional na decada de 1990 facilitou a popularizacao de aplicacoes
dos metodos bayesianos.
3.1.1 Amostrador de Gibbs
O amostrador de Gibbs foi proposto por Geman e Geman (1984), sendo popularizado por
Gelfand e Smith (1990). O amostrador de Gibbs e um esquema iterativo de amostragem
de uma cadeia de Markov, utilizando tal esquema para amostrar uma distribuicao a
posteriori p(x) do vetor parametrico x = (x1,x2, · · · ,xd)′, desde que as distribuicoes
23
condicionais completas destes parametros sejam conhecidas para amostragem. A distribuicao
condicional completa e a distribuicao da j−esima componente de x condicionada a informacao
de todos os outros parametros, podendo ser denotada por p(xj | x−j) onde x−j =
(x1, · · · ,xj−1,xj+1, · · · ,xd). A distribuicao condicional completa e obtida a partir da
distribuicao conjunta.
Em muitos casos, a geracao de uma amostra diretamente de p(x) pode ser custosa ou
impossıvel, mas se as distribuicoes condionais completas forem completamente conhecidas,
entao o algoritmo do amostrador de Gibbs pode descrito atraves dos seguintes passos:
(i) Especifique-se um valor arbitrario inicial para x(0) para x e inicialize-se o contador
j = 1;
(ii) Obtem-se um novo valor x(j) a partir de x(j−1) atraves de geracoes sucessivas dos
valores
x(j)1 ∼ p(x1 | x(j−1)
2 , . . . ,x(j−1)d )
x(j)2 ∼ p(x2 | x(j)
1 ,x(j−1)3 , . . . ,x
(j−1)d )
x(j)3 ∼ p(x3 | x(j)
1 ,x(j)2 ,x
(j−1)4 , . . . ,x
(j−1)d )
...
x(j)d ∼ p(xd | x(j)
1 ,x(j)2 , . . . ,x
(j)d−1).
(iii) Atualiza-se o contador de j para j+ 1 e retorna-se ao passo (ii) ate que a convergencia
seja obtida.
Carlin et al. (1992) introduziram o uso do amostrador Gibbs para realizar inferencia para
modelos de espaco de estados nao-normais e nao-lineares, com perturbacoes tanto na equacao
observacao e na equacao do sistema. Embora Carlin et al. (1992) nao trabalhem com MDLG,
a ideia motivou algumas solucoes baseadas em MCMC para MDLG.
Fahrmeir et al. (1992) propos o amostrador de Gibbs para analise de MDLG utilizando
amostragem por rejeicao para gerar a partir das condicionais completas de xt. No contexto
de MDLGD a equacao do sistema tem perturbacoes normais, e assim existem mt e Ct tal
24
que p(xt | x−t, Dt,W ) ∝ p(yt | xt)N(mt,Ct) = f(xt),∀t. Fahrmeir et al. (1992) entao
propuseram g(xt) = N(mt, Ct) o que parece natural, mas pode levar a baixas taxas de
aceitacao.
A desvantagem e que este abordagem nao e aplicavel a famılia exponencial em geral,
pois nao e possıvel amostrar as distribuicoes condicionais completas de todos os parametros
de interesse, como descritos no paso 2 acima. Nestes casos, passos de Metropolis Hastings
podem ser inseridos no algoritmo de Gibbs
3.1.2 Algoritmo de Metropolis Hastings
O metodo de Metropolis-Hastings foi apresentado por Metropolis et al. (1953) e
posteriormente estendido por Hastings (1970), resultando o algoritmo de Metropolis-Hastings.
O metodo e usado geralmente para a geracao de amostras da distribuicao de parametros de
interesse cujas condicionais completas nao tenham forma analıtica fechada. Nesse caso,
sao gerados valores para cada parametro a partir de uma distribuicao auxiliar, chamada de
distribuicao proposta q(·), tais valores sao aceitos ou nao com uma certa probabilidade.
Para superar a dificuldade de amostragem das distribuicoes condicionais completas de
todos os parametros de interesse, alguns autores vem apresentando abordagens com base no
algoritmo de Metropolis-Hastings. A ideia consiste em usar o amostrador de Gibbs com passos
de Metropolis-Hastings para conseguir amostrar as densidades dos parametros de estado. No
caso da amostragem dos parametros de estado, um de cada vez, este passo de Metropolis-
Hastings e composto de dois passos. Sejam p(·) a funcao densidade de probabilidade de
interesse e q(·) a funcao de probabilidade ou funcao densidade de probabilidade proposta.
Entao:
(i) Gera-se um valor ξ de uma distribuicao proposta q(x(j−1), ξ).
(ii) Aceita-se o valor gerado em (ii) com probabilidade
α (x, ξ) = min
{1,
p(ξ)/q(x(j−1), ξ)
p(x(j−1))/q(ξ,x(j−1))
}.
Se o valor for aceito, x(j) = ξ. Caso contrario, a cadeia nao se move e x(j) = x(j−1).
25
A diferenca entre os metodos apresentados por diferentes autores, e a amostragem
individual ou em blocos dos estadoe e a escolha da distribuicao proposta q(x(j−1), ξ). Embora
a distribuicao proposta possa ser escolhida arbitrariamente, na pratica e crucial, e deve-se
tomar alguns cuidados para garantir rapida convergencia para a distribuicao a posteriori e
eficiencia do algoritmo MCMC.
Fruhwirth Schnatter (1994) e Carter e Kohn (1994) propoem amostrar todos os parametros
de estado em forma simultanea, usando o amostrador de Gibbs, acelerando a convergencia
das cadeias em comparacao ao caso de movimentos individuais.
Shephard e Pitt (1997) e Knorr-Held (1999) sugeriram o uso de blocos aleatorios (cujos
tamanhos podem variar ), contendo alguns dos componentes de xt. Suas propostas foram
baseadas em aproximacoes normais para a verossimilhanca e movimentos de passeio aleatorios,
respectivamente, utilizando a metodologia via algoritmo de Metropolis Hastings como uma
ferramenta simples, porem flexıvel e eficiente para inferencia bayesiana via MCMC em modelos
dinamicos.
Gamerman (1998) apresenta um abordagem metodologica para a realizacao de inferencia
bayesiana sobre modelos dinamicos para observacoes da famılia exponencial, baseando-se em
simulacao, utilizando a tecnica MCMC, considerando passos de Metropolis-Hastings dentro
do amostrador de Gibbs. Nesse trabalho o autor sugeriu o uso de um modelo dinamico
normal ajustado com o objetivo de construir uma boa densidade proposta, q(.), baseada no
modelo do trabalho de Singh e Roberts (1992) definido por (2.11), (2.12) e (2.13). Para
organizar esta ideia, implementa tres esquemas de amostragem diferentes: a amostragem
individual ou em blocos dos estados, ou amostragem individual dos termos do erro da equacao
de evolucao; concluindo que os movimentos individuais sao preferıveis aos movimentos em
blocos e que amostrar os erros e mais eficiente, porque as cadeias resultantes sao menos
autocorrelacionadas, o que ajuda a acelerar a convergencia do algoritmo. Para maiores
detalhes sob o algoritmo, veja Tierney (1994), Gamerman (1997) e Ferreira e Gamerman
(2000) Gamerman e Lopes (2006).
Migon et al. (2013) Propoem um esquema de amostragem em blocos para os parametros
de estados de modelos dinamicos nao-gaussianos e nao-linear de series temporais univariadas.
Este esquema de amostragem combina a abordagem de Conjugate Updating para modelos
26
lineares dinamicos generalizados, com amostragem retrospectiva (Backward Sampling) dos
parametros de estado usados em modelos lineares dinamicos normais. Denominado algoritmo
CUBS, abreviacao do ingles Conjugate updating backward sampling. As vantagens do
algoritmo sao a facilidade de implementacao e a reducao significativa do tempo computacional
necessario para se atingir a convergencia das cadeias.
Para a avaliacao da convergencia, existem metodos formais para diagnostico de
convergencia das cadeias, como os propostos por Geweke (1992) e a estatıstica de Gelman e
Rubin (1992). Entretanto, nenhum destes metodos e conclusivo, fornecendo apenas indıcios
de convergencia.
Um modo informal simples de verificacao da convergencia, adotado neste trabalho, e a
inspecao visual das cadeias de Markov, que consiste em analisar a trajetoria de pelo menos
duas cadeias independentes (definidas por diferentes valores iniciais) dos proprios parametros
que oscilam na mesma regiao, mesclando-se, a medida em que o numero de iteracoes cresce
e verificar se todas convergem para a mesma regiao de estabilidade.
Vale ressaltar que tecnicas informais de verificacao da convergencia devem ser utilizado
com cautela e sempre acompanhadas de alguma fundamentacao teorica.
Supondo que sejam realizadas N iteracoes do algoritmo MCMC e que a partir da iteracao J
a cadeia tenha atingido a convergencia para a distribuicao de equilıbrio, os valores xJ , . . . ,xN
constituem uma amostra aproximada da distribuicao a posteriori de x.
3.2 O Metodo de Aproximacao INLA
Modelos de espaco de estados pressupoem uma estrutura markoviana de evolucao entre
seus estados e, ainda, que condicionalmente a tais estados, as observacoes sejam serialmente
independentes Petris et al. (2009). Modelos dinamicos lineares gaussianos, como descritos
em West e Harrison (1997) sao, portanto, modelos de espaco de estados em que se pressupoe
normalidade, tanto da variavel resposta, quanto dos erros de evolucao dos estados. O metodo
INLA e adequado a modelos mais gerais que os MDLs, ja que nao impoe a restricao de
normalidade sobre a variavel reposta. Ao inves disso, assume-se resposta pertencente a famılia
exponencial, explicada por preditores aditivos (e, portanto, nao necessariamente lineares)
27
e, ainda, que seu vetor de estados x = (x1, . . . , xn)′ possa ser escrito como um campo
latente gaussiano markoviano com matriz de precisao Q(θ) esparsa, definida por um vetor
de hiperparametros θ de dimensao m baixa (m ≤ 6). θ pode incluir, ainda, parametros
definidores de p(yi|x). Segundo Rue et al. (2009), usualmente a eficiencia do metodo depende
da baixa dimensao do vetor de hiperparametros e da esparsidade da matriz de precisao, mas
ha excecoes. A esparsidade de Q viabiliza a adocao de metodos numericos muito mais rapidos
que seus equivalentes para matrizes densas.
Sejam yi, i = 1, . . . , T os dados observados, supostos condicionalmente
independentes,dado o campo latente x, ou seja,
p(y|x,θ) =∏i
p(yi|xi,θ), (3.1)
onde cada observacao esta conectada a somente um elemento do campo latente. Apos
atribuir uma priori para o vetor de hiperparametros θ, temos a distribuicao a posteriori de
interesse, dada por:
p(x,θ|y) ∝ p(θ)p(x|θ)∏i
p(yi|xi,θ) (3.2)
∝ p(θ)|Q(θ)|1/2exp
{−1
2x′Q(θ)x +
∑i
log(p(yi|xi,θ))
}. (3.3)
O principal objetivo da abordagem INLA e aproximar as distribuicoes a posteriori marginais
para as quantidades latentes p(xi|y), ∀i = 1, . . . , n e para os hiperparametros do modelo
gaussiano latente p(θj|y), ∀j = 1, . . . ,m, dadas por:
p(xi|y) =
∫p(θ|y)p(xi|θ,y)dθ (3.4)
p(θj|y) =
∫p(θ|y)dθ−j, (3.5)
Como as integracoes em (3.4) e (3.5) nao podem ser resolvidas analiticamente, sao
necessarias aproximacoes para p(xi|θ,y) e p(θ|y) denotadas respectivamente por p(xi|θ,y)
e p(θ|y). Ainda, as integrais sao possıveis devido a baixa dimensao suposta para θ.
28
p(xi|y) =
∫p(θ|y)p(xi|θ,y)dθ (3.6)
p(θj|y) =
∫p(θ|y)dθ−j. (3.7)
Estas aproximacoes baseiam-se em uma combinacao eficiente de aproximacoes de Laplace
(proposta por Tierney e Kadane (1986)) para as condicionais completas p(θ|y) e p(xi|θ,y),
i = 1, · · · , n, e rotinas de integracao numerica e obtidas atraves de
p(xi|y) ≈∑j
p(xi|θj,y)p(θj|y)∆j (3.8)
p(θj|y) ≈∑i
p(θi|y)∆ji (3.9)
Para resolver (3.8) e (3.9) e necessario obter as aproximacoes de p(θ|y) e p(xi|θ,y) e
avaliar p(θ|y) em uma grade, descrito a seguir (secoes 3.2.1 e 3.2.1).
3.2.1 Parametrizacao adequada do vetor parametrico e exploracao
da grade
Quando se aplica metodos numericos para a obtencao das distribuicoes de interesse, em
muitos casos e necessario reparametrizar o vetor parametrico, ja que boas reparametrizacoes
simplificam integracoes numericas e facilitam a construcao de grades otimizadas.
Parametrizacao adequada do vetor parametrico
Segundo Rue et al. (2009), e conveniente fazer uma reparametrizacao do vetor parametrico
θ = (θ1, · · · ,θm) ∈ <m, sendo m a dimensao de θ, de modo a obter um novo vetor,
por exemplo, z, cuja distribuicao seja aproximadamente uma distribuicao normal padrao
multivariada. Essa transformacao ira facilitar a construcao de grades otimizadas, em que
as regioes de maior massa de probabilidade tenham mais pontos em relacao as demais areas.
Na literatura sao conhecidas varias propostas para achar a melhor reparametrizacao. Em
particular, a abordagem de Smith et al. (1987), fornece bons resulatdos. Para maiores
29
referencias sugerimos Martins (2010). Segundo Smith et al. (1987), procede-se em duas
etapas:
1. Reparametriza-se os parametros individuais de tal forma que tenham suporte na reta
real. Por exemplo, parametros com suporte em intervalos da forma (0,∞) sao
transformados utilizando θ∗ = log(θ) e, quando o intervalo e do tipo (a, b), a
transformacao utilizada e θ∗ = [log(θ − a) − log(b − θ)]. Denote-se θ∗ como o
novo vetor parametrico, tal que θ∗ ∈ <m.
2. Deseja-se, agora, um conjunto centrado, padronizado e mais ortogonal. Para isso,
calcular a moda da densidade a posteriori de θ∗, denotada por θ∗Mo e encontrar a
matriz Hessiana, H avaliada na moda θ∗Mo. Considerar a decomposicao espectral, Σ,
dada por Σ = VΛV′ e que Σ = −H−1. O novo conjunto de parametros z e tal que
θ(z) = θ∗Mo + VΛ1/2z ou z = VTΛ−1/2(θ − θ∗Mo)
Exploracao da Grade
A ideia e explorar log[p(θ|y)] suficientemente bem, a fim de selecionar bons pontos, o
que e crucial para a precisao da integracao numerica e assim localizar as regioes com maior
massa de probabilidade. Portanto, nao e necessario representar p(θ|y) parametricamente.
Rue et al. (2009) sugere fazer tal exploracao, utilizando a parametrizacao z, detalhados da
seguinte forma:
1. Iniciar pela moda z = 0 (θ = θ∗Mo), e ir na direcao positiva de z1 com espacamento de
tamanho δz, por exemplo, δz = 1; ate quando a seguinte condicao for satisfeita
log{p(0|y)} − log{p(θ(z)|y)} < δπ (3.10)
sendo por exemplo δπ = 2.5.
2. Repetir o procedimento na direcao negativa de z1. Todas as demais coordenadas de z
seguem o mesmo procedimento analogo. Os pontos assim obtidos sao denominados de
pontos pretos.
30
3. Completar a construcao da grade com todas as combinacoes dos pontos pretos,
verificando para cada combinacao se a condicao (3.10) e satisfeita. Esses novos pontos
sao chamados de pontos cinzas.
O conjunto de pontos pretos e cinzas e utilizado para obter as distribuicoes marginais de
θ|y via integracao numerica.
3.2.2 Aproximacao para p(θ|y)
Seja a identidade
p(θ|y) =p(x,θ|y)
p(x|θ,y), ∀x. (3.11)
Rue e Martino (2007) propoem aproximar p(θ|y), baseadas em aproximacoes de Laplace
para a distribuicao a posteriori marginal de θ, dada por
p(θ|y) ∝ p(x,θ,y)
pG(x|θ,y)
∣∣∣∣x=x∗(θ)
, (3.12)
em que pG(x|θ,y) do denominador de (3.12) e uma aproximacao Gaussiana para a
distribuicao condicional completa de x, e x∗(θ) e a moda da distribuicao condicional completa
de x para um dado θ,
p(x|θ,y) ∝ exp
{−1
2x′Q(θ)x +
∑i
log(p(yi|xi,θ))
}. (3.13)
A moda, x∗(θ) e obtida atraves de algum metodo de otimizacao, por exemplo,
utilizando o metodo de Newton-Raphson descritos por Martins (2010) e Resende (2011).
A proporcionalidade de (3.12) e devido ao fato de nao conhecer a constante normalizadora
de p(x,θ|y).
Rue e Martino (2007) utilizaram a aproximacao (3.12) em diversos modelos latentes
Gaussianos e concluıram que tal aproximacao e muito precisa, de tal forma que nem longas
cadeias de MCMC foram capazes de detectar erros na aproximacao.
A distribuicao a posteriori p(θ|y) sera usada posteriormente para integrar a incerteza com
relacao a θ ao aproximar a marginal a posterior de xi.
31
3.2.3 Aproximacao para p(xi|θ,y)
Rue et al. (2009) sugerem tres abordagens de aproximacao para p(xi|θ,y): a aproximacao
Gaussiana, aproximacao de Laplace, aproximacao de Laplace simplificada, que variam em
termos de custo computacional e/ou precisao.
Duas de tais aproximacoes propostas para a obtencao das condicionais completas de um
modo mais simples sao aproximacoes Gaussiana e Laplace simplificada.
A aproximacao Gaussiana fornece resultados razoaveis em tempo computacional curto.
No entanto, de acordo com a Rue e Martino (2007), sua precisao pode ser afetada por erros
na localizacao e/ou erros devidos a falta de simetria. Estas deficiencias podem ser corrigidas
com um custo adicional razoavel, usando uma versao simplificada da aproximacao de Laplace,
definida como a expansao de Taylor da serie de pLA(xi|θ,y) em torno de xi = µi(θ), a media
de p(xi|y,θ) (para detalhes ver Rue et al. (2009), Ruiz-Cardenas et al. (2010), Martins
(2010), Resende (2011)).
Metodo de Laplace
Uma forma natural para melhorar a aproximacao Gaussiana e pela aproximacao de Laplace
para p(xi|θ,y), dada por:
pLA(xi|θ,y) ∝ p(x,θ,y)
pGG(x−i|xi,θ,y)
∣∣∣∣x−i=x∗−i(xi,θ)
, (3.14)
em que pGG(x−i|xi,θ,y) e a aproximacao Gaussiana para (x−i|xi,θ,y) e x∗−i(xi,θ) e a
moda de p(x−i|xi,θ,y).
E necessario avaliar a expressao (3.14) em diversos valores de xi para obter a densidade,
como por exemplo, {x(1)i , . . . , x
(k)i }. Estes pontos sao selecionados, considera-se x
(j)i =
µi(θ)+σi(θ)x(j)ab , em que µi(θ) e σi(θ) sao a media e desvio-padrao da aproximacao Gaussiana
pG(xi|θ,y) dada em (3.12) e x(j)ab e um ponto da abcissa dado pela quadratura de Gauss-
Hermite (Resende (2011)). Assim, a aproximacao de Laplace pLA(xi|θ,y) e representada
por:
pLA(xi|θ,y) ∝ N{µi(θ), σ2i (θ)}exp{h(xi)}, (3.15)
32
sendo a funcao h(x) uma funcao spline cubica ajustada a diferenca log {pLA(xi|θ,y)} −
{pG(xi|θ,y)}, avaliada nos pontos selecionados. Em seguida, a densidade (3.15) e
normalizada utilizando integracao por quadratura.
A aproximacao p(xi|θ,y) usando a expressao (3.14) pode levar a um alto custo
computacional, uma vez que pGG(x−i|xi,θ,y) deve ser recalculado para cada valor de xi
e θ uma vez que sua matriz de precisao depende de xi e θ.
Rue et al. (2009) propoem duas alteracoes a expressao (3.14) que a tornam viavel
computacionalmente. A primeira modificacao consiste em evitar o passo de otimizacao e
calcular pGG(x−i|xi,θ,y) atraves de uma aproximacao para a moda, dada por
x∗−i(xi,θ) ≈ EpG(x−i|xi). (3.16)
A segunda modificacao materializa a intuicao de que, devido a propriedades de Markov
do campo latente de x, apenas aqueles elementos do campo latente xj que sao de uma certa
forma ”proximos”de xi devem ter um efeito (alguma influencia) sobre a distribuicao marginal
de xi. Denote-se por Ri(θ) a ”regiao de interesse”com relacao a distribuicao marginal de xi.
A esperanca condicionada na approximacao acima implica que
EpG(xj|xi)− µj(θ)
σj(θ)= aij(θ)
xi − µi(θ)
σi(θ)(3.17)
para algum aij(θ), j 6= i. Portanto, considere-se uma simples regra para a construcao de
Ri(θ),
Ri(θ) = {j : |aij(θ)| > 0.001}. (3.18)
O mais importante e que agora so sera necessario fatorizar matrizes esparsas de dimensao
|Ri(θ)| × |Ri(θ)|.
3.2.4 Algoritmo INLA
O metodo de aproximacao INLA, pode ser resumido em tres passos:
33
1. Explorar p(θ|y).
• Localizar a moda θ∗ de p(θ|y), otimizando log{p(θ|y)} com respeito a θ.
• Usar a Hessiana para construir a nova variavel z, isto e, fazendo a transformacao
adequada para θ, descrito na secao (3.2.1).
• Construir a grade, explorando log{p(θ|y)} utilizando a nova parametrizacao z.
• Aproximar p(θj|y) atraves de integracao numerica. A ideia principal e usar os
pontos calculados no item anterior e construir um interpolante para log{p(θ|y)}
e fazer a integracao numerica atraves deste interpolante.
2. Obter p(xi|θ,y).
Para cada θj da grade, encontrar a aproximacao de Laplace para cada elemento do
campo latente x.
3. Obter p(xi|y), conbinando os dois passos anteriores, por meio de integracao numerica,
isto e∑
j p(xi|θj,y)p(θj|y)∆j.
3.3 Modelo Poisson Dinamico com dados artificiais,
exemplo
Na presente secao, apresentamos o ajuste do modelo Poisson dinamico com dados
artificiais. Para a comparacao do modelo utilizamos os metodos de aproximacao INLA, MCMC
e Linear Bayes (LB), descrito na secao anterior (2.4.2).
Seja Y1, . . . , YT variaveis aleatorias condicionalmente independentes seguindo uma
distribuicao Poisson com media λt > 0, com funcao de probabilidae dada por:
p(yt|λt) = eλtλtyt 1
yt!
34
Assume-se que a media λt varia no tempo atraves de uma estrutura dinamica,
Yt | λt ∼ Poisson(λt)
ηt = log(λt) (3.19a)
ηt = F′txt = αt + βtzt (3.19b)
αt = αt−1 + vt, vt ∼ N [0,V], (3.19c)
βt = βt−1 + ωt, ωt ∼ N [0,W] (3.19d)
em que F′t = (1,wt) e um vetor de dimensao p = 2 que descreve os valores das covariaveis
observadas para cada indivıduo t e xt = (αt, βt) o vetor de parametros associados ao preditor
linear ηt.
Na estimacao bayesiana, devemos considerar as prioris para os parametros (αt, βt) no
instante inicial. Assume-se que
α1 ∼ N [0, R1] (3.20a)
β1 ∼ N [0, R2]. (3.20b)
sendo R1 e R2 valores conhecidos refletindo a incerteza a priori sobre o processo no
instante inicial.
Assim, as quantidades a serem estimadas sao Ψ = (V,W, α1, α1, . . . , αT ; β1, β2, . . . βT ),
em que θ = (V,W ) e o parametro fixo e (α1, α1, . . . , αT ; β1, β2, . . . βT ) parametros do campo
latente.
3.3.1 Prioris para os parametros fixos
Assume-se uma especificacao para as priori dos hiperparametros θ = (V,W ) que
sao independentes, e a distribuicao a priori para V ∼ GI(av, bv) e W ∼ GI(aw, bw)
respectivamente, sendo av, bv, aw e bw valores conhecidos.
Agora vamos reparametrizar os parametros fixos de modo que eles possuam suporte na
reta real.
35
V ∗ = log(1/V ) = log(τ1), sendo τ1 ∼ G(av, bv)
W ∗ = log(1/W ) = log(τ2), sendo τ2 ∼ G(aw, bw),
Sendo, τ1 = 1/V e τ2 = 1/W . Assim os calculos das prioris apos a nossa reparametrizacao
e,
p(V ∗) ∝ exp(V ∗)ave−bvexp(V∗) (3.21)
p(W ∗) ∝ exp(W ∗)awe−bwexp(W∗) (3.22)
3.3.2 Prioris para variaveis gaussianas latentes
A distribuicao conjunta dos parametros do campo latente xt = (α1, · · · , αT, β1, · · · , βT)
condicional a θ = (V,W ), e dada por:
p(x|θ) = p(α1)T∏t=2
p(αt|αt−1)p(β1)T∏t=2
p(βt|βt−1) (3.23a)
∝ exp
{−1
2R−1
1 α21
} T∏t=2
exp
{−1
2V −1 (αt − αt−1)2
}× (3.23b)
× exp{−1
2R−1
2 β21
} T∏t=2
exp
{−1
2W−1 (βt − βt−1)2
}
∝ exp
{−1
2
[R−1
1 α21 + V −1
T∑t=2
(αt − αt−1)2 +R−12 β2
1 +W−1
T∑t=2
(βt − βt−1)2
]}(3.23c)
mais
T∑t=2
(αt − αt−1)2 =∑T
t=2 α2t +
∑Tt=2 α
2t−1 − 2
∑Tt=2 α
2tα
2t−1 (3.24a)
= α21 + α2
T + 2∑T−1
t=2 α2t − 2
∑Tt=2 αtαt−1 (3.24b)
36
rescrevendo a distribuicao conjunta dos parametros do campo latente por
p(x|θ) ∝ exp
{−1
2
[(V −1 +R−1
1 )α21 + V −1α2
T + 2V −1
T∑t=2
α2t − 2V −1
T∑t=2
αtαt−1
(3.25a)
+(W−1 +R−12 )β2
1 +W−1β2T + 2V −1
T∑t=2
α2t − 2V −1
T∑t=2
βtβt−1
]}
p(x|θ) ∝ exp
{−1
2
[αTQ(V )α+ βTQ(W )β
]}(3.25b)
p(x|θ) ∝ exp
{−1
2xTQ(θ)x
}(3.25c)
Com isso temos (x|θ) ∼ N [0,Q(θ)], onde Q(θ) e a matriz de precisao, em que
Q(θ) =
Q(V) 0
0 Q(W)
2T×2T
(3.26)
em que,
Q(V) =
R−11 + V −1 −V −1 0 · · · 0 0
−V −1 2V −1 −V −1 · · · 0 0
0 −V −1 2V −1 · · · 0 0
0 0 −V −1 · · · 0 0...
......
. . ....
...
0 0 0 · · · 2V −1 −V −1
0 0 0 · · · −V −1 V −1
T×T
(3.27)
e
Q(W) =
R−12 +W−1 −W−1 0 · · · 0 0
−W−1 2W−1 −W−1 · · · 0 0
0 −W−1 2W−1 · · · 0 0
0 0 −W−1 · · · 0 0...
......
. . ....
...
0 0 0 · · · 2W−1 −W−1
0 0 0 · · · −W−1 W−1
T×T
(3.28)
37
Neste caso a verossimilhanca associada a η = (η1, · · · , ηT )′, sendo ηt = log λt temos:
p(yt|ηt, θ) ∝ e−λtλtyt = exp(−eηt)exp(ηt)yt (3.29)
log p(yt|ηt) ∝ −eηt + ytηt (3.30)
3.3.3 Aproximacao Gaussiana para a distribuicao Condicional
Completa xt
A aproximacao Gaussiana para a distribuicao condicional completa de xt, dado por
p(x|θ,y) ∝ p(y|x,θ)p(x|θ) =∏T
t=1 p(yt|xt,θ)p(xt|θ) (3.31)
∝ exp
−12xTQ(θ)x +
∑log p(yt|xt,θ)︸ ︷︷ ︸g(xt)=g(ηt)
(3.32)
Denote gt(xt) = log p(yt|xt,θ) como a verossimilhanca descrita em (3.30) e que xt =
ηt = (αt,βt) com t = 1, · · · , T entao g(ηt) = g(xt) e θ = (V,W ), sendo αt e βt
independentes.
gt(ηt) = log p(yt|ηt, θ) ∝ −eηt + ytηt (3.33)
ηt = αt + βtzt (3.34)
Em que o gradiente (∇g) e a matriz Hessiana (Hg) e dado por,
∇g =
(∂g(ηt)
∂αt,∂g(ηt)
∂βt
)= (−eηt + yt,−eηtzt + ytzt) (3.35)
Hg =
∂2g(ηt)
∂α2t
0
0 ∂2g(ηt)
∂β2t
=
−eηt 0
0 −eηt
(3.36)
Seja η(0) = (α(0),β(0)) um vetor de valores iniciais utilizados para a inicializacao do
algoritmo em busca da moda de p(θ | y). Expandindo gt(ηt) em torno de η(0)t , t = 1, · · · , T
segue:
38
gt(ηt) ≈ gt(η(0))︸ ︷︷ ︸
cte
+b(η(0))η − 1
2η′C(η(0))η
em que
b(η(0)) = ∇g(η(0))−Hg(η
(0))η(0)
C(η(0)) = −Hg(η(0))
fazendo as devidas operacoes algebricas,
bt(η(0)t ) =
(yt − eη
(0)t + eη
(0)t α
(0)t , ytzt − eη
(0)t zt + z2
teη(0)t β
(0)t
)(3.37)
Ct(η(0)t ) = diag
(eη
(0)t , eη
(0)t z2
t
)(3.38)
Logo
gt(ηt) ∝(eη
(0)t(α0t − 1
)+ yt, zte
η(0)t(ztβ
0t − 1
)+ ztyt
)ηt − ηTdiag
(eη
(0)t , eη
(0)t z2
t
)ηt
aplicando somatorios∑gt(ηt) =
∑t∈I log p(yt|ηt,θ) ≈ −1
2ηTC(η(0))η + b(η(0))ηt (3.39)
em que,
b(η(0)) =
eη
(0)1
(α
(0)1 − 1
)+ y1 z1e
η(0)1
(z1β
(0)1 − 1
)+ z1y1
eη(0)2
(α
(0)2 − 1
)+ y2 z2e
η(0)2
(z2β
(0)2 − 1
)+ z2y2
......
eη(0)T
(α
(0)T − 1
)+ yT zT e
η(0)T
(zTβ
(0)T − 1
)+ zTyT
(3.40)
e,
C(η(0)) = diag((eη
(0)1 , eη
(0)2 , · · · , eη
(0)T
),(eη
(0)1 z2
1, eη(0)2 z2
2, · · · , eη(0)T z2
T
))(3.41)
39
3.3.4 Resultados
Geramos uma serie temporal de dados artificiais, segundo um modelo Poisson dinamico
descrito acima, de tamanho T = 1444, sendo 1443 para estimacao do modelo e o ponto
restante, para a avalicao preditiva um passo a frente, isto devido ao MCMC levar horas para
rodar. A regressora z utilizada na simulacao de dados e uma das regressoras de valores reais
(monoxido de carbono) em que se tem interesse na aplicacao do Capıtulo 4. Portanto, a serie
temporal do modelo Poisson foi simulado para y condicionado a w, fixando-se V = 0.001
(1/V = 1000), W = 0.0005 (1/W = 2000). A figura 3.1 mostra a serie da regressora zt.
0 500 1000 1500
02
46
time
zt
Figura 3.1: Regressora zt utilizada na simulacao para o resposta Poisson.
A figura 3.2 mostra a serie temporal dos dados artificiais para o modelo poisson, gerados
com V = 0.001, W = 0.0005, e valores α1 e β1 sorteados, respectivamente, de α1 ∼
N(0, R1 = 10) e β1 ∼ N(0, R2 = 10).
Considerando a matriz de precisao a priori para o campo latente descrita em (3.26),
assumindo-se distribuicoes a priori para os parametros α1 ∼ N(0, R1 = 10) e β1 ∼ N(0, R2 =
10), e prioris nao informativas para os hiperparametros V e W com av = 2, bv = 0.0018,
aw = 2, bw = 0.0018, a estimacao foi feita usando os metodos INLA, MCMC e Linear Bayes.
40
0 500 1000 1500
05
1015
20
time
Y
Figura 3.2: Dados artificialmente gerados de um modelo Poisson dinamico, com covariaveis
de tamanho T = 1444.
A implementacao computacional do INLA e LB foram programadas utilizando o software R
(R Development Core Team (2013)) e para a abordagem do MCMC, o WinBUGS. Os codigos
sao apresentados no apendice A. No caso do INLA, foi tomado como base inicial o relatorio
de Albi et al. (2009).
No caso MCMC, para lograr convergencia das cadeias de Markov, foram simulados
duas cadeias de 27000 iteracoes cada uma, sendo as primeiras 5000 iteracoes consideradas
como perıodo de aquecimento, sendo as amostras restantes tomadas a espacamentos de 30
iteracoes, de forma a reduzir a autocorrelacao serial, com isso obtendo-se uma amostra final
de tamanho 1466.
A tabela (3.1) e a figura 3.3 mostram os valores estimados das distribuicoes a posteriori
41
dos hiperparametros V e de W , assim como os intervalos de credibilidade de 95% para cada
um deles, obtidas utilizando os metodos INLA e MCMC. Observe-se que ambos metodos
fornecem estimativas muito proximas e foram capazes de recuperar os valores teoricos dos
hiperparametros.
Tabela 3.1: Estimativas a posteriori com intervalos de credibilidade a 95% e valor verdadeiro
dos hiperparametros.
INLA MCMC
ICI media ICS ICI media ICS
V (0.001) 0.000549 0.001012 0.001728 0.000568 0.001059 0.001828
W (0.0005) 0.000247 0.000507 0.000961 0.000275 0.000553 0.000971
V
Den
sida
de
0.0005 0.0015 0.0025
020
060
010
0014
00
simulatedINLAMCMC
W
Den
sida
de
0.0000 0.0005 0.0010 0.0015
050
010
0015
0020
0025
00
simulatedINLAMCMC
Figura 3.3: Aproximacao da distribuicao a posteriori para V e W , pelos metodos INLA (linha
azul) e MCMC (histograma, linha preta), com intervalos de credibilidade ao 95% para INLA
(linha vertical tracejada azul) e MCMC (pontos pretos). A linha vertical vermelha refere-se
ao valor verdadeiro.
As Figuras 3.4 e 3.5 mostram as distribuicoes a posteriori marginais para os coeficientes
de regressao αt e o coeficiente βt obtidas utilizando os metodos INLA, MCMC e LB. Os
intervalos de credibilidade obtidos pelos metodos INLA e MCMC tem nıvel de credibilidade
de 95%. Ja os intervalos do LB tem nıvel de credibilidade indeterminado e foram obtidos
42
0 500 1000 1500
0.0
0.5
1.0
1.5
2.0
Tempo
alph
a_t
ArtificialINLAMCMCLB
Figura 3.4: Media a posteriori e intervalos de credibilidade de 95% para (α1, · · · , αT ) usando
os metodo INLA (linha azul), MCMC (linha vermelha) e LB (linha cinza). A linha preta
refere-se aos valores reais.
0 500 1000 1500
−0.
6−
0.2
0.2
0.4
0.6
Tempo
beta
_t
ArtificialINLAMCMCLB
Figura 3.5: Media a posteriori e intervalos de credibilidade de 95% para (β1, · · · , βT ) usando
os metodo INLA (linha azul), MCMC (linha vermelha) e LB (linha cinza). A linha preta
refere-se aos valores reais.
43
somando-se e subtraindo-se dois desvios padrao as medias a posteriori das componentes αt
e βt. Ainda, observe-se que de fato o modelo ajustado pelo metodo LB difere do modelo
ajustado por MCMC e INLA, uma vez que, neste dois ultimos metodos, pressupos-se V (αt) =
V , V (βt) = W (fixas no tempo) e no metodo LB em contrapartida, tem-se V (αt) = Vt,
V (βt) = Wt, variando no tempo e especificadas por meio de fator de desconto δ1 = 0.95
e δ2 = 0.97, respectivamente. Uma outra diferenca essencial entre os metodos deve-se ao
fato de que o LB processa informacao em tempo real, produzindo inferencia a cada tempo t
condicionalmente a informacao disponıvel ate esse instante, Dt. Ja INLA e MCMC produzem
inferencia condicional a toda informacao disponıvel, Dt. Estes dois ultimos metodos produzem
inferencia completa para os estados (diferentemente do LB, em que tal inferencia resume-se
a primeiros e segundos momentos).
Tabela 3.2: Tempo computacionais necessarios para a convergencia dos diferentes metodos,
em minutos.
LB INLA MCMC
1 15 300
Os tempos computacionais necessarios para a convergencia do metodo MCMC podem ser
bastantes elevados, devido a autocorrelacao dos pararametros com variacao temporal, como se
muestra na tabela (3.2). A adocao do metodo INLA mostra-se, entao, computacionalmente
tao rapido quanto o metodo Linear Bayes, produzindo porem inferencias arbitrariamente
proximas daquelas obtidas via MCMC.
A figura 3.6 e a figura 3.7 apresenta as trajetorias dos valores observados yt (simulados)
junto a resposta media, λt e intervalos de credibilidade a posteriori de 95%, estimados
pelos metodos INLA e MCMC. A figura 3.7 apresenta as trajetorias dos metodos INLA e
MCMC respectivamente. Pode-se observar que a resposta media estimada acompanha bem
as flutuacoes da serie da variavel resposta por ambos metodos e que os valores teoricos da
funcao resposta media parecem ser capturados pelos intervalos de credibilidade a 95%.
A figura 3.8 mostra o histograma das amostras das distribuicoes preditivas um passo a
frente pelos metodos INLA e MCMC, a linha vertical vermelha corresponde ao valor verdadeiro
da variavel resposta. Ja a figura 3.9 apresenta a distribuicao preditiva um passo a frente com
44
0 500 1000 1500
05
10
15
20
25
Tempo
La
mb
da
_t
ArtificialMCMCINLA
Figura 3.6: Media a posteriori e intervalo de credibilidade de 95% para a resposta media λt
usando os metodo INLA (linha azul) e MCMC (linha vermelha) . A linha cinza refere-se aos
valores reais.
0 500 1000 1500
05
1015
2025
Tempo
Lam
bda_
t
INLAICArtificial
0 500 1000 1500
05
1015
2025
Tempo
Lam
bda_
t
MCMCICArtificial
Figura 3.7: Media a posteriori e intervalo de credibilidade de 95% para a resposta media
λt usando os metodos INLA (figura (a)) e MCMC (figura (b)). A linha preta refere-se aos
valores reais.
45
Den
sity
0 5 10 15
0.00
0.05
0.10
0.15
INLAvalor real
Den
sity
0 5 10 15
0.00
0.05
0.10
0.15 MCMC
valor real
Figura 3.8: Histogramas das amostras das distribuicoes preditivas um passo a frente usando
o metodo INLA (esquerda) e MCMC (direita). A linha vertical vermelha refere-se ao valor
real da variavel resposta.
Tempo
0 500 1000 1500
−5
05
1015
20
Figura 3.9: Media a posteriori e intervalo de credibildade para a predicao um passo a frente
junto a serie observada (linha preta) pelo metodo LB.
46
intervalos de credibilidade a 95%, junto aos valores observados para a serie artificial, pelo
metodo LB. Pode-se ver que o modelo foi capaz de capturar bastante bem as flutuacoes da
variavel resposta. O esquema de atualizacao em tempo real do metodo linear bayes permite
incorporacao de nova informacao e previsao a cada instante. Ja os metodos MCMC e INLA
produzem inferencia a posteriori com respeito a toda a serie observada. Os elevados tempos
computacionais em esquemas MCMC dificultam, a atualizacao de informacao para realizacao
de inferencia em tempo real. Adotando-se metodo INLA, entretanto, com baixıssimo tempo
computacional e possıvel incorporar novas informacoes, atualizando-se predicoes e obtendo-se
resultados muito proximos daquele fornecido pelo metodo MCMC, como se ve na figura 3.8.
47
Capıtulo 4
Aplicacoes com Dados Reais
Nesta secao, apresentamos duas aplicacoes de interesse no contexto de epidemiologia
ambiental. Na primeira aplicacao, secao 4.1, tratamos dos efeitos de poluicao atmosferica, em
particular do poluente monoxido de carbono, sobre o numero de obitos de criancas menores
de cinco anos, por doencas respiratorias, em Sao Paulo. Para modelagem do desfecho,
utiliza-se modelo Poisson e Poisson com inflacao de zeros. Como covariaveis, adota-se,
alem do poluente, variaveis climaticas (temperatura, umidade) e de calendario. Uma segunda
aplicacao apresentada na secao 4.2, trata dos nıveis de material particulado (PM10) na cidade
do Rio de Janeiro, usando como regressores variaveis climaticas, como a temperatura mınima
diaria e umidade relativa, com particular atencao centrada sobre o efeito do volume da chuva
sobre tais poluentes. Inicialmente consideramos como resposta o nıvel de material particulado,
descrito por um modelo Gama e em seguida, considera-se resposta dicotomica, com interesse
em analisar fatores associados a ultrapassagem de um limiar de seguranca no nıvel de material
particulado.
4.1 Efeito de Monoxido de Carbono sobre Obitos de
Criancas em Sao Paulo
Apresentamos aqui uma analise da relacao entre o numero diario (contagens) de obitos
de criancas menores de 5 anos por doencas respiratorias, e os nıveis do poluente atmosferico
48
Monoxido de Carbono, CO (ppm) e variaveis climaticas, temperatura mınima (oC) e umidade
relativa media (%), registrados na cidade de Sao Paulo, desde o primeiro de janeiro 1994 a 31
de dezembro 1997, medidos diariamente. Apresenta-se como valor diario a media aritmetica
sobre as oito estacoes de mensuracao disponıveis no banco de dados. Estas informacoes estao
bem documentadas em Alves (2006).
A figura 4.1 mostra a serie temporal de contagens diarias de obitos por doencas
respiratorias para criancas com idades inferiores a 5 anos, na cidade de Sao Paulo, no perıodo
de 1994 a 1997. Observe-se um claro comportamento sazonal com periodo anual e possui
grandes picos nos meses de julho, implicando o incremento de obitos durante o inverno, alem
de uma queda de tendencia do numero de obitos ao longo dos anos.
Tempo
Núm
ero
de Ó
bito
s
0 500 1000 1500
02
46
810
12
Figura 4.1: Contagem diaria de obitos de criancas menores de 5 anos em Sao Paulo nos anos
de 1994 a 1997.
A figura 4.2(a) apresenta a serie temporal de registros do poluente monoxido de carbono
(CO) na cidade de Sao Paulo, no perıodo de 1994 a 1997, observe-se que os meses de
inverno possuem altos picos o que significa alta concentracao de poluentes, eventualmente
ultrapassando o limite aceitavel (Alves 2006), ao igual que na figura 4.1 possui uma queda
de tendencia ao longo dos anos.
49
Tempo
CO
0 500 1000 1500
51
01
5
(a)
Tempo
Te
mp
era
tura
0 500 1000 1500
05
10
15
20
(b)
Tempo
Um
ida
de
0 500 1000 1500
50
60
70
80
90
(c)
Figura 4.2: Registros diario de concentracao de CO, temperatura mınima e umidade media
em Sao Paulo no perıodo 1994 a 1997.
50
A figura 4.2(b) e 4.2(c), exibe as series temporais da temperatura mınima e umidade
media na cidade de Sao Paulo. Note-se que no painel da esquerda (temperatura mımima)
que os meses de inverno possuem picos baixos, associados a maior concentracao da poluicao.
4.1.1 Formulacao do Modelo Proposto
A variavel de interesse e contagem diaria de obitos de criancas menores de 5 anos por
doencas respiratorias na cidade de Sao Paulo, sendo esta discreta com numero medio de
contagens baixo. E usual, no contexto epidemiologico adotar-se o modelo poisson para a
descricao do comportamento probabilıstico desse tipo de variavel. Nosso interesse e analisar
a associacao entre o numero diario de obitos e os nıveis do poluentes atmosfericos monoxido
de carbono (CO), ajustada pelas variaveis climaticas, como a temperatura e umidade, onde
o logaritmo da taxa media da distribuicao e relacionada com cada uma das covariaveis, nıvel
e efeitos sazonais anual (cos( 2πt365
) e sen( 2πt365
)), por meio de Modelos Dinamicos Lineares
Generalizados.
Alem da relacao do numero de obitos com os nıveis dos poluentes atmosfericos no dia
da ocorrencia, considera-se tambem a relacao atraves dos dias passados. A fim de propor
os modelos mais proximos da realidade, optamos por modelos com efeitos das covariaveis
defasados a dias anteriores a ocorrencia do numero de obitos. Portanto considere-se modelos
com efeitos de cada uma das covariaveis desfasados por dois, tres e dois dias anteriores a
observacao do numero de obitos para CO, temperatura e umidade, respectivamente.
Denote-se yt como o numero de obitos de criancas menores de 5 anos, no dia
t, t = 1, · · · , T com distribuicao poisson, yt ∼ Poisson(λt), e as covariaveis nıvel, CO,
temperatura, umidade e sasonalidade. Assim, a forma geral do modelo basico e representado,
Yt | λt ∼ Poisson(λt)
ηt = log(λt), (4.1)
em que o preditor linear ηt associado ao vetor de parametros, assume quinze diferentes
especificacoes, determinadas pela natureza dos parametros que variam no tempo atraves
51
de uma estrutura dinamica e estatica. A evolucao temporal dos parametros e um passeio
aleatorio, o que torna flexıveis os modelos. Em seguida, descrevemos os diferentes modelos.
Modelo 1: Modelo com nıvel dinamico.
log(λt) = ηt = αt + β1COt−2 + β2Tempt−3 + β3Umidt−2 + (4.2a)
β4Cos(2πt
365) + β5Sen(
2πt
365)
αt = αt−1 + vt, vt ∼ N [0, V ] (4.2b)
Modelo 2: Modelo com dinamica no coeficiente de CO.
log(λt) = ηt = α + β1tCOt−2 + β2Tempt−3 + β3Umidt−2 (4.3a)
+β4Cos(2πt
365) + β5Sen(
2πt
365)
β1t = β1,t−1 + ω1t, ω1t ∼ N [0,W1] (4.3b)
Modelo 3: Modelo com dinamica no coeficiente da temperatura.
log(λt) = ηt = α + β1COt−2 + β2tTempt−3 + β3Umidt−2 (4.4a)
+β4Cos(2πt
365) + β5Sen(
2πt
365)
β2t = β2,t−1 + ω2t, ω2t ∼ N [0,W2] (4.4b)
Modelo 4: Modelo com dinamica no coeficiente da umidade.
log(λt) = ηt = α + β1COt−2 + β2Tempt−3 + β3tUmidt−2 (4.5a)
+β4Cos(2πt
365) + β5Sen(
2πt
365)
β3t = β3,t−1 + ω3t, ω3t ∼ N [0,W3] (4.5b)
Modelo 5: Modelo com dinamica no nıvel e o coeficientes de CO.
log(λt) = ηt = αt + β1tCOt−2 + β2Tempt−3 + β3Umidt−2 (4.6a)
+β4Cos(2πt
365) + β5Sen(
2πt
365)
αt = αt−1 + vt, vt ∼ N [0, V ] (4.6b)
β1t = β1,t−1 + ω1t, ω1t ∼ N [0,W1] (4.6c)
52
Modelo 6: Modelo com dinamica no nıvel e o coefieciente da temperatura.
log(λt) = ηt = αt + β1COt−2 + β2tTempt−3 + β3Umidt−2 (4.7a)
+β4Cos(2πt
365) + β5Sen(
2πt
365)
αt = αt−1 + vt, vt ∼ N [0, V ] (4.7b)
β2t = β2,t−1 + ω2t, ω2t ∼ N [0,W2] (4.7c)
Modelo 7: Modelo com dinamica no nıvel e no coeficiente da umidade.
log(λt) = ηt = αt + β1COt−2 + β2Tempt−3 + β3tUmidt−2 (4.8a)
+β4Cos(2πt
365) + β5Sen(
2πt
365)
αt = αt−1 + vt, vt ∼ N [0, V ] (4.8b)
β3t = β3,t−1 + ω3t, ω3t ∼ N [0,W3] (4.8c)
Modelo 8: Modelo com dinamica nos coeficientes do CO e temperatura.
log(λt) = ηt = α + β1tCOt−2 + β2tTempt−3 + β3Umidt−2 (4.9a)
+β4Cos(2πt
365) + β5Sen(
2πt
365)
β1t = β1,t−1 + ω1t, ω1t ∼ N [0,W1] (4.9b)
β2t = β2,t−1 + ω2t, ω2t ∼ N [0,W2] (4.9c)
Modelo 9: Modelo com dinamica nos coeficientes do CO e umidade.
log(λt) = ηt = α + β1tCOt−2 + β2Tempt−3 + β3tUmidt−2 (4.10a)
+β4Cos(2πt
365) + β5Sen(
2πt
365)
β1t = β1,t−1 + ω1t, ω1t ∼ N [0,W1] (4.10b)
β3t = β3,t−1 + ω3t, ω3t ∼ N [0,W3] (4.10c)
Modelo 10: Modelo com dinamica nos coeficientes do temperatura e umidade.
53
log(λt) = ηt = α + βtCOt−2 + β2tTempt−3 + β3tUmidt−2 (4.11a)
+β4Cos(2πt
365) + β5Sen(
2πt
365)
β2t = β2,t−1 + ω2t, ω2t ∼ N [0,W2] (4.11b)
β3t = β3,t−1 + ω3t, ω3t ∼ N [0,W3] (4.11c)
Modelo 11: Modelo com dinamica no nıvel e nos coeficientes de CO e temperatura.
log(λt) = ηt = αt + β1tCOt−2 + β2tTempt−3 + β3Umidt−2 (4.12a)
+β4Cos(2πt
365) + β5Sen(
2πt
365)
αt = αt−1 + vt, vt ∼ N [0, V ] (4.12b)
β1t = β1,t−1 + ω1t, ω1t ∼ N [0,W1] (4.12c)
β2t = β2,t−1 + ω2t, ω2t ∼ N [0,W2] (4.12d)
Modelo 12: Modelo com dinamica no nıvel e nos coeficientes de CO e umidade.
log(λt) = ηt = αt + β1tCOt−2 + β2Tempt−3 + β3tUmidt−2 (4.13a)
+β4Cos(2πt
365) + β5Sen(
2πt
365)
αt = αt−1 + vt, vt ∼ N [0, V ] (4.13b)
β1t = β1,t−1 + ω1t, ω1t ∼ N [0,W1] (4.13c)
β3t = β3,t−1 + ω3t, ω3t ∼ N [0,W3] (4.13d)
Modelo 13: Modelo com dinamica no nıvel e nos coeficientes de temperatura e umidade.
log(λt) = ηt = αt + βtCOt−2 + β2tTempt−3 + β3tUmidt−2 (4.14a)
+β4Cos(2πt
365) + β5Sen(
2πt
365)
αt = αt−1 + vt, vt ∼ N [0, V ] (4.14b)
β2t = β2,t−1 + ω2t, ω2t ∼ N [0,W2] (4.14c)
β3t = β3,t−1 + ω3t, ω3t ∼ N [0,W3] (4.14d)
54
Modelo 14: Modelo com dinamica nos coeficientes de CO, temperatura e umidade.
log(λt) = ηt = α + β1tCOt−2 + β2tTempt−3 + β3tUmidt−2 (4.15a)
+β4Cos(2πt
365) + β5Sen(
2πt
365)
β1t = β1,t−1 + ω1t, ω1t ∼ N [0,W1] (4.15b)
β2t = β2,t−1 + ω2t, ω2t ∼ N [0,W2] (4.15c)
β3t = β3,t−1 + ω3t, ω3t ∼ N [0,W3] (4.15d)
Modelo 15: Modelo com dinamica no nıvel e nos coeficientes de CO, temperatura e
umidade.
log(λt) = ηt = αt + β1tCOt−2 + β2tTempt−3 + β3tUmidt−2 (4.16a)
+β4Cos(2πt
365) + β5Sen(
2πt
365)
αt = αt−1 + vt, vt ∼ N [0, V ] (4.16b)
β1t = β1,t−1 + ω1t, ω1t ∼ N [0,W1] (4.16c)
β2t = β2,t−1 + ω2t, ω2t ∼ N [0,W2] (4.16d)
β3t = β3,t−1 + ω3t, ω3t ∼ N [0,W3] (4.16e)
Com base na estrutura acima, composto de 1458 dados, dos quais 1443 foram utilizados
para estimacao do modelo e os outros 15 dados ficam para a avalicao preditiva.
4.1.2 Inferencia Bayesiana Utilizando INLA
Os modelos acima foram implementados usando o metodo INLA como descrito no capıtulo
3. Para as distribuicoes a priori dos hiperparametros foram usadas prioris nao informativas.
A estrutura dos quinze modelos tem a forma requerida da estrutura do metodo INLA. Nesta
secao descrevemos os modelos 1 e 11. Os outros modelos podem ser feitos por analogia.
Modelo 1: Nıvel Dinamico
O preditor linear ηt descrito pela equacao (4.2a) e com dinamica no nıvel ao longo do
tempo, dado pela equacao (4.2b). Neste caso, as quantidades a serem estimadas sao Ψ =
55
(α1, · · · , αT , β1, · · · , β5, V ), em que o parametro fixo e θ = V e o campo latente e definido
como x = (α1, · · · , αT , β1, · · · , β5). Assuma-se a priori para o hiperparametro θ = V ∼
GI(av, bv), sendo av, bv valores conhecidos. Usando as equacoes (4.2a) e (4.2b) temos que
a distribuicao a priori p(x|θ) e dado por:
p(x|θ) = p(α1)T∏t=2
p(αt|αt−1)5∏i=1
p(βi) (4.17a)
p(x|θ) ∝ exp
{−1
2xTQ(θ)x
}(4.17b)
Com isso temos (x|θ) ∼ N [0,Q(θ)−1], onde Q(θ) e a matriz de precisao, dada por
Q(θ) =
Q(V) 0 0 0 0 0
0 Qβ1β1 0 0 0 0
0 0 Qβ2β2 0 0 0
0 0 0 Qβ3β3 0 0
0 0 0 0 Qβ4β4 0
0 0 0 0 0 Qβ5β5
(T+5)×(T+5)
(4.18)
em que,
Q(V) =
R−11 + V −1 −V −1 0 · · · 0 0
−V −1 2V −1 −V −1 · · · 0 0
0 −V −1 2V −1 · · · 0 0
0 0 −V −1 · · · 0 0...
......
. . ....
...
0 0 0 · · · 2V −1 −V −1
0 0 0 · · · −V −1 V −1
T×T
(4.19)
Neste caso o log da funcao verossimilhanca associada ηt = log λt e:
log p(yt|ηt) ∝ −eηt + ytηt. (4.20)
56
Portanto, a aproximacao Gaussiana para a distribuicao condicional completa xt e dada
por
p(x|θ,y) ∝ p(y|x,θ)p(x|θ) =∏T
t=1 p(yt|xt,θ)p(xt|θ) (4.21)
∝ exp
−12xTQ(θ)x +
∑log p(yt|x,θ)︸ ︷︷ ︸g(xt)=g(ηt)
(4.22)
Modelo 11: Dinamica no nıvel e nos coeficientes de CO e temperatura
O preditor linear ηt descrito pela equacao (4.12a) e com dinamica descrita
pelas equacoes (4.12b) ate (4.12d). Neste caso, as quantidades a serem
estimadas sao Ψ = (α1, · · · , αT , β1,1, · · · , β1,T , β2,1, · · · , β2,T , β3, β4, β5, V,W1,W2), e
o parametro fixo e θ = (V,W1,W2) e o campo latente e definido como x =
(α1, · · · , αT , β1,1, · · · , β1,T , β2,1, · · · , β2,T , β3, β4, β5). Assuma-se que os hiperparametros θ
sao independentes e a distribuicao a priori para V,W1,W2 sera uma GI(av, bv), sendo av, bv
valores conhecidos. Usando as equacoes (4.12a) a (4.12d) temos que a distribuicao a priori
p(x|θ) e dada por:
p(x|θ) =p(α1)T∏t=2
p(αt|αt−1)p(β1,t)T∏t=2
p(β1,t|β1,t−1) (4.23a)
p(β2,t)T∏t=2
p(β2,t|β2,t−1)5∏i=3
p(βi)
p(x|θ) ∝ exp
{−1
2xTQ(θ)x
}(4.23b)
Com isso temos (x|θ) ∼ N [0,Q(θ)−1], onde Q(θ) e a matriz de precisao, dada por
57
Q(θ) =
Q(V) 0 0 0 0 0
0 Q(W1) 0 0 0 0
0 0 Q(W2) 0 0 0
0 0 0 Q(β3β3) 0 0
0 0 0 0 Q(β4β4) 0
0 0 0 0 0 Q(β5β5)
(3T+3)×(3T+3)
(4.24)
em que,
Q(V) =
R−11 + V −1 −V −1 0 · · · 0 0
−V −1 2V −1 −V −1 · · · 0 0
0 −V −1 2V −1 · · · 0 0
0 0 −V −1 · · · 0 0...
......
. . ....
...
0 0 0 · · · 2V −1 −V −1
0 0 0 · · · −V −1 V −1
T×T
(4.25)
e submatriz para para a presicao Wi com i = 1, 2, e
Q(Wi) =
R−12 +W−1
i −W−1i 0 · · · 0 0
−W−1i 2W−1
i −W−1i · · · 0 0
0 −W−1i 2W−1
i · · · 0 0
0 0 −W−1i · · · 0 0
......
.... . .
......
0 0 0 · · · 2W−1i −W−1
i
0 0 0 · · · −W−1i W−1
i
T×T
(4.26)
Neste caso a verossimilhanca associada ηt = log λt e:
log p(yt|ηt) ∝ −eηt + ytηt (4.27)
58
Portanto, a aproximacao Gaussiana para a distribuicao condicional completa xt e dado
por
p(x|θ,y) ∝ p(y|x,θ)p(x|θ) =∏T
t=1 p(yt|xt,θ)p(xt|θ) (4.28)
∝ exp
−12xTQ(θ)x +
∑log p(yt|x,θ)︸ ︷︷ ︸g(xt)=g(ηt)
(4.29)
4.1.3 Escolha do melhor Modelo
A implementacao computacional dos 15 modelos acima, foi feita usando a biblioteca
INLA, disponıvel no pacote R-INLA, na pagina http://www.r-inla.org, que deve ser instalado
no software R.
Para a escolha do melhor modelo foi usado o criterio DIC, utilizados em aplicacoes do
metodo INLA. O modelo que apresenta melhor ajuste ao dados sao aqueles que presentam
o menor valor. Na tabela (4.1), sao apresentados os valores obtidos para o criterio de
comparacao, DIC de cada um dos modelos estimados.
O melhor modelo, segundo o criterio DIC, e o modelo 1, com dinamica no nıvel.
59
Tabela 4.1: Comparacao de cada um dos modelos (especificacoes do preditor) ajustados para
o numero de obitos
Modelo Dinamica DIC
1 Nıvel 5331.4
2 Poluente 5441.6
3 Temperatura 5415.6
4 Umidade 5467.1
5 Nivel e Poluente 5337.1
6 Nivel e Temperatura 5338.9
7 Nivel e Umidade 5341.1
8 Poluente e Temperatura 5393.8
9 Poluente e Umidade 5434.7
10 Temperatura e Umidade 5415.1
11 Nivel, Poluente e temperatura 5344.4
12 Nivel, Poluente e Umidade 5347.6
13 Nivel, Temperatura e Umidade 5348.2
14 Poluente, Temperatura e Umidade 5408.4
15 Nivel,Poluente, Temperatura e Umidade 5354.0
60
4.1.4 Resultados para o Modelo 1: Dinamica no Nıvel
Como ja foi mencionado a implementacao deste modelo, preditor linear com dinamica no
nıvel, foi feita pela biblioteca INLA e apresentada no Apendice B.1.
Para a comparacao do metodo INLA com MCMC, o modelo com dinamica no nıvel,
tambem foi estimado via MCMC, como descrito no capıtulo (2) utilizando o software
WinBugs. Amostras a posterior no MCMC foram simuladas a partir de duas cadeias de
33000 iteracoes, sendo descartadas as primeiras 6000, denominado perıodo de aquecimento,
logo foi considerado um espacamento de 25 iteracoes. Assim, 1080 amostras foram deixados
para cada cadeia. Portanto, a amostra final e de 2160 para a obtencao das estimativas da
distribuicao a posteriori.
Tabela 4.2: Estimativas das distribuicoes a posteriori com intervalos de credibilidade a 95%
dos parametros estaticos.
INLA MCMC
ICI media ICS ICI media ICS
β1 0.0142154 0.0525784 0.0905837 0.013688 0.0531 0.091782
β2 -0.1296084 -0.0746466 -0.0191963 -0.1323025 -0.074886 -0.017695
β3 -0.0297995 0.0077015 0.0453201 -0.0299445 0.0077052 0.0447915
β4 -0.2833013 -0.1713995 -0.0596559 -0.2949075 -0.1741129 -0.0689595
β5 -0.043985 0.0648599 0.1731649 -0.0278118 0.0780079 0.1728275
A tabela 4.2 mostra os resultados das estatısticas das distribuicoes a posteriori para os
parametros estaticos para o modelo Poisson com nıvel variando no tempo (modelo 1), e
figura 4.3 mostra os graficos das distribuicoes a posteriori para os parametros estaticos,
estimados pelos metodos INLA e MCMC, mostrando resultado muito semelhante com as duas
abordagens, mas o tempo computacional no INLA foi muito rapido (segundos) em compracao
a MCMC. Os coeficientes β1, β2, β4 sao significativamente diferentes de zero ao nıvel del 95%
de credibilidade. A estimativa para o coeficiente (β1) do monoxido de carbono (CO) e
significativo com efeito positivo sobre a taxa de obitos. O coeficiente (β2) da temperatura
mımima e significativa com efeito negativo sobre a variavel resposta, isto e, temperaturas
61
Bet
a1
−0.02 0.02 0.04 0.06 0.08 0.10 0.12
05
1015
20
INLAMCMC
Bet
a2
−0.15 −0.10 −0.05 0.00
02
46
810
1214 INLA
MCMCB
eta3
−0.06 −0.02 0.02 0.04 0.06 0.08
05
1015
20
INLAMCMC
Bet
a4
−0.35 −0.25 −0.15 −0.05
01
23
45
67 INLA
MCMC
Bet
a5
−0.1 0.0 0.1 0.2
02
46
INLAMCMC
Figura 4.3: Aproximacoes das distribuicoes a posteriori e intervalos de credibilidade a 95% para
os parametros estaticos do modelo Poisson com nıvel variando no tempo estimados pelos metodos
INLA (linha azul) e MCMC (histograma, cor preto). A linha azul vertical pontilhada e os pontos
pretos sao os intervalo de credibilidade do INLA e MCMC, respectivamente.
62
baixas estao associadas a maior numero de obitos, como ja esperado. Apesar de a estimativa
do coeficiente (β3) da umidade nao ser significativo, pode-se considerar com efeito positivo
com probabilidade a posteriori de aproximadamente de 91.2%.
Para se obter o risco relativo da estimativa obtida para o monoxido de carbono, basta
exponenciar os valores da a posteriori para β1, obtendo-se 1.054 o que siginifica que o
risco relativo associado a uma elevacao de um desvio padrao no nıvel de CO provocaria
um acrescimo de 5.4% na taxa de obitos; sendo o intervalo de credibilidade de 95%,
[1.014317, 1.094813].
0 500 1000 1500
0.2
0.4
0.6
0.8
1.0
1.2
1.4
Tempo
Niv
el_
t
INLAMCMC
Figura 4.4: Media a posteriori e intervalos de credibilidade de 95% para a componente αt (Nıvel).
A linha azul refera-se ao metodo INLA e a vermelha a MCMC.
A Figura 4.4 apresenta o nıvel do preditor linear, pelos metodos INLA e MCMC,
acompanham os valores observados da resposta. Pode-se observar que o nıvel capta a
tendencia da variavel resposta (numero de obitos) ao longo dos anos.
A figura 4.5 mostra a aproximacao da distribuicao a posteriori da marginal V −1, e os
resultados sao apresentados na tabela 4.3, observa-se que ambos metodos fornecem resultados
equivalentes.
63
1/V
Den
sity
500 1000 1500 2000
0.00
000.
0005
0.00
100.
0015
INLAMCMC
Figura 4.5: Distribuicao a posteriori marginal para o hiperparametro V −1 e intervalos de
credibilidade de 95%. A linha azul refera-se ao metodo INLA e o Histograma ao metodo MCMC.
Tabela 4.3: Estimativas das distribuicoes a posteriori com intervalos de credibilidade a 95% do
hiperparametro V.
INLA MCMC
ICI media ICS ICI media ICS
V −1 556.4787 1012.709 1679.086 585.185 998.0053 1559
A figura 4.6 mostra a serie do numero de obitos ao lado da serie da resposta media, λt
(taxa de obitos) pelos dois metodos INLA e MCMC. Observa-se que o modelo acompanha o
comportamento do numero de obitos.
A distribuicao preditiva nos fornece as previsoes para horizontes futuros. A figura 4.7
64
0 500 1000 1500
02
46
810
12
Tempo
Lam
bda_
t
ObservadoMCMC e ICINLA e IC
Figura 4.6: Media a posteriori para a resposta media, λt com intervalo de credibilidade de 95%.
Numero de obitos (linha preta). Refera-se ao metodo INLA (linha azul) e MCMC (linha vermelha).
mostra o grafico de barras das amostras das distribuicoes preditiva dos 15 dias finais que
foram reservados, obtidas pelo metodo INLA. As linhas verticais vermelhas de cada grafico
de barras indicam o numero de obitos observados em cada um desses dias. Intervalos de
predicao a 95% capturam os reais valores observados.
Por outra lado, observou-se que os dados contem varios dias com contagens nulas
e a adocao de um modelo Poisson pode ser questionada. Portanto, uma alternativa e
realizar um ajuste que contemple o excesso de zeros, sendo comum utilizar o modelo
com distribuicao Poisson inflacionada de zeros (ZIP) (Fernandes 2006). Na abordagem
INLA temos duas propostas da distribuicao Poisson inflacionada de zeros (ZIP): inflacao
de zeros de tipo 0 e inflacao de zeros de tipo 1 (http://www.math.ntnu.no/inla/r-
inla.org/doc/likelihood/zeroinflated.pdf.), dada por
tipo 0: p(y|...) = p× 1y=0 + (1− p)Poisson (y|y > 0).
tipo 1: p(y|...) = p× 1y=0 + (1− p)Poisson(y).
65
0 2 4 6
1
0100
200300
0 2 4 6
2
0100
200300
0 2 4 6
3
050
100150
200250
300350
0 2 4 6
4
050
100150
200250
300350
0 2 4 6
5
050
100150
200250
300350
0 2 4 6
6
050
100150
200250
300350
0 2 4 6
7
050
100150
200250
300350
0 2 4 6
8
050
100150
200250
300
0 2 4 6
9
050
100150
200250
300350
0 2 4 7
10
050
100150
200250
300350
0 2 4 6
11
050
100150
200250
300350
0 2 4 6
12
050
100150
200250
300350
0 2 4 6
13
050
100150
200250
300350
0 2 4 6
14
050
100150
200250
300
0 2 4 6
15
050
100150
200250
300350
Figura 4.7: Grafico de barra das amostras das ditribuicoes preditivas para o numero de obitos, para
dados reais de acordo com o modelo Poisson, com horizontes variando de 1 a 15
66
A implementacao destas duas propostas de inflacao de zeros para o melhor modelo do
contagens de obitos (modelo 1: Nıvel mudando no tempo), sao apresentadas no Apendice
B.2 e B.3.
Para comparar estas tres famılias (Poisson tradicional, distribuicao Poisson inflacionada
de zeros (ZIP): tipo 0 e tipo 1) usamos o erro quadratico medio (EQM) de predicao, para
os ultimos 15 dias disponıveis. Os resultados do erro quadratico medio associado ao modelo
Poisson tradicional e de 18.67, correspondente a um erro relativo medio de 27.8%. O erro
quadratico medio do modelo com inflacao de zeros do tipo 0 resultou em 19.22 e erro relativo
medio de 30%. Ja erro quadratico medio para inflacao de zeros do tipo 1 e = 18.70 e
erro relativo medio 27.98%. Como se pode ver os EQM preditivos para as tres abordagens
apontam uma ligeira vantagem do modelo Poisson tradicional, em relacao aos demais.
67
4.2 Efeito de Chuva sobre Nıveis de material Particulado
no Rio de Janeiro
Pode haver interesse, no contexto de epidemiologia ambiental, em modelos para explicacao
de nıveis de poluicao, ja que frequentemente ha observacoes faltantes nesse tipo de dado. Tais
modelos podem ser utilizados como um nıvel hierarquico para modelos que objetivem explicar
contagens epidemiologicas. Ainda, pode-se ter interesse na modelagem de desfechos contınuos
e estritamente positivo (por exemplo, nıvel de um poluente) ou binarios (ultrapassagem ou nao
de certo limite de seguranca para a saude). Todas essas aplicacoes podem ser acomodadas
pela classe de modelos dinamicos lineares generalizados e estas sob a forma de um modelo
de campo gaussiano markoviano latente.
Portanto, a metodologia INLA pode-se mostrar bastante util a implementacao de
inferencia bayessiana nesse contexto. Apresentamos a seguir, a tıtulo de ilustracao dois
modelos: uma em que se modela o nıvel de certo poluente atmosferico em funcao de variaveis
climaticas e de calendario e, em seguida, um segundo modelo em que adotamos como resposta
a variavel indicadora de que o nıvel desse poluente ultrapasse um limite de seguranca para a
saude.
4.2.1 Descricao dos Dados
Apresentamos aqui um analise da relacao entre o nıveis do poluente material particulado
(PM10) e variaveis climaticas, temperatura mınima e umidade media e ainda utilizam-se
variaveis indicadoras de dia da semana, tomando como base o domingo, e tendo especial
atencao a propagacao do efeito de volumes de chuva sobre os nıveis do poluente, registrados
na cidade do Rio de Janeiro, desde o primeiro de setembro de 2000 a 31 de agosto de 2002,
medidos diariamente. Apresenta-se como valor diario a media aritmetica sobre as diferentes
estacoes de mensuracao disponıveis no banco de dados. Segundo informacoes colhidas no sıtio
da Companhia de Tecnologia de Saneamento Ambiental (CETESB), o material particulado e
um conjunto de poluentes constituıdos de poeiras, fumacas e todo tipo de material solido e
lıquido que se mantem suspenso na atmosfera por causa de seu pequeno tamanho. O tamanho
68
das partıculas esta diretamente associado ao seu potencial para causar problemas a saude,
sendo que quanto menores, maiores os efeitos provocados. As partıculas inalaveis (PM10)
sao definidas como aquelas cujo diametro aerodinamico e menor que 10µm. Dependendo
da distribuicao de tamanho (0 a 10 µm), podem ficar retidas na parte superior do sistema
respiratorio ou penetrar mais profundamente, alcancando os alveolos pulmonares. Nesta
aplicacao consideram-se o material particulado, PM10.
Na figura 4.8 exibe-se a serie temporal das concentracoes medias diarias de PM10 (µg/m3)
e as outras tres series mostram volume de chuva, temperatura mınima e umidade media, na
cidade de Rio de Janeiro, no perıodo de primeiro de setembro de 2000 ao 31 de agosto de
2002.
69
0 100 200 300 400 500 600 700
2040
6080
100
120
140
Tempo
PM
0 100 200 300 400 500 600 700
020
4060
80
Tempo
Chuv
a
0 100 200 300 400 500 600 700
1618
2022
2426
Tempo
Temp
eratur
a míni
ma
0 100 200 300 400 500 600 700
6065
7075
8085
9095
Tempo
Umida
de m
édia
Figura 4.8: Serie de Nıveis diarios de PM10(µg/m3) e variaveis meteorologicas no Rio de Janeiro,
no perıodo de Setembro de 2000 a Agosto de 2002.70
4.2.2 Modelo Gama
Denote-se por Yt o nıvel de material particulado (PM10) como a variavel a ser explicada
no dia t, t = 1, · · · , T com distribuicao Gama, Yt ∼ Gama(ϕ, λt) e como regressoras,
nıvel padronizado de chuva, nıvel padronizado de temperatura mınima, nıvel padronizado de
umidade e indicadoras de dia da semama e feriado.
O preditor linear ηt associado ao vetor de parametros com nıvel dinamico e parametros
estaticos. Para evolucao temporal do nıvel, assume-se um passeio aleatorio. Portanto, o
modelo ajustado segue a seguinte estrutura:
Yt ∼ Gama(ϕ, λt) (4.30a)
log(ϕ
λt) = ηt = αt + β1Chuva+ β2Temp+ β3Umid+ β4Seg + β5Ter +
β6Qua+ β7Qui+ β8Sex+ β9Sab+ β10Fer (4.30b)
em que,
αt = αt−1 + vt, vt ∼ N [0,V] (4.31)
4.2.3 Resultados do modelo
A implementacao deste modelo, preditor linear com dinamica no nıvel, foi feita pela
biblioteca INLA e apresentada no Apendice C.1.
A figura 4.9 mostra as distribuicoes a posterioris marginais dos parametros estaticos para o
modelo Gama com nıvel variando no tempo, e tabela 4.4 mostra os resultados das estatısticas
dos parametros estaticos.
Observe-se que todas as estimativas para os coeficientes sao significativamente diferentes
de zero. As estimativas dos dias da semana crescem a medida em que se aproxima a sexta-
feira e decresce no sabado, como se exibe na figura 4.10, isto pode ser pelo trafego de veıculos
e atividades industriais nos dias da semana. Ja feriados tem efeito negativo, provavelmente
devido a reducao de atividades industriais e trafego de veıculos. O efeito da chuva (β1), e
negativo, o que indica decrescimo nos nıveis de PM10 em funcao de aumento nos volumes
de chuva. A chuva ”limpa a atmosfera”, ao depositar partıculas no solo, entao em dias de
chuva e dias subsequentes, espera-se menos material particulado no ar. A temperatura media
apresenta efeito positivo (β2) e a estimativa do efeito da umidade media (β3) e negativo.
71
−0.06 −0.04 −0.02 0.00
010
2030
40
Beta.Chuva
0.05 0.10 0.15 0.20
05
1015
20
Beta.Temp
−0.12 −0.10 −0.08 −0.06 −0.04 −0.02
010
2030
40
Beta.Umid
0.00 0.05 0.10 0.15 0.20
05
1015
Beta.Seg
0.05 0.10 0.15 0.20 0.25
05
1015
Beta.Ter
0.10 0.15 0.20 0.25 0.30
05
1015
Beta.Qua
0.10 0.15 0.20 0.25 0.30
05
1015
Beta.Qui
0.10 0.15 0.20 0.25 0.30 0.35
05
1015
Beta.Sex
0.00 0.05 0.10 0.15 0.20
05
1015
Beta.Sab
−0.25 −0.15 −0.05 0.05
02
46
810
Beta.Feriado
Figura 4.9: Aproximacoes das distribuicoes a posteriori e intervalos de credibilidade a 95% para os
parametros estaticos do modelo Gama com nıvel variando no tempo. A linha vertical pontilhada
sao os intervalo de credibilidade.
72
Tabela 4.4: Estimativas das distribuicoes a posteriori com intervalos de credibilidade a 95%
dos parametros estaticos.
Media e.p ICI Mediana ICS
ϕ 35.9606 4.29025 28.0649 35.7944 44.9050
β1 -0.0257 0.0083 -0.0419 -0.0257 -0.0093
β2 0.1321 0.0204 0.0918 0.1321 0.1720
β3 -0.0717 0.0101 -0.0916 -0.0717 -0.0519
β4 0.1048 0.0241 0.0574 0.1048 0.1522
β5 0.1591 0.0254 0.1092 0.1591 0.2090
β6 0.1922 0.0258 0.1415 0.1922 0.2429
β7 0.1980 0.0258 0.1473 0.1980 0.2486
β8 0.2208 0.0254 0.1708 0.2208 0.2707
β9 0.1053 0.0242 0.0579 0.1053 0.1527
β10 -0.0898 0.0360 -0.1601 -0.0899 -0.0189
−0
.10
.00
.10
.20
.3
Se
g
Ter
Qu
a
Qu
i
Sex
Sa
b
Fe
r
Figura 4.10: Media a posteriori e intervalos de credibilidade de 95% para os parametros estaticos
dos dias da semana comparados ao domingo do modelo Gama
73
O risco relativo associado a variacoes no nıvel da chuva e 0.98 o que siginifica que o risco
relativo associado a um aumento de um desvio padrao no nıvel de chuva provocaria uma
reducao de 2% no nıvel de material particulado; sendo o intervalo de credibilidade de 95%,
[0.9589657, 0.9907431].
0 100 200 300 400 500 600 700
3.0
3.5
4.0
4.5
Tempo
Niv
el(a
lpha
_t)
Figura 4.11: Media a posteriori e intervalos de credibilidade de 95% para a componente αt (Nıvel).
100 200 300 400
0.000
0.004
0.008
0.012
1/V
Figura 4.12: Distribuicao a posteriori marginal para o hiperparametro V −1 e intervalos de
credibilidade de 95%.
74
A Figura 4.11 apresenta o nıvel do preditor linear. Pode-se observar que a serie capta a
tendencia da variavel resposta ao longos dos anos. A figura 4.12 mostra a aproximacao
da distribuicao a posteriori da marginal V −1. A tabela 4.5 apresenta as estatısticas do
hiperparametro.
Tabela 4.5: Estimativas das distribuicoes a posteriori com intervalos de credibilidade a 95% do
hiperparametro V −1.
Media e.p. ICI Mediana ICS
V −1 117.78621 32.704912 68.14643 112.81178 195.49330
A figura 4.13 mostra a serie do material particulado ao lado da serie da resposta
media µt, com intervalo de credibilidade de 95%. Observe-se que o modelo acompanha
o comportamento da serie do PM10.
0 100 200 300 400 500 600 700
20
40
60
80
10
01
20
14
0
Tempo
y vs
mu
_t
Figura 4.13: Nıveis do material particulado (lihna cinza) e a funcao da resposta media, µt (linha
preta) com intervalo de credibilidade de 95% (linha tracejada).
A figura 4.14 mostra os histogramas das amostras das distribuicoes preditiva para os nıveis
de material particulado dos 15 dias finais do perıodo. As linhas verticais vermelhas de cada
histograma indicam os valores reais observados em cada um desses dias. Pode-se observar
que quase todas as distribuicoes preditivas contemplam os valores observados.
75
1
Density
40 80 120
0.0000.005
0.0100.015
0.0200.025
2
Density
20 60 100
0.0000.005
0.0100.015
0.0200.025
3
Density
40 80 140
0.0000.005
0.0100.015
0.0200.025
4
Density
50 150
0.0000.005
0.0100.015
5
Density
50 150
0.0000.005
0.0100.015
6
Density
50 150
0.0000.005
0.0100.015
7
Density
50 150
0.0000.005
0.0100.015
8
Density
50 150
0.0000.005
0.0100.015
9
Density
50 150
0.0000.005
0.0100.015
0.020
10
Density
0 100 200
0.0000.005
0.0100.015
11
Density
0 100 200
0.0000.005
0.0100.015
12
Density
50 150 250
0.0000.002
0.0040.006
0.0080.010
0.012
13
Density
50 150 250
0.0000.002
0.0040.006
0.0080.010
0.012
14
Density
50 200
0.0000.002
0.0040.006
0.0080.010
0.012
15
Density
0 100 200
0.0000.005
0.0100.015
Figura 4.14: Histogramas das amostras das ditribuicoes preditivas para nıvel de material particulado
do modelo Gama, com horizontes variando de 1 a 15. As linhas verticais sao os valores reais
observado.
76
4.2.4 Modelo Bernoulli
Denote-se Yt a variavel a ser explicada como a indicadora de que o nıvel de material
particulado ultrapassou o limite recomendado pela CETESB como seguro a saude (PM10 >
50µg/m3); isto e que a qualidade do ar e regular, o que significa que pessoas de grupos
sensıveis (criancas, idosos e pessoas com doencas respiratorias e cardıacas) podem apresentar
sintomas como tosse seca e cansaco. A populacao, em geral, nao e afetada.
Portanto, Yt ∼ Bernoulli(pt), com t, t = 1, · · · , T e como regressoras: nıvel padronizado
de chuva, nıvel padronizado de temperatura mınima, nıvel padronizado de umidade e
indicadoras de dia da semama e feriado. O preditor linear ηt esta associado ao vetor de
parametros com nıvel dinamico e demais parametros estaticos.
Portanto o modelo ajustado segue a seguinte estrutura:
Yt ∼ Binomial(1, pt) (4.32a)
log(pt
1− pt) = ηt = αt + β1Chuva+ β2Temp+ β3Umid+ β4Seg + β5Ter +
β6Qua+ β7Qui+ β8Sex+ β9Sab+ β10Fer (4.32b)
em que a evolucao temporal do nıvel e um passeio aleatorio.
αt = αt−1 + vt, vt ∼ N [0,V] (4.33)
4.2.5 Resultados do modelo
Para a implementacao do modelo, usamos a biblioteca INLA e apresentada no Apendice
C.2. A figura 4.15 mostra as distribuicoes a posteriori marginais dos parametros estaticos
para o modelo Bernoulli com nıvel variando no tempo, e tabela 4.6 mostra os resultados das
estatısticas dos parametros estaticos.
Observa-se que os efeitos dos parametros das distribuicoes a posteriori sao
significativamente diferentes de zero; exceto as segundas-feira (β4) e sabados (β9).
Novamente as estimativas dos dias da semana crescem a medida em que se aproxima a
sexta-feira e decrescem no sabado, como se exibe na figura 4.16, isto pode ser pelas diferentes
atividades industriais, trafego de veıculos nos dias da semana, ja o feriado tem efeito negativo,
77
−1.0 −0.5 0.0 0.5
0.00.5
1.01.5
2.02.5
3.0
Beta.Chuva
0.5 1.0 1.5 2.0
0.00.5
1.01.5
2.02.5
Beta.Temp
−1.4 −1.2 −1.0 −0.8 −0.6 −0.4 −0.2
01
23
4
Beta.Umid
−1 0 1 2
0.00.2
0.40.6
0.81.0
Beta.Seg
−1 0 1 2 3
0.00.2
0.40.6
0.81.0
Beta.Ter
−1 0 1 2 3
0.00.2
0.40.6
0.81.0
Beta.Qua
0 1 2 3
0.00.2
0.40.6
0.81.0
Beta.Qui
0 1 2 3 4
0.00.2
0.40.6
0.81.0
Beta.Sex
−1 0 1 2
0.00.2
0.40.6
0.81.0
1.2
Beta.Sab
−4 −3 −2 −1 0 1 2
0.00.2
0.40.6
0.8
Beta.Feriado
Figura 4.15: Aproximacoes das distribuicoes a posteriori e intervalos de credibilidade a 95% para os
parametros estaticos do modelo Bernoulli com nıvel variando no tempo. A linha vertical pontilhada
sao os intervalo de credibilidade.
78
Tabela 4.6: Estimativas das distribuicoes a posteriori com intervalos de credibilidade a 95%
dos parametros estaticos para o modelo Bernoulli.
Media Sd ICI Mediana ICS
β1 -0.2888031 0.1310572 -0.55602824 -0.2853235 -0.04103419
β2 1.1992542 0.1471949 0.92054531 1.1957534 1.49772529
β3 -0.8431380 0.1032854 -1.05254980 -0.8407065 -0.64733609
β4 0.5775461 0.3468859 -0.09876125 0.5760566 1.26239545
β5 1.2286691 0.3588336 0.53377010 1.2253947 1.94213943
β6 1.2314627 0.3617086 0.53084082 1.2282191 1.95048542
β7 1.3654385 0.3575843 0.67317938 1.3621016 2.07661565
β8 1.5402010 0.3644277 0.83574698 1.5364285 2.26601823
β9 0.6257250 0.3406822 -0.03849647 0.6242737 1.29827086
β10 -1.0257493 0.4570002 -1.93256933 -1.0221242 -0.13972123
pois tais atividades diminuem. O efeito da chuva (β1) e negativo, o que indica decrescimo na
probabilidade de ultrapassagem do nıvel de seguranca de PM10 em funcao de aumentos dos
volumenes de chuva.
A estimativa da razao de chances associado a um incremento c sobre os nıveis da chuva,
(β1), segue a razao de chances associado a um desvio padrao (c = 1) no nıvel de chuva
provocaria uma reducao de 24.34% no nıvel na chance; com intervalo de credibilidade de
95%, [0.5734823, 0.9597963].
A Figura 4.17 apresenta o nıvel do preditor linear. Pode ser observado que a serie
acompanha a tendencia da funcao de resposta media ao longos dos anos. A figura 4.18
mostra a aproximacao da distribuicao a posteriori da marginal do hiperparametro (V −1). A
tabela 4.7 apresentam as estatısticas do hiperparametros.
A figura 4.19 apresenta a resposta media do material particulado com intervalo de
credibilidade de 95%. Observa-se que o modelo acompanha o comportamento da serie do
PM10. Adotando-se o modelo ajustado e classificando-se como dia de ultrapassagem do
limite de seguranca aquele em que pt > 0.50, pt a media a posteriori de pt. Comparando-se
tal classificacao a real situacao de cada dia, pode-se avaliar a sensibilidade e a especificidade
79
−2−1
01
23
Seg
Ter
Qua
Qui
Sex
Sab
Fer
Figura 4.16: Media a posteriori e intervalos de credibilidade de 95% para os parametros estaticos
dos dias da semana comparados ao domingo do modelo Bernoulli
Tabela 4.7: Estimativas das distribuicoes a posteriori com intervalos de credibilidade a 95% do
hiperparametro V −1
Media e.p. ICI Mediana ICS
V −1 11.15328 6.020395 3.637271 9.812335 26.44925
80
0 100 200 300 400 500 600 700
−4−2
02
4
Tempo
Niv
el
Figura 4.17: Media a posteriori e intervalos de credibilidade de 95% para a componente αt (Nıvel).
0 10 20 30 40
0.00
0.02
0.04
0.06
0.08
1/V
Figura 4.18: Distribuicao a posteriori marginal para o hiperparametro V −1 e intervalos de
credibilidade de 95%.
81
do procedimento de classificacao. Os erros de classificacao sao exibidos na figura 4.20. A
sensibilidade e a probabilidade de classificar como ultrapassagem, dado que ultrapassou e foi
igual a 93% e especifidade e definida como classificacao como nao ultrapassagem, dado que
nao ultrapassou com 87.5%. A taxa glovbal de acertos foi de 80%.
0 100 200 300 400 500 600 700
0.0
0.2
0.4
0.6
0.8
1.0
Tempo
p_t
Figura 4.19: Estimativa da media a posteriori (linha preta) com intervalo de credibilidade de 95%
(linha tracejada) para o modelo Bernoulli.
A figura 4.22 mostra o grafico de barras das amostras das distribuicoes preditiva para os
nıveis de material particulado dos 20 dias finais do perıodo de analises. As linhas verticais
vermelhas de cada uma das barras indicam os valores reais observados em cada um desses
dias, percebe-se que quase todas as distribuicoes preditivas estam centradas em torno dos
valores observados. A figura 4.21 mostra a media a posteriori e intervalos de credibilidade de
95% da distribuicao preditiva.
82
0 100 200 300 400 500 600 700
−1.
0−
0.5
0.0
0.5
1.0
Tempo
erro
Figura 4.20: Erro de classificacao do modelo Bernoulli, para cada dia do perıodo de observacao
690 700 710 720 730
0.2
0.4
0.6
0.8
1.0
Tempo
p_t
Figura 4.21: Intervalo de credibilidade de 95% (linha tracejada vermelha) e media a posteriori da
distribuicao preditiva (linha vermelha). A linha cinza refere-se aos valores observados.
83
0 1
1
0200
400600
800
0 1
2
0200
400600
8001000
0 1
3
0200
400600
8001000
0 1
4
0200
400600
800
0 1
5
0200
400600
800
0 1
6
0200
400600
800
0 1
7
0200
400600
800
0 1
8
0200
400600
800
0 1
9
0200
400600
8001000
0 1
10
0200
400600
800
0 1
11
0200
400600
800
0 1
12
0200
400600
800
0 1
13
0200
400600
800
0 1
14
0200
400600
0 1
15
0200
400600
800
0 1
16
0200
400600
800
0 1
17
0200
400600
800
0 1
18
0200
400600
8001000
0 1
19
0200
400600
800
0 1
20
0100
200300
400500
600700
Figura 4.22: Modelo Bernoulli: Grafico de barras das amostras das distribuicoes preditivas do
ultrapassagem do limiar de seguranca de PM10, com horizontes variando de 1 a 20. As linhas
verticais vermelhas sao os valores reais observado.
84
Capıtulo 5
Conclusoes
Esta dissertacao comeca descrevendo metodos de aproximacao para realizar inferencia
Bayesiana de modelos dinamicos com estrutura observacional pertencente a familia
exponencial. Dentre os metodos de aproximacao de inferencia, estao os metodos de
aproximacao LB, MCMC e INLA, descrito no Capıtulo 2 e Capıtulo 3. Esta dissertacao
foca na comparacao destes metodos.
Uma ilustracao de aplicacao dos metodos de aproximacao: Linear Bayes, MCMC e INLA
em um modelo Poisson dinamico com dados simulados foi apresentada. Uma questao
importante da comparacao entre os tres metodos e que de fato o modelo ajustado pelo
metodo LB difere do modelo ajustado por MCMC e INLA, uma vez que, neste dois ultimos
metodos, pressupos-se variancia fixas no tempo e no metodo LB em contrapartida, tem-se,
variando no tempo e especificadas por meio de fator de desconto. Outra diferenca essencial
entre os metodos deve-se ao fato de que o LB processa informacao em tempo real, produzindo
inferencia a cada tempo t condicionalmente a informacao disponıvel ate esse instante, Dt. Ja
INLA e MCMC produzem inferencia condicional a toda informacao disponıvel, Dt. Embora
estes dois ultimos metodos produzam inferencia completa para os estados (diferentemente do
LB, em que tal inferencia resume-se a primeiros e segundos momentos).
Os tempos computacionais necessarios para a convergencia do metodo MCMC podem
ser bastante elevados, devido a autocorrelacao dos pararametros com variacao temporal. A
adocao do metodo INLA mostra-se, entao, computacionalmente rapido quanto o metodo
Linear Bayes, produzindo porem inferencias arbitrariamente proximas daquelas obtidas via
85
MCMC. Este ultimo metodo continua sendo um abordagem mais flexıvel, praticamente
qualquer modelo pode ser ajustado usando MCMC, alem disso a documentacao e exemplos
e mais extensa para os modelos.
O Capıtulo 4 apresenta duas aplicacoes com dados reais usando os metodos MCMC e
INLA, a primeira aplicacao, trata do efeito do poluentes sobre o numero de obitos de criancas
menores de cinco anos por doencas respiratorias, na cidade de Sao Paulo e na segunda trata
dos nıveis de material particulado em funcao de variaveis climaticas e de calendario, na cidade
do Rio de Janeiro, adotando-se os modelos Gama e Bernoulli. Esta ultima com interesse em
analisar fatores associados a ultrapassagem de um limiar de seguranca no nıvel de material
particulado.
Como trabalhos futuros, poderıamos considerar a modelagem da propagacao de efeitos de
uma regressora sobre a resposta esperada ou, ainda, considerar estruturas espaco-temporais
no preditor.
86
Apendice A
Codigos Usados para Dados Artificiais
A seguir, tem-se os codigo utilizados pelos metodos de aproximacao para dados artificiais
para ajustar o modelo com dois covariaveis para dados simulados no capıtulo 3.
A.1 Codigo usado pelo metodo Linear Bayes
#########################################################
#modelo :
# y_t|lambda_t ~ poisson(lambda_t)
# eta_t=log(lambda_t)=alpha_t + beta_t*CO, F_t=(1,CO)
# alpha_t = alpha_{t-1} + v_t, v_t ~ (0, V_t)
# beta_t = beta_{t-1} + w_t,w_t ~ (0, W_t)
#########################################################
phi<-1 #parametro de escala da familia exponencial
d=2 #ncol(matrizdados)-1 # 5 dim. do vetor parametrico
T=n #nrow(matrizdados) #1458 comprimento da serie temporal
yt<-y
G = diag(d)
m0 = matrix(ncol=1, rep(0,d))
C0 =2*diag(d) # 2*diag(d)
vetdes<-c(0.95,0.97) # c(0.92,0.97)
87
B<-diag(1/sqrt(vetdes))
#Inicializando a dimens~ao das matrizes
mt=matrix(ncol=T,nrow=d)
Rt=array(rep(diag(d),T),dim=c(d,d,T))
ft=matrix(nrow=T,ncol=1)
pred=icl.pred=icu.pred=ft
pred.aj=hh=icl=icu=qt=rt=rtp=gtp=zt=ztp=gt=pt=ft
at=St=mt
Ct=Rt
# et=ft
at[,1]=G%*%m0
Rt[,,1]=B%*%G*C0*t(G)%*%B
ft[1]=t(F[1,])%*%at[,1]
St[,1]=Rt[,,1]%*%F[1,]
qt[1]=t(F[1,])%*%St[,1]
rt[1]=1/qt[1]
zt[1]=exp(-ft[1])/qt[1]
# a posteriori
rtp[1]=rt[1]+yt[1]
ztp[1]=zt[1]+1
gt[1]=log(rtp[1]/ztp[1])+1/(2*rtp[1])
pt[1]=(2*rtp[1]-1)/(2*rtp[1]^2)
# atualizac~ao
mt[,1]=at[,1]+St[,1]*(gt[1]-ft[1])/qt[1]
Ct[,,1]=Rt[,,1]-St[,1]%*%t(St[,1])*(1-pt[1]/qt[1])/qt[1]
pred[1]<-rt[1]/ zt[1]
var.pred<-rt[1]*(zt[1]+1)/(zt[1])^2
icl.pred[1]<-pred[1]-2*sqrt(var.pred)
icu.pred[1]<-pred[1]+2*sqrt(var.pred)
for(t in 2:T){
88
at[,t]=G%*%mt[,t-1]
Rt[,,t]=B%*%G*Ct[,,t-1]*t(G)%*%B
ft[t]=t(F[t,])%*%at[,t]
St[,t]=Rt[,,t]%*%F[t,]
qt[t]=t(F[t,])%*%St[,t]
rt[t]=1/qt[t]
zt[t]=exp(-ft[t])/qt[t]
rtp[t]=rt[t]+yt[t]
ztp[t]=zt[t]+1
gt[t]=log(rtp[t]/ztp[t])+1/(2*rtp[t])
pt[t]=(2*rtp[t]-1)/(2*rtp[t]^2)
mt[,t]=at[,t]+St[,t]*(gt[t]-ft[t])/qt[t]
Ct[,,t]= Rt[,,t]- St[,t]%*%t(St[,t]) *
(1 -pt[t]/qt[t])/qt[t]
pred[t]<-rt[t]/ zt[t]
var.pred<-rt[t]*(zt[t]+1)/(zt[t])^2
icl.pred[t]<-pred[t]-2*sqrt(var.pred)
icu.pred[t]<-pred[t]+2*sqrt(var.pred)
icl.ct1<-sqrt(Ct[1,1,])
list(mt[,t],diag(Ct[,,t]))
}
89
A.2 Codigo usado em WinBUGS
#########################################################
#modelo :
# y_t|lambda_t ~ poisson(lambda_t)
# eta_t=log(lambda_t)=alpha_t + beta_t*CO, F_t=(1,CO)
# alpha_t = alpha_{t-1} + v_t, v_t ~ (0, V_t)
# beta_t = beta_{t-1} + w_t,w_t ~ (0, W_t)
#########################################################
model
alfa[1]~dnorm(b1,tau1)
beta[1]~dnorm(b2,tau2)
lnlambda[1]<- alfa[1]+ beta[1]*x[1]
lambda[1] <- exp(lnlambda[1])
y[1] ~ dpois(lambda[1])
y.rep[1] ~ dpois(lambda[1])
for(i in 2:1443){
alfa[i]~dnorm(alfa[i-1], tau1)
beta[i]~dnorm(beta[i-1], tau2)
lnlambda[i]<- alfa[i] + beta[i]*x[i]
lambda[i] <- exp(lnlambda[i])
y[i] ~ dpois(lambda[i])
y.rep[i]~dpois(lambda[i])}
sigma1<- 1/tau1
sigma2<- 1/tau2
b1~dnorm(0, 0.7)
b2~dnorm(0, 0.7)
tau1 ~dgamma(2,0.0018)
tau2 ~dgamma(2,0.0018)
90
A.3 Codigo usando a Biblioteca INLA
#########################################################
#modelo :
# y_t|lambda_t ~ poisson(lambda_t)
# eta_t=log(lambda_t)=alpha_t + beta_t*CO, F_t=(1,CO)
# alpha_t = alpha_{t-1} + v_t, v_t ~ (0, V_t)
# beta_t = beta_{t-1} + w_t,w_t ~ (0, W_t)
#########################################################
require(INLA)
raizdados = "C:\\Users\\Techy\\Desktop\\Cap_INLAMCMCLB_2regresoras\\
MDLD\\Dados"
setwd(raizdados)
dados1 = matrix(scan("dadosim3.txt",dec=".",skip=0,
what=double(0)),1458,5,byrow=T)
dados1[1:3,]
n.tot=nrow(dados1) # 1458
n=n.tot-14 # 1444
dados=dados1[1:n,]
y=dados[,1]
x1=dados[,3]
b0=dados[,4] #valor verdadeiro do nivel
b1=dados[,5] #valor verdadeiro do coef da regressora
# defining indices for rw1 coefficients
# -------------------------------------
id <- id1 <- 1:n
# formulating the model
# ---------------------
formula <- y ~ f(id, model="rw1",param=c(2,0.0018),initial=5, constr=F) +
91
f(id1, x1, model="rw1",param=c(2,0.0018),initial=5,
constr=F) -1
r = inla(formula, family="poisson", data = data.frame(id,id1,y),
control.predictor=list(compute=TRUE))
#plot(r)
#summary(r)
##################################################################
## graph for y
win.graph()
rang <- range(r[[18]][1:n, 3:5], y)
plot(y, type="l", ylim=rang, xlim=c(1,n),ylab="y",xlab="time")
lines(r[[18]][1:n,1], col=2, lty=3)
lines(r[[18]][1:n,3], col="blue", lty=3)
lines(r[[18]][1:n,5], col="blue", lty=3)
legend("topleft", legend=c("simulated","posterior mean","95% CI"),
col=c("black", "red","blue"),lty=c(1,1,2),bty="n")
## graph for states (b0_t)
win.graph()
rang <- range(r[[13]][[1]][1:n,4:6], b0[1:n])
plot(r[[13]][[1]][1:n,2], type="l",
ylim=rang, col="red", xlim=c(1,n),ylab="b0",xlab="time")
lines(r[[13]][[1]][1:n,4], col="blue", lty=3)
lines(r[[13]][[1]][1:n,6], col="blue", lty=3)
lines(b0[1:n])
legend("topleft", legend=c("simulated","posterior mean","95% CI"),
col=c("black", "red","blue"),lty=c(1,1,2),bty="n")
## graph for states (b1_t)
win.graph()
92
rang <- range(r[[13]][[2]][1:n,4:6], b1[1:n])
plot(r[[13]][[2]][1:n,2], type="l",
ylim=rang, col="red", xlim=c(1,n),ylab="b1",xlab="time")
lines(r[[13]][[2]][1:n,4], col="blue", lty=3)
lines(r[[13]][[2]][1:n,6], col="blue", lty=3)
lines(b1[1:n])
legend("bottomleft", legend=c("simulated","posterior mean","95% CI"),
col=c("black", "red","blue"),lty=c(1,1,2),bty="n")
# hyperparamerter estimation values
# ---------------------------------
r$summary.hyperpar
# posterior densities for precision parameters
# --------------------------------------------
W0=0.001
W1=0.0005
# summary of posterior hyperparameters
# ------------------------------------
r[[21]] #r$summary.hyperpar, s~ao os hiperparametro dinamico (nivel)
# mean sd 0.025quant 0.5quant 0.975quant
#Precision for id 748.0065 208.3632 418.2954 721.6714 1230.515
#Precision for id1 2506.0592 896.1714 1168.1072 2369.5505 4639.932
# posterior density for precision parameter
# -----------------------------------------
# precision for V_t
win.graph()
plot(r[[22]][[1]], type="l", xlab="1/V", main="precision for V_t", ylab="",
xlim=c(0,2000),lwd=2)
abline(v=r[[21]][1,1],col=4, lty=2,lwd=2)#mean
93
points(x=r[[21]][1,3],y=0,cex=1.5,pch=19,col=4) #q25
points(x=r[[21]][1,5],y=0,cex=1.5,pch=19,col=4) #q75
#legend(2000,0.0015,legend=c("simulated","INLA","MCMC"),
col=c(2,4,1),lty=c(1,2,1),bty="n")
abline(v=1/W0,col=2, lwd=3)
# precision for W_t
win.graph()
plot(r[[22]][[2]], type="l", xlab="1/V", main="precision for W_t",
ylab="", xlim=c(0,6000),lwd=2)
abline(v=r[[21]][2,1],col=4, lty=2,lwd=2)#mean
points(x=r[[21]][2,3],y=0,cex=1.5,pch=19,col=4) #q25
points(x=r[[21]][2,5],y=0,cex=1.5,pch=19,col=4) #q75
#legend(2000,0.0015,legend=c("simulated","INLA","MCMC"),
col=c(2,4,1),lty=c(1,2,1),bty="n")
abline(v=1/W1,col=2, lwd=3)
#####################################################################
#Distribuic~ao preditiva para o preditor linear \eta_T|\eta_{-T}
n.pred=1 ##forescat
n1=n-n.pred
yf = c(y[1:n1],rep(NA,n.pred))
id <- id1 <- 1:n
formula1 <- yf ~ f(id, model="rw1",param=c(2,0.0018),initial=5,constr=F) +
f(id1, x1, model="rw1",param=c(2,0.0018),initial=5,constr=F) -1
r1 = inla(formula1, family="poisson", data = data.frame(id,id1,y),
control.predictor=list(compute=TRUE))
#letura das marginais do preditor linear
PredH1<-r1[[17]][[(n1+1)]] #r1$marginals.linear.predictor[[(n1+1)]]
94
win.graph()
plot(PredH1,xlab="",ylab="Density")
lines(inla.smarginal(PredH1))
#Grafico da taxa marginal
win.graph()
plot(inla.tmarginal(exp,PredH1)) # taxa.marg3)
summary(inla.tmarginal(exp,PredH1))
amostra.taxa1<-inla.tmarginal(exp,PredH1)[,1]
amostra.y1<-rpois(length(amostra.taxa1),amostra.taxa1)
#----------------------------------
# Grafics da distribuic~ao preditiva
#----------------------------------
win.graph()
hist(amostra.y1,freq=F,breaks = 17,main="",xlab="")
abline(v=y[(n1+1)],col=2,lwd=2) #valor dos dados reais
legend("topright", legend=c("INLA","valor real"),
col=c("black", "red"),lwd=c(2,2), bty="n")
95
Apendice B
Codigo Usados aos Dados de
Contagem de Obitos
A seguir, tem-se os codigos utilizados pela biblioteca INLA, para ajustar o melhor modelo,
para dados de contagem de obitos de criancas menores de 5 anos. O melhor modelo e dado
por:
#########################################################
#Biblioteca INLA com parametros dinamicos e estaticos #
#########################################################
## regression Poisson dynamic: nivel variando no tempo
##-------------------------------------------------------
## Observational equation:
## y_t \sim Poisson(\lambda_t)
## log(\lambda_t) = \eta_t = \alpha_t + \beta*F_t , t=1,...,n
## System equations:
## \alpha_t = \alpha_{t-1} + w_{1,t} , t=2,...,n
## F_t =cbind(pollag2,templag3,umidlag2,cost,sent)
## \beta = c(b1,b2,b3,b4,b5)
##-------------------------------------------------------
96
B.1 Codigo do modelo Poisson Tradicional
require(INLA)
#Leitura de dados para inla (dados reais)
dados1 = matrix(scan("dadosreais_P2T3U2CS.txt")
n1=nrow(dados1)
n.pred=15 ##forescat
n=n1-n.pred # 1458-15=1443
dados=dados1[1:n,]
y=dados[,1] #obitos
nivel=dados[,2]
pollag2=dados[,3]
templag3=dados[,4]
umidlag2=dados[,5]
cost=dados[,6]
sent=dados[,7]
# defining indices for rw1 coefficients
# -------------------------------------
id <- 1:n
# formulating the model
# ---------------------
formula <- y ~ pollag2+templag3+umidlag2+cost + sent +
f(id, model="rw1", param=c(2,0.007), initial=5, constr=F) -1
# call to fit the model
# ---------------------
r = inla(formula, family="poisson",
data = data.frame(pollag2, templag3, umidlag2, cost, sent, id,y),
control.predictor=list(compute=TRUE))
plot(r)
summary(r)
97
B.2 Codigo do modelo Poisson inflacao-zeros do tipo 0
require(INLA)
# leitura dod dados do modelo completo
dados1 = matrix(scan("dadosreais_P2T3U2CS.txt")
# formulating the model
# ---------------------
id <- 1:n
formula <- y ~ pollag2 +templag3 +umidlag2 +cost +sent +
f(id, model="rw1",param=c(2,0.007), initial=5, constr=F) -1
# call to fit the model
# ---------------------
r = inla(formula, family="zeroinflatedpoisson0",
data = data.frame(pollag2,templag3,umidlag2,cost, sent,id,y),
control.predictor=list(compute=TRUE))
summary(r)
plot(r)
B.3 Codigo do modelo Poisson inflacao-zeros do tipo 1
require(INLA)
# leitura dod dados do modelo completo
dados1 = matrix(scan("dadosreais_P2T3U2CS.txt")
# formulating the model
# ---------------------
id <- 1:n
formula <- y ~ pollag2 +templag3 +umidlag2 +cost +sent +
f(id, model="rw1",param=c(2,0.007), initial=5, constr=F) -1
# call to fit the model
# ---------------------
98
r = inla(formula, family="zeroinflatedpoisson1",
data = data.frame(pollag2,templag3,umidlag2,cost, sent,id,y),
control.predictor=list(compute=TRUE))
summary(r)
plot(r)
99
Apendice C
Codigo Referentes aos Dados de
Material Particulado
A seguir, tem-se os codigos utilizados pela biblioteca INLA, para ajustar os modelos Gama
e Binimial, para dados de material particulado. O melhor modelo e dado por:
#########################################################
#Modelo Gama e Binomial: nivel variando no tempo #
#########################################################
## Observational equation:
## log(phi/\lambda_t) = \eta_t = \alpha_t + \beta*F_t , t=1,...,n
## System equations:
## \alpha_t = \alpha_{t-1} + w_{1,t} , t=2,...,n
## F_t =cbind(rainpad,tmpminpad,wetpad,MON,TUE,WED, THU, FRI,SAT,FERIADO)
## \beta = c(b1,...,b10)
##-------------------------------------------------------
C.1 Codigo do modelo Gama
require(INLA)
# leitura dod dados do modelo completo
100
dados1 = matrix(scan("dados_Riof.CSV")
n1=nrow(dados1) # 730
n=n1-15
dados=dados1[1:n,]
y=PM=dados[,1]
rainpad=dados[,4] #chuva padronizada
wetpad=dados[,8] # umindade padronizada
tmpminpad=dados[,9] #temp minima padronizada
MON=dados[,10]
TUE=dados[,11]
WED=dados[,12]
THU=dados[,13]
FRI=dados[,14]
SAT=dados[,15]
FERIADO=dados[,16]
id <- 1:n
# formulating the model
# ---------------------
formula <- y ~ rainpad+tmpminpad+wetpad+MON+TUE+WED+THU+ FRI+
SAT+FERIADO+ f(id, model="rw1",param=c(2,0.007), initial=5,constr=F) -1
# call to fit the model
# ---------------------
r = inla(formula, family="gamma",
data = data.frame(rainpad,tmpminpad,wetpad,MON, TUE,WED,THU,FRI,
SAT,FERIADO,id,y),
scale =1.1,control.predictor=list(compute=TRUE))
summary(r)
plot(r)
101
C.2 Codigo do modelo Bernoulli
require(INLA)
# leitura dod dados do modelo completo
dados1 = matrix(scan("dados_Riof.CSV"))
n1=nrow(dados1) # 730
n=n1-20 #pred horizonte 20
id <- 1:n
# formulating the model
# ---------------------
formula <- y ~ rainpad+tmpminpad+wetpad+MON+TUE+WED+THU+
FRI + SAT + FERIADO + f(id, model="rw1", hyper=list(list(prior="loggamma",
param = c(2,0.007), initial = 2)))-1
r = inla(formula, family="binomial", Ntrials=rep(1, n),
data = data.frame(rainpad,tmpminpad,wetpad,MON, TUE,WED,THU,FRI,SAT,
FERIADO,id,y),
control.family = list(link = "logit"),
control.predictor=list(compute=TRUE))
summary(r)
plot(r)
102
Referencias Bibliograficas
Abramowitz, M. e Stegun, I, A. (1964) Handbook of Mathematical Functions. New York:
Dover.
Albi, M. O. S., Migon, H. S. e da Silva, C. Q. (2009) Bayesian inference for the beta
distribution. Technical Report-RTPE07-Dep. Estatıstica-UnB.
Alves, M. B. (2006) Funcoes de Transferencia em Modelos Dinamicos Lineares Generalizados.
Tese de Doutorado, Universidade Federal do Rio de Janeiro, Rio de Janeiro.
Alves, M. B., Gamerman, D. e Ferreira, M. A. (2010) Transfer functions in dynamic
generalized linear models. Statistical Modelling, 10(1), 03–40.
Carlin, B. P., Polson, N. G. e Stoffer, D. S. (1992) A Monte Carlo approach to nonnormal and
nonlinear state-space modeling. Journal of the American Statistical Association, 87(418),
493–500.
Carter, C. K. e Kohn, R. (1994) On Gibbs sampling for state space models. Biometrika,
81(3), 541–553.
Dominici F, Daniels M, Z. S. e J., S. (2002) Air pollution and mortality: estimating regional
and national dose-response relationships. Journal of the American Statistical Association,
97, 100–111.
Fahrmeir, L. (1992) Posterior mode estimation by extended Kalman filtering for multivariate
dynamic generalized linear models. Journal of the Americal Statistical Association,
87(418), 501–509.
103
— (1997) Penalized likelihood estiamtion and iterative Kalman smoothing for non-Guassian
dynamic regression models. Computational Statistics and Data Analysis, 24, 295–320.
Fahrmeir, L., Hennevogl, W. e Klemme, K. (1992) Smoothing in dynamic generalized linear
models by Gibbs sampling. In advances in glim and Statistical Modelling, Lecture Notes in
Statistics. Springer, New York, 78, 85–90.
Fernandes, M. V. M. (2006) Modelos para Processos EspaA§os-Temporais Inflacionados de
Zeros. Dissertacao de Mestrado, Universidade Federal do Rio de Janeiro, Rio de Janeiro.
Ferreira, M. A. R. e Gamerman, D. (2000) Dynamic generalized linear models. In Dey. D. K.,
Ghosh, S. K., e Mallick, B. K. Generalized Linear Models, 85, 57–72.
Franco, G. C., Gamerman, D. e Santos, T. R. (2009) Modelos de Espacos de Estados:
Abordagens Classica e Bayesiana. Associacao Brasileira de Estatıstica.
Fruhwirth Schnatter, S. (1994) Data augmentation and dynamic linear models. Journal of
Time Series Analysis, 15, 183–202.
Gamerman, D. (1997) Markov Chain Monte Carlo Stochastic Simulation for Bayesian
Inference. Chapman & Hall.
— (1998) Markov chain Monte Carlo for dynamic generalized linear models. Biometrika, 85,
215–227.
Gamerman, D. e Lopes, H. (2006) Markov Chain Monte Carlo: Stochastic Simulation for
Bayesian Inference. New York: Chapman & Hall / CRC.
Gelfand, A. E. e Smith, A. F. M. (1990) Sampling-based approaches to calculating marginal
densities. Journal of the American Statistical Association, 85, 398–409.
Gelman, A. e Rubin, D. R. (1992) Inference from iterative simulation using multiple sequences
(with discussion). Statistical Science, 7, 457–511.
Geman, S. e Geman, D. (1984) Stochastic relaxation, Gibbs distributions and the Bayesian
restoration of images. IEEE Transactions on Pattern Analysis and Machine Intelligence, 6,
721–741.
104
Geweke, J. (1992) Evaluating the accuracy of sampling-based approaches to the calculation
of the posterior moments (with discussion). Em Bernardo, J.M. , Berger, J. O., Dawid, A.
P. e Smith, A. F. M. Em bayesian statistics 4. Oxford University Press, 4, 169–193.
Harrison, P. e Stevens, C. (1976) Bayesian forecasting (with discussion). J. Roy. Statist. Soc.,
B38, 205–247.
Hastings, H. (1970) Monte Carlo sampling methods using Markov chains and their
applications. Biometrika, 57, 97–109.
Junger, W. L. (2002) Imputacao de dados faltando em series temporais multivariadas via
algoritmo EM. Monografia de graduacao, Universidade Federal do Rio de Janeiro.
Knorr-Held, L. (1999) Conditional prior proposals in dynamic models. Scandinavian Journal
of Statistics, 26, 129–144.
Martins, T. (2010) Aproximacoes Analıticas de Distribuicoes a Posteriori. Dissertacao de
Mestrado, Universidade Federal do Rio de Janeiro, Rio de Janeiro.
Metropolis, N., Rosenbulth, A., Rosenbulth, M., Teller, A. e Teller, E. (1953) Equation of
state calculations by fast computing machine. Journal of Chemical Physics, 21, 1087–1091.
Migon, H. S. e Gamerman, D. (1999) Statistical Inference: an Integrated Approach. Arnold.
Migon, H. S., Gamerman, D., Lopes, H. F. e Ferreira, M. A. R. (2005) Bayesian dynamic
models. Em Dey, D. e Rao, C.R. Handbooks of Statistics, 25, 553–588.
Migon, H. S. e Harrison, P. (1985) An Application of Non-Linear Bayesian Forecasting to
Television Advertising, in Bayesian Statistic 2. Amsterdam: North Holland.
Migon, H. S., Schmidt, A. M., Ravines, R. e Pereira, J. (2013) An efficient sampling scheme
for generalized dynamic models. Computational Statistics, 28, 2267–2293.
Nelder, J. A. e Wedderburn, R. W. M. (1972) Generalized linear models. Journal of the Royal
Statistical Society A, 135, 370–384.
105
Petris, G., Campagnoli, P. e Petrone, S. (2009) Dynamic Linear Models with R. New York:
Springer.
R Development Core Team (2013) R: A Language and Environment for Statistical Computing.
R Foundation for Statistical Computing, Vienna, Austria. URLhttp://www.R-project.
org/.
Resende, C. M. C. (2011) Inferencia Bayesiana Aproximada em Modelos de Espaco de Estados.
Dissertacao de Mestrado, Universidade Federal do Rio de Janeiro, Rio de Janeiro.
Rue, H. e Martino, S. (2007) Approximate bayesian inference for hierarchical gaussian markov
random field models. Journal of statistical planning and inference, 137(10), 3177–3192.
Rue, H., Martino, S. e Chopin, N. (2009) Approximate Bayesian inference for latent Gaussian
models by using integrated nested Laplace approximations. Journal of the Royal Statistical
Society B, 71, 319–392.
Ruiz-Cardenas, R., Krainski, E. T. e Rue, H. (2010) Fitting dynamic models using
integrated nested laplace approximations - inla. Preprint Statistics No 12/2010,
Department of Mathematical Sciences, Norwegian University of Science and Technology,
Trondheim, Norway. URLhttp://www.math.ntnu.no/inla/r-inla.org/papers/
Ruiz_Krainski_Rue_TR_23112010.pdf.
Shephard, N. e Pitt, M. K. (1997) Likelihood analysis of non-Gaussian measurement time
series. Biometrika, 84, 653–667.
Singh, A. C. e Roberts, G. R. (1992) State space modelling of cross-classified time series of
counts. International Statistics Review, 60, 321–336.
Smith, A. F. M., Skene, A. M., Shaw, J. E. H. e Naylor, J. C. (1987) Progress with numerical
and graphical methods for practical bayesian statistics. Journal of the Royal Statistical
Society. Series D (The Statistician), 36, 75–82.
Tierney, L. (1994) Markov chains for exploring posterior distribuitions. Annals of Statistics,
22, 1701–1728.
106
Tierney, L. e Kadane, J. B. (1986) Accurate approximations for posterior moments and
marginal densities. Journal of the American Statistics Association, 81, 82–86.
Vedal S, Brauer M, W. R. e J., P. (2003) Air pollution and daily mortality in a city with low
levels of pollution. Environmental Health Perpespectives, 111, 45–51.
West, M. e Harrison, J. (1997) Bayesian Forecasting and Dynamic Models. New York:
Springer-Verlag, 2a edn.
West, M., Harrison, J. e Migon, H. S. (1985) Dynamic generalized linear models and bayesian
forecasting. Journal of the Americal Statistical Association, 80, 389, 73–83.
107
