Modelos de Previsão para Populações Raras
e Agrupadas sob Amostragem Adaptativa
TESE DE DOUTORADO
por
Kelly Cristina Mota Gonçalves
Universidade Federal do Rio de Janeiro
Instituto de Matemática
Departamento de Métodos Estatı́sticos
Modelos de Previsão para Populações Raras
e Agrupadas sob Amostragem Adaptativa
Kelly Cristina Mota Gonçalves
Tese de Doutorado submetida ao Corpo Docente do Instituto de Matemática Departamento de Métodos Estatı́sticos da Universidade Federal do Rio de Janeiro UFRJ, como parte dos requisitos necessários à obtenção do grau de Doutor em Estatı́stica.
Aprovada por:
Prof. Fernando A. S. Moura
PhD - UFRJ - Presidente.
Prof. Alexandra Mello Schmidt
PhD - UFRJ.
Prof. Mariane Branco Alves
PhD - UFRJ.
Prof. Heleno Bolfarine
PhD - USP.
Prof. Josemar Rodrigues
PhD - UFSCAR.
Rio de Janeiro, RJ - Brasil
2014
ii
CIP - Catalogação na Publicação
G635m
Gonçalves, Kelly Cristina Mota
Modelos de Previsão para Populações Raras e
Agrupadas sob Amostragem Adaptativa / Kelly
Cristina Mota Gonçalves. -- Rio de Janeiro, 2014.
143 f.
Orientador: Fernando Antônio da Silva Moura.
Tese (doutorado) - Universidade Federal do Rio
de Janeiro, Instituto de Matemática, Programa de
Pós-Graduação em Estatística, 2014.
1. Modelos de superpopulação. 2. Amostragem
informativa. 3. Modelos de mistura. 4. Inferência
bayesiana. I. Moura, Fernando Antônio da Silva,
orient. II. Título.
Elaborado pelo Sistema de Geração Automática da UFRJ com os
dados fornecidos pelo(a) autor(a).
iii
À minha mãezinha e ao meu paizinho (in memorian),
meus orgulhos.
iv
“Quero falar de uma coisa
Adivinha onde ela anda
Deve estar dentro do peito
Ou caminha pelo ar
Pode estar aqui do lado
Bem mais perto que pensamos
A folha da juventude
É o nome certo desse amor
Já podaram seus momentos
Desviaram seu destino
Seu sorriso de menino
Quantas vezes se escondeu
Mas renova-se a esperança
Nova aurora, cada dia
E há que se cuidar do broto
Pra que a vida nos dê
Flor, flor, e fruto
Coração de estudante
Há que se cuidar da vida
Há que se cuidar do mundo
Tomar conta da amizade
Alegria e muito sonho
Espalhados no caminho
Verdes, planta e sentimento
Folhas, coração,
Juventude e fé.”
Coração de estudante - Milton Nascimento.
v
Agradecimentos
Agradeço sempre em primeiro lugar a Deus pelo dom da vida e por iluminar meus
caminhos. Por estar ao meu lado em todos os momentos me protegendo e provendo várias
bençãos em minha vida. Sem Ele nada disso seria possı́vel.
À minha mãezinha Tereza por estar sempre ao meu lado cuidando de mim e torcendo
pelo meu sucesso. Agradeço por ser minha melhor companheira e por ter ajudado no diaa-dia para que eu pudesse dedicar-me exclusivamente a minha formação acadêmica nestes
anos. Ao meu paizinho Juarez (in memorian) pelo seu carinho e por ter se esforçado o
máximo para me dar educação. Sei que no céu o senhor está em festa e como sempre
cheio de orgulho da sua Kellynha. Meus pais amados, essa vitória também é de vocês!
Agradeço também aos tios e primos pela torcida e por terem estado sempre ao meu
lado, principalmente nos momentos em que mais precisei.
Ao meu orientador Fernando Moura, por acreditar em mim e estar sempre disponı́vel
para me ajudar. Meu crescimento durante estes 6 anos de trabalho juntos (entre mestrado
e doutorado) também se deve a você.
Ao meu amorzinho Andrés por sempre me apoiar em tudo e me dar o amor que muitas
vezes curou o meu estresse nestes anos. Obrigada por ser o anjinho que tornou meus dias
mais felizes nestes anos de muito estudo!
Aos professores do DME-UFRJ que passaram pela minha formação acadêmica
nestes anos.
Em especial ao professor Hélio Migon pela força e oportunidade de
trabalhar juntos em outros assuntos, e à professora Alexandra Schmidt pela torcida de
sempre e por ter incentivado a minha entrada neste programa de pós-graduação. Aos
vi
inesquecı́veis professores do IM-UFRJ que ajudaram a formar minha base matemática
nesta instituição.
Aos amigos que fiz durante estes anos de pós-graduação no DME-UFRJ. Em especial,
à Panela Camila, João, Larissa e Renata pela torcida e amizade verdadeira. À minha
turma Gustavo, João, Jony e Larissa pelo companheirismo nas disciplinas cursadas. Aos
demais amigos Patrı́cia, Mariana, Josiane, Vera (in memorian) e Felipe, veteranos que
estiveram sempre por perto. Agradeço a todos vocês pelos inesquecı́veis momentos que
passamos juntos. Grandes amizades que espero levar para toda a vida.
Agradeço também aos professores Alexandra Schmidt, Mariane Branco, Heleno
Bolfarine e Josemar Rodrigues por aceitarem participar desta banca.
Agradeço à CAPES pelo apoio financeiro, sem o qual não seria possı́vel realizar este
sonho. Ao GET-UFF pela flexibilidade, que me ajudou a exercer esta dupla jornada.
Agradeço também pelas experiências acadêmicas que tive no GET ao longo desses anos
e que me ajudaram a amadurecer em diversos aspectos.
Finalmente, agradeço à UFRJ, que tornou-se minha segunda casa nestes anos.
Quando entrei nesta instituição era uma menina de 17 anos ainda em dúvida sobre
sua carreira. Ao longo desses 9 anos aqui me graduei, encontrei uma área pela qual
me apaixonei, me tornei uma profissional e amadureci como pessoa. Sou profundamente
grata a esta instituição por hoje ser quem eu sou.
Ao escrever estes Agradecimentos a emoção algumas vezes tomou conta de mim, isso
mostra a importância desta conquista em minha vida. É um filme que passa na cabeça
neste momento. Obrigada a todos pela realização deste sonho!
vii
Resumo
Populações raras, como animais em extinção, pessoas infectadas por doenças raras,
usuários de drogas, entre outros, tendem a distribuir-se de forma agrupada em regiões.
Em levantamentos estatı́sticos com populações deste tipo, em que o principal interesse
é estimar o total populacional, este comportamento dificulta o processo de obtenção de
informação por meio de uma amostra aleatória simples, tornando-se necessários métodos
de amostragem complexos. Thompson (1990) propôs um método eficiente para estas
situações, denominado amostragem adaptativa por conglomerados.
Por outro lado, Rapley e Welsh (2008) propuseram uma abordagem para inferência em
populações deste tipo baseada em modelos. Sob o enfoque Bayesiano, o modelo proposto
é construı́do no nı́vel agregado dos grupos e incorpora o planejamento da amostragem
adaptativa por conglomerados à verossimilhança. Além disso, supõe homogeneidade entre
todas as unidades, mesmo as pertencentes a grupos distintos, o que resulta na frequência
esperada do total do fenômeno dentro de um grupo proporcional ao seu tamanho.
O objetivo deste trabalho é criar modelos alternativos para a previsão do total
populacional em uma determinada região. Inicialmente, o modelo agregado é estendido
para populações que evoluam dinamicamente.
Em particular, o interesse está em
populações raras que apresentam crescimento ou decrescimento dentro dos grupos até
a estabilização com a evolução do tempo.
Em seguida, o interesse é propôr um modelo de mistura alternativo ao modelo
agregado, que contemple situações mais gerais. A proposta é formulada em um nı́vel
desagregado da população, o que possibilita a inserção de estruturas com suposições
mais realistas, como a heterogeneidade entre grupos. O modelo é avaliado sob diversos
estudos de simulação e, finalmente, aplicado ao plano amostral adaptativo duplo, o qual
é um plano que permite a extração de mais informações acerca da população, mas sem
exceder os custos.
Palavras-chave: Amostragem informativa; modelos de mistura Poisson; RJMCMC.
viii
Abstract
Rare populations, such as endangered species, individuals infected by rare diseases and
drug users tend to cluster in regions. In many research studies with those populations,
where the main interest is to predict the population total, this behavior makes it difficult
the selection of a representative sample, making necessary complex sampling methods.
Thompson (1990) introduced an efficient method for these situations, called adaptive
cluster sampling.
On the other hand, Rapley e Welsh (2008) proposed a model-based approach to
make inference in those populations. From the Bayesian point of view, the proposed
model is built on the aggregated level of groups and takes into account the inclusion
probability of the adaptive sampling in the model likelihood. Furthermore, their model
supposes homogeneity between all units, even those belonging to different networks,
which is equivalent to assuming that the expected total in a group is proportional to its
size.
The aim of this work is to propose alternative models in order to predict the
population total in a region. Initially, the agregated model is extended to populations
that dinamically evolve. In particular, the interest is in rare populations which present
an increase or decrease within the groups, but stabilizes after some time.
Then, the interest is to propose a mixture model for more general situations,
alternative to the agregated model. The formulation of the model is done in the unit level,
what allows incorporating more realistic structures, such as the heterogeneity among units
belonging to different groups. The model is evaluated by carrying out some simulation
studies and finally applied to the adaptive cluster double sampling, which extracts more
informations about the population, without exceeding the costs.
Keywords: Informative sampling; Poisson mixture model; RJMCMC.
ix
Sumário
1 Introdução
1
1.1
Contribuições da tese . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4
1.2
Organização da tese . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5
2 Inferência em população finita
7
2.1
Introdução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7
2.2
Amostragem adaptativa por conglomerados . . . . . . . . . . . . . . . . .
9
2.3
2.4
2.2.1
Estimador do tipo Horvitz-Thompson modificado . . . . . . . . .
13
2.2.2
Amostragem estratificada adaptativa por conglomerados . . . . .
15
2.2.3
Amostragem adaptativa por conglomerados em dois estágios . . .
16
2.2.4
Custo operacional do plano amostral . . . . . . . . . . . . . . . .
16
2.2.5
Eficiência do plano amostral . . . . . . . . . . . . . . . . . . . . .
18
Modelos de superpopulação . . . . . . . . . . . . . . . . . . . . . . . . .
20
2.3.1
Desenho amostral informativo . . . . . . . . . . . . . . . . . . . .
21
Conclusões . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
23
3 Amostragem adaptativa por conglomerados baseada em modelos
3.1
3.2
25
Um modelo agregado . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
25
3.1.1
Possı́veis cenários gerados pelo modelo . . . . . . . . . . . . . . .
29
3.1.2
Estudo simulado para alguns cenários . . . . . . . . . . . . . . . .
30
3.1.3
Estudo simulado com população real . . . . . . . . . . . . . . . .
37
Um modelo para populações móveis, em crescimento ou decrescimento . .
40
3.2.1
41
Amostragem adaptativa para populações móveis . . . . . . . . . .
x
3.3
3.2.2
Incorporando estrutura de crescimento e decrescimento ao modelo
43
3.2.3
Modelo de crescimento exponencial . . . . . . . . . . . . . . . . .
45
3.2.4
Estudo simulado . . . . . . . . . . . . . . . . . . . . . . . . . . .
51
3.2.5
Comparação do modelo de crescimento com outras abordagens . .
55
Conclusões . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
58
4 Modelo de mistura para populações raras e agrupadas sob amostragem
adaptativa
60
4.1
Uma revisão sobre modelos de mistura de distribuições . . . . . . . . . .
62
4.1.1
Inferência Bayesiana em modelos de mistura . . . . . . . . . . . .
64
Modelo de mistura Poisson proposto . . . . . . . . . . . . . . . . . . . .
68
4.2.1
Distribuição a priori para λ . . . . . . . . . . . . . . . . . . . . .
72
4.2.2
Inferência para o modelo . . . . . . . . . . . . . . . . . . . . . . .
74
Estudo simulado . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
80
4.3.1
Considerando diferentes configurações . . . . . . . . . . . . . . . .
81
4.3.2
Considerando diferentes nı́veis de heterogeneidade . . . . . . . . .
84
4.3.3
Análise de sensibilidade da distribuição a priori . . . . . . . . . .
88
Comparação com o modelo agregado . . . . . . . . . . . . . . . . . . . .
91
4.4.1
Simulação baseada no desenho amostral . . . . . . . . . . . . . .
92
4.4.2
Simulação baseada no modelo . . . . . . . . . . . . . . . . . . . .
95
Modelo de mistura sob amostragem adaptativa dupla . . . . . . . . . . .
97
4.5.1
Amostragem adaptativa dupla . . . . . . . . . . . . . . . . . . . .
98
4.5.2
Modelo proposto sob amostragem dupla com variável auxiliar
4.2
4.3
4.4
4.5
indicadora de presença . . . . . . . . . . . . . . . . . . . . . . . .
4.5.3
99
Avaliação do modelo proposto sob amostragem adaptativa e
adaptativa dupla . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
4.6
Conclusões . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
5 Conclusões e trabalhos futuros
5.1
108
Trabalhos futuros . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
5.1.1
Planejamento amostral ótimo . . . . . . . . . . . . . . . . . . . . 110
xi
A Resultados dos modelos ajustados no Capı́tulo 3
112
A.1 Modelo (3.1) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112
A.2 Modelo de crescimento (3.4) . . . . . . . . . . . . . . . . . . . . . . . . . 116
B Cálculos envolvidos na inferência para o modelo proposto
118
B.1 Distribuições condicionais completas . . . . . . . . . . . . . . . . . . . . 118
B.2 Probabilidade de aceitação do algoritmo RJMCMC . . . . . . . . . . . . 121
xii
Lista de Tabelas
3.1
RaEQM e RaVAR dos estimadores para α, β, γ e T , entre os valores
obtidos no ajuste usando a probabilidade de seleção da amostra na função
de verossimilhança (3.3) e sem usá-la, sob 100 amostras artificiais. . . .
3.2
36
Estudo simulado com a população de marrecos da asa azul: eficiência
relativa para o estimador do total populacional com base no desenho
amostral adaptativo (estimador de Horvitz-Thompson modificado) e no
ajuste do modelo (3.1), com relação à amostragem aleatória simples de
tamanho n. A eficiência do estimador Bayesiano com relação ao estimador
de Horvitz-Thompson também é apresentada na última coluna. . . . . . .
3.3
40
Sumário da distribuição a posteriori dos parâmetros do modelo de
crescimento proposto: são apresentados o EQM e EAM, a amplitude média
dos intervalos HPD de 95% e a probabilidade de cobertura para as 100
populações geradas. Os resultados estão separadas para as populações em
crescimento e decrescimento. . . . . . . . . . . . . . . . . . . . . . . . . .
4.1
55
Análise da convergência das cadeias a posteriori dos parâmetros do modelo
proposto supondo distribuição a priori independente e dependente para λ
para uma população artificial. . . . . . . . . . . . . . . . . . . . . . . . .
4.2
78
Sumário a posteriori da estimação pontual e intervalar dos parâmetros do
modelo proposto e de T sob as 500 simulações, para diferentes valores de
α, β e N = 200. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
xiii
82
4.3
Sumário a posteriori da estimação pontual e intervalar dos parâmetros do
modelo proposto e de T sob as 500 simulações, para diferentes valores de
α, β e N = 400. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.4
83
Sumário a posteriori da estimação pontual e intervalar dos parâmetros do
modelo proposto e de T sob as 500 simulações, para diferentes valores de
α, β e N = 600. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.5
84
Sumário para a estimação pontual e intervalar dos parâmetros do
modelo e o total populacional para as 500 populações, variando o nı́vel
de homogeneidade nas redes, a partir do valor do CV fixado para a
distribuição de λ, para N = 400. . . . . . . . . . . . . . . . . . . . . . .
4.6
Análise da convergência das cadeias com a distribuição a posteriori dos
parâmetros dos modelos de mistura e agregado para a população real. . .
4.7
94
Sumário da estimação pontual e intervalar do total populacional obtido do
ajuste do modelo de mistura e do modelo agregado.
4.8
87
. . . . . . . . . . . .
95
Sumário a posteriori para a estimação pontual e intervalar dos parâmetros
dos modelos sob as 500 simulações onde λ foi gerado de uma distribuição
Gama com CV=50% e CV=25%, para N = 400 e (α, β) = (0.15, 0.10). .
4.9
96
Sumário a posteriori do total populacional T para os quatro planejamentos
considerados com base nas 500 amostras simuladas.
. . . . . . . . . . . 104
4.10 Resumo das principais conclusões acerca dos estudos simulados realizados
com o modelo de mistura proposto em (4.4). . . . . . . . . . . . . . . . . 107
xiv
Lista de Figuras
2.1
Ilustração do procedimento de amostragem adaptativa por conglomerados
para uma população rara e agrupada distribuı́da em uma região com 400
unidades.
No painel à esquerda temos uma amostra inicial de n1 =
10 unidades representadas pelos quadrados em cinza.
A partir desta
amostra, vizinhos são adicionados à amostra sempre que há pelo menos
uma observação (pontos em preto) na unidade selecionada, configurando
finalmente o plano amostral da direita. . . . . . . . . . . . . . . . . . . .
2.2
11
Ilustração dos conceitos importantes na amostragem adaptativa por
conglomerados: os quadrados com borda em negrito correspondem ao
conglomerado observado, os quadrados em cinza são as unidades da
rede e a parte hachurada as unidades da borda. A unidade selecionada
inicialmente está em cinza mais escuro. . . . . . . . . . . . . . . . . . . .
3.1
13
Populações artificiais geradas a partir do modelo proposto por Rapley e
Welsh (2008), para alguns valores fixos para os parâmetros α e β e para
γ = 10, numa grade regular de tamanho N = 400. . . . . . . . . . . . . .
3.2
31
População real de marrecos da asa azul na região da Flórida, nos Estados
Unidos, no ano de 1992, disposta numa grade regular de tamanho N = 200. 38
3.3
Ilustração da evolução dinâmica de interesse de uma população rara e
agrupada numa região sobreposta a uma grade regular com N = 400
unidades. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
xv
44
3.4
Curvas de crescimento e decrescimento de interesse para αt , t = 1, . . . , 50.
Em (a) fixou-se a = −1.73, b = −1.41 e c = −0.15, e em (b) a = −2.20,
b = 0.94 e c = −0.15, o que resulta no parâmetro αt variando de 0.05 e
0.15 e de 0.2 a 0.1, respectivamente. . . . . . . . . . . . . . . . . . . . . .
47
3.5
Distribuição a priori conjunta para o vetor (a, b)0 . . . . . . . . . . . . . .
49
3.6
Sumário da distribuição a posteriori de αt e do total populacional para
uma população em crescimento e decrescimento ao longo do tempo. Em
preto está a média a posteriori de αt e total populacional Tt , t = 1, . . . , 50,
com intervalo HPD de 95% em cinza e valor verdadeiro em azul. . . . . .
3.7
54
Sumário da distribuição a posteriori do total populacional a cada instante
de tempo T para 100 populações em crescimento e outras 100 em
decrescimento geradas. São apresentados os EQMR, EAR, probabilidade
de cobertura e amplitude média dos intervalos HPD de 95%. . . . . . . .
3.8
56
Comparação do modelo proposto de crescimento exponencial (3.4) com o
ajuste independente ao longo do tempo do modelo (3.1). Em (a) estão
as probabilidades de cobertura dos intervalos HPD de 95%, em (b) a
amplitude média destes intervalos, em (c) está a REQMR para cada
abordagem utilizada e em (d) as REQMR para todos os tempos incluindo
na comparação o estimador de Horvitz-Thompson. . . . . . . . . . . . . .
4.1
Comparação das médias da distribuição de Poisson e Poisson truncada no
zero. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2
58
75
Densidade a posteriori para alguns parâmetros do modelo proposto e para o
total populacional T com base em um dado artificial supondo distribuição
a priori para λ independente. A linha vertical cheia representa o valor
verdadeiro e a linha pontilhada o intervalo HPD de 95%. . . . . . . . . .
4.3
79
Densidade a posteriori para alguns parâmetros do modelo proposto e para o
total populacional T com base em um dado artificial supondo distribuição
a priori para λ dependente. A linha vertical cheia representa o valor
verdadeiro e a linha pontilhada o intervalo HPD de 95%. . . . . . . . . .
xvi
80
4.4
Erro relativo para λs e λs̄ ao longo de 500 simulações, para N = 400 e
diferentes configurações de α e β. . . . . . . . . . . . . . . . . . . . . . .
4.5
Distribuição a priori para λj usada nas simulações variando o valor do
CV da distribuição. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.6
85
85
Sumário da distribuição a posteriori de R assumindo diferentes
distribuições a priori para λ.
As cruzes representam a mediana da
distribuição a posteriori, o cı́rculo o valor verdadeiro de R e a linha o
intervalo HPD de 95%. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.7
89
EMQR para cada λj assumindo diferentes distribuições a priori para λ. Os
resultados com a distribuição a priori independente são representados pelos
cı́rculos vazios e a linha cheia, os resultados para a distribuição dependente
com τ = 5 são representados pelos triângulos e a linha tracejada, as cruzes
com a linha pontilhada representam os resultados quando τ = 10 e τ = 20
são os cı́rculos cheios e a linha traço e ponto. . . . . . . . . . . . . . . .
4.8
90
EQMR, probabilidade de cobertura e amplitude média do intervalo HPD
de 95% para o total populacional T sob cada distribuição a priori assumida
para λ e para cada valor de R fixado.
Os cı́rculos vazios e a linha
representam os resultados para R = 5, os triângulos com a linha tracejada
quando R = 6 e as cruzes com a linha pontilhada para R = 7. . . . . . .
4.9
91
Traço das cadeias com a distribuição a posteriori para α, β e T obtida do
ajuste do modelo de mistura (a) e do modelo agregado (b). A linha em
cinza representa o valor verdadeiro de T . . . . . . . . . . . . . . . . . . .
93
4.10 ER para T para as 500 amostras obtidos a partir do ajuste do modelo de
mistura e do modelo agregado. . . . . . . . . . . . . . . . . . . . . . . . .
95
4.11 Boxplot com o ER para T , a partir do modelo de mistura e do modelo
agregado para as 500 populações, tal que λ foi gerado de uma distribuição
Gama com CV=50% e CV=25%. . . . . . . . . . . . . . . . . . . . . . .
97
4.12 Sumário a posteriori de λs2 para os planejamentos (i) e (ii-a) com base
nas 500 amostras simuladas. . . . . . . . . . . . . . . . . . . . . . . . . . 105
xvii
1.1
Traços das cadeias dos parâmetros α, β, γ e total populacional T para um
dado artificial gerado fixando α = 0.05 e β ∈ {0.05, 0.1, 0.15, 0.2}, com
respectivos valores verdadeiros em cinza. . . . . . . . . . . . . . . . . . . 113
1.2
Traços das cadeias dos parâmetros α, β, γ e total populacional T para
um dado artificial gerado fixando α = 0.1 e β ∈ {0.05, 0.1, 0.15, 0.2}, com
respectivos valores verdadeiros em cinza. . . . . . . . . . . . . . . . . . . 113
1.3
Traços das cadeias dos parâmetros α, β, γ e total populacional T para um
dado artificial gerado fixando α = 0.15 e β ∈ {0.05, 0.1, 0.15, 0.2}, com
respectivos valores verdadeiros em cinza. . . . . . . . . . . . . . . . . . . 114
1.4
Traços das cadeias dos parâmetros α, β, γ e total populacional T para
um dado artificial gerado fixando α = 0.2 e β ∈ {0.05, 0.1, 0.15, 0.2}, com
respectivos valores verdadeiros em cinza. . . . . . . . . . . . . . . . . . . 114
1.5
Sumário da distribuição a posteriori dos parâmetros α, β, γ e T para
100 populações em 16 cenários com amostra inicial de 5%N e 10%N . Em
(a) os triângulos representam as probabilidades de cobertura dos intervalos
HPD de 95% para a amostra de 5%, os cı́rculos cheios para a amostra de
10% e a linha tracejada em vermelho o nı́vel nominal de 95%. Em (b)
estão o EQM para cada parâmetro e o EQMR para T . . . . . . . . . . . . 115
1.6
Sumário da distribuição a posteriori de Θ e do total populacional para
uma população em crescimento ao longo do tempo. Em (a)-(e) estão os
traços das cadeias da distribuição a posteriori dos parâmetros a, b, c, β e
γ. De (f )-(j) estão os traços das cadeias para os totais em alguns tempos.
A linha em cinza representa o valor verdadeiro usado na geração dos dados
artificiais. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116
1.7
Sumário da distribuição a posteriori de Θ e do total populacional para
uma população em decrescimento ao longo do tempo. Em (a)-(e) estão
os traços das cadeias da distribuição a posteriori dos parâmetros a, b, c,
β e γ. De (f )-(j) estão os traços das cadeias para os totais em alguns
tempos. A linha em cinza representa o valor verdadeiro usado na geração
dos dados artificiais. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
xviii
Capı́tulo 1
Introdução
Em diversos levantamentos estatı́sticos é possı́vel deparar-se com dificuldades na
coleta de dados, devido ao objeto de estudo ser difı́cil de ser observado. Isto pode ocorrer
simplesmente por ser um subconjunto pequeno da população toda, exibir um padrão de
grupos esparsamente distribuı́dos numa região, ou ainda por apresentar uma mobilidade
ao longo do tempo. São alguns exemplos de populações com estas caracterı́sticas: animais
e plantas em extinção, minorias étnicas, usuários de drogas, indivı́duos com doenças
raras e imigrantes recentes numa região. Problemas de monitoramento de populações
raras tornaram-se uma prioridade para muitos órgãos públicos, como por exemplo o
monitoramento de espécies ameaçadas de extinção para as agências de conservação.
Em geral detectar e estimar a abundância ou distribuição de populações com estas
caracterı́sticas é uma tarefa difı́cil.
Kalton e Anderson (1986) afirmam que populações raras são definidas basicamente
como uma pequena fração da população total, como por exemplo em estudos de
doenças raras, em que o interesse se concentra em grupos especı́ficos de sexo e idade.
No entanto, McDonald (2004) afirma que populações raras não são necessariamente
aquelas que possuem poucos indivı́duos, e sim aquelas em que os indivı́duos apresentam
comportamento elusivo ou estão esparsamente distribuı́dos em grandes espaços. Nesta
abordagem estão as populações raras e agrupadas, as quais apresentam um padrão de
distribuição espacial altamente concentrado, com grupos esparsos em uma região. Assim,
uma população com comportamento em forma de grupos espalhados em um espaço
1
geográfico grande tem uma raridade geográfica maior do que uma população de mesmo
tamanho confinada em um espaço geográfico menor.
A amostragem de populações raras é uma tarefa árdua, porque os custos de localização
de tais populações são substanciais e podem exceder os recursos disponı́veis. Além disso,
em geral, a densidade populacional média é pequena com relação à área total, mas quando
uma abundância substancial em alguns pontos é localizada, concentrações em vizinhanças
tendem a ser detectadas, e ao aplicar-se um planejamento amostral tradicional, muitas
unidades podem apresentar zeros na contagem, enquanto a maior parte das unidades
com contagens diferentes de zero se mantém concentrada em alguns locais que não foram
amostrados. Este fenômeno resulta em estimadores altamente imprecisos. Por esses
motivos, métodos especı́ficos têm sido desenvolvidos para a amostragem de populações
raras e agrupadas.
Em meio ao surgimento de diversas técnicas de amostragem para populações raras,
como as revisadas em Sudman e Kalton (1986), Kalton e Anderson (1986) e Kalton
(2001), a amostragem proposta por Thompson (1990) ganhou destaque na literatura
como uma técnica eficiente para levantamentos estatı́sticos em populações deste tipo.
Denominada como amostragem adaptativa por conglomerados, a técnica aproveita a ideia
intuitiva de que se os elementos da população foram encontrados em uma área, as áreas
vizinhas têm maior probabilidade de possuı́rem elementos com as mesmas caracterı́sticas.
Extensões desta técnica de amostragem podem ser vistas em Thompson e Seber (1996)
e Turk e Borkowski (2005).
Por outro lado, a biosfera está constituı́da de sistemas que mudam com o passar
do tempo, dependendo da organização do sistema e dos recursos disponı́veis. Kalton
(1991) revisa métodos de amostragem para populações móveis. O estudo da dinâmica
das populações naturais é importante para compreender o que ocorre nos ecossistemas em
equilı́brio. Da mesma forma, populações raras e agrupadas também podem apresentar
uma dinâmica populacional ao longo do tempo e tal fator pode ser gerador de dificuldades
maiores ainda nos levantamentos estatı́sticos.
McDonald (2004) apresenta estudos
por amostragem que produzem estimativas inadequadas simplesmente pelo fato do
pesquisador perder a população-alvo em um curto intervalo de tempo, devido ao grande
2
poder de deslocamento, mortes, entre outros fatores. Estudos acerca de populações de
animais selvagens constituem um campo de aplicação que em muitos aspectos difere de
levantamentos com uma população de árvores, por exemplo. Os animais podem circular
e se esconder naturalmente, e além disso o próprio processo de amostragem em si pode
induzir a esta mobilidade. Assim, um planejamento amostral eficiente pode não existir
e a probabilidade de inclusão de um animal na amostra é calculada depois da amostra
ter sido planejada. Por isso, a probabilidade de obter erros amostrais é também maior
em pesquisas com uma população de animais ou outra com esta mesma caracterı́stica.
Para estes e outros casos, um levantamento estatı́stico por amostragem, que considera
esta dinâmica da população e trabalha com coletas de amostras ao longo de um perı́odo
de tempo, pode produzir resultados mais precisos que planejamentos que não levem tal
dinâmica em consideração.
Todas as técnicas citadas acima fundamentam-se na teoria de amostragem baseada
na aleatorização do desenho amostral, ou seja, o mecanismo probabilı́stico de seleção da
amostra define um procedimento predeterminado de aleatorização, denominado desenho
amostral. Como apontado por Skinner et al. (1989), a principal razão desta abordagem
é sua caracterização como livre de distribuição.
Em algumas situações especı́ficas, como em estimação em pequenos domı́nios, esta
abordagem, baseada no desenho amostral, pode mostrar-se ineficiente, fornecendo
preditores inadequados. Isto porque neste caso, o tamanho da amostra resultante de
uma pesquisa é muito pequeno para que estimadores baseados somente no desenho
amostral apresentem precisão aceitável. Além disso, em termos de estimação intervalar,
é necessário recorrer ao Teorema Central do Limite, o qual não pode ser aplicado em
muitas situações práticas, em que o tamanho da amostra não é suficientemente grande
e/ ou no caso em que suposições de independência das variáveis aleatórias envolvidas
não são realistas. Uma possı́vel solução para estes casos é a utilização de modelos de
superpopulação. Nesta abordagem são usadas suposições explı́citas, buscando realizar
inferência sobre a parte desconhecida, que não seja baseada apenas na parte observada,
mas na distribuição conjunta das variáveis de interesse.
3
Com base nestas ideias, é possı́vel também fazer inferência em populações raras e
agrupadas usando as técnicas de amostragem citadas, mas sob a abordagem baseada
em modelos, em particular sob o enfoque Bayesiano. Nestes problemas a perspectiva
Bayesiana pode ter grandes vantagens sobre abordagens baseadas em desenho amostral
ou em modelos frequentistas, tais como: (i) podem-se obter estimativas para quantidades
para as quais a amostra coletada é pequena, incorporando informações a priori do
comportamento da população; (ii) a incerteza inerente ao procedimento de estimação
é levada em consideração na previsão, pois seguindo o paradigma de Bayes, é possı́vel
obter uma distribuição preditiva, entre outras.
Neste contexto, Rapley e Welsh (2008) propõem, de forma pioneira, um modelo,
sob o enfoque Bayesiano, que incorpora o planejamento da amostragem adaptativa por
conglomerados, a fim de inferir sobre o total populacional em uma região de interesse.
Uma caracterı́stica importante de tal modelo é que a unidade de análise é dada por um
nı́vel agregado de unidades menores, dessa forma trata-se de uma alternativa à introdução
das localizações espaciais, a fim de facilitar a inferência. No entanto, não incorporar
efeitos espaciais e estimar parâmetros populacionais em nı́veis agregados pode trazer
perdas de informações de interesse em nı́veis menores e na precisão das estimativas. Além
disso, duas suposições fortes deste modelo são que em média as unidades da população são
homogêneas com relação ao fenômeno de interesse e que o total esperado de ocorrências
do fenômeno em um determinado grupo é proporcional ao tamanho deste grupo na região.
1.1
Contribuições da tese
O objetivo deste trabalho é fazer previsões em populações raras, agrupadas e móveis
usando amostragem adaptativa por conglomerados, sob uma abordagem baseada em
modelos de superpopulação, sob o enfoque Bayesiano.
Primeiramente, o interesse está em estender o modelo de Rapley e Welsh (2008) com
o objetivo de fazer inferências sobre populações dinâmicas. Em particular, o interesse
está em populações em crescimento ou decrescimento que atingem a uma estabilização
com a evolução do tempo.
4
Em seguida, sem considerar evolução no tempo, é proposto um modelo para
populações raras e agrupadas, alternativo ao de Rapley e Welsh (2008). baseado em
misturas de distribuições.
Tal modelagem possibilita fazer inferência em um nı́vel
desagregado da população e suposições mais realistas, como por exemplo heterogeneidade
entre unidades que compõem grupos distintos.
Finalmente, esta proposta é estendida para problemas em que a amostragem
adaptativa por conglomerados torna-se muito custosa e faz-se necessário o uso de um
planejamento alternativo. Em particular, será considerada a amostragem adaptativa
dupla por conglomerados proposta por Felix-Medina e Thompson (2004). Neste contexto,
é considerada também a inserção de variáveis auxiliares que podem ajudar na estimação.
O software livre R (R Core Team, 2013) foi utilizado tanto para programar os
algoritmos quanto para a construção dos gráficos apresentados.
1.2
Organização da tese
No Capı́tulo 2 é introduzida a notação de amostragem de população finita, a qual
será utilizada ao longo do texto, e é feita uma ampla revisão de literatura sobre
planos amostrais informativos, modelos de superpopulação e amostragem adaptativa por
conglomerados.
No Capı́tulo 3 é apresentado o modelo proposto por Rapley e Welsh (2008), descrito
acima, o qual serviu-nos de inspiração para as propostas deste trabalho. Um estudo
simulado é apresentado, a fim de verificar o desempenho do modelo para alguns cenários.
Além disso, é apresentada uma população real, a qual é utilizada ao longo deste trabalho,
e em particular neste capı́tulo, esta é usada em uma avaliação do desempenho do modelo
em questão. Finalmente, é proposta uma extensão deste modelo para uma classe de
populações móveis e, em crescimento ou decrescimento, ao longo do tempo.
No Capı́tulo 4 é proposto um novo modelo de mistura de probabilidades para previsão
em populações deste tipo. Este modelo é mais geral que o proposto por Rapley e Welsh
(2008) pois modela as unidades desagregadas, o que permite prever neste nı́vel menor
e incorporar estruturas que acomodem suposições mais complexas para a população.
5
Alguns estudos simulados são apresentados a fim de avaliar o desempenho do modelo
proposto.
Experimentos baseados em modelos e desenho são feitos com o objetivo
de comparar o modelo proposto neste trabalho com o modelo de Rapley e Welsh
(2008).
Finalmente, é feita uma aplicação do modelo de mistura ao planejamento
amostral apresentado em Felix-Medina e Thompson (2004), o qual permite a realização
de pesquisas com um custo mais controlado e o uso de variáveis auxiliares.
Finalmente, o Capı́tulo 5 conclui o trabalho, resumindo o que foi desenvolvido e
apresentando propostas futuras.
6
Capı́tulo 2
Inferência em população finita
Neste capı́tulo são apresentados a notação e definições importantes na teoria de
amostragem de população finita que serão utilizadas ao longo deste trabalho. Neste
contexto, existem duas possı́veis abordagens: (i) a baseada na aleatorização do desenho
amostral, com a população fixa, e (ii) modelos de superpopulação (detalhes em Bolfarine
e Zacks (1992)). Na Seção 2.1 a primeira abordagem é apresentada. Em particular,
a Seção 2.2 apresenta um plano amostral utilizado para populações raras e agrupadas
proposto por Thompson (1990) e algumas extensões. Finalmente, na Seção 2.3 a segunda
abordagem é apresentada, com ênfase a modelos, para os quais o planejamento amostral
é relevante para a análise Bayesiana do modelo.
2.1
Introdução
Segundo Cassel et al. (1977), uma população finita é uma coleção de N unidades
denotada pelo conjunto de ı́ndices P = {1, . . . , N }, para a qual temos interesse numa
caracterı́stica y, para N supostamente conhecido. Associada à unidade i, i = 1, . . . , N ,
tem-se o valor yi . Se a unidade i é observada, não é somente o valor de yi que é registrado
mas, também, o fato de que foi exatamente a unidade i que gerou essa medida. Denote
a observação completa pelo par (i, yi ) e, portanto, existem N pares, (i, yi ), i = 1, . . . , N ,
para a população toda.
7
Defina y = (y1 , . . . , yN )0 como o parâmetro populacional da população finita. Por
exemplo, o número de pessoas com alguma doença em N bairros, ou o número de animais
de uma determinada uma espécie em N localizações. No contexto de populações finitas,
em geral o objetivo é estimar funções de y, como por exemplo o total populacional
P
0
T = N
i=1 yi = 1N y, onde 1N é o vetor unitário de dimensão N × 1, a média populacional
P
2
µ = T /N e a variância populacional σ 2 = N
i=1 (yi − µ) /N . Em particular, o interesse
neste trabalho concentrar-se-á em estimar o total populacional.
A inferência sobre estes parâmetros é feita com base em informações obtidas sobre
o vetor y por meio de uma amostra ordenada s ⊂ P , de tamanho n, dada por s =
{i1 , . . . , in }. A amostragem de população finita baseada na aleatorização do desenho
amostral distingue-se de outras partes da estatı́stica, pois trata a população de forma
fixa. Nesta abordagem, o mecanismo probabilı́stico de seleção da amostra define um
procedimento predeterminado de aleatorização, denominado desenho amostral. Este é
representado por uma função de probabilidade, conhecida como planejamento amostral,
definida no conjunto S de todas as possı́veis amostras s, onde [s] fornece a probabilidade
de selecionar a amostra s. Um desenho amostral [.] é chamado não informativo se, e
somente se, [.] é uma função que não depende dos valores de y associados a s. Denote
um planejamento amostral informativo por [s | y].
Uma vez que s é selecionada, o resultado observado pode ser especificado como o
conjunto de pares d = {(i, yi ) : i ∈ s}. Em alguns casos, o interesse está apenas nos
valores de y e não no par completo, por isso defina ys = {yi : i ∈ s}. Sejam s̄ = P − s e
portanto ys̄ = {yi : i ∈ P − s}, os valores de y que não pertencem à amostra.
Neste contexto, um conceito importante que virá a facilitar expressões mais à frente
é o conceito de consistência. De acordo com Cassel et al. (1977), uma amostra s é
0 0
dita consistente com uma particular população y0 = (y10 , . . . , yN
) se, e somente se,
yi = yi0 para todo i ∈ s. Em outras palavras uma amostra é consistente com uma
particular população se, e somente, se os valores de y das unidades amostradas coincidem
com os valores de y das mesmas unidades na população. Dessa forma, para qualquer
planejamento amostral dado por [.] e, qualquer vetor populacional y, tem-se que a
8
probabilidade de uma quantidade aleatória D tomar um valor d é dada por: [s], se
s é consistente com y e 0, caso contrário.
Analogamente, pode-se definir I como o vetor de dimensão N indicador de inclusão
na amostra s ⊂ S, de cada unidade da população, isto é Ii = 1 se i ∈ s e Ii = 0 se
i ∈
/ s. Note que Ii segue uma distribuição de Bernoulli com probabilidade de sucesso
πi , i = 1, . . . , N, tal que πi é a probabilidade de inclusão da unidade i na amostra.
Assim, por exemplo, o estimador de Horvitz-Thompson (Horvitz e Thompson (1952))
para o total T e sua variância podem ser escritos como:
T̂HT =
N
X
yi Ii
i=1
πi
,
V (T̂HT ) =
N
X
1 − πi
πi
i=1
yi2
N X
X
πij − πi πj
+2
yi yj ,
πi πj
i=1 j>i
(2.1)
tal que πij representa a probabilidade de inclusão das unidades i e j conjuntamente na
amostra.
A outra técnica usada na inferência em populações finitas é a baseada em modelos
de superpopulação, na qual a amostra permanece fixa, e as observações populacionais
são representadas por realizações de variáveis aleatórias, e a inferência se refere a uma
superpopulação hipotética, na qual uma lei de probabilidade governa as variáveis de
interesse. Esta metodologia também será vista com detalhes na Seção 2.3.
Na próxima seção é apresentado um planejamento amostral especı́fico, voltado para
levantamentos em populações raras e agrupadas.
2.2
Amostragem adaptativa por conglomerados
Em pesquisas dentro de regiões pode-se sobrepor uma grade regular e a seleção da
amostra envolve a seleção de um subconjunto de células da grade. Para populações
esparsas e agrupadas, a maioria das amostras de tamanho pequeno consistem de células
vazias, resultando em muitas amostras que geram estimativas imprecisas da quantidade
de interesse. A amostragem adaptativa por conglomerados é uma alternativa para esta
dificuldade pois trata-se de um planejamento voltado para populações raras e agrupadas.
Proposta inicialmente por Thompson (1990), o método mostrou-se eficiente em pesquisas
epidemiológicas, sobre doenças raras, com animais, plantas e de caráter social.
9
A técnica utiliza informações dos valores observados para ter mais êxito na coleta
de unidades da população, aumentado assim a eficiência dos estimadores. Isso se deve
ao fato de que se espera ser mais provável encontrar um elemento com caracterı́sticas
semelhantes a outro na sua vizinhança, quando a população é agrupada. Dessa forma,
este desenho caracteriza-se como informativo, pois a probabilidade de seleção da amostra
depende dos valores de y.
Na Figura 2.1 o método é ilustrado para uma população distribuı́da em uma região
particionada em uma grade regular no plano com N = 400 quadrados. Assim como
em Thompson (1990), defina os quadrados como unidades de observação primária e a
vizinhança de um quadrado como o conjunto de quadrados que apresentam um lado
contı́guo a este. Daqui em diante no lugar do termo quadrado será utilizado unidade. O
procedimento de amostragem inicia-se com a amostragem aleatória simples sem reposição
de n1 = 10 unidades, as quais estão dispostas em cinza na grade. Suponha que uma
unidade é classificada como de interesse se pelo menos uma observação é encontrada
nesta. Note que das 10 unidades selecionadas, apenas 2 satisfazem esta condição. Em
seguida, as unidades vizinhas a estas 2 unidades são também incluı́das na amostra.
O processo continua até que todas as unidades vizinhas com observações de interesse
sejam adicionadas à amostra e finaliza nas unidades vizinhas que não apresentem tais
observações. Observe na Figura 2.1 à direita o processo finalizado com n = 45 unidades
amostrais, representados pelas unidades em destaque.
Ainda que no exemplo descrito na Figura 2.1, a vizinhança tenha sido definida dessa
forma, outros tipos de vizinhanças podem ser consideradas, como por exemplo uma grade
sistemática em torno da unidade inicial, ligações genéticas e sociais no caso de populações
humanas, entre outras.
A condição para adição de vizinhos à amostra pode ser também definida de forma
mais geral como ter mais observações que um número mı́nimo fixado.
Além disso, note que à medida que as unidades vizinhas são agregadas à amostra,
em torno da primeira unidade selecionada é formado um grupo de unidades amostrais,
estes grupos formados são denominados conglomerados. Tal conglomerado só tem sua
fronteira finalizada até que vizinhos observados não satisfaçam a condição de interesse,
10
●
●
●
●●●●●
● ●●●
●
●
●
●● ●● ●● ●
●
●
●
●
●
●
●
●
●●●●●
● ●●●
●
●
●
●● ●● ●● ●
●
●
●
●
●
●
●
●●
●●
● ● ●●●
● ●●
●
●
●
●
●● ● ●
●● ●
● ●●●
●●●●●
●
●
●
●
●
●●●●
●● ●●●● ●
●●
●●
●
● ●
●
●
●
●●
●
●
●●
●●
● ● ●●●
● ●●
●
●
●
●
●● ● ●
●● ●
● ●●●
●●●●●
●
●
●
●
●
●●●●
●● ●●●● ●
●●
●●
●
● ●
●
●
●
●●
●
●
● ●●
●
●
● ●●●●●
●
● ● ●●● ●● ●
●
●● ●●
●
● ●
●
● ●●
●
●
● ●●●●●
●
● ● ●●● ●● ●
●
●● ●●
●
● ●
●
Figura 2.1: Ilustração do procedimento de amostragem adaptativa por conglomerados para
uma população rara e agrupada distribuı́da em uma região com 400 unidades. No painel à
esquerda temos uma amostra inicial de n1 = 10 unidades representadas pelos quadrados
em cinza. A partir desta amostra, vizinhos são adicionados à amostra sempre que há
pelo menos uma observação (pontos em preto) na unidade selecionada, configurando
finalmente o plano amostral da direita.
portanto todo conglomerado é formado por unidades na fronteira que não satisfazem tal
condição. Estas unidades são chamadas unidades de borda. Se uma unidade selecionada
na amostra inicial não é de interesse, não há acréscimos de vizinhos na amostra a partir
desta unidade.
Um conglomerado, descontadas as unidades de borda, é denominado rede. Note
que neste planejamento uma rede é sempre a mesma, independente da unidade da rede
selecionada na amostragem inicial.
Embora as unidades da amostra inicial selecionadas via amostragem aleatória simples
sem reposição sejam distintas, seleções repetidas podem ocorrer na amostra final quando
um conglomerado inclui mais de uma unidade na amostra inicial. Ou seja, se duas
unidades que não sejam de borda no mesmo conglomerado são selecionadas inicialmente,
então este conglomerado pode ocorrer duas vezes na amostra final. Uma unidade i da
11
população pode ser incluı́da na amostra tanto se qualquer unidade da rede a qual i
pertence é selecionada na amostra inicial, ou se qualquer unidade da rede a qual i é
uma unidade de borda é selecionada. Por definição as unidades que não satisfazem a
condição de interesse, assim como as unidades de borda, são também redes de tamanho
1. Portanto, uma amostra adaptativa por conglomerados, que se inicia com a seleção sem
reposição de n1 unidades iniciais, tem no final um número de redes não vazias distintas
sempre menor ou igual a n1 , mas note que o tamanho final da amostra é uma variável
aleatória e, portanto, não pode ser fixado.
A fim de ilustrar os conceitos de conglomerado, de rede e unidades de borda descritos,
na Figura 2.2 está uma parte da amostra vista na Figura 2.1. Os quadrados com borda
em negrito correspondem ao conglomerado observado, os quadrados em cinza compõem
a rede não vazia e a parte hachurada são as unidades da borda. A unidade selecionada
inicialmente está em cinza mais escuro.
Em geral, as redes é que são usadas como unidades de análise no lugar das células
da grade, pois as células da grade dentro de redes têm uma estrutura de dependência
e trabalhar no nı́vel de rede permite-nos evitar fazer esta estrutura de dependência de
forma explı́cita.
Segundo Cassel et al. (1977) um desenho amostral é chamado não informativo ou
ignorável se, e só se, a função planejamento amostral [.] não depende dos valores de y
associados aos ı́ndices em s. Desenhos informativos podem afetar as inferências quando
são erroneamente ignorados.
Note que o desenho adaptativo é informativo, pois a
probabilidade de seleção de uma amostra depende dos valores da variável de interesse.
Este tipo de planejamento será descrito com mais detalhes na Seção 2.3.
Estimadores convencionais sob este planejamento amostral tendem a ser viesados,
pois as unidades com observação de interesse são amostradas desproporcionalmente. Com
base nesta ideia, Thompson (1990) obteve um estimador não viesado sob este desenho
amostral para a média populacional, o qual está brevemente descrito a seguir.
12
●
●
●
●
●
●
●
●
●
●
●
●●
●●● ●
●
● ● ●●
●
●
●
●
●●
●
●
● ●●
●
●●●
●●
●
● ● ●●
●
●
●
●
●
● ●● ●
●
●●
●●
●
●
●
●
●
● ●
●
●
●
●
●●
●●
Figura 2.2:
Ilustração dos conceitos importantes na amostragem adaptativa por
conglomerados: os quadrados com borda em negrito correspondem ao conglomerado
observado, os quadrados em cinza são as unidades da rede e a parte hachurada as unidades
da borda. A unidade selecionada inicialmente está em cinza mais escuro.
2.2.1
Estimador do tipo Horvitz-Thompson modificado
Thompson (1990) apresentou um estimador não viesado para a média populacional
que corresponde a uma modificação do estimador de Horvitz-Thompson, no qual cada
observação yi na unidade amostral é dividida pela sua probabilidade de inclusão. Em
particular, será descrito a seguir o estimador do total populacional, que é uma simples
transformação da média.
Nesse caso uma unidade i é incluı́da na amostra se qualquer unidade da rede a qual
i pertence (incluindo ela mesma) é observada na amostra inicial, ou se qualquer unidade
da rede a qual i é uma unidade de borda é selecionada. Dessa forma, defina ai como o
número de unidades na rede para os conglomerados em que i é uma unidade de borda e
ci como o número de unidades na rede que contém i. Note que se i satisfaz a condição
de interesse, ou seja se i é uma unidade em cinza na Figura 2.2, tem-se ai = 0 e ci = 10.
Mas se i não satisfaz a condição de interesse, ou seja se i é uma unidade hachurada na
Figura 2.2, ci = 1 e ai = 10.
13
A probabilidade de inclusão da unidade i para qualquer uma das n1 seleções é dada
por
N − c i − ai
N
πi = 1 −
/
.
n1
n1
(2.2)
Note que, ao final do processo de amostragem, ci é uma quantidade conhecida para
as unidades amostradas, enquanto que ai pode ser maior do que o observado na amostra,
pois não temos o conhecimento se existe outra rede na qual i seja unidade de borda,
i = 1, . . . , N , tal que N é o número de unidades da grade. Portanto, o estimador de
Horvitz-Thompson para o total populacional em (2.1), com probabilidade de inclusão πi
dado por (2.2) não deve ser usado sob este desenho amostral.
Um estimador não-viesado para este caso pode ser obtido como uma modificação
do estimador de Horvitz-Thompson, apresentado em (2.1). O estimador faz uso das
observações que não satisfazem a condição de interesse só quando estas são observadas na
amostra inicial. Assim, a probabilidade de que uma unidade seja utilizada no estimador
pode ser calculada, mesmo se sua verdadeira probabilidade de inclusão seja desconhecida.
Portanto, defina a probabilidade de inclusão neste caso por:
N − ck
N
∗
πk = 1 −
/
,
n1
n1
em que ck é o número de unidades na rede que inclui a unidade k.
Seja a variável indicadora Ik∗ que assume o valor 0 se a unidade k na amostra s não
satisfaz a condição de interesse ou se k não foi selecionada na amostra inicial, e caso
contrário assume o valor 1. O estimador modificado portanto é dado por:
T̂HT ∗ =
ν
X
yk I ∗
k=1
k
∗
πk
,
(2.3)
em que ν é o tamanho efetivo da amostra final, ou seja o número de unidades distintas.
Para obter a expressão da variância do estimador é mais conveniente formulá-lo em
termo das redes do que das unidades individuais. Denote por N ∗ o número de redes na
população. Note que para toda unidade k da rede j, j = 1, . . . , N ∗ , Ik∗ é sempre a mesma,
portanto Ij∗ seria uma variável indicadora que assume o valor 0 se a rede j é vazia ou
se não foi observada na amostra, caso contrário assume o valor 1. A probabilidade de
14
inclusão πk∗ de uma unidade k é igual para todas as unidades na mesma rede j. Denote a
probabilidade de inclusão de uma rede j na amostra por αj . O total na rede j é definido
X
como yj∗ =
yk , em que Uj é o conjunto de unidades que compõem a rede j.
{k:k∈Uj }
Dessa forma, (2.3) pode ser reescrito como:
∗
T̂HT ∗ =
N
X
yj∗ Ij∗
j=1
αj
.
(2.4)
Note que como as redes são as unidades de análise neste caso, a fim de compatibilizar
∗
0
a notação com a Seção 2.1, o vetor populacional agora seria dado por y∗ = (y1∗ , . . . , yN
∗)
e o tamanho da população de interesse então deixaria de ser N um número conhecido e
passaria a ser N ∗ , um número desconhecido.
Para calcular a variância do estimador é necessário calcular a probabilidade αjl de
se selecionar duas redes simultaneamente, e dessa forma tem-se (detalhes em Thompson
(1990)):
∗
V (T̂HT ∗ ) =
∗
N X
N
X
yj∗ yl∗
j=1 l=1
em que αjl = 1 −
N −cj
n1
/
N
n1
−
N −cl
n1
/
N
n1
αj αl
−
(αjl − αj αl ),
N −cj −cl
n1
/
N
n1
.
A partir do trabalho de Thompson (1990), algumas extensões deste planejamento
amostral, além da seleção inicial baseada na amostragem aleatória simples, surgiram na
literatura e serão apresentadas a seguir.
2.2.2
Amostragem estratificada adaptativa por conglomerados
Uma das extensões naturais desta técnica de amostragem seria considerar o primeiro
estágio de amostragem não como uma amostra aleatória simples, mas como amostragem
estratificada. Tal extensão foi proposta em Thompson (1991). A amostragem adaptativa
tira vantagens de tendências de agrupamento da população, quando a localização e forma
dos conglomerados não podem ser previstos a priori. Enquanto a tradicional amostragem
estratificada (detalhes em Bolfarine e Zacks (1992)) é usada a fim de agrupar unidades
mais homogêneas entre si, baseada em informação a priori sobre a população ou na
15
simples proximidade das unidades. O planejamento amostral proposto combina estes
dois métodos.
Nesta abordagem a população é divida na grade em estratos e unidades dentro destes
estratos são selecionadas por amostragem aleatória simples. Se a unidade selecionada
satisfaz a condição, todas as unidades na sua vizinhança são observadas e a amostragem
adaptativa é realizada.
2.2.3
Amostragem adaptativa por conglomerados em dois
estágios
Proposta por Salehi e Seber (1997), esta é uma extensão do método introduzido em
Thompson (1991). Neste caso, a grade de tamanho N é particionada em M (M < N )
unidades primárias. Num primeiro estágio uma amostra de m das M unidades primárias é
selecionada sem reposição, num segundo estágio, observa-se nas m unidades maiores uma
amostra de unidades sem reposição. A partir destas unidades secundárias observadas,
a amostragem nas m unidades segue usando a técnica de amostragem adaptativa por
conglomerados. Note que quando m = M voltamos à metodologia de amostragem
estratificada adaptativa por conglomerados, pois todas as partições teriam amostras
coletadas.
2.2.4
Custo operacional do plano amostral
Assim como a amostragem por conglomerados convencional, a amostragem adaptativa
por conglomerados possui a vantagem de agrupar as unidades de análise em
conglomerados, o que minimiza o tempo e os custos de deslocamento. Mas se muitas
unidades na vizinhança satisfazem a condição de interesse, a amostra pode consistir da
maioria das unidades na população e, portanto, ser muito custosa. Logo, o esforço na
obtenção da amostra está associado à estrutura da população, e por isso é importante
que a população seja rara.
Algumas sugestões para a limitação do esforço na amostragem adaptativa são descritas
em Thompson e Seber (1996). Além disso, Brown e Manly (1998) propõem um método
16
chamado de amostragem adaptativa restrita por conglomerados, o qual limita o esforço
na obtenção da amostra e permite que uma aproximação para o tamanho da amostra final
seja obtida previamente. Na proposta, uma amostra inicial de tamanho fixo é selecionada
e amostragem adaptativa por conglomerados é feita. Se o tamanho da amostra final é
menor que um limite pré-definido, então outra unidade “inicial” é selecionada. Se incluir
esta unidade e sua vizinhança, caso a condição de interesse seja cumprida, resultar numa
amostra de tamanho maior que o limite pré-definido, então o conglomerado é incluı́do
na amostra mas nenhuma outra unidade é observada. Logo, esta metodologia exige
uma redução do tamanho da amostra inicial, para que esta produza uma amostra final
com tamanho próximo do limite desejado. Dessa forma, a variação no tamanho final é
reduzida e o planejamento dos esforços envolvidos na coleta de observações pode ser feito
com menos incerteza.
Por outro lado, também com o objetivo principal de controlar o número de medidas
da variável de interesse, Felix-Medina e Thompson (2004) introduziram a técnica de
amostragem adaptativa dupla por conglomerados, a qual combina ideias de amostragem
em dois estágios e amostragem adaptativa por conglomerados e exige a disponibilidade
de uma variável auxiliar mais fácil de medir. Na primeira fase a variável auxiliar é
usada para selecionar uma amostra adaptativa por conglomerados. Com a rede obtida
nesta primeira fase, são selecionadas subamostras subsequentes, as quais são obtidas
usando planos amostrais convencionais. Apenas nesta última fase os valores da variável
de interesse são registrados e estimativas para a média populacional, por exemplo, são
obtidas usando um estimador do tipo regressão.
Este plano amostral proposto permite ao pesquisador controlar o número de medições
da variável de interesse, alocar a subamostra na fase final próximo a lugares interessantes,
iniciar a coleta da segunda fase antes da primeira estar concluı́da e usar a variável auxiliar
na estimação.
Note que podem ser usados diferentes tipos de variáveis auxiliares neste caso, como
as de avaliação rápida que levam o pesquisador para as áreas mais promissoras, onde
observações exatas da variável podem ser feitas. Por exemplo, numa pesquisa sobre
mexilhões de água doce, a amostragem é feita a partir de mergulho para observar a
17
abundância de mexilhões. Assim, a variável auxiliar pode ser uma avaliação preliminar
da presença ou ausência de mexilhões, e a variável de interesse o número de mexilhões,
a qual é uma variável difı́cil de ser medida porque alguns mexilhões são parcialmente
escondidos pela areia e pedras no fundo do rio.
Note que este procedimento não controla o número de observações da variável auxiliar
e sim da variável de interesse. No entanto, em geral, procura-se escolher variáveis
auxiliares correlacionadas com a variável de pesquisa mas que sejam mais fáceis de serem
observadas e que produzam menos custos.
2.2.5
Eficiência do plano amostral
Ao comparar a eficiência da amostragem adaptativa por conglomerados com a
amostragem aleatória simples, por exemplo, Thompson e Seber (1996) notam que um
fator decisivo para uma maior eficiência relativa é a variabilidade dentro da rede.
Os estimadores sob o desenho da amostragem adaptativa por conglomerados, como o
apresentado em (2.4), não levam em conta a variabilidade dentro das redes pois a variável
resposta é dada pelos valores agregados dentro destas. Quanto maior essa variabilidade,
maior a vantagem, em termos de eficiência relativa, em usar amostragem adaptativa por
conglomerados do que a aleatória simples.
Portanto, conclui-se que, para que a amostragem adaptativa por conglomerados seja
um plano amostral eficiente em termos de precisão e custos é necessário que a população
de estudo exiba de fato um comportamento raro e agrupado. Logo, antes de propor
um planejamento amostral complexo como este, é importante conhecimentos a priori da
população em análise. Neste contexto, supondo que a variável y seja uma variável de
contagem do número de elementos que apresentam o atributo de interesse, para avaliar
a raridade da população pode ser utilizada a proporção de unidades contendo ao menos
um elemento da população rara, definida como:
PR =
N
1 X
I(yi > 0),
N i=1
18
(2.5)
onde I(.) é a função indicadora que assume o valor 1, se a unidade i apresenta ao menos
um elemento de interesse, e 0 caso contrário. Para avaliar a variabilidade dentro das
redes defina
PN ∗ P
V IR =
j=1
2
{i:i∈Uj } (yi − µj(i) )
,
PN
2
i=1 (yi − µ)
(2.6)
em que µj(i) é a média dos valores de yi nas unidades da rede que contém a unidade i e µ
é a média global da população. Note que se não há redes de tamanho maior que 1, tem-se
que V IR = 0, mas caso todas as unidades estejam numa única rede, V IR = 1. Dessa
forma, V IR pode ser considerada uma medida relacionada ao grau de agrupamento da
população.
Apresentamos portanto o método de amostragem adaptativa por conglomerados e
suas extensões propostas na literatura.
Vimos que o método é flexı́vel e pode ser
aplicado a diversos problemas estatı́sticos reais. No entanto, é importante ressaltar que a
eficiência do método depende da raridade e agrupamento espacial da população, portanto
é interessante o conhecimento prévio da população em estudo, dada a complexidade desta
metodologia. Smith et al. (2004) apresentam estas e outras questões práticas que devem
ser tratadas antes da proposta de tal planejamento num estudo por amostragem.
Alguns trabalhos na literatura mostram a eficiência deste tipo de amostragem
comparado a outros planos convencionais em aplicações a problemas reais, entre eles
podemos citar Thompson e Collins (2002), Danaher e King (1994), Smith et al. (1995),
Roesch (1993) e Conners e Schwager (2002).
A amostragem adaptativa por conglomerados fornece uma forma de lidar com
populações agrupadas sob o paradigma baseado no desenho amostral. Entretanto, sob
a abordagem baseada em modelo a metodologia de Rapley e Welsh (2008) é até então
a única proposta na literatura para este cenário. Na próxima seção é apresentada a
abordagem de modelos de superpopulação para um contexto geral.
19
2.3
Modelos de superpopulação
Outra abordagem de inferência, amplamente utilizada na literatura, para populações
finitas é a baseada em modelos de superpopulação.
Basicamente, o processo de
inferência estatı́stica a partir de uma amostra compreende um conjunto de princı́pios
e procedimentos que podem envolver, por exemplo, o conhecimento de algum processo
aleatório que possa ter gerado o verdadeiro valor desconhecido da caracterı́stica de
interesse para cada unidade da população. Esse processo é representado por um modelo
que é utilizado como base para se fazer inferência.
Enquanto na teoria convencional de amostragem as unidades da população são
tratadas como constantes fixas, não expressando nenhuma relação entre as unidades da
amostra e as unidades não amostradas, sob o enfoque de modelos de superpopulação, os
valores das caracterı́sticas de interesse são considerados realizações de variáveis aleatórias,
para os quais existe uma distribuição conjunta de todos os valores da população, a qual
é uma forma de expressar uma relação entre as unidades amostradas e não amostradas.
Logo, este enfoque complementa o planejamento amostral não informativo em relação às
unidades não amostradas. O vetor populacional y = (y1 , . . . , yN )0 é, portanto, tratado
como uma realização do vetor aleatório Y = (Y1 , . . . , YN )0 . A inferência clássica sobre
uma função do vetor populacional de interesse y procede com respeito à distribuição
amostral de uma estatı́stica, sob repetidas realizações geradas pelo modelo, com a amostra
selecionada permanecendo fixa. Esta forma de inferência em populações finitas pode ser
vista com maiores detalhes em Cassel et al. (1977).
Segundo o modelo, suponha que Y dado θ ∈ Θ segue uma distribuição de
probabilidades dada por [Y | θ]. Seja y = (y1 , . . . , yN )0 o vetor populacional gerado
segundo a distribuição [Y | θ]. Pode-se definir uma matriz H = (H1 , . . . , HN ) de
dimensão N × k, tal que Hi = (Hi1 , . . . , Hik )0 representa variáveis adicionais associadas
com a estrutura da população. Suponha que a distribuição conjunta de H, a qual depende
de um parâmetro φ ∈ Φ ∈ Rk , é dada por [H | φ].
20
2.3.1
Desenho amostral informativo
De forma mais complexa, o mecanismo de seleção amostral pode depender dos valores
das variáveis de interesse na população, ou seja, as probabilidades de inclusão das
unidades na amostra estariam relacionadas com as variáveis respostas. Tal situação
caracteriza um plano amostral informativo. Um exemplo tı́pico são os estudos de casocontrole, em que a amostra é selecionada de tal forma que haja casos (unidades com
determinada condição de interesse) e controles (unidades sem essa condição), sendo de
interesse a modelagem do indicador de presença ou ausência da condição em função de
variáveis preditoras. Esse indicador é uma das variáveis de pesquisa e é considerado no
mecanismo de seleção da amostra.
Sob a abordagem de modelos de superpopulação, é importante antes de propor
o modelo, analisar se as probabilidades de seleção dos elementos da população estão
relacionadas com as variáveis respostas, mesmo condicionado a covariáveis do modelo.
Neste caso, é relevante para inferência levar em consideração o plano amostral, seja na
definição do modelo ou na construção da função de verossimilhança.
Segundo, Gelman et al. (1995) é natural nestes casos expandir o espaço amostral e
incluir na verossimilhança o planejamento amostral. A verossimilhança completa, da
amostra s, do vetor Y, e das variáveis H pode ser escrita como:
[s, Y, H | θ, φ] = [s | Y, H][Y | H, θ][H | φ].
(2.7)
A expressão em (2.7) é avaliada em todos os valores da variável, mas na verdade a
real informação que tem-se a partir de uma amostra é (s, Ys , Hs ). A verossimilhança dos
dados observados, supondo continuidade, é dada por:
Z Z
[s, Ys , Hs | θ, φ] =
[s, Y, H | θ, φ]dYs̄ dHs̄
Z Z
=
[s | Y, H][Y | H, θ][H | φ]dYs̄ dHs̄ .
(2.8)
Já no caso discreto tem-se:
[s, Ys , Hs | θ, φ] =
X
X
{Yi :i∈s̄} {Hi1 :i∈s̄}
X
···
{Hik :i∈s̄}
21
[s | Y, H][Y | H, θ][H | φ].
(2.9)
Em particular, escolheu-se apresentar os demais resultados supondo variáveis
contı́nuas. Sob o enfoque Bayesiano, o interesse está na obtenção da distribuição a
posteriori do vetor paramétrico. Neste caso, a distribuição conjunta a posteriori dos
parâmetros (θ, φ), é dada por:
[θ, φ | s, Ys , Hs ] ∝ [θ, φ][s, Ys , Hs | θ, φ]
Z Z
= [θ, φ]
[s, Y, H | θ, φ]dYs̄ dHs̄
Z Z
= [θ, φ]
[s | Y, H][Y | H, θ][H | φ]dYs̄ dHs̄ .
A distribuição a posteriori de θ, em geral é a de maior interesse, e é obtida integrando
a expressão acima em φ, da seguinte forma:
Z Z Z
[θ | s, Ys , Hs ] ∝ [θ]
[φ | θ][s | Y, H][Y | H, θ][H | φ]dYs̄ dHs̄ dφ.
(2.10)
No caso de optar-se por ignorar o mecanismo de seleção da amostra, a distribuição a
posteriori de θ é dada por:
[θ | Ys , Hs ] ∝ [θ][Ys | Hs , θ][Hs | φ]
Z Z
= [θ]
[Y | H, θ][H | φ]dYs̄ dHs̄ .
(2.11)
Quando os dados não observados não fornecem informação adicional, ou seja, quando
[θ | Ys , Hs ] dada em (2.11) se iguala a [θ | s, Ys , Hs ] dada em (2.10), diz-se que o
desenho amostral é ignorável, por exemplo no caso da amostragem aleatória simples
com reposição. Entretanto, esquemas amostrais desse tipo são raramente empregados
na prática, por razões de eficiência e custo. Em vez disso, são geralmente empregados
planos amostrais que envolvem algum conhecimento da estrutura da população, como
a estratificação, conglomeração e probabilidades desiguais de seleção (amostragem
complexa).
Duas condições neste caso são suficientes para garantir ignorabilidade do desenho: (i)
[s | Y, H] = [s | Ys , Hs ]; (ii) [φ | θ] = [φ]. A importante consequência destas definições
é que, de (2.10), segue que, de fato, se o plano amostral é ignorável com respeito ao
parâmetro de interesse θ, [θ | s, Ys , Hs ] = [θ | Ys , Hs ]. Logo, a informação adicional
trazida por s pode ser descartada quando se deseja fazer inferência sobre θ, caso contrário
22
não pode ser eliminada. Ignorar erroneamente o plano amostral informativo na inferência
pode trazer consequências na estimação dos parâmetros.
Como consequência ainda se tem os seguintes resultados:
(i) se s é consistente com y então [s | Y] = [s | Ys ], e assim [s | Y] = [s] se, e somente
se, [s | Ys ] = [s];
(ii) se s é consistente com y, [s | Y, H] = [s | Ys , H] e diz-se que o planejamento
amostral é não informativo em relação a Ys̄ ;
(iii) se em (2.7) [s, Y, H | θ, φ] = [s | H][Y | H, θ][H | φ], diz-se que o planejamento é
informativo para H, mas não informativo para Y. Neste caso, se H é conhecido a
expressão em (2.8) pode ser reescrita da forma:
Z
[s, Ys , H | θ, φ] = [s | H][H | φ]
[Y | H, θ]dYs̄ .
Neste trabalho será amplamente utilizada a abordagem baseada em modelos de
superpopulação, discutindo a inferência sobre os parâmetros do modelo e previsão de
ys̄ a partir de dados obtidos por amostragem adaptativa por conglomerados, o qual é um
plano amostral informativo.
Como visto, a inferência para populações raras e agrupadas é usualmente abordada
com base no desenho amostral. De forma alternativa, Rapley e Welsh (2008) propõem
uma inferência neste contexto baseada em modelos usando a amostragem adaptativa.
Este plano amostral é informativo e, portanto, as ideias discutidas na Seção 2.3.1
são aplicadas a este modelo. Esta metodologia será apresentada no próximo capı́tulo,
juntamente com uma proposta de extensão do modelo para populações dinâmicas.
2.4
Conclusões
Neste capı́tulo foi feita uma revisão das duas possı́veis abordagens de inferência em
população finita. Como o objetivo deste trabalho é inferir acerca de populações raras
e agrupadas, o foco deste capı́tulo foi apresentar o plano amostral adaptativo e suas
extensões na literatura, por ser um plano amostral cabı́vel a este tipo de população. A
23
eficiência e o custo desta metodologia estão relacionados diretamente com a estrutura da
população em questão, portanto um conhecimento a priori pode auxiliar na construção
do planejamento amostral. Em particular, com relação ao custo operacional do método,
existem propostas na literatura, e algumas destas foram apresentadas neste capı́tulo.
Por outro lado, como o interesse deste trabalho é propor um modelo de
superpopulação para este contexto, fez-se necessário apresentar o conceito de plano
amostral informativo, pois este deverá ser relevante na construção da função de
verossimilhança do modelo neste caso.
24
Capı́tulo 3
Amostragem adaptativa por
conglomerados baseada em modelos
Como uma alternativa à inferência sobre o total populacional baseada nos planos
amostrais descritos anteriormente, Rapley e Welsh (2008) tratam tal problema sob uma
perspectiva baseada em modelos.
A inferência para este modelo fundamenta-se no
paradigma Bayesiano e leva em consideração o fato de que as unidades foram amostradas
de forma adaptativa por conglomerados, um plano informativo.
Na Seção 3.1 esta
metodologia é apresentada, o ajuste do modelo é estudado em alguns cenários e sua
eficácia é ilustrada para uma população real.
Na Seção 3.2 é proposta uma extensão deste modelo para populações em crescimento
ou decrescimento ao longo do tempo. Tal proposta é comparada com o ajuste do modelo
de Rapley e Welsh (2008) de forma independente ao longo do tempo.
3.1
Um modelo agregado
Rapley e Welsh (2008) propõem um modelo complexo, que usa as redes como unidades
de análise, de forma a não ter que introduzir componentes espaciais no modelo, o que
pode vir a facilitar a inferência. Portanto, por este motivo, nos referimos a este modelo
como um modelo agregado. O uso da abordagem Bayesiana é uma extensão natural da
ideia da amostragem adaptativa por conglomerados, pois incorpora o conhecimento a
25
priori de que a população é rara e agrupada tanto para a inferência como para o desenho
amostral. A fim de ilustrar a eficiência de sua proposta, Rapley e Welsh (2008) comparam
seus estimadores com os estimadores desenvolvidos em Thompson (1990) por meio de
um estudo de simulação, mostrando ser mais eficiente, principalmente num contexto de
conhecimento a priori. O modelo está descrito a seguir.
Seja Ω uma região que contém uma população esparsa e agrupada, na qual sobrepõese uma grade regular com N unidades. Uma unidade é dita não vazia se esta contém pelo
menos uma observação, e vazia caso contrário. Seja X ≤ N o número de unidades não
vazias em Ω. Seja R ≤ X o número de redes não vazias em Ω, Ci o número de unidades
não vazias dentro da rede i não vazia e portanto C = (C1 , . . . , CR )0 é o vetor com o número
P
de unidades não vazias dentro de cada rede não vazia. Logo X = R
i=1 Ci . Como existem
N − X unidades vazias, as quais são definidas como redes vazias de tamanho 1, então
há N − X + R redes em Ω. Dessa forma, pode-se estender o vetor de dimensão R para
Z = (C0 , 10N −X )0 em que 10N −X é um vetor de 1’s de dimensão N − X, logo Zi = Ci , se i
é uma rede não vazia e Zi = 1, caso contrário, para i = 1, . . . , N − X + R.
Seja Yi∗ o total observado na rede não vazia i e, portanto, Y∗ = (Y1∗ , . . . , YR∗ )0 denota
o vetor com o total populacional em cada uma das R redes não vazias. Também podemos
estender neste caso o vetor de dimensão R para um de dimensão N − X + R da forma
(Y∗ 0 , 00N −X )0 , em que 00N −X é um vetor de 0’s de dimensão N − X, o qual representa o
número de observações em cada rede vazia. O objetivo é fazer inferência sobre o total da
P
∗
população de interesse T = R
i=1 Yi .
Fazendo uma analogia com a notação definida na Seção 2.3 do Capı́tulo 2, note que
é possı́vel obter a seguinte relação: N ∗ = N − X + R, Hi1 = Ci e Hi2 = X, θ = γ,
φ = (α, β)0 e n = m. Note que apesar do tamanho da grade N ser conhecido, o tamanho
da população de interesse (redes não vazias), a qual está sendo modelada, ou seja, R,
é desconhecido e precisa ser estimado, portanto também pode ser interpretado como
Hi3 = R.
Isto é feito especificando a distribuição conjunta de X, R, C e Y∗ para a população
toda e o mecanismo de amostragem que fornece uma particular amostra s = {i1 , . . . , im }
de m redes das N − X + R redes na população. Um aspecto importante desta proposta
26
é que a estrutura da rede é totalmente determinada por X, R e C e não se faz necessário
modelar as localizações espaciais das redes.
Primeiramente modela-se a estrutura de rede vazia/ não vazia e então, condicional à
estrutura de rede, modela-se a contagem nas redes não vazias. Como o modelo aplica-se a
unidades não vazias, para evitar problemas de degeneração assume-se que há pelo menos
uma célula não vazia em Ω e, portanto uma rede não vazia, logo as distribuições são
truncadas à esquerda no valor igual a 1. Dessa forma, o modelo é dado por:
Yi∗ | Ci , R, γ ∼ Poisson Truncada independente (γCi ), Yi∗ ≥ Ci , i = 1, . . . , R,
R
X
1
C | X, R ∼ 1R + Multinomial X − R, 1R , Ci = 1, . . . , X − R + 1,
Ci = X
R
i=1
R | X, β ∼ Binomial Truncada (X, β), R = 1, . . . , X,
(3.1)
X | α ∼ Binomial Truncada (N, α), X = 1, . . . , N.
O truncamento na distribuição de Poisson também faz-se necessário para levar em
conta o fato de que cada unidade em uma rede não vazia deve conter ao menos uma
observação de interesse, logo Yi∗ ≥ Ci , i = 1, . . . , R.
Note que o parâmetro γ é
interpretado como o número médio de observações em cada célula não vazia, dentro
de cada rede não vazia na população. Vale ressaltar que a distribuição de Poisson pode
ser trocada por outro modelo, mas Rapley e Welsh (2008) mantiveram-se nesta proposta.
Além disso, um modelo log-linear comum não foi adotado para a variável resposta por
questões de custo computacional e problemas numéricos no ajuste, mas o uso de técnicas
mais eficientes de aproximação, tais como em Gilks e Wild (1992), poderia facilitar a
implementação deste modelo.
Este modelo é aplicado a amostras coletadas segundo o método adaptativo descrito
na Seção 2.2. Lembrando que o procedimento de amostragem consiste em observar Yi
para i ∈ s e seu delineamento depende da estrutura da população, a qual é desconhecida,
portanto este plano amostral caracteriza-se como informativo e deve ser incorporado à
função de verossimilhança do modelo para realização de inferência.
Logo, o próximo passo é definir a probabilidade de selecionar uma amostra s =
{i1 , . . . , im }, ou seja, [s]. Já vimos que tal mecanismo utiliza o argumento de que se
27
uma célula dentro de uma rede é amostrada, então toda a rede deve ser observada e,
portanto, a probabilidade de selecionar uma rede é proporcional ao seu tamanho. Para
motivar a construção da probabilidade de seleção de uma amostra, considere o seguinte
exemplo: seja uma população com 8 redes de tamanhos {5, 5, 1, 1, 1, 3, 3, 1} dos quais
obtemos a amostra {5, 1, 5, 3}. A probabilidade de selecionar a primeira unidade é igual
à probabilidade de selecionar uma unidade de tamanho 5, que é igual a 5 × 2/20, a
probabilidade de selecionar uma unidade de tamanho 1 no segundo passo, dado o anterior
é de 1 × 4/15 e, assim a probabilidade de seleção da particular amostra é igual a
1×4
× 20−5
×
5×1
20−5−1
×
5×2
20
×
3×2
.
20−5−1−5
Portanto, a probabilidade de seleção de uma particular amostra pode ser generalizada
da forma:
[s | C, R, X] =
m
Y
Zij × gij ,j
,
P
PN −X+R
Zi − j−1
k=0 Zik
i=1
j=1
(3.2)
onde gij ,j é o número de redes de tamanho Zij que restam após j − 1 redes terem sido
selecionadas e Zi0 = 0. Note que a probabilidade da seleção de s depende apenas das
variáveis associadas com a estrutura da população e não diretamente com Y∗ , logo, o
resultado (iii) da Subseção 2.3.1 se aplicaria neste caso e diz-se que o plano amostral é
informativo com relação a H.
Incorporando esta probabilidade de seleção da amostra ao modelo, tem-se por (2.7)
com [s | Y∗ , H] = [s | H], a seguinte função de verossimilhança global:
[s,Y∗ , C, R, X | α, β, γ] = [s | C, R, X][Y∗ | C, R, X, γ][C, R, X | α, β, γ]


