Seleção de Preditores em Modelos de Regressão

Prévia do material em texto

Seleção de Preditores em Modelos de
Regressão
Carolina Marques Bastos
Universidade Federal do Rio de Janeiro
Instituto de Matemática
Departamento de Métodos Estat́ısticos
2011
Seleção de Preditores em Modelos de Regressão
Carolina Marques Bastos
Dissertação submetida ao Corpo Docente do Instituto de Matemática - Departamento
de Métodos Estat́ısticos da Universidade Federal do Rio de Janeiro - UFRJ, como parte
dos requisitos necessários à obtenção do grau de Mestre em Estat́ıstica.
Aprovada por:
Profo. Helio S. Migon
DME - UFRJ - Orientador
Profo. Thais C. O. Fonseca
DME - UFRJ
Profo. Marco A. Rodŕıguez
Université du Québec à Trois-Rivières
Rio de Janeiro
Agosto de 2011
ii
Agradecimentos
Agradeço a Deus por sempre colocar na minha vida ótimas oportunidades, por ter
me capacitado e guiado para a conclusão de mais uma etapa da minha vida.
A minha famı́lia, agradeço por sempre estar ao meu lado, dando apoio em todos os
passos da minha vida e comemorando cada vitória alcançada. Agradeço aos meus pais,
Katia e Heloy, por me apoiarem em todos os momentos, por todo o esforço para que
eu tivesse as melhores condições de estudo e por darem muito valor a cada conquista.
Obrigada por tudo!
Ao meu noivo Luiz, que acompanha de perto todas as vitórias ao longo do tempo
que estamos juntos. Gostaria de agradecer a sua compreensão e paciência em todos os
meus surtos diante desta dissertação. Sem o seu apoio, carinho e incentivo, eu não teria
chegado até aqui.
Agradeço às companheiras da FGV, Lúısa e Samanta, pela força e apoio de sempre.
Ao Marcelo Neri, por me incentivar desde o ińıcio e mostrar a importância do mestrado
na minha formação. Agradeço pela licença que me foi concedida durante esse peŕıodo,
por todos os conselhos e conversas.
As minhas amigas de turma: Carol, Camila, Renata e Priguete, que ajudaram muito
com estudos, trabalhos e momentos de lazer. Compartilhamos ótimos momentos. Obri-
gada por sempre estarem presentes ajudando! A amiga Renata um agradecimento es-
pecial por ter me apresentado minha best Luana (hehe)! Agradeço por compreender o
que ela chama de desapego, por segurar a vontade de fofocar durante horas enquanto eu
terminava a minha dissertação e, ela já tinha terminado a dela! Por me proporcionar
muitas risadas, momentos de reflexão, explicações sobre estat́ıstica bayesiana, mcmc e R.
Inexplicável o quanto você me incentivou... Obrigada por toda sua ”bestice”! Aos demais
amigos do LPGE, principalmente Joãozinho, Larisson, Kelly Dance, Sheila, Patylene e
iii
Mari que ajudaram com contas, disciplinas, provas e etc. Também proporcionaram mo-
mentos de risadas e muita descontração, tornando a vida mais alegre em dias de Fundão.
Ao Vini, por ter me ensinado a rodar o WinBugs no R, deixando a parte da programação
da dissertação muito mais prática!!!
Agradeço ao Migon pela orientação, incentivo e paciência. Obrigada por todo conhe-
cimento passado durante o peŕıodo em que estivemos envolvidos na dissertação.
Agradeço ao Marco Rodŕıguez e a Thais Fonseca, por aceitarem fazer parte da minha
banca. A Thais agradeço também por sua disponibilidade e boa vontade ao me passar
seus conhecimentos e sugestões referentes a elaboração desta dissertação.
Agradeço a todos aqueles que de alguma forma torceram por mim e contribúıram
para que eu conclúısse esse curso de mestrado.
Finalmente, agradeço ao CNPQ pelo financiamento dos meus estudos durante este
peŕıodo.
iv
Resumo
O estudo de técnicas que selecionam os preditores de um modelo estimado de forma
criteriosa, é parte fundamental do processo de construção do modelo estat́ıstico. Nesta
dissertação, a proposta é fazer a comparação de métodos de estimação de modelos que
utilizam técnicas para a seleção de preditores. A comparação será feita por meio de
critérios de seleção de modelos já conhecidos na literatura.
A primeira técnica a ser utilizada para a estimação do modelo e seleção de predi-
tores, se baseia na imersão do modelo de regressão em uma estrutura hierárquica de
mistura de normais, onde uma variável latente irá sinalizar quais preditores devem ser
inclúıdos no modelo ou não. Dessa forma, esta técnica não só estima o modelo, como
também tem a capacidade de selecionar os preditores mais relevantes para o mesmo. A
segunda técnica para a estimação de modelos consiste no uso de distribuições de con-
tração para o vetor paramétrico. As distribuições de contração são obtidas via mistura
de funções de distribuições cont́ınuas. Estamos interessados em duas formas particula-
res de obtenção de funções de distribuição de contração: a primeira delas consiste na
mistura do parâmetro de escala da distribuição normal com uma função de distribuição
que seja exponencialmente distribúıda. Também estamos interessados em outra forma de
obtenção de distribuições de contração, em que é feita a mistura do parâmetro de escala
da distribuição normal com uma função de distribuição Cauchy, truncada nos valores
reais positivos.
Verificaremos as vantagens e desvantagens associadas a estas propostas para a es-
timação de modelos, que também tem o objetivo de fazer seleção ou contração dos pre-
ditores. Todo o procedimento de inferência será feito sob o enfoque bayesiano, isto é:
atribúıremos uma distribuição a priori para os parâmetros de interesse do modelo, a fim
de obtermos a distribuição a posteriori que, em nosso caso, não é conhecida. Métodos
v
de Monte Carlo via Cadeias de Markov (MCMC, sigla em inglês) serão utilizados para
obter amostras dessa distribuição.
As técnicas para a estimação do modelo serão aplicadas a um conjunto de dados
gerados de maneira artificial. Para esse conjunto de dados, a quantidade de variáveis
preditoras, a correlação entre elas e o tamanho da amostra, serão variados. Dessa ma-
neira, iremos avaliar qual técnica de estimação de modelos foi a mais eficiente. Um ńıvel
de esparsidade será atribúıdo ao vetor paramétrico, fazendo com que muitos de seus
elementos sejam nulos. Exerćıcios de simulação nos permitem avaliar qual dos métodos
capta melhor a estrutura de esparsidade associada ao vetor paramétrico e calibrar de ma-
neira adequada a implementação das propostas para estimação de modelos. Finalmente,
as técnicas de estimação propostas e avaliadas serão aplicadas a exemplos que utilizam
dados reais.
Palavras-Chaves: Estimação de modelos, seleção de preditores, misturas cont́ınuas,
distribuição de contração.
vi
Abstract
In statistics, a crucial problem in building a multiple regression model is the selection
of predictors to include. In this work, we will compare methods for model estimation
that use techniques that select the predictors. The comparison will be made using known
criteria for model selection.
The first technique to be used for model estimation and selection of predictors, entails
embedding the regression setup in a hierarchical normal mixture model where latent
variables are used to identify which predictors should be included in the model. This
technique can estimate and select the most relevants predictors for this. The second
technique for model estimation, is based on shrinkage priors obtained by normal scale
mixtures. We are interested in two particular ways of obtaining shrinkage distributions:
the first one is obtained by normal scale mixtures with exponential distributions. Also
we are interested in another way of obtaining shrinkage distributions, by normal scale
mixtures with a standard half-Cauchy distribution on the positive reals.
We examine the proposal’s advantages and disadvantages. These proposals for model
estimation also objectively select or shrink predictors. All the inference procedure follows
the Bayesian approach, that is, we attribute a priordistribution for the parameters of
interest of each model to obtain the posterior distribution which, in our case, is not
known. Markov chain Monte Carlo methods (MCMC) are used to obtain samples of this
distribution
The proposed techniques for model estimation will be applied to data sets having
different numbers of predictors, correlation among predictors and sample size. We analyze
which technique for model estimation is more efficient. The parametric vector has a
sparsity level, such that many of its elements are null. A simulation exercise allows
us to evaluate which method better captures the sparsity level and standardizes the
vii
implementation of proposals for model estimation. Finally, the proposed estimation
techniques will be applied in a example based on a real data set.
Keywords: Model estimation, predictor selection, continous mixtures, shrinkage dis-
tributions.
viii
Sumário
1 Introdução 1
2 Modelos Bayesianos e Métodos de Estimação 6
2.1 Inferência Bayesiana e Métodos de Estimação . . . . . . . . . . . . . . . 6
2.1.1 Estimação Pontual . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.1.2 Estimação por Intervalo . . . . . . . . . . . . . . . . . . . . . . . 8
2.2 Inferência Via Simulação Estocástica . . . . . . . . . . . . . . . . . . . . 8
2.2.1 Inferência Via MCMC . . . . . . . . . . . . . . . . . . . . . . . . 8
2.2.2 WinBugs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
3 Seleção de Variáveis 14
3.1 Seleção de Variáveis via Busca Estocástica . . . . . . . . . . . . . . . . . 16
3.2 Operador de Seleção e Contração com Penalidade em Valor Absoluto . . 20
3.2.1 Formulação Hierárquica do Modelo Lasso Bayesiano . . . . . . . . 22
3.2.2 Função de Contração . . . . . . . . . . . . . . . . . . . . . . . . . 24
3.2.3 Função de Influência . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.3 Mistura de normais na Escala Usando Distribuições de Cauchy . . . . . . 31
3.3.1 Formulação Hierárquica do Modelo . . . . . . . . . . . . . . . . . 31
3.3.2 Função de Contração . . . . . . . . . . . . . . . . . . . . . . . . . 33
3.3.3 Função de Influência . . . . . . . . . . . . . . . . . . . . . . . . . 36
4 Critérios de Seleção de Modelos 39
4.1 Critérios Baseados na Função de Verossimilhança Marginal . . . . . . . . 40
ix
4.1.1 Fator de Bayes . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
4.1.2 Escores Logaŕıtmicos . . . . . . . . . . . . . . . . . . . . . . . . . 42
4.2 Critério de Informação Baseado no Desvio - DIC . . . . . . . . . . . . . . 43
4.3 Critério Baseado na Função de Perda Canônica . . . . . . . . . . . . . . 44
4.4 Critérios Baseados em Função de Perda Quadrática e Valor Absoluto . . 46
5 Métodos de Estimação de Modelos Aplicados em Modelos de Regressão
Linear 48
5.1 Aplicação a Dados Simulados . . . . . . . . . . . . . . . . . . . . . . . . 48
5.2 Estimação de Modelos - Aplicação a Dados Reais . . . . . . . . . . . . . 63
6 Conclusão 73
Cálculo da Distribuição de Contração 77
Cálculo dos Estimadores da Média Harmônica 79
Análises Referentes a Aplicação dos Dados sobre Câncer de Próstata 81
x
Caṕıtulo 1
Introdução
Um dos grandes objetivos na estat́ıstica, é o desenvolvimento e a estimação de um
modelo estocástico, descrevendo as variáveis de interesse para certo estudo. Modelos
estocásticos podem ser usados em diversas áreas do conhecimento. Com efeito, em eco-
nomia podemos fazer uso de um modelo estátistico para estimar o lucro de uma empresa,
sujeito a determinadas caracateŕısticas. Na demografia, um modelo pode ser elaborado
para estimar qual será a população do páıs daqui a 10 anos. Já na área da saúde, um
modelo estat́ıstico pode ser capaz de associar fatores relativos ao estilo de vida de um
paciente (prática de atividades f́ısicas, boa alimentação, ingestão de cálcio) com a chance
dele adquirir uma doença, como a osteoporose, daqui a alguns anos. Assim, com poucos
exemplos, já foi posśıvel perceber o quão importante e abrangente um modelo estat́ıstico
pode ser.
A inferência em modelos estat́ısticos pode não ser uma tarefa simples, mas é um
conceito amplamente difundido. Além disso, a grande maioria dos softwares estat́ısticos
possuem pacotes que são capazes de estimar modelos sem maiores dificuldades. Os pro-
blemas podem começar a surgir quando a quantidade de vaŕıaveis candidatas a predizer
uma determinada caracteŕıstica for muito elevada, uma vez que modelos com muitas
variáveis explicativas tendem a ser complicados para interpretar.
É cada vez mais frequentente na estat́ıstica moderna, estudos aplicados em que surge a
necessidade de lidar com bases de dados muito grandes. Para a manipulação de tais bases,
sem lançar mão de nenhuma informação, é necessário o conhecimento de ferramentas que
1
tenham a capacidade de lidar com problemas de dimensão muito elevadas.
Um estudo apresentado em Chandulaka et al. (2010), que será denominado como o
problema de marketing, lida com a estimação de um modelo cuja variável de interesse é
a intenção dos consumidores na compra de um determinado produto. Nesse problema a
intenção de compra de um produto, pode ser influenciada por variáveis relacionadas às
atitudes dos consumidores, suas crenças e a publicidade do produto. É importante notar
que, estamos lidando com um número elevado de variáveis explicativas, uma vez que
temos distintas quantidades de variáveis associadas às atitudes dos consumidores, tais
como: recomendação do produto a amigos, ”test-drive”do produto, estudo do produto
antes da compra, entre outras. Algumas variáveis associadas às crenças dos consumidores
são: durabilidade, segurança, qualidade do produto e outras. Por fim, variáveis associ-
adas a publicidade do produto são: publicidade na internet, publicidade na televisão e
outras diversas formas de publicidade de um produto ou marca.
Visando o manuseio de tantas variáveis explicativas e a estimação de um modelo que
explique o problema de maketing, Chandulaka et al. (2010) desenvolvem um modelo de
efeitos hierárquicos nas variáveis. Eles caracterizam os denominados modelo de efeitos
hierárquicos nas variáveis, isto é, a relação entre certas variáveis aleatórias são media-
das por outras variáveis. Assim, modelos de efeitos hierárquicos nas variáveis podem
ser analisados fatorando-se a distribuição conjunta das variáveis como um produto de
distribuições condicionais e marginais.
Para a melhor compreensão do conceito de efeitos hierárquicos nas variáveis e do
modelo proposto, suponha o grupo das três variáveis aleatórias associadas ao problema de
marketing: (x, b, z), onde x corresponde a um grupo de variáveis aleatórias relacionadas a
atitudes dos consumidores mediante um determinado produto, b é um grupo de variáveis
relacionadas às crenças dos consumidores a respeito do produto em questão e z são
variáveis relacionadas a publicidade do produto. O objetivo é analisar a influência desse
grupo de variáveis (x, b, z) na quantidade aleatória y, que é a intenção do consumidor na
compra do produto, e que deve ser entendida como a variável de interesse. Obviamente a
distribuição conjunta dessas variáveis (y, x, b, z) admite distintas fatorações. No contexto
de efeitos hierárquicos nas variáveis, cada uma dessas distintas fatorações descreverá um
2
modelo a ser analisado.
Uma posśıvel fatoração da distribuição conjunta das variáveis em questão, indica que
as atitudes e as crenças dos consumidores, x e b respectivamente, influenciam diretamente
a intenção de compra do produto. A Figura a seguir ilustra os efeitos hierárquico nas
variáveis para esta particular fatoração:
Para este caso particular, o modelo a ser estudado é dado por [y|x, b], [x] e [b], onde os
colchetes representamdistribuições condicionais e marginais. Assim, a intenção de com-
pra do produto pode ser explicada pelas atitudes e crenças dos consumidores a respeito
do mesmo.
Uma fatoração distinta, indica que z e x, as variáveis relacionadas a publicidade
do produto e atitudes dos consumidores, influenciam diretamente a intenção de compra
do produto. Adicionalmente, temos as variáveis associadas a publicidade do produto
afetando diretamente as variáveis relacionadas as atitudes dos consumidores, x. A Figura
a seguir representa os efeitos hierárquicos nas variáveis para esta distinta fatoração:
A ilustração acima exibe a representação do efeito indireto da publicidade na intenção
de compra do produto, uma vez que a publicidade está afetando diretamente as atitudes
dos consumidores. Posteriormente, as atitudes dos consumidores x irão ter um efeito
direto na intenção de compra do produto. Portanto, nesse contexto, variáveis associadas
3
a publicidade afetam a intenção de compra de maneira direta e indireta, o que pode
ser melhor verificado na Figura acima. Para esta particular fatoração, o modelo a ser
analisado é dado por [y|x, z], [x|z] e [z], onde os colchetes representam distribuições
condicionais e marginais.
Exitem outras fatorações posśıveis, porém, a descrição de cada uma destas fatorações
não é relevante para o contexto. O importante é notar que que cada uma dessas distintas
fatorações, descreverá diferentes modelos a serem analisados. Além disso, podemos ter
muitas variáveis explicativas associadas a variável de interesse, o que dificultará a es-
timação e interpretação dos modelos. Ainda podemos lidar com o caso em que algumas
variáveis explicativas tem pouca influência na intenção de compra do produto, ou uma
influência não significativa. Nesse contexto, o uso de técnicas que façam o procedimento
de seleção das variáveis relevantes pode auxiliar muito, fazendo com que tenhamos um
modelo mais parcimonioso e de fácil interpretação. Foi diante deste contexto, que surgiu
pela primeira vez, a necessidade da implementação de um método que selecione preditores
relevantes para o modelo de maneira eficiente.
Nesta dissertação, iremos revisar e discutir a aplicação de algumas técnicas de seleção
de variáveis em modelos lineares de regressão. Todo o procedimento de inferência será
feito sob o enfoque bayesiano, isto é, atribuiremos uma distribuição a priori para os
parâmetros de interesse de cada modelo a fim de obter a distribuição a posteriori, e a
partir dela, realizar todo o processo de estimação. Obteremos amostras desta distribuição
a posteriori por meio de métodos de simulação estocástica, particularmente utilizaremos
os métodos de Monte Carlo via cadeias de Markov (MCMC na sigla em inglês).
Procuraremos verificar as vantagens e desvantagens entre cada uma das técnicas pro-
postas para a seleção de preditores. Avaliaremos tais vantagens e desvantagens sob o
contexto teórico e aplicado associado a cada uma das técnicas, e também utilizaremos
critérios de seleção de modelos conhecidos na literatura, como o fator de bayes e o DIC,
visando a comparação dos modelos estimados por cada uma das diferentes técnicas.
A relevância deste tema pode ser notada quando analisamos o contexto da modela-
gem estat́ıstica moderna, onde é cada vez mais frequente nos depararmos com situações
em que uma grande quantidade de variáveis regressoras estão associadas a uma determi-
4
nada variável de interesse. Este problema acaba por trazer dificuldades na estimação do
modelo. Por exemplo, preditores relacionados de forma exata ou aproximada geram di-
ficuldades de estimação. Também podemos citar o problema de obtenção de estimativas
imprecisas ou até mesmo não significativas para o modelo. Técnicas de seleção de predi-
tores são necessárias, pois dão a possibilidade de estimar um modelo mais parcimonioso,
com menos variáveis preditoras, facilitando assim, o processo de estimação do modelo e
a sua interpretação.
Esta dissertação está organizada em 6 Caṕıtulos e 3 Apêndices. No Caṕıtulo 2,
faremos uma breve revisão sobre estimação e sobre procedimentos de inferência sob o
enfoque bayesiano. Discutiremos também métodos de simulação estocástica, particu-
larmente métodos MCMC. Também apresentaremos alguns pontos relevantes sobre o
pacote estat́ıstico WinBUGS, utilizado para a estimação dos modelos presentes nesta
dissertação. No Caṕıtulo 3, revisaremos técnicas de interesse para a estimação de mode-
los com seleção de preditores. Aqui serão feitas comparações teóricas entre as diferentes
técnicas propostas, procurando avaliar as vantagens e desvantagens associadas a cada
um dos métodos. No Caṕıtulo 4, apresentaremos as técnicas de comparação de modelos,
que serão utilizadas para a avaliar os métodos de estimação de modelos em um contexto
aplicado. Tal contexto aplicado, será apresentado no Caṕıtulo 5, onde faremos um estudo
simulado para avaliar em quais aspectos as técnicas de estimação de modelo com seleção
de preditores possuem um melhor desempenho. Finalmente, um exemplo com dados será
trabalhado na Seção 5.2, onde as técnicas de estimação de modelos propostas serão apli-
cadas. Por fim, no Caṕıtulo 6, apresentaremos as conclusões e posśıveis extensões desta
dissertação.
5
Caṕıtulo 2
Modelos Bayesianos e Métodos de
Estimação
2.1 Inferência Bayesiana e Métodos de Estimação
Este caṕıtulo tem por objetivo revisar os principais conceitos do procedimento de
inferência sob o enfoque bayesiano. Considere y, uma variável aleatória ou vetor aleatório
com função de probabilidade ou função de densidade de probabilidade p(y|θ) em que θ
é um parâmetro ou vetor paramétrico que caracteriza a distribuição de probabilidade de
y. O valor de θ é desconhecido e queremos estimá-lo. Sob o ponto de vista da inferência
bayesiana, podemos incorporar nossa própria incerteza na estimação de θ, assumindo
uma distribuição de probabilidade para este parâmetro, p(θ), a distribuição a priori.
Esta distribuição é atribúıda antes da observação dos dados e mede a nossa incerteza a
priori a respeito de θ. Uma vez que os dados são observados, os quais denotaremos por y,
podemos encontrar a distribuição a posteriori de θ, p(θ|y), obtida a partir da combinação
da função de verossimilhança p(y|θ) com a distribuição a priori de θ, p(θ), via teorema
de Bayes, da forma:
p(θ|y) = p(y|θ)p(θ)
p(y)
. (2.1)
A quantidade p(y) =
∫
Θ
p(y, θ)dθ =
∫
Θ
p(y|θ)p(θ)dθ, em que Θ é o espaço paramétrico
6
de θ. Note que p(y) não depende de θ, logo o denominador da equação acima pode ser
considerado constante com relação a θ. Portanto, podemos rescrever a equação 2.1 como:
p(θ|y) ∝ p(y|θ)p(θ) (2.2)
O procedimento de inferência bayesiano é baseado fundamentalmente na distribuição
a posteriori de θ. Esta distribuição contém toda informação probabiĺıstica a respeito do
parâmetro de interesse. No entanto, em algumas situações torna-se necessário resumir
a informação contida na distribuição a posteriori. O caso mais simples é a estimação
pontual, descrita na próxima subseção:
2.1.1 Estimação Pontual
Na estimação pontual, nosso objetivo é a minimização de uma função perda L(δ(Y ), θ)
para algum estimador δ(Y ) de θ. Observe que o valor de θ é estimado a partir de
elementos da amostra. Para cada valor de θ e cada posśıvel estimativa d pertencente ao
espaço paramétrico Θ, associamos uma função de perda L(d, θ). Neste caso, podemos
calcular a perda esperada a posteriori ou risco a posteriori, da forma:
r(d, θ) = E(L(d, θ)|y) =
∫
Θ
L(d, θ)p(θ|y)dθ (2.3)
A regra de Bayes consiste em escolher o valor de d ótimo, ou seja, o valor de d que
minimiza a perda esperada E(L(d, θ)|y). Os estimadoresd(Y), obtidos minimizando esta
perda esperada, serão chamados estimadores de Bayes. As funções perda mais utilizadas
na literatura e seus respectivos estimadores são:
• Perda Quadrática: L(δ, θ) = (δ− θ)2. Neste caso, o estimador resultante é a média
a posteriori de θ, isto é, δ = θ̂ = Eθ|y(θ);
• Perda Absoluta: L(δ, θ) = |δ − θ|. O estimador associado a perda quadrática
absoluta é a mediana a posteriori de θ, δ = med(θ);
• Perda 0-1: L(δ, θ) = lim�→0 I|θ−δ|([�,∞)), onde Ix(A) = 1 se x ∈ A e 0 caso
contrário. O estimador de θ é a moda da distribuição a posteriori de θ.
7
2.1.2 Estimação por Intervalo
Resumir a informação contida na distribuição a posteriori através de um único valor,
θ̂, resulta numa sumarização extrema da informação dispońıvel. É interessante obter pelo
menos uma medida sobre quão precisa é a estimativa de θ̂. Uma maneira de fazer isso é
fornercer uma região de valores θ ∈ Θ, que têm associados a eles uma grande massa de
probabilidade a posteriori. Idealmente, gostaŕıamos de descrever uma região de valores
de θ que é tão pequena quanto posśıvel, mas contém o máximo de probabilidade a
posteriori. Assim define-se o intervalo de credibilidade a posteriori de θ, uma quantidade
desconhecida definida em Θ, como sendo:
uma região C ∈ Θ é uma região de 100(1 − α)% de credibilidade para θ se P (θ ∈
C) ≥ 1 − α. Neste caso, 1 − α é chamado ńıvel de credibilidade. No caso escalar, C é
usualmente dado por um intervalo, por exemplo, [c1, c2].
2.2 Inferência Via Simulação Estocástica
No contexto da inferência estat́ıstica, a simulação estocástica tem o objetivo de esti-
mar caracteŕısticas probabiĺısticas de modelos ou distribuições de interesse, as quais não
poderiam ser obtidas analiticamente.
Métodos de simulação estocástica são comumente utilizados ao fazer inferência sob
a abordagem bayesiana. Eles são uma alternativa razoável para, por exemplo, simular
pontos de forma indireta da distribuição a posteriori, quando esta não possui forma
anaĺıtica fechada ou quando a avaliação por métodos numéricos é inviável, devido a
grandes dimensões paramétricas.
Em geral utilizam-se métodos de Monte Carlo via Cadeias de Markov (MCMC) para
realizar-se o processo inferencial.
2.2.1 Inferência Via MCMC
A inferência sob o paradigma bayesiano parte do pressuposto que a incerteza sobre
uma quantidade desconhecida, digamos θ, pode ser representada por modelos proba-
8
biĺısticos.
Por vezes, o denominador da equação 2.2 não possui forma anaĺıtica fechada e a
avaliação por métodos numéricos quando a dimensão é grande é inviável. Dáı surge
a necessidade de métodos de simulação estocástica, tais como os de Monte Carlo via
Cadeias de Markov (MCMC).
Se uma cadeia de Markov homogênea é irredut́ıvel, recorrente positiva e aperiódica,
então possui distribuição limite, a qual depende apenas da matriz de transição da ca-
deia. Além disso, uma vez que a cadeia atinja a distribuição limite, todos os estados
subsequentes seguirão tal distribuição.
Os métodos MCMC consistem na construção de uma cadeia de Markov que, por meio
de escolhas adequadas de núcleos de transição, tenha como distribuição estacionária a
distribuição de interesse. No contexto de estimação bayesiana, a distribuição a posteriori.
Uma vez que a convergência da cadeia tenha sido atingida, as amostras estarão sendo
geradas da distribuição estacionária. Para aproximar a distribuição a posteriori, utilizam-
se amostras suficientemente grandes dessa distribuição.
Os algoritmos MCMC mais utilizados no contexto de inferência bayesiana são o amos-
trador de Gibbs e o algoritmo de Metropolis-Hastings, que serão descritos a seguir.
(i) Amostrador de Gibbs
O amostrador de Gibbs é um método de MCMC em que o núcleo de transição é
formado pelas distribuições condicionais completas do vetor paramétrico. Assuma
que a distribuição de interesse é π(θ) em que θ = (θ1, . . . , θd)
′. Considere também
que as distribuições condicionais completas πi(θi) = π(θi|θ−i), i = 1, . . . , d são
conhecidas e dispońıveis para a amostragem.
Quando há necessidade de amostrar de π, mas a sua geração direta é complicada,
custosa, ou simplesmente inviável, o amostrador de Gibbs permite um processo de
geração alternativo baseada em gerações sucessivas das distribuições condicionais
completas.
Para construção de uma cadeia de Markov cujas transições sejam definidas pelas
9
condicionais completas, o amostrador de Gibbs procede da seguinte forma:
• Inicialize o contador da cadeia em j = 1 e assuma valores iniciais θ(0) =
(θ
(0)
1 , . . . , θ
(0)
d )
′
• Obtenha o novo valor θ(j) = (θ(j)1 , . . . , θ
(j)
d )
′ de θ(j−1) a partir de gerações
sucessivas:
θ
(j)
1 ∼ π(θ1|θ
(j−1)
2 , . . . , θ
(j−1)
d ) (2.4)
θ
(j)
2 ∼ π(θ2|θ
(j−1)
1 , θ
(j−1)
3 , . . . , θ
(j−1)
d )
...
θ
(j)
d ∼ π(θd|θ
(j−1)
1 , . . . , θ
(j−1)
d−1 )
• Faça j = j + 1 e volte ao passo anterior até obter convergência.
À medida que o número de iterações cresce, a cadeia de Markov simulada aproxima-
se de sua distribuição de equiĺıbrio. Sendo assim, θ(i) = (θ
(i)
1 , . . . , θ
(i)
d )
′ pode ser
considerado um ponto amostrado de π(θ).
(ii) Metropolis-Hastings
Assuma que a distribuição de interesse é π(θ) em que θ = (θ1, . . . , θd)
′. O algoritmo
Metropolis-Hastings é útil para a geração de valores de parâmetros cujas distri-
buições condicionais completas não tenham forma anaĺıtica fechada e, portanto,
não estejam dispońıveis para amostragem, diferentemente do caso do amostrador
de Gibbs.
Neste caso, gera-se valores do parâmetro a partir de uma distribuição proposta
arbitrária e este é aceito ou não com uma certa probabilidade de aceitação, que
depende da qualidade do movimento proposto, avaliado com base na distribuição
proposta e da distribuição de interesse π(θ).
O algoritmo de Metropolis-Hastings procede da seguinte forma:
• Inicialize o contador da cadeia em j = 1 e assuma valores iniciais θ(0) =
(θ
(0)
1 , . . . , θ
(0)
d )
′
10
• Obtenha um valor proposto θ∗ da distribuição proposta q(θ∗|θ(j−1))
• Aceite o valor proposto com probabilidade α(θ∗|θ(j−1)) = min
[
1, π(θ
∗)q(θ(j−1)|θ∗)
π(θ(j−1))q(θ∗|θ(j−1))
]
,
ou seja, θ(j) = θ∗. Caso o valor proposto não seja aceito, faça θ(j) = θ(j−1).
• Faça j = j + 1 e volte ao passo segundo passo até obter convergência.
A escolha da distribuição proposta é uma questão importante ao se utilizar métodos
MCMC com base no algoritmo de Metropolis-Hastings. Sob o ponto de vista
prático, tal escolha é crucial para a sua convergência para a distribuição a pos-
teriori.
Uma das propostas mais comuns são chamadas de cadeias simétricas, quando a dis-
tribuição proposta é simétrica em torno da iteração anterior, isto é, q(θ(j)|θ(j−1)) =
q(θ(j−1)|θ(j)). Neste caso, a probabilidade de aceitação se reduz à razão da distri-
buição de interesse, ou seja, α(θ∗|θ(j−1)) = min
[
1, π(θ
∗)
π(θ(j−1))
]
.
É importante ressaltar que a eficiência do método está diretamente ligada à escala
da distribuição proposta. Caso a variância da distribuição proposta seja muito pe-
quena, a cadeia de Markov irá convergir lentamente, uma vez que seus incrementos
serão pequenos. Se a variância for grande, a taxa de rejeição dos valores propostos
será alta e a cadeia tenderá a não se mover. Muitos autores sugerem que a taxa de
aceitação do algoritmo deve estar entre 20% e 50%, ver Gamerman e Lopes (2006).
Estando decidido o método a ser utilizado, e obtida uma simulação da cadeia, deve-se
verificar se a convergência foi obtida, para assim poder formar a amostra da distribuição
a posteriori das quantidades desconhecidas do modelo.
Existem várias formas de se realizar uma análise a respeito da convergência da cadeia.
Uma das abordagensmais informais é a inspeção gráfica, onde analisa-se a trajetória de
uma ou mais cadeias, com valores iniciais distintos e considera-se que a convergência é
alcançada quando todas as cadeias monitoradas permanecem em torno de um mesmo
ponto. Outros critérios, mais formais, também podem ser utilizados, como os métodos
propostos por Gelman (1992) e Geweke (1992). Neste estudo será utilizado este último
11
critério mencionado. Geweke (1992) sugere um procedimento para teste de convergência
a partir da avaliação de médias ergódicas de uma única cadeia gerada, com base na
idéia de que, após convergência, diferentes intervalos da cadeia gerada devam apresentar
comportamentos semelhantes. Seja uma cadeia gerada com um número de iterações n
suficientemente grande. A idéia é testar a igualdade das médias x̄1 e x̄2, calculadas,
respectivamente, a partir da fração 0.1n inicial e 0.5n final da amostra. Considerando os
respectivos estimadores das variâncias assintóticas de x̄1 e x̄2, dados por V (x̄1) e V (x̄2),
tem-se que, quando n→∞,
Gk =
x̄1 − x̄2√
V (x̄1)/0.1n+ V (x̄2)/0.5n
→ N(0, 1). (2.5)
Desta maneira, valores extremos de Gk indicam falta de convergência. A técnica de
Geweke está implementada no pacote CODA (Best et al. (1995)), executável no software
R (R Development Core Team (2006)).
Após a obtenção da amostra, deve-se analisar a autocorrelação existente entre θ(j)
e θ(j−1). A amostra obtida a partir de uma cadeia de Markov é aleatória, mas não é
independente. Isso não afeta as estimativas dos parâmetros, mas tem influência sobre as
variâncias das estimativas resultantes desse procedimento de amostragem Gamerman e
Lopes (2006). Assim, nos casos em que for constatada uma forte correlação serial na ca-
deia, após verificada a convergência, recomenda-se a retirada de uma amostra sistemática
de seus valores para compor uma nova amostra. A forma como a amostragem sistemática
será realizada pode ser baseada em um gráfico contendo a função de autocorrelação da
cadeia.
2.2.2 WinBugs
O pacote estat́ıstico WinBUGS é uma versão em ambiente Windows do pacote BUGS
(Bayesian Inference Using Gibbs Sampling). É utilizado para análise bayesiana de mo-
delos estat́ısticos simples ou complexos, tendo a capacidade de estimar seus parâmetros
via MCMC. O WinBUGS consiste em um conjunto de funções que permitem a especi-
ficação do modelo e das distribuições de probabilidade para todos os seus componentes
12
aleatórios. Foi implementado por Thomas et al. (1992) e amplamente discutido em Lunn
et al. (2000).
O WinBUGS possui a capacidade de reconhecer formas de distribuições conjuga-
das, distribuições log-côncavas, distribuições com amplitudes restritas e etc. Com base
nesta informação, o algoritmo de amostragem mais eficiente é selecionado para simulação.
Quando nenhuma destas propriedades é identificada, uma mensagem avisa a incapacidade
na escolha do método de atualização.
Dentro do WinBugs existe uma ordenação dos métodos de amostragem dispońıveis
para serem utilizados, que depende da forma da distribuição de interesse. Primeiramente,
métodos de amostragem utilizando algoritmos padrões serão utilizados caso a distribuição
condicional seja conjugada. Caso essa condição não seja satisfeita, o amostrador de
Gibbs passa a ser utilizado: a ARS (Adaptive Rejection Sampling) é usada para amostrar
eficientemente qualquer distribuição condicional com função densidade log-côncava e a
ARMS (Adaptive Rejection Metropolis Sampling) generaliza a rotina ARS para o caso de
funções que não são log-côncavas, mas que possuem amplitudes restritas. Para o caso de
funções que não são log-côncavas e que não possuem amplitudes restritas, são utilizados
passos de Metropolis.
Para o algoritmo Metropolis-Hastings, o pacote usa como densidade de transição
q(θ(j), .) uma distribuição gaussiana centrada no valor atual do parâmetro θ(j).
Todo o processo inferencial utilizado neste trabalho foi implementado no software
WinBUGS versão 1.4.
13
Caṕıtulo 3
Seleção de Variáveis
Com frequência, em estudos aplicados, a modelagem estat́ıstica envolve um grande
número de regressores. Este problema acaba por trazer dificuldades na estimação do
modelo. Por exemplo, preditores relacionados de forma exata ou aproximada geram
dificuldades de estimação. Também pode-se citar o problema de obtenção de estimativas
imprecisas ou até mesmo não significativas para o modelo.
Eventualmente, também pode-se lidar com aplicações em que a quantidade de regres-
sores p é maior que n, número de observações. Um exemplo deste tipo, apresentado em
West (1993), consiste em prever o teor de gordura da massa de um determinado biscoito.
As caracteŕısticas desta massa são medidas por uma técnica chamada NIR (near infrared
spectroscopy). Os preditores são p = 300 ńıveis de reflectância mensurados, obtidos pela
técnica NIR, com uma amostra de 39 massas de biscoito.
Para esse exemplo, o método de mı́nimos quadrados não tem a capacidade de fazer a
estimação eficiente do modelo. Devido ao grande número de variáveis independentes, são
grandes as chances delas possúırem relações lineares exatas ou aproximadamente exatas
entre si, gerando o problema de multicolinearidade. Além disso, com tantas variáveis
regressoras, a variância associada aos parâmetros regressores pode ser muito alta e a
matriz X’X−1 intratável. Uma solução para a estimação de modelos em que p > n, seria
a obtenção de um pequeno número de combinações lineares do conjunto de variáveis
independentes, que retenham o máximo da informação contida nessas variáveis. Essa
técnica é conhecida como componentes principais, e frequentemente usada para cuidar
14
de multicolinearidade. Em geral, esses procedimentos são feitos em duas etapas: primeiro
obtem-se as componentes principais e depois a regressão estimada. Em West (1993), um
método integrado é apresentado sob a ótica bayesiana.
Em contextos onde p é uma quantidade muito grande, com o objetivo de evitar a
estimação de modelos complexos, alguma forma de redução de dimensionalidade, no que
diz respeito a quantidade de regressores p, é necessária. Com efeito, suponha o seguinte
modelo de regressão:
y = Xβ + �, (3.1)
onde β = (β1, . . . , βp)
T é o vetor paramétrico, y é um vetor n × 1 da variável resposta,
X é a matriz n× p dos regressores, e � é o vetor de erros de dimensão n× 1; as hipóteses
do modelo de regressão estabelecem que esses erros seguem uma distribuição normal, são
independentes e identicamente distribúıdos, com média 0 e variância desconhecida σ2.
Buscar soluções esparsas para o modelo de regressão em questão, é o mesmo que
identificar de maneira eficiente os coeficientes βp que são iguais a zero ou muito próximos
de zero. Logo, o regressor correspondente ao parâmetro βp = 0 ficará fora do modelo de
regressão, levando a redução da dimensão de p.
A partir de uma perspectiva bayesiana, existem duas principais abordagens para
a estimação da esparsidade associada aos regressores: misturas discretas e prioris de
contração (shrinkage). A primeira abordagem associa a cada βp uma distribuição a
priori que possui um ponto de massa no valor βp = 0 e uma alternativa absolutamente
cont́ınua; a segunda abordagem, que será utilizada nesta dissertação, modela cada βp
com distribuições a priori de contração, centradas em zero. Essas prioris são obtidas a
partir de misturas cont́ınuas.
Na seção 3.1, será apresentada uma técnica que introduz uma variável latente do tipo
Bernoulli na distribuição a priori de β, tal variável sinalizará os preditores que deverão
ser inclúıdos ou não no modelo. Enquanto que nas seções 3.2 e 3.3, serão apresentadas
técnicas para a estimação do modelo que utilizam distribuições de contração,obtidas
via misturas cont́ınuas, para o vetor paramétrico β. Todo procedimento de inferência
necessário nas técnicas a serem apresentadas será feito sob o enfoque bayesiano, isto é,
serão atribúıdas distribuições a priori para os parâmetros de interesse a fim de obter
15
a distribuição a posteriori, que em nosso caso não é conhecida. Técnicas de simulação
estocástica (MCMC) serão utilizadas para a obtenção de amostras desta distribuição.
3.1 Seleção de Variáveis via Busca Estocástica
O SSVS (Seleção de Variáveis via Busca Estocástica em inglês), é a técnica de seleção
de variáveis proposta por George e Robert (1993), a qual baseia-se na incorporação da
regressão em um modelo hierárquico de mistura de normais, onde um vetor de variáveis
latentes é capaz de sinalizar quais os melhores subconjuntos deX1, . . . , Xp. Cabe ressaltar
que os p regressores associados a y, fazem com que tenhamos 2p posśıveis modelos a serem
estimados.
Um fato interessante associado a este método de seleção de variáveis é que ele ”vi-
sita”mais vezes os modelos mais relevantes, no sentido de possúırem os regressores mais
apropriados para explicar a quantidade y. A estimação do modelo é posśıvel a partir da
seguinte estrutura hierárquica a priori para os parâmetros da regressão:
y|X,β, σ2 ∼ Nn(Xβ, σ2In)
βp|γp ∼ (1− γp)N(0, τ 2p ) + γpN(0, c2pτ 2p ) (3.2)
γp ∼ Bern(πp)
σ2 ∼ IG (ν, λ) ,
onde 0 ≤ πp ≤ 1, τp > 0, cp ∈ <, ν > 0 e λ > 0 são quantidades de ”sintonia”, isto é,
quantidades que precisam ser determinadas pelo pesquisador.
A quantidade πp pode ser interpretada como a probabilidade a priori de inclusão do
regressor Xp no modelo. Logo, πp = 0 indica que, a priori, o pesquisador assume que o
regressor Xp deve ser exclúıdo no modelo; de maneira análoga, quando πp = 1 assume-se
que o respectivo regressor deve ser inclúıdo no modelo estimado.
Um dos recursos da técnica Ssvs é que cada componente do vetor β é modelada
como mistura de distribuições normais com diferentes variâncias, conforme apresentado
na equação 3.2. Quando γp = 0, βp ∼ N(0, τ 2p ), indicando que esta componente βp deve
16
ser exclúıda do modelo. Logo a quantidade τp, que é o desvio-padrão da componente
βp neste caso particular onde γp = 0, deve ser determinada de tal maneira que o valor
estimado para esse parâmetro βp possa ser substitúıdo por 0. Desta maneira, o regressor
Xp correspondente estará exclúıdo do modelo. Analogamente, se βp ∼ N(0, c2pτ 2p ), então
a respectiva componente γp = 1. Nesse caso, estamos interessados na escolha de valor
para cpτp que nos leve a uma estimativa não nula para βp, fazendo com que o regressor
Xp seja inclúıdo no modelo estimado.
Portanto, pode ser observado que o ajuste dos parâmetros de sintonia τ e c não é
tarefa fácil. Diferentes escolhas para essas quantidades devem ser testadas. A Figura
3.1 ilustra como que distintas escolhas para τp e cp afetam a distribuição a priori de βp
quando γp = 0 e 1, respectivamente. De acordo com a Figura 3.1(a), podemos observar
que a distribuição de βp quando γp = 1 concentra uma grande massa de probabilidade
em torno do valor zero. Isso não é o desejável, uma vez que valores de γp = 1 sugerem
que a estimativa para o parâmetro βp seja não nula. A análise da distribuição de βp
quando γp = 0, representada pela linha cheia, concentra menos massa de probabilidade
em torno do valor zero, quando comparada com a linha tracejada, que é a distribuição do
parâmetro quando γp = 1. Neste caso, o ideal é a estimativa do parâmetro ser zero, mas
o ajuste do valor τp tal que a distribuição fique muito concentrada em torno do valor zero
deve ser evitado, afim de evitar que o parâmetro tenha uma distribuição muito restritiva.
Na Figura 3.1(b), ilustramos o comportamento para a distribuição de βp quando
γp = 0 e a estimativa para βp deve ser zero. Observamos que a função densidade da
distribuição a priori concentra uma grande massa de probabilidade em torno do valor
zero. Analogamente, quando γp = 1 e a estimativa para o parâmetro deve ser não
nula, observamos que a função de distribuição para βp é mais vaga e atribúı massa de
probablidade a valores mais distantes de zero.
A terceira combinação que apresentamos para a distribuição a priori para βp é ilus-
trada na Figura 3.1(c), onde observamos que quando γp = 0 e a estimativa para β deve
ser zero, vemos que a função densidade da distribuição a priori concentra uma massa
de probabilidade elevada em torno do valor zero. Em contrapartida, a linha tracejada
mostra que, quando γp = 1 e e a estimativa para βp deve ser não nula, a função de den-
17
sidade do parâmetro atribúı massa de probabilidade a valores mais afastados de zero. A
situação ilustrada por esta Figura consiste em um caso similar ao apresentado na Figura
3.1(b), a diferença está em prioris para βp quando γp = 0 ou 1 que atribuem massa de
probabilidade em intervalos com maiores amplitudes, ou seja, têm um comportamento
mais vago.
18
β
D
en
si
da
de
γ = 0
γ = 1
−10 −5 0 5 10
0.
0
0.
1
0.
2
0.
3
0.
4
0.
5
(a) τ = 2, c = 0.5
β
D
en
si
da
de
γ = 0
γ = 1
−10 −5 0 5 10
0.
0
0.
1
0.
2
0.
3
0.
4
0.
5
(b) τ = 2, c = 5
β
D
en
si
da
de
γ = 0
γ = 1
−10 −5 0 5 10
0.
0
0.
1
0.
2
0.
3
0.
4
0.
5
(c) τ = 10, c = 5
Figura 3.1: Diferentes configurações da distribuição a priori para β.
Como veremos nas aplicações referentes a esta dissertação, utilizamos a distribuição
a priori para o parâmetro βp que possui o comportamento mais vago em torno do valor
zero, quando a respectiva componente γp = 1. Esta configuração para o parâmetro βp foi
utilizada afim de evitar a concentração da alta massa de probabilidade em torno do valor
zero, quando a estimativa para a componente βp deverá ser não-nula. Especificamente,
as configurações apresentadas na Figura 3.1(b) e 3.1(c) foram utilizadas em distintas
aplicações.
Por fim, para as quantidades ν e λ, George e Robert (1993) ressaltam que a escolha de
ν próximos de 0 e qualquer valor de λ podem ser utilizadas para representar ignorância
a respeito do parâmetro σ2. Usaremos a configuração em que os parâmetros ν e λ são
iguais a 0.001. Note que, utilizando tais valores, consideramos a distribuição a priori
vaga, permitindo que os dados tenham maior influência na distribuição a posteriori.
Dentre as vantagens do Ssvs, podemos citar a possibilidade do usuário determinar
a importância prática de alguma variável regressora. Com efeito, suponha que estudos
teóricos mostrem que Xp é extremamente relevante para explicar y. Com o Ssvs, o usuário
pode levar em conta tal relevância, ao determinar que o elemento γp é Bernoulli com seu
parâmetro πp próximo a 1.
Uma das desvantagens do método está associada a grande quantidade de parâmetros
de ”sintonia”presentes na equação 3.2. O ajuste adequado dos parâmetros πp, τp, cp, ν
e λ dificulta o processo de estimação, uma vez que diferentes valores fixados podem
19
alterar drasticamente as estimativas obtidas. Nesta dissertação, o processo de ajuste dos
parâmetros foi simplificado quando fixamos que πp = π, τp = τ e cp = c, para todos p
regressores do modelo.
Quando a quantidade de regressores p > n, o Ssvs não é capaz de estimar o mo-
delo. Tal limitação é provavelmente decorrente, devido a um passo dentro do amostrador
de Gibbs, onde são necessárias as estimativas para β obtidas via mı́nimos quadrados.
Essa limitação é uma grande desvantagem comparativa aos métodos mais modernos de
estimação de modelos com uso de técnicas de seleção de preditores, que se baseiam na
obtenção de distribuições de contração para o vetor paramétrico β, e que serão apresen-
tados nas próximas seções. Tais métodos produzem estimativas coerentes até mesmo nos
casos multidimensionais.
Finalmente, muitasvezes a interpretação dos resultados obtidos pelo Ssvs é restringida
somente ao conhecimento dos melhores subconjuntos de X1, . . . , Xp, analisando-se apenas
a contagem dos modelos mais frequentementes visitados, e não as estimativas de β obtidas
pelo método. Neste trabalho não estaremos interessados em avaliar qual o modelo foi mais
frequentemente ”visitado”, e sim, nas estimativas para β, permitindo assim a comparação
do Ssvs com os demais métodos que serão apresentados. Além disso, pode ser observado
que o Ssvs indica e seleciona os modelos mais frequentes, enquanto que os métodos que
serão apresentados a seguir, naturalmente fazem uma mistura de modelos.
3.2 Operador de Seleção e Contração com Penali-
dade em Valor Absoluto
Dentre os métodos que fazem tanto a contração cont́ınua quanto a seleção de variáveis,
uma técnica promissora que utiliza o operador de seleção e contração com penalidade em
valor absoluto, foi proposta por Tibshirani (1996). Essa técnica será denominada como
Lasso, que sintetiza least absolute shrinkage and selection operator, isto é, operador de
seleção e contração mı́nimo absoluto, em português.
O Lasso é um penalizador do procedimento de mı́nimos quadrados, que minimiza a
20
soma dos quadrados dos reśıduos com uma restrição na norma L1 dos coeficientes β’s.
Assim, a estimativa de β sob o método do Lasso é dada por:
β̂ = arg min
β
(ỹ −Xβ)′(ỹ −Xβ) + λ
p∑
j=1
|βj|. (3.3)
Observe que X é a matriz dos regressores padronizados, a quantidade ỹ = y − ȳ1n e
λ é um parâmetro de ”sintonia”.
Uma maneira de ilustrar o funcionamento do estimador Lasso, é no caso onde temos
somente dois preditores. O losângulo da figura 3.2 caracteriza a restrição imposta pelo
método Lasso na estimação de β, enquanto que as elipses são as curvas de ńıveis das
estimativas de verossimilhança de βp. As curvas de ńıvel poderão interceptar o losângulo
em um de seus quatro vértices. A solução para o estimador Lasso corresponde a inter-
ceptação de uma dessas elipses com o losângulo. Se isto ocorrer no vértice (como na
Figura 3.2) a estimativa de um dos parâmetros será nula, caso contrário representará um
contração dessa estimativa em torno do valor zero.
21
ββi
ββ k
 0.
