ads:
NAYARA FRANCINE DE MOURA GONÇALVES
BOOTSTRAP EM MODELOS AUTO-REGRESSIVOS
ADITIVOS GENERALIZADOS
UNIVERSIDADE FEDERAL DE MINAS GERAIS
MAIO 2009
ads:
Livros Grátis
http://www.livrosgratis.com.br
Milhares de livros grátis para download.
II
BOOTSTRAP EM MODELOS AUTO-REGRESSIVOS
ADITIVOS GENERALIZADOS
NAYARA FRANCINE DE MOURA GONÇALVES
Orientadora: GLAURA DA CONCEIÇÃO FRANCO - UFMG
Co-Orientador: VALDÉRIO ANSELMO REISEN - UFES
Dissertação apresentada ao Departamento de Estatística do Instituto de
Ciências Exatas da Universidade Federal de Minas Gerais como requisito
parcial à obtenção do título de MESTRE em ESTATÍSTICA.
Maio 2009
ads:
III
Ah, se o mundo inteiro me
pudesse ouvir! Tenho muito pra contar,
dizer que aprendi...
Azul da cor do mar - Tim Maia
IV
Agradecimentos
Aos meus pais e irmão
Desde o começo vocês estiveram comigo. Inúmeras foram as vezes que sacrificaram
seus objetivos em favor dos meus. Também não foram raros os momentos em que
buscaram meu sorriso e me encontraram tão cheia de pressa. Eu agradeço por todo
amor que recebi, pois em todas as lições de vida vocês estavam presentes e sempre
que foi preciso decidir vocês acreditaram em mim!
À Glaura
Você me convidou a voar na sua sabedoria e o que eu aprendi foi que este voar sempre
dependeu das minhas próprias asas. Agradeço pela amizade, paciência e pela tão
dedicada orientação.
Ao Prof. Valdério e à minha amiga Fabiana
Pelas fundamentais contribuições.
Ao Prof. Paulo Sérgio Lúcio (Departamento de Estatística da UFRN)
Por gentilmente ceder o banco de dados reais utilizado neste trabalho.
V
Resumo
A classe dos Modelos Aditivos Generalizados (MAG), considerados uma extensão
dos Modelos Lineares Generalizados, vem atraindo a atenção de pesquisadores
principalmente em função de sua flexibilidade. Apesar de construído sob a hipótese de
independência dos dados, os MAG’s são muito aplicados em estudos de séries
temporais, sobretudo como alternativa para modelagem de variáveis de confusão tais
como tendência e sazonalidade. Recentemente, modelos mais gerais, que consideram a
estrutura de correlação entre os dados, como os modelos GLARMA (
autoregressive
moving average generalized linear models
), têm sido utilizados. Este trabalho estende
os modelos GLARMA para uma classe de modelos auto-regressivos aditivos
generalizados para séries de contagem cuja distribuição condicional, dadas as
observações passadas e as variáveis explicativas, segue uma distribuição de Poisson.
Além de apresentar uma conceituação desses modelos bem como procedimentos
de ajustes, este trabalho emprega, em um estudo empírico, o procedimento
bootstrap
em três formas (
bootstrap
nas observações,
bootstrap
condicional e
bootstrap
nos
resíduos)
na inferência pontual dos parâmetros do modelo e compara dois métodos de
construção de intervalos de confiança
bootstrap
-
bootstrap
percentílico e
bootstrap
com correção do vício – na estimação intervalar.
Os resultados mostram que, em geral, os procedimentos e os intervalos de
confiança
bootstrap
apresentam um bom desempenho quando utilizados na classe de
MAG’s que por sua vez, quando auxiliados pela modelagem GLARMA, modelam bem
dados de contagem com estrutura auto-regressiva de ordem 1, apresentando
estimativas próximas dos valores verdadeiros dos parâmetros.
VI
Abstract
The class of Generalized Additive Models (GAM), considered an extension of the
Generalized Linear Models (GLM), is attracting the attention of researchers mainly due
to the flexibility of these procedures. In spite of being built under the hypothesis of
independency of the data, the GAM is widely applied to time series data, as an
alternative to model variables such as trend and seasonality. Recently, more general
models, which consider the correlation structure among the data, like the GLARMA
models (autoregressive moving average generalized linear models), are being used.
This work extends the GLARMA models to a class of autoregressive generalized additive
models of count series whose conditional distribution, given the past observations and
the independent variables, follows a Poisson distribution.
Besides presenting the definition of the model, as well as the fitting procedures,
this work employs, in a empirical study, the bootstrap procedure in three different ways
(bootstrap in the observations, conditional bootstrap and the bootstrap in the residuals)
in the interval inference of the parameters, comparing two bootstrap methods of
building confidence intervals – percentile bootstrap and bootstrap with bias correction.
The results show that, in general, the procedures and the bootstrap confidence
intervals present a satisfactory performance when used in the GAM models with the
GLARMA structure, modeling count data with an autoregressive structure of order 1,
and presenting estimates close to the true values of the parameters.
VII
Sumário
1. INTRODUÇÃO ...................................................................................................... 1
1.1 Revisão de Literatura ................................................................................... 1
1.2 Objetivos ..................................................................................................... 4
1.3 Organização do Trabalho .............................................................................. 5
2. TÉCNICAS DE SUAVIZAÇÃO .................................................................................. 7
2.1 Definição e propriedades dos suavizadores .................................................... 7
2.2 Seleção do parâmetro de suavização ........................................................... 10
2.3 O SUAVIZADOR
loess
................................................................................. 13
3. MODELOS ADITIVOS GENERALIZADOS................................................................ 17
3.1 Modelos Lineares Generalizados .................................................................. 17
3.2 Modelos Aditivos Generalizados .................................................................. 19
3.2.1 Ajuste dos modelos não paramétricos ................................................ 20
3.2.2 Ajuste dos modelos semiparamétricos ............................................... 24
3.2.3 Função Desvio .................................................................................. 26
3.2.4 Seleção do parâmetro de suavização ................................................. 27
4. MODELOS AUTO-REGRESSIVOS GENERALIZADOS ............................................... 29
4.1 Modelos Poisson auto-regressivo média móvel linear generalizados ............... 29
4.2 Modelos Poisson auto-regressivos aditivos generalizados .............................. 32
5. TÉCNICA BOOTSTRAP ........................................................................................ 35
5.1
Bootstrap
nas observações ......................................................................... 35
5.2
Bootstrap
não paramétrico nos resíduos ...................................................... 36
5.3
Bootstrap
condicional ................................................................................. 37
5.4 Intervalos de Confiança
bootstrap
............................................................... 38
5.4.1 Intervalos de confiança
bootstrap
percentílico .................................... 38
5.4.2 Intervalos de confiança
bootstrap
com correção do vício ..................... 39
6. ANÁLISE DOS DADOS SIMULADOS ...................................................................... 41
6.1 Resultados das simulações de dados independentes .................................... 42
VIII
6.2 Resultados das simulações de séries temporais ............................................ 43
7. APLICAÇÃO A SÉRIES REAIS ............................................................................... 47
7.1 Análise descritiva ....................................................................................... 48
7.2 Modelagem MAG ........................................................................................ 52
7.3 Modelagem MAG-AR(1) .............................................................................. 53
8. CONCLUSÕES RELEVANTES ................................................................................ 56
9. REFERÊNCIAS BIBLIOGRÁFICAS.......................................................................... 58
1
1. INTRODUÇÃO
1.1 Revisão de Literatura
Durante muitos anos os modelos lineares, sob a suposição de normalidade,
foram utilizados para descrever fenômenos aleatórios. Se o fenômeno sob estudo não
apresentasse uma resposta para a qual fosse razoável supor a distribuição Gaussiana,
algum tipo de transformação da variável era utilizado, com o propósito de se alcançar a
normalidade.
Mais tarde, pesquisadores abriram o leque de opções para a distribuição da
variável resposta permitindo que a mesma pertencesse à família exponencial de
distribuição. Em 1972, Nelder e Wedderburn unificaram estes procedimentos
introduzindo a classe dos Modelos Lineares Generalizados (MLG). Muitas distribuições
conhecidas pertencem à família exponencial, como a Normal, a Poisson, a Binomial e a
Gama. Uma característica dos MLG’s é que a forma da relação funcional entre a média
da variável resposta e as variáveis preditoras é completamente linear e especificada por
termos paramétricos.
Como a suposição de linearidade pode ser irrealista em situações práticas, Hastie
e Tibshirani (1990) propuseram os Modelos Aditivos Generalizados (MAG), cuja principal
diferença com os MLG’s é a substituição da usual forma linear das covariáveis por
funções suavizadoras não paramétricas que sumarizam a associação entre a variável
resposta e as variáveis explicativas.
Existem na literatura várias técnicas de suavização – Hastie e Tibshirani (1990)
apresentam algumas delas. Em qualquer uma das técnicas, a suavização é obtida
ajustando-se uma curva aos dados de tal forma que a cada ponto, a curva dependa
somente das observações naquele ponto (
0
x
) e em uma vizinhança.
Entre as formas mais simples de suavização estão as técnicas
running-mean
e a
running-line
(Hastie e Tibshirani, 1896). A primeira, também conhecida como
moving-
average
, é popular por ser utilizada em dados temporais igualmente espaçados,
entretanto, apesar da simplicidade dos cálculos, esta técnica tende a ser viciada por
suavizar a tendência nos pontos iniciais e finais da série.
2
O suavizador r
unning-line
, uma simples generalização do suavizador
running-
mean
, reduz o problema do vício utilizando os mínimos quadrados lineares ao invés da
média em cada vizinhança. Por outro lado, a técnica tende a gerar curvas bastante
irregulares então um segundo estágio de suavização geralmente é necessário. Uma
maneira de melhorar os resultados em um suavizador
running-line
é usar o ajuste de
mínimos quadrados ponderados (MQP) em cada vizinhança. O suavizador
running-line
pode produzir resultados irregulares porque pontos em uma determinada vizinhança
têm pesos iguais (não nulos), enquanto pontos fora desta região têm peso zero.
Em 1977, Cleveland abranda este problema propondo uma técnica mais robusta
de suavização,
locally weighted running line smoother
(
loess
). O procedimento, que
pode ser usado para suavizar dados com configurações mais gerais que não sejam
necessariamente para dados de séries temporais igualmente espaçados, é uma
adaptação do ajuste sucessivo de modelos de regressão pelo método de mínimos
quadrados ponderados (Beaton and Tukey, 1974; Andrews, 1974). A proposta é dar
maior peso para
0
x
e para os pontos da vizinhança e pesos que decrescem suavemente
para pontos mais afastados desta região. Visando tornar o ajuste
loess
ainda mais
robusto, Cleveland propõe, em 1979, o suavizador
robust locally weighted
running line
smoother
introduzindo um novo conjunto de pesos a serem utilizados no ajuste de
MQP.
Outras alternativas de técnicas de suavização são propostas na literatura.
Introduzida inicialmente por Whittaker (1923), a suavização por
spline
é uma técnica
bastante utilizada no ajuste do MAG; autores como Dominici
et al.
(2002) e Ramsay
et
al.
(2003), descrevem e utilizam as funções splines e algumas de suas extensões –
penalized splines e parametric splines – como forma de suavização; a técnica Kernel é
destacada, por exemplo, em Silverman (1986), Buja
et al.
(1989) e Härdle(1990).
Nenhuma das técnicas de suavização faz alguma suposição paramétrica sobre a
curva a ser estimada; cada suavizador tem um parâmetro que determina o “quanto” a
função estará suavizada. Para a função
loess
, por exemplo, este parâmetro é chamado
span
. Para
splines
e
parametric splines
o grau de suavização é especificado através do
número de graus de liberdade.
Também fazem parte da classe de MAG’s os modelos semiparamétricos –
segundo Buja
et al.
(1989) pode-se entender estes modelos como uma generalização
dos MAG’s – , constituídos pela soma de termos paramétricos de algumas variáveis
3
preditoras e funções suavizadoras de outras. Discutido por Stone (1985), Hastie e
Tibshirani (1990) e Lee (1990), este modelo tem sido bastante explorado na literatura
devido a esta peculiaridade de compor a flexibilidade do modelo aditivo não
paramétrico com um componente paramétrico. Uma contribuição de interesse é feita
por Buja
et al.
(1989) ao apresentar diversas técnicas de estimação aplicáveis a este
tipo de modelo. Além desses trabalhos destacam-se outros como Schick (1986,
1993,1996) e Bhattacharya e Zhao (1997) que enriqueceram a teoria assintótica acerca
do modelo semiparamétrico apresentando assim uma ponte entre os resultados
assintóticos para o modelo puramente não paramétrico já discutidos por Stone (1977),
Cox (1983) ou Rice e Rosenblatt (1983).
Um grande número de pesquisadores vem utilizando os MAG’s e sua extensão
semiparamétrica em séries temporais, principalmente em estudos cujo objetivo é
avaliar os efeitos da poluição atmosférica sobre a saúde de seres humanos (Schwartz
et
al.
, 1993). Alguns fatores como condições meteorológicas e os dias da semana
influenciam os dados (Díez, 1999) e confundem a associação entre a exposição de
interesse e o desfecho. Além disso, ainda existem as componentes da própria série
temporal como tendência, sazonalidade e autocorrelação. O MAG tem sido uma técnica
alternativa facilitadora no controle desses fatores de confusão já que esta modelagem
elimina a necessidade de especificar uma forma paramétrica para a associação entre
covariáveis e preditor.
Nestes estudos epidemiológicos a variável resposta é, geralmente, alguma
contagem de eventos que representam danos à saúde, como o número de óbitos ou o
número de internações por determinada causa respiratória. Na maioria das vezes, estes
desfechos são modelados com o pressuposto de que as contagens dos eventos seguem
uma distribuição Poisson – como exemplos de aplicações desta técnica ver estudos de
Schwartz
et al.
(1992), Conceição
et al.
(2001) e Lima
et al.
(2001).
Em 1999, Davis
et al.
faz uma revisão dos modelos já propostos na literatura
para séries temporais de contagem que seguem a distribuição Poisson. Em particular,
uma nova classe de modelo, GLARMA (
autoregressive moving average generalized
linear models
), é introduzida e suas propriedades desenvolvidas em parte. Por serem
capazes de capturar uma gama de estruturas de dependência entre as observações das
séries temporais, os modelos GLARMA também vêm sendo utilizados nestes estudos
epidemiológicos nos casos, diferentemente do MAG, em que a relação entre a variável
4
resposta e as covariáveis assume a forma linear e os dados apresentam uma estrutura
de dependência, por exemplo a auto-regressiva.
Um problema comumente encontrado na estimação dos termos dos modelos
semiparamétricos é a baixa freqüência de dados. Pequenas amostras e a dificuldade em
se determinar a distribuição assintótica dos dados fazem com que os métodos de
estimação percam um pouco de sua eficiência. Nestes casos o
bootstrap
, sugerido por
Efron (1979), pode ser utilizado para a melhoria das inferências intervalares.
O método
bootstrap
consiste em uma técnica de reamostragem que permite
aproximar a distribuição de uma função das observações pela distribuição empírica dos
dados, baseando-se em uma amostra finita (Efron e Tibshirani, 1993). Em série
temporal essa técnica tem sido bastante empregada, porém, o fato dos dados em série
não serem independentes torna a aplicação do método bastante criteriosa – as
observações não devem ser reamostradas diretamente, pois sua estrutura original pode
ser perdida. A utilização do
bootstrap
em modelos MAG ainda é pouco discutida na
literatura estatística, mas já considerada em alguns trabalhos recentes – Härdle
et al.
(2004) mostram como o procedimento pode ser utilizado na correção do vício das
estimativas paramétricas e não-paramétricas, em testes de hipótese e na construção de
bandas de confiança; Figueiras
et al.
(2005), utilizam um
bootstrap
condicional para
corrigir problemas de concurvidade em modelos MAG. Entretanto, não foram
encontrados trabalhos utilizando o
bootstrap
em modelos GLARMA.
Existem vários métodos
bootstrap
de construção de intervalos de confiança e
entre eles estão os intervalos
bootstrap
percentílico (Efron e Tibshirani, 1986) e o
bootstrap
com correção do vício (Efron e Tibshirani, 1986 e Hall, 1988) utilizados neste
estudo.
1.2 Objetivos
Este trabalho tem por objetivos (a) propor uma modelagem para dados de
contagem que substitua a usual forma linear do GLARMA pela estrutura
semiparamétrica dos MAG’s – neste caso, além dos termos lineares serão consideradas
também funções não lineares de variáveis explicativas na construção dos modelos; (b)
comparar a estimação dos parâmetros lineares da modelagem proposta com a
estimação do MAG em estudos simulados de dados independentes e de séries
5
temporais, ambos com respostas que seguem a distribuição Poisson e (c) utilizar a
técnica
bootstrap
com o intuito de fazer inferências sobre os parâmetros lineares dos
modelos.
1.3 Organização do Trabalho
O abundante uso da técnica MAG em séries temporais, apesar de não existirem
resultados teóricos sobre sua utilização em dados dependentes, e a presença de
estruturas de correlação nas séries, motivaram a proposta deste trabalho que é
estender os modelos GLARMA para uma classe de modelos auto-regressivos aditivos
generalizados, MAG-AR, ainda não abordada na literatura.
No decorrer do texto, são apresentados os modelos generalizados MAG e
GLARMA bem como algumas de suas propriedades. Alguns resultados gerais sobre
suavizadores lineares também serão apresentados, no entanto o foco maior é dado ao
suavizador
loess
– escolhido por apresentar diferentes propriedades estatísticas e
também por ser bastante utilizado em estudos temporais sobre os efeitos da poluição
atmosférica na saúde dos seres humanos; a citada técnica será utilizada no ajuste da
parte não-paramétrica do MAG e MAG-AR. Três diferentes abordagens
bootstrap
–
bootstrap
nas observações,
bootstrap
condicional e
bootstrap
nos resíduos – são
avaliadas.
Para atender aos objetivos acima, estudos de simulação Monte Carlo serão
realizados para que se compare, segundo o vício e o erro quadrático médio, o
comportamento dos estimadores do parâmetro linear dos modelos. Além disto,
intervalos de confiança obtidos através da técnica
bootstrap
, como o intervalo
bootstrap
percentílico e o
bootstrap
com correção do vício, serão comparados ao intervalo de
confiança assintótico quanto ao percentual de cobertura e tamanho dos intervalos.
Este trabalho está organizado da seguinte maneira: o Capítulo 2 apresenta as
técnicas de suavização descrevendo com maiores detalhes o suavizador
loess
. No
Capítulo 3 apresenta-se o modelo MAG dissertando sobre sua estimação e
propriedades. A descrição do modelo GLARMA e a proposta do novo modelo MAG-AR
são feitas no Capítulo 4. No Capítulo 5 o procedimento
bootstrap
é abordado.
6
Resultados das simulações e da análise de dados reais são discutidos, respectivamente,
nos Capítulos 6 e 7 e o Capítulo 8 conclui o trabalho.
7
2. TÉCNICAS DE SUAVIZAÇÃO
Um suavizador é uma ferramenta utilizada para descrever a dependência de uma
variável resposta
Y
como uma função de uma ou mais variáveis preditoras
X
. Uma
importante propriedade de um suavizador é sua natureza não paramétrica: a forma da
curva estimada é determinada pelos próprios dados e, de fato, não é necessário nem
mesmo conhecer previamente a forma dessa relação para estimá-la, por esse motivo os
procedimentos de suavização são também denominados técnicas de regressão não
paramétrica. Os valores destas funções devem ser mais “suaves” do que os valores de
Y
, isto porque a estimativa obtida de um procedimento de suavização tem
variabilidade menor que a de
Y
, daí a razão do nome suavizador ou alisador.
2.1 Definição e propriedades dos suavizadores
Formalmente, um suavizador é definido como uma função de
),...,(
1 n
xxx
=
e
),...,(
1 n
yyy
=
,
)
|
(
x
y
f
δ
=
, com mesmo domínio de
x
. Para alguns suavizadores, a
função )|( xyf
δ
=
calculada em
0
x
, )(
0
xf , é definida para todo
0
x
. Outras vezes ela é
definida apenas para os valores observados
n
xx ,...,
1
e, neste caso, algum tipo de
interpolação é necessário para obter estimativas para outros valores de
X
.
As curvas suavizadas podem ser utilizadas com diferentes objetivos; tipicamente
são empregadas no ajuste do modelo
)()|(
ii
xfxXYE
=
=
, (2.1)
uma generalização dos modelos de regressão linear simples, onde
f
é uma função
arbitrária desconhecida não especificada a priori.
Os suavizadores são classificados como lineares ou não lineares. Um suavizador
é linear se
)|()|()|(
2121
xybxyaxbyay
δ
δ
δ
+
=
+
(Hastie e Tibshirani, 1990).
8
Focando no vetor de valores estimados
fxfxfyyy
nn
ˆ
))'(
ˆ
),...,(
ˆ
()'
ˆ
,...,
ˆ
(
ˆ
11
===
, um
suavizador linear pode ser escrito como nas equações (2.2) e (2.3)
=f
ˆ
Sy
(2.2)
nnxxxi
ysysysxf
iii
,22,11,
...)(
ˆ
+++=
,
n
i
,...,
1
=
‘(2.3)
onde
}{}{
ijjx
ssS
i
=
=
é uma matriz de dimensão
n
n
×
chamada matriz suavizadora
que não depende de
y
, mas apenas de
X
e da técnica de suavização adotada, e
nxxx
iii
sss
,2,1,
,...,,
é a i-ésima linha da matriz
S
. Pelo fato de suas matrizes não
dependerem da variável resposta, a análise por meio destes suavizadores é
relativamente simples (Hastie e Tibshirani 1989).
Dado um algoritmo de um suavizador linear pode-se produzir a correspondente
matriz suavizadora
S
suavizando o vetor de base unitária: o resultado de suavização
do i-ésimo vetor unitário é a i-ésima coluna de
S
(Buja
et al.
1989).
Running-mean
,
running-line
,
spline
e
loess
são exemplos de suavizadores lineares.
Um simples exemplo de suavizador não linear para o qual a matriz suavizadora
não pode ser construída é o suavizador
running median
. A diferença entre esta técnica
e o suavizador
running mean
é que o valor ajustado através do
running median
para
0
x
e vizinhança é dado pela mediana destas observações. Neste caso, para variáveis
aleatórias
X
e
Y
, mediana(
Y
X
+
)
≠
mediana (
X
) + mediana (
Y
), assim as
estimativas dependem de
Y
de um modo não linear.
A matriz suavizadora desempenha papel semelhante ao da matriz
hat
1
no
método de estimação de mínimos quadrados e algumas de suas propriedades são
demonstradas por Hoaglin e Welsch (1978).
.
i
10
≤
≤
ii
s
;
ii
11
≤
≤
−
ij
s
para
j
i
≠
;
1
No modelo de regressão normal linear, a estimativa do vetor de médias
βµ
ˆ
ˆ
X=
, onde
yXXX ')'(
ˆ
1
−
=
β
, pode ser
reescrita como
Hy=
µ
ˆ
, com
')'(
1
XXXXH
−
=
. A matriz
H
é a matriz de projeção ortogonal dos vetores de
n
R
no subespaço gerado pelas colunas da matriz
X
.
9
iii
1
=
ii
s
se e somente se
0
=
ij
s
para todo
j
i
≠
;
iv
∑
=
=
n
j
ij
s
1
1
.
Porém, para algumas técnicas de suavização como
running line
e
loess
S
não é
simétrica e nem idempotente, ou seja, não é um operador de projeção como ocorre
com a matriz
hat
. Para os suavizadores
Bin
,
least-square line
,
polynomial regression
e
splines
, por exemplo, a matriz suavizadora é simétrica e seus autovalores são sempre
reais.
Se o suavizador é linear, o estimador
f
ˆ
para uma função arbitrária
f
é sempre
viciado
SffSYEffEfffE −=−=−=− )()
ˆ
()
ˆ
(
(2.5)
e a matriz de covariância de
f
ˆ
é dada por
')()()
ˆ
( SYSVarSYVarfVar ==
(2.6)
e, sob a suposição de que os
i
Y
’s são independentes com
2
)(
σ
=
i
YVar
, a expressão
(2.6) pode ser reescrita
2
')
ˆ
(
σ
SSfVar =
. (2.7)
Considerando estas informações, o erro quadrático médio
n
xfxfE
EQM
i
n
i
i
2
1
)](
ˆ
)([
∑
=
−
=
(2.8)
assume a forma
(
)
∑∑
==
+=
n
i
i
n
i
i
b
n
xfVar
n
EQM
1
2
1
1
)(
ˆ
1
10
(2.9)
n
bb
n
SStr ')'(
2
+=
σ
,
onde
b
é o vetor de vício definido em (2.5). O parâmetro
2
σ
em (2.7) e em (2.9)
geralmente é desconhecido. Um estimador para este parâmetro, assumindo que
f
é
não viciado é dado por
(
)
)'2(
ˆ
'
ˆ
)'2(
)(
ˆ
ˆ
2
2
SSStrn
ee
SSStrn
xfy
ii
−−
=
−−
−
=
∑
σ
, (2.10)
onde
ySISyyfye )(
ˆ
ˆ
−=−=−=
, com
I
representando uma matriz identidade
n
n
×
. No caso em que
S
é idempotente,
)
(
)
(
)
'
(
S
rank
S
tr
SS
tr
=
=
e
∑
=
−=−
n
i
ii
SSStr
1
2
)2()'2(
θθ
, sendo
i
θ
,
n
i
,...,
1
=
, os autovalores de
S
.
O resultado (2.10) é apresentado sob a suposição de ausência de viés de
f
ˆ
.
Porém o viés é nulo apenas para uma classe restrita de funções. O suavizador
loess
,
por exemplo, fornece funções
f
ˆ
viciadas para uma
f
arbitrária. Uma solução para
esse problema é considerar o comportamento assintótico discutidos amplamente por
Stone (1977), Cox (1983) ou Rice e Rosenblatt (1983).
2.2 Seleção do parâmetro de suavização
Na maioria das técnicas lineares de suavização, o valor alisado é obtido segundo
o comportamento de uma vizinhança. Diferentes formas de cálculos nesta vizinhança
definem as diferentes técnicas de suavização. A escolha do tamanho da vizinhança está
associada a um parâmetro
λ
, denominado parâmetro de suavização. Sendo mais
importante até do que a escolha da técnica de suavização, a definição dos valores para
λ
é um passo importante no processo, isto porque o parâmetro está diretamente
relacionado à relação de ganho e perda entre o viés e a variância da curva estimada:
aumentar
λ
implica aumentar a suavização da curva, logo, diminui-se a variância, por
11
outro lado, perde-se informação no ajuste implicando no aumento do viés. A Figura 2.1
ilustra este efeito.
É interessante observar também a influência de
λ
nos termos do
EQM
(2.9);
de forma geral, aumentando
λ
, o
)
'
(
SS
tr
tende a diminuir e o vício tende a aumentar.
-2 -1 0 1 2 3
-1.5 -0.5 0.0 0.5 1.0
X
Y
Ajuste com lambda=0.7
-2 -1 0 1 2 3
-1.5 -0.5 0.0 0.5 1.0
X
Y
Ajuste com lambda=0.5
-2 -1 0 1 2 3
-1.5 -0.5 0.0 0.5 1.0
X
Y
Ajuste com lambda=0.3
-2 -1 0 1 2 3
-1.5 -0.5 0.0 0.5 1.0
X
Y
Ajuste com lambda=0.1
Figura 2.1 – Diagrama de dispersão de
X
e
Y
com curva suavizada pelo
método
loess
com diferentes valores de
λ
.
Não existe um critério rígido para a escolha do valor de
λ
. Na prática estes
valores são escolhidos a
priori
através da inspeção visual da curva ou através de um
método automático que geralmente testa vários valores de
λ
para um mesmo conjunto
de dados. A Validação Cruzada (
cross-validation
), por exemplo, é uma técnica
automática de seleção do parâmetro suavizador; o método consiste em retirar o ponto
12
),(
ii
yx
da base de dados e calcular
)(
ˆ
i
xf
apenas com os
1
−
n
pontos restantes,
n
i
,...,
1
=
. A estatística de validação cruzada é dada por
∑
=
−
−=
n
i
i
i
i
xfy
n
VC
1
2
))(
ˆ
(
1
)(
λ
λ
,
n
i
,...,
1
=
(2.4)
onde
)(
ˆ
i
i
xf
−
λ
indica o valor estimado para
Y
quando o ponto
),(
ii
yx
é eliminado. A
expressão (2.4) é calculada considerando-se um conjunto de valores de
λ
fixados a
priori
, sendo selecionado o valor de
λ
que minimize esta expressão. Alguns detalhes
podem ser encontrados em Silverman (1985) e Craven e Wahba (1979).
Apesar da validação cruzada e de outras formas de seleção automática para o
parâmetro de suavização parecerem bem fundamentadas, suas performances são
questionáveis. Hastie e Tibshirani (1990) mostraram, em um estudo simulado, que os
valores de
λ
assim obtidos apresentam grande variabilidade, mesmo para dados
gerados a partir de modelos simples, com pequena variância para os erros. Os autores
sugerem que a escolha deste parâmetro seja feita com métodos gráficos auxiliados por
medidas dos graus de liberdade dos suavizadores.
Dada uma matriz suavizadora
S
de um suavizador linear com um
λ
fixado, o
número dos graus de liberdade pode ser definido como
)
(
S
tr
gl
=
. (2.12)
Quanto maior o número de graus de liberdade, menor o valor de
λ
e, por
conseqüência, menor a quantidade de suavização. Existem ainda pelo menos outras
duas definições de graus de liberdade – ver em Hastie e Tibshirani (1990). Estas
definições, assim como (2.12) são derivadas da analogia dos modelos de regressão
linear e podem ser utilizadas com vários propósitos, entre eles comparar os
suavizadores descritos nas seções anteriores levando em conta à “quantidade” de
suavização ou comparar duas técnicas com base no valor esperado da soma de
quadrado residual (ver maiores detalhes em Buja
et al.
1989); além disso, qualquer
uma das expressões para graus de liberdade pode ser usada para auxiliar na escolha de
um valor para o parâmetro de suavização.
13
2.3 O SUAVIZADOR
loess
Proposto por Cleveland em 1977, o
loess
(l
ocally weighted running line
smoother
) é um método de suavização que se baseia no ajuste sucessivo de
n
modelos de regressão pelo método de mínimos quadrados ponderados (MQP).
Para cada ponto (
ii
yx ,
) define-se uma vizinhança e aos pontos (
kk
yx ,
) nessa
vizinhança são atribuídos pesos através de uma função
U
. Esses pesos são utilizados
no ajuste de um polinômio por MQP.
A vizinhança de cada (
ii
yx ,
) é constituída por
l
pares de observações (
kk
yx ,
)
que possuem as coordenadas
k
x
mais próximas a
i
x
. A quantidade
l
a ser
considerada é dada por
n
l
λ
=
(2.15)
onde
λ
, pertencente ao intervalo
]1,0(
, é o parâmetro de suavização correspondente
à proporção do número total de observações a ser utilizado em cada ajuste local.
O grau do polinômio deve ser fixado com base no padrão apresentado pelos
dados num diagrama de dispersão. De uma forma geral, se a nuvem de pontos sugere
uma tendência sem máximos ou mínimos locais, então um ajuste linear é adequado. Se
existem regiões com máximos ou mínimos locais, então um ajuste quadrático
geralmente produz uma curva que descreve localmente melhor o padrão dos dados.
A função
U
que atribui pesos em cada ajuste local do polinômio, tendo (
ii
yx ,
)
como ponto alvo tem a forma
))((
1
, kiikx
xxhUu
i
−=
−
,
n
i
,...,
1
=
,
n
k
,...,
1
=
(2.16)
onde
i
h
é a distância entre
i
x
e seu
l
-ésimo vizinho mais próximo, com
l
definido em
(2.15), ou seja, colocando em ordem crescente as distâncias
||
ki
xx
−
,
n
i
,...,
1
=
,
n
k
,...,
1
=
,
i
h
é a distância que ocupa a
l
-ésima posição nesta seqüência. A função
U
deve ser especificada de forma a possuir as seguintes propriedades:
14
.
i
0
)
(
>
g
U
para
1
1
<
<
−
g
;
ii
)
(
g
U
é uma função par, isto é,
)
(
)
(
g
U
g
U
−
=
;
iii
)
(
g
U
é uma função decrescente para
0
≥
g
;
iii
0
)
(
=
g
U
para
1
|
|
=
g
A função tricúbica dada por
(
)
<−
=
contráriocaso
gparag
gU
0
1||||1
)(
3
3
(2.17)
apresenta as propriedades descritas acima e, de acordo com Cleveland (1979), fornece
uma suavização adequada na maioria dos casos.
Com base na função (2.17) obtém-se a matriz de pesos referente ao ponto alvo
(
ii
yx ,
)
=
i
x
U
diagonal
}u,...,u{
n,x1,x
ii
(2.18)
com elementos dados por
(
)
<−−−
=
−−
contráriocaso
xxhparaxxh
u
kiikii
kx
i
0
1|)(||)(|1
1
3
31
,
. (2.19)
Assim, por (2.18), em um ajuste local tendo como ponto alvo (
ii
yx ,
), este ponto fica
associado a um peso 1; os pesos diminuem à medida que os pontos se afastam de
(
ii
yx ,
) e pontos fora da vizinhança de
i
x
ficam associados a pesos nulos. A reta de
regressão assim ajustada
xy
βα
ˆ
ˆ
ˆ
+=
fornece o valor previsto
ii
xxf
βα
ˆ
ˆ
)(
ˆ
+=
.
Refazendo todos os passos considerando cada uma das
n
observações (
y
x
,
)
como ponto alvo, obtêm-se os pontos (
)(
ˆ
, xfx
) que formam a curva suavizada.
15
O
loess
é um suavizador linear; neste caso, os valores previstos de
Y
obtidos no
procedimento de suavização podem ser escritos da forma dada em (2.2). Os elementos
da matriz suavizadora
S
são denotados por
=
=
nnnn
x
x
x
nxxx
nxxx
nxxx
s
s
s
sss
sss
sss
S
'
'
'
...
...
...
2
1
222
111
,2,1,
,2,1,
,2,1,
MMOMM
O vetor
i
x
s'
, referente à i-ésima linha da matriz
S
, corresponde também à i-ésima
linha da matriz
iii
xxx
UXXUXXS ')'(
1−
=
, (2.20)
construída no ajuste da regressão ponderada local que tem (
ii
yx ,
) como ponto alvo e
matriz de pesos
i
x
U
definida em (2.18),
n
i
,...,
1
=
.
Dessa forma o valor previsto correspondente a
i
x
pode ser dado como em (2.3)
e ser reescrito como
ys
i
x
'
,
n
i
,...,
1
=
. (2.21)
De forma geral, pode-se mostrar que o elemento
ij
da matriz suavizadora do
loess
é dado por
,
)(
2
1
,
1
,
2
1
,
1 1
,
1
,
2
,,
1
,
2
,
,
−
++−
=
∑∑∑
∑ ∑∑∑
===
= ===
n
j
jxj
n
j
jxj
n
j
jx
n
j
n
j
jx
n
j
jxjjijxjjijx
n
j
jxjjx
jx
iii
iiiiii
i
uxuxu
uuxxxuxxxuuxu
s
(2.22)
onde
jx
i
u
,
é definido de acordo com (2.16).
16
A expressão 2.22 mostra claramente que os elementos de
S
dependem apenas
das covariáveis e do parâmetro de suavização. Quanto menor o
λ
, maior o número de
elementos iguais a zero na diagonal de
i
x
U
e maior o número de elementos nulos nas
linhas de
S
.
17
3. MODELOS ADITIVOS GENERALIZADOS
Os suavizadores podem ser utilizados em modelos de regressão com o objetivo
de descrever a relação entre a média da variável resposta e as variáveis preditoras.
Neste capítulo é feita uma breve introdução dos Modelos Lineares Generalizados (MLG)
e descritos os Modelos Aditivos Generalizados (MAG) – modelos de regressão que vêm
ganhado destaque na literatura, principalmente por sua flexibilidade – mostrando o
emprego dos suavizadores nesta classe.
3.1 Modelos lineares generalizados
Os Modelos Lineares Generalizados (MLG) são formados por um componente
aleatório, um componente sistemático e uma função de ligação que “liga” os dois
componentes. A resposta
Y
, componente aleatória do modelo, tem função de
densidade de probabilidade (ou função de probabilidade) dada por
{
}
),()]([exp);;(
Φ
+
Θ
−
Θ
=
Φ
Θ
ycbyy
Y
φ
ρ
(3.1)
onde
Θ
é chamado parâmetro natural e
Φ
parâmetro de dispersão,
)
(
⋅
b
e
)
(
⋅
c
são
funções especificadas,
)
(
⋅
b
é duas vezes diferenciável e
0
1
>
Φ
−
. Assume-se também
que a esperança de
Y
, denotada por
µ
, está relacionada às covariáveis
d
XX
,...,
1
por
η
µ
=
)
(
g
, onde
dv
XX
β
β
α
η
+
+
+
=
...
11
;
η
é a componente sistemática do modelo
linear generalizado chamada preditor linear e
)
(
⋅
g
é a função de ligação. Um caso
particular ocorre quando o preditor linear coincide com o parâmetro
Θ
, isto é,
Θ
η
=
;
neste caso, a função de ligação é chamada ligação canônica. Estas ligações
desempenham papel muito importante na teoria dos MLG’s e muitas vezes são
escolhidas por possuírem propriedades estatísticas e matemáticas convenientes (ver
Paula, 2000).
18
McCullagh e Nelder (1989) mostraram que se
Y
tem função densidade de
probabilidade dada por (3.1) seu valor esperado e sua variância estão relacionados ao
parâmetro natural da seguinte forma
Θ
Θ
µ
∂
∂
=
)(b
e (3.2)
VYVar
11
)(
−−
=
∂
∂
=
Φ
Θ
µ
Φ
,
onde
V
é chamada função de variância.
A estimação do vetor dos
1
+
d
parâmetros dos MLG,
β
ββ
β
=
)',...,,(
1 d
β
β
α
é feita
por máxima verossimilhança e a estimativa deste vetor é calculada resolvendo-se as
seguintes equações escore
0)(
1
1
=−
∂
∂
−
=
∑
iii
n
i
i
i
ij
yVx
µ
η
µ
,
d
j
,...,
1
=
(3.3)
onde
)(
ii
YVarV
=
e
1
0
=
i
x
.
As equações (3.3) podem ser resolvidas pelo método
scoring
de Fisher
(McCullagh e Nelder, 1989); um procedimento equivalente e também conveniente para
a resolução destas equações é o procedimento de mínimos quadrados reponderados
iterativamente (MQRI). Dado um vetor inicial
β
ββ
β
(0)
, calcula-se a resposta modificada
)(m
i
z
(3.4) e os pesos
)(m
i
w
(3.5 ) no passo
m
.
)(
)()(
)(
m
i
i
m
ii
m
i
m
i
yz
∂
∂
−+=
µ
η
µη
(3.4)
1)(
2
)(
)(
−
∂
∂
=
m
i
i
i
m
i
Vw
η
µ
,
n
i
,...,
1
=
(3.5)
19
onde
∑
=
−−
+=
d
j
ij
m
j
mm
i
x
1
)1()1()(
βαη
,
)(
)(1)( m
i
m
i
g
η
µ
−
=
e
)(
)(
m
i
i
m
i
V
∂
∂
=
Θ
µ
.
O vetor
)()(1)()(
')'(
mmmm
zWXXWX
−
=
β
ββ
β
,
...,
2
,
1
=
m
é obtido da regressão de
)(m
i
z
em
i
x
com peso
)(
m
i
w
,
n
i
,...,
1
=
. A matriz
X
é a matriz de planejamento de dimensão
)
1
(
+
×
v
n
, cuja i-ésima linha corresponde ao vetor
),1(
i
x
,
),...,(
)()(
1
)(
m
n
mm
zzz
=
e
=
)(m
W
diagonal
},...,{
)()(
1
m
n
m
ww
. Considerando agora o novo vetor
β
ββ
β
, o critério de
parada no
scoring
de Fisher se baseia numa medida de proximidade das estimativas;
assim o processo é repetido até que
δ
β
ββ
≤
−
∑
∑
=
−
=
−
d
j
m
j
d
j
m
j
m
j
1
)1(
1
)1()(
||||
||||
, (3.6)
para um valor
0
>
δ
pré-estabelecido.
3.2 Modelos aditivos generalizados
Na Seção 3.1 foi visto que o preditor linear é uma função linear de cada uma das
variáveis preditoras
d
XX ,...,
1
. No entanto, uma relação menos rígida pode ser adotada
substituindo o termo linear correspondente a cada covariável por uma função não
especificada dessa variável, obtendo-se o preditor aditivo
)(...)()(
1
d
XfXfg
+
+
+
=
=
α
µ
η
. (3.7)
A classe de modelos assim obtida é denominada Modelos Aditivos Generalizados e pode
ser vista como uma generalização dos MLG’s.
O preditor (3.7) corresponde a um modelo totalmente não paramétrico. Porém,
também fazem parte dos MAG’s os modelos semiparamétricos, cujo preditor combina
20
formas paramétricas de algumas das
r
variáveis preditoras com termos não
paramétricos das outras
)
(
r
d
−
variáveis. Nestes casos o preditor pode ser escrito
como
)(...)(...)(
111 drrr
XfXfXXg
+
+
+
+
+
+
=
=
+
β
β
α
µ
η
. (3.8)
Considerando
n
observações
),,...,,(
21 idiii
YXXX
do modelo Poisson
semiparamétrico
)(~
ii
PoissonY
µ
n
i
,...,
1
=
, o preditor
η
assume a forma
).(...)(...)log(
)1(11 diirririii
XfXfXX
+
+
+
+
+
+
=
=
+
β
β
α
µ
η
(3.9)
A estimação dos MAG’s e testes de hipóteses sobre os componentes do modelo
foram desenvolvidos em analogia a procedimentos utilizados com esses objetivos nos
MLG’s, modificando-os de forma que as funções
f
em (3.7) sejam estimadas por meio
da utilização de suavizadores.
O processo de ajuste dos MAG’s baseia-se na combinação de dois procedimentos
iterativos: o procedimento de ponderação local – PPL – (
local scoring
) e o retroajuste
(
backfitting
).
O PPL é um procedimento similar ao procedimento MQIR utilizado no ajuste dos
MLG’s e corresponde a um ciclo externo no processo e estimação necessário para o
ajuste de um modelo com estrutura semelhante a um MLG; o retroajuste, um ciclo
interno ao PPL, é o algoritmo responsável pela estimação de cada função
f
por meio
da utilização de suavizadores ponderados.
3.2.1 Ajuste dos modelos não paramétricos
O ajuste de um MAG pode ser efetuado nas três etapas do algoritmo PPL
esquematizado a seguir. O passo 2 é o retroajuste. Sejam
)(m
α
o valor estimado de
α
e
)(m
j
f
a estimativa de
j
f
,
d
j
,...,
1
=
, no passo
m
do procedimento iterativo, com
0
=
m
denotando o passo inicial.
21
Dados valores iniciais para
=
∑
=
n
i
i
n
y
g
1
)0(
α
e
0...
)0(
)0(
1
===
d
ff
, os valores da
variável modificada
)(m
i
z
e dos pesos
)(m
i
w
são calculados da mesma forma que em
(3.4) e (3.5), porém, neste caso,
∑
=
−−
+=
d
j
ij
m
j
mm
i
xf
1
)1()1()(
αη
. O algoritmo consiste em
iterar as seguintes etapas para
,...
2
,
1
=
m
Passo 1 – Para
n
i
,...,
1
=
calcular a variável resposta modificada e os pesos:
)(
)()(
)(
m
i
i
m
ii
m
i
m
i
yz
∂
∂
−+=
µ
η
µη
(3.10)
1)(
2
)(
)(
−
∂
∂
=
m
i
i
i
m
i
Vw
η
µ
, (3.11)
onde
∑
=
−−
+=
d
j
ij
m
j
mm
i
xf
1
)1()1()(
αη
,
)(
)(1)( m
i
m
i
g
η
µ
−
=
e
)(
)(
m
i
i
m
i
V
∂
∂
=
Θ
µ
.
Passo 2 – Os valores
)()(
1
,...,
m
d
m
ff
são obtidos com o uso do algoritmo de
retroajuste:
Inicializa-se
)(m
α
com a média amostral da variável modificada no passo
m
,
∑
=
==
n
i
m
i
mm
n
z
z
1
)(
)()(
α
e
)1(0 −
=
m
jj
ff
,
d
j
,...,
1
=
e calcula-se
)(
)(
)()(
)(
m
vj
m
j
m
vj
rSf =
, (3.12)
até que
ε
≤−
−
||||
)(
)1(
)(
)(
m
vj
m
vj
ff
, para um valor
0
>
ε
pré-estabelecido. Em (3.12),
),...,(
)()()(
m
vij
m
vij
m
vj
rrr =
é o vetor de resíduos parciais com elementos dados por
∑∑
+=
−
−
=
−−−=
d
jk
ik
m
vk
j
k
ik
m
vk
mm
i
m
ij
xfxfzzr
1
)(
)1(
1
1
)(
)(
)()()(
)()(
,
22
para
,...
2
,
1
=
v
, e
)(m
j
S
, de dimensão
n
n
×
, é a matriz suavizadora ponderada relativa
à j-ésima covariável. No caso do suavizador
loess
, a i-ésima linha de
)(m
j
S
corresponde
à i-ésima linha da matriz suavizadora ponderada
iii
xx
m
x
AXXAXXS ')'(
1)( −
=
, onde
=
X
),1(
j
X
,
1
é um vetor (
1
×
n
) de valores unitários, e
i
x
A
=diagonal
},...,{
)(
,
)(
11,
m
nnx
m
x
wuwu
ii
com
jx
i
u
,
e
)(m
j
w
definidos, respectivamente, em
(2.19) e (3.11).
Passo 3 – Os passos 2 e 3 são repetidos até que
δ
≤
−
∑
∑
=
−
=
−
d
j
m
j
d
j
m
j
m
j
f
ff
1
)1(
1
)1()(
||||
||||
,
para um valor
0
>
δ
pré-estabelecido.
O algoritmo de retroajuste corresponde ao método de Gauss-Seidel para resolver
o sistema de equações lineares:
=
)()(
)()(
2
)()(
1
)(
)(
2
)(
1
)()(
)(
2
)(
2
)(
1
)(
1
)(
1
)(
2
mm
d
mm
mm
m
d
m
m
m
d
m
d
mm
mmm
d
m
zS
zS
zS
f
f
f
ISS
SSI
SSS
S
S
I
MM
L
MOMM
L
L
M
. (3.13)
O retroajuste é um método eficiente para resolvê-lo, principalmente quando o
número de parâmetros é grande. Motivações para a utilização dessas equações na
obtenção de
)()(
1
,...,
m
d
m
ff
podem ser encontradas em Hastie e Tibshirani (1990).
Se a variável resposta segue uma distribuição Normal, a função de ligação é a
identidade, então
Y
Z
=
,
I
W
=
e o procedimento MQIR é substituído por um método
direto, ou seja, apenas o ciclo interno, correspondente ao retroajuste, é necessário. No
caso de um MAG com apenas uma função não especificada, isto é,
1
=
d
em (3.7), o
23
algoritmo de retroajuste não é necessário, pois
)(m
f
pode ser obtido diretamente com
a utilização de um alisador ponderado aplicado aos resíduos aplicado aos resíduos
)()()( mm
i
m
i
zzr
−
=
em função de
i
x
,
n
i
,...,
1
=
com matriz de pesos
)(m
W
.
Embora o retroajuste seja um algoritmo eficiente para resolver (3.13), pelo
menos conceitualmente, a solução direta pode ser utilizada. Estimativas para
d
ff
,...,
1
podem ser obtidas pela relação
CzMf
1
ˆ
−
=
com
=f
ˆ
d
f
f
f
ˆ
ˆ
ˆ
2
1
M
,
=
M
ISS
SSI
SSS
S
S
I
dd
d
L
MOMM
L
L
M
22
111
2
e
=
C
d
S
S
S
M
2
1
,
se a inversa de
M
existe.
Em particular, escrevendo
CMER
jj
1
−
=
,
d
j
,...,
1
=
,
onde
j
E
denota uma matriz de dimensão (
nd
n
×
) composta por
d
”blocos” de
dimensão (
n
n
×
), com todos os blocos nulos à exceção do j-ésimo, que é uma matriz
identidade de tal maneira que
zRf
j
=
ˆ
.
Seja
d
fff
+
+
=
...
1
, então
zRzRzRf
NDd
=++= ...
ˆ
1
,
24
onde
dND
RRR
+
+
=
...
1
é a matriz suavizadora ponderada que produz
f
ˆ
a partir de
z
.
Para modelos que envolvem apenas duas matrizes suavizadoras em seu ajuste,
2
=
d
, Hastie e Tibshirani (1990) fornecem expressões mais simples para
1
R
e
2
R
:
)()(
1
1
211
SISSIIR
−
−
−
=
−
(3.14)
)()(
2
1
122
SISSIIR
−
−
−
=
−
.
Nesta caso,
)())(()(
1
1
21221
SISSISIIRRR
ND
−
−
−
−
=
+
=
−
.
Expressões recursivas para modelos envolvendo mais de dois suavizadores foram
deduzidas por Opsomer (2000). O custo computacional para obter
j
R
,
d
j
,...,
1
=
, a
partir destas expressões é, entretanto, elevado.
O algoritmo do retroajuste é mais eficiente do ponto de vista computacional, mas
a solução direta pode ser utilizada como uma alternativa para obter expressões para
j
f
ˆ
e
η
ˆ
que tornem mais simples o estudo de suas propriedades estatísticas.
A convergência do procedimento de ajuste dos MAG’s está condicionada à
convergência do retroajuste, uma vez que o PPL não apresenta, em geral, problemas
dessa ordem (Hastie e Tibshirani, 1990). Resultados sobre a convergência desse
procedimento podem ser encontrados em Buja
et al.
. (1989) e Opsomer (2000).
3.2.2 Ajuste dos modelos semiparamétricos
Considere o modelo semiparamétrico
25
)(...)(...)(
111 drrr
xfxfXXg
+
+
+
+
+
+
=
=
+
β
β
α
η
µ
. (3.15)
Os parâmetros
r
β
β
α
,...,,
1
e as funções
dr
ff ,...,
1+
podem ser estimados através
do PPL e do retroajuste. Dado os valores iniciais
)0(
β
ββ
β
=
)',...,,(
)0()0(
1
)0(
r
β
β
α
e
)0()0(
1
,...,
dr
ff
+
estimativas para
β
ββ
β
e
dr
ff ,...,
1+
são obtidas resolvendo-se, iterativamente,
as seguintes equações
−=
∑
+=
−
d
rj
m
j
mmmm
fzWXXWX
1
)()()(1)()(
')'(
β
(3.16)
e
−−=
∑
≠
+=
d
ji
ri
m
i
mmm
j
m
j
fXzSf
1
)()()()()(
β
,
d
r
j
,...,
1
+
=
(3.17)
onde
),...,,1(
1 r
XXX
=
é a matriz de especificação correspondente aos termos
paramétricos do modelo com
j
X
,
r
j
,...,
1
=
, denotando o vetor dos valores
observados da j-ésima covariável e
)(m
j
S
,
d
r
j
,...,
1
+
=
, a matriz suavizadora
ponderada relativa à j-ésima covariável no m-ésimo passo do procedimento iterativo
PPL. Após obter as estimativas
)(m
β
ββ
β
e
)(m
j
f
,
d
r
j
,...,
1
+
=
, pelo retroajuste, valores de
)1(
+
m
η
,
)1(
+
m
µ
,
)1( +m
z
e
)1( +m
w
são calculados pelo PPL e o processo é repetido até a
convergência.
Quando existe apenas uma função não especificada, isto é,
1
+
=
r
d
em (3.8)
as expressões (3.16) e (3.17) se reduzem a
(
)
)()()(1)()(
')'(
mmmmm
fzWXXWX
−=
−
β
(3.18)
e
(
)
)()()()( mmmm
XzSf
β
−=
. (3.19)
Substituindo-se (3.19) em (3.18) obtém-se
26
(
)
)()()(1)()()(
'])('[
mmmmmm
zSIWXXSIWX −−=
−
β
(3.20)
e dessa forma o retroajuste pode ser evitado. Após a obtenção de uma estimativa
)(m
β
segundo (3.20), uma estimativa
)(m
f
é calculada segundo (3.19). Por intermédio do
PPL estimam-se novos valores
)1( +m
η
,
)1( +m
µ
,
)1( +m
z
e
)1( +m
w
e o processo é repetido até
a convergência.
Um fato pouco evidenciado na literatura é que, em geral, os estimadores dos
modelos semiparamétricos não são identificáveis quando incluem o intercepto
(Opsomer e Ruppert, 1999). Neste caso, quando a soma dos elementos das linhas de
S
é igual a 1, o que acontece no caso
loess
,
XSIWX
mm
)('
)()(
−
e
XRIWX
m
ND
m
)('
)()(
−
são singulares e uma solução simples para esse problema é
substituir as matrizes
)(m
j
S
por matrizes centradas, da forma
)(
)/'11(
m
j
SnI −
; este
procedimento faz com que a média de
f
ˆ
seja igual a zero em cada passo e o modelo
torne-se identificável.
3.2.3 Função desvio
Nos MLG’s o desvio para o modelo ajustado
µ
ˆ
,
)
ˆ
;(
µ
yD
e o desvio parcial para
dois modelos ajustados
1
ˆ
µ
e
2
ˆ
µ
,
)
ˆ
;
ˆ
(
21
µ
µ
D
, são quantidades bem conhecidas,
utilizadas, respectivamente, para avaliar a qualidade do ajuste e comparar dois modelos
ajustados. Sem perda de generalidade, pode-se escrever
)
ˆ
;(
η
yD
como sendo o desvio
para o modelo ajustado, uma vez que
µ
ˆ
está relacionado com
η
ˆ
por meio da função
de ligação
η
µ
=
)
(
g
.
No caso dos MAG’s, o desvio também pode ser usado como uma medida de
ajuste e a comparação entre modelos. No modelo Poisson, o desvio é dado por
∑
=
−−=
n
i
iiiii
yyyyD
1
)
ˆ
()
ˆ
/log((2)
ˆ
;(
µµµ
.
Embora as distribuições assintóticas dessas estatísticas não tenham sido
determinadas, Hastie e Tibshirani (1990) mostram, por simulações, que a distribuição
27
2
χ
é uma boa aproximação. Uma medida aproximada para os graus de liberdade de
)
ˆ
;(
η
yD
é dada por
∑
=
−−−=
d
j
j
Strngl
1
]1)([1
. (3.21)
3.2.4 Seleção do parâmetro de suavização
Um possível critério para selecionar parâmetros de suavização em um MAG com
p
termos não paramétricos é baseado na estatística
VC
∑
=
−
=
n
i
i
yD
n
VC
1
1
)
ˆ
;(
1
λ
µ
. (3.22)
A idéia é minimizar esta quantidade sob
d
λ
λ
,...,
1
, os parâmetros de suavização para
cada curva ajustada; no entanto o custo computacional é muito grande já que são
necessárias
n
aplicações completas do procedimento PPL para cada valor pré-fixado
dos parâmetros de suavização.
Uma opção para a seleção destes parâmetros é baseado na estatística
∑
=
+=
n
i
ii
Rtr
n
yD
n
AIC
1
)(
2
)
ˆ
;(
1
Φµ
(3.23)
inspirada no critério de informação de Akaike (ver Hastie e Tibshirani, 1990 e
Hastie,1992), valores pequenos desta estatística indicam um bom ajuste do modelo. O
valor
Φ
, definido como parâmetro de dispersão na equação (3.1), está associado à
distribuição de
Y
: se a variável resposta tem distribuição Normal, com variância
2
σ
,
2
1
σ
Φ
=
; se a variável resposta tem distribuição Poisson,
1
=
Φ
(veja resultados para
outras distribuições em Paula, 2000). Embora muito empregada na prática, não existem
resultados teóricos sobre a adequação de sua utilização como um critério para seleção
28
do parâmetro de suavização. O ganho computacional é dado pelo fato de que a
estatística
AIC
requer somente uma aplicação do PPL para cada valor
d
λ
λ
,...,
1
.
29
4. MODELOS AUTO-REGRESSIVOS GENERALIZADOS
Introduzido por Davis
et al.
(1999), o GLARMA –
autoregressive moving average
generalized linear models
– é um modelo utilizado em estudos de séries temporais
sendo capaz de capturar uma gama de estruturas de dependência nas observações.
Neste capítulo o GLARMA é estendido para uma classe de modelos auto-regressivos
aditivos generalizados para séries de contagem cuja distribuição condicional dada as
observações passadas e variáveis explicativas segue uma distribuição de Poisson.
4.1 Modelos Poisson auto-regressivo média móvel linear generalizados
(Poisson-GLARMA)
A classe GLARMA é uma classe de modelos que estende o processo ARMA (auto-
regressivo médias móveis) Gausssiano de séries temporais para um modelo mais
flexível para séries de contagem não-Gaussianas. A variável dependente é suposta ter
uma distribuição condicional na família exponencial dado todo o passado do processo.
Para introduzir o modelo Poisson-GLARMA, assuma que a observação
i
Y
dado o
passado histórico
1−i
F
tem distribuição Poisson com média
i
µ
,
)(~|
1 iii
PoissonFY
µ
−
,
n
i
,...,
1
=
.
A função de ligação
η
segue a forma
ididiii
TXX
+
+
+
+
=
=
β
β
α
µ
η
...)log(
11
(4.1)
onde
i
T
, responsável pela estrutura de correlação de
i
Y
, é dado por
∑
∞
=
−
=
1j
jiji
T
επ
, e
i
ε
assume a forma dos resíduos de Pearson
30
i
ii
i
Y
µ
µ
ε
−
=
. (4.2)
A estrutura média móvel infinita de
i
T
pode ser especificada em termos de um
número finito de parâmetros. Uma forma de parametrizar os pesos médias móveis
j
π
é
e expressá-los como coeficientes de um filtro auto-regressivo médias móveis (Box e
Jenkins, 1976)
1
)(
)(
)(
1
−==
∑
∞
=
B
B
BB
i
i
i
φ
θ
ππ
onde
p
p
BBBB
φφφφ
−−−−= ...1)(
2
2
1
1
q
q
BBBB
θθθθ
++++= ...1)(
2
2
1
1
são os polinômio auto-regressivo e polinômio média móvel com todas raízes fora do
círculo unitário. Dessa forma
i
T
pode ser expresso por
∑∑
=
−
=
−
++=
q
j
jiji
p
j
jiji
TT
11
)(
εθεφ
. (4.3)
A estimação dos parâmetros do GLARMA,
)'
'
,
'
(
ξ
β
ϕ
=
, onde
)'
'
,
'
(
θ
φ
ξ
=
, é feita
conjuntamente através da função de verossimilhança, maximizada pelo método
numérico Newton-Raphson (Davis
et al.
, 2003).
Considere
l
a densidade Poisson condicional de
i
Y
dado
1−
i
F
e defina
)|(log)(
1−
=
iii
FyL l
ϕ
. A log-verossimilhança pode ser escrita como
∑
=
n
i
i
L
1
)(
ϕ
que, ignorando termos que não envolvem os parâmetros, se torna
31
)
(
ϕ
L
∑
=
=
n
i
ii
Y
1
(
η
−
)
(
ϕ
i
η
ε
)
)(
ϕ
(4.4)
onde
ii
η
µ
=
)log(
)
(
ϕ
∑
∞
=
++++=
1
11
(...
j
ididi
XX
πββα
)
ξ
ji−
ε
)
(
ϕ
e
i
ε
=
)
(
ϕ
iii
Y
µµ
)( −
.
Para facilitar a compreensão dos cálculos, a dependência de
i
ε
em
ϕ
foi
desconsiderada. A primeira e a segunda derivadas de
L
são dadas pelas expressões
(4.5) e (4.6)
∑∑
==
∂
∂
=
∂
∂
−=
∂
∂
n
i
i
ii
n
i
i
ii
Y
L
11
)(
ϕ
η
µε
ϕ
η
µ
ϕ
(4.5)
∑
=
∂
∂
∂
∂
−
∂∂
∂
−=
∂∂
∂
n
i
ii
i
i
ii
Y
L
1
2
2
''
)(
'
ϕ
η
ϕ
η
µ
ϕϕ
η
µ
ϕϕ
(4.6)
∑
=
∂
∂
∂
∂
−
∂∂
∂
=
n
i
ii
i
i
ii
1
2
''
ϕ
η
ϕ
η
µ
ϕϕ
η
µε
.
Maiores detalhes sobre expressões úteis no cálculo dessas derivadas e
resultados assintóticos das estimativas dos parâmetros podem ser encontrados em
Davis
et al.
(2003).
Para inicializar o método recursivo de Newton Raphson na maximização
numérica da log-verossimilhança
)(
ϕ
i
L
, Davis
et al.
(2003) sugerem que os valores
obtidos das estimativas do GLARMA sem os termos auto-regressivos média móveis
sejam utilizados como valores iniciais. A convergência, na maioria dos casos, ocorre
após 10 iterações. A matriz de covariância dos estimadores é estimada por
32
1
2
'
)
ˆ
(
ˆ
−
∂∂
∂
−=Ω
ϕϕ
θ
L
. 4.7
Maiores detalhes sobre as condições de estacionaridade, propriedades,
estimação e inferência dos modelos GLARMA podem ser vistos em Benjamin
et al.
(2003) e Drescher (2005).
4.2 Modelos Poisson auto-regressivos aditivos generalizados
(Poisson MAG-AR)
Considere que
i
Y
dado o passado histórico
1−i
F
tem distribuição Poisson com
média
i
µ
como na Seção 4.1. Utilizando a construção MAG semiparamétrico e
adicionando uma estrutura de correlação entre os dados proposta no modelo GLARMA,
a equação (4.1) pode ser reescrita como
idiirririii
TXfXfXX
+
+
+
+
+
+
+
=
=
+
)(...)(...)log(
)1(11
β
β
α
µ
η
, (4.8)
com
i
T
sendo um processo auto-regressivo de ordem
p
, AR(
p
),
i
p
j
jiji
TT
εφ
+=
∑
=
−
)(
1
,
i
ε
definido em (4.2) e
n
i
,...,
1
=
.
Muitos pesquisadores vêm utilizando o MAG na modelagem de séries temporais
ignorando a dependência entre as observações. O modelo MAG-AR além de contar com
a vantagem que o MAG oferece em eliminar a necessidade de especificar uma forma
paramétrica para a associação de algumas covariáveis com a variável dependente ainda
é capaz de capturar estruturas de correlação entre os dados.
A estimação do MAG-AR foi desenvolvida em analogia aos procedimentos de
estimação do GLARMA, descritos na Seção 4.1, modificando-os de forma que as
33
Tempo
Y
0 20 40 60 80 100
0 5 10 15
Série Resposta
Tempo
Resíduos
0 20 40 60 80 100
-2 0 2 4 6 8
Resíduos
0 5 10 15 20
-0.2 0.0 0.2 0.4 0.6 0.8 1.0
Lag
ACF
Series Resíduos
5 10 15 20
-0.2 -0.1 0.0 0.1 0.2 0.3
Lag
Partial ACF
Series Resíduos
funções
f
em (4.8) sejam estimadas por meio da utilização de suavizadores como nos
procedimentos de estimação dos MAG’s, descritos na Seção 3.2.1.
Para verificar se o procedimento proposto realmente funciona em processos que
exibem uma estrutura auto-regressiva, um exercício de simulação simples foi
implementado para
i
T
seguindo um processo AR(1). Assim, foram geradas uma
variável resposta
Y
seguindo um processo de Poisson e duas séries temporais como
variáveis explicativas,
1
X
– relacionada linearmente com a variável resposta – e
2
X
–
que apresenta uma relação não-linear com
Y
. A Variável
1
X
é um modelo auto-
regressivo (AR) de ordem 1 da forma
i
i
i
eXX
+
=
− )1
(11
4,0
,
)1,0(~ Ne
i
e
i
viX
+
=
)12/2cos(
2
π
,
)01,0;0(~ Nv
i
,
100
,...,
1
=
i
. Considere
i
Y
gerado a partir de
uma Poisson com média
}exp{
211 iiii
Txx
+
+
+
=
β
α
µ
, com
α
fixado em 0,01,
1
β
fixado em 0,08 e
i
T
um AR de ordem 1 com
6
,
0
=
φ
. Os dados simulados foram
ajustados através do MAG semiparamétrico, ignorando a correlação existente nas
observações, e também através do MAG-AR(1), com o objetivo de se comparar as duas
modelagens.
A Figura 4.1 apresenta a série
i
Y
gerada com a estrutura descrita acima, assim
Figura 4.1 – Avaliação residual da série
Y
modelada via MAG
34
como gráficos de resíduos obtidos através do ajuste de um MAG a esta série. Dos
gráficos de autocorrelação (ACF) e autocorrelação parcial (
partial
ACF) é clara a
presença de um AR(1) nos resíduos, estrutura não capturada pelo modelo.
Para corrigir este problema, os dados são então ajustados através do MAG-AR(1)
– os gráficos de resíduos são apresentados na Figura 4.2. Quando a nova modelagem é
aplicada, a estrutura auto-regressiva é “captada” pelo modelo e os resíduos se tornam
um ruído branco, logo parece que o modelo proposto é capaz de incorporar a estrutura
auto-regressiva, fazendo com que se obtenha um melhor ajuste para as séries
envolvidas.
Tempo
Y
0 20 40 60 80 100
0 5 10 15
Série Resposta
Tempo
Resíduos
0 20 40 60 80 100
-2 0 2 4 6
Resíduos
0 5 10 15 20
-0.2 0.0 0.2 0.4 0.6 0.8 1.0
Lag
ACF
Series Resíduos
5 10 15 20
-0.2 -0.1 0.0 0.1 0.2
Lag
Partial ACF
Series Resíduos
Figura 4.2 – Avaliação residual da série
Y
, modelada via MAG-AR
35
5. TÉCNICA BOOTSTRAP
Quando se deseja medir a precisão de estimadores para os parâmetros
desconhecidos de uma dada distribuição, geralmente calcula-se uma medida que
expresse a variabilidade dos mesmos. Mas se a distribuição exata do estimador é
desconhecida ou se o pesquisador tem acesso apenas à sua distribuição assintótica,
este cálculo pode ser complicado. Há mais de duas décadas surgiu um procedimento
computacional –
Bootstrap
– uma técnica de reamostragem que pode ser utilizada para
aproximar a distribuição teórica pela distribuição empírica de uma amostra finita de
observações (Efron, 1979). Porém, sendo um método numérico, a sua operacionalidade
somente se tornou viável com o advento dos computadores.
Em séries temporais, devido ao fato das observações serem correlacionadas, a
aplicação desta técnica requer vários cuidados e a reamostragem direta das
observações não pode ser feita. Nestes casos pode-se aplicar o
Bootstrap
amostrando
os resíduos diretamente de sua distribuição (
Bootstrap
paramétrico) ou reamostrando
os resíduos do modelo ajustado (
Bootstrap
não-paramétrico).
Em MAG’s e GLARMA’s o uso da técnica é ainda pouco aplicada e discutida. Em
estudo recente, Härdle
et al.
(2004) mostram como o procedimento pode ser utilizado
na correção do vício das estimativas paramétricas e não-paramétricas dos MAG’s, em
testes de hipótese e na construção de bandas de confiança.
5.1
Bootstrap
nas observações
A proposta original da técnica
bootstrap
é a reamostragem direta, com
reposição, das observações. Para ilustrar a técnica na modelagem aditiva generalizada,
considere
i
Y
,
di
i
XX ,...
1
,
n
i
,...,
1
=
, vetores de dados independentes. Considere ainda
que
i
Y
tenha distribuição Poisson e que sua média possa ser modelada por um MAG
semiparamétrico através da relação
)(...)(...)log()()]|([
)1(11 diirririi
XfXfXXgXYEg
+
+
+
+
+
+
=
=
=
+
β
β
α
µ
µ
,
36
n
i
,...,
1
=
. (5.1)
A técnica
bootstrap
consiste em reamostrar os pontos (
diii
xxy ,...,,
1
), com reposição,
obtendo vetores bootstrap
*
i
Y
,
**
1
,...,
dii
XX
.
5.2
Bootstrap
não paramétrico nos resíduos
Com o objetivo de não reamostrar diretamente uma série temporal devido a não
independência das observações, uma das alternativas é reamostrar os resíduos
utilizando o método
bootstrap
(Efron e Tibshirani, 1993): inicialmente ajusta-se um
modelo aos dados e reamostra-se os resíduos (que devem ser independentes e
identicamente distribuídos).
A abordagem não-paramétrica é assim classificada por não utilizar nenhuma
suposição quanto à distribuição dos resíduos ao reamostrá-los; neste caso usa-se uma
distribuição empírica. Os procedimentos
bootstrap
não-paramétrico para o MAG e para
o MAG-AR utilizando a distribuição de Poisson são descritas a seguir.
Seja
i
Y
,
n
i
,...,
1
=
, um vetor de observações independentes com distribuição
Poisson e esperança modelada por um MAG,
)(...)(...)log()()]|([
)1(11 diirririi
XfXfXXgXYEg
+
+
+
+
+
+
=
=
=
+
β
β
α
µ
µ
,
n
i
,...,
1
=
(5.2)
onde
d
XX ,...,
1
são covariáveis também independentes.
Após estimar os parâmetros
β
’s e as funções arbitrárias
f
’s os resíduos de
Pearson podem ser obtidos através da expressão
i
ii
i
Y
µ
µ
ε
ˆ
ˆ
−
=
,
n
i
,...,
1
=
, (5.3)
37
onde
)}(
ˆ
...)(
ˆ
ˆ
...
ˆ
ˆ
exp{
ˆ
)1(11 diirririi
xfxfxx ++++++=
+
ββαµ
. Em seguida reamostra-
se, com reposição,
i
ε
atribuindo-se a cada um uma massa de probabilidade igual a
n/1
. Dessa forma obtêm-se os resíduos
bootstrap
*
i
ε
. A partir daí, conforme expressão
(5.4), constrói-se, recursivamente, a série
bootstrap
*
i
Y
iiii
Y
µµε
ˆˆ
**
+=
,
n
i
,...,
1
=
, (5.4)
Considerando
i
Y
,
dii
XX
,...,
1
,
n
i
,...,
1
=
, séries temporais
com distribuição
Poisson e esperança modelada por um MAG-AR,
idiirririi
TXfXfXXgXYEg
+
+
+
+
+
+
+
=
=
=
+
)(...)(...)log()()]|([
)1(11
β
β
α
µ
µ
n
i
,...,
1
=
(5.5)
e os resíduos de Pearson definidos como em (5.3) com
})(
ˆ
...)(
ˆ
ˆ
...
ˆ
ˆ
exp{
ˆ
)1(11 idiirririi
Txfxfxx +++++++=
+
ββαµ
reamostra-se, com
reposição,
i
ε
, cada um com a mesma massa de probabilidade, obtendo os resíduos
bootstrap
*
i
ε
do MAG-AR. A série
bootstrap
*
i
Y
é construída recursivamente da mesma
forma que em (5.4).
5.3
Bootstrap
condicional
O
bootstrap
condicional, sugerido por Figueiras
et al.
(2005), é um método que
considera dados do tipo (
ii
yx ,
),
n
i
,...,
1
=
ou, de forma mais geral,
,...),,(
21
iii
xxy
,
assumindo que a distribuição de
i
Y
é conhecida e que seus valores são condicionais aos
valores
,...),(
21
ii
xx
. Os vetores das variáveis
Y
e
j
X
,
d
j
,...,
1
=
, podem ser dados
correlacionados ou não.
38
Para ilustrar a técnica na modelagem MAG, suponha que
PoissonY
i
~
com
esperança
i
µ
e considere
d
XX ,...,
1
relacionadas à
i
Y
pela expressão
)}(
ˆ
...)(
ˆ
ˆ
...
ˆ
ˆ
exp{
ˆ
)1(11 diirririi
xfxfxx ++++++=
+
ββαµ
.
A técnica
bootstrap
condicional consiste em gerar um
*
i
y
,
*
i
Y
com distribuição
Poisson de média
)}(
ˆ
...)(
ˆ
ˆ
...
ˆ
ˆ
exp{
ˆ
)1(11 diirririi
xfxfxx ++++++=
+
ββαµ
,
para cada
ponto
),...,(
1 dii
xx
. O vetor
*
i
Y
é a série
bootstrap
.
Considerando
d
XX ,...,
1
, séries temporais, relacionadas à série
PoissonY
i
~
pela expressão
})(
ˆ
...)(
ˆ
ˆ
...
ˆ
ˆ
exp{
ˆ
)1(11 idiirririi
Txfxfxx +++++++=
+
ββαµ
a técnica
segue a mesma forma.
5.4 Intervalos de confiança
bootstrap
Existem vários métodos para calcular intervalos de confiança para um
parâmetro, mas devido a algumas restrições tal como a imprecisão causada por
aproximações feitas através da distribuição assintótica, Efron & Tibshirani (1986)
propuseram métodos que fazem uso da técnica
bootstrap
na construção destes
intervalos. Neste capítulo são apresentadas
duas técnicas
bootstrap
de construção de
intervalos –
bootstrap
percentílico e
bootstrap
com correção do vício – utilizadas na
inferência da parte linear dos modelos MAG e MAG-AR.
5.4.1 Intervalos de confiança
bootstrap
percentílico
O intervalo percentílico
)%
1
(
γ
−
para o parâmetro
β
é definido pela expressão
(
)
*
)2/1(
*
)2/(
ˆ
;
ˆ
γγ
ββ
−
(5.5)
onde, por definição,
*
)(
ˆ
a
β
é o
)
100
(
a
⋅
-ésimo percentil empírico da distribuição
bootstrap (Efron e Tibshirani, 1993). Na prática, se são geradas
B
amostras
bootstrap
independentes,
B
xxx
*2*1*
,...,,
, e estima-se
*
ˆ
β
para cada uma delas tem-se que
*
)2/(
ˆ
γ
β
39
é o
)
2
/
(
γ
B
º valor ordenado das replicações
*
ˆ
β
; a mesma interpretação é dada para
*
)2/1(
ˆ
γ
β
−
.
5.4.2 Intervalos de confiança
bootstrap
com correção do vício
Um dos principais objetivos da teoria
bootstrap
é produzir intervalos de confiança
que realmente fornecem coberturas probabilísticas confiáveis para o parâmetro de
interesse. O método
bootstrap
com correção do vício, apesar de ser um método mais
complicado de se definir, é tão simples de ser usado quanto o método percentílico.
Além disso, os intervalos de confiança correção do vício levam em consideração o vício
do parâmetro estimado (Efron, 1982).
Os limites do intervalo de confiança correção do vício são encontrados através de
percentis empíricos da distribuição
bootstrap
, mas não são necessariamente
2
/
γ
e
2
/
1
γ
−
tal como no método percentílico. Os percentis usados para o cálculo dos
limites inferiores e superiores dos intervalos de confiança correção do vício dependem
do número
0
k chamado
Bias-corrected
, ou corretor do vício, que é definido pela
expressão
≤=
∑
=
−
)
ˆˆ
(
1
1
*1
0
B
t
t
I
B
k
ββ
Φ
ΦΦ
Φ
,
em que
I
é uma função indicadora que recebe valor 1 se
)
ˆ
ˆ
(
*
ββ
≤
t
e valor 0 caso
contrário;
B
é o número de amostras
bootstrap
independentes;
β
ˆ
é a estimativa do
parâmetro para os dados observados;
*
ˆ
β
é a estimativa do parâmetro para cada uma
amostra
bootstrap
e )(
⋅
Φ
é a função de distribuição da Normal padronizada.
Os limites do intervalo
bootstrap
com correção do vício são dados por
(
)
*
)(
*
)(
21
ˆ
,
ˆ
pp
ββ
40
sendo
+=
2
01
2
γ
zkp Φ
ΦΦ
Φ
e
+=
−
2
1
02
2
γ
zkp Φ
ΦΦ
Φ
– em que
x
z é o )100( x
⋅
-ésimo ponto
percentil da distribuição Normal padronizada – e
*
)(
1
ˆ
p
β
igual ao
)(
1
pB
O
valor ordenado
das replicações
*
ˆ
β
.
Caso a distribuição do estimador
*
ˆ
a
β
seja simétrica, tem-se que 0
0
=k ,
2/
1
γ
=
p
e
2/1
2
γ
−
=
p
. Logo, nesta situação o intervalo de confiança obtido é o
intervalo percentílico.
41
6. ANÁLISE DOS DADOS SIMULADOS
O desempenho da classe semiparamétrica dos MAG’s no ajuste de dados
independentes e no ajuste de séries temporais bem como a performance dos MAG-AR’s
no ajuste de séries temporais foi averiguada segundo a estimação do parâmetro linear
1
β
, via simulação de dados com variável resposta inteira e não negativa modelada pela
distribuição Poisson
)(~
ii
poissonY
µ
.
Para isto, foram consideradas duas covariáveis, sendo a primeira,
1
X
,
relacionada linearmente com a variável resposta e a segunda,
2
X
, relacionada a
Y
de
uma forma não linear através de uma função desconhecida. Dessa forma, a média da
variável
Y
– vetor de dados independentes ou de séries temporais – quando modelada
por um MAG é expressa pela relação
)()log(
211 i
i
ii
XfX
+
+
=
=
β
α
µ
η
,
n
i
,...,
1
=
. (6.1)
Considerando que a média da variável
Y
– agora uma série temporal – seja modelada
por um MAG-AR o preditor aditivo assume a forma
iiiiii
TXfXg
+
+
+
=
=
=
)()log()(
211
β
α
µ
µ
η
,
n
i
,...,
1
=
(6.2)
onde
i
T é um AR de ordem 1 da forma,
iii
TT
ε
φ
+
=
−1
.
As estatísticas que viabilizaram a comparação das performances do MAG, do
MAG-AR e dos procedimentos
bootstrap
foram o vício e o erro quadrático médio (EQM)
das estimativas. Os intervalos de confiança foram comparados via tamanho e
probabilidade de cobertura – expressa pela razão entre o número de intervalos que
contem o valor verdadeiro do parâmetro e o número total de intervalos construídos –
com o nível nominal fixado em 95%, isto é,
05
,
0
=
γ
. Para fins de comparação o
42
intervalo assintótico também foi construído, utilizando-se a distribuição assintótica
normal para o estimador de
1
β
.
Dois valores foram definidos para
n
: 100 e 500. Para cada um deles, o número
de simulações Monte Carlo (MC) e de replicações
bootstrap
foi fixado em 500. A técnica
loess
de suavização foi utilizada na estimação de
f
; o valor do parâmetro de
suavização foi escolhido com base na estatística AIC dado os valores pré-fixados 0,5;
0,6; 0,7 e 0,8.
Os dados foram simulados e modelados através da linguagem de programação
do
software
R; os algoritmos da modelagem MAG já estão implementados no software.
6.1 Resultados das simulações de dados independentes
O vetor resposta
i
Y
,
n
i
,...,
1
=
, foi gerado a partir de uma Poisson com média
}exp{
211 iiii
xx
τ
β
α
µ
+
+
+
=
,
onde
α
e
1
β
foram fixados, respectivamente, em 0,08 e 0,02. A covariável
1
X
foi
gerada a partir da distribuição Normal com média 3 e variância 4 e a variável
2
X
gerada a partir da equação
)(001.0
3
2
aaX
+
×
=
, com
a
sendo uma variável
aleatória de distribuição
)
4
,
10
(
N
. O ruído
τ
,
)
1
,
0
(
~
Normal
τ
, foi incluído no modelo
para introduzir aleatoriedade no cálculo de
Y
.
As Tabelas 1 e 2 mostram os resultados sobre a estimação pontual e intervalar,
respectivamente, do parâmetro
1
β
obtidos através da modelagem MAG. Da Tabela 1,
nota-se que a média das estimativas pontuais obtidas no Monte Carlo (MC)
superestimaram o verdadeiro valor do parâmetro e, apesar do vício ser praticamente o
mesmo tanto para as simulações de tamanho 100 quanto para as simulações de
tamanho 500, os valores do EQM neste último cenário têm uma ordem de grandeza
menor (da ordem de
4
10
−
para
100
=
n
e
5
10
−
para
500
=
n
.
Os resultados do procedimento
bootstrap
, estão bem próximos à média das
estimativas obtidas no MC, sendo que o melhor desempenho corresponde ao
bootstrap
nas observações. O
bootstrap
nos resíduos apresenta um EQM muito maior que os
43
outros dois tipos do
bootstrap
, do que se pode concluir que, se os dados são
independentes, o melhor é utilizar o
bootstrap
direto nas observações, pois como o
bootstrap
nos resíduos é dependente do modelo ajustado, ele pode carregar o vício das
estimativas calculadas para os parâmetros.
Tabela 1. Médias das estimativas, médias de vício e EQM na estimação de
1
β
no MAG para dados
independentes
n
MC
Bootstrap nas
observações
Bootstrap
condicional
Bootstrap nos
resíduos
100 0,0217 0,0211 0,0220 0,0219
0,0017 0,0011 0,0020 0,0019
(5,81e-04) (5,94e-04) (5,87e-04) (7,70e-04)
500 0,0216 0,0218 0,0219 0,0219
0,0016 0,0018 0,0019 0,0019
(3,45e-05) (3,82e-05) (3,95e-05) (5,10e-05)
Nota: valores em negrito são as médias de vício, valores entre parênteses são os EQM’s.
Da Tabela 2, pode-se concluir que os intervalos de confiança apresentam, em
geral, probabilidade de cobertura bastante próxima ao valor nominal (0,95) e que os
intervalos
bootstrap
com correção do vício têm melhor desempenho se comparados aos
intervalos
bootstrap
percentílico, com a única exceção para o procedimento
bootstrap
nas observações. Em geral, os intervalos de confiança construídos através da técnica
bootstrap
apresentam-se mais próximos a 0,95 que o intervalo assintótico, tanto para
100
=
n
quanto para
500
=
n
.
Tabela 2. Intervalos de confiança para
1
β
no MAG para dados independentes
n
Assintótico
Bootstrap nas observações Bootstrap condicional Bootstrap nos resíduos
Percentílico BC Percentílico BC Percentílico BC
100
0,928 0,946 0,938 0,936 0,938 0,928 0,948
(0,090) (0,095) (0,094) (0,091) (0,091) (0,099) (0,097)
500
0,932 0,944 0,942 0,940 0,942 0,934 0,950
(0,094) (0,097) (0,094) (0,095) (0,091) (0,107) (0,101)
Nota: valores em negrito são as probabilidades de cobertura, valores entre parênteses são as amplitudes dos
intervalos.
6.2 Resultados das simulações de séries temporais
Para os dados que representam séries temporais, a variável resposta
i
Y
,
n
i
,...,
1
=
, foi gerada a partir de uma Poisson com média
44
}exp{
211 iiii
Txx
+
+
+
=
β
α
µ
, (6.3)
onde
α
foi fixado em 0,08 e
1
β
em 0,02. A covariável
1
X
é um modelo auto-
regressivo da forma
iii
eXX
+
=
−1
4,0
,
)1,0(~ Ne
i
,
i
viX
+
=
)12/2cos(
2
π
,
)01,0;0(~ Nv
i
e
i
Z
um modelo auto-regressivo de ordem 1,
iii
TT
ε
φ
+
=
−1
, com
4
,
0
=
φ
e
6
,
0
e
i
ε
definido em (4.2). O procedimento
bootstrap
nas observações,
conforme justificado no Capítulo 5, não é realizado para as séries temporais simuladas.
Nesta seção, o objetivo é avaliar as estimativas do parâmetro linear
1
β
quando
os dados são modelados através do MAG, ignorando a dependência entre as
observações das séries temporais, e compará-las com as estimativas encontradas na
modelagem MAG-AR, classe capaz de capturar estruturas de dependência entre os
dados.
As Tabelas 3 e 4 apresentam as estimativas pontuais das séries temporais
geradas, a primeira para o MAG e a segunda para o MAG-AR(1). Ao contrário das
estimativas apresentadas na Seção 6.1 de dados independentes, todas as estimativas
apresentadas nas Tabelas 3 e 4 subestimaram o verdadeiro valor do parâmetro linear,
1
β
.
Tabela 3. Médias das estimativas, médias de vício e EQM na estimação de
1
β
no MAG para séries
temporais
n
φ
MC Bootstrap condicional Bootstrap nos resíduos
100 0,4 -0,0043 -0,0043 -0,0028
-0,0243 -0,0243 -0,0228
(2,97e-02) (2,97e-02) (2,75e-02)
0,6 0,0013 0,0015 0,0014
-0,0187 -0,0185 -0,0186
(3,93e-02) (3,93e-02) (3,47e-02)
500 0,4 -0,0031 -0,0033 -0,0063
-0,0231 -0,0233 -0,0263
(2,32e-02) (2,32e-02) (1,95e-02)
0,6 0,0023 0,0024 0,0023
-0,0177 -0,0176 -0,0177
(3,44e-02) (3,45e-02) (2,77e-02)
Nota: valores em negrito são as médias de vício, valores entre parênteses são os EQM’s.
Comparando o ajuste das séries de contagem utilizando o MAG e o MAG-AR(1)
pode-se verificar que as estimativas apresentam maior vício e EQM para o MAG (Tabela
45
3) que para o MAG-AR(1) (Tabela 4), evidenciando a necessidade de se utilizar este
último modelo para dados de séries temporais. Das Tabelas 3 e 4 verifica-se também
que são mais viciadas e possuem EQM’s maiores os resultados dos dados gerados com
φ
igual a
4
,
0
.
Novamente, o desempenho dos procedimentos
bootstrap
é muito
semelhante ao do MC. Neste caso, percebe-se que o EQM das estimativas utilizando o
bootstrap
nos resultados é ligeiramente menor que a do
bootstrap
condicional, com o
vício das estimativas para o modelo MAG-AR(1) também menor para o
bootstrap
nos
resíduos.
Tabela 4. Médias das estimativas, médias de vício e EQM na estimação de
1
β
no MAG-AR(1) para
séries temporais
n
φ
MC Bootstrap condicional Bootstrap nos resíduos
100 0,4 0,0169 0,0165 0,017
-0,0031 -0,0035 -0,0030
(7,80e-03) (7,86e-03) (7,78e-03)
0,6 0,0179 0,0174 0,0176
-0,0021 -0,0026 -0,0024
(3,66e-03) (3,91e-03) (3,47e-03)
500 0,4 0,0178 0,0177 0,0181
-0,0022 -0,0023 -0,0019
(7,44e-03) (7,54e-03) (7,25e-03)
0,6 0,0186 0,0186 0,0189
-0,0014 -0,0014 -0,0011
(2,14e-03) (2,32e-03) (1,93e-03)
Nota: valores em negrito são as médias de vício, valores entre parênteses são os EQM’s.
Na Tabela 5 estão os resultados das estimações intervalares do parâmetro linear
1
β
das séries temporais de contagem modeladas através do MAG e na Tabela 6 os
resultados das séries com ajuste MAG-AR. Em todos os casos nota-se que a cobertura
dos intervalos está mais próxima ao nível de 95% no modelo MAG-AR(1), com
intervalos praticamente do mesmo tamanho.
Comparando os procedimentos
bootstrap
, o
bootstrap
nos resíduos apresenta
melhor desempenho quanto à cobertura, apesar da amplitude dos mesmos ser um
pouco maior. Os intervalos BC apresentam-se levemente superiores ao percentílico,
mas ambos têm melhor desempenho se comparados ao intervalo assintótico. Nestes
casos, os resultados foram melhores para
4
,
0
=
φ
.
46
Tabela 5. Estimação intervalar do parâmetro
1
β
e comparação dos procedimentos MC e
bootstrap
’s –
Séries temporais modeladas através do MAG
n
φ
Assintótico
Bootstrap condicional Bootstrap nos resíduos
Percentílico BC Percentílico BC
100
0,4
0,870 0,884 0,886 0,922 0,924
(0,115) (0,118) (0,118) (0,132) (0,133)
0,6
0,866 0,880 0,882 0,914 0,916
(0,113) (0,114) (0,116) (0,134) (0,134)
500
0,4
0,874 0,906 0,908 0,930 0,936
(0,117) (0,121) (0,122) (0,138) (0,138)
0,6
0,872 0,888 0,888 0,926 0,926
(0,117) (0,119) (0,119) (0,136) (0,137)
Nota: Valores em negrito são as probabilidades de cobertura, valores entre parênteses são as amplitudes dos
intervalos
Tabela 6. Estimação intervalar do parâmetro
1
β
e comparação dos procedimentos MC e
bootstrap
’s –
Séries temporais modeladas através do MAG-AR
n
φ
Assintótico
Bootstrap condicional Bootstrap nos resíduos
Percentílico BC Percentílico BC
100
0,4
0,898 0,902 0,914 0,938 0,938
(0,116) (0,118) (0,119) (0,135) (0,133)
0,6
0,884 0,888 0,902 0,938 0,940
(0,110) (0,109) (0,113) (0,134) (0,135)
500
0,4
0,901 0,922 0,926 0,942 0,940
(0,121) (0,119) (0,119) (0,137) (0,138)
0,6
0,898 0,916 0,920 0,940 0,938
(0,116) (0,118) (0,118) (0,134) (0,134)
Nota: Valores em negrito são as probabilidades de cobertura, valores entre parênteses são as amplitudes dos
intervalos
47
7. APLICAÇÃO A SÉRIES REAIS
Estudos epidemiológicos realizados em diferentes centros de pesquisa têm
detectado associações significativas entre morbi-mortalidade por causas respiratórias e
poluição atmosférica em populações urbanas – Schwartz (1994) e Saldiva
et al.
(1995),
por exemplo. A maior parte desses estudos é do tipo ecológico, isto é, de base
populacional e um grupo, ao invés de um indivíduo, constitui a unidade de observação
seguida ao longo do tempo. Consistem, em geral, da observação de eventos tais como
mortalidade, admissões hospitalares ou sintomas respiratórios. Esse tipo de
planejamento é menos suscetível a variáveis de confusão individuais como fumo,
pressão arterial e fatores sócio-econômicos (Rothman
et al.
, 1998), pois esses fatores
não variam de dia para dia com a poluição atmosférica.
A cidade de São Paulo, o segundo maior centro urbano da América Latina, possui
um cenário apropriado para o desenvolvimento de estudos dos efeitos da poluição
atmosférica sobre a saúde de seres humanos; um dos principais motivos é a grande
oferta de transporte coletivo e uma malha viária de uma frota de pouco mais de
6.100.000 veículos leves em toda região metropolitana (DETRAN-SP) que constituem a
principal fonte da poluição do ar. A poluição atmosférica na cidade é
predominantemente gerada por fontes poluidoras móveis, além disso, suas condições
geográficas e meteorológicas desfavorecem a dispersão dos poluentes, principalmente
durante os meses de inverno, quando, com freqüência, ocorrem inversões térmicas.
Neste Capítulo, os modelos MAG e MAG-AR(1) são ajustados a dados mensais do
número de pacientes internados com afecção das vias aéreas superiores – AVAS – na
cidade de São Paulo, nos anos de 1997 a 2000, para se verificar o desempenho da
modelagem proposta, MAG-AR(1), em relação aos MAG’s numa aplicação a dados reais,
segundo a estimação do parâmetro linear
1
β
. As covariáveis consideradas são as séries
dióxido de enxofre - SO
2
/(
µ
g/m
3
) - e monóxido de carbono - CO/ppm. A variável
tempo
i
(
n
i
,...,
1
=
) e as variáveis harmônicas
)2( nisen
π
e
)nicos(
π
2
foram
inseridas no modelo com o objetivo de ajustar o efeito da tendência e da sazonalidade
(Moretin
et al
., 2004).
48
MESES DE OBSERVAÇÃO
AVAS
50403020100
12
11
10
9
8
7
6
5
SÉRIE AVAS
MESES DE OBSERVAÇÃO
SO2
50403020100
30
25
20
15
10
SÉRIE SO2
Ano
SO2
2000199919981997
30
25
20
15
10
BOXPLOT SO2
Ano
CO
2000199919981997
7
6
5
4
3
2
BOXPLOT CO
ANO
AVAS
2000199919981997
12
11
10
9
8
7
6
5
BOXPLOT AVAS
7.1 Análise descritiva
A seguir são apresentadas tabelas e gráficos construídos com o objetivo de
resumir os dados mensais das variáveis de interesse: paciente com afecção das vias
aéreas superiores – AVAS – e os poluentes SO
2
e CO, nos anos de 1997 a 2000.
Figura 7.1 – Representação gráfica das séries AVAS, SO
2
e CO nos anos de 1997 a 2000.
Figura 7.2 – Gráficos do tipo
Box-plot
para AVAS, SO
2
e CO nos anos de 1997 a 2000.
MESES DE OBSERVAÇÃO
CO
50403020100
7
6
5
4
3
2
SÉRIE CO
49
A Figura 7.1 apresenta as séries cronológicas e a Figura A.1, no Apêndice A,
apresenta gráficos do tipo
box-plot
. Na Figura A.1 não se nota valores discrepantes da
AVAS nos anos avaliados. Na Tabela 7, observa-se que o número mensal médio de
pacientes com afecção das vias aéreas superiores foi maior no ano de 1997 (desvio
padrão igual a 2,02). Nos anos de 1998, 1999 e 2000 essa média foi, respectivamente,
igual a 6,92 (desvio padrão igual a 1,44), 6,42 (desvio padrão igual a 0,51) e 7,42
(desvio padrão igual a 1,08). Observando a representação gráfica da série AVAS (Figura
7.1) nota-se que em 1997, ano de maior ocorrência de pacientes com afecção das vias
aéreas superiores os picos ocorrem nos meses mais frios do ano, em 1998 e 1999 os
valores mais altos ocorrem nos últimos meses e em 2000 os picos ocorrem em fevereiro
e julho.
As concentrações dos poluentes SO
2
e CO também apresentam valor médio mais
alto em 1997: 19,39 (desvio padrão igual a 6,94) e 5,11 (desvio padrão igual a 1,08),
respectivamente – veja Tabela 7.
Tabela 7 – Média (
±
desvio padrão) para os dados AVAS e poluentes – SO
2
e CO.
Variáveis
N 1997 1998 1999 2000 Total
AVAS
48
8,50 (
±
2,02) 6,92 (
±
1,44) 6,42 (
±
0,51) 7,42 (
±
1,08) 7,31 (
±
1,55)
SO
2
48
19,39 (
±
6,94) 13,12(
±
4,02) 14,88 (
±
4,75) 17,24 (
±
4,84) 16,16 (
±
5,62)
CO
48
5,11 (
±
1,08) 4,21 (
±
0,64) 3,73 (
±
0,67) 3,54 (
±
0,84) 4,15 (
±
1,01)
Os
box-plots
não apontam concentrações aberrantes de nenhum dos dois
poluentes e evidenciam que os níveis de concentração de CO na cidade de São Paulo
diminuíram, ano após ano, no período avaliado.
A relação entre a variável AVAS e os poluentes SO
2
e CO estão representadas na
Figura 7.3. Tanto SO
2
quanto CO possuem correlação positiva com o número de
pacientes com afecções das vias aéreas superiores – correlação de Person iguais a
0,614 (p-valor aproximadamente 0) e 0,417 (p-valor igual a 0,003) – porém, da
avaliação da Figura 7.3 nota-se o fato do tipo de relação existente entre a resposta e o
poluente CO poder não ser estritamente linear. Um polinômio envolvendo um termo de
maior ordem ou uma função suavizadora podem ser, por exemplo, mais adequados
para descrever a relação entre a resposta e esse poluente.
50
Figura 7.3 – Representação gráfica da relação entre as séries AVAS e CO, AVAS e SO
2
.
A série de dióxido de enxofre apresenta uma estrutura auto-regressiva de ordem
1, conforme visto na Figura 7.4 e na Tabela 8 – note que o coeficiente AR(1) é
significativo tal como a constante; o teste Ljung-Box, Tabela 9, não apresenta indícios
que levem à rejeição de uma estrutura AR(1) na série SO
2
e a análise gráfica dos
resíduos (Figura 7.5) conduzem à conclusão de normalidade e homocedasticidade dos
mesmos.
Figura 7.4 – Gráficos de autocorrelação e autocorrelação parcial de SO
2
.
Tabela 8 – Modelagem AR(1) da série SO
2
.
Coeficiente Erro padrão T P
AR(1) 0,75 0,09 7,77 ~0
Constante 3,98 0,54 7,37 ~0
Média 16,23 2,20
SO2
AVAS
3025201510
12
11
10
9
8
7
6
5
RELAÇÃO AVAS E SO2
CO
AVAS
765432
12
11
10
9
8
7
6
5
RELAÇÃO AVAS E CO
PASSOS
AUTOCORRELAÇÃO
121110987654321
1,0
0,8
0,6
0,4
0,2
0,0
-0,2
-0,4
-0,6
-0,8
-1,0
FUNÇÃO DE AUTOCORRELAÇÃO DE SO2
(limites com 5% de significância para as autocorrelações)
PASSOS
AUTOCORRELAÇÃO PARCIAL
121110987654321
1,0
0,8
0,6
0,4
0,2
0,0
-0,2
-0,4
-0,6
-0,8
-1,0
FUNÇÃO DE AUTOCORRELAÇÃO PARCIAL DE SO2
(limites com 5% de significância para as autocorrelações parciais)
51
Tabela 9 – Teste Ljung-Box para modelagem AR(1) da série SO
2
.
Passos
12 24 36
Qui-quadrado 16,10 28,9 52,8
Graus de liberdade 10 22 34
p-valor 0,10 0,15 0,02
Figura 7.5 – Análise residual da modelagem AR(1) de SO
2
.
PASSOS
AUTOCORRELAÇÃO
121110987654321
1,0
0,8
0,6
0,4
0,2
0,0
-0,2
-0,4
-0,6
-0,8
-1,0
ACF DOS RESÍDUOS DA MODELAGEM AR(1) DE SO2
(limites com 5% de significância para as autocorrelações)
PASSOS
AUTOCORRELAÇÃO PARCIAL
121110987654321
1,0
0,8
0,6
0,4
0,2
0,0
-0,2
-0,4
-0,6
-0,8
-1,0
PACF DOS RESÍDUOS DA MODELAGEM AR(1) DE SO2
(limites com 5% de significância para as autocorrelações parciais)
RESÍDUO
PERCENTUAL
1050-5-10
99
95
90
80
70
60
50
40
30
20
10
5
1
NORMALIDADE DOS RESÍDUOS
Resposta SO2
VALORES AJUSTADOS
RESÍDUOS
27,525,022,520,017,515,012,510,0
10
5
0
-5
-10
RESÍDUOS x VALORES AJUSTADOS
Resposta SO2
MESES DE OBSERVAÇÃO
RESÍDUOS
454035302520151051
10
5
0
-5
-10
RESÍDUOS x MESES DE OBSERVAÇÃO
Resposta SO2
RESÍDUOS
FREQUÊNCIA
840-4-8
12
10
8
6
4
2
0
HISTOGRAMA DOS RESÍDUOS
RESPOSTA SO2
52
7.2 Modelagem MAG
A seguir são apresentados resultados quando um MAG é adotado para relacionar
a série AVAS aos poluentes SO
2
e CO, na forma
AVAS
)(~
i
poisson
µ
+
+
=
21i
SO)log(
β
α
µ
f
)
CO
(
432
)2(
β
π
β
β
+
+
+
niseni )nicos(
π
2
, (7.1)
n
i
,...,
1
=
. O método de suavização utilizado para a estimação do termo não
paramétrico foi o
loess
. A escolha do parâmetro de suavização,
λ
, foi feita a partir de
valores pré-fixados em 0,3, 0,5, 0,7 e 0,8 de forma que o modelo final apresentasse
valor mínimo para a estatística AIC. A Tabela 10 apresenta os resultados da estimação
pontual e intervalar do parâmetro linear e os resultados dos procedimentos
bootstrap
condicional e
bootstrap
nos resíduos – médias das 500 replicações
bootstrap
. Os
intervalos de confiança foram calculados com nível de 95% de confiança.
Tabela 10. Estimativas de
1
β
na modelagem MAG do AVAS
MAG
Bootstrap
condicional
Bootstrap nos
resíduos
1
ˆ
β
0,024
(0,0139)
0,024
(0,0141)
0,026
(0,0146)
AIC 208,10 - -
IC assintótico [-0,003;0,051] - -
IC bootstrap percentílico - [-0,006;0,054] [0,016;0,037]
IC bootstrap com correção
do vício
- [-0,004;0,057] [0,012;0,034]
Nota: Valores em negrito são o desvio padrão da estimativa de
1
β
.
A estimativa inicial obtida para o parâmetro linear
1
β
é igual a 0,024. Quando a
abordagem
bootstrap
é utilizada a média das estimativas de cada replicação apresenta
valor bastante próximo ao valor encontrado na modelagem inicial – 0,024 no
procedimento
bootstrap
condicional e 0,026 no procedimento
bootstrap
nos resíduos.
Quanto aos intervalos de confiança
bootstrap
, os percentílicos têm as menores
amplitudes, porém, no estudo simulado apresentado no Capítulo 6, foi mostrado que,
53
apesar de menos acurado, o IC’s
bootstrap
’s com correção do vício possuem
probabilidade de cobertura mais próxima ao nível nominal e apresentam-se levemente
superiores aos percentílicos, mas ambos com melhor desempenho que os intervalos
assintóticos.
As estimativas para os demais parâmetros envolvidos no modelo,
2
β
,
3
β
e
4
β
são todas significativas a um nível de 10%:
=
2
ˆ
β
-0,004 (desvio padrão igual a 0,007),
=
3
ˆ
β
0,029 (desvio padrão igual a 0,014) e
=
4
ˆ
β
-0,025 (desvio padrão igual a 0,089).
A análise gráfica dos resíduos da modelagem da série AVAS através do MAG está
representada na Figura 7.6; dos ACF’s e PACF’s é clara a presença de uma estrutura
auto-regressiva de ordem 1 nos resíduos indicando que o MAG não foi capaz de
capturar a estrutura de correlação presente no desfecho.
Figura 7.6 – Avaliação residual da série AVAS, modelada via MAG.
7.3 Modelagem MAG-AR(1)
Para corrigir a estrutura auto-regressiva presente nos resíduos da modelagem
MAG, a relação entre a série AVAS e os poluentes SO
2
e CO foram ajustados através do
MAG-AR(1) proposto neste trabalho; o modelo é da forma
0 5 10 15
-0.2 0.0 0.2 0.4 0.6 0.8 1.0
Lag
ACF
Series residuos
5 10 15
-0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4
Lag
Partial ACF
Series residuos
54
AVAS
)(~
i
poisson
µ
+
+
=
21i
SO)log(
β
α
µ
f
)
CO
(
432
)2(
β
π
β
β
+
+
+
niseni )nicos(
π
2
T
+
, (7.2)
onde
T
representa a estrutura de correlação presente nas observações e
n
i
,...,
1
=
.
O método de suavização utilizado na estimação do termo não paramétrico, os
valores pré-fixados de
λ
e o nível de confiança utilizado no computo dos intervalos de
confiança são os mesmo utilizados no ajuste do MAG na Seção 7.2 deste Capítulo. A
Tabela 11 apresenta os resultados da estimação pontual e intervalar do parâmetro
linear dos dados ajustados através do MAG-AR(1) conforme expressão (7.2).
As estimativas obtidas para o parâmetro linear
1
β
, tanto nos casos em
que a abordagem
bootstrap
é utilizada quanto no caso em que não há qualquer
reamostragem, apresentam valores muito próximos (iguais até pelo menos a terceira
casa decimal) às estimativas encontradas no ajuste do MAG.
As estimativas intervalares do parâmetro linear quando AVAS é modelada através
do MAG-AR(1) são bastantes semelhantes aos intervalos calculados na Seção 7.2, onde
os dados são ajustados via MAG. Da mesma forma, as conclusões da comparação entre
Tabela 11. Estimativas de
1
β
na modelagem MAG-AR(1) do AVAS
MAG-AR
Bootstrap
condicional
Bootstrap nos
resíduos
1
ˆ
β
0,023
(0,0139)
0,024
(0,0141)
0,026
(0,0147)
AIC 209,62 - -
IC assintótico [-0,004;0,050] - -
IC bootstrap percentílico - [-0,005;0,054] [0,015;0,037]
IC bootstrap com correção
do vício
- [-0,004;0,055] [0,012;0,035]
Nota: Valores em negrito são o desvio padrão da estimativa de
1
β
.
as duas técnicas
bootstrap
de construção de intervalo –
bootstrap
percentílico e
bootstrap
com correção dos vícios – são as mesmas consideradas na seção anterior.
As estimativas para os demais parâmetros envolvidos no modelo,
2
β
,
3
β
e
4
β
são bastante parecidas às estimativas calculadas na modelagem MAG além de serem
55
significativas a um nível de 10%:
=
2
ˆ
β
-0,005 (desvio padrão igual a 0,007),
=
3
ˆ
β
0,019
(desvio padrão igual a 0,014) e
=
4
ˆ
β
-0,022 (desvio padrão igual a 0,089).
A principal diferença entre a modelagem MAG e MAG-AR(1) está na capacidade
deste último de capturar estruturas de correlação entre as observações. Percebe-se, da
análise gráfica dos resíduos – Figura 7.7 –, que a estrutura auto-regressiva dos dados
foi captada pelo modelo e que os resíduos se tornaram um ruído branco.
Figura 7.7 – Avaliação residual da série AVAS, modelada via MAG-AR(1).
0 5 10 15
-0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0
Lag
ACF
Series residuos
5 10 15
-0.3 -0.2 -0.1 0.0 0.1 0.2 0.3
Lag
Partial ACF
Series residuos
56
8. CONCLUSÕES RELEVANTES
O Modelo Aditivo Generalizado e suas extensões constituem uma ampla classe
de modelos de regressão, na qual o efeito das variáveis preditoras na variável resposta
pode ser modelado de forma bastante flexível por meio de uma função não
especificada. Apesar de bastante utilizados em estudos de séries temporais,
principalmente em casos em que a variável resposta é uma contagem de eventos, os
MAG’s têm toda sua teoria de estimação e inferência construída sobre a hipótese de
independência dos dados. Apesar de não apresentar a flexibilidade do MAG, o GLARMA,
recentemente proposto na literatura estatística, tem a vantagem de ser um modelo
capaz de capturar a estrutura de dependência existente entre as observações de séries
temporais.
Neste trabalho é proposta uma nova classe de modelos, MAG-AR, que estende a
estrutura linear do GLARMA para a estrutura semiparamétrica do MAG acomodando
variáveis que tem relação não linear com a variável resposta em dados com estrutura
AR de dependência. Além disto, a técnica
bootstrap
foi utilizada para se fazer
inferências sobre o estimador dos parâmetros da parte linear do modelo.
O desempenho do MAG e do MAG-AR no ajuste de séries temporais foi,
empiricamente, comparado através de alguns experimentos Monte Carlo com o cálculo
do vício e do erro quadrático médio das estimativas do parâmetro linear dos modelos.
Os resultados mostraram estimativas mais consistentes e menos viciadas para o MAG-
AR(1) aplicados aos dados de séries temporais, quando comparados ao ajuste
utilizando o MAG. A performance de algumas abordagens
bootstrap
’s –
bootstrap
condicional e
bootstrap
nos resíduos – nas séries temporais simuladas foi avaliada e os
resultados indicaram que o
bootstrap
pode ser utilizado neste caso para se fazer
inferências intervalares sobre os parâmetros lineares do modelo por apresentar
resultados bastante semelhantes aos dos experimentos Monte Carlo.
Estimativas intervalares também foram calculadas e os intervalos de confiança
bootstrap
percentílico e
bootstrap
com correção do vício foram comparados ao intervalo
assintótico quanto à probabilidade de cobertura e o tamanho dos mesmos. Em geral os
intervalos
bootstrap
com correção do vício apresentaram resultados mais próximos à
57
cobertura nominal fixada (0,95) e os intervalos de confiança
bootstrap
tiveram melhor
desempenho que o intervalo assintótico.
As metodologias MAG e MAG-AR foram utilizadas no ajuste da relação entre a
série AVAS (número de pacientes com afecção das vias aéreas superiores) e os
poluentes dióxido de enxofre (SO
2
) e monóxido de carbono (CO). As análises
mostraram as mesmas conclusões sobre o efeito dos poluentes na variável resposta, no
entanto a estrutura de correlação presente entre as observações só foi capturada na
modelagem MAG-AR(1).
Como continuidade do trabalho, podem ser sugeridas pesquisas futuras que
incluam a extensão do modelo MAG-AR para estruturas de correlação mais complexas,
como a adição de termos médias móveis, no caso um MAG-ARMA, ou modelos de longa
dependência, MAG-ARFIMA.
58
REFERÊNCIAS BIBLIOGRÁFICAS
[1] Andrews, .F. (1974), A Robust Method for Multiple Linear Regression, Technometrics, 16, 523 – 531.
[2] Beaton, A.E.; Tukey, J.W. (1974), The fitting of power Series, Meaning Polynomials, Illustrated on Band-
Spectroscopic Data, Technometrics, 16, 147 – 185.
[3] Benjamin, R.A.; Rigby, M.A.; Stasinopoulos, M.D. (2003), Generalized autoregressive moving average
models. Journal of the American Statistical Association, 98(461), 214-223.
[4] Bhattacharya, P.K.; Zhao, P.L. (1997), Semiparametric inference in a partial linear models. The annals of
statistics, 25, 244-262.
[5] Box, G.E.P.; Jenkins, G.M. (1976), Time series analysis: forecasting and control. San Francisco: Holden-Day.
[6] Buja, A.; Hastie, T.J.; Tibshirani, R.J. (1989), Linear Smoothers and Additive Models. The Annals of Statistics,
16, 136-146.
[7] Cleveland, W.S. (1977), Locally Weighted Regression and Smoothing Scatterplots, Bell Laboratories
memorandum.
[8] Cleveland, W.S. (1979), Robust locally-weighted regression and smoothing scatterplots. J. Am. Statis. Assoc.,
74, 829 – 836.
[9] Conceição, G.M.S; Saldiva, P.H.N; Singer, J.M. (2001), Modelos MLG e MAG para análise da associação
entre poluição atmosférica e marcadores de morbi-mortalidade: uma introdução baseada em dados da cidade de
São Paulo. Revista Brasileira de Epidemiologia, 4.
[10] Cox D.D. (1983), Asymptotics for M-type smoothing splines. Ann. Statist. 11, 530-51.
[11] Craven, P.; Wahba, G. (1979), Smoothing noisy data with spline functions. Nunmer. Math., 31, 377-403.
[12] Davis, R.A.; Dunsmuir, W.T.M.; Wang,Y. (1999), Modelling time series of count data. Asymptotics,
Nonparametrics and Time Series, 63-114, New York.
[13] Davis, R.A.; Dunsmuir, W.T.M.; Street, S.B. (2003), Observation-driven Models for Poisson Counts,
Biometrika, 90, 4, 777-790.
[14] Deaton, A.; Muellbauer, J. (1980), Economics and Consumer Behavior. Cambridge University Press.
[15] Díez F.B.; Tenías, J.N.; Pérez-Hoyos, S. (1999), Efecto de la contaminación atmosférica sobre a salud: una
introducción. Rev Esp Salud Pública, 73, 109-121.
[16] Dominici, F.; McDermott, A.; Zeger, S.L.; Samet, J.M. (2002), On the use of generalized additive models in
time-series studies of air pollution and health. Am. J. Epidemiol. 156 (3), 193-203.
[17] Drescher, D. (2005). Alternative distributions for observation driven count series models. Economics Working
Paper, 11, Christian-Albrechts-Universitat Kiel.
[18] Efron, B. (1979) Bootstrap methods: Another look at the Jackknife. The Annals of Statistics, 7, 1-26.
[19] Efron, B.; Tibshirani, R. (1986), Bootstrap methods for standard errors, confidence intervals, and other
measures of statistical accuracy. Satatistical Science 1, 54-77.
[20] Efron, B.; Tibshirani, R. (1993), An introduction to the Bootstrap, New York: Chapman and Hall.
[21] Figueiras, A.; Roca-Pardiñas, J.; Suárez, C.C. (2005), A bootstrap method to avoid the effect of concurvity in
generalised additive models in time series studies of air pollution. J Epidemiol Community Health, 59, 881-
884.
[22] Goldberger, A.S. (1964), Econometric Theory. Wiley.
[23] Hall, P. (1988), Theoretical comparison of bootstrap confidence intervals (with discussion), Annals Statistical,
16, 927-953.
[24] Härdle, W. (1990), Applied Nonparametric Regression. Cambridge University Press. New York.
[25] Härdle, W.; Huet, S.; Mammen, E.; Sperlich, S. (2004), Bootstrap inference in semiparametric generalized
additive models, Econometric Theory, 20, 265-300.
[26] Hastie, T.J.; Tibshirani, R.J. (1990) Generalized additive models, London:Chapman and Hall
[27] Hastie, T. J. (1992), Generalized additive models. Capítulo 7 de Statistical Models in S eds J. M. Chambers
and T. J. Hastie, Wadsworth & Brooks/Cole.
[28] Hoaglin, D.C.; Welsch, R.E. (1978), The hat matrix in regression and ANOVA, Amer. Statist., 32, 17-22.
[29] Horowitz, J. (1998), Semiparametric methods in econometrics. Lecture Notes in Statistics 131, Springer,
Heidelberg, Berlin, New York.
[30] Lee, D.K.C. (1990), Cross-Validation in Semiparametric Models: Some Monte Carlo Results. Journal of
Statistical Computation and Simulation, 37, 171-187.
[31] Leontief, W. (1947), Introduction to a theory of the internal structure of functional relationships,
Econometrica, 15, 361-373.
[32] Lima, L.P.; André, C.D.S.; Singer, J.M. (2001), Modelos Aditivos Generalizados: metodologia e prática, R.
Bras. Estat., Rio de Janeiro, 62, 37-69
[33] McCullagh, P.; Nelder, J.A. (1989), Generalized linear models. 2.ed. London, Chapman and Hall.
[34] Moretin, P.A.; Toloi, C.M. (2004), Análise de séries temporais, São Paulo, Editora Blücher.
59
[35] Nelder, J.A.; Wedderburn, R.W.M.; (1972) Generalized linear models. J.R. Statis. Soc. A 135, 370-84.
[36] Opsomer, J.D.; Ruppert, D. (1999), A root-n consistent backfitting estimator for semiparametric additive
modeling. J. Comput. Graph. Statist., 4, 715-732.
[37] Opsomer, J.D. (2000). Asymptotic properties of backfitting estimators. J. Multivariate Anal., 73, 166-179.
[38] Paula, G.A. (2004), Modelos de regressão com apoio computacional, Instituto de Matemática e Estatística,
USP.
[39] Ramsay, T.O.; Burnett, R.T; Krewski, D. (2003), The Effect of Concurvity in Generalized Additive Models
Linking Mortality to Ambient Particulate Matter, Epidemiology,14.
[40] Rice, J.A.; Rosenblatt, M. (1983), Smoothing splines, regression derivatives and convolution. Ann. Statist, 11,
141-56.
[41] Rothman, K.J.; Greenland, S. (1998). Modern epidemiology. 2.ed., Philadelphia, Lippincott-Raven, 737.
[42] Saldiva, P.H.N.; Pope III, C.A.; Schwartz, j.; Dockery, D.W.; Lichtenfels, A.J.F.C.; Salge, J.M.; Barone, I.A.;
Böhm, G.M. (1995), Air pollution and mortality in elderly people: a time series study in São Paulo, Brazil.
Arch. Environmental Health, 50, 150-163.
[43] Schick, A. (1986), On Asymptotically Efficient Estimation in Semiparametric Models. The Annals of Statistics,
14, 1139 – 1151.
[44] Schick, A. (1993), On efficient estimation in regression models. The Annals of Statistics, 21, 1486-1521.
[45] Schick, A. (1996), Root-n-consistent and Efficient estimation in semiparametric additive regression models.
Statistical & probability Letters, 30, 45-51.
[46] Schwartz, J.; Slater, D.; Larson, T.V.; Pierson, W.E.; Koenig, J.Q. (1993), Particulate air pollution and hospital
emergency room visits for asthma in Seattle, Am Rev RespirDis, 147, 826-831.
[47] Schwartz, J. (1994), Nonparametric smoothing in the analysis of air pollution and respiratory illness. Canad. J.
Statist., 22, 4, 471-487
[48] Silverman, B.W.(1985), Some aspects of the spline smoothing approach to nonparametric regression curve
fitting (with discussion). J. Roy. Statist. Soc. Ser, 47, 1-52.
[49] Stone, C.J. (1977), Consistent nonparametric regression, Ann. Statist. 5, 549-645.
[50] Stone, C.J. (1985), Additive Regression and Other Nonparametric Models. The Annals of Statistics, 13, 689-
705.
[51] Whittaker, E. (1923), On a New Method of Graduation. Proceedings of the Edinburgh Mathematical Society,
41, 63-75.
[52] DETRAN-SP – Departamento Estadual de Trânsito de São Paulo. Frota de veículos. Disponível em:
<http://www.detran.sp.gov.br>. Acesso em 21 abr. 2009.
Livros Grátis
( http://www.livrosgratis.com.br )
Milhares de Livros para Download:
Baixar livros de Administração
Baixar livros de Agronomia
Baixar livros de Arquitetura
Baixar livros de Artes
Baixar livros de Astronomia
Baixar livros de Biologia Geral
Baixar livros de Ciência da Computação
Baixar livros de Ciência da Informação
Baixar livros de Ciência Política
Baixar livros de Ciências da Saúde
Baixar livros de Comunicação
Baixar livros do Conselho Nacional de Educação - CNE
Baixar livros de Defesa civil
Baixar livros de Direito
Baixar livros de Direitos humanos
Baixar livros de Economia
Baixar livros de Economia Doméstica
Baixar livros de Educação
Baixar livros de Educação - Trânsito
Baixar livros de Educação Física
Baixar livros de Engenharia Aeroespacial
Baixar livros de Farmácia
Baixar livros de Filosofia
Baixar livros de Física
Baixar livros de Geociências
Baixar livros de Geografia
Baixar livros de História
Baixar livros de Línguas
Baixar livros de Literatura
Baixar livros de Literatura de Cordel
Baixar livros de Literatura Infantil
Baixar livros de Matemática
Baixar livros de Medicina
Baixar livros de Medicina Veterinária
Baixar livros de Meio Ambiente
Baixar livros de Meteorologia
Baixar Monografias e TCC
Baixar livros Multidisciplinar
Baixar livros de Música
Baixar livros de Psicologia
Baixar livros de Química
Baixar livros de Saúde Coletiva
Baixar livros de Serviço Social
Baixar livros de Sociologia
Baixar livros de Teologia
Baixar livros de Trabalho
Baixar livros de Turismo