m
Y
N
Zij × gij ,j
αX (1 − α)N −X


=
×
PN −X+R
Pj−1
1 − (1 − α)N
Zi − k=0
Zik
X
i=1
j=1


Ci −1
R
Y
X
β R (1 − β)X−R
1
1


×
× (X − R)!
X
1 − (1 − β)
(Ci − 1)! R
R
i=1
×
R
Y
i=1
(3.3)
exp{−γCi + Yi∗ log(γCi )}
.
P
i −1
Yi∗ ![1 − C
exp{−γC
+
j
log(γC
)
−
log(j!)}]
i
i
j=0
Com a amostra coletada, parte das variáveis do modelo é conhecida. Usando o ı́ndice s
para identificar a parte observada e s̄ a parte não observada, os vetores são particionados
da seguinte forma: Y∗ = (Ys∗ 0 , Ys̄∗ 0 )0 , C = (C0s , C0s̄ )0 , R = Rs + Rs̄ e X = Xs + Xs̄ .
28
A função de verossimilhança marginal dos dados observados é obtida somando a
expressão acima sob todas as quantidades desconhecidas, como visto em (2.9).
3.1.1
Possı́veis cenários gerados pelo modelo
A distribuição espacial da população ao longo da região é caracterizada no modelo
pelos parâmetros α e β. O parâmetro α controla o número esperado de unidades não
vazias, pois E(X | α) = N α/{1 − (1 − α)N } e β o número esperado condicional de redes
não vazias pois, E(R | X, β) = Xβ/{1 − (1 − β)X }. Note que se α se aproxima de 0 então
E(X | α) se aproxima de 1, que é o menor valor que X pode assumir segundo o modelo
proposto, mas se α está próximo de 1 então E(X | α) tende a N . De forma análoga
temos que, condicional a X, se β está próximo de 0 então E(R | β) está próximo de 1,
mas para valores de β perto de 1, E(R | β) tende a X, o número total de unidades não
vazias.
Como tratamos de populações esparsas, ambos os parâmetros são pequenos em
geral, e combinados, controlam a raridade e agrupamento destas. Populações raras são
caracterizadas pelo modelo para valores pequenos de α, enquanto populações agrupadas
estão caracterizadas para valores pequenos de β, mas este nı́vel de agrupamento depende
também do valor de X, o qual depende de α devido à estrutura condicional do modelo.
Além disso, as probabilidades da distribuição multinomial são tratadas como
conhecidas e iguais. Sob o modelo, o tamanho esperado, condicional a X e R, da rede é
1 + (X − R)/R = X/R.
Para ilustrar o impacto dos parâmetros no modelo, na Figura 3.1 temos alguns dados
artificiais gerados a partir do modelo para alguns valores fixos de α e β, γ = 10 e uma
grade regular de tamanho N = 400.
Observe que para α e β iguais a 0.05 tem-se uma população altamente rara, portanto,
intuitivamente, espera-se dificuldades de estimação numa população deste tipo, mesmo
utilizando a técnica de amostragem adaptativa por conglomerados.
Em contrapartida, para α e β iguais a 0.20 terı́amos uma população altamente dispersa
na região, o que estaria descaracterizando a raridade e agrupamento geográfico. Logo, o
uso deste modelo complexo não seria justificável.
29
Note também que fixando α igual a 0.05 e aumentando β, isto reflete uma população
com poucas unidades com observações, porém mais espalhada que o primeiro caso.
Finalmente, aumentando o valor de α e fixando β igual a 0.05, há um maior número de
unidades não vazias, o que ainda assim resulta em mais redes que o primeiro caso devido
à estrutura de condicionamento do modelo, diminuindo o grau de raridade espacial, mas
sem destruir o comportamento agrupado da população.
Note que como a partir do modelo não temos informação sobre a localização das redes,
na Figura 3.1 a localização destas foi feita de forma arbitrária e sem perda de generalidade,
sem comprometer a ilustração. Além disso, como estas populações foram geradas sob o
modelo agregado, não é possı́vel verificar o agrupamento da população usando a medida
em (2.6), pois nesta necessita-se da contagem em cada unidade da grade, o que não é
obtido na geração dos dados. Portanto, esta ilustração do comportamento do modelo
será feita apenas de forma visual.
A partir desta ilustração espera-se que populações raras e agrupadas possam ser
geradas a partir deste modelo para valores controlados de α e β. Lembre-se que temos
particular interesse em populações deste tipo, pois o interesse é explorar cenários em que,
com um custo controlado, a amostragem adaptativa possa ser mais eficiente, em termos
de precisão, que qualquer plano amostral não informativo e mais comumente utilizado.
3.1.2
Estudo simulado para alguns cenários
Como o procedimento de inferência baseia-se na metodologia Bayesiana, a fim de
avaliar o modelo apresentado por Rapley e Welsh (2008), foram analisadas amostras das
distribuição a posteriori dos parâmetros do modelo e do total populacional T . Para
isso o modelo proposto deve ser completado com uma distribuição a priori para o vetor
(α, β, γ). Supondo independência a priori entre estes, assume-se:
α ∼ Beta(aα , bα ), β ∼ Beta(aβ , bβ ) e γ ∼ Gama(aγ , bγ ),
em que Beta(a, b) representa a distribuição Beta parametrizada com média igual a
variância
a
a
b
ab
(a+b+1)(a+b)2
e variância
a
a+b
e
e Gama(a, b) a distribuição Gama parametrizada com média igual
a
.
b2
30
●
●
● ●
● ●
●
●
●●●
●
●● ●●●
●●
● ●
●
●●●
●
●●
●●
●
● ●●
●● ●
●
● ● ●●
● ● ●●
●
●
●
●●●●
● ●●
●
●
●
●
●
●
●
●
●
● ●
●
●
● ●● ●
●
●
●
●● ● ●
●●
●●●
●
● ● ● ●
●●
●
●● ●●
●●
●●
●
●●●● ●
●
● ●● ●●●
●●
●
●●
● ●●
●
●●●
●
●●●
● ●
●
●
●●●●●
●●
●●
●●●●
●●
●
●●
●
●
●●
●
●
●
●●
●
●●
● ●
●
●
●
●
●
●
●
●
●
●
● ●●
●
●●
● ● ● ●●
● ●
●
●
●
●
●
●
●●
●
●●●
●
●
●● ●
●
●●●
●
●
●
●
●●
●
● ● ●
●● ●
● ●
● ●
●●●●● ●
●
●
●
● ●●
●●●
●
●
●
●
●●
●
● ● ●
●●●
●●
● ● ● ●●
●
●● ●● ●
●
● ● ●●●
●
●
● ●●
●
●●
●
●●
●
● ●●●
●
●●
●
●●
●
●●● ●●
(α,β)=(0.05,0.05)
●
●
(α,β)=(0.05,0.20)
●● ●
● ●
● ●
●
●●
●
●
●●
●
●●
●
●
●●●
●●
●
●●●
●●
●
●
●●
●
●
●
●●●
●
●
●●
●
●●●
● ●●
●
●●
●●● ●
●
●
●
●
●
●
●●
● ●●
●
●● ●
●
●
●●
●●●
●●
●
●
●
●
●
● ●●
●
●
●●●
●
●
●
●
● ●●
●
●
●
●
●
●
● ●●
●
●
●●
●
●
●
●●
●●
●
●●
●
●
●●
●
●●
●
● ●●●
●
●●
●
●
●
●
●
●
●
●●● ●
●●
●●
●● ●
●
●●
●●●
●
●●●
●●
●
●●
●●
●●●
●
●
●●
●●
●
●● ●
●●●
●● ●● ●
●
●
●●
●
●
●● ●
●●● ●
●●
●
● ●●
●
●
●
●
●
●
●
●
●●
●● ● ●
●
●● ●
●
●
●●
●●●
●
●● ●
●
●●●●●
●●
●●●
●
●●●
●● ●
●●
●
●
●
●
●
●
● ●
●●
●
●
●
●
●
●
● ●●
●
●
●●
●●●●
●
●●● ●
●
●
●
●● ● ●
●● ●●
●
●
●●
● ●
●●
●
●
●
●●
●
●
●
●●●
●
●
●
●
●●
● ● ●
●●
●
●
●
●
●
●●
●
●
● ●
●●
●
●
●●●
●
●
●●●
●
●
●●●●
●
●
●● ●
●
●●
●
●
●●●
●
●
●● ● ●
●
● ● ●●
●●
●●●
●●
●
●●● ●
●●
●
●
●
●
●
●
●●
●●
●
●
●● ●●
● ●●●
● ●●
●
●
●
●
●
●● ●
●
●●●● ● ●●
● ●●
●
●● ●
● ●●
●
●
●
●
● ●
● ●
●
●
●●●●
●●●●
●
● ●●
●
●● ●
●
●
●
●
●
●●
●● ● ● ●
●
●
●
●
●● ●
●
●●
● ●
●● ●
● ●●
●
●
●
● ●●
●● ●
●●●●
●
●
●
●
●
●● ●●●●● ●
●
●●
●
● ●●
● ●●
●
●●
● ●●●
● ● ● ●● ●
● ●●
●
●●
●●● ●
●
●●
● ●
●●
●
●
●
●
●
●
●
●
●
●
●●
●
●●
●
●
●
●
●
●
● ● ●
● ● ● ●●●●
●
●
●
●●
●
●●● ●●
●
●●
●●
●
●●●
●●●●●●
●
●
●
●
●
●
● ●
●
●
● ●
●●● ●
●
●
● ● ●●
● ●●
●● ●
● ●●
●● ●
● ● ●
●●
●● ●
●
●●● ●
●● ●
●● ●
●●
●
●●● ● ●●●
● ●● ●●
●
●
●
●
●
●
●
● ● ●●●●
●●
●●
●● ● ●●●
●●
●
● ●
●●● ●●● ● ●
●
●
●
● ●● ●●
●● ●●
●●
● ●● ●
●
●●
●●
● ●●
●
●●
● ●● ● ● ● ●●
●
●
●●
●
●
●
●
●
● ● ●
● ● ●●
●
●●
●
●
●●
●●●
●● ●
● ●
●
● ●●
●
●● ●
●
●
●
●●●● ●●
●
●●●●
● ●
●
●
●
●
●● ●
●●
●●
●
●
●
●
●●
●●
●
●●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ● ● ●●
● ●●
● ●
●●
●
●●●
●
●●
●
●
●
●
●
●
●●
●●
●●
●
●
●
●●●
●
●
●
●
●
●●
●●●●
●
●
●
●
●
●
●
●
●
●●
●●
●
● ●●
●
●●
●
●
●
●● ●●
●
●●
●
●
● ●● ●
●
●
●●
●
● ●●
●●
●
●
● ●●
●●
●
●
●
●
●
● ●●● ●
●
●
●● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●● ●
●●
● ●●
●
● ●● ●
● ●●
●● ●
● ● ●●
●●
●●
●
● ●● ●
●
●
●
●●
●
● ●
●
● ● ●
● ●●● ● ● ●
●
●●
●● ● ●●
●
●●●
● ●
●● ●
●●
●
● ●●●
●●
● ●
●●
●
●
●
●● ●
●
● ●●
●●●
●●
●●
●●●●●
●●● ●
●
●●
●
●●●●●
●●
●● ●
●
●●
●●●
●●
●●●
● ●●●
●
●●●●
●
●● ●
●●
●
●●
●
●●
●●
●
●
●
●
●
● ●●●
●
●
●
●●●
●●
●
●●
●
●
●
●
●
●● ●
●●●
●●
●
●
●●●
●
●
●
●
●●
●
●●●
●
●●
●
●●
●
●●
●●
●● ●
●●
●●
●
●●● ●●●
●● ●●
●
●●●● ●●
●●●● ●
●●
●● ●
●●
● ● ●●●
●
●
●
●
●
●●
●
●
●
●●
●
●
● ● ●●
●●●●● ●●
●● ●●
●●
●
●
● ●
● ●
●● ●●● ●
●
●● ● ●
● ●●
●
●
●
●
● ●
●
● ●●
●● ●
● ● ● ●● ●●●
●●●●●
●
● ● ●●
●●●● ●
●● ●
● ●●
●● ●
●
●
●
●
●
● ●
●
●●
●●●● ●●
●
●
● ●
●●●
● ● ● ●●
●●●● ●
●
● ●
●
● ●●
● ●
●
●
●
●
●
●
●
●
●
●
●
●
●● ●● ●● ●
●
●
●●● ●●●
● ●●●●
●●
●
●● ●●
●●
●
● ● ●●
●
●
● ●
●
●
(α,β)=(0.20,0.05)
●
●
●
●
●
●
●
●● ● ●
●
●●●●
● ●●
● ●
●●●
● ●●
●●
●
●●
●
●
●
● ● ●● ●
● ●● ●
●
●
●
● ●●●
●●● ●●
●
●
● ●
●
●
●
● ●●
● ●
●●●● ●
●● ●
●● ●●
●
● ●
●● ●●●
●
●●●● ●
● ●
●●
● ●● ●●
●
●
●
● ● ●
● ●
●● ● ●
●
●
●
●
(α,β)=(0.20,0.20)
Figura 3.1: Populações artificiais geradas a partir do modelo proposto por Rapley e Welsh
(2008), para alguns valores fixos para os parâmetros α e β e para γ = 10, numa grade
regular de tamanho N = 400.
Rapley e Welsh (2008) fazem um estudo de elicitação da distribuição a priori para
estes parâmetros, avaliando a sensibilidade dos estimadores. Vale ressaltar que, ainda
sob distribuições a priori não informativas, o modelo fornece estimativas razoáveis para
os parâmetros e para o total populacional. No entanto, visto que o modelo é voltado
para aplicações a populações raras e agrupadas e dada a análise ilustrativa feita na
Figura 3.1, foram utilizados os seguintes valores: aα = aβ = 2 e bα = bβ = 9,
caracterizando distribuições a priori para α e β informativas. No entanto, neste contexto,
31
esta distribuição com alta probabilidade centrada em um intervalo apenas reflete a
priori a estrutura rara e agrupada da população, o que é o mı́nimo de conhecimento
para justificar o uso de tal modelo complexo. Para γ utilizou-se aγ = 1 e bγ = 0.1,
caracterizando assim uma distribuição a priori pouco informativa para γ, mas com
mais massa de probabilidade no valor médio de unidades por rede com base na amostra
selecionada, ou seja pela média do vetor Ys∗ /Cs .
Como a distribuição a posteriori do vetor paramétrico Θ = (Xs̄ , Rs̄ , Cs̄ , Ys̄∗ , α, β, γ)
não possui forma analı́tica fechada faz-se necessário o uso de métodos de simulação
estocástica, como o método de Monte Carlo via Cadeias de Markov (MCMC). Em
particular, o amostrador de Gibbs com passos de Metropolis-Hastings foi utilizado. Além
disso, o preditor do total populacional T é dado por:
T̂ = 10Rs Ys∗ + 10R̂s̄ Ŷs̄∗ ,
cuja amostra da distribuição a posteriori também pode ser obtida via MCMC.
Os passos da amostragem são descritos por:
(1) faça j = 1 e especifique valores iniciais para Xs̄ , Rs̄ , Cs̄ e Ys̄∗ ;
(2) sorteie α da distribuição condicional completa [α | X, R, C, Y∗ , β, γ] = [α | X];
(3) sorteie β de [β | X, R, C, Y∗ , α, γ] = [β | X, R];
(4) sorteie γ de [γ | X, R, C, Y∗ , α, β] = [γ | R, C, Y∗ ];
(5) sorteie (Xs̄ , Rs̄ , Cs̄ , Ys̄∗ ) de [Xs̄ , Rs̄ , Cs̄ , Ys̄∗ | Xs , Rs , Cs , Ys∗ , α, β, γ];
(6) faça j = j + 1 e volte ao passo (2).
As condicionais completas e as distribuições propostas podem ser vistas com
detalhes em Rapley e Welsh (2008). A fim de mostrar a eficiência do modelo para a
previsão do total populacional, foram geradas algumas populações raras e agrupadas
artificiais, para alguns valores fixos dos parâmetros, e o modelo (3.1) foi ajustado a tais
dados. Dessa forma é possı́vel comparar a estimativa do total com o valor verdadeiro
gerado. Cada população foi simulada numa grade regular com N = 400 unidades.
32
Populações foram geradas para 16 cenários diferentes a partir das combinações de
α, β ∈ {0.05, 0.10, 0.15, 0.20} e γ = 10. Para cada valor dos parâmetros gerou-se 100
populações, e de cada uma selecionou-se uma amostra adaptativa com dois tamanhos
iniciais distintos de 5%N e 10%N . Vale ressaltar que, apesar da amostra aleatória simples
inicial ser de 20 ou 40 unidades, o número de redes observadas ao final da amostragem
adaptativa era menor ou igual a esse número, pois em alguns casos duas ou mais unidades
selecionadas faziam parte da mesma rede na população.
As Figuras 1.1, 1.2, 1.3 e 1.4 no Apêndice A apresentam as trajetórias das cadeias
obtidas para cada parâmetro e para o total populacional T com o respectivo valor
verdadeiro em cinza, para uma das 100 populações geradas com amostra inicial de
tamanho 10%N . Para todas as cadeias foram geradas 200.000 iterações, sendo as 10.000
primeiras descartadas como aquecimento e foram tomadas amostras de 190 em 190, a
fim de obter-se 1.000 amostras independentes. Há indı́cios de convergência para todos
os 16 cenários simulados, visto que as cadeias são estacionárias e movem-se em torno do
valor verdadeiro fixado na geração dos dados. O mesmo ocorre quando seleciona-se uma
amostra adaptativa de tamanho inicial n1 = 5%N .
Na Figura 1.5 no Apêndice A estão um sumário da distribuição a posteriori dos
parâmetros α, β, γ e de T para as 100 populações artificiais para cada um dos 16 cenários
gerados a partir do modelo e para os dois tamanhos de amostra distintos. Tais cenários
estão na seguinte ordem na figura: fixa-se um valor de α e depois varia-se β. A Figura
1.5 (a) apresenta uma análise de propriedades frequentistas dos estimadores. Nela estão
as probabilidades de cobertura dos intervalos HPD de 95% para as amostra de 5%N e
10%N , o erro quadrático médio (EQM) para cada parâmetro e o erro quadrático médio
relativo (EQMR) para o total populacional. Os intervalos HPD apresentados ao longo
deste trabalho foram obtidos usando o comando emp.hpd do pacote TeachingDemos do
software R.
Note que em termos da cobertura média dos intervalos, enquanto os parâmetros β e γ
apresentam resultados próximos do desejado para todos os cenários, o parâmetro α tem
maior variabilidade e resultados mais satisfatórios são obtidos no geral à medida que o
valor de α aumenta, para β não muito pequeno. O mesmo se passa com a estimação do
33
total populacional T . Isto ocorre pois quanto maior α, mais unidades com observações de
interesse, o que traz mais informações que auxiliam na estimação e previsão. Por outro
lado, analisando o EQMR de T , que é o nosso maior interesse, observa-se também que,
fixado α, no geral os valores do EQM e EQMR diminuem à medida que β aumenta. Isto
ocorre pois o parâmetro β está associado ao número de redes e à medida que β aumenta,
cresce o número de redes, fazendo com que os grupos na população se espalhem mais, o
que também facilita o procedimento de inferência com base numa amostra. Neste mesmo
caso, observe que, mesmo com uma amostra de 5%N o modelo já se ajusta bem aos
dados e as conclusões são análogas.
Uma alternativa para melhorar o ajuste deste modelo sob cenários em que α e β são
extremamente pequenos é elicitar outras distribuições a priori independentes para α e β.
Rapley (2004) apresenta uma lista de distribuições a priori utilizadas e que resultaram
num melhor desempenho do modelo em populações geradas para diferentes valores de α e
β. Neste trabalho, foi utilizada apenas uma distribuição a priori informativa para todos
os cenários, com o único interesse de garantir que o desenho amostral seja razoável ao
problema e a robustez desta para diferentes valores de α e β num intervalo. Portanto, é
uma possı́vel distribuição a priori a ser utilizada quando o único conhecimento prévio que
se tem a respeito da população é que esta é rara e agrupada. No entanto, se informações
mais precisas sobre o tipo ou estrutura da população estão disponı́veis, resultados mais
vantajosos podem ser obtidos para alguns casos especı́ficos.
É importante mencionar que este estudo simulado foi feito sob todas as possı́veis
amostras. Por exemplo, nos casos em que α = β = 0.05 a população é extremamente rara
e agrupada, portanto é alta a probabilidade de selecionar uma amostra que não contenha
unidade alguma com observação ou que contenha todas as unidades não vazias da
população. Isso prejudica a qualidade das estimativas. Esta é mais uma explicação para
o fato de que os resultados são mais próximos do desejado para populações menos raras
e agrupadas. Uma possibilidade para este caso é repetir o estudo simulado descartando
estas amostras não representativas, no entanto, como elas têm alta chance de ocorrer
em alguns casos optou-se por mantê-las, a fim de não mascarar estes problemas nos
resultados.
34
Conclui-se desta forma que, em termos de estimativas pontuais e intervalares, a
eficiência do modelo aumenta à medida que os valores de α e β aumentam. Entretanto,
é importante lembrar que a amostragem adaptativa pode ser custosa, portanto esta é
razoável em cenários de raridade e agrupamento da população. Logo, recomenda-se o
uso de tal modelo complexo nestes cenários, mas com um número esperado controlado
de unidades e redes com a caracterı́stica de interesse, de forma que a amostra adaptativa
coletada seja a mais representativa possı́vel sem altos custos.
Por outro lado, já foi visto que como o plano amostral adaptativo por conglomerados é
não-ignorável, a probabilidade de seleção deve ser incluı́da na função de verossimilhança,
pois esta também traz informações para a estimação dos parâmetros do modelo. O
objetivo agora é simplesmente verificar o ajuste do modelo para o caso em que o plano
amostral é erroneamente considerado ignorável, ou seja, quando a probabilidade de
seleção é descartada da função de verossimilhança completa em (3.3).
Para isso, o modelo em (3.1) foi ajustado para as mesmas 100 amostras do estudo
anterior, mas, agora, desconsiderando a probabilidade de seleção da amostra em (3.3).
Na Tabela 3.1 é apresentada uma comparação entre as duas abordagens usando
a razão dos EQM (RaEQM) e das variâncias (RaVAR) entre os estimadores obtidos
considerando a probabilidade de seleção e sem considerá-la, para n1 = 10%N .
Portanto, valores menores que 1 indicam que considerar o plano amostral na função de
verossimilhança produz resultados mais vantajosos sob ambos os critérios. Vale informar
que as probabilidades de cobertura para os intervalos HPD de 95% gerados para os dois
métodos apresentam-se próximo do nı́vel nominal desejado, logo não seriam um critério
relevante na comparação e, por isso, não foram apresentadas.
Observando a Tabela 3.1 é possı́vel verificar que desconsiderar esta parcela na função
de verossimilhança completa, gera na grande maioria das vezes, estimativas viesadas
e com maior variância, principalmente para o parâmetro α e para o total T . Apenas
para dados artificiais gerados a partir do modelo fixando α = β = 0.20 esta conclusão
é diferente em termos do EQM para todos os parâmetros. Contudo, a variância ainda
permanece menor quando incluı́da a probabilidade de seleção. Isso ocorre pois, este
cenário gera uma população mais esparsa, e menos rara que os outros cenários estudados.
35
Logo, fazer uma amostragem não informativa, como a aleatória simples por exemplo, ou
adaptativa, teria o mesmo efeito, e não justificaria assim o uso do modelo complexo.
Tabela 3.1: RaEQM e RaVAR dos estimadores para α, β, γ e T , entre os valores obtidos
no ajuste usando a probabilidade de seleção da amostra na função de verossimilhança
(3.3) e sem usá-la, sob 100 amostras artificiais.
(α, β) fixos
α
β
γ
T
RaEQM RaVAR
RaEQM RaVAR
RaEQM RaVAR
RaEQM RaVAR
(0.05, 0.05)
0.26
0.10
1.16
1.32
0.69
0.97
0.25
0.07
(0.05, 0.10)
0.23
0.10
1.08
1.31
1.17
3.22
0.23
0.12
(0.05, 0.15)
0.24
0.11
1.12
1.26
1.07
2.28
0.23
0.10
(0.05, 0.20)
0.21
0.13
1.08
1.25
0.72
1.26
0.19
0.10
(0.10, 0.05)
0.38
0.12
1.13
1.38
1.21
3.41
0.37
0.10
(0.10, 0.10)
0.30
0.13
1.03
1.31
0.84
1.04
0.32
0.10
(0.10, 0.15)
0.21
0.15
0.83
1.15
0.80
3.12
0.27
0.14
(0.10, 0.20)
0.25
0.17
1.23
1.28
0.93
3.47
0.30
0.16
(0.15, 0.05)
0.45
0.15
0.88
1.35
0.99
0.99
0.51
0.09
(0.15, 0.10)
0.38
0.16
1.21
1.21
0.91
1.02
0.45
0.12
(0.15, 0.15)
0.42
0.16
0.89
1.29
1.18
1.04
0.51
0.13
(0.15, 0.20)
0.63
0.21
1.13
1.21
0.89
1.02
0.75
0.20
(0.20, 0.05)
0.52
0.17
1.13
1.29
1.11
0.98
0.53
0.10
(0.20, 0.10)
0.49
0.19
0.83
1.10
0.83
0.96
0.55
0.15
(0.20, 0.15)
0.83
0.28
1.10
1.15
0.97
0.99
0.81
0.24
(0.20, 0.20)
1.48
0.40
1.25
1.08
1.23
1.02
1.17
0.39
Esta conclusão pode ser vista na forma analı́tica da expressão (3.2). Por exemplo,
numa situação extrema, suponha que a população esteja totalmente espalhada numa
região, dessa forma é razoável supor que todas as redes existentes (vazias e não vazias)
sejam de tamanho 1, ou seja, Z1 = · · · = ZN −X+R = 1. Neste caso, o número de redes
36
não vazias passa a ser o número de unidades não vazias na população, ou seja, R = X.
PN −X+R
Portanto, para todo j = 1, . . . , m, Zij = 1, gij ,j = N − (j − 1),
Zi = N e
i=1
Pj−1
k=0 Zik = m − [m − (j − 1)]. Portanto, a probabilidade de seleção em (3.2) se reduz a:
[s | C, R, X] =
1×N
1 × (N − 1)
1 × [N − (m − 1)]
×
× ··· ×
= 1,
N −0
N −1
N − (m − 1)
para qualquer amostra s sorteada.
Logo, a probabilidade de inclusão da amostra
permanece inalterada para qualquer amostra s selecionada desta população.
3.1.3
Estudo simulado com população real
A fim de ilustrar a eficiência do modelo em (3.1), será feita a seguir uma comparação
do estimador obtido do ajuste de tal modelo com o estimador de Horvitz-Thompson
modificado, dado em (2.4), obtido com base no desenho amostral adaptativo por
conglomerados. Além disso, ambos serão comparados à amostragem aleatória simples
sem reposição. Esta ilustração será feita a partir de sorteios de repetidas amostras de
uma população verdadeira. Tal população constitui-se de marrecos da asa azul na região
da Flórida, nos Estados Unidos, no ano de 1992. Em particular, esta é uma espécie
rara de aves aquáticas com um comportamento agrupado. Esta mesma população e
outras duas espécies, as quais apresentam diferentes graus de agrupamento, foram usadas
para comparação da eficiência da amostragem adaptativa com relação a outros planos
amostrais em Smith et al. (1995).
A Figura 3.2 corresponde à área de estudo, dada em Smith et al. (1995), a qual foi
subdividida em N = 200 unidades de uma grade regular, tal que cada unidade apresenta
o número de indivı́duos da população de marrecos da asa azul naquela região. Observe
que esta população caracteriza-se com um aspecto raro e extremamente agrupado.
Além disso, usando as expressões em (2.5) e (2.6) para avaliar numericamente estas
propriedades na população, obteve-se P R = 0.11 e V IR = 0.71, o que também indica
que a população em estudo tem estas caracterı́sticas, justificando assim o uso do plano
amostral adaptativo.
Para avaliar a eficiência dos métodos de amostragem citados, para esta particular
população, foram sorteadas 100 amostras e para cada amostra obtivemos uma estimativa
37
1
7144 6399
103 150
6
10
2
60
122 114
14
3
2
3
12
2
4
20
5
2
3
Figura 3.2: População real de marrecos da asa azul na região da Flórida, nos Estados
Unidos, no ano de 1992, disposta numa grade regular de tamanho N = 200.
do total populacional T . Tal estimativa foi obtida com base no estimador não viesado
para o total sob os planos adaptativo e aleatória simples, e no caso do ajuste do modelo
Bayesiano em (3.1) são obtidas amostras da distribuição a posteriori, e tal estimativa
pontual é dada pela média a posteriori de T .
Em cada uma das 100 amostras, sorteia-se aleatoriamente e sem reposição n1 unidades
iniciais na grade e, se pelo menos um marreco da asa azul é observado, as unidades
vizinhas, ou seja, as de lado contı́guo, são incluı́das na amostra, e o procedimento é
repetido até o momento em que uma unidade de borda, ou seja, sem qualquer marreco
de asa azul, é obtida. Dessa forma, cada amostra adaptativa possui n unidades divididas
em m redes (m ≤ n1 ). E com base nestas n unidades, estimamos o total populacional
a partir do estimador em (2.4) e no modelo (3.1). Além disso, também foram obtidas
estimativas para T considerando amostras aleatórias simples de tamanho n, com base no
estimador T̂AAS = N ȳ.
A mesma distribuição a priori descrita anteriormente para o modelo (3.1) foi utilizada
neste estudo, exceto a distribuição de γ, para o qual foram usados aγ = 5 e bγ =
2, como recomendado em Rapley e Welsh (2008) para a maioria dos casos. Notouse que ao atribuir distribuições para γ com alta massa de probabilidade em valores
38
maiores, surgiram problemas de superestimação do total populacional, devido às amostras
coletadas conterem na sua maioria a rede de maior tamanho, a qual apresenta maiores
valores de Y , diferente dos dados artificiais que eram gerados de um modelo que supõe
homogeneidade entre as unidades.
Na Tabela 3.2 temos a eficiência de cada estimador para alguns tamanhos de amostra
iniciais. A eficiência de um estimador é dada pela razão entre as variâncias para cada
estimador em questão, logo se esta razão é maior que 1 significa que, em termos de
precisão, o estimador do denominador é mais eficiente do que o outro. Em particular,
AAS
defina, ef(T̂HT
∗ ) a eficiência do estimador da amostragem aleatória simples com relação ao
estimador de Horvitz-Thompson modificado descrito pela expressão em (2.4), ef(T̂BAAS )
a eficiência do estimador da amostragem aleatória simples com relação ao estimador
∗
Bayesiano e ef(T̂BHT ) denota a eficiência do estimador de Horvitz-Thompson modificado
com relação ao estimador Bayesiano. Além disso, E(n) denota o valor esperado do
tamanho final da amostra adaptativa utilizando as 100 amostras geradas, portanto, é o
tamanho médio das amostras aleatórias simples selecionadas para a comparação.
Observe que, para qualquer tamanho de amostra, as duas abordagens que usam o
plano amostral adaptativo são mais eficientes que a amostragem aleatória simples. Exceto
para n1 = 4, em que a conclusão se inverte quando compara-se T̂AAS com relação a
T̂HT ∗ . Quando comparados entre si, o modelo em (3.1) apresenta maior eficiência que a
estimação com base no desenho amostral adaptativo.
Portanto, conclui-se que o modelo (3.1) é eficiente e apresenta vantagens quando
comparado com as outras metodologias. Com base nesta conclusão, o interesse agora
é estender este modelo para outros contextos usuais. Na próxima seção é proposta
uma extensão do modelo (3.1) para populações que apresentam constante mobilidade,
incorporando esta caracterı́stica ao próprio modelo.
Vale ressaltar que um modelo inflacionado de zeros poderia ser uma alternativa para
previsão nestas populações raras, devido ao excesso de zeros. Esta classe de modelos
ganhou destaque com Lambert (1992). A ideia geral desta classe de modelos é baseada
na inclusão de massa de probabilidade no ponto zero, inflacionando suas possibilidades
de existir no modelo, por meio de uma mistura de distribuições. No entanto, neste
39
Tabela 3.2: Estudo simulado com a população de marrecos da asa azul: eficiência
relativa para o estimador do total populacional com base no desenho amostral adaptativo
(estimador de Horvitz-Thompson modificado) e no ajuste do modelo (3.1), com relação à
amostragem aleatória simples de tamanho n. A eficiência do estimador Bayesiano com
relação ao estimador de Horvitz-Thompson também é apresentada na última coluna.
∗
AAS
ef(T̂HT
ef(T̂BAAS ) ef(T̂BHT )
∗ )
n1
E(n)
4
16.74
0.44
14.37
33.33
10
25.23
1.68
12.36
7.14
20
39.91
2.60
7.12
2.70
40
66.63
3.19
4.30
1.35
trabalho o objetivo é fazer previsão acerca de uma população dividida em redes, as quais
por definição são unidades não vazias, portanto não é contemplada a possibilidade de ser
zero. A amostragem adaptativa por conglomerados é portanto uma abordage, totalmente
cabı́vel a esta situação e não fornece informações sobre as unidades vazias, apenas sobre
as não vazias. Por isso o modelo de Rapley e Welsh (2008) é formulado apenas para as
redes não vazias.
3.2
Um
modelo
para
populações
móveis,
em
crescimento ou decrescimento
A biosfera está constituı́da de sistemas que mudam com o passar do tempo. O modo
pelo qual o sistema muda depende de sua organização e dos recursos disponı́veis a ele. Por
exemplo, alguns ecossistemas aumentam em tamanho e complexidade, enquanto outros
detêm seu crescimento. O estudo da dinâmica das populações naturais é importante para
compreender o que ocorre nos ecossistemas em equilı́brio. Este tipo de comportamento,
em geral, é observado em populações de animais, habitats ou outra espécie sensı́vel a
mudanças.
40
Neste caso, o mais comum é trabalhar com modelos espaço-temporais, mas quando
trata-se de populações raras e agrupadas podemos ter grandes dificuldades em ajustar tais
modelos comumente vistos na literatura, principalmente se a elaboração do planejamento
amostral não levar este fator em consideração na coleta dos dados.
McDonald
(2004) apresenta estudos por amostragem que resultaram em estimativas altamente
imprecisas simplesmente pelo fato do pesquisador em curto intervalo de tempo “perder”
a população-alvo, devido ao grande poder de deslocamento, mortes, entre outros fatores.
Inclusive, o próprio procedimento de coleta dos dados pode ser um fator gerador de
dispersão da população de interesse. Uma opção para estes cenários é a replicação
da coleta de dados ao longo de um perı́odo de tempo, com o objetivo de ganhar mais
informações sobre este comportamento móvel, difı́cil de ser estudado. Dessa forma, além
de gerar estimativas mais precisas, tal abordagem pode ser altamente relevante para
possibilitar possı́veis intervenções mais precisas no futuro neste tipo de população, em
casos de epidemia, por exemplo.
O objetivo desta seção é propor para situações como as descritas acima, modelos
de previsão que incorporem o plano de amostragem adaptativa, mas que leve em conta
não só a raridade e esparsidade geográfica, como o modelo proposto por Rapley e Welsh
(2008), mas que também levem em conta a mobilidade da população ao longo de um
perı́odo de tempo.
3.2.1
Amostragem adaptativa para populações móveis
Um comportamento de mobilidade, crescimento ou decrescimento em um espaço ao
longo de um perı́odo de tempo é comumente visto em populações biológicas. Em geral,
esta caracterı́stica é algo natural da espécie em estudo, ou pode simplesmente surgir
num estudo por levantamentos estatı́sticos, pelo fato do método de amostragem utilizado
alterar seu habitat natural, incentivando esta dinâmica populacional.
Por outro lado, estas populações biológicas, por exemplo, também em geral são uma
fração pequena da população e estão distribuı́das numa região em grupos. Já foi visto que,
para populações com tais comportamentos, a amostragem adaptativa por conglomerados
pode ser bastante eficiente quando comparada a outros planos mais comuns e menos
41
custosos. Mas, segundo McDonald (2004), se a população, além destas caracterı́sticas,
tem alta mobilidade por fatores naturais, ou se move ou se destrói na coleta dos dados,
adaptações neste planejamento devem ser realizadas. O mesmo ocorre se a população de
interesse tende a crescer, indicando situações de alastramento.
McDonald (2004) apresenta algumas alternativas para o problema da mobilidade,
tais como: redefinir a vizinhança de forma que não inclua somente unidades de lado
contı́guo e a criação de um ı́ndice de presença de espécies, que não seja a observação
direta. Este último recai na amostragem adaptativa dupla, proposta por Felix-Medina
e Thompson (2004) e apresentada na Seção 2.2.4. Um exemplo desta é um estudo de
monitoramento da abundância de gambás na Nova Zelândia. Para detectar a região
de interesse, são colocados de forma adaptativa blocos de cera com algum atrator e a
frequência de mordidas neste bloco é um indicador da distribuição de gambás na região.
Em seguida, uma subamostra desta amostra adaptativa é observada nestes locais a fim
de obter uma estimativa do total de gambás na região.
Por outro lado, sob o ponto de vista de inferência baseada em modelos de
superpopulação, o modelo (3.1), proposto por Rapley e Welsh (2008) e apresentado
na seção anterior, não se ajusta explicitamente a populações com esta dinâmica. A
princı́pio, para inferência num único instante de tempo, as alternativas descritas acima
e apresentadas por McDonald (2004) podem ser facilmente inseridas na função de
verossimilhança do modelo, com mudanças somente na definição da vizinhança e redes.
Uma outra alternativa, que pode gerar estimativas ainda mais confiáveis é a coleta de
dados ao longo de um perı́odo de tempo e uso destas amostras repetidas para inferir sobre
os parâmetros populacionais. Esta abordagem pode ser útil também para o entendimento
do comportamento elusivo da população em perı́odos de tempo, além de previsão para
tempos futuros. Neste caso, para cada tempo terı́amos uma amostra coletada, e para cada
tempo terı́amos uma estimativa calculada com base no estimador de Horvitz-Thompson
dado em (2.4), por exemplo. No caso da abordagem baseada em modelo, o modelo (3.1)
seria ajustado para cada tempo de forma independente. E poucas são as alternativas na
literatura para dados deste tipo. Em particular, temos interesse em estender o modelo
42
(3.1), proposto por Rapley e Welsh (2008), incorporando este comportamento móvel para
que se ajuste a populações deste tipo.
3.2.2
Incorporando estrutura de crescimento e decrescimento
ao modelo
Como o objetivo é propor um modelo para previsão em populações que evoluem
dinamicamente, de forma que a amostragem adaptativa por conglomerados ainda seja
um plano amostral eficiente, serão tratadas apenas situações em que este crescimento
se dá em sua maior parte dentro das redes, de forma a não descaracterizar a raridade
e agrupamento da população, os quais são os principais motivos para o uso deste plano
amostral.
Na Figura 3.3 temos uma ilustração da dinâmica de uma população artificial
sobreposta a uma grade regular com N = 400 unidades. Para gerar esta população
foi utilizado o processo pontual conglomerado de Poisson (ver Diggle et al. (1983)), o
qual gera configurações de eventos agregados, onde os conglomerados são interpretados
como grafos e, portanto, formados por pais e filhos. Em particular, fixou-se o número de
redes e de observações em cada rede na geração. Dessa forma, dado o número R de redes
não-vazias, as coordenadas dos centroides (pais) destas R redes (grafos) são sorteadas de
uma distribuição Uniforme definida neste espaço. A partir destas R localizações, com
o número de observações Yi (filhos), para cada rede i, i = 1, . . . , R, as localizações dos
Yi − 1 filhos são gerados para cada rede de uma distribuição Normal com média nas
coordenadas dos pais e variância fixada. O número Yi para cada rede i foi gerado de uma
distribuição Poisson. Como o objetivo era apenas ilustrar uma população dinâmica de
interesse, observe que para tal ilustração não foram necessárias as variáveis número de
células não-vazias e número de células em cada rede que fazem parte do modelo (3.1), pois
o processo utilizado na geração é um processo pontual, e processos deste tipo independem
da divisão da área, no caso da grade regular, o que não comprometeu de forma alguma
a ilustração.
43
Observe na Figura 3.3 que ao longo do tempo o número de unidades com observações
aumenta e o número de redes varia de forma estável.
●
●
●
●
● ●
●
● ●
●
●● ● ● ●
●
●●
●
●●
●●●
●
●
●●●●●
●
●
●●
● ●
● ● ●●
● ● ● ●
●
●● ●
●●
●●
●
● ●
● ●
●
●
●
●
●
●
●
●
● ●
●
●
●
●●
●
●●
●●
●
●●
●
●
● ●● ●
●
●●● ● ●
●●
●
●
● ● ●●
●● ●●
●
●
●
●
●● ●●
●
●● ●●
●●●●
●●● ●●
●
●●
●
●●
●● ●
● ● ●● ● ●
● ●● ●
●
●●●
●
●●● ●●
●●
●●
●
●
●●●
●●
●● ●●●● ●
● ●● ●
●●
●
●
● ●
●● ●
●
●●●●●●● ● ●
● ●
●● ●● ●
●
●
●● ● ●● ● ●
●
● ●
●
● ●●
●
●
●
●
●●
●
●
●● ●● ●● ● ● ●
●
●●
●●
●
● ●
● ●●
●
●
●●
●
●
●
●
●
●
●
●
●
●
● ●
●
●● ● ●● ●● ●
●● ●
● ●● ●●
● ● ●
● ●●
●
●●●● ●●● ●
● ●●●
●●
●
●
●
●● ●
●●●
●
●
●
●●
● ●
●●
●● ●
●● ●
●●
●
●
●●
●
●
●●●●
● ●●
● ●●●
●●
●●● ●
●
● ●● ●
●●
●
●
●
●●
●
●
●
●
●
●
●
●
●●
● ●
●
●
●
●
●
●
● ●●
●●
●
●
●●
●
●
●● ● ●
●● ● ●
● ●●
●
●
●
● ● ●● ● ●
●
●
● ●
●●
● ●●
●
●● ●
● ●
●
t=1
●
●
●
● ●
●
●
●
● ●●● ●
●
●● ● ●
●● ● ●
● ● ● ●●
●●● ●
●
●●●
●
● ●●
●●
●●● ●
● ●●
●●
● ●
●
●●●
●●●● ●●● ●
●
● ●
●
●
●
●● ●
●
●
●●
●
●
●
●
● ●● ● ● ●
●
●● ● ● ●
●
●●
● ●
●●● ● ●
●●
● ● ●
● ●●
● ● ●
● ●● ● ●●
●●● ● ●●
●●
● ●●
●
● ● ● ● ●●
●●●●
●
●
●
●●
● ●●● ●
●
●
● ●●
● ●
●● ●
●● ●
● ● ●
●●● ●● ●
●
●
●●
●● ●●● ●
●●
●
●
●
●●
● ●
●
●
●
●
●●
●
●
●
●
●
t=2
●
●
●
●
●
● ● ●
●
●
●
●
●●
● ●●
● ●
●●
●●
●●
● ●●
●
● ●
●● ●●
●● ●
●● ●●●
●● ●●
● ●● ●●●●
●●● ● ●
●
●●
●
●●●●●
● ● ●
●●●●
●
●● ●●
● ●
●●●●●●
●
●●●
●●
●●
●●
●● ●●
●
● ●
●●●
●●
●●●
● ●● ●
●● ●
●
●●
●
●
●
●●
●
●
●
● ●
●●
●●●
●●
●●
●
●●
●
●
●
●
●●
●
●●
●
●●●
●●
●
●
● ●●●
●●
● ● ●● ●●●●
●
●
●
●●
●
●●●
●
●
●●
●●
●
●●
●●
●
● ●●
●
●●
●●
●●
●●
●●
●
●
●●●●
●
●●
●
●●
● ●●●
●● ●●●
●
●●
●
●
●●
●●●● ●
●● ● ●●
●● ●
●
●●
●
●●●
●●●
●
● ●●●● ●
●● ●
●●
● ●●
●● ●● ●
●
● ●
●
●
●●
●
●
●● ●
●
● ● ●
●● ●
●
●
●
●
●● ● ● ●
●
●●
●● ●
●
●
●● ● ● ●
● ●●
●● ● ● ●
● ●● ● ● ●
●
●
●
●
●
● ●
●● ●
● ●
●●
●●
●●
●
●
●
●
●● ●●
●● ●
●●●
●● ●
● ● ● ●●●
●●● ● ● ●
●
●
●●●●
●
●●
●● ●
●●●
●
●●● ●●
●●
●
●●●
●
●
●●
●
● ● ●
●● ● ● ●
● ●
●●
● ●● ●●
●●●●●
●●
●
●●●
●●
●
●●●
●
●●
●●
● ● ●●●●●
●
● ●
●●●
●●
●●●
●
●
●
●
●
●
●
●●
●●
●
●
●
● ●●●
● ●●
●●
●
●● ●●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●● ● ● ●
●
● ●
●●●
●
● ●
●
●●●●●
●
●
● ●●●
● ● ●●●
●
●●
●
●●
●●
●●
●
●
● ●●●
●●
●
● ●
●
●
● ●●
●●● ●●
●●
●
●●
●
●
●
●●
●
●
●● ●●●●
●●
●●
●
●
●
●●
●
●
●●● ●●● ●
●
●
●
●● ● ●
●
●●● ● ●
●
● ●●
● ●●
●● ● ●
●
●
●● ● ●● ● ●
●
●
●
●●● ●
●●
●
●
●
●
●
●
●
● ●●
● ● ●
●
●
●
●
●
●
●● ● ●
●
●
● ●
●
●
●● ●
●
●
●● ● ●
●
●
●●
●●
●
● ●
●
●
●● ●● ●●
● ●
●
●
●
●●●●● ●
●
●
●
●
● ●●●● ●
●
●●
●●
●
●
●
● ●
●●
●● ●
●●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
●
●● ● ●● ●● ●
●● ●
● ●● ●●
● ● ●
● ●●
●
●●●● ●●● ●
● ●●●
●●
●
●
●
●● ●
●●●
●
●
●
●●
● ●
●●
●● ●
●● ●
●●
●
●
●●
●
●
●●●●
● ●●
● ●●●
●●
●●● ●
●
● ●● ●
●●
●
●
●
●●
●
●
●
●
●
●
●
● ●● ●●
●
●
●
●
●● ● ● ● ●
●●
●●●● ●● ●
● ●
●
●●●
●●
●●
●
●●
●●●● ●
●
●
●●
●
●
●●
● ●●●
●
●
●
●
●
●●
●●
●
●
●●
●
●
●
●●
●
● ●● ●
●●
●
●●
● ●● ●● ●
●
●
●
● ●●
●● ●
●●
●
●
●
●
●●●●
●
●
●●●●●●●●
●
●
●●
●
●
●
●●●●
●●
●
●
●
●●
●●
●
●
●●
●
●
●
●
●● ● ●
●
●●
● ●●
●
●
●
●●
●
●
●
●●
●
●
●
●●●
●● ● ●
●
●●●●●●
●●
●
●●
●●
●
●
●●●
●
●
●
●●
●●●●● ●
●
●
●●
●●●●
● ●●
●
●●
●●
●
●●● ●●
●●● ●●
●●● ●●
●
●
●●●
●●●
●
● ●●
●●●●
●
●
●
●
●
●●
●
●●
● ●
●
●
●
● ●
●
● ● ●
●
●
●
●
●
●
●
●
●●
●●
●
●
●
● ●
● ● ● ●● ●
●
●● ●
● ●
●
●
● ● ●●
●
●
● ●
● ● ●●
●
● ●●● ●● ●
● ●
● ● ●●●
● ● ●● ●
●
● ● ● ●●
●●
●●
●
●
● ●●●●
●
●●●
●
●●● ●●● ●
●
●
●
●●● ●
● ●
● ●● ● ●●
●●●
●●
●
●
●●● ●
●●
●
●
●●●
●
●
● ●●
● ●
●● ●
●
●
●
●
●●
●● ●
●●●●
●
● ●
●●
●
● ●●●●●
● ●
●
●
●
●●
●
●
●●●●● ●
●
●
●
●
●●● ●
●● ● ●
●
● ● ●●
● ● ●
●
●
●
●
●
●
●
●
●
●
●
● ●●
●● ●● ●
●●● ● ●
●●
● ● ● ●●●
● ●
●●
● ●●● ● ●●●
●
●
●●● ● ●●
●●
● ●● ●
●
●
●
●● ●●
● ●
●
●●
●●●
●
●●● ●● ●
●●●
●
●
●
●
● ●●● ●
●
●
●●
●●
●●
●
●● ●
●
●●● ●
●
●●
● ● ●
●●● ●
● ● ●●●
●
● ●
● ●●●
●●
●●
● ●● ●● ● ●
●●●
●●
●
●●● ●● ● ●
●●
● ●●● ● ●
● ●●
● ●
●●
● ●
● ●
●
●
●● ●
● ●
●
●
●●●
●
●
●
● ●●
●
●
●
●
●
●
● ●
●
●
● ●
● ●● ●●
●
●●●
●●●
●●
●
●●
●●● ● ●● ●●
●●●● ●
●
●●
●
● ●● ●●
●
●●
●
●●
●
●
●
●
● ●●
●
●●
●
●
●●●●
●●●●
●
●
●
●●
●
●● ●●●
●●
●
●
● ●●
●
●●
●
●
●
●
●
●
●
●
●●
●●
●● ●●●
●
●
●
●●
●
●●
●
●
●
●●
●●
●●
●●
●
●
●
●●
●●
●
●
●
●●
●
●●●●●●●
●●
●
●●
●
●
●
●
●
●
●●●
●●
●
●
●
●
●
●
●●
●●●
●●
●●
●●
●
● ●
● ●● ●
●
●
●
●●
●
●
●
●
●●
●
●
●● ●
●●
●●
●●
●●
●
●●
●
●
●●●
●
●
●
●
●
●
●
●
● ● ●●
●
●● ●
●
●● ●●
●● ●
● ●●●
● ●●
●
● ●●
● ●●
● ● ●
●●
● ●
● ●●
●●
●● ●
● ●
●
●●
●
●
●
●
●
●
t=3
t=4
●
●
● ●
●●
●● ●
●
●
● ●●
●●
●
●
●
●
●
● ● ●
●●
● ●●
● ●●●●
●● ●
●
●
●●● ●●●
●●
●
●● ●
●
●
●●
●
●●●
●
● ●
●●
●●
●
●
●●
●●●
●
●●
●●
●●
●
●
●●●●●
●
●
●●● ●
●●
●●
●
● ●●
●
●
● ●
●
●●
●
●
●
●●
●●● ●
● ● ●
●●
●●●
●● ●
●
●●
●
●
●
● ●
●●●
●
●●●●●● ●
●●
● ●
●●
● ●
● ●●
●
●
●
●
●
● ●
●
●
●●
● ●
●
● ●●
●
● ●●
● ●●● ● ●
●
● ● ●●
●
● ● ●●
● ●
●
● ●
●●●
●●
●●● ●●●●
●●
●●●●●●
●
●●●
●
●
●
●●●
●
●●●●●●●●●
●●●
●
●
●●
●
●
●
●
●
●
●
●
●
●●
●●
●
●
●●●
●●●
●●
●
●
●●
●
●
●
●
●●
●
●●
●
●
●
●
●
●
●●
●●
●●●
●
●●●
●●
●
●
●
●
●●
●
●● ●
●●● ●
●●
●●
●
● ●
●●
● ●●●● ●
●●●
●● ● ●
●
●
●●●●●● ●●●
●
●
● ●●
●●
●
●
●
●
●
● ●
●
●
●
●
●
● ●●
●
● ●● ●●
● ●●
●● ●
●
●
●●
●●
● ●●●●
●● ●
●●
●
●●
●
●
●●
●
●
●●
● ●●
●
●●
●●
●
●
●
●
●●●●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●●
●
●● ● ●● ●●
●●
●
●
●
●●● ●
●
●
●
●
●
●●
●
●
●
●
●
●● ●● ●
●●
●
●
●
●
●
●
●
●●
●
●●
●
● ●●
●● ●
●
●
●
●
●●●●●●
●●
●
●
●
●
●
●●●
● ●●
●●
●
●
●
●
●●
●
● ●● ●
●
●
●●
●
●
●
●●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●● ●●
●
●
●
●
●
●●
●
●
●●
●
●
●
●●
●●
●●
● ●● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●●
●
● ●●
● ●●●
●●
●●
●
●
●
●●●
●●●●
●●● ●●●
●●●● ●●
●●
●
●●
●●
●
● ●●●
●●●●
●
●
●
● ●
●
● ●
●●
●
●
●
●
●● ●
●
●●
●●
●●
●
●
●●
●
● ●
● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●● ●
●●
● ●●
●●
●
●● ●
● ●●
●
●●
●
●
●●
●●
●
●
●
●
●●
●● ●● ●
●
●
●●
●● ●●
●
●●●
●●
●
●
●
● ●●
●
●
●
●
●●
●
● ●●
●●
●
●
●
●
●
●
●
●
●●●
●●
●
●
●
●
●
●
●
●●
●
●
●●●●
●
●●
●
●●
●
●
●●●
●
●
●
●
●
●●
●●
● ●
●●●
●
●
● ●●●●
●
●
●
●
●
●
●
●
●
●
●●●
●●●
●
●
●
●●
●
●
●
●
●
●
●●
●●
●●
● ● ●●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●●
●●
●●
●
●
●●
●●
●
●●
●
●
●
●
●●
●●●
●●
●
●
●
●
●
●
●
●
●
●●
●
●
●●●
●
●●
●
●
●●
●
●
●
●●
●
●
●●
●
●●
●
●●
● ●●
●
●
●
●
●
● ●●
●
●
●●
●
●
●
●●
●● ● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●●
●
● ●● ●
●●
● ●●
●●
●
● ●
●
●
●
●
●
●● ●
● ●
● ●
●● ● ●●
● ●
● ● ●
●
● ●●
●
●
●
t=5
Figura 3.3: Ilustração da evolução dinâmica de interesse de uma população rara e
agrupada numa região sobreposta a uma grade regular com N = 400 unidades.
44
Note que a partir do modelo em (3.1) é possı́vel incorporar esta dinâmica populacional
acrescentando alguma estrutura temporal aos parâmetros α e β.
Se tornarmos o
parâmetro α dinâmico, deixando β fixo, o número de unidades não-vazias na população
se altera ao longo do tempo e, portanto, pela estrutura de condicionamento do modelo,
o número de redes também pode se alterar.
Se for feito o contrário, ou seja tornar o parâmetro β dinâmico, deixando α fixo,
teremos uma população cujo número médio de unidades não-vazias não se altera ao longo
do tempo, mas sua disposição dentro das redes sim, o que pode criar novas redes com
o número de unidades reduzido, ou ainda desaparecer redes com o número de unidades
crescendo dentro de algumas redes.
Uma outra possibilidade intuitiva é a incorporação de dinâmica nos dois parâmetros α
e β ao mesmo tempo, isso geraria uma população menos estável que em qualquer um dos
dois cenários citados anteriormente. Isto porque estarı́amos alterando diretamente tanto
o número de unidades na população, quanto o número de redes. Note que a estrutura
dinâmica a ser imposta deve ser de forma controlada, a fim de que a população-alvo rara
e agrupada não se descaracterize ao longo do tempo.
Neste trabalho temos particular interesse na primeira extensão, onde cresce o número
total de unidades na população ao longo do tempo e, por conta do condicionamento do
modelo, o número de redes varia, mas de forma mais estável que as outras duas opções.
Dessa forma, serão contemplados comportamentos de mobilidade caracterizados pelo
surgimento e desaparecimento de redes, e ainda pelo crescimento ou decrescimento do
número de observações nas redes que permanecem, mas com uma estabilização no final
do tempo. Portanto, não serão considerados cenários com alastramento desordenado
ou desaparecimento global da observação de interesse, como uma epidemia, no caso de
doença, por exemplo.
3.2.3
Modelo de crescimento exponencial
Como o interesse está em modelar populações que apresentam um crescimento ou
decrescimento médio de observações dentro das redes, mas com uma estabilização ao longo
do tempo, em particular, modelos de crescimento exponencial podem gerar populações
45
com esta estrutura, além de serem amplamente utilizados em problemas reais em diversas
áreas, como na ecologia. Considere que as observações obtidas a partir de um processo Yt
ao longo de perı́odos de tempo t = 1, . . . , L são modeladas a partir de uma distribuição
de probabilidade na famı́lia exponencial, tal que E(Yt | θ t ) = λt , onde θ t é um vetor de
parâmetros. Modelos caracterizados pela parametrização θ t = (a, b, c)0 e por uma função
de ligação h tal que
h(λt ) = a + b exp(ct) e