02
 
 0.
04
 
 0
.0
6 
 0.08 
 0.
1 
 0
.1
2 
 0.14 
 0.1
6 
 0.
18
 
 0.2 
 0.
22
 
0.24
 
−2 −1 0 1 2 3 4 5
−
2
0
2
4
6
Figura 3.2: Restrição na estimação dos β’s imposta pelo Lasso bayesiano.
Com a presença do termo penalizador λ
∑p
j=1 |βj| na equação 3.3, Tibshirani (1996)
nota que a estimativa do Lasso pode ser interpretada como a moda a posteriori es-
timada quando os parâmetros β’s da regressão possuem distribuição a priori Laplace
(exponencial dupla) independentes. Uma vantagem desta distribuição é que ela pode
ser expressada como uma mistura na escala de distribuições normais com variâncias que
seguem distribuições exponenciais independentes.
3.2.1 Formulação Hierárquica do Modelo Lasso Bayesiano
Recentemente, Park e Casella (2008) propuseram o amostrador de Gibbs para o Lasso,
a partir da seguinte estrutura hierárquica para o modelo:
22
y|X,β, σ2 ∼ Nn(Xβ, σ2In)
p(β|σ2, τ 21 , . . . , τ 2p ) ∼ N(0p, σ2Dτ ) (3.4)
Dτ = diag(τ
2
1 , . . . , τ
2
p )
σ2, τ 21 , . . . , τ
2
p ∼ π(σ2)dσ2
p∏
j=1
λ2
2
e−λ
2τ2j /2dτ 2j
Especificamente, foi considerada uma análise bayesiana completa a partir do uso da
distribuição a priori Laplace, condicionada a σ2, no modelo hierárquico. Com efeito, tal
distribuição é da seguinte forma:
π(β|σ2) =
p∏
j=1
λ
2
√
σ2
e−λ|βj |/
√
σ2 (3.5)
Esta especificação condicional a σ2 é particularmente importante, uma vez que garante
que a distribuição conjunta π(β, σ2|ỹ) seja unimodal, segundo Park e Casella (2008). Uma
consequência da não unimodalidade é a dificuldade de convergência do amostrador de
Gibbs.
Como dito anteriormente, esta distribuição Laplace pode ser expressada como uma
mistura na escala de distribuições normais com variâncias que seguem distribuições ex-
ponenciais independentes, isto é:
a
2
e−a|z| =
∫ ∞
0
1√
2πs
e−z
2/(2s)a
2
2
e−a
2s/2ds, a > 0, (3.6)
onde temos a variável aleatória Z|s ∼ N(0, s) combinada com S ∼ Exp(a2/2). Essa
representação foi exemplificada em Andrews e Mallows (1974) e é exatamente a mistura de
normais na escala obtida a partir do Lasso bayesiano. Desta maneira, o Lasso bayesiano é
uma metodologia de estimação que utiliza distribuições de contração, obtidas via misturas
cont́ınuas, conforme pode ser constatado na equação 3.6.
Note que a representação da distribuição Laplace como uma mistura de normais na
escala exibida na equação 3.6 é facilmente obtida:
23
f(βp|σ2) =
∫ ∞
0
f(βp|σ2, τ 2p )f(τ 2p )dτ 2p
f(βp|σ2) =
∫ ∞
0
1√
2πσ2τ 2p
e
−1
2σ2τ2
i
β2p λ2
2
e
−λ2τ2p
2 dτ 2p (3.7)
f(βp|σ2) =
λ
2
√
σ2
e
−λ|βp|√
σ2
Como os βp’s e τp’s são independentes, a distribuição de β|σ2 é obtida pelo produtório
de cada uma das f(βp|σ2), chegando ao resultado descrito em (3.5).
Para a especificação do modelo hierárquico associado ao Lasso bayesiano, é necessária
a especificação das distribuições a priori associadas aos parâmetros σ2 e λ, presentes na
equação 3.4. Nesta dissertação, utilizaremos a distribuição a priori Inversa Gama para
o parâmetro σ2, como recomendado em Park e Casella (2008).
Sob a perspectiva bayesiana, o parâmetro λ pode ser estimado através do procedi-
mento bayesiano emṕırico ou pelo uso de uma distribuição a priori apropriada. Nesta
dissertação o parâmetro será estimado a partir da especificação de uma distribuição a
priori, embora o procedimento bayesiano emṕırico tenha sido utilizado em trabalhos
anteriores. Park e Casella (2008) consideram o uso de uma distribuição gamma a pri-
ori para λ2, uma vez que a conjugação resultante permite que o amostrador de Gibbs
seja mais facilmente especificado. Deve-se evitar a especificação de prioris vagas, como
(p(λ2) ∝ 1/λ2), uma vez que a distribuição a posteriori resultante será imprópria. O
ideal é que p(λ2) se aproxime de 0 suficientemente rápido quando λ2 → ∞, sendo ao
mesmo tempo relativamente vaga.
Nas aplicações, recomenda-se a padronização da matriz de covariáveis X.
3.2.2 Função de Contração
Para uma melhor compreensão a respeito do método do Lasso, podemos definir uma
quantidade, função dos parâmetros, denomida parâmetro de contração. Este parâmetro,
será representado pela quantidade κi = 1/(1 + τ
2
i ).
Com efeito, suponha que o nosso objetivo seja estimação do seguinte modelo:
24
yi|βi, σ2 ∼ N(βi, σ2) (3.8)
βi|τ 2i , σ2 ∼ N(0, τ 2i σ2)
Quando fixamos a quantidade σ2 = 1, o valor esperado a posteriori do parâmetro βi
fica definido por:
E(βi|yi, τ 2i ) =
1
1 + τ 2i
0 +
τ 2i
1 + τ 2i
yi = (1− κi)yi (3.9)
É importante ressaltar que no modelo proposto em (3.8), para cada observação yi
temos uma estimativa βi associada ao valor. Assim, a quantidade de parâmetros β’s a
serem estimados é exatamente igual ao tamanho da amostra n. Esse caso é diferente do
modelo apresentado na equação 3.4, onde existem variáveis independentes associadas aos
p preditores β’s.
Voltando a Equação 3.9, observamos que o parâmetro de contração κi pode ser inter-
pretado como a quantidade de peso que a média a posteriori de β concentra no ponto 0.
Note que, valores de κi próximos a zero fazem com que a média a posteriori de β seja
o próprio valor observado yi, indicando que não houve contração do parâmetro β. Por
outro lado, valores de κi próximos a um, fazem com que a média a posteriori de β seja
o valor zero, representando a contração total do parâmetro estimado.
Uma vez que a quantidade κi ∈ [0, 1], podemos eliminar a condicionalidade associada
ao parâmetro τ 2i da seguinte forma:
E(βi|y) =
∫ 1
0
(1− κi)yip(κi|y)dκi = [1− E(κi|yi)]y. (3.10)O núcleo da função de densidade do parâmetro κi, associado ao Lasso bayesiano é
apresentada a seguir:
p(κ) ∝ exp
(
−1
2κ
)
κ−2 (3.11)
O cálculo desta função densidade está apresentado no Anexo desta dissertação. O
gráfico desta função de densidade associada ao Lasso bayesinao é apresentado na Figura
3.3:
25
κi
D
en
si
da
de
 a
 m