λφt , se φ 6= 0
h(λt ) =
 log(λ ), se φ ≈ 0
t
são chamados modelos de crescimento exponencial generalizados e podem ser vistos com
detalhes em Migon e Gamerman (2006).
O parâmetro c está relacionado com a velocidade de crescimento/decrescimento (ou
curvatura), o parâmetro b com a intensidade do crescimento/decrescimento e a com a
localização da curva. Derivando a expressão a+b exp(ct) em relação a t, podemos concluir
que a curva será crescente se b e c tiverem o mesmo sinal e decrescente caso contrário.
Pela derivada segunda, podemos concluir que a curva tem concavidade voltada para cima
se b > 0 e para baixo se b < 0. Vale notar ainda que se c < 0 então a curva tem um
comportamento não explosivo, convergindo para a quando t → ∞.
A principal vantagem em utilizar estes modelos é a possibilidade de manter as
medições de Yt na escala original, transformando apenas a trajetória de Yt , o que torna
a interpretação dos resultados mais simples. Além disso, os intervalos de tempo não
precisam ser igualmente espaçados, permitindo que se trabalhe com dados provenientes
de pesquisas com datas de referência distintas através de uma codificação do ı́ndice t de
tempo.
A proposta, portanto, é modelar o parâmetro α do modelo (3.1) a partir de uma
curva de crescimento exponencial. Como este parâmetro é uma probabilidade, é natural
usar-se na modelagem uma função de ligação logı́stica, portanto o modelo apresenta-se
da seguinte forma:
logit(αt ) = log
αt
1 − αt
= a + b exp(ct), t = 1, . . . , L.
46
Em particular, ao modelar o parâmetro α desta forma, os possı́veis valores que os
parâmetros da curva exponencial a, b e c podem assumir devem estar compatı́veis com
o contexto de populações raras e agrupadas e, portanto, serão usados os resultados
apresentados na Seção 3.1 para efetuar esta escolha.
O interesse neste trabalho
concentrar-se-á em dois casos: (i) crescimento da população e estabilização com a
evolução do tempo; (ii) população decrescente ao longo do tempo e estabilização, de
forma que a população não desapareça. Na Figura 3.4 estão as duas curvas de crescimento
que caracterizam os dois cenários de interesse neste trabalho, para L = 50. Para obter
a curva crescente (a) assumiu-se a = −1.73, b = −1.41 e c = −0.15, o que resulta no
parâmetro αt iniciando em 0.05 e estabilizando em 0.15, os quais são valores razoáveis
para este parâmetro no modelo (3.1) no contexto de populações raras e agrupadas. Já
a curva (b) é obtida para a = −2.20, b = 0.94 e c = −0.15, produzindo uma curva
decrescente para αt iniciando em 0.20 e estabilizando em 0.10.
0.18
●
●
●
●
0.14
αt
0.10
αt
0.14
●●●●●●●●●●●●●●●●
●●●●●●●●●
●●●●
●●●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●●
●●●
●●●●
●●●●●●●●●●●●
●●●●●●●●●●●●●●●●●
0.10
0.06
●
●
●
●
0
10
20
30
40
50
0
t
10
20
30
40
50
t
(a) Crescimento
(b) Decrescimeno
Figura 3.4: Curvas de crescimento e decrescimento de interesse para αt , t = 1, . . . , 50.
Em (a) fixou-se a = −1.73, b = −1.41 e c = −0.15, e em (b) a = −2.20, b = 0.94
e c = −0.15, o que resulta no parâmetro αt variando de 0.05 e 0.15 e de 0.2 a 0.1,
respectivamente.
Desta forma, uma extensão do modelo (3.1) para populações que evoluem ao longo do
tempo com uma dinâmica semelhante à descrita anteriormente é dada, para t = 1, . . . , L,
por:
47
Yit∗ | Cit , Rt , γ ∼ Poisson Truncada independente (γCit ), Yit∗ ≥ Cit , i = 1, . . . , Rt ,
Rt
X
1
Ct | Xt , Rt ∼ 1Rt + Multi Xt − Rt , 1Rt , Cit = 1, . . . , Xt − Rt + 1,
Cit = Xt ,
Rt
i=1
Rt | Xt , β ∼ Binomial Truncada (Xt , β), Rt = 1, . . . , Xt ,
(3.4)
Xt | αt ∼ Binomial Truncada (N, αt ), Xt = 1, . . . , N,
logit(αt ) = a + b exp(ct), t = 1, . . . , L,
em que Xt é o número de células não vazias no tempo t, Rt é o número de redes não vazias
no tempo t, Ct é o vetor com o número de células não vazias em cada uma das Rt redes
não vazias no tempo t, Yit∗ é o número de observações na rede não vazia i no tempo t,
P t ∗
t = 1, . . . , L. O maior interesse está em prever T = (T1 , . . . , TL )0 , em que Tt = R
i=1 Yit .
Além disso, a variável resposta neste caso segue uma distribuição de Poisson, pois o
interesse concentra-se em dados de contagem, embora seja possı́vel estender esta ideia
para outras distribuições na famı́lia exponencial, assim como (3.1).
Como trata-se de inferência Bayesiana, o modelo (3.4) deve ser completado com a
distribuição a priori para o vetor paramétrico (a, b, c, β, γ). Neste caso, ao ajustar o
modelo, sob distribuições a priori não-informativas para os parâmetros a, b e c, surgiram
problemas de identificabilidade, pois estes parâmetros devem estar restritos a valores de
αt pequenos. Portanto, conclui-se que para o ajuste razoável do modelo (3.4) é necessário
atribuir uma distribuição a priori informativa para estes parâmetros. Por outro lado,
há interesse em permitir ao modelo que verifique se os dados fornecem um cenário de
crescimento ou decrescimento da população, o que é fornecido pelos valores de a e b, para
c negativo, afim de garantir cenários de estabilização com o passar do tempo. Portanto,
supondo independência entre os parâmetros a priori, será atribuı́da a priori para c uma
distribuição Normal truncada nos reais negativos, denotada por c ∼ N(−∞,0) (µc , σc2 ).
Para (a, b) será atribuı́da a priori uma mistura de distribuições normais bivariadas que
contemplam os dois possı́veis cenários, da seguinte forma:
(a, b) ∼ w1 N(µ1 , Σ1 ) + (1 − w1 )N(µ2 , Σ2 ),
48
em que µ1 = (−2.20, 0.94)0 e µ2 = (−1.73, −1.41)0 são os vetores de média para cada
Normal e Σ1 e Σ2 são as matrizes de covariância de cada componente da mistura e
caracterizam o quanto esta distribuição é informativa a priori. Note que a primeira
distribuição da mistura caracteriza a curva de decrescimento e a segunda de crescimento.
Logo, o valor de w1 reflete o peso que será dado às duas distribuições. Neste caso, se não
há informação a priori para dar mais probabilidade de ocorrência a uma das situações,
fixar w1 em 0.5 é uma forma de ser não-informativo e permitir ao modelo que recupere o
comportamento da população ao longo do tempo. Na Figura 3.5 está a distribuição para
os parâmetros (a, b), fixando Σ1 = Σ2 = diag(0.01, 0.01).
Densidade
a
b
Figura 3.5: Distribuição a priori conjunta para o vetor (a, b)0 .
Além disso, supõe-se a priori que β ∼ Beta(aβ , bβ ) e γ ∼ Gama(aγ , bγ ), como no
modelo em (3.1).
Note que os parâmetros a, b e c controlam os valores que αt assumem ao longo
do tempo e, usando o argumento de que as populações-alvo são raras e agrupadas, é
necessário atribuir distribuições a priori para estes informando que αt assume valores
pequenos. O conhecimento a priori mı́nimo para aplicabilidade das técnicas propostas
neste trabalho é que as populações em estudo são raras e agrupadas, portanto a
distribuição a priori deve ser ao menos informativa sobre este comportamento.
Vale lembrar que as variáveis do modelo são compostas por uma parte conhecida e
outra desconhecida como descrito na Seção 3.1, no entanto agora devemos definir esta
partição para cada tempo t, da seguinte forma: Yt∗ = (Ys∗t 0 , Ys̄∗t 0 )0 , Ct = (C0st , C0s̄t )0 ,
49
Zt = (C0t , 10N −Xt )0 , Rt = Rst + Rs̄t e Xt = Xst + Xs̄t , para t = 1, . . . , L. Sejam então
X = (X1 , . . . , XL )0 , R = (R1 , . . . , RL )0 , C = (C01 , . . . , C0L )0 e Y∗ = (Y1∗ 0 , . . . , YL∗ 0 )0 .
Além disso, a cada tempo t uma amostra é selecionada de forma adaptativa e
independente dos outros tempos. Note que o modelo escrito desta forma agregada apenas
nos fornece informações numéricas sobre as redes, portanto não seria possı́vel incorporar
a este modelo um planejamento amostral que não fosse aplicado independentemente
ao longo do tempo, caso contrário necessitarı́amos de informações adicionais, como
localização das unidades que foram selecionadas. Desta forma, a probabilidade de seleção
de uma amostra st = {i1t , . . . , imt } de mt redes no tempo t, t = 1, . . . , L, é dada por:
[st | Xt , Rt , Ct ] =
mt
Y
Zijt × gijt ,jt
,
PN −Xst −Xs̄t +Rst +Rs̄t
P
Z
Zit − jktt−1
i
jt =1
it =1
k
=0
t
onde mt é o número de redes na amostra no tempo t, gijt ,jt é o número de redes de
tamanho Zijt que restam após jt − 1 redes terem sido selecionadas e Zi0 = 0.
Portanto, a função de verossimilhança completa é dada por:
L Y
mt
Y
Zij gij ,jt
PN −Xt +Rt t t Pjt −1
Zit − kt =0 Zikt
it =1
t=1 jt =1




N −Xt
Xt
Xt −Rt
Rt
N
X
 αt (1 − αt )
 t  β (1 − β)
×
×
1 − (1 − αt )N
1 − (1 − β)X
t
Xt
Rt
Cit −1
Rt
Y
1
1
× (Xt − Rt )!
(Ci − 1)! Rt
i=1
[s,X, R, C, Y∗ | a, b, c, β, γ] =
×
Rt
Y
i=1
para logit(αt ) = log
exp{−γCit + Yit∗ log(γCit )}
,
P
it −1
Yit∗ ![1 − C
j=0 exp{−γCit + j log(γCit ) − log(j!)}]
αt
1−αt
= a + b exp(ct), t = 1, . . . , L e s = (s1 , . . . , sL ).
Como a distribuição a posteriori do vetor paramétrico Θ = (X, R, C, Y∗ , a, b, c, β, γ)
não possui forma analı́tica fechada, faz-se necessário o uso de métodos de simulação
estocástica, como o MCMC. Em particular, o amostrador de Gibbs com passos de
Metropolis-Hastings foi utilizado. Além disso, o preditor de T = (T1 , . . . , TL )0 é obtido
da seguinte forma:
T̂t = 10Rst Ys∗t + 10R̂s̄ Ŷs̄∗t ,
t
50
para t = 1, . . . , L cuja amostra da distribuição pode ser obtida via MCMC.
Os passos do algoritmo MCMC são descritos como:
(1) Faça j = 1 e especifique valores iniciais para X1t , Rs̄t , Cs̄t e Ys̄∗t , para t = 1, . . . , L;
(2) Gere a da distribuição [a | X, R, C, Y∗ , b, c, β, γ] = [a | X, b, c];
(3) Gere b da distribuição [b | X, R, C, Y∗ , a, c, β, γ] = [b | X, a, c];
(4) Gere c da distribuição [c | X, R, C, Y∗ , a, b, β, γ] = [c | X, a, b];
(5) Gere β da distribuição [β | X, R, C, Y∗ , a, b, c, γ] = [β | X, R];
(6) Gere γ da distribuição [γ | X, R, C, Y∗ , a, b, c] = [γ | R, C, Y∗ ];
(7) Gere (Xs̄t , Rs̄t , Cs̄t , Ys̄∗t ) da distribuição [Xs̄t , Rs̄t , Cs̄t , Ys̄∗t | Xst , Rst , Cst , Ys∗t ,
a, b, c, β, γ], para t = 1, . . . , L;
(8) Faça j = j + 1 e volte ao passo (2).
Note que este modelo é usado para estimação do total populacional após a observação
de todos os tempos. Se optássemos por uma estimação sequencial dos parâmetros do
modelo e de T à medida que as observações fossem coletadas, algumas dificuldades
poderiam surgir. Os parâmetros do modelo de crescimento exponencial estão associados
ao crescimento e estabilização da população, portanto, só serão realmente bem estimados
após coletadas todas as observações ao longo dos tempos. Logo, recomenda-se que, se o
ajuste de tal modelo for feito de forma sequencial, já tenham sido observados um número
razoável de instantes de tempo, para que o limite do crescimento tenha sido ao menos
atingido.
3.2.4
Estudo simulado
A fim de examinar o desempenho do modelo proposto, foram geradas 100 populações
artificiais com N = 400 unidades ao longo de L = 50 tempos, a partir do modelo
(3.4) com parâmetros (c, β, γ) fixados em (−0.15, 0.10, 10). Para contemplar os dois
tipos de evolução ao longo do tempo de interesse para algumas populações, fixou-se
51
(a, b) = (−1.73, −1.41), o que caracteriza um cenário de crescimento ao longo do tempo,
e para outras (a, b) = (−2.20, 0.94), o que caracteriza um decrescimento. Estes valores
foram escolhidos de modo a representar as mesmas situações descritas na Figura 3.4 e
com cenários de raridade e agrupamento semelhantes aos apresentados na Subseção 3.1.2.
Para cada tempo t, inicialmente 5% das unidades amostrais foram selecionadas
por amostragem aleatória simples sem reposição e uma amostra adaptativa por
conglomerados para cada tempo foi selecionada de forma independente.
A distribuição a priori utilizada para (a, b, c, β, γ) é descrita a seguir. Supondo
independência a priori para cada parâmetro, para β e γ foram utilizadas as mesmas
distribuições descritas na Subseção 3.1.2. Os parâmetros a, b e c controlam os valores
que αt assume ao longo do tempo e, usando o argumento de que as populações-alvo são
raras e agrupadas, será atribuı́da uma distribuição a priori para estes que informe que αt
é um parâmetro com valor pequeno e que a está relacionado com a convergência da curva
e b e c com o crescimento. Além disso, note que a + b é igual ao valor inicial da curva
de crescimento. Portanto, foram atribuı́das distribuições a priori informativas para estes
parâmetros. No entanto, na distribuição a priori para (a, b), usou-se w1 = 0.5, logo a
priori não está sendo informado se, com a evolução do tempo, a população passou por
um crescimento ou decrescimento.
Para verificar o ajuste do modelo foram geradas 200.000 iterações, sendo as 10.000
primeiras descartadas como aquecimento e tomamos amostras de 190 em 190, a fim de
obtermos 1.000 amostras independentes. Nas Figuras 1.6 e 1.7, apresentadas no Apêndice
A, está um sumário da distribuição a posteriori para duas das 100 populações geradas.
Na Figura 1.6 está o resultado a posteriori para uma população que cresce ao londo tempo
e na Figura 1.7 para uma população que decresce. Na Figura 1.6 (a)-(e) e na Figura 1.7
(a)-(e) estão os traços das cadeias da distribuição a posteriori dos parâmetros a, b, c, β e
γ. As Figuras 1.6 e 1.7 (f)-(j) mostram as cadeias para o total em alguns instantes de
tempo arbitrários. A linha em cinza representa o valor verdadeiro usado na geração dos
dados artificiais. Observe que há indı́cios de convergência para todos os parâmetros do
modelo.
52
Como temos amostras das distribuições a posteriori de a, b e c, temos uma amostra
da distribuição de αt , t = 1, . . . , L. Nas Figuras 3.6 (a) e (b) estão em preto a média a
posteriori de αt , em azul o valor verdadeiro e em cinza o intervalo HPD de 95%. Pela
proximidade das linhas azuis e preta e pela linha azul estar contemplada pelo intervalo de
95% é possı́vel concluir que este parâmetro é bem estimado. Finalmente, as Figuras 3.6
(c) e (d) apresentam os valores estimados para o total populacional para os 50 instantes
de tempo. Em preto está a média a posteriori do total para cada tempo, as cruzes em
azul representam os valores verdadeiros do total populacional para cada tempo e em
cinza o intervalo HPD de 95%. Em sua grande maioria os pontos em azul pertencem ao
intervalo, portanto podemos concluir que os totais estão sendo bem estimados para cada
tempo. Note que em ambos os casos o modelo recupera a estrutura ao longo do tempo
dos dados.
Na Tabela 3.3 temos um sumário da distribuição a posteriori dos parâmetros do
modelo de crescimento proposto para as 100 populações geradas. São apresentadas o
EQMR, erro médio absoluto relativo (EAR), a amplitude média dos intervalos HPD de
95% relativizada com relação ao verdadeiro valor e respectiva probabilidade de cobertura.
Note que todos os parâmetros do modelo são bem estimados, pois os respectivos EQMR
e EAR são pequenos. As probabilidades de cobertura dos intervalos na maioria das vezes
se apresentam abaixo do nı́vel nominal desejado de 95%, com exceção dos parâmetros γ,
β e c no caso de população em decrescimento. Uma explicação plausı́vel para este caso é
o fato de que as populações artificiais em decrescimento foram geradas de forma que em
todos os 50 instantes de tempo o número de unidades não vazias estivesse em torno de
10% e 20%, portanto de uma maneira geral estas populações apresentam-se menos raras
e agrupadas que as populações geradas para o cenário de crescimento.
Na Figura 3.7 estão os EQMR, EAR, cobertura e amplitude média relativizada com
relação ao valor verdadeiro dos intervalos HPD de 95% para os totais populacionais para
cada tempo para as populações simuladas de crescimento e decrescimento. Por questões
de melhor visualização gráfica, a amplitude média apresentada é dada pela média dos
valores da amplitude obtida, para cada simulação, dividida pelo valor verdadeiro do total
53
800
600
0.14
400
T
αt
0.10
200
0.06
0
10
20
30
40
50
0
10
20
t
30
40
50
t
(b) T - crescimento
200
0.10
400
αt
0.14
T
600
0.18
800
(a) αt - crescimento
0
10
20
30
40
50
0
t
10
20
30
40
50
t
(c) αt - decrescimento
(d) T - decrescimento
Figura 3.6: Sumário da distribuição a posteriori de αt e do total populacional para uma
população em crescimento e decrescimento ao longo do tempo. Em preto está a média
a posteriori de αt e total populacional Tt , t = 1, . . . , 50, com intervalo HPD de 95% em
cinza e valor verdadeiro em azul.
populacional. Desta forma, é possı́vel uniformizar a escala destes valores para os dois
modelos.
Note que para ambos os cenários, o EQMR e EAR diminuem à medida que a
população torna-se menos rara e agrupada. Portanto, no caso de uma população em
crescimento os erros diminuem com a evolução do tempo e no caso de decrescimento
os erros tendem a ter um ligeiro aumento com o passar do tempo. Com relação à
probabilidade de cobertura dos intervalos HPD de 95% para T nota-se uma subestimação
do nı́vel de 95% para ambos os casos e observando a amplitude média dos intervalos
nota-se um aumento na precisão do intervalo com a evolução do tempo, no caso de
uma população em crescimento, e uma diminuição da mesma, no caso de decrescimento.
54
Tabela 3.3:
Sumário da distribuição a posteriori dos parâmetros do modelo de
crescimento proposto: são apresentados o EQM e EAM, a amplitude média dos intervalos
HPD de 95% e a probabilidade de cobertura para as 100 populações geradas. Os resultados
estão separadas para as populações em crescimento e decrescimento.
parâm
EQMR
EAR
cob
ampl EQMR
Pop. em cresc.
EAR
cob
ampl
Pop. em decresc.
a
0.02
0.10
0.83
0.13
0.02
0.11
0.82
0.10
b
0.01
0.06
0.86
0.21
0.04
0.16
0.85
0.29
c
0.02
0.14
0.80
0.43
0.01
0.08
0.96
0.52
β
0.01
0.06
0.88
0.25
0.01
0.07
0.94
0.30
γ
0
0.01
0.96
0.04
0
0.01
0.94
0.04
Este fato já era esperado e ocorre devido ao aumento e diminuição, respectivamente, do
número de observações com a caracterı́stica de interesse, com a evolução do tempo.
Logo, dada a dificuldade de previsão em populações raras, agrupadas e móveis, no
geral o modelo (3.4) proposto parece ser eficiente para estimar o total populacional para
este tipo de população.
3.2.5
Comparação do modelo de crescimento com outras
abordagens
Outras duas possı́veis abordagens para previsão do total populacional neste cenário
ao longo de instantes de tempo, cujos dados são obtidos por amostragem adaptativa
por conglomerados, são: o simples ajuste de forma independente ao longo do tempo do
modelo (3.1); estimação para cada tempo com base no estimador de Horvitz-Thompson
modificado (2.4).
O objetivo desta seção é comparar o modelo de crescimento com as abordagens citadas
acima. Desta forma, para cada tempo t uma amostra adaptativa é selecionada e: (i)
com base em todas as amostras observadas o modelo (3.4) é ajustado; (ii) à medida
que uma amostra é coletada, o modelo (3.1) é ajustado com base nestes dados e assim
55
0.4
0.3
Erro para T
EAR
EQMR
0.1
0.2
0.4
0.3
0.2
0.1
Erro para T
EAR
EQMR
0
10
20
30
40
50
0
10
20
t
20
50
30
40
0.7
0.6
Amplitude média (T)
50
Decrescimento
Crescimento
0.5
0.90
0.80
0.70
Cobertura de T
Decrescimento
Crescimento
10
40
(b) EQMR e EAR para T (decrescimento)
0.8
(a) EQMR e EAR para T (crescimento)
0
30
t
0
t
10
20
30
40
50
t
(c) Cobertura T (cresc. e decresc.)
(d) Amplitude T (cresc. e decresc.)
Figura 3.7: Sumário da distribuição a posteriori do total populacional a cada instante de
tempo T para 100 populações em crescimento e outras 100 em decrescimento geradas.
São apresentados os EQMR, EAR, probabilidade de cobertura e amplitude média dos
intervalos HPD de 95%.
sucessivamente para todo t = 1, . . . , L; (iii) estima-se o total com base no estimador de
Horvitz-Thompsom para cada tempo separadamente. Note que, exceto para a primeira
metodologia, não estarı́amos incluindo a estrutura dinâmica à estimação em nenhuma
das outras abordagens e, portanto, não seria possı́vel usar todos os dados coletados ao
longo de L tempos para estimar o total populacional, ao contrário da proposta em (3.4).
Utilizando as mesmas 100 populações geradas no estudo anterior foram analisados
o modelo de crescimento proposto (3.4) e as duas abordagens descritas acima. Para
os modelos considerados foram usadas as mesmas distribuições a priori descritas nas
Subseções 3.2.4 e 3.1.2.
56
Dado que a convergência foi obtida para todos os parâmetros e como nosso maior
interesse está na previsão do total populacional com base em modelos, na Figura (3.8)
estão em (a) as probabilidades de cobertura dos intervalos HPD de 95% para T e em (b)
as amplitudes médias relativizadas destes intervalos para todos os 50 tempos, usando o
modelo proposto em (3.4) e as replicações ao longo do tempo do modelo (3.1). Chamamos
de “Crescimento”o modelo proposto em (3.4), “Estático”o modelo estático em (3.1) e “HT”o estimador de Horvitz-Thompson.
Note que apesar do ajuste independente fornecer uma probabilidade de cobertura
em média mais próxima do nı́vel desejado de 95% que a extensão proposta, temos
uma incerteza significativamente maior. Por outro lado, em (c) temos um gráfico de
dispersão com a raiz quadrada do erro quadratico médio relativo (REQMR) sob as
mesmas abordagens. Como todos os pontos encontram-se abaixo da reta, conclui-se que
a extensão proposta em (3.4) produz erros menores, logo em termos de estimação pontual
parece ser mais vantajoso. Finalmente, em (d) estão os diagramas boxplot com os REQMR
para todos os tempos sob as duas abordagens baseadas em modelos e adicionalmente
para o estimador de Horvitz-Thompson modificado. É possı́vel concluir que o ajuste
baseado em modelos de superpopulação produz erros menores que a abordagem usando
a aleatorização do plano amostral e que o modelo proposto é o mais eficiente neste caso,
já que usamos as observações coletadas até o tempo L na estimação dos parâmetros.
Portanto, conclui-se que a extensão em (3.4) proposta para populações dinâmicas é
vantajosa em relação ao estimador de Horvitz-Thompson e ao ajuste repetido ao longo
do tempo do modelo (3.1).
57
5
●
4
●
●
3
●
●
1
2
Amplitude média
6
0.95
0.90
0.85
0.80
Probabilidade de cobertura
●
Estático
Crescimento
●
●
●
Estático
(b)
●
●
●
●
2
3
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
0
●
●
●● ●●
●
●
●
●
●
●●
●●
●
●
●
●
●
●
●
●
●
1
1.0
REQMR para T
4
1.5
●
0.5
REQMR para T (Crescimento)
(a)
Crescimento
0.5
1.0
1.5
H_T
Estático Crescimento
REQMR para T (Estático)
(c)
(d)
Figura 3.8: Comparação do modelo proposto de crescimento exponencial (3.4) com o
ajuste independente ao longo do tempo do modelo (3.1). Em (a) estão as probabilidades
de cobertura dos intervalos HPD de 95%, em (b) a amplitude média destes intervalos,
em (c) está a REQMR para cada abordagem utilizada e em (d) as REQMR para todos
os tempos incluindo na comparação o estimador de Horvitz-Thompson.
3.3
Conclusões
Neste capı́tulo foi apresentada a proposta de inferência baseada em modelos
introduzida por Rapley e Welsh (2008) para populações raras e agrupadas cujas amostras
são selecionadas de forma adaptativa por conglomerados. O modelo foi avaliado sob
58
diversos estudos simulados, variando o grau de esparsidade e agrupamento da população,
o tamanho da amostra selecionada e ainda quando na inferência o plano amostral
não é incorporado erroneamente na expressão da verossimilhança. O modelo teve um
bom desempenho em todos os casos, principalmente para valores de α e β que gerem
populações raras com um número mı́nimo de redes não vazias.
O modelo de Rapley e Welsh (2008) é construı́do em um nı́vel agregado da população,
o que produz algumas facilidades na inferência, no entanto, esta estrutura agregada induz
a algumas restrições nas hipóteses do modelo, como a homogeneidade entre grupos e a
hipótese de que o número de observações numa rede está relacionado diretamente com
seu tamanho.
É comum que as populações consideradas neste trabalho, além de terem um
comportamento raro e agrupado, apresentem uma constante mobilidade dentro de uma
região ao longo de um perı́odo de tempo. Visando a este problema, foi introduzida uma
extensão dinâmica do modelo agregado. O modelo se ajusta a amostras adaptativas
coletadas de forma independente ao longo do tempo e supõe uma dinâmica populacional
de crescimento e decrescimento ao longo do tempo. Na inferência sob o modelo de
crescimento é preciso observar as amostras ao longo de todos os tempos, ou de pelo
menos uma quantidade razoável destes, o que exige maior custo operacional. Por outro
lado, o modelo em (3.1) só permite o uso de amostras independentes ao longo do tempo
e previsão para cada tempo separadamente, isto porque o modelo escrito desta forma
agregada apenas nos fornece informações numéricas sobre as redes, portanto não seria
possı́vel incorporar a este modelo um planejamento amostral diferente, caso contrário
necessitarı́amos de informações adicionais, como a localização das unidades que foram
selecionadas. Por esta e outras razões, um modelo desagregado que informe as localizações
pode ser de interesse. No próximo capı́tulo será proposto um modelo que atenda a estas
necessidades.
59
Capı́tulo 4
Modelo de mistura para populações
raras e agrupadas sob amostragem
adaptativa
O modelo proposto por Rapley e Welsh (2008) usa as redes como unidades de análise,
de forma a não ter que introduzir componentes espaciais no modelo, o que pode vir a
facilitar a inferência.
Neste caso, a modelagem é feita supondo que em média a distribuição do total
das redes seja proporcional ao tamanho das redes. Isto equivale a tratar as unidades
populacionais como sendo homogêneas e que a intensidade do fenômeno em uma rede
depende do seu tamanho, ou seja, redes maiores apresentam em média sempre maior total.
No entanto, esta suposição nem sempre é válida. Por exemplo, é comum que a intensidade
de um fenômeno em uma unidade varie de acordo com a rede à qual ela pertence devido
à influência da vizinhança. Ou ainda, dentro de uma mesma rede é possı́vel que as
unidades de borda apresentem menor taxa de ocorrência do que as unidades no centro da
rede. Além disso, uma rede pode ter maior incidência do fenômeno em suas unidades não
somente por ser maior, mas por outros fatores externos que influenciem na sua disposição.
Mas o fato de agregar a informação para todas as unidades dentro de uma mesma
rede, além de não permitir previsão do total populacional em cada unidade da grade
regular, impossibilita a incorporação de estruturas mais complexas úteis na inserção das
60
suposições descritas acima. Ou ainda, num contexto de populações móveis, como visto
na Seção 3.2, o modelo não possibilita a inserção de um planejamento amostral que ao
longo do tempo use informações de tempos anteriores, pois o modelo não apresenta uma
variável de identificação das unidades i da rede j que pertencem à amostra, e sim uma
variável que agrega numericamente estas informações para cada rede. Por estas razões,
a proposta de um modelo desagregado pode ser interessante em muitos contextos com
populações raras e agrupadas.
Por outro lado, a modelagem de eventos raros usando distribuição de Poisson algumas
vezes revela uma significante sobredispersão, a qual pode ser diminuı́da usando modelos
mistos hierárquicos. Esta abordagem vem sendo usada em muitas aplicações, como por
exemplo na modelagem de doenças raras, em que o número de casos por área é pequeno,
como pode ser visto em Clayton e Bernardinelli (1992). Viallefont et al. (2002) sugerem
para a modelagem de eventos deste tipo um modelo de mistura Poisson, cujo número de
componentes de mistura é desconhecido.
O objetivo deste capı́tulo é propor um modelo de mistura desagregado que suponha
heterogeneidade entre unidades pertencentes a redes distintas e, portanto, o total em
cada rede não dependeria somente do tamanho desta. A proposta é que este modelo se
aplique a populações raras e agrupadas, que são amostradas usando o desenho adaptativo
por conglomerados, portanto a probabilidade de seleção deve ser incorporada à função
de verossimilhança dos parâmetros do modelo. Note que o fato de modelar cada unidade
da grade permite também construir um modelo com suposição de heterogeneidade entre
as células de uma mesma rede, o qual será nosso interesse futuro.
Na Seção 4.1 é definida a classe de modelos de mistura de distribuições de
probabilidades e uma forma de fazer inferência sob o enfoque Bayesiano para modelos
deste tipo. Na Seção 4.2 é apresentado o modelo proposto neste trabalho, discutindo
pontos como elicitação de distribuição a priori, inferência e convergência das cadeias
com as amostras da distribuição a posteriori do vetor paramétrico. Vários estudos
simulados são realizados com o objetivo de avaliar o desempenho do modelo sob diferentes
configurações da população. Finalmente, na Seção 4.4 o modelo proposto neste capı́tulo
é comparado ao modelo de Rapley e Welsh (2008), a partir de experimentos baseado em
61
modelos e no desenho amostral usando a população real de marrecos da asa azul, descrita
na Seção 3.1.3 do Capı́tulo 3. Daqui em diante, o modelo de Rapley e Welsh (2008) será
referido como “modelo agregado”.
4.1
Uma revisão sobre modelos de mistura de
distribuições
Modelos com mistura de distribuições são frequentemente utilizados para modelar
fenômenos cujas observações são provenientes de uma população composta por k
subpopulações, onde k pode ser conhecido ou desconhecido. Um modelo com mistura é
dado por uma soma ponderada de distribuições de probabilidades. Vejamos a definição
a seguir.
Definição 4.1.1 Qualquer combinação linear convexa
k
X
wj f (Y | φj ), com 0 < wj < 1 e
j=1
k
X
wj = 1,
(4.1)
j=1
das distribuições f (· | φj ) pertencentes a uma famı́lia de distribuições paramétricas
indexadas pelo vetor paramétrico φj , é denominada uma mistura de distribuições com
k componentes, tal que wj , j = 1, . . . , k são os pesos da mistura.
O modelo em (4.1) assume que temos uma população heterogênea, com k
subpopulações, de tamanhos proporcionais aos pesos wj , j = 1, . . . , k.
Em Marin et al. (2005) pode-se ver uma revisão desta modelagem e exemplos de
dados aos quais se aplica esta abordagem, como os dados de peso de militares recrutados
na França, que apresentam um comportamento bimodal de acordo com seu lugar de
origem. Assim, cada peso yi é proveniente a priori da densidade f1 ou f2 , em que f1 está
modelando os pesos dos homens das planı́cies e f2 os pesos dos homens das montanhas,
com probabilidades w1 e w2 = 1 − w1 . Note que a diferença fundamental entre uma
regressão simples usando o lugar de origem como covariável, para um modelo de mistura
deste tipo, é que as observações são coletadas indiscriminadamente para toda a população,
ou seja, não se conhece o lugar de origem de todos os militares, apenas supõe-se a priori
62
que a variável resposta (peso) é influenciada por fatores (origem) que podem ou não
ser conhecidos. Logo, a estrutura de mistura é cabı́vel neste caso devido à perda de
informação sobre a origem de cada homem.
A função de verossimilhança de um modelo de mistura com k componentes, para uma
amostra y = (y1 , . . . , yn )0 é dada por:
[y | w, φ] =
n X
k
Y
wj f (yi | φj ),
(4.2)
i=1 j=1
onde φ = (φ1 , . . . , φk )0 e w = (w1 , . . . , wk )0 . Esta função tem uma forma complexa
pois envolve uma expansão em k n termos, o que torna computacionalmente custoso o
desenvolvimento de estimadores de máxima verossimilhança, ou no contexto Bayesiano,
a obtenção da distribuição a posteriori do vetor paramétrico.
Como uma alternativa a este problema, a estrutura oculta é explorada para facilitar o
procedimento de estimação dos parâmetros. Utilizando o fato de que, para todo o vetor
aleatório Y, proveniente de um modelo com mistura de distribuições com k componentes,
é possı́vel associar uma variável latente Z, de dimensão n, que indica a componente da
qual a observação Yi é proveniente, isto é Zi = j, se a unidade i é proveniente da
componente j, i = 1, . . . , n, j = 1, . . . , k.
O vetor de dados observados em conjunto com as variáveis latentes produz mais
informação para o modelo e passa a ser chamado de dados aumentados. Segundo Tanner
(1993), o objetivo de aumentar os dados é simplificar a forma analı́tica da função de
verossimilhança, condicionando-os à variável latente. Condicional a Z = (Z1 , . . . , Zn )0 ,
a função de verossimilhança em (4.2) passa a ser escrita em termos de produtórios de
densidades simples, da forma:
[y | Z, w, φ] =
k
Y
Y
f (yi | φj ),
j=1 {i:Zi =j}
tal que Zi ’s são supostamente independentes, com função de probabilidade dada por
P (Zi = j) = wj ,
para j = 1, . . . , k e i = 1, . . . , n. Ao integrar sobre as variáveis latentes Z, retorna-se à
expressão em (4.1).
63
Finalmente, estimar o número de componentes de mistura é uma questão importante
e complexa. Ao atribuir um número menor que o necessário, o modelo não consegue
capturar a verdadeira estrutura dos dados. Por outro lado, se esse número for superior
ao ideal, o modelo torna-se menos parcimonioso e atribui massa de probabilidade
desnecessária em algumas regiões do espaço e, consequentemente, a densidade fica
subestimada e não identificável.
4.1.1
Inferência Bayesiana em modelos de mistura
Abordagens Bayesianas para inferência em modelos de mistura têm despertado grande
interesse entre pesquisadores. Além de permitir a inclusão de conhecimento a priori sobre
os parâmetros do modelo na análise, diminui a complexidade do modelo decompondo-o
em estruturas mais simples. Segundo Richardson e Green (1997) o paradigma Bayesiano
é o mais adequado ao contexto de misturas, principalmente quando o número de
componentes é desconhecido e deve ser estimado.
No contexto Bayesiano o modelo de mistura deve ser completado com uma distribuição
a priori para o vetor paramétrico Θ = (k, w, φ). Como pode ser visto em Richardson
e Green (1997), supondo independência a priori e que [φ | Z, w, k] = [φ | k] e [y |
φ, Z, w, k] = [y | φ, Z] a distribuição conjunta das variáveis do modelo é dada por:
[k, w, Z, φ, y] = [k][w | k][Z | w, k][φ | k][y | φ, Z].
Para completa flexibilidade, atribuem-se aos hiperparâmetros da distribuição a priori
também distribuições a priori independentes.
4.1.1.1
Identificabilidade do modelo
Uma caracterı́stica importante de um modelo de mistura é que este apresenta-se
invariante sob permutações dos ı́ndices dos seus componentes.
Isto implica que os
parâmetros φi não são marginalmente identificáveis, ou seja a partir da função de
verossimilhança não é possı́vel distinguir por exemplo φi de φj , para i 6= j. Nesta classe
de modelos, esta identificabilidade é tratada na realização da inferência Bayesiana.
64
Primeiramente, note que se (φ1 , . . . , φk )0 é máximo local, logo qualquer permutação
dentro deste vetor também o é, portanto existem “k!” modas. Além disso, se é utilizada
uma distribuição a priori para φ permutável, todas as condicionais completas são
idênticas para todas as componentes, então a distribuição a posteriori para φ se apresenta
multimodal o que dificulta a análise e interpretação dos resultados.
Considere por
exemplo, uma população constituı́da de duas distribuições Normais, inequivocamente
rotuladas. A distribuição a posteriori das duas médias irão se sobrepor, mas a extensão
desta sobreposição depende da separação entre elas e do tamanho da amostra. Quando as
médias são bem diferentes, o rótulo na distribuição a posteriori ao ordenar suas médias
geralmente coincidem com a rotulagem da população. Mas, se a diferença diminui, o
fenômeno conhecido como label switching tende a ocorrer.
Portanto, uma alternativa usada neste problema é a imposição de um único tipo
de rótulo para as componentes. Por exemplo, num caso de mistura de distribuições
normais, em que φj = (µj , σj2 ), tal que µj e σj2 são respectivamente a média e a variância
da distribuição para a componente j, pode-se identificar as componentes de acordo
com a ordem crescente da média. Desta forma, a distribuição a priori conjunta é “k!
multiplicado pela distribuição original com a restrição para identificabilidade imposta”.
Vale ressaltar que em alguns casos a melhor alternativa será impor tal restrição na média,
outras vezes na variância, outras ainda no peso, como pode ser visto em Richardson e
Green (1997).
4.1.1.2
Algoritmo MCMC para modelos de mistura supondo k desconhecido
Em geral, os métodos de MCMC são utilizados para amostrar da distribuição a
posteriori do vetor paramétrico. Os métodos de MCMC inicialmente eram utilizados
apenas para problemas em que a distribuição a posteriori tivesse uma densidade com
respeito a uma medida subjacente fixa e, portanto, não podiam ser utilizados em casos,
como na mistura de distribuições, em que o tamanho do espaço paramétrico é também
um parâmetro desconhecido. Alguns trabalhos na literatura surgiram propondo métodos
de MCMC para problemas de dimensão variante, entre elas a abordagem denominada
MCMC com saltos reversı́veis (do inglês, Reversible Jump Markov Chain Monte Carlo,
65
RJMCMC) ganhou destaque. Proposto em Green (1995), o algoritmo RJMCMC é como
um algoritmo de Metropolis-Hastings que permite a movimentação entre modelos que
possuem espaços paramétricos de diferentes dimensões.
Em particular, Richardson e Green (1997) desenvolveram uma metodologia Bayesiana
para estimação em modelos de mistura com número de componentes desconhecido usando
métodos de RJMCMC. O método é brevemente descrito a seguir.
Considere que temos os modelos M1 , . . . ,Mm , em que o modelo Mj , j = 1, . . . , m
é indexado pelo vetor paramétrico Θj pertencente ao espaço paramétrico Φj . Suponha
que a distribuição a priori para (Θj ,Mj ) é dada pelo produto entre [Θj |Mj ] e [Mj ].
Para este estado corrente, propõe-se um movimento do modelo Mj para o modelo Ml ,
tal que l = 1, . . . , m, j 6= l com probabilidade pl|j . Como os espaços paramétricos Θj e Θl
possuem dimensões diferentes, é preciso completar um dos espaços com espaços artificiais
adequados, para criar assim uma bijeção entre eles. Isto é feito, em geral, aumentando
o espaço paramétrico do modelo com menor dimensão. Ou seja, se dim(Θj ) < dim(Θl ),
gera-se um vetor aleatório u, independente de Θj , de uma distribuição q(u) de dimensão
igual à diferença dim(Θl ) − dim(Θj ) e obtém-se Θl usando uma transformação ϕ : j → l
definida por T (Θj , u) = Θl . Este movimento proposto é aceito com probabilidade igual
a min(1, A), tal que
[y | Θl , Ml ] [Θl | Ml ] [Ml ] pj|l ∂Θl A=
,
[y | Θj , Mj ] [Θj | Mj ] [Mj ] pl|j q(u) ∂(Θj , u) (4.3)
onde o último termo é o Jacobiano da transformação ϕ : j → l.
O movimento contrário, ou seja de Θl para Θj , supondo dim(Θj ) < dim(Θl ), é feito
usando a transformação inversa, logo o valor proposto para u é determinı́stico. Assim, o
movimento inverso é aceito com probabilidade min(1, A−1 ).
Em geral, as probabilidades pj|l são construı́das de forma que qualquer um dos
movimentos tenham a mesma probabilidade de serem propostos, a não ser que no passo
corrente o valor de k seja 1 ou kmax , em que kmax é o maior valor que k pode assumir a
priori. Se k = 1 não seria possı́vel excluir alguma componente e diminuir a dimensão do
66
espaço paramétrico, e se k = kmax não seria possı́vel incluir alguma. Nesses casos, então
só é possı́vel propor um dos dois movimentos.
Richardson e Green (1997) apresentam o método de inferência RJMCMC para um
modelo de mistura normal. Restritos a uma vizinhança de modelos Ml e Mp , tais que
dim(Θl ) = k + 1 e dim(Θp ) = k − 1, Richardson e Green (1997) utilizam o algoritmo
RJMCMC com os movimentos de “divisão”/ “combinação” e “inclusão”/ “exclusão”
para estimar o valor de k. A parte de “inclusão”/ “exclusão” é inicialmente usada no
método, com o objetivo de tratar de amostras finitas que não apresentem observações de
todos os grupos da população.
A proposta de “inclusão” consiste em adicionar uma nova componente j ∗ vazia na
mistura, com novos parâmetros (wj ∗ , φj ∗ ) gerados de distribuições propostas. Além disso,
os pesos devem ser reescalados para somar 1, fazendo wj0 = wj (1 − wj ∗ ), para todo j 0 .
A proposta de “exclusão” consiste em remover uma componente vazia da mistura e
reescalar os pesos para somar 1, fazendo wj0 =
wj
.
1−wj ∗
Estes passos são aceitos ou não
com probabilidade dadas respectivamente por min(1, A) e min(1, A−1 ), para A como
descrito em (4.3).
Para evitar problemas de alta taxa de rejeição da proposta “inclusão” sob distribuição
a priori não informativas, Richardson e Green (1997) propõem ainda o uso de movimentos
do tipo “divisão”/ “combinação” de componentes existentes.
No movimento de
“divisão”propõe-se a passagem de um modelo Mj de dimensão k para um modelo Ml
de dimensão k + 1 da seguinte forma: escolhe-se aleatoriamente uma componente de
mistura indexada por j ∗ e propõe-se a divisão desta única em um par (j1 , j2 ). Neste caso,
é necessário além de realocar as observações Yi tais que zi = j ∗ em zi = j1 e zi = j2 ,
definir os novos valores de (wj1 , wj2 , φj1 , φj2 )0 . Estas componentes devem ser adjacentes
em relação aos valores de φ, por conta da identificabilidade do modelo, como já discutido
anteriormente. Por outro lado, ao propor o movimento de Mj para Mp , precisa-se
diminuir a dimensão do espaço paramétrico correspondente em 1 unidade. De forma
análoga, isto é feito a partir da escolha aleatória de um par de componentes de mistura
indexados por (j1 , j2 ) e propondo a combinação deste par em uma única componente j ∗ .
Este movimento é chamado de “combinação”.
67
A probabilidade de aceitação das propostas “divisão” e “combinação” também são
dadas, respectivamente por min(1, A) e min(1, A−1 ), para A dado em (4.3).
A inferência prossegue usando MCMC para amostrar da distribuição a posteriori dos
parâmetros. Em Richardson e Green (1997) podem ser vistos maiores detalhes sobre esta
metodologia, incluindo a elicitação da distribuição a priori para o vetor paramétrico, a
eficiência dos passos do RJMCMC para modelos de mistura e uma aplicação a dados
reais.
Dentre os aspectos observados no uso do RJMCMC para estimação em modelos de
mistura normais, em Richardson e Green (1997) destaca-se o fenômeno conhecido como
label switching. Este comportamento pode ocorrer mesmo que a restrição no rótulo dos
parâmetros seja imposta na distribuição a priori e caracteriza-se pela invariância da
função de verossimilhança sob nova rotulagem das componentes da mistura, conduzindo
a uma distribuição a posteriori dos parâmetros sendo altamente simétrica e multimodal,
dificultando assim sua sumarização.
O modelo proposto neste trabalho é um modelo de mistura de distribuições de Poisson.
Na literatura, o trabalho de Viallefont et al. (2002) é um dos que usam esta classe de
modelos. No entanto, vale ressaltar que em todos os modelos citados supõe-se que o plano
amostral é não informativo e, portanto, a forma de seleção da amostra não contribui
para a função de verossimilhança do modelo. Com esse objetivo, na próxima seção
será proposto um modelo de superpopulação de mistura, aplicável a populações raras e
agrupadas, para dados coletados por amostragem adaptativa por conglomerados.
4.2
Modelo de mistura Poisson proposto
O modelo proposto a seguir pode ser aplicado a populações raras e agrupadas,
dispostas sobre uma grade regular, sob o enfoque Bayesiano. Tal modelo apresentase como uma alternativa mais flexı́vel ao modelo agregado, pois usa como unidade de
análise as próprias unidades da grade e supõe heterogeneidade entre as redes, e assim
também não necessariamente redes maiores devem ter mais observações, ou vice-versa.
Além disso, tratar a modelagem no nı́vel da unidade primária, permitiria incorporar
68
estruturas na média que supõem heterogeneidade para unidades dentro de uma mesma
rede. Por enquanto, o interesse está apenas no primeiro caso.
Diferentemente dos modelos de mistura anteriormente apresentados, como o de
Viallefont et al. (2002), o modelo proposto é ajustado a dados obtidos sob um plano
amostral informativo, logo deve-se incluir a probabilidade de seleção da amostra na função
de verossimilhança e, além disso, as interpretações das variáveis mudam com relação aos
problemas comuns de mistura. Em geral, o objetivo, ao ajustar um modelo de mistura,
é fazer inferência acerca de: k, φj , wj , para j = 1, . . . , k, no entanto, neste caso, como
se trata de um modelo de superpopulação, o vetor paramétrico do modelo é composto
também por partes das variáveis que não foram observadas e o principal objetivo é fazer
previsão. O modelo proposto é descrito a seguir.
Considere uma população rara com N unidades, das quais X apresentam uma
caracterı́stica de interesse, ou seja são unidades não-vazias e estão divididas em R redes
não-vazias. Logo, tem-se N − X redes vazias. Seja Yi a contagem deste determinado
fenômeno de interesse na unidade não-vazia i, i = 1, . . . , X, logo Yi ≥ 1. Como se tratam
de populações raras, ou seja, cujo número de unidades vazias é extremamente alto, assim
como Rapley e Welsh (2008), vamos modelar apenas as unidades não-vazias da grade.
Suponha que a rede j não-vazia é composta por Cj unidades primárias, j = 1, . . . , R. Para
facilitar o procedimento de inferência, é natural definir uma variável aleatória latente de
alocação i , supostamente independentes para todo i e tais que P (i = j) = wj = Cj /X,
j = 1, . . . , R. Dado o valor da variável i , as contagens Yi nas redes não-vazias seguem
uma distribuição de Poisson truncada em 0 independente, cuja média se altera de acordo
com a rede à qual pertence. O modelo, para j = 1, . . . , R, pode ser escrito da seguinte
forma:
69
Yi | i = j, λj , X ∼
Poisson Truncada independente(λj ), Yi ≥ 1,
(4.4a)
P (i = j) = wj = Cj /X,
(4.4b)
λ | θ ∼ [. | θ, R],
(4.4c)
R
X
1
Ci = X,
C | X, R ∼ 1R + Multinomial (X − R, 1R ),
R
i=1
(4.4d)
R | X, β ∼
Binomial Truncada (X, β), R = 1, . . . , X,
(4.4e)
X|α ∼
Binomial Truncada (N, α), X = 1, . . . , N,
(4.4f)
em que [. | θ, R] representa a distribuição a priori de λ = (λ1 , . . . , λR )0 , a qual
depende do número de grupos R e de um vetor de hiperparâmetros θ. Lembrando que
esta distribuição deve satisfazer alguma restrição no ı́ndice, devido à identificabilidade
do modelo. Além disso, os parâmetros α, β e θ podem ser desconhecidos e, portanto,
são atribuı́das distribuições a priori independentes a estes também. Denote por [.] cada
uma destas distribuições.
Este modelo, somente aplicável às X unidades não-vazias de uma população,
apresenta a estrutura de uma mistura de probabilidades, cujas componentes são as R
redes não-vazias, supostamente heterogêneas, e seus pesos são proporcionais ao número de
unidades nas redes, Cj , j = 1, . . . , R. No entanto, o fato do modelo proposto ser ajustado
a dados provenientes de amostragem adaptativa por conglomerados cria uma grande
diferença nas variáveis entre um modelo de mistura comum e este modelo proposto. Num
modelo de mistura, a amostra selecionada somente traz informações acerca da variável
Y ; por exemplo, não é possı́vel saber a qual grupo a unidade i observada pertence, pois
isto é uma divisão artificial. No entanto, neste modelo, a amostragem por conglomerados
adaptativos traz informações acerca de todas as variáveis. Ao coletar-se uma amostra
adaptativa, além de observar a variável Y para as unidades amostradas, sabe-se as redes
às quais elas pertencem e o tamanho desta rede, pois de acordo com este plano amostral
se uma unidade não-vazia é selecionada aleatoriamente, toda a rede é observada.
Dessa forma, é adequado dividir as variáveis indicando pelo ı́ndice s a parte observada
e por s̄ a parte não observada. Sejam então X = Xs + Xs̄ , R = Rs + Rs̄ , = (0s , 0s̄ )0 ,
70
C = (C0s , C0s̄ )0 , Y = (Ys0 , Ys̄0 )0 . Neste caso, o objetivo está em não somente estimar os
parâmetros do modelo com base numa amostra, mas também fazer previsão das partes
X
X
não-observadas. O objetivo final é prever o total populacional T =
Yi .
i=1
Finalmente, como o modelo aplica-se a dados coletados de forma adaptativa e este
planejamento amostral é não-ignorável, a probabilidade de seleção deve ser acrescentada
à função de verossimilhança completa, de forma a trazer mais informações ao processo
de estimação dos parâmetros. Como no modelo (3.1), a probabilidade de seleção de uma
particular amostra s = {i1 , . . . , im }, composta por m redes, é dada por:
[s | X, R, C] =
m
Y
l=1
Zi × gil ,l
,
PN −X+R l
P
Zi − j−1
i=1
k=0 Zik
(4.5)
onde gil ,l é o número de redes de tamanho Zil que restam após a seleção de j − 1 redes e
Zi0 = 0. O vetor Z é construı́do de forma a ter os tamanhos de todas as redes vazias e
não-vazias, ou seja Z = (C0 , 10X−R )0 .
A função de verossimilhança completa é dada por:
[s, X, R, , C, Y | λ, α, β] = [s | X, R, C][Y | , λ, X][ | C, R, X]
× [C | R, X][R | X, β][X | α]
=
m
Y
l=1
RY
s +Rs̄
Y
λYj i exp(−λj )
zil × gil ,l
×
Pj−1
PN −X+R
Y ![1 − exp(−λj )]
Zik
zi − k=0
i=1
j=1 {i:i =j} i
Cj −1
RY
RY
s +Rs̄
s +Rs̄
1
1
1
Cj
×
C × (Xs + Xs̄ − Rs − Rs̄ )!
(Xs + Xs̄ )Xs +Xs̄ j=1 j
(Cj − 1)! Rs + Rs̄
j=1