en
os
 d
e 
co
ns
ta
nt
es
0.0 0.2 0.4 0.6 0.8 1.0
0.
0
0.
5
1.
0
1.
5
2.
0
2.
5
Figura 3.3: Densidade de κi ∈ [0, 1] associado ao método de estimação Lasso bayesiano.
Com a inspeção da função de distribuição a priori do parâmetro de contração κi, é
posśıvel analisar de maneira mais clara como o método faz o discernimento entre as ob-
servações associadas aos rúıdos e aquelas que são associadas a valores espúrios. A Figura
3.3 ilustra a função de densidade do parâmetro de contração para o Lasso bayesiano; é
posśıvel notar que a massa de probabilidade concentrada em torno do valor zero é pe-
quena, isso indica que a probabilidade desse parâmetro de contração ser igual a zero é
pequena. Como vimos anteriormente que a esperança a posteriori de β é igual a (1−κi)yi,
é posśıvel concluir que raramente o Lasso associa ao valor de β o próprio valor observado
yi.
Analisando o comportamento da função de distribuição do parâmetro de contração
onde κi é próximo ao valor um, é posśıvel notar que a distribuição é limitada. Em outras
palavras, a probabilidade do parâmetro κi ser igual a um não é tão alta. Ressaltando
que altas probabilidades de κi = 1 indicam a capacidade de contração do parâmetro β.
Conclúımos que o Lasso possui uma capacidade limitada de contração da estimativa do
parâmetro β.
26
3.2.3 Função de Influência
Outro instrumento que permite um melhor entendimento sobre o comportamento do
Lasso bayesiano como método de estimação de modelos e seleção de preditores, é a função
de influência. Como seu próprio nome sugere, o estudo desta função permite a análise da
maneira como os dados serão tratados pelo método de estimação. Nesta análise, temos
o particular interesse em analisar o comportamento dessa função em valores associados
as observações espúrias.
Um resultado básico e necessário para o cálculo da função de influência foi apresentado
em Pericchi e Smith (1992) e merece ser revisitado:
Suponha que x1, . . . , xn seja uma amostra aleatória de uma distribuição normal com
média β e variância σ2. Logo y =
∑
xi/n ∼ N(β, σ2/n) tem distribuição de densidade
p(y|β). Defina a quantidade m(y), dada por:
m(y) =
∫
p(y − β)π(β)dβ. (3.12)
Este resultado é aplicável para qualquer função de distribuição a priori para o vetor β
que obedeça a condição π(β) ≥ 0 em valores de β pertencentes ao conjunto dos números
reais. Também defina as seguintes quantidades:
s(y) = −∂{log(m(y)}
∂y
e S(y) = −∂{log(s(y)}
∂y
(3.13)
A função de influência é dada pela quantidade s(y). Como dito anteriormente, o
estudo desta função irá auxiliar a compreender a maneira como os dados são tratados
pelo método de estimação.
Com efeito, suponha que y ∼ N(β, 1), com β = 0. Neste caso, é fácil observar que
a função de influência associada a este modelo é dada por s(y) = y. A Figura 3.4, que
será apresentada mais adiante, possui o gráfico relativo a esta função. Sua análise mostra
que valores pequenos, tem uma pequena influência no modelo, ao contrário de valores
grandes, que associam uma grande influência ao modelo, indicando que este modelo não
é robusto com relação as observações espúrias.
27
Pericchi e Smith (1992) mostram que tanto a esperança quanto a variância a posteriori
de β podem ser escritas como função das quantidades apresentadas na equação 3.13,
portanto:
E(β|y) = y + σ
2
n
s(y) e V ar(β|y) = σ
2
n
+
σ2
n2
S(y). (3.14)
Ainda neste estudo, Pericchi e Smith (1992) apresentam os valores das quantidades de
interesse: s(y), a função de influência, e E(β|y), a esperança a posteriori do parâmetro,
quando a distribuição a priori de β é exponencial dupla. Este é exatamente o caso do
método de estimação do Lasso bayesiano, que associa ao parâmetro β a distribuição
exponencial dupla, da seguinte forma:
p(β) =
1√
2σ2
exp
[
− |β|√
σ2
]
.
Essa é a função de distribuição exponencial dupla associada ao Lasso bayesiano quando
λ = 1. Neste contexto, os valores das quantidade de interesse são dados por:
s(y) = − a
σ2
[F (y)−G(y)] onde,
F (y) = exp[c(y)]Φ
[√
(n)
σ
(−y − b)
]
, (3.15)
G(y) = exp[−c(y)]Φ
[
−
√
(n)
σ
(−y + b)
]
,
a = exp[
1
n
] , b =
√
2
n
, c(y) =
√
2y
σ2
Em que Φ(.) denota a função de distribuição acumulada normal padrão. A média a
posteriori de β pode ser obtida a partir da seguinte expressão
E(β|y) = w(y)(y + b) + [1− w(y)](y − b) , onde (3.16)
w(y) =
F (y)
F (y) +G(y)
(3.17)
28
Na Figura 3.4 é posśıvel observar o comportamento da função de influência associada
ao Lasso bayesiano. Tal função tem a caracteŕıstica de truncagem da influência determi-
nadas observações, assim, observações associadas as observações espúrias possuem uma
influência constante na estimação do modelo nesta metodologia. Observe que a com-
paração da função de influência do Lasso bayesiano com a função associada ao modelo
normal, exalta a diferença dos modelos no que diz respeito ao tratamento das observações
espúrias.
29
Y
F
un
çã
o 
de
 In
flu
ên
ci
a
−10 −5 0 5 10
−
2
−
1
0
1
2
Figura 3.4: Funções de influência associadas ao modelo normal e Lasso, linha tracejada
e cheia, respectivamente.
Dentre as vantagens da utilização deste método na estimação de um modelo, citamos
a presença do termo penalizador λ. A restrição imposta por esse termo se mostra como
uma qualidade interessante na proposta de contração das estimativas de β. Tal proce-
dimento só é posśıvel a partir da idéia do Lasso. Comparado com o Ssvs, notamos que
o tempo computacional para a estimação de um mesmo modelo é bem menor. Ainda
verificamos a presença de menos termos de sintonia. Em contrapartida, o ajuste de uma
distribuição a priori adequada para λ2 pode ser uma tarefa delicada, uma vez que é ideal
que essa distribuição não seja muito vaga, para não haver o risco do amostrador de Gibbs
fornecer estimativas imprecisas. Uma desvantagem do método é a sensibilidade associada
a distribuição de λ2.
30
3.3 Mistura de normais na Escala Usando Distri-
buições de Cauchy
A estimação de modelos via mistura na escala de distribuições normais com a distri-
buição de Cauchy é umas das técnicas mais recentemente apresentadas no contexto do
uso de distribuições de contração, obtidas via misturas cont́ınuas. Ressaltando que, mis-
turas cont́ınuas para a obtenção de distribuições de contração também foram utilizadas
para a estimação do modelo pelo Lasso bayesiano.
3.3.1 Formulação Hierárquica do Modelo
O estimador de modelos via mistura na escala de normais com distribuições Cauchy,
será definido daqui em diante como estimador Horseshoe. A estimação do modelo via
Horseshoe é um método eficiente não só de estimação, mas também de seleção de predi-
tores no modelo proposto. Sua metodologia foi proposta em Carvalho et al. (2010). O
método de estimação via Horseshoe assume que cada um dos parâmetros βp’s possuem
distribuições condicionalmente independentes dado λ, o parâmetro de contração global.
Dessa maneira, o modelo de estimação é definido com a seguinte mistura na escala de
normais:
y|X,β, τ, λ, σ2 ∼ Nn(Xβ, σ2In)
βp|τp ∼ N(0, τ 2p ) (3.18)
τp|λ ∼ C+(0, λ)
λ|σ ∼ C+(0, σ)
onde C+(0, a) é uma distribuição Cauchy padrão truncada nos reais positivos, com
parâmetro de escala a.
As quantidades τp’s podem ser interpretadas como parâmetros de contração local,
no sentido de estarem associadas a cada βp.Observe que na estrutura hierárquica im-
posta pelo método de estimação via Horseshoe, só precisamos fixar os valores dos hi-
31
perparâmetros associados a variância do modelo, σ2. Os demais parâmetros são devida-
mente estimados a partir da estrutura hierárquica imposta. A distribuição a priori para
o parâmetro de variância σ2 será a priori de Jeffrey’s. Assim p(σ2) ∝ 1/σ2, tendo sua
distribuição relativamente vaga e permitindo que os dados tenham maior influência na
distribuição a posteriori do parâmetro de variância.
O estimador Horseshoe tem a liberdade de fazer a contração dos elementos de βp de
maneira global, através do parâmetro λ, e de maneira local através das quantidades τp.
O parâmetro λ estima o ńıvel de esparsidade associado ao vetor paramétrico, enquanto
que os parâmetros de contração locais são capazes de reduzir os valores associados ao
vetor paramátrico β. Essa caracteŕıstica é uma vantagem do método Horseshoe quando
comparado aos demais métodos de seleção de preditores, já que nenhum outro tem essa
mesma capacidade.
A Figura 3.5 ilustra o comportamento da distribuição a priori para o parâmetro β.
Tal distribuição é obtida a partir da mistura no parâmetro de escala da distribuição
normal associada a β|τ , com a distribuição Cauchy truncada nos valores reais positivos.
A função de distribuição para o parâmetro β é limitada da seguinte maneira:
(2π2)−1/2
2
log
(
1 +
4
β2
)
< p(β) < (2π2)(−1/2)log
(
1 +
2
β2
)
,
conforme demonstrado em Carvalho et al. (2010).
32
β
D
en
si
da
de
−3 −2 −1 0 1 2 3
0.
0
0.
1
0.
2
0.
3
0.
4
0.
5
(a)
β
D
en
si
da
de
3 4 5 6 7
0.
00
0.
01
0.
02
0.
03
0.
04
(b)
Figura 3.5: Comparação entre as distribuição a priori para β. As linhas tracejada e cheia
representam a distribuição associada ao método Lasso e Horseshoe, respectivamente.
A Figura 3.5 ilustra o comportamento das distribuições a priori associadas aos métodos
de estimação de modelo Horseshoe e Lasso. Podemos observar no gráfico 3.5(a) que a
distribuição Horseshoe possui um alongamento nos valores onde β é próximo de zero. Tal
comportamento é a chave para a boa performance do método de estimação com relação
aos rúıdos associados ao vetor paramétrico β. Este comportamento é menos percebido
quando estamos analisando a distribuição exponencial dupla, associada ao método de
estimação via Lasso bayesiano. Na Figura 3.5(b) podemos observar que a distribuição
Horseshoe apresenta sua cauda mais pesada, quando comparada com a cauda da distri-
buição exponencial dupla. É exatamente essa caracteŕıstica das caudas pesadas, que faz
com que o método de estimação através do Horseshoe lide melhor com os valores espúrios
associados ao modelo.
3.3.2 Função de Contração
Na seção 3.2.2 foi apresentado o parâmetro de contração, definido pela quantidade
κi = 1/(1 + τ
2
i ). A função deste parâmetro associa uma regra de contração local ao
33
método de estimação de modelo em análise, uma vez que o parâmetro τi está diretamente
relacionado com a variabilidade a priori do vetor paramétrico β.
Assim como feito no Lasso bayesiano, para o método de estimação pelo Horseshoe
também calculamos a função de densidade do parâmetro de contração. O cálculo desta
função é apresentado no Anexo deste trabalho. A função de distribuição do parâmetro
κ é dada pela seguinte expressão:
p(κ) ∝ (κ)−0.5(1− κ)−0.5 (3.19)
A Figura a seguir ilustra o comportamento da função de distribuição do parâmetro
de contração associado ao método de estimação de modelo via Horseshoe:
34
κ
D
en
si
da
de
0.0 0.2 0.4 0.6 0.8 1.0
0
0.
5
1
1.
5
2
2.
5
Figura 3.6: Densidade de κi ∈ [0, 1] associado ao método horseshoe para a estimação.
Diferentemente da função de contração associada ao método de estimação pelo Lasso
bayesiano, neste caso observamos que a densidade da função p(κi) é ilimitada em valores
de κi próximos a 0 e 1, indicando que probabilidades elevadas são associadas a κi ≈ 0 e
1.
Devemos relembrar que o parâmetro de contração pode ser interpretado como a quan-
tidade de peso que a média a posteriori de β concentra no ponto 0. Além disso, como
vimos anteriormente que E(β|y) = (1−κi)yi, então, valores de κi próximos a zero, fazem
com que esta média a posteriori seja igual ao próprio valor observado yi, indicando que
não houve a contração do parâmetro β. Enquanto que valores de κi próximos a um,
fazem com que a média a posteriori do parâmetro β seja o valor zero, indicando a total
contração do parâmetro estimado.
A fato da função de contração associada ao método de estimação via Horseshoe ser
ilimitada em valores de κi ≈ 0, indica a alta probabilidade do parâmetro κi ser igual a
zero. Isso sugere que na estimação do modelo pelo método Horseshoe, a capacidade que o
modelo tem de avaliar a importância da observação, sem fazer a contração do parâmetro
relacionado a esta observação, é elevada. Assim, o método Horseshoe de estimação capta
35
bem as observações importantes do modelo, e não contráı desnecessariamente o parâmetro
associado a esta observação.
A análise do comportamento da função de contração em valores de κi próximos a
um, indica que há uma probabilidade alta associada a chance do parâmetro κi = 1.
Relembrando que, quando o parâmetro de contração é igual a um, então a esperança a
posteriori de β será igual a zero. Portanto, o fato da função de contração ser ilimitada
em valores próximos a um, indica que o método de estimação via Horseshoe possui uma
grande capacidade de contrair a estimativa do parâmetro β.
3.3.3 Função de Influência
Também podemos analisar o comportamento da função de influência. A função de
influência foi um conceito introduzido anteriormente e seu estudo permite compreen-
der melhor como o método de estimação irá tratar os dados. O estudo dessa função é
particularmente interessante para o dados associados a rúıdos e valores espúrios.
A média a posteriori do parâmetro β é dada por:
E(β|y) = y
[
1− 2Φ(1/2, 1, 3/2, y
2/2, 1− 1/τ 2)
3Φ(1/2, 1, 5/2, y2/2, 1− 1/τ 2)
]
, (3.20)
onde Φ(α, β, γ, x, y) é a função hipergeométrica degenerada de duas variáveis. A
obtenção do valor da média a posteriori do parâmetrico β é apresentada de maneira
detalhada em Carvalho et al. (2010).
A função de influência para o modelo Horseshoe foi calculada a partir da média a
posteriori apresentada anteriormente, e seu gráfico é exibido na Figura a seguir:
36
Figura 3.7: Função de Influência para o modelo Horseshoe.
O comportamento da função de influência associada ao método de estimação pelo
Horseshoe é apresentada na Figura 3.7. Podemos observar que esta função de influência
rendescende para o valor zero. Este comportamento indica que este método de estimação
tem a capacidade de diminuir a influência dos valores espúrios associados ao modelo.
Ainda na Figura 3.7, a linha tracejada apresenta a função de influência associada ao
modelo normal. É posśıvel notar que a influência da observação no modelo é seu próprio
valor observado. Enquanto que a função de influência associada ao método de estimação
do Lasso bayesiano, tende a fixar o valores da influência a partir de um determinado valor
constante observado. Carvalho et al. (2010) ressaltam que nenhum método de estimação
conhecido, tem a capacidade de reduzir a influência de observações associadas a valores
espúrios no modelo, assim como o Horseshoe faz. Portanto, essa caracteŕıstica é uma
grande vantagem deste método de estimação.
Dentre as principais vantagens associadas a este método, podemos citar que a im-
plementação do método Horseshoe para a estimação de um modelo, não exige que hi-
37
perparâmetros sejam ajustados.Diferentemente do Lasso e do método Ssvs, que podem
ter suas estimativas severamente comprometidas caso seus hiperparâmetros sejam mal
ajustados, o Horseshoe não corre esse risco.
Além disso, o método Horseshoe apresenta grandes vantagens teóricas quando compa-
rado diretamente com o Lasso. A distribuição a priori de β possui caudas bem pesadas,
que permitem um tratamento mais apropriado as posśıveis observações espúrias associa-
das ao modelo.
Uma análise com foco no contexto de seleção de preditores nos mostra que este método
também possui propriedades mais desejáveis, visto que ele possui uma capacidade mais
apurada de encolher os valores estimados de β. A contração total de valores associados
ao rúıdo é a palavra-chave no contexto de seleção de preditores, já que esta caracteŕıstica
é a que permite a redução da dimensão do vetor paramétrico β.
38
Caṕıtulo 4
Critérios de Seleção de Modelos
No caṕıtulo 3 foram apresentados métodos distintos que estimam e selecionam de
maneira eficiente os parâmetros para um modelo linear. Porém, a utilização de distintas
técnicas para a estimação de um mesmo modelo, acaba por gerar dúvidas sobre qual
deles obteve o melhor desempenho, e consequentemente, qual a técnica de estimação é a
mais adequada para a estimação. Neste contexto, insere-se a aplicação de critérios para
seleção de modelos, baseados em medidas que tentam quantificar sua qualidade de ajuste
e aplicam alguma penalidade a modelos complexos, com muitos parâmetros, apontando
para a escolha daqueles que sejam mais parcimoniosos.
Apresentaremos a seguir critérios para a seleção de modelos baseados no cálculo da
verossimilhança marginal, como o fator de Bayes. Também utilizaremos o critério de
informação baseado no desvio; adicionalmente, apresentaremos uma proposta baseada na
função de perda canônica. Finalmente, apresentaremos os critérios baseados em função
de perda quadrática e valor absoluto.
39
4.1 Critérios Baseados na Função de Verossimilhança
Marginal
4.1.1 Fator de Bayes
O fator de bayes é um critério que se baseia na comparação de modelos via função
de verossimilhança marginal. Seguindo a definição de Kass e Raftery (1995), os dados
y são melhores representados sob determinadas especificações, que aqui serão represen-
tadas pelas hipóteses H0 e H1, desta maneira, a função de verossimilhança marginal é
representada pela densidade p(y|H0) ou p(y|H1). Sabemos pelo teorema de bayes, que:
p(Hk|y) =
p(y|Hk)p(Hk)
p(y|H0)p(H0) + p(y|H1)p(H1)
, k = 0, 1. (4.1)
As densidades p(H0) e p(H1) são as distribuições a priori de cada uma das hipóteses.
A partir da equação 4.1, as seguintes relações são obtidas:
(i) FB(H1, H0) =
p(y|H1)
p(y|H0)
.
(ii)
p(H1|y)
p(H0|y)
=
p(y|H1)p(H1)
p(y|H0)p(H0)
, isto é, a razão de chances a posteriori é igual a razão de
chances a priori versus o fator de Bayes.
Analisando a igualdade apresentada em (i), é posśıvel concluir que, quando as hipóteses
H0 e H1 são igualmente prováveis a priori, teremos p(H0) = p(H1) = 0.5. Desta maneira,
o fator de bayes será igual a razão de chances a posteriori contra a hipótese H0. Além
disso, sabemos que p(H0|y) = 1− p(H1|y), então, uma nova relação, baseada no fator de
bayes, é obtida:
p(H1|y) =
FB(H1, H0)
1 + FB(H1, H0)
(4.2)
Portanto, a partir da equação 4.2, é posśıvel obter a probabilidade da hipótese H1 ser
a especificação adequada para os dados y.
40
Nos casos em que θ é um vetor paramétrico, o fator de bayes ainda será a razão das
verossimilhanças, porém, a função p(y|H) será obtida através da integração no espaço
paramétrico, desta maneira:
p(y|Hk) =
∫
p(y|θk, Hk)π(θk|Hk)dθk, (4.3)
onde θk é o vetor que contempla todos os parâmetros envolvidos na hipótese Hk.
Quando utilizamos métodos MCMC para a estimação dos parâmetros, algumas alter-
nativas podem ser utilizadas para estimar a equação 4.3. Algumas delas são a média
harmônica e a verossimilhança ponderada calculada via Bootstrap, ambas porpostas por
Newton e Raftery (1994), além da média harmônica generalizada apresentada por Gel-
fand e Dey (1993). Nesta dissertação utilizaremos as médias harmônicas para a estimação
dessa equação. Tais métodos serão apresentados mais a frente.
Com base na proposta de Jeffreys (1961), Kass e Raftery (1995) sugerem a seguinte
interpretação para o fator de bayes:
Tabela 4.1: Interpretação do fator de Bayes
FB(H1, H0) p(H1|y) Evidência contra H0
1 a 3 0.5 a 0.75 Nenhuma Evidência
3 a 20 0.75 a 0.95 Postiva
20 a 150 0.95 a 0.99 Forte
> 150 > 0.99 Muito Forte
A seguir, discutiremos alternativas para a obtenção da equação 4.3.
Média Harmônica
O estimador da média harmônica foi apresentado por Newton e Raftery (1994) e é
obtido através da seguinte expressão:
41
p̂mh(y) =
[
1
G
G∑
g=1
1
f(y|θ(g))
]−1
, (4.4)
onde θ(g) é o vetor paramétrico simulado da distribuição a posteriori, obtido na g-
ésima iteração do método MCMC.
Uma das vantagens do estimador da média harmônica é a facilidade com a qual ele
pode ser obtido: só precisamos conhecer a função verossimilhança associada ao modelo e
obter os valores simulados pelas G iterações do MCMC. Por outro lado, uma desvanta-
gem associada ao estimador é sua instabilidade, uma vez que é senśıvel a valores muito
pequenos da função de verossimilhança.
Média Harmônica Generalizada
O estimador da média harmônica generalizada foi proposto por Gelfand e Dey (1993)
e é uma extensão do estimador da média harmônica. Seu valor é obtido por:
p̂mhg(y) ==
[
1
G
G∑
g=1
g(θ(g))
f(y|θ(g))f(θ(g))
]−1
, (4.5)
onde g(θ(g)) é uma função de importância que deve ser escolhida cuidadosamente. Nt-
zoufras (2009) pondera que uma uma distribuição normal p-variada ou uma distribuição
t-student p-variada, com média e variância iguais as médias e variâncias a posteriori
estimadas pelo MCMC usualmente fornecem um bom estimador da média harmônica
generalizada.
Para esta dissertação, a função de importância g(θ(g)) é distribuição t-student p-
variada com poucos graus de liberdade, com média e variância iguais as médias e variâncias
estimadas a posteriori.
4.1.2 Escores Logaŕıtmicos
O escore logaŕıtmico é um critério de seleção de modelo que é facilmente obtido através
das sáıdas do MCMC. Seu valor é calculado a partir de regras de escore estritamente
42
próprias, como descrito em Gneiting e Raftery (2004). Eles consideram regras escore para
obter bondades de ajuste que num contexto bayesiano estão relacionadas com particulares
funções de utilidades.
Assuma que o escore médio seja dado por:
Sn(θ) =
1
n
n∑
i=1
S(p(y|θ)), (4.6)
onde S(.) é uma regra escore própria, assim, o maior escore será obtido para o
verdadeiro modelo. Gneiting e Raftery (2004) apresentam o escore logaŕıtmico, que é
dado por:
S(p(y|θ)) = log[p(yrep = yi|y)],
onde p(yrep = yi|y) denota a função de densidade preditiva a posteriori no ponto yrep = yi
para o modelo que está sendo avaliado. Quando a amostra das iterações obtidas via
MCMC para o vetor paramétrico θ está dispońıvel, uma aproximação para a quantidade
log[p(yrep = yi|y)] para a i-ésima observação pode ser obtida:
ˆlog[p(yrep = yi|y)] =
1
G
G∑
g=1
log[p(yi|θ̂
(g)
)], (4.7)
no qual θ(g) é a g-ésima sáıda do MCMC para o vetor paramétrico.
Maiores valores de escores logaŕıtmicos indicam um melhor ajuste do modelo.
4.2 Critério de Informação Baseado no Desvio - DIC
O critério DIC (Deviance Information Criterion) foi apresentado por Spiegelhalter
et al. (2002) como uma medida de adequalibilidade e comparação de modelos. O DIC é
dado pela seguinte expressão:
DIC = D̄ + pD, (4.8)
43
onde D̄ é a deviance, que é dada por: D̄ = 1
G
∑G
g=1−2logL(θ(g)), sendo que G é a
quantidade de amostras