Xs +Xs̄
Rs +Rs̄
Xs + Xs̄
N
(1 − β)Xs +Xs̄ −Rs −Rs̄ 
(1 − α)N −Xs −Xs̄
β
α
.
×
×
1 − (1 − β)Xs +Xs̄
1 − (1 − α)N
Rs + Rs̄
Xs + Xs̄
A função de verossimilhança marginal é obtida da seguinte forma:
[Xs , Rs , s , Cs , Ys | λ, α, β] =
X
[s, X, R, , C, Y | λ, α, β]
{Ys̄ ,Cs̄ ,s̄ ,Rs̄ ,Xs̄ }
=
X
{Ys̄ ,Cs̄ ,s̄ ,Rs̄ ,Xs̄ }
71
[s, X, R, C | α, β][Y | , λ, X].
4.2.1
Distribuição a priori para λ
Segundo Richardson e Green (1997) usar distribuição a priori completamente não
informativa e gerar distribuição a posteriori própria não é possı́vel em modelos de mistura.
Como existem componentes da mistura que não apresentam observações na amostra,
distribuições a priori independentes impróprias e não informativas não podem ser usadas.
A alternativa neste caso é manter-se com a estrutura de independência a priori usando
distribuições pouco informativas, as quais podem ou não depender dos dados observados,
o que pode ser feito, por exemplo, inserindo estruturas a priori para os hiperparâmetros.
Por outro lado, existe uma relação direta entre a distribuição a priori de λ e a
distribuição a posteriori de R. Uma sugestão neste caso é considerar distribuições a
priori dependentes para λ, de forma a modelar a distância entre λj s consecutivos. Esta
distribuição a priori foi introduzida por Roeder e Wasserman (1997) para misturas de
normais e é muito utilizada quando deseja-se ser não informativo.
Note que, neste caso, o vetor λ pode ser definido como λ = (λ0s , λ0s̄ )0 , em que λs
refere-se à parte associada às redes observadas na amostra, para o qual espera-se obter
melhores resultados, e λs̄ refere-se à parte associada às variáveis não observadas. A fim
de garantir a identificabilidade, será imposta sobre a distribuição a priori de λ alguma
restrição sobre o ı́ndice dos parâmetros. Mas esta restrição é necessária apenas aos
elementos de λ que estão associadas às redes não amostradas, ou seja à λs̄ .
Com base nestas ideias serão utilizados dois tipos de distribuições a priori para o
parâmetro λ em (4.4c), as quais estão descritas a seguir.
4.2.1.1
Distribuição a priori independente
Primeiramente, será considerada a independência entre os λj ’s, tal que a distribuição
conjunta de λ é dada por:
[λ | θ, R] = Rs̄ ![λ1 | θ] . . . [λR | θ], tal que λj < λj+1 , para todo j ∈ [Rs + 1, Rs + Rs̄ ).
Em particular, será considerado, que
λj ∼ Gama(d, ν), j = 1, . . . , R, para θ = (d, ν)
72
e introduz-se um nı́vel hierárquico adicional assumindo que ν ∼ Gama(e, f ).
Gelman (2006) apresenta formas de elicitar a priori esta distribuição Gama. Uma
forma usual de ser não informativo é escolher valores pequenos para seus dois parâmetros,
como 0.01.
No entanto, deve-se evitar distribuições que tenham altas massas de
probabilidade no zero, o que pode incluir componentes com médias pequenas, tornando
difı́cil estimar o modelo de mistura.
Viallefont et al. (2002) relatam uma sensibilidade da distribuição a posteriori de
outros parâmetros do modelo de mistura Poisson de acordo com a escolha dos parâmetros
da Gama. Foi usada então uma distribuição a priori pouco informativa descrita em
Viallefont et al. (2002). Para d escolhe-se um valor maior que 1, por exemplo 1.1, pois
isso permite evitar a forma exponencial da distribuição sem reduzir muito o coeficiente
de variação (CV). Para o parâmetro ν escolhe-se e e f a priori tal que a aproximação
a média de λj , d/(e/f ) seja igual ao ponto médio das observações com variância e/f 2
controlada.
4.2.1.2
Distribuição a priori dependente
Esta distribuição leva em conta a informação da distância entre dois parâmetros da
Poisson para duas componentes que são consecutivas em termos dos valores de λj . Neste
caso o hiperparâmetro θ em (4.4c) é igual a uma constante positiva τ e a distribuição a
priori conjunta para λ é dada por:
[λ | τ 2 , R] = [λR | λR−1 , τ 2 ][λR−1 | λR−2 , τ 2 ] . . . [λ1 ],
onde [λj | λj−1 , τ 2 ] é N(λj−1 ,∞) (λj−1 , τ 2 ), que denota a densidade de uma Normal centrada
em λj−1 com variância τ 2 , truncada para ser maior que λj−1 e [λ1 ] ∝ 1, o que garante a
identificabilidade do modelo. Esta distribuição indica baixa probabilidade a priori que
duas redes vizinhas sejam mais distantes que τ desvio padrões.
Segundo Viallefont et al. (2002) uma vantagem deste modelo é que o hiperparâmetro
τ 2 , o qual controla a distância entre as médias de duas componentes e sua variabilidade,
é explı́cito, e controla o número de grupos. Eles discutem as dificuldades de elicitar τ 2 e
a influência deste hiperparâmetro na distribuição a posteriori do vetor paramétrico, em
73
especial na de R. Por exemplo, se τ 2 é pequeno quando comparado à verdadeira distância
entre dois λj s consecutivos, há uma tendência em ajustar componentes intermediários
entre os verdadeiros e assim obter uma distribuição a posteriori favorecendo valores
mais altos para R. Baseado num estudo de simulação, Roeder e Wasserman (1997)
recomendam assumir τ = 5, pois esta escolha resulta em resultados razoáveis.
O modelo proposto (4.4) é um modelo de superpopulação e seu ajuste depende
da estimação de parâmetros e previsão de quantidades populacionais que não foram
observadas na amostra. Além disso, tal modelo aplica-se a cenários com populações
que podem ser extremamente raras e agrupadas. Logo, para um tamanho de amostra
relativamente pequeno, podem ser selecionadas amostras com poucas unidades nãovazias e pouco representativas, o que deve produzir estimativas inadequadas para o total
populacional. Para estes casos, recomenda-se elicitar distribuições a priori informativas.
4.2.2
Inferência para o modelo
Como o modelo é descrito por um vetor paramétrico Θ = (Xs̄ , Rs̄ , s̄ , Cs̄ , Ys̄ , α, β, λ)
de dimensão desconhecida, o algoritmo RJMCMC será também utilizado neste caso,
como apresentado na Seção 4.1.1.2 para modelos normais. Basicamente o procedimento
de estimação consiste dos seguintes passos:
(1) atualização de α, β, θ e λ;
(2) atualização das variáveis não observadas Xs̄ e Ys̄ ;
(3) atualização da alocação s̄ e diretamente Cs̄ é atualizado;
(4) proposta de “divisão” de uma rede em duas ou “combinação” de duas redes em
uma.
As distribuições condicionais completas podem ser vistas no Apêndice B. Será descrito
a seguir com detalhes o passo (4).
De forma análoga ao procedimento de inferência descrito em Viallefont et al. (2002),
serão utilizados os momentos de ordem zero e primeira ordem na proposta de “divisão”
de uma componente da mistura j ∗ em duas novas j1 e j2 , mas como a distribuição da
74
variável Yi , i = 1, . . . , X é Poisson Truncada, diferentemente de Viallefont et al. (2002),
o momento de primeira ordem não é λj , j = 1, . . . , R, e sim λ0j = λj /{1 − exp(−λj )}.
Logo, os parâmetros propostos satisfazem as seguintes equações:
wj ∗ = wj1 + wj2 ,
wj ∗ λ0j ∗ = wj1 λ0j1 + wj2 λ0j2 ,
tal que λ0j−1 < λ0j1 < λ0j2 < λ0j+1 , devido a questões de identificabilidade do modelo.
Mas, para valores de λj razoavelmente grandes, λj e λ0j se aproximam, como podemos
ver na Figura 4.1, logo nestes casos as equações acima podem ser escritas em função
dos λj ’s.
Por outro lado, para os casos em que esta aproximação não é válida, a
solução é estimar λ0j , e quando necessário expressar λj em função de λ0j , como na
função de verossimilhança, uma aproximação numérica faz-se útil como, por exemplo,
a aproximação de Taylor. Isto porque esta função, apesar de ser inversı́vel, envolve um
polinômio com uma exponencial, para a qual, em geral, é impossı́vel obter uma solução
analı́tica exata. Neste trabalho, utilizamos a aproximação pela própria função identidade
7
em todos os exemplos, pois são estudados casos em que λj é razoavelmente grande.
λ
6
λ
0
1
2
3
4
5
1 − exp(− λ)
0
1
2
3
4
5
6
7
λ
Figura 4.1: Comparação das médias da distribuição de Poisson e Poisson truncada no
zero.
Para determinar os parâmetros associados a estas novas componentes, basta resolver
o sistema de equações anterior. Mas, como tem-se um sistema com 4 incógnitas e 2
equações, para resolvê-lo é preciso completá-lo gerando um vetor aleatório u = (u1 , u2 ).
Viallefont et al. (2002) consideram 3 formas diferentes de fazer isso, as quais baseiam-se
75
em diferentes intuições de como induzir a positividade dos parâmetros da Poisson. Neste
trabalho, será utilizada apenas uma destas, a qual baseia-se em adição de vizinhos de λj ∗
dependentes e está descrita a seguir.
São geradas duas variáveis auxiliares u1 ∼ U(0, 1) e u2 ∼ U(0, 1) e então define-se as
seguintes transformações determinı́sticas:
wj2 = wj ∗ (1 − u1 ),
w j1 = w j ∗ u 1 ,
λj1 = λj ∗ − ρu2 (1 − u1 ), λj2 = λj ∗ + ρu2 u1 ,
onde



min{(λj ∗ − λj ∗ −1 )/(1 − u1 ), (λj ∗ +1 − λj ∗ )/u1 },




 min{λ /(1 − u ), (λ − λ )/u },
1
1
2
1
1
ρ=


(λj ∗ − λj ∗ −1 )/(1 − u1 ),




 λ /(1 − u ),
1
1
1 < j ∗ < Rs̄ ,
1 = j ∗ < Rs̄ ,
1 < j ∗ = Rs̄ ,
1 = j ∗ = Rs̄ .
No Apêndice B é apresentada a expressão da probabilidade de aceitação do movimento
descrito.
4.2.2.1
Diagnóstico de convergência
Para verificar que a convergência é atingida no ajuste do modelo serão apresentados
histogramas com a distribuição a posteriori dos parâmetros e medidas que avaliam a
convergência propostas por Geweke (1992) e Raftery e Lewis (1992). A primeira medida
baseia-se em um teste de igualdade das médias da primeira e última partes da cadeia de
Markov. Se as amostras são resultantes de uma distribuição estacionária, as duas médias
devem ser iguais e a estatı́stica de teste tem assintoticamente uma distribuição Normal
padrão. A outra medida verifica a independência entre os valores gerados para a cadeia
baseado num fator de dependência, se este for maior que 5 pode-se dizer que existe forte
autocorrelação entre os valores da cadeia.
A fim de examinar o desempenho do modelo em (4.4), foram analisadas amostras da
distribuição a posteriori dos parâmetros. Para isso, gerou-se uma população artificial em
uma grade regular de tamanho N = 400 para α = 0.15 e β = 0.1 fixados. Os valores
76
de λ foram gerados aleatoriamente de uma distribuição Gama centrada em 8.5 com CV
√
fixado em 95%. Como o CV de uma distribuição Gama(d,ν) é dado por 1/ d, sob
estas condições tem-se d = 1.1 e ν = 0.13. Foi selecionada uma amostra adaptativa por
conglomerados com tamanho inicial 5%N e, em particular, a população gerada apresenta
R = 8 redes e as redes observadas na amostra são s = {2, 4, 7}, de acordo com a ordem
crescente de λ.
Assumindo que os parâmetros α, β e λ são independentes a priori considerouse a priori que α ∼ Beta(3, 15) e β ∼ Beta(1, 9).
Estas distribuições, apesar de
informativas, neste caso estão trazendo o mı́nimo de informação necessária para aplicação
deste modelo complexo, ou seja que a população é rara e agrupada. Isto porque, como
visto no Capı́tulo 3, α e β são parâmetros relacionados ao número de células não-vazias
e redes não-vazias. Para λ foram consideradas as duas distribuições a priori citadas
na Seção 4.2.1, i.e. : (i) λj | ν ∼ Gama(d, ν), para j = 1, . . . , R, independentes; (ii)
λj | λj−1 ∼ N(λj−1 ,∞) (λj−1 , τ 2 ), para j = 1, . . . , R. Para a segunda distribuição de λ
assumiu-se τ = 5, que é uma das sugestões de Roeder e Wasserman (1997).
Para a obtenção de amostras da distribuição a posteriori do vetor paramétrico
Θ = (Xs̄ , Rs̄ , s̄ , Cs̄ , Ys̄ , α, β, λ, ν) é necessário o uso de métodos de simulação estocástica,
em particular como a dimensão de Θ é também um parâmetro, utiliza-se o método
de RJMCMC, como descrito na Seção 4.2.2, com passos de Metropolis-Hastings e
Amostrador de Gibbs. Foram geradas 200.000 amostras, sendo as 10.000 primeiras
descartadas como aquecimento e amostras de 190 em 190 foram tomadas, a fim de obterse 1.000 amostras independentes resultantes.
A Tabela 4.1 apresenta o valor da estatı́stica de teste de Geweke e do fator de
dependência do critério de Raftery-Lewis. Todos os resultados mostram a convergência
das cadeias e a ausência de forte autocorrelação.
Nas Figuras 4.2 e 4.3 estão os histogramas com as densidades a posteriori para
os parâmetros α, β, ν, λ e o total populacional T , supondo distribuição a priori
para λ independente e dependente, respetivamente.
O respectivo valor verdadeiro
está representado pela linha cheia e intervalo HPD de 95% pela linha pontilhada. A
77
Tabela 4.1: Análise da convergência das cadeias a posteriori dos parâmetros do modelo
proposto supondo distribuição a priori independente e dependente para λ para uma
população artificial.
α
β
ν
T
λ1
λ2
λ3
λ4
λ5
λ6
λ7
λ8
indep.
0.7
-0.4
-1.6
0.4
1.4
-1.3
1.4
-0.4
1.5
1.5
1.2
1.5
dep.
-1.1
0.4
-
-0.8
-0.3
-1.0
1.2
1.0
0.6
-0.4
-1.0
1.2
indep.
1.3
1.1
1.1
1.8
0.9
1.0
1.0
1.0
0.9
1.0
1.1
1.1
dep.
2.5
1.1
-
3.2
0.9
1.0
1.0
1.1
1.0
0.9
1.4
1.0
Geweke
R-L
distribuição a posteriori de λj para j ∈ s̄ apresentada é condicional às amostras em que
R é estimado como o valor verdadeiro.
Note que a maioria dos parâmetros são bem estimados sob as duas distribuições a
priori, com maior densidade a posteriori em torno do valor verdadeiro e o mesmo contido
no intervalo HPD de 95%. O parâmetro populacional β apresenta um pequeno viés,
mas em todos os casos ainda contido no intervalo HPD de 95%. Alguns λj ’s para j ∈ s̄
apresentam um comportamento bimodal e baixa precisão, um comportamento que pode
ser esperado em modelos de mistura. Neste caso esta bimodalidade não influenciou na
convergência dos outros parâmetros e principalmente do total T , portanto não afetou o
desempenho do modelo.
Por outro lado, λj ’s para j ∈ s apresentam estimativas melhores, o que também era
esperado já que existem informações adicionais com respeito às redes amostradas.
78
0.3
4
5
6
7
8
200
5
10
15
20
800
0.4
0.2
Densidade
5
10
15
20
25
4
5
6
7
8
9 10
λ4
Densidade
0.00
8
λ6
600
0.0
0
Densidade
0
400
T
0.0 0.1 0.2 0.3
0.20
0.10
Densidade
40
0.4
λ3
0.00
30
λ5
0.3
0.10
Densidade
3
0.002
Densidade
0.2
0.00
2
λ2
0.12
20
0.1
ν
0.4
20
0.06
10
0.000
0.0
0.0
15
0.00
0
4
0.4
0.2
Densidade
0.20
Densidade
0.10
10
λ1
Densidade
0.2
0.10
0.1
β
0.00
5
2
Densidade
0
0.0
α
0
6
8
6
4
0
2
Densidade
12
8
4
Densidade
0
0.05 0.10 0.15 0.20
9
11
λ7
13
0
5
10
15
20
25
λ8
Figura 4.2: Densidade a posteriori para alguns parâmetros do modelo proposto e para o
total populacional T com base em um dado artificial supondo distribuição a priori para
λ independente. A linha vertical cheia representa o valor verdadeiro e a linha pontilhada
o intervalo HPD de 95%.
79
0.25
0.0
0.2
0.6
5
10
15
0.20
0.10
15
4
6
8
15
10
14
5
10
Densidade
15
20
λ5
λ4
0.00
0.00
10
λ6
0.10
Densidade
20
0.30
Densidade
0.20
0.10
0.00
Densidade
10
λ1
0.20
0
5
5
0.00
Densidade
12
λ3
0
0
0.20
10
800
0.10
8
λ2
0.15
6
600
0.0 0.1 0.2 0.3
0.20
0.00
4
400
T
0.10
Densidade
0.30
0.15
0.00
2
Densidade
200
β
α
Densidade
0.4
0.00
0.000
0
0.15
0.002
Densidade
3
2
1
Densidade
4
12
Densidade
0 2 4 6 8
0.05
8
10
12
λ7
14
16
0
5
10
15
20
λ8
Figura 4.3: Densidade a posteriori para alguns parâmetros do modelo proposto e para o
total populacional T com base em um dado artificial supondo distribuição a priori para
λ dependente. A linha vertical cheia representa o valor verdadeiro e a linha pontilhada
o intervalo HPD de 95%.
4.3
Estudo simulado
Para examinar o desempenho da metodologia proposta e a influência da distribuição a
priori nos resultados, foram feitos alguns estudos de simulação sob repetidas populações.
O objetivo é verificar o desempenho do modelo para diferentes cenários que possam
existir.
80
4.3.1
Considerando diferentes configurações
Foram geradas 500 populações considerando diferentes configurações para alguns
parâmetros variando os valores de N , R e X, assim como o nı́vel de homogeneidade/
heterogeneidade da população. Primeiramente, N foi fixado em 200, 400 e 600, e para
cada um destes valores, os valores de α e β foram fixados, com o objetivo de criar
diferentes populações raras e agrupadas e variar R e X na simulação. Em particular,
as populações foram simuladas para 4 pares de (α, β) com α, β ∈ {0.1, 0.15}. Portanto
neste estudo são apresentados resultados para 12 diferentes configurações. Neste caso,
foi considerada para λ apenas a distribuição a priori independente, portanto, para cada
população λ foi gerado a partir de uma distribuição Gama com d = 1.1 e ν = 0.13, o
que produz um CV igual a 95%. Desta maneira, estes valores fixados permitem gerar
populações raras e agrupadas com redes heterogêneas entre si. Finalmente, uma amostra
adaptativa foi selecionada de cada população, com primeiro estágio caracterizado por
uma amostra aleatória simples sem reposição de tamanho 5%N .
As Tabelas 4.2, 4.3 e 4.4 mostram um sumário com algumas propriedades frequentistas
da distribuição a posteriori dos parâmetros do modelo proposto após a convergência,
para cada configuração testada. São apresentados o EQMR, o EAR, a probabilidade de
cobertura (em porcentagem) dos intervalos HPD de 95%, com sua respectiva amplitude
média ao longo das 500 simulações. Em particular, as amplitudes dos intervalos para
T e para λs e λs̄ estão relativizadas com relação ao valor verdadeiro. Os resultados
para λj ’s estão sumarizados em relação a λs e λs̄ , pois na simulação o valor de R não
foi fixado, R foi gerado de sua distribuição, condicional ao valor de β, portanto para
cada população foi simulado um valor de R distinto, o que impede a apresentação de
propriedades frequentistas para cada λj separadamente.
No geral, é possı́vel observar que os parâmetros são bem estimados. A cobertura
dos intervalos de 95% é próxima do nı́vel nominal desejado. O EQMR e o EAR são
pequenos para a maior parte dos parâmetros, exceto para β em alguns casos especı́ficos.
Entretanto, este fato não tem um impacto significante na previsão de T , o qual é o maior
81
interesse deste trabalho. Como esperado, os resultados para λj , para j ∈ s mostram
erros menores e maior precisão do que para para j ∈ s̄.
À medida que o valor de N cresce, o EQMR e o EAR da maioria dos parâmetros
diminui. Isto ocorre porque nestes casos existe um maior número de redes não-vazias
do que para um valor menor de N , melhorando assim as estimativas de α e β, e
consequentemente de outros parâmetros. Uma sugestão de melhoria nestes casos é o
aumento do tamanho da amostra. Por outro lado, pela mesma razão, para um valor fixo
de N , os erros diminuem para valores maiores de α e β.
Tabela 4.2: Sumário a posteriori da estimação pontual e intervalar dos parâmetros do
modelo proposto e de T sob as 500 simulações, para diferentes valores de α, β e N = 200.
(α, β) = (0.10, 0.10)
(α, β) = (0.10, 0.15)
T
α
β
ν
λs
λs̄
T
α
β
ν
λs
λs̄
EQMR
0.21
0.38
0.53
0.56
0.03
0.29
0.22
0.29
0.29
0.39
0.03
0.28
RAE
0.35
0.17
0.25
0.60
0.12
0.46
0.36
0.16
0.35
0.47
0.13
0.45
Cob.
95.0
91.1
96.7
89.5
91.7
87.8
93.8
93.7
98.1
89.7
90.3
87.7
Ampl.
1.60
0.20
0.31
0.28
0.58
1.23
1.60
0.19
0.31
0.28
0.57
1.26
(α, β) = (0.15, 0.1)
(α, β) = (0.15, 0.15)
EQMR
0.09
0.20
0.50
0.22
0.02
0.31
0.06
0.10
0.19
0.32
0.02
0.27
RAE
0.24
0.31
0.45
0.40
0.11
0.46
0.21
0.27
0.21
0.47
0.10
0.41
Cob.
94.6
90.9
97.1
90.2
93.6
89.1
97.3
97.0
98.5
90.5
94.1
89.8
Ampl.
1.22
0.19
0.21
0.22
0.50
1.33
1.24
0.20
0.23
0.21
0.56
1.51
82
Tabela 4.3: Sumário a posteriori da estimação pontual e intervalar dos parâmetros do
modelo proposto e de T sob as 500 simulações, para diferentes valores de α, β e N = 400.
(α, β) = (0.10, 0.10)
(α, β) = (0.10, 0.15)
T
α
β
ν
λs
λs̄
T
α
β
ν
λs
λs̄
EQMR
0.06
0.15
0.42
0.14
0.02
0.29
0.05
0.08
0.15
0.10
0.02
0.31
RAE
0.21
0.32
0.35
0.28
0.10
0.43
0.20
0.23
0.29
0.21
0.12
0.43
Cob.
96.7
91.1
96.0
90.8
94.2
91.0
96.8
95.1
98.1
90.5
94.3
91.8
Ampl.
1.04
0.09
0.20
0.19
0.47
1.38
1.05
0.10
0.21
0.18
0.55
1.64
(α, β) = (0.15, 0.1)
(α, β) = (0.15, 0.15)
EQMR
0.04
0.06
0.35
0.04
0.02
0.30
0.05
0.03
0.15
0.03
0.02
0.36
RAE
0.18
0.18
0.39
0.18
0.09
0.42
0.20
0.15
0.21
0.15
0.10
0.43
Cob.
93.4
91.2
96.9
96.7
94.2
93.9
92.4
97.0
98.7
96.5
93.5
95.6
Ampl.
0.79
0.11
0.15
0.14
0.45
1.43
0.77
0.11
0.16
0.13
0.51
1.77
83
Tabela 4.4: Sumário a posteriori da estimação pontual e intervalar dos parâmetros do
modelo proposto e de T sob as 500 simulações, para diferentes valores de α, β e N = 600.
(α, β) = (0.10, 0.10)
(α, β) = (0.10, 0.15)
T
α
β
ν
λs
λs̄
T
α
β
ν
λs
λs̄
EQMR
0.04
0.05
0.25
0.10
0.02
0.32
0.05
0.03
0.11
0.09
0.02
0.35
RAE
0.17
0.17
0.28
0.12
0.09
0.42
0.20
0.14
0.26
0.11
0.11
0.42
Cob.
96.3
91.8
98.1
98.0
93.5
93.1
92.8
97.5
98.3
97.0
93.8
96.1
Ampl.
0.79
0.08
0.22
0.20
0.46
1.40
0.78
0.08
0.23
0.19
0.52
1.70
(α, β) = (0.15, 0.10)
(α, β) = (0.15, 0.15)
EQMR
0.05
0.04
0.21
0.06
0.01
0.37
0.09
0.08
0.06
0.05
0.02
0.35
RAE
0.19
0.17
0.30
0.09
0.09
0.44
0.29
0.24
0.18
0.09
0.10
0.43
Cob.
90.4
91.1
98.7
98.9
95.3
96.0
90.0
90.5
98.8
98.4
95.5
96.8
Ampl.
0.78
0.08
0.17
0.18
0.43
1.49
0.53
0.08
0.20
0.17
0.53
1.79
Como mencionado acima, não foi possı́vel apresentar os resultados para cada λj pois
o valor de R não foi fixado para as simulações e portanto a dimensão de λ varia em
cada simulação. Portanto na Figura 4.4 é apresentado um diagrama boxplot com o erro
relativo (ER) para λs e λs̄ para todas as redes e todas as populações, para diferentes
valores de α e β e para N = 400. Note que em todos os casos o ER está em torno de
zero e o ER para λs é menor que para λs̄ , como esperado. Além disso λs̄ é ligeiramente
subestimado com respeito à mediana da distribuição a posteriori.
4.3.2
Considerando diferentes nı́veis de heterogeneidade
As 500 populações usadas neste estudo de simulação foram geradas para alguns valores
dos parâmetros fixados, em particular foi assumido que λj segue uma distribuição Gama
com hiperparâmeros d = 1.1 e ν = 0.13. Como mencionado acima, com esses valores, esta
distribuição Gama tem um CV de aproximadamente 95%, o que geraria populações com
redes heterogêneas, com respeito à média do número de observações dentro das unidades
que as compõem. A partir de agora o interesse é avaliar o desempenho do modelo proposto
com respeito ao nı́vel de homogeneidade e heterogeneidade da população. Para realizar
84
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
ER
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
(α,β)=(0.1,0.1)
(α,β)=(0.1,0.15)
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
1
●
●
●
●
●
●
●
●
●
●
●
2
●
●
●
●
●
−0.4
0
ER
0.0 0.2 0.4 0.6
●
●
●
●
●
●
●
●
●
(α,β)=(0.1,0.1)
(α,β)=(0.1,0.15)
●
●
●
●
●
●
●
●
●
●
●
●
●
(α,β)=(0.15,0.1) (α,β)=(0.15,0.15)
(a) RE - λs
(α,β)=(0.15,0.1) (α,β)=(0.15,0.15)
(b) RE - λs̄
Figura 4.4: Erro relativo para λs e λs̄ ao longo de 500 simulações, para N = 400 e
diferentes configurações de α e β.
esta análise geramos outros dois cenários fixando o CV da distribuição Gama de λ em
50% e 25% com média fixada em 8.5. Ao fixar o valor do CV em 50% obtém-se d = 4 e
ν = 0.47 e para CV igual a 25%, d = 16 e ν = 1.89.
A Figura 4.5 apresenta as curvas das distribuições de λj para cada um dos valores de
CV fixado nesta análise. Note que à medida que o CV diminui a distribuição a priori
para λj se torna simétrica e mais concentrada em torno da média da distribuição, e
portanto as redes se tornam mais homogêneas com respeito ao total de observações em
0.10
CV=95%
CV=50%
CV=25%
0.00
Densidade
0.20
cada unidade.
0
5
10
20
30
λj
Figura 4.5: Distribuição a priori para λj usada nas simulações variando o valor do CV
da distribuição.
85
Desta forma, a análise apresentada a seguir é feita a partir de outras 500 populações
geradas fixando o CV da distribuição de λj em 50% e outras 500 com CV fixado em 25%.
Em particular, esta simulação foi feita apenas para o valor de N = 400 pois o interesse
neste era apenas verificar o desempenho do modelo variando o nı́vel de homogeneidade.
Além disso, de acordo com a Tabela 4.3 já foi visto que este valor para N apresentou
resultados razoáveis segundo as propriedades frequentistas analisadas.
Na Tabela 4.5 é apresentado novamente um sumário com algumas propriedades
frequentistas dos estimadores obtidos a partir da distribuição a posteriori para as 500
populações geradas com CV da distribuição de λj fixada em 50% e 25%. Note que ainda
para casos mais homogêneos o modelo proposto em (4.4) tem um bom desempenho,
resultando em estimadores para todos os parâmetros com pequenos EQMR e EAR
e intervalos HPD de 95% com probabilidade de cobertura próxima do nı́vel nominal
desejado.
Em particular, observe que os EQMR e EAR para T são muito semelhantes aos valores
apresentados na Tabela 4.3 quando o CV da distribuição de λj foi fixado em 95%, exceto
para o caso em que (α, β) = (0.10, 0.10), para o qual existe um número pequeno de redes
não vazias na população, e portanto um número de redes ainda menor na amostra de
5%. O EQMR e o EAR para λs̄ são menores que os observados na Tabela 4.3, apesar
dos mesmos para ν serem maiores. Além disso, à medida que o CV diminui a cobertura
empı́rica dos intervalos de 95% é subestimada, principalmente as obtidas para ν e λ. Uma
possı́vel explicação para estes resultado é que a distribuição Gama com um CV pequeno
se aproxima de uma distribuição Normal, o que parece complicar a inferência para seus
hiperparâmetros, portanto para λ e consequentemente os outros parâmetros do modelo.
Logo, uma alternativa neste caso pode ser assumir uma outra distribuição a priori para
λ. No entanto, vale destacar que ainda para esses casos estudados os EQMR e EAR são
pequenos para a maioria dos parâmetros e principalmente para o total populacional, o
qual é o maior interesse neste trabalho.
86
Tabela 4.5:
Sumário para a estimação pontual e intervalar dos parâmetros do modelo
e o total populacional para as 500 populações, variando o nı́vel de homogeneidade nas
redes, a partir do valor do CV fixado para a distribuição de λ, para N = 400.
CV = 50%
(α, β) = (0.10, 0.10)
(α, β) = (0.10, 0.15)
T
α
β
ν
λs
λs̄
T
α
β
ν
λs
λs̄
EQMR
0.13
0.15
0.52
0.16
0.02
0.04
0.06
0.09
0.18
0.10
0.02
0.03
EAR
0.26
0.32
0.27
0.30
0.10
0.15
0.18
0.24
0.36
0.23
0.11
0.15
Cob.
95.3
87.2
97.0
95.3
94.7
97.0
96.7
95.0
98.2
95.0
94.5
97.6
Ampl.
1.38
0.11
0.26
0.91
0.51
1.27
1.24
0.11
0.27
0.82
0.55
1.31
(α, β) = (0.15, 0.1)
(α, β) = (0.15, 0.15)
EQMR
0.03
0.04
0.40
0.08
0.02
0.03
0.03
0.03
0.10
0.06
0.02
0.03
EAR
0.15
0.15
0.50
0.21
0.10
0.12
0.16
0.14
0.26
0.18
0.10
0.13
Cob.
96.5
94.7
97.3
97.8
95.6
98.0
95.8
97.3
98.0
97.5
95.8
97.9
Ampl.
0.95
0.11
0.23
0.75
0.48
1.28
0.92
0.11
0.24
0.70
0.53
1.36
CV = 25%
(α, β) = (0.10, 0.10)
(α, β) = (0.10, 0.15)
EQMR
0.09
0.30
0.50
0.36
0.03
0.08
0.05
0.18
0.12
0.34
0.03
0.08
EAR
0.23
0.48
0.37
0.47
0.13
0.24
0.19
0.37
0.29
0.44
0.14
0.26
Cob.
89.7
86.8
98.0
75.0
79.7
90.1
94.7
90.1
98.2
74.9
74.5
90.0
Ampl.
0.96
0.12
0.25
3.01
0.47
0.70
0.91
0.12
0.27
2.83
0.51
0.75
(α, β) = (0.15, 0.1)
(α, β) = (0.15, 0.15)
EQMR
0.03
0.08
0.41
0.25
0.02
0.03
0.04
0.05
0.07
0.19
0.02
0.04
EAR
0.14
0.22
0.49
0.34
0.10
0.15
0.17
0.15
0.21
0.24
0.11
0.17
Cob.
96.6
91.7
97.5
80.8
84.6
94.4
91.9
92.5
98.3
83.2
83.8
93.9
Ampl.
0.70
0.12
0.22
2.48
0.46
0.74
0.70
0.12
0.23
2.25
0.50
0.79
87
4.3.3
Análise de sensibilidade da distribuição a priori
O interesse agora é comparar o desempenho do modelo sob uma outra alternativa de
distribuição a priori para λ usada na literatura, que é a distribuição a priori dependente.
Neste caso, o estudo de análise de sensibilidade é feita com a geração de 500 populações
e uma amostra adaptativa de tamanho inicial de 5% selecionada de cada uma. Como
o maior interesse está na comparação da influência de ambas as distribuições a priori
para λ nos resultados, foi escolhido efetuar a análise para somente alguns valores fixos de
R, a fim de viabilizar a apresentação dos resultados para cada λj separadamente e não
somente para para os parâmetros sumarizados em λs e λs̄ . Em particular, para gerar as
populações usadas neste estudo, fixou-se N = 400 e buscou-se uma configuração para α e
β que produz uma população rara e agrupada que tenha apresentado um bom desempenho
no estudo simulado realizado na Subseção 4.3.1. Portanto, fixou-se (α, β) = (0.15, 0.10)
e gerou-se um grande número de populações até obter 500 populações com R = 5, outras
500 com R = 6 e 500 com R = 7. Estes valores de R foram escolhidos porque, de acordo
com sua distribuição no modelo, dada pela equação em (4.4e), estes são valores com alta
probabilidade para (α, β) = (0.15, 0.10). Finalmente, como estão sendo especificadas
neste estudo duas distribuições a priori para λ, na geração dos dados fixou-se para todas
as populações λ em um valor arbitrário gerado de uma distribuição Uniforme definida
no intervalo (3,15).
Todos os resultados apresentados a seguir correspondem a 1.000 amostras
independentes da distribuição a posteriori do vetor paramétrico, geradas de 200.000
iterações do RJMCMC, com um aquecimento de 10.000 e um espaçamento de 190. As
mesmas distribuições a priori usadas para α e β no estudo de simulação apresentado na
Subseção 4.3.1. Para λ foi considerada então a distribuição a priori Gama e a Normal
truncada dependente com desvio padrão τ ∈ {1, 5, 10, 20}, ambas descritas na Seção 4.2.1
Primeiramente, para uma única população gerada ajustamos o modelo proposto com
as distribuições a priori a fim de visualizarmos de forma preliminar o desempenho das
distribuições a priori consideradas e se os valores de τ considerados eram razoáveis. Dessa
forma, foi visto que o maior impacto da distribuição a priori de λ está na distribuição
88
a posteriori de R. A Figura 4.6 apresenta o intervalo HPD de 95% obtido para R para
cada distribuição a priori de λ considerada. Note que a distribuição a posteriori de R
é altamente imprecisa quando τ = 1, no entanto à medida que o valor de τ assumindo
aumenta este comportamento melhora, e com τ = 20 tem-se inclusive um comportamento
similar ao obtido no caso da distribuição a priori independente. Logo, deste momento
em diante decidiu-se descartar a distribuição a priori Normal truncada dependente com
τ = 1.
_
15
R
20
25
_
_
_
5
10
_
●
●
●
●
●
_
_
_
_
_
Indep
τ=1
τ=5
τ = 10
τ = 20
Distribuições a priori
Figura 4.6: Sumário da distribuição a posteriori de R assumindo diferentes distribuições
a priori para λ. As cruzes representam a mediana da distribuição a posteriori, o cı́rculo
o valor verdadeiro de R e a linha o intervalo HPD de 95%.
Para as 500 populações geradas foi feita uma análise de sensibilidade com respeito
à distribuição a posteriori de cada λj . A Figura 4.7 apresenta o EQMR para cada
λj em ordem crescente, mas separados para as amostras em que a rede j é observada
(a) e quando não é (b). Os resultados com a distribuição a priori independente são
representados pelos cı́rculos vazios e pela linha cheia, já os resultados para a distribuição
dependente com τ = 5 são representados pelos triângulos e a linha tracejada, já as cruzes
e linha pontilhada representam τ = 10 e τ = 20 é representado pelos cı́rculos cheios e a
linha traço e ponto. Pela Figura 4.7 (a) é possı́vel concluir que a distribuição a priori
independente produz na maioria dos casos EQMR menor que a distribuição dependente,
principalmente para os λj ’s com valor absoluto menor. Os resultados mostram-se muito
similares para os diferentes valores de τ . Para λj para o caso em que a rede j não pertence
à amostra o EQMR é maior do que quando j pertence à amostra, como esperado, e neste
89
caso os resultados sob cada distribuição considerada tornam-se mais similares entre si
0.00
●
●
λ1
λ2
λ3
●
●
●
●
λ4
λ5
●
λ2
●
●
●
●
λ3
λ4
●
●
●
●
λ5
λ6
●
0.04
0.06
●
●
λ1
EQMR
●
0.00
EQMR
●
0.00
0.02
EQMR
●
0.02
0.04
●
●
0.02
0.04
para maiores valores de R.
●
●
●
●
●
●
●
●
●
●
●
●
λ1 λ2 λ3 λ4 λ5 λ6 λ7
●
●
●
●
0.0
0.0
●
●
●
λ1
λ2
λ3
λ4
λ5
λ1
λ2
●
●
●
λ3
λ4
●
●
0.4
●
●
●
EQMR
●
●
●
●
●
0.0
●
0.8
0.6
EQMR
●
0.2
●
0.4
●
0.4
EQMR
0.8
(a) RMSE for λj , j ∈ s
λ5
λ6
●
●
●
●
●
●
●
λ1 λ2 λ3 λ4 λ5 λ6 λ7
(b) RMSE for λj , j ∈ s̄
Figura 4.7: EMQR para cada λj assumindo diferentes distribuições a priori para λ. Os
resultados com a distribuição a priori independente são representados pelos cı́rculos vazios
e a linha cheia, os resultados para a distribuição dependente com τ = 5 são representados
pelos triângulos e a linha tracejada, as cruzes com a linha pontilhada representam os
resultados quando τ = 10 e τ = 20 são os cı́rculos cheios e a linha traço e ponto.
Finalmente, como prever o total populacional é o maior interesse neste trabalho, foi
avaliado também o impacto destas distribuições a priori na distribuição a posteriori de
T . Na Figura 4.8 estão os EQMR de T para cada valor de R considerado, a cobertura
dos intervalos HPD de 95% e sua respectiva amplitude média relativa. Os cı́rculos vazios
e a linha representam neste caso os resultados para R = 5, os triângulos com a linha
tracejada para R = 6 e as cruzes com a linha pontilhada para R = 7. Observe que o
EQMR obtido no caso em que se assume uma distribuição a priori Gama independente
para λ é sempre maior do que quando se assume a distribuição dependente. No geral, os
90
intervalos de 95% apresentam maior probabilidade de cobertura do que o nı́vel desejado e
no caso da hipótese de independência a priori estes são mais precisos, quando comparados
aos obtidos sob hipótese de dependência a priori. Note que condicional ao valor de R,
●
96
Amplitude
●
●
92
●
●
94
●
Cobertura
0.04
●
0.02
EQMR
●
Indep
τ=5
τ = 10
τ = 20
Distribuições a priori
Indep
τ=5
τ = 10
τ = 20
Distribuições a priori
0.6 0.8 1.0 1.2 1.4
98
0.06
os resultados para a distribuição dependente são bastante similares quando varia-se τ .
●
●
τ=5
τ = 10
●
●
Indep
τ = 20
Distribuições a priori
Figura 4.8: EQMR, probabilidade de cobertura e amplitude média do intervalo HPD de
95% para o total populacional T sob cada distribuição a priori assumida para λ e para
cada valor de R fixado. Os cı́rculos vazios e a linha representam os resultados para R = 5,
os triângulos com a linha tracejada quando R = 6 e as cruzes com a linha pontilhada
para R = 7.
Portanto, sob alguns critérios considerar uma distribuição a priori independente
parece ser mais eficiente que a distribuição a priori dependente e vice-versa. No entanto,
vale destacar que a distribuição dependente é mais fácil de interpretar, o que torna sua
elicitação mais intuitiva em muitos casos, em que não há conhecimento a priori adequado
sobre a população.
4.4
Comparação com o modelo agregado
O modelo de mistura em (4.4) foi proposto neste trabalho como uma alternativa ao
modelo agregado, principalmente quando não é adequada a suposição de homogeneidade
entre redes com respeito ao número de observações dentro destas e quando o número
esperado de observações não é proporcional à sua área. O objetivo deste modelo proposto
é, portanto, aprimorar as estimativas populacionais obtidas com o ajuste do modelo
91
agregado através do uso de um modelo que leve em conta na sua formulação a suposição
de heterogeneidade entre redes. Isto é realizado na proposta através da modelagem no
nı́vel das unidades primárias, no lugar das redes.
Para acessar a eficiência da metodologia proposta, nesta seção é feita uma comparação
do desempenho do modelo de mistura com o modelo agregado em duas situações. A
primeira comparação consiste de um experimento de simulação baseado no desenho
amostral com uma população real, já o outro estudo é baseado em simulações sob o
modelo.
Para ajustar ambos os modelos, foram assumidas as mesmas distribuições a priori
usadas na Subseção 4.3.
Na execução dos métodos de MCMC e RJMCMC foram
realizadas 200.000 iterações cada, 10.000 foram descartadas como aquecimento da cadeia
e as amostras finais foram tomadas de 190 em 190, a fim de obter 1.000 amostras
independentes.
4.4.1
Simulação baseada no desenho amostral
O estudo apresentado a seguir baseia-se em verificar propriedades frequentistas dos
estimadores obtidos do ajuste de cada modelo, a partir da seleção de várias amostras
de uma população real. Tal população está descrita na Seção 3.1.3 e é composta por
marrecos da asa azul na região da Flórida, Estados Unidos, no ano de 1992.
O estudo consiste em selecionar 500 amostras adaptativas com tamanho inicial de
10%N desta população real. Note que esta população, a qual pode ser vista na Figura
3.2, é composta por 3 principais redes, as quais apresentam no geral um número médio de
marrecos diferente para cada rede. E o total em cada rede não é proporcional ao número
de unidades em cada uma. Logo, as hipóteses do modelo agregado não seriam adequadas
a este conjunto de dados. Por outro lado, o modelo de mistura assume heterogeneidade
entre redes, o que parece mais razoável ao observar a Figura 3.2.
Além disso, observe que existem duas unidades com um número discrepante de
marrecos da asa azul, logo se as amostras selecionadas não contivessem estas unidades,
seria extremamente difı́cil estimar o total populacional próximo do valor verdadeiro.
92
Portanto, optou-se por fixar esta rede na amostra, de modo que a probabilidade de
seleção desta fosse igual a 1.
A Figura 4.9 apresenta os traços das cadeias com a distribuição a posteriori de α,
β e T , partindo de dois pontos iniciais distintos, sob ambos os modelos para uma das
amostras selecionadas. A linha cinza representa o valor verdadeiro do total. Note que o
modelo agregado tende a sobreestimar o total populacional, um comportamento esperado
neste caso devido a heterogeneidade presente nos dados. Além disso, observe que o
parâmetro α é estimado num valor mais alto quando ajustado o modelo agregado que
quando ajustado o modelo de mistura. Como este parâmetro está relacionado com o
número de unidades não-vazias, o modelo agregado estima um número maior de unidades
0
200
600
T
1000
0
200
iterações
600
14200 14400 14600
0.0
0.02
0.2
β
α
0.4
0.06
0.6
não vazias na população que o modelo de mistura.
1000
0
200
iterações
600
1000
iterações
T
14200
0.0
0
200
600
1000
0
iterações
14800
0.4
0.2
β
0.15
0.05
α
15400
(a) Modelo de mistura
200
600
1000
iterações
0
200
600
1000
iterações
(b) Modelo agregado
Figura 4.9: Traço das cadeias com a distribuição a posteriori para α, β e T obtida do
ajuste do modelo de mistura (a) e do modelo agregado (b). A linha em cinza representa
o valor verdadeiro de T .
A Tabela 4.6 apresenta os valores da estatı́stica de Geweke e do fator de dependência
do diagnóstico de Raftery-Lewis.
Sob ambos os critérios é possı́vel observar que
93
a convergência foi alcançada.
Esta mesma conclusão vale para todas as amostras
selecionadas.
Tabela 4.6:
Análise da convergência das cadeias com a distribuição a posteriori dos
parâmetros dos modelos de mistura e agregado para a população real.
Geweke
Mistura
Raftery-Lewis
Agregado Mistura
Agregado
α
0.54
-0.12
1.05
0.95
β
-0.75
-0.93
1.03
1.01
T
0.32
0.27
1.55
1.08
Na Tabela 4.7 apresenta-se uma comparação com base em propriedades frequentistas
dos estimadores obtidos para o total populacional T , sob os dois modelos. Como temos
a população inteira de marrecos da asa azul e selecionamos amostras desta população, é
possı́vel usar critérios de comparação entre os modelos que usam o valor verdadeiro. Logo
são apresentados o EQMR, o EAR, as probabilidades de cobertura dos intervalos HPD
de 95%, a média das amplitudes relativas destes intervalos, expressa pela razão entre
seu valor e o valor verdadeiro de T . É apresentado também a eficiência do estimador do
total obtido do ajuste do modelo de mistura com relação ao modelo agregado, sob as 500
amostras. Observe que o modelo proposto apresenta menores valores para os EQMR e
EAR, além de probabilidade de cobertura dos intervalos mais próxima do nı́vel nominal
desejado de 95%, mesmo tendo uma amplitude menor, portanto os intervalos gerados
sob esta abordagem são mais precisos. Além disso, como a eficiência é menor que 1,
isto indica que a variância do estimador para o total sob o modelo (4.4) proposto neste
trabalho é menor que sob o modelo agregado.
Finalmente, na Figura 4.10 é apresentado um diagrama boxplot com os ER para o
total populacional para as 500 amostras sorteadas sob os dois modelos em questão. Note
que os ER obtidos com base no modelo de mistura são inferiores aos obtidos com o
modelo agregado, apesar de ambos no geral sobrestimarem o valor verdadeiro de T .
94
Tabela 4.7: Sumário da estimação pontual e intervalar do total populacional obtido do
ajuste do modelo de mistura e do modelo agregado.
EQMR
EAR
Cobertura Amplitude ef(T̂ )
Modelo proposto
0.02
0.05
97.7
0.45
Modelo agregado
0.05
0.18
84.5
0.63
0.15
−0.05
0.05
ER
0.25
0.78
Modelo de mistura Modelo agregado
Figura 4.10: ER para T para as 500 amostras obtidos a partir do ajuste do modelo de
mistura e do modelo agregado.
4.4.2
Simulação baseada no modelo
O interesse neste estudo é avaliar o desempenho do modelo agregado sob populações
mais homogêneas simuladas do modelo de mistura em (4.4). Portanto, para as mesmas
500 populações usadas no estudo simulado apresentado na Subseção 4.3.1 ajustou-se
o modelo agregado.
Em particular este estudo destina-se as simulações realizadas
assumindo (α, β) = (0.15, 0.10) e a distribuição Gama para λ com CV=50% e 25%,
que caracterizam populações com redes mais homogêneas.
Na Tabela 4.8 são apresentadas propriedades frequentistas dos estimadores obtidos
do ajuste do modelo agregado. A fim de facilitar a comparação, os resultados obtidos
do ajuste do modelo de mistura com as mesmas 500 populações são apresentados em
parênteses na tabela. Os resultados para T indicam um maior EAR e EQMR no ajuste
do modelo agregado, mas essa diferença parece diminuir à medida que o CV diminui. Por
outro lado, o estimador para β produzido no ajuste do modelo agregado apresenta para
95
todos os casos um menor EAR e EQMR do que o estimador obtido do ajuste do modelo
de mistura. Portanto, de acordo com este critério, é possı́vel concluir que à medida que o
nı́vel de heterogeneidade aumenta os resultados são favoráveis ao modelo de mistura com
relação a previsão de T , e para populações mais homogêneas os resultados tendem a se
tornar mais semelhantes. Entretanto, com relação a estimação de β o modelo agregado
apresenta um melhor desempenho.
Tabela 4.8: Sumário a posteriori para a estimação pontual e intervalar dos parâmetros
dos modelos sob as 500 simulações onde λ foi gerado de uma distribuição Gama com
CV=50% e CV=25%, para N = 400 e (α, β) = (0.15, 0.10).
CV=50%
CV=25%
T
α
β
T
α
β
EQMR
0.05 (0.03)
0.04 (0.04)
0.10 (0.40)
0.03 (0.03)
0.05 (0.08)
0.18 (0.41)
EAR
0.21 (0.15)
0.19 (0.15)
0.37 (0.50)
0.17 (0.14)
0.16 (0.22)
0.32 (0.49)
Cob.
95.6 (96.5)
98.1 (94.7)
97.4 (97.3)
96.8 (96.6)
97.1 (91.7)
95.6 (97.5)
Ampl.
0.85 (0.95)
0.16 (0.11)
0.18 (0.23)
0.86 (0.70)
0.16 (0.12)
0.19 (0.22)
Finalmente, na Figura 4.11 são apresentados os diagramas boxplot com o ER para
T sob ambos os modelos. Note que um maior ER é obtido quando ajustado o modelo
agregado, em particular T é subestimado se a mediana dos ER é observada. Entretanto,
este comportamento tende a ser atenuado à medida que o grau de homogeneidade
aumenta.
Portanto, a partir destes resultados pode-se concluir que à medida que o nı́vel de
heterogeneidade entre as redes diminui, o desempenho dos modelos torna-se similar, com
relação a previsão de T , o qual é o maior interesse neste trabalho. A principal diferença
seria o número de parâmetros a estimar e o custo computacional na implementação dos
métodos de aproximação necessários no ajuste de cada modelo.
96
0.3
0.4
0.1
0.2
−0.1
ER
0.0
ER
−0.3
−0.2
−0.4
Modelo de mistura Modelo agregado
Modelo de mistura Modelo agregado
(a) CV = 50%
(b) CV = 25%
Figura 4.11: Boxplot com o ER para T , a partir do modelo de mistura e do modelo
agregado para as 500 populações, tal que λ foi gerado de uma distribuição Gama com
CV=50% e CV=25%.
4.5
Modelo de mistura sob amostragem adaptativa
dupla
Apesar do planejamento amostral adaptativo por conglomerados mostrar-se
apropriado em levantamentos cuja população-alvo se comporta de forma rara e agrupada,
uma de suas principais desvantagens é a impossibilidade de controlar o tamanho da
amostra final. Neste sentido, algumas alternativas surgiram na literatura visando impor
um limite a este tamanho final para amostras coletadas de forma adaptativa. Neste
trabalho, temos particular interesse na abordagem de Felix-Medina e Thompson (2004),
chamada amostragem adaptativa dupla por conglomerados.
O interesse agora está em aplicar o modelo de mistura (4.4) à populações raras e
agrupadas, cujas amostras são provenientes do planejamento amostral elaborado por
Felix-Medina e Thompson (2004). Com essa mudança, algumas adaptações devem ser
feitas no modelo proposto. A probabilidade de seleção dada em (4.5) deve ser recalculada
e o método de inferência reescrito. Veremos também que, sob algumas condições, a
amostragem adaptativa por conglomerados pode ser tratada como um caso particular da
amostragem dupla. Além disso, como com este desenho é possı́vel aumentar o tamanho
97
da amostra e usar informações auxiliares, espera-se uma melhora na qualidade das
estimativas dos parâmetros do modelo e do total populacional sem exceder abusivamente
os custos disponı́veis.
4.5.1
Amostragem adaptativa dupla
Proposto por Felix-Medina e Thompson (2004), este plano amostral trata-se de uma
variação com múltiplos estágios da amostragem adaptativa por conglomerados. Chamado
amostragem adaptativa dupla, o método permite ao pesquisador atingir aos seguintes
objetivos: controlar o número de observações da variável de interesse; alocar a amostra
final próxima a locais interessantes; e utilizar uma variável auxiliar na estimação do
parâmetro populacional de interesse.
A metodologia pode ser decomposta em três estágios e está descrita a seguir. Seja
H uma variável auxiliar menos custosa que a variável de interesse e mais fácil de medir.
Suponha que nada se conhece sobre os valores desta variável auxiliar antes do inı́cio da
coleta da amostra.
A primeira fase do método consiste em selecionar uma amostra adaptativa por
conglomerados s1 baseada nos valores da variável auxiliar H, gerando m1 diferentes
redes, vazias e não-vazias.
A segunda fase consiste em selecionar uma subamostra s2 de m2 redes das m1
diferentes redes que estão na amostra s1 . Esta seleção pode ser feita segundo planos
amostrais probabilı́sticos convencionais.
Finalmente, a terceira fase consiste em selecionar uma subamostra de unidades
primárias dentro de cada uma das redes em s2 e observar o valor da variável de interesse
Y associada em cada uma destas. Denote por s3i (i = 1, . . . , m2 ) a amostra de unidades
S 2
observada na rede i, cujo tamanho é dado por n3i , e portanto, s3 = m
i=1 s3i .
Segundo Felix-Medina e Thompson (2004), existem várias possibilidades de variações
dentro destas três fases.
Uma destas é omitir a segunda fase e subamostrar toda
rede em s1 . Cada rede pode ser subamostrada antes mesmo do pesquisador terminar
o planejamento s1 .
Neste último caso, há necessidade em controlar o tamanho da
amostra antes de iniciar as outras fases. Outra possibilidade é combinar diferentes planos
98
probabilı́sticos para selecionar s2 e s3i (i = 1, . . . , m2 ). A maioria das combinações
permite ao pesquisador um controle sobre custos e número de medidas da variável de
interesse.
Além disso, com relação às variáveis auxiliares, podem ser usadas variáveis quaisquer
correlacionadas com a variável de interesse e mais fáceis de medir, ou ainda por exemplo
variáveis de avaliação rápida, as quais conduzem o pesquisador para as áreas mais
promissoras, onde observações exatas da variável podem ser feitas posteriormente. Por
exemplo, numa pesquisa sobre mexilhões de água doce, cujo interesse é estimar o total
de mexilhões numa região, a variável de interesse, ou seja o número de mexilhões, é uma
variável difı́cil de ser medida porque alguns mexilhões são parcialmente escondidos pela
areia e pedras no fundo do rio. Desta forma pode-se recorrer a amostragem adaptativa
dupla, com primeiro estágio caracterizada por uma amostra adaptativa somente para
detectar a presença ou ausência de mexilhões, e esta ser usada como uma variável auxiliar
no método.
4.5.2
Modelo proposto sob amostragem dupla com variável
auxiliar indicadora de presença
O modelo de mistura em (4.4) deve ser ajustado a populações raras e agrupadas, as
quais são amostradas de forma adaptativa. Por outro lado, como o desenho amostral
adaptativo por conglomerados é informativo, à verossimilhança completa do modelo
(4.4) acrescenta-se a probabilidade de inclusão da amostra, dada em (4.5).
Neste
momento a ideia é substituir este desenho amostral, pelo proposto por Felix-Medina
e Thompson (2004). Esta pequena mudança traz adaptações na verossimilhança, por
conta da probabilidade de inclusão, e em alguns aspectos do procedimento de inferência,
os quais serão descritos a seguir.
Assim como no exemplo do mexilhão, há particular interesse em uma variável auxiliar
H binária, que assume o valor 1 se há ao menos uma observação de interesse, ou seja se
Yi > 0, e 0 caso contrário. Além disso, suponha que s2 e s3i , (i = 1, . . . , m2 ) são sorteadas
99
segundo um desenho amostral aleatório simples. Este estudo será restrito a um plano
amostral adaptativo duplo com estas caracterı́sticas.
Desta forma, a amostra final s é composta pelas unidades que compõem s1 e s3 .
Ou seja, pelas m1 redes amostradas de forma adaptativa na primeira fase e pelas n3i ,
i = 1, . . . , m2 unidades selecionadas dentro das m2 redes amostradas no segundo estágio.
Note que de s1 só se extrai informações acerca da estrutura das redes, sem observar Y
dentro destas. Enquanto que de s3i , para i = 1, . . . , m2 , se extrai informações acerca da
variável de interesse Y dentro das unidades primárias selecionadas. Por esse motivo s é
caracterizada pela união de s1 e s3 .
Portanto, ao selecionar uma amostra adaptativa dupla as informações observadas
surgem em etapas. Na primeira fase, a amostragem adaptativa com a variável auxiliar do
tipo presença/ ausência, fornece informações acerca das variáveis X, R e C. Portanto,
de s1 tem-se Xs , Rs e Cs no modelo (4.4). O segundo estágio não fornece nenhuma
informação a mais sobre as variáveis do modelo. Finalmente, na terceira fase uma parte
da variável de interesse Y é observada, ou seja Ys , o qual neste caso indica os totais
observados em uma subamostra de unidades de uma subamostra de redes não-vazias.
Portanto, ao aplicar este planejamento amostral ao modelo proposto, este continua
com a mesma estrutura descrita em (4.4). Entretanto, a probabilidade de seleção de uma
amostra s deve ser revista, pois o planejamento amostral foi alterado. Em particular,
à probabilidade de inclusão dada em (4.5), devem ser acrescentadas a probabilidade de
inclusão de s2 e s3 . Em particular, neste caso, em que consideramos s2 e s3 selecionadas
aleatoriamente, esta probabilidade é obtida da seguinte forma:
m1
Y
m
Y2
zi × gil ,l
1
[s | X, R, C] =
×
×
PN −X+Rl
Pj−1
m
−
(h
−
1)
z
−
z
1
i
i
k
i=1
k=0
l=1
h=1
n3h
m2 Y
Y
1
.
×
C
h − (i − 1)
i=1
h=1
(4.6)
O segundo termo da multiplicação na equação em (4.6) refere-se justamente à amostra
s2 , e é a probabilidade de seleção de m2 redes dentre m1 sob amostragem aleatória simples
sem reposição. O terceiro fator refere-se à amostra s3 , ou seja é a probabilidade de seleção
de n3h unidades, h = 1, . . . , m2 , dentro das m2 redes observadas na segunda fase. Observe
100
que como os planos amostrais da segunda e terceira fases constituem-se de amostragem
aleatória simples, os quais são desenhos ignoráveis, estes não fornecem informação a mais
para a previsão das variáveis não observadas. A única parcela que depende das variáveis
não observadas vem da expressão em (4.5), logo as outras parcelas são constantes na
distribuição a posteriori.
4.5.2.1
Inferência
O procedimento de inferência baseia-se na obtenção da distribuição a posteriori para
o vetor paramétrico Θ = (Xs̄ , Rs̄ , s̄ , Cs̄ , Ys̄ , Ys∩s̄3 , α, β, λ). Note que, à primeira vista,
a diferença entre aplicar o modelo a este planejamento ou ao anterior está na inserção de
Ys∩s̄3 . Pois neste caso, além da previsão de Yi para as unidades i ∈ s̄, também devem
ser preditos Yi para as unidades i que apesar de fazerem parte da amostra s, não foram
observadas em s3 e portanto são desconhecidas, ou seja, para i ∈ s ∩ s̄3 . Uma vantagem
é que, com este plano amostral menos custoso, a amostra s pode aumentar, portanto s̄
diminui e, portanto a dimensão do vetor paramétrico diminui. Esta e outras diferenças
serão apresentadas a seguir.
Note que, diferente da amostragem adaptativa por conglomerados, o atual
planejamento induz uma nova partição, de Y, tal que Y = (Ys3 , Ys∩s̄3 , Ys̄ )0 . Note
que apesar de usarmos a notação de s para as unidades que pertencem a amostra, como
a amostra é formada pela união de subamostras e apenas em s3 é que valores de Y
são observados, Ys3 é a única parte conhecida de Y e portanto Ys∩s̄3 e Ys̄ devem ser
preditos. A diferença entre estes dos últimos é que existem informações adicionais sobre
a estrutura das redes que contém as unidades em s ∩ s̄3 , o que auxilia na previsão de
Ys∩s̄3 , melhorando assim a qualidade das previsões dos totais nestas unidades, quando
comparado a s̄. Portanto, no processo de inferência com base na obtenção da distribuição
a posteriori, é necessário incluir às distribuições condicionais completas do Apêndice B a
distribuição de Ys∩s̄3 . Dessa forma a expressão em (2.2) dada no Apêndice B é reescrita
da seguinte maneira:
101
[Ys∩s̄3 , Ys̄ | ·] ∝

 Y


Y
{j:j∈Λ} {i:i =j}
λYj i  
Yi !  

Y
Y
λYj i 
{j:j∈s2 } {i∈s̄3 :i =j}
Yi ! 
,
tal que Λ = s̄ ∪ {s1 ∩ s̄2 }.
Com relação a estimação de λ também existe uma diferença.
O atual desenho
amostral induz a uma partição deste parâmetro um pouco diferente da obtida quando
se realiza somente a amostragem adaptativa por conglomerados em um único estágio.
No caso da amostragem dupla terı́amos uma partição da forma λ = (λs2 , λs1 ∩s̄2 , λs̄ )0 ,
onde λs2 está associado às redes que foram amostradas em s2 e portanto apresentam
informação adicional Y para algumas unidades que as compõem, λs1 ∩s̄2 às redes que foram
amostradas em s1 , mas que não fazem parte de s2 , e λs̄ continua se referindo a parte de
λ associada às redes não amostradas, sequer no primeiro estágio. Observe a distribuição
condicional completa de λ na equação (2.1) no Apêndice B, esta depende das variáveis
Y e C, logo quanto maior o conhecimento acerca destas variáveis, melhor a estimação
deste parâmetro. Portanto, espera-se que λs2 seja o parâmetro melhor estimado, pois
além do conhecimento de uma parte de C proveniente de s1 , s3 fornece adicionalmente
informações sobre Y para as redes selecionadas em s2 . Por outro lado, λs1 ∩s̄2 deve ser o
segundo melhor estimado pois para as redes em s1 ∩ s̄2 há apenas o conhecimento de uma
parte de C. Finalmente, o subvetor λs̄ continua sendo o mais difı́cil de ser estimado, por
falta de informação.
Portanto, como este planejamento amostral permite aumentar o número de
observações com um custo controlado, espera-se melhorar a estimação de parâmetros e
previsão de quantidades populacionais que apresentaram alguma dificuldade. Isso porque
com este método é possı́vel diminuir o número de redes não-vazias para as quais não se tem
nenhum conhecimento. Com o desenho amostral construı́do em 3 estágios, é possı́vel ao
menos conhecer para algumas redes o tamanho destas, mesmo sem observar diretamente
a variável de interesse Y . Inclusive esta foi a maior motivação para estendermos o modelo
(4.4) para um plano amostral alternativo que extraı́sse maiores informações da população,
sem extrapolar os custos operacionais. Neste caso escolheu-se a amostragem adaptativa
dupla, com variável auxiliar do tipo ausência/ presença da caracterı́stica de interesse.
102
4.5.3
Avaliação
do
modelo
proposto
sob
amostragem
adaptativa e adaptativa dupla
Este estudo baseia-se na avaliação do modelo de mistura proposto em (4.4) quando
se considera os dois planejamentos amostrais estudados neste trabalho: amostragem
adaptativa por conglomerados e a amostragem dupla. Note que neste particular estudo
optou-se por não utilizar a população real de marrecos da asa azul, descrito na Subseção
3.1.3, pois seu tamanho é relativamente pequeno para fins desta comparação. Portanto,
foram geradas 500 populações com N = 600 unidades, X = 15%N unidades não-vazias
e R = 10%X = 9 redes não-vazias, e de cada uma destas foram simuladas as seguintes
amostras:
(i) adaptativa por conglomerados com tamanho inicial n1 = 10%N produzindo m1
redes na amostra;
(ii) adaptativa dupla por conglomerados com tamanho inicial n1 = 10%N produzindo
m1 redes na amostra e
(a) m2 = 100%m1 e n3i = 70%Ci , i = 1, . . . , m2 ;
(b) m2 = 70%m1 e n3i = 100%Ci , i = 1, . . . , m2 ;
(c) m2 = 70%m1 e n3i = 70%Ci , i = 1, . . . , m2 .
O interesse é comparar o ajuste do modelo sob estes quatro planejamentos. Os
cenários (ii-a), (ii-b), (ii-c) tratam-se de variações do plano amostral duplo. Observe que,
apesar do cenário (ii-b) estar caracterizado como uma amostragem adaptativa dupla,
este também pode ser tratado como o planejamento (i), porém com um menor tamanho
inicial de amostra.
Para este estudo, foi utilizada a mesma distribuição a priori usada na Subseçção 4.3,
supondo a distribuição a priori para λ independente. Após 200000 iterações, com um
burn-in de 10000 e espaçamento de 190, foram obtidas 1000 amostras independentes da
distribuição a posteriori do vetor paramétrico Θ. Para todos os parâmetros observou-se
a convergência.
103
Na Tabela 4.9 estão os EQMR, EAR, probabilidade de cobertura do intervalo
HPD de 95% e sua respectiva amplitude média relativizada para a previsão do total
populacional T . Note que para todos os planejamentos temos erros pequenos e intervalos
HPD com probabilidade de cobertura próxima do nı́vel desejado de 95%. Mesmo no
planejamento (ii-c), em que se reduz de forma mais significante o tamanho da amostra
quando comparado aos demais, tem-se resultados que mostram boas previsões neste caso.
Portanto, mesmo com um número menor de observações da variável de interesse é possı́vel
obter resultados tão eficientes quanto os obtidos usando a amostragem adaptativa em um
estágio.
Tabela 4.9: Sumário a posteriori do total populacional T para os quatro planejamentos
considerados com base nas 500 amostras simuladas.
Amostra
EQMR
EAR
Cobertura (%)
Amplitude relativa
(i)
0.02
0.12
96.0
0.61
(ii-a)
0.03
0.14
95.9
0.62
(ii-b)
0.02
0.12
93.3
0.69
(ii-c)
0.03
0.13
95.8
0.62
Com relação aos planos amostrais (i) e (ii-a) a diferença está no número de unidades
que são observadas dentro das redes amostradas. O segundo observa um número menor
de unidades com relação a variável de interesse, portanto em contextos em que observar
Y é altamente custoso, pode-se preferir o plano (ii-a). Desta forma, o interesse agora
concentrar-se-á em comparar a performance do modelo de mistura sob estes dois planos
em particular. Quando comparados ambos os planejamentos com relação a previsão do
total populacional T , não foram observadas grandes diferenças, portanto com base neste
critério ambos mostraram-se eficientes. Portanto, será feita uma comparação de ambas
as metodologias a partir da estimação do parâmetro λs2 . A ideia em usar λ como critério
de comparação dá-se pois este é um parâmetro importante para a previsão do total e está
relacionado diretamente com as informações extraı́das dentro das redes.
104
Na Figura 4.12 está um sumário da distribuição a posteriori de λ ao longo das
500 simulações.
Nesta são apresentados o EAR, a probabilidade de cobertura do
intervalo HPD de 95% e sua respectiva amplitude média relativizada com relação ao valor
verdadeiro. O triângulo com linha cheia representa o plano amostral (i) e o cı́rculo cheio
com linha pontilhada o plano amostral (ii-a). Note que o plano (i) produz erros relativos
ligeiramente menores que o plano (ii-a), o que era de se esperar pois este apresenta
um maior tamanho de amostra final. Além disso, os intervalos HPD de 95% são mais
precisos para todos os λj s sob o plano amostral (i). Com relação as probabilidades de
cobertura não há nada conclusivo sobre qual plano é mais eficiente, ora um se apresenta
mais próximo do nı́vel desejado, ora outro se apresenta. Observe que λ6 apresenta
uma subestimação da probabilidade de cobertura, mas este fato ocorre para os dois
planejamentos em questão.
●
●
●
λ1 λ2 λ3 λ4 λ5 λ6 λ7 λ8 λ9
●
●
●
0.9
0.7
●
●
●
●
●
0.5
●
●
●
●
●
●
λ1 λ2 λ3 λ4 λ5 λ6 λ7 λ8 λ9
0.3
●
●
●
Amplitude
●
0.92
●
●
Cobertura
●
0.96
●
0.88
0.03
0.01
EAR
0.05
●
●
λ1 λ2 λ3 λ4 λ5 λ6 λ7 λ8 λ9
Figura 4.12: Sumário a posteriori de λs2 para os planejamentos (i) e (ii-a) com base nas
500 amostras simuladas.
4.6
Conclusões
Neste capı́tulo apresentou-se a principal contribuição deste trabalho, que foi a
proposta de um modelo desagregado que se ajuste a amostras adaptativas selecionadas
de populações raras e agrupadas. O modelo é construı́do no nı́vel das unidades da grade,
o que permitiu a inserção da suposição de heterogeneidade entre redes distintas. A
inferência Bayesiana para o modelo é feita usando o método RJMCMC, pois neste caso o
tamanho do espaço paramétrico é desconhecido. Portanto, o ajuste do modelo proposto
105
necessita de métodos mais custosos computacionalmente do que o modelo agregado, onde
apenas o MCMC é necessário.
No geral, o modelo apresentou uma boa performance nos estudos de simulação
realizados e ao ajustá-lo com a população real do marreco da asa azul, resultados
mais satisfatórios foram obtidos quando comparado com o modelo agregado. Por outro
lado, foi possı́vel observar que ao diminuir o grau de heterogeneidade da população
o desempenho do modelo agregado com relação a estimação de T , o qual é o maior
interesse neste trabalho, tende a melhorar e a tornar-se mais próximo ao obtido quando
ajustado o modelo de mistura. Portanto, recomenda-se o uso do modelo proposto quando
de fato a heterogeneidade é um comportamento presente nos dados, visto que o custo
computacional é maior neste caso.
Um sumário das conclusões mais relevantes extraı́das dos estudos de simulação
realizados neste capı́tulo é apresentado na Tabela 4.10.
Finalmente, com o propósito de melhorar a previsão e estimação do modelo de
mistura, foi apresentada uma aplicação do modelo de mistura ao plano amostral
adaptativo duplo.
Este planejamento tende a fornecer mais informações sobre a
população de pesquisa, com um custo operacional controlado. Nesta extensão verificou-se
que é possı́vel obter resultados eficientes ainda que com um número menor de observações
da variável de interesse e usando uma variável auxiliar indicadora de presença da
caracterı́stica de interesse.
106
Tabela 4.10: Resumo das principais conclusões acerca dos estudos simulados realizados
com o modelo de mistura proposto em (4.4).
Variando N , α e β
(1) Melhores resultados à medida que os valores de N , α e β aumentam.
(2) Maiores dificuldades de estimação de λs̄ que λs̄ .
Distribuição a priori de λ
(1) Distribuição a posteriori de R sensı́vel à escolha de τ .
(2) Escolha de τ não afeta a distribuição a posteriori de T .
(3) Os EQMR obtidos na previsão de T são menores quando assume-se distribuição
a priori dependente para λ.
Nı́vel de heterogeneidade
(1) Mesmo sob nı́veis mais intensos de homogeneidade bons resultados são
atingidos na previsão de T , mas surgem problemas na estimação de ν e β.
(2) Comparando com o modelo agregado, percebe-se que o modelo proposto é adequado
principalmente para populações heterogêneas. Sob maiores nı́veis de homogeneidade,
o desempenho dos modelos torna-se similar.
107
Capı́tulo 5
Conclusões e trabalhos futuros
Ao longo deste trabalho foram revisadas duas possı́veis formas de fazer previsão em
populações raras e agrupadas: a inferência baseada na aleatorização do plano amostral e
a abordagem baseada em modelos de superpopulação. No primeiro caso, apresentouse o planejamento amostral adaptativo por conglomerados e, no segundo, o modelo
proposto por Rapley e Welsh (2008), o qual é ajustado sob o enfoque Bayesiano. Estudos
simulados com base em populações artificiais e real foram apresentados e ambas as
abordagens foram comparadas principalmente em nı́veis de eficiência da previsão do
total populacional. Tendo em vista um bom desempenho do modelo de Rapley e Welsh
(2008), as metodologias propostas neste trabalho permanecem no contexto de inferência
em população finita baseada em modelos.
Realizar pesquisas em populações raras e agrupadas é uma tarefa árdua e necessita em
geral de metodologias especı́ficas que usem na sua formulação a estrutura da população.
No entanto, estas populações podem ser ainda mais problemáticas se apresentarem
uma dinâmica populacional, o que é uma caracterı́stica também comum neste contexto.
Buscando tratar situações como esta, foi apresentada uma extensão do modelo de Rapley
e Welsh (2008). Em particular, a extensão é voltada principalmente para populações em
crescimento ou decrescimento e final estabilização com a evolução do tempo.
Por outro lado, questões como a modelagem no nı́vel agregado das redes, suposições
de homogeneidade entre as redes e de relação direta entre a frequência esperada de
um fenômeno e o tamanho de uma rede no qual ele é observado, restringem o modelo
108
de Rapley e Welsh (2008) a algumas especı́ficas populações com estas caracterı́sticas.
Com o objetivo de tratar destas questões, foi proposto um modelo de mistura a nı́vel
desagregado que supõe heterogeneidade entre as redes, e consequentemente que o número
de ocorrências de um fenômeno em uma rede não depende necessariamente apenas do
tamanho desta. Como foi visto, para fazer inferência para este modelo fez-se necessário
técnicas mais sofisticadas, pois a dimensão do vetor paramétrico é também um parâmetro.
Em particular, foi utilizado o método de RJMCMC. O modelo mostrou-se mais eficiente
que o modelo agregado em casos de heterogeneidade. Por outro lado, à medida que o
nı́vel de heterogeneidade diminui a performance dos modelos torna-se semelhante.
Finalmente, a metodologia proposta foi aplicada ao plano amostral adaptativo duplo
por conglomerados, com o objetivo de adquirir mais informações que auxiliem a estimar
os parâmetros do modelo de mistura proposto em (4.4) associados às unidades que não
foram observadas. Em particular, a variável auxiliar utilizada nesta extensão caracterizase como uma indicadora da ausência ou presença da observação de interesse, ou seja, está
totalmente relacionada com a variável de pesquisa.
5.1
Trabalhos futuros
Na extensão apresentada na Seção 3.2 do Capı́tulo 3 supor uma amostra independente
a cada instante de tempo pode não ser viável em algumas situações práticas. No entanto,
como o modelo é formulado de forma agregada, isto traz dificuldades a incorporar outros
planejamentos mais viáveis. Com isso, há interesse em aplicar o modelo de mistura
proposto a planos amostrais que apresentem dependência temporal.
Com relação ao desenho amostral adaptativo duplo, seria interessante investigar
um tamanho de amostra ótimo na primeira e/ou na segunda fase, de modo a ser
eficiente e minimizar o custo operacional. Além disso, há interesse também em aplicar a
metodologia, supondo outras variáveis auxiliares relacionadas com a variável de interesse
que não somente indicadoras de presença da caracterı́stica de interesse.
Além disso, dentro de uma rede é comum que unidades tenham frequência de
observações que varia de acordo com a distância ao centroide da rede. Por exemplo,
109
espera-se que unidades dentro de uma rede tenham frequência de observações que varia de
acordo com a distância ao centroide da rede. O processo pontual conglomerado de Poisson
(ver Diggle et al. (1983)) é um exemplo de população com este comportamento. Dessa
forma, uma ideia futura para o modelo de mistura proposto é a inserção de componentes
espaciais na média da distribuição da variável resposta que dependam da distância. Um
importante aspecto a ser considerado nesta proposta futura é a definição do centroide,
visto que uma rede em geral não é regular. Além disso, a proposta seria incorporar esta
estrutura espacial na parte do modelo que se ajusta à amostra coletada, pois para a parte
não amostrada não há conhecimento da localização e nem das unidades que compõem as
redes, o que inviabilizaria a ideia nestas unidades.
5.1.1
Planejamento amostral ótimo
Como o desenho amostral adaptativo caracteriza-se pela seleção da amostra em fases,
seria razoável estudar a incorporação de um planejamento amostral ótimo, a fim de buscar
unidades amostrais que possam ser mais promissoras para a estimação do parâmetro
populacional de interesse.
Em desenhos amostrais convencionais a amostra completa é planejada de uma vez,
antes mesmo da seleção. Um exemplo de planejamento em duas fases é aquele em
que a amostra inicial de n1 unidades é selecionada e os valores de Y são observados
e, posteriormente, uma amostra adicional de n2 unidades é selecionada, cujo tamanho
depende dos valores observados na primeira amostra. A amostragem adaptativa seria
uma classe de desenhos com L fases, em que L é uma variável aleatória.
De forma geral, um planejamento ótimo é uma tarefa que costuma envolver
metodologias para obtenção de máximos e mı́nimos de funções objetivo. Estas funções
objetivo quantificam os ganhos e perdas associados às possı́veis decisões a serem tomadas.
A ideia de um planejamento ótimo com duas fases é descrito a seguir e pode ser visto
com maiores detalhes em Thompson e Seber (1996).
Suponha um desenho com tamanho amostral fixo em n unidades e suponha que
dessas unidades n1 foram selecionadas e observadas. Seja ys1 os valores de Y associados
a esta amostra inicial. A amostra restante a ser observada é s2 e de forma análoga
110
defina ys2 . Logo, a amostra completa é dada por s = (s1 , s2 ), com respectivos ys =
(ys1 , ys2 ). O objetivo é prever uma função populacional qualquer W = w(Y), como
o total populacional por exemplo, a partir de uma função da amostra H(d), tal que
d = (s, ys ). Deseja-se que H seja não viesado de acordo com o modelo. A função que
minimiza o erro quadrático médio de previsão E ((H − W )2 | s) é a esperança condicional
H(d) = E (W | d).
Finalmente, a questão é se a seleção das n2 unidades restantes deve depender dos
valores de ys1 . Neste caso, a ideia seria selecionar uma amostra s2 que minimize a função
objetivo:
gs2 (s1 , ys1 ) = E (h(s1 , ys1 , s2 , Ys2 ) − w(Y))2 | s1 , ys1
Z
= (h − w)2 [ys̄1 | s1 , ys1 ]dys̄1 .
Este mesmo argumento pode ser estendido para desenhos com múltiplas fases.
Portanto, um interesse futuro seria incorporar o planejamento amostral ótimo nos
modelos estudados neste trabalho. Por exemplo, ao selecionar uma amostra adaptativa
inicial e verificar alguns locais mais informativos na região, é possı́vel continuar o processo
de seleção da amostra propondo outros locais mais eficientes do que uma amostra aleatória
simples inicial. Ou ainda no plano adaptativo duplo, onde a segunda fase depende da
primeira. Nessa mesma proposta é possı́vel também avaliar um tamanho de amostra
ótimo.
111
Apêndice A
Resultados dos modelos ajustados
no Capı́tulo 3
A.1
Modelo (3.1)
Neste apêndice são apresentados os traços das cadeias dos parâmetros para o modelo
(3.1) e uma ilustração do sumário da distribuição a posteriori dos parâmetros obtido para
100 populações em 16 cenários gerados.
112
400
1000
3000
T
2000
T
1000
0
400
1000
iteração
T
400
500
7 9
0
1000
1000
0
iteração
1000
0
T
γ
400
iteração
8 12
0.4
0.1
β
0.3
1000
1000
12
400
iteração
0.0
400
400
iteração
γ
0
iteração
0
400
iteração
500
0
1000
0.1 0.3
β
400
0
14
γ
0
iteração
0.05 0.30
0
1000
7 10
β
0.1
1000
400
iteração
0.4
0.4
α
400
iteração
α
0
1000
iteração
0.1
0
α
400
0
γ
0
1500
1000
iteração
0 20
β
0.0
400
0.4
0.3
α
0.0
0
0
iteração
400
0
1000
400
1000
0
iteração
iteração
400
1000
iteração
Figura 1.1: Traços das cadeias dos parâmetros α, β, γ e total populacional T para um
dado artificial gerado fixando α = 0.05 e β ∈ {0.05, 0.1, 0.15, 0.2}, com respectivos valores
T
400
0
T
9
400
0
1000
400
1000
0
iteração
1000
T
γ
8
0.1
400
iteração
12
β
0.3
0
0
400
1000
0
400
iteração
iteração
1000
400
1000
iteração
0.4
iteração
400
iteração
500 1500
γ
β
0
1500
1000
iteração
0.05
1000
iteração
1000
0
0
1000
400
iteração
T
γ
400
iteração
0.25
0.4
α
400
0
6
0
0.1
0
1000
10
β
1000
400
iteração
0.1
400
iteração
0.0
0
1000
0.4
0.4
α
0.1
0
α
400
iteração
1000
500 2000
0
500 1500
γ
1000
iteração
12
400
8
0.05
0.1
0
10 12
β
α
0.4
0.25
verdadeiros em cinza.
0
400
1000
iteração
Figura 1.2: Traços das cadeias dos parâmetros α, β, γ e total populacional T para um
dado artificial gerado fixando α = 0.1 e β ∈ {0.05, 0.1, 0.15, 0.2}, com respectivos valores
verdadeiros em cinza.
113
400
2000
1000
0
T
400
1000
0
iteração
iteração
400
1000
iteração
11
0
1000
1000
T
400
iteração
γ
0
iteração
400
iteração
500
γ
0
1000
8
β
0.05
1000
0
8
400
0.30
0.4
α
400
1000
10 12
0.25
0.05
0
iteração
0.1
0
400
iteração
β
0.4
1000
iteração
1000
0
0
1000
400
iteração
T
γ
400
iteração
0.1
400
0
6
0
iteração
0
1000
12
0.5
1000
400
iteração
0.1
400
T
0
1000
β
0.3
α
0.0
0
α
400
iteração
1500
0
500 1500
1000
600 1600
γ
β
0.05
400
iteração
9.0 10.5
0.4
α
0.1
0
400
1000
iteração
Figura 1.3: Traços das cadeias dos parâmetros α, β, γ e total populacional T para um
dado artificial gerado fixando α = 0.15 e β ∈ {0.05, 0.1, 0.15, 0.2}, com respectivos valores
400
400
2000
0
400
1000
0
iteração
1000
T
γ
10
7
0.1
400
iteração
0
400
1000
0
400
iteração
iteração
1000
400
1000
iteração
0.4
β
α
0.1
0
400
iteração
2500
0
1000
iteração
0.4
1000
T
γ
0
1000
500
0.30
0.05
1000
400
iteração
β
0.4
400
iteração
400
iteração
T
0
1000
iteração
0.1
0
0
500
γ
0
1000
8.5 10.5
0.20
1000
iteração
400
iteração
0.05
400
T
0
1000
β
0.4
α
0.1
0
α
400
iteração
2000
0
1000
500 2000
1000
9 11
400
iteração
500
γ
0.05
0
8.0 10.0
0.20
β
0.4
0.1
α
verdadeiros em cinza.
0
400
1000
iteração
Figura 1.4: Traços das cadeias dos parâmetros α, β, γ e total populacional T para um
dado artificial gerado fixando α = 0.2 e β ∈ {0.05, 0.1, 0.15, 0.2}, com respectivos valores
verdadeiros em cinza.
114
●
●
●
●
●
EQM − α
●
●
●
0.050
●
●
●
●
●
●
0.020
●
●
●
●
●
●
●
●
●
●
●
●
●
●
α=0.10
●
●
●
●
●
●
●
α=0.05
●
0.8
●
●
●
●
α=0.20
●
0.6
EQM − γ
●
●
●
●
●
●
●
●
●
●
α=0.05
●
●
●
α=0.15
●
α=0.10
●
α=0.15
●
●
α=0.20
1.5
●
●
●
●
●
0.7
EQMR − T
●
●
●
1.0
●
α=0.10
●
●
●
α=0.20
●
0.5
0.9
●
●
α=0.15
●
●
●
0.4
0.7
0.5
α=0.10
α=0.20
●
●
●
●
α=0.05
α=0.15
●
●
α=0.20
●
●
●
●
●
0.02
0.9
●
●
1.0
●
●
α=0.15
●
●
●
●
α=0.10
●
●
0.08
●
●
EQM − β
●
α=0.05
0.5
●
α=0.05
0.7
0.9
●
●
●
α=0.20
0.5
Cobertura − β
●
α=0.15
0.06
α=0.10
●
Cobertura média − γ
●
0.04
α=0.05
●
●
●
●
●
0.0
0.3
Cobertura − T
●
0.035
1.0
0.8
0.6
●
●
●
0.4
Cobertura − α
●
●
α=0.05
α=0.10
(a)
α=0.15
α=0.20
α=0.05
●
●
α=0.10
●
●
●
α=0.15
●
●
●
●
●
α=0.20
(b)
Figura 1.5: Sumário da distribuição a posteriori dos parâmetros α, β, γ e T para 100
populações em 16 cenários com amostra inicial de 5%N e 10%N . Em (a) os triângulos
representam as probabilidades de cobertura dos intervalos HPD de 95% para a amostra
de 5%, os cı́rculos cheios para a amostra de 10% e a linha tracejada em vermelho o nı́vel
nominal de 95%. Em (b) estão o EQM para cada parâmetro e o EQMR para T .
115
A.2
Modelo de crescimento (3.4)
Nas Figuras 1.6 e 1.7 estão os resultados do ajuste do modelo (3.4) para duas das
populações artificiais geradas. A primeira é para uma população em crescimento ao longo
1000
1000
1000
400
0
1000
400
1000
iteração
T49
800
(e) γ
T37
800
0
10.2
9.8
400
(d) β
400
400
0
iteração
T25
800
700
0
γ
β
0.09
1000
(c) c
T13
(f) T1
1000
400
iteração
300
400
iteração
0
(b) b
T1
300
100
400
iteração
(a) a
0
0.12
−0.10
c
−1.3
0
0
400
1000
400
1000
400
400
iteração
−0.20
0
−1.7
b
−1.6
−1.9
a
do tempo e a segunda para uma que decresce.
0
400
iteração
iteração
iteração
iteração
(g) T13
(h) T25
(i) T37
(j) T49
1000
Figura 1.6: Sumário da distribuição a posteriori de Θ e do total populacional para uma
população em crescimento ao longo do tempo. Em (a)-(e) estão os traços das cadeias
da distribuição a posteriori dos parâmetros a, b, c, β e γ. De (f )-(j) estão os traços das
cadeias para os totais em alguns tempos. A linha em cinza representa o valor verdadeiro
usado na geração dos dados artificiais.
116
1000
400
10.2
9.6
0.08
γ
β
1000
400
1000
400
1000
400
1000
iteração
(e) γ
T37
450
0
0
(d) β
T25
500
0
0
iteração
200
200
1000
(c) c
T13
500
900
T1
(f) T1
1000
400
iteração
(b) b
500
400
iteração
0
T49
500
400
iteração
(a) a
0
0.13
−0.10
0
0
400
1000
200
1000
200
400
iteração
−0.20
c
1.0
0.7
b
−2.1
a
−2.4
0
0
400
iteração
iteração
iteração
iteração
(g) T13
(h) T25
(i) T37
(j) T49
1000
Figura 1.7: Sumário da distribuição a posteriori de Θ e do total populacional para uma
população em decrescimento ao longo do tempo. Em (a)-(e) estão os traços das cadeias
da distribuição a posteriori dos parâmetros a, b, c, β e γ. De (f )-(j) estão os traços das
cadeias para os totais em alguns tempos. A linha em cinza representa o valor verdadeiro
usado na geração dos dados artificiais.
117
Apêndice B
Cálculos envolvidos na inferência
para o modelo proposto
Neste apêndice são apresentadas expressões importantes envolvidas no algoritmo
RJMCMC, utilizado para inferência a posteriori para o modelo de mistura proposto (4.4).
Primeiramente estão as distribuições condicionais completas para o vetor paramétrico
Θ = (Xs̄ , Rs̄ , s̄ , Cs̄ , Ys̄ , α, β, λ, ν). Dessa forma, a variável resposta Ys̄ , por exemplo,
é também considerada um parâmetro e portanto é estimada da mesma maneira que
as demais quantidades. Além disso, será apresentada a probabilidade de aceitação do
algoritmo RJMCMC, passando por alguns cálculos importantes.
B.1
Distribuições condicionais completas
Para as distribuições condicionais completas que apresentam forma analı́tica
conhecida, o Amostrador de Gibbs pode ser utilizado. Para as que não apresentam
forma fechada um método indireto de amostragem é necessário, em particular, passos
de Metropolis-Hastings podem caracterizar a obtenção dessas amostras a posteriori. As
118
distribuições apresentadas a seguir são obtidas ao assumir-se as seguintes distribuições a
priori independentes para o modelo (4.4):
λj ∼ Gama(d, ν), j = 1, . . . , R,
ν ∼ Gama(e, f ),
α ∼ Beta(aα , bα ),
β ∼ Beta(aβ , bβ ).
A seguir estão as distribuições condicionais completas.
• De α:
[α | ·] ∝
αXs +Xs̄ (1 − α)N −Xs −Xs̄ aα −1
α
(1 − α)bα −1 .
N
1 − (1 − α)
Para gerar amostras desta distribuição deve-se utilizar passos de Metropolis-Hastings,
visto que esta não apresenta forma analı́tica conhecida.
• De β:
[β | ·] ∝
β Rs +Rs̄ (1 − β)Xs +Xs̄ −Rs −Rs̄ aβ −1
β
(1 − β)bβ −1 .
1 − (1 − β)Xs +Xs̄
Como [β | ·] também não possui forma analı́tica fechada, deve-se utilizar passos de
Metropolis-Hastings, para amostrar desta distribuição de probabilidade.
• De λ: Para j = 1, . . . , Rs + Rs̄ ,
P
[λj | ·] ∝
λj
{i:i =j}
Yi +d−1
exp{−λj (ν + Cj )}
.
1 − exp(−λj )
(2.1)
Observe que [λj | ·] não possui forma fechada conhecida. Para gerar amostras de sua
distribuição a posteriori é necessário utilizar um passo de Metropolis-Hastings.
• De s̄ : Para i, j ∈ s̄,
119
λYj i exp(−λj )
Cj
[i = j | ·] ∝
.
Xs + Xs̄ Yi ![1 − exp(−λj )]
Neste caso, i é amostrado diretamente dos possı́veis valores, com a probabilidade
acima. Note que o modelo proposto é aplicável a populações divididas em redes nãovazias, logo toda rede deve ter pelo menos uma observação. Portanto, na condicional
completa de i ainda é incluı́do uma indicadora de que todas as Rs̄ redes tenham pelo
menos uma unidade alocada.
• De (Xs̄ , Cs̄ ):
m
Y
αXs̄ (1 − α)−Xs̄
Zi × gil ,l
(1 − β)Xs̄
Pl−1
PN −X+R l
Zi − k=0 Zik (N − Xs − Xs̄ )! (1 − (1 − β)Xs +Xs̄ )
i=1
l=1
Y C Y
1
R−(Xs̄ −Rs̄ )
× (Xs + Xs̄ )−(Xs +Xs̄ )
Cj j
(Cj − 1)!
[Xs̄ , Cs̄ | ·] ∝
{j:j∈s̄}
×
Y
{j:j∈s̄}
{j:j∈s̄}
exp{−λj Cj }
.
[1 − exp(−λj )]Cj
A amostragem de (Xs̄ , Cs̄ ) é feita de forma conjunta, e como a distribuição condicional
completa não tem forma analı́tica fechada, o algoritmo de Metropolis-Hastings é utilizado.
A proposta de Xs̄ é baseada num passeio aleatório em torno do valor corrente de Xs̄ e a
proposta de Cs̄ baseia-se na Multinomial(Xs̄ − Rs̄ , R1s̄ 1Rs̄ ).
• De Ys̄ :
P
Y λj
Q
[Ys̄ | ·] ∝
{j:j∈s̄}
{i:i =j}
{i:i =j}
Yi
Yi !
.
(2.2)
Portanto, Ys̄i ∼ Poisson truncada(λs̄i ), j ∈ s̄. Logo, para amostrar desta distribuição
podemos utilizar o Amostrador de Gibbs.
• De ν:
[ν | ·] ∝ ν
(Rs +Rs̄ )d+e−1
exp{−ν(f +
RX
s +Rs̄
λj )}.
j=1
Logo, ν ∼ Gamma ((Rs + Rs̄ )d + e, f +
PRs +Rs̄
j=1
podemos utilizar o Amostrador de Gibbs.
120
λj ) e para amostrar desta distribuição
B.2
Probabilidade
de
aceitação
do
algoritmo
RJMCMC
Se é proposto um movimento de “divisão”, ou seja que leva de (Cj ∗ , j ∗ , λj ∗ ) a
(Cj1 , Cj2 , j1 , j2 , λj1 , λj2 ), para j ∗ , j1 e j2 pertencentes a s̄, o movimento é aceito com
probabilidade dada por min(1, A), tal que A é dada por (4.3), e para este modelo tem a
seguinte forma:
P
A=
exp{−(Cj1 λj1 + Cj2 λj2 )}λj1
{i:i =j1 }
Yi
P
λj 2
P
exp{−Cj ∗ λj ∗ }λj ∗
{i:i =j2 }
{i:i =j ∗ }
Yi
Yi
(1 − exp(−λj1 ))−Cj1 (1 − exp(−λj2 ))−Cj2
(1 − exp(−λj ∗ ))−Cj∗
C
C
Cj1j1 Cj2j2
[{ij1 , ij2 }] p(Rs̄ + 1)
(Cj ∗ − 1)!
−(Cj1 +Cj2 −Cj ∗ )
×
×
(Rs + Rs̄ )
×
× (Rs̄ + 1)
C ∗
[{ij ∗ }]
p(Rs̄ )
(Cj1 − 1)!(Cj2 − 1)!
Cj ∗j
d−1
νd
λj1 λj2
exp{−ν(λj1 + λj2 − λj ∗ )}
×
Γ(d)
λj ∗
pk|k+1
Cj ∗
×ρ
,
×
pk+1|k Palloc q(u1 )q(u2 )
Xs + Xs̄
onde a primeira linha consiste da razão das verossimilhança avaliada nestes pontos, a
segunda e terceira linha apresentam a distribuição a priori dos parâmetros. No final da
segunda linha, o termo Rs̄ +1 vem da razão (Rs̄ +1)!/Rs̄ !, devido a ordem dos parâmetros.
A última linha apresenta a razão das probabilidades de transição entre os espaços, onde
Palloc é a probabilidade desta particular alocação ser feita, e o último termo é o jacobiano
da transformação.
121
Referências Bibliográficas
Besag, J. (1974) Spatial interaction and the statistical analysis of lattice systems. Journal
of the Royal Statistical Society. Series B (Methodological), 36, 192–236.
Bolfarine, H. e Zacks, S. (1992) Prediction theory for finite populations. Springer-Verlag
New York:.
Brown, J. A. e Manly, B. J. F. (1998) Restricted adaptive cluster sampling.
Environmental and Ecological Statistics, 5, 49–63.
Cassel, C.-M., Särndal, C.-E. e Wretman, J. H. (1977) Foundations of inference in survey
sampling. Wiley New York.
Clayton, D. e Bernardinelli, L. (1992) Bayesian methods for mapping disease risk.
Geographical and environmental epidemiology: methods for small area studies, 205–
220.
Conners, M. e Schwager, S. (2002) The use of adaptive cluster sampling for hydroacoustic
surveys. ICES Journal of Marine Science: Journal du Conseil, 59, 1314–1325.
Danaher, P. e King, M. (1994) Estimating rare household characteristics using adaptive
sampling. NZ Stat, 29, 14–23.
Diggle, P. J. et al. (1983) Statistical analysis of spatial point patterns. Academic Press.
Felix-Medina, M. H. e Thompson, S. K. (2004) Adaptive cluster double sampling.
Biometrika, 91, 877.
122
Gelman, A. (2006) Prior distributions for variance parameters in hierarchical models
(comment on article by browne and draper). Bayesian analysis, 1, 515–534.
Gelman, A., Carlin, J. B., Stern, H. S. e Rubin, D. B. (1995) Bayesian data analysis.
Chapman & Hall.
Geweke, J. (1992) Evaluating the accuracy of sampling-based approaches to the
calculations of posterior moments. Em Bayesian Statistics (eds. A. D. J. Bernardo,
J. Berger e A. Smith). Oxford University Press, New York.
Gilks, W. R. e Wild, P. (1992) Adaptive rejection sampling for Gibbs sampling. Applied
Statistics, 337–348.
Green, P. (1995) Reversible jump markov chain monte carlo computation and bayesian
model determination. Biometrika, 82, 711–732.
Horvitz, D. e Thompson, D. (1952) A generalization of sampling without replacement
from a finite universe. Journal of the American Statistical Association, 47, 663–685.
Kalton, G. (1991) Sampling flows of mobile human populations. Survey Methodology,
17, 183–194.
— (2001) Practical methods for sampling rare and mobile populations. Em Proceedings
of the Annual Meeting of the American Statistical Association, 5–9.
Kalton, G. e Anderson, D. (1986) Sampling rare populations. Journal of the Royal
Statistical Society. Series A (General), 149, 65–82.
Lambert, D. (1992) Zero-inflated poisson regression, with an application to defects in
manufacturing. Technometrics, 34, 1–14.
Marin, J.-M., Mengersen, K. e Robert, C. P. (2005) Bayesian modelling and inference on
mixtures of distributions. Handbook of statistics, 25, 459–507.
McDonald, L. L. (2004) Sampling rare populations.
Em Sampling rare or elusive
species: concepts, designs, and techniques for estimating population parameters (ed.
W. Thompson), cap. 4, 11–42. Island Press Washington, DC, USA.
123
Migon, H. e Gamerman, D. (2006) Generalized exponential growth models a bayesian
approach. Journal of Forecasting, 12, 573–584.
Neyman, J. e Scott, E. (1958) Statistical approach to problems of cosmology. Journal of
the Royal Statistical Society. Series B (Methodological), 20, 1–43.
R Core Team (2013) R: A Language and Environment for Statistical Computing. R
Foundation for Statistical Computing, Vienna, Austria. URLhttp://www.R-project.
org.
Raftery, A. E. e Lewis, S. M. (1992) One long run with diagnostics: Implementation
strategies for markov chain monte carlo. Statistical Science, 7, 493–497.
Rapley, V. (2004) Model-Based Adaptive Cluster Sampling. Tese de Doutorado, University
of Southampton.
Rapley, V. e Welsh, A. (2008) Model-based inferences from adaptive cluster sampling.
Bayesian Analysis, 3, 717–736.
Richardson, S. e Green, P. (1997) On bayesian analysis of mixtures with an unknown
number of components. Journal of the Royal Statistical Society, Series B, 59, 731–792.
Roeder, K. e Wasserman, L. (1997) Practical bayesian density estimation using mixtures
of normals. Journal of the American Statistical Association, 92, 894–902.
Roesch, F. (1993) Adaptive cluster sampling for forest inventories. Forest Science, 39,
655–669.
Salehi, M. M. e Seber, G. (1997) Two-stage adaptive cluster sampling. Biometrics, 53,
959–970.
Skinner, C., Holt, D. e Smith, T. (1989) Analysis of complex surveys. John Wiley &
Sons.
Smith, D., Brown, J. e Lo, N. (2004) Application of adaptive sampling to biological
populations. Em Sampling rare or elusive species: concepts, designs, and techniques for
124
estimating population parameters. Island, Washington, DC, USA (ed. W. Thompson),
cap. 5, 77–122. Island Press Washington, DC, USA.
Smith, D., Conroy, M. e Brakhage, D. (1995) Efficiency of adaptive cluster sampling for
estimating density of wintering waterfowl. Biometrics, 51, 777–788.
Spiegelhalter, D. J., Best, N. G., Carlin, B. P. e Van Der Linde, A. (2002) Bayesian
measures of model complexity and fit. Journal of the Royal Statistical Society: Series
B (Statistical Methodology), 64, 583–639.
Sudman, S. e Kalton, G. (1986) New developments in the sampling of special populations.
Annual Review of Sociology, 12, 401–429.
Tanner, M. A. (1993) Tools for Statistical Inference: Methods for the Exploration of
Posterior Distributions and Likelhood Functions. Springer-Verlag.
Thompson, S. e Collins, L. (2002) Adaptive sampling in research on risk-related
behaviors. Drug and Alcohol Dependence, 68, 57–67.
Thompson, S. K. (1990) Adaptive cluster sampling. Journal of the American Statistical
Association, 85, 1050–1059.
— (1991) Stratified adaptive cluster sampling. Biometrika, 78, 389–397.
Thompson, S. K. e Seber, G. A. F. (1996) Adaptive sampling. Wiley New York.
Turk, P. e Borkowski, J. (2005) A review of adaptive cluster sampling: 1990–2003.
Environmental and Ecological Statistics, 12, 55–94.
Viallefont, V., Richardson, S. e Green, P. J. (2002) Bayesian analysis of Poisson mixtures.
Journal of Nonparametric Statistics, 14, 181–202.
125