MODELOS DE REGRESSAO COM RESPOSTAS~ PARCIAIS Heleno Bolfarine Jorge Bazan - SBM · 2014. 10. 20. · esta ultima sendo conhecida como raz~ao de Mills. Os resultados para modelos de

MODELOS DE REGRESSÃO COM RESPOSTASPARCIAIS

Heleno Bolfarine

Jorge Bazan

MODELOS DE REGRESSÃO COM RESPOSTASLIMITADAS E CENSURADAS

H. Bolfarine

Universidade de São PauloInstituto de Matemática e Estat́ıstica

Departamento de Estat́ıstica

J. Bazan

Universidade de São PauloInstituto de Ciências Matemáticas e de ComputaçãoDepartamento de Matemática Aplicada e Estat́ıstica

Prefácio

Neste trabalho, desenvolvemos análises Bayesiana e clássica para modelosde regressão com respostas limitadas ou censuradas.São consideradas ex-tensões do modelo tobit usual normalmente distribúıdo em duas direções. Aprimeira considera modelos mais gerais que o modelo normal proporcionadapelo modelo potência-normal, o qual pode ajustar dados com certo grau deassimetria e bimodalidade. Uma outra direção em que estendemos o modeloestá voltada para situações onde temos excesso (inflação) de zeros. No casoem que as observações são proporções (no intervalo (0, 1)), podemos ter da-dos com excesso de zeros e uns. Discute-se especificacação de prioris poucoinformativas e algoritimos tipo MCMC para estimação dos parâmetros domodelo. Procedimentos de estimação alternativos são desenvolvidos usandoo método de máxima verossimilhança. Aplicações a vários conjunto de da-dos são apresentadas. Um conjunto de dados, em especial, é o conjuntode dados sobre a resposta sorológica em um programa de vacinação contrasarampo no Haiti. Além disso, são estudadas aplicações a outros conjuntosde dados relacionados com os modelos considerados.

Este manuscrito, direcionado a extensões do modelo tobit, está organi-zado da seguinte forma: o Caṕıtulo 1 enfoca resultados básicos de modelospara dados censurados e truncados. No Caṕıtulo 2 apresentamos uma breverevisão do modelo tobit com sugestões de extensões que podem ser consid-eradas substituindo-se a distribuição normal por modelos mais robustos eflexiveis como os modelos potência-normal (Pewsey et al., 2012) e t-Student.Aplicações a dados reais mostram bom desempenho dos modelos propostos.O Caṕıtulo 3 está dedicado ao modelo tobit com excesso de zeros em queduas extensões são consideradas. Análise de dados reais são apresentadasilustrando o bom desempenho dos modelos estudados. O Caṕıtulo 4 discutemodelos α-potência para dados duplamente censurados com ênfase nos casos(0, 1), com posśıveis excessos de zeros e uns. O Caṕıtulo 5 estuda modelosbimodais censurados. Este texto está direcionado a alunos do último anodo bacharelado e ińıcio do mestrado em Estat́ıstica.

v

vi

Heleno Bolfarine [email protected] Bazan [email protected]

São Carlos, SP, janeiro de 2013

Sumário

1 Dados limitados 1

1.1 Truncamento . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

1.2 Censura . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

1.3 Resultados básicos sobre truncamento e censura . . . . . . . . 2

1.3.1 Distribuição normal truncada . . . . . . . . . . . . . . 2

1.3.2 Distribuição normal censurada . . . . . . . . . . . . . 3

1.4 Alguns conjuntos de dados . . . . . . . . . . . . . . . . . . . . 3

1.4.1 Vacinação no Haiti . . . . . . . . . . . . . . . . . . . . 4

1.4.2 Horas trabalhadas por ”donas”de casas . . . . . . . . 4

2 O modelo tobit 7

2.1 O modelo tobit normal . . . . . . . . . . . . . . . . . . . . . . 7

2.2 Extensões robustas do modelo tobit . . . . . . . . . . . . . . 11

2.3 Aplicações . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

2.3.1 Aplicação do modelo tobit-normal . . . . . . . . . . . 13

2.4 Aplicação do modelo tobit potência-normal . . . . . . . . . . 13

3 O modelo tobit com excesso de zeros 15

3.1 Modelos com excesso de zeros . . . . . . . . . . . . . . . . . . 15

3.2 A distribuição log-α-potência . . . . . . . . . . . . . . . . . . 16

3.3 O modelo bernoulli/log-α-potência . . . . . . . . . . . . . . . 18

3.4 Aplicação: dados do Haiti . . . . . . . . . . . . . . . . . . . . 22

3.5 Aplicação: dados de Mroz . . . . . . . . . . . . . . . . . . . . 25

4 Modelo α-potência inflacionado de zeros e/ou uns 27

4.1 Modelos duplamente censurados . . . . . . . . . . . . . . . . 27

4.2 Distribuições PN para dados censurados . . . . . . . . . . . . 29

4.3 Modelo potência-normal duplamente censurado . . . . . . . . 29

4.4 A transformação logaŕıtmica . . . . . . . . . . . . . . . . . . 32

vii

viii SUMÁRIO

4.5 O modelo Bernoulli duplamente censurado com mistura potência-normal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

4.6 Estimação . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 334.7 Matriz de informação observada . . . . . . . . . . . . . . . . . 344.8 Modelos censurados para inflação de zeros e uns . . . . . . . 374.9 Mistura Bernoulli/LPN . . . . . . . . . . . . . . . . . . . . . 394.10 Ilustração com dados reais . . . . . . . . . . . . . . . . . . . . 394.11 Testando modelos disjuntos . . . . . . . . . . . . . . . . . . . 404.12 Conclusões . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

5 Modelos bimodais censurados 435.1 Modelos assimétricos bimodais . . . . . . . . . . . . . . . . . 445.2 Extensões bimodais para modelos simétricos . . . . . . . . . . 44

5.2.1 Aplicação: Dados de poluição. . . . . . . . . . . . . . 465.3 Modelo flex́ıvel normal censurado . . . . . . . . . . . . . . . . 47

5.3.1 Momentos . . . . . . . . . . . . . . . . . . . . . . . . . 485.3.2 Extensão para localização-escala . . . . . . . . . . . . 485.3.3 Estimação . . . . . . . . . . . . . . . . . . . . . . . . . 495.3.4 Matriz de informação . . . . . . . . . . . . . . . . . . 50

5.4 O modelo bimodal simétrico normal censurado . . . . . . . . 535.4.1 Estimação por máxima verossimilhança . . . . . . . . 545.4.2 Matriz de informação esperada . . . . . . . . . . . . . 55

5.5 Modelo bimodal normal-assimétrico . . . . . . . . . . . . . . 555.5.1 A função log-verossimilhança . . . . . . . . . . . . . . 56

5.6 Analizando um conjunto de dados reais. Concentração de HIV. 575.7 Discussão final . . . . . . . . . . . . . . . . . . . . . . . . . . 59

Referências bibliográficas 61

Caṕıtulo 1

Dados limitados

A principal causa da ocorrência de dados incompletos é devido a (i) trunca-mento e (ii) censura.

1.1 Truncamento

Truncamento ocorre quando algumas observações tanto na variável respostacomo indepedentes (covariáveis, regressores) não estão dispońıveis. Por ex-emplo, a variável resposta (dependente) pode ser renda e somente pessoascom baixa (propriamente definida) renda são inclúıdadas na pesquisa. Por-tanto, truncamento ocorre quando a amostra é esolhida somente em parteda população.

1.2 Censura

Censura ocorre quando dados sobre a variável dependente não estão dispońıveispara algumas unidades da amostra. Mas para estas unidades, os dadospara a variáveis independentes (regressores) estão dispońıveis. por exem-plo, pessoas de todos os ńıveis de renda são incluidas na amostra mas, poralguma razão, pessoas com alto ńıvel de renda tem a mesma codificadaem R$100.000. Censura pode ser visto como um defeito na amostra - nãohavendo censura, amostra seria representativa. Truncamento em geral pro-duz maior perda de informação.

1

2 DADOS LIMITADOS 1.3

1.3 Resultados básicos sobre truncamento e cen-sura

É comum considerar que a variável resposta (Y ) é normalmente distribúıdacom média µ e variância σ2, que denotamos por N(µ, σ2), de tal forma que

E[Y ] = µ e V ar[Y ] = σ2.

O caso particular em que µ = 0 e σ = 1, ou seja, Z ∼ N(0, 1), temos afunção de densidade de probabilidade (fdp)

f(z) = φ(z) =1√2π

e−z2/2, z ∈ R.

A fdp de Y ∼ N(µ, σ2) segue da tansformação Y = µ+ σZ.A função de distribuição acumulada pode ser escrita como

Φ(y) = P [Y ≤ y] = Φ((y − µ)/σ),

de modo que

P [Y ≥ y] = 1− Φ((y − µ)/σ).

1.3.1 Distribuição normal truncada

Para truncamento pela esquerda (s.p.g.), com ponto de trunamento ”c”,temos

f(y|y > c) = f(y)1− F (c)

,

de modo que para Y ∼ N(µ, σ),

f(y|y > c) =1σφ(

y−µσ )

1− Φ( c−µσ ).

Assim, a função de verossimilhança para uma amostra de tamanho n dadistribuição normal truncada pode ser escrita como

L(µ, σ) =n∏i=1

1σφ(

yi−µσ )

1− Φ( c−µσ ).

1.4 ALGUNS CONJUNTOS DE DADOS 3

Temos também

E[y|y > c] = µ+ σλ(αc)

e

V ar[y|y > c] = σ2[1− δ(αc)],

onde αc = (c− µ)/σ,

δ(αc) = λ(αc)[λ(αc)− αc] e λ(αc) =φ(αc)

1− Φ(αc),

esta última sendo conhecida como razão de Mills.

Os resultados para modelos de regressão seguem dos resultados acimasubstituindo µ = x′β.

1.3.2 Distribuição normal censurada

Quando a distribuição é censurada à esquerda no ponto ”c”, observaçõescom valores menores ou iguais a c são substituidas por c ou seja,

y =

{y∗i , se y

∗i > c

c, se y∗i ≤ c.

Se uma variável continua Y com fdp f(.), e c é uma constante, entãopara variáveis censuradas à esquerda

f(y) = [f(yi)]Ii [F (c)]1−Ii ,

onde

Ii =

{1, se y∗i > c0, se y∗i ≤ c,

i = 1, . . . , n. O caso particular do modelo normal censurado segue tomandof = φ.

1.4 Alguns conjuntos de dados

Alguns conjuntos de dados muito utilizados na literatura são descritos aseguir.


1.4.1 Vacinação no Haiti

Dados contém informações sobre concentração de anticorpos em um grupode 330 crianças de até um ano no Haiti após serem vacinadas contra osarampo. As medições das concentrações são feitas por laboratórios comlimite de deteção mı́nimo (LDM) de 0.1 mm/l (ou -2.16 na escala logarit-mica). Isto significa que valores de concentrações iguais ou menores que 0.1são reportadas como sendo 0.1. Temos informação sobre a concentração (Y- variável resposta) , tipo de vacina (X1: Edmonton-Zagreb (1) e Schwarz(0)), dose (X2: alta (1) e médio (1)) e sexo (X3: masculino (O) e feminino(1)). O total de criana̧s no (ou abaixo do) limite de deteção é de 86. Umresumo dos dados é apresentado na tabela abaixo. Este conjunto de dadosesta disponibilizado em Moulton and Halsey (1995). Da Tabela 1.1. temosentão que a primeira criança tem concentração 0.1, tomou a vacina tipo 0(Schwarz) com dose média e é do sexo masculino.

Tabela 1.1: Dados sobre vacinação no HaitiCriança Concentração (Y ) Tipo (X1) Dose (X2) Sexo (X3)

1 0.1 0 0 02 0.1. 0 0 03 0.1 0 0 0... ... ... ... ...

316 15.475 1 0 0

Moulton and Halsey (1995) consideram uma distribuição log-normal paraobservações acima do LDM, e modelam o excesso de zeros com um modelologito, extendendo a proposta de Cragg (1971). Como se depreende dohistograma, a concentração de observações no LDM é bastante alta.

1.4.2 Horas trabalhadas por ”donas”de casas

Este conjunto de dados (Mroz, 1987) foi tomado do estudo da dinâmica derenda de 1975 com 753 observações das quais 428 correspondem a mulheres(casadas) com Y horas trabalhadas (não nulas) e as 325 remanescentes, cor-respondem a mulheres que não trabalharam (Y=0). O conjunto de dadoscompreende um total de 19 variáveis das quais consideramos

1.4 ALGUNS CONJUNTOS DE DADOS 5

1. LPF: variável ”dummy”= 1 se esposa trabalhou em 1975; =0, casocontrário;

2. WHRS: horas trabalhadas pela esposa em 1975;3. KL6: Número de criança com crianças menores que 6 anos no domićılio;4. K618: Número de crianças com idade entre 6 e 18 anos no domićılio;5. WA: Idade da esposa;6. WE: Escolaridade da esposa, em anos;7. WW: Salário da esposa em 1975.

Tabela 1.2: Dados sobre horas trabalhadashoras kids5 kids618 age educ nwifeinc exper

1610 -10.5 -9.2 0.8 -1.4 -0.3 3.41656 -12.5 -0.6 -0.2 0.6 -0.3 -5.61980 -7.5 -8.1 0.8 1.6 -0.3 4.4456 -8.5 -13.3 -0.2 1.6 -0.3 -4.61568 -11.5 0.0 0.8 0.6 1.7 -3.62032 11.5 -10.3 -0.2 -1.4 -0.3 22.41440 -5.5 -11.0 -0.2 0.6 3.7 0.41020 11.5 -9.2 -0.2 -1.4 -0.3 24.4

... ... ... ... ... ... ...0 0 3 39 9 28.3 12

Caṕıtulo 2

O modelo tobit

Neste caṕıtulo discutimos alguns resultados básicos sobre o modelo tobit.Apresentamos inicialmente o modelo tobit normal, a função de verossimi-lhança e as equações de estimação correspondentes. Mencionamos tambémextensões robustas com a substituição da distribuição normal pelos modelost-Student e potência-normal.

2.1 O modelo tobit normal

Pesquisadores são frequentemente confrontados com dados para os quais avariável resposta tem um limite inferior (que pode ser considerado comozero, sem perda de generalidade) e toma este valor para uma parte con-siderável das unidades amostrais. Este é o caso, por exemplo, dos dadossobre horas trabalhadas por donas de casa (Mroz, 1987).

Uma outra maneira é tratar os zeros como observações latentes (não ob-servadas) cont́ınuas. Esta idéia é popularizada em Tobin (1956) e o modeloresultante é chamado modelo tobit.

Formalmente, dada a variável de interesse Y , o modelo tobit pode serformulado como

yi =

{0, se wi ≤ 0wi, se wi > 0,

7

8 O MODELO TOBIT 2.1

onde a variável latente é wi = x′iβ + �i, com �i ∼ N(0, σ2), i = 1, . . . , n.

Consequentemente, denotamos as respostas observadas por yi, o valor das kvariáveis explanatórias para a i-ésima observação por xi ∈ Rk, os parâmetrosde regressão por β = (β0, β1, . . . , βk)

′ e o i-ésimo termo residual por �i.

Pode-se escrever o modelo acima como

yoi = Iiyi, yi = x′iβ + �i,

onde Ii = I(yi > 0), com xi = (xi1, . . . , xik)′, i = 1, . . . , n.

Com as suposições temos

yiind.∼ N(x′iβ, σ2),

i = 1, . . . , n. Note que, sendo yi ∼ N(µi, σ2), temos que

P [y0i = 0] = P [yi ≤ 0] = 1− Φ(µi/σ).

Por outro lado, sendo

yoi > 0, temos yoid= yi,

de modo que a função de verossimilhança pode ser escrita como

LN (β, σ2) =

n∏i=1

[1− Φ( 1σ

x′iβ)]1−Ii [(

1

σφ(

1

σ(yi − x′iβ)]Ii ,

com φ e Φ sendo a fdp e a fda da N(0,1).

Derivando a log-verossimilhança, temos as equações de verossimilhança

σ2 =1

n1y′D(y −Xβ),

X′(In −D)η = X′D(y −Xβ),

onde n1 =∑n

i=1 Ii, D = diag(I1, . . . , In),

η = (σr(−x1β′/σ), . . . , σr(−x′nβ/σ))′, r(z) =φ(z)

Φ(z).

2.1 O MODELO TOBIT NORMAL 9

Como as equações acima são não lineares, métodos numéricos são necessáriospara a sua solução. Uma alternativa seria a maximização direta da funçãolog-verossimilhança, o que pode ser implementado no aplicativo R.

A partir da derivada da função escore (avaliada no estimador de máximaverossimilhança (EMV)), podemos obter a matriz de informação observada(MIO). Invertendo a MIO, temos estimativas consistentes para a matriz decovariâncias assintóticas dos parâmetros do modelo.

Considerando a reparametrização

γ = β/σ, τ = 1/σ,

pode-se mostrar que as derivadas segundas podem ser escritas como

∂2 logLN∂γ∂γ′

= −n∑i=1

(1− Ii)r(−ci)(r(−ci)− ci)xix′i −n∑i=1

Iixix′i,

∂2 logLN∂γ∂τ

=n∑i=1

Iixix′i,

∂2 logLN∂γ∂γ′

=n1τ2−

n∑i=1

Iiy2i ,

onde ci = x′iγ.

A matriz de informação de Fisher (MIF) pode ser calculada a partir dasderivadas segundas acima usando os seguintes resultados (Arellano-Valle etal., 2012):

E[Ii] = P [Yi > 0] = Φ(ci),

E[IiYi] = E[Ii]E[Yi|Yi > 0] = (1/τ)(ciΦ(ci) + φ(ci)),

E[IiY2i ] =

1

τ2[(1 + ci)Φ(ci) + c

2iφ(ci)].

Para implementar o enfoque Bayesiano para o modelo tobit normal,pode-se usar o programa OpenBugs de duas maneiras diferentes. Uma dasmaneiras é entrar diretamente no OpenBugs usando

dummy[i] ∼ loglik(logLike[i]),


onde logLike[i] é o logaritimo da função de verossimilhança.

Uma maneira alternativa é entrar com o modelo de regressão normalcensurado, isto é

yi|β, σ2 ∼ NT (x′iβ, σ2, A),

com A = [0,∞).

Em geral,

Y ∼ TN(x′iβ, σ2, A)

se

fTN (y|µ, σ2, A) = c−1fN (y|µ, σ2)I(y ∈ A),

com

c =

∫AfN (y|µ, σ2)dx.

A função de log-verossimilhança para o modelo tobit (para T=c) para asituação onde o erro �i segue uma função de distribuição F , pode ser escritacomo

`(θ;Y ) =∑i

(1− Ii) ln[F (c− µσ

)] +∑i

Ii{− ln(σ) + ln(f(yi − µσ

)}

onde f = F ′, e

Ii =

{1, se y∗i > c0, se y∗i ≤ c,

A distribuição comumente usada com o modelo acima é a distribuiçãonormal, isto é, X ∼ N(µ, σ2),

F ′(x) = f(x) =1√

2πσ2e−

(x−µ)2

2σ2 .

2.2 EXTENSÕES ROBUSTAS DO MODELO TOBIT 11

2.2 Extensões robustas do modelo tobit

Uma distribuição que pode ser empregada no lugar da distribuição normalé a distribuição t-Student com fdp

F ′(z) = f(z) =Γ(ν+12 )√νπσΓ(ν2 )

(1 +(z − µ)2

νσ2)−

ν+12 , z ∈ R,

onde Γ(.) é a função gamma. A distribuição t-Student vem sendo bas-tante utilizada na presença de observações extremas (”outliers”) e uma re-visão bastante completa de sua mais importantes propriedades aparece emArellano-Valle e Bolfarine (1995). Uma versão assimétrica do modelo t-Student é considerada em Gomea et al. (2007). Veja também Arellano-Valleet al. (2012).

O modelo de regressão t-Student pode ser escrito através da hierarquia

Yi|Vi = viind∼ N(x′iβ, v−1i σ

2),

Viiid∼ Gama(ν/2, ν/2),

i = 1, . . . , n.

Para implementar o enfoque Bayesiano e algoritmo EM, pode-se usar averossimilhança completa

Lc(θ) =n∏i=1

(1√

2πσ2)e−

12σ2

vi(yi−x′iβ)2 .(ν2 )

ν2

Γ(ν2 )vν2−1

i e−vi ν2 .

Uma outra distribuição que pode ser usada é a distribuição α-potência(Pewsey et al., 2012). Uma variável aleatória Y segue a distribuição α-potência com parâmetros α, µ, σ, que denotamos por Y ∼ AP (µ, σ, α) sesua fdp é dada por

f(y|µ, σ, µ) = 1σf(y − µσ

){F (y − µσ

)}α−1,

onde α ∈ R. Para o caso normal, isto é F = Φ, temos o modelo potêncianormal. Usamos a notação Y ∼ PN(µ, σ, α). Este modelo é proposto comouma alternativa ao modelo ”skew-normal”, com fdp

f(y|λ) = 2φ(y)Φ(λy),


que apresenta algumas dificuldades como a singularidade da MIF (Azza-lini, 1985). Isto implica, por exemplo, que a distribuição da estat́ıstica darazão de verossimilhanças (ERV) não é assintoticamente distribúıda comdistribuição χ2.

Para o modelo potência-normal pode-se mostrar que a matriz de in-formação de Fisher para θ = (µ, σ, α = 1)′ é dada por

IF (θ) =

1σ2 0 0.9031920 2σ2

−0.5956360.903197 −0.595636 1

.Pode-se verificar que

|IF (θ)| = 0.013688/σ4.Então, para este modelo a matriz de informação de Fisher não é singularno ponto de simetria. Por outro lado, Azzalini (1985) mostra que o modelo”skew normal”apresenta matriz de informação de Fisher singular. Isso im-plica que as condições usuais de regularidade (Bolfarine e Sandoval, 2005)não estão satisfeitas neste caso.

Para o caso do modelo tobit (com T=c), a log-verossimilhança para omodelo α-potência pode ser escrita como

`(θ;Y ) = α∑i

(1− Ii) log[F

(c− µσ

)]+∑i

Ii

{log(α)− log(σ) + log

(f

(yi − µσ

))+ (α− 1) ln

(F

(yi − µσ

))},

onde

Ii =

{1, se yi > c0, se yi ≤ c.

2.3 Aplicações

Nesta seção consideramos aplicações ao conjunto de dados usando o modelotobit e o modelo tobit potência-normal. Os resultados a seguir aparecemem Martinez et al. (2013).

2.4 APLICAÇÃO DO MODELO TOBIT POTÊNCIA-NORMAL 13

2.3.1 Aplicação do modelo tobit-normal

Vamos ilustrar uma aplicação do modelo tobit-normal para parte dos dadosem Fair (1978). Para uma amostra de 601 homens e mulheres casados pelaprimeira vez, temos como variável resposta (Y), o número de casos extra-conjugais. Parte das variáveis usadas no estudo foram as seguintes:

Y : número de casos extraconjugais no ano anterior

X1: anos de casado

X2: idade

X3 : religosidade (escala de 1 (ateu) a 5 (frequenta regularmente)

X4: avaliação casamento (escala de 1 (muito infeliz) a 5 (muito feliz)

Dos 601 entrevistados, 451 não tiveram casos. Temos, portanto, dadoscom censura em zero.

Tabela 2.1: Estat́ısticas descritivas para dados de Fairn Média Variância Assimetria curtose

601 7.45 17.11 0.1553 3.7

Note que existe indicação de assimetria e curtose acima do esperadopara a distribuição normal. Temos também as estimativas: β̂0 = 9.08 (2.66),β̂1 = −0.16 (0.077), β̂2 = 0.54 (0.13), β̂3 = 1.72 (0.41), β̂4 = −2.26 (0.41),σ = 8.27 (0.55). Além disso, Log − lik = −706.4. Portanto as variáveisinfluenciam significativamente no número de casos extraconjugais.

2.4 Aplicação do modelo tobit potência-normal

Para os dados de Fair (1978), usando o modelo tobit potência-normal (to-bit/PN) temos α̂ = 10.26 (0.56), com Log − lik = −581.22, indicando forteevidência de que o modelo tobit/PN apresenta melhor ajuste para os dados


de Fair. Não existem disferenças significativas nas estimativas dos outrosparâmetros.

Caṕıtulo 3

O modelo tobit com excessode zeros

Neste caṕıtulo consideramos modelos para a situação onde temos ajuste domodelo tobit com excesso de zeros. Basicamente, consideramos os mode-los propostos em Moulton and Halsey (1995) e Cragg (1971). Discutimosestimação por métodos clássicos e Bayesianos.

3.1 Modelos com excesso de zeros

Existem situações reais onde a quantidade de zeros é maior que o esper-ado com o modelo tobit-normal. Uma possibilidade é considerar que partedos zeros observados vem de uma massa pontual concentrada no limite dedetecção mı́nimo (LDM) não explicada pela distribuição correspondente àresposta não nula.

O modelo tobit com excesso de zeros pode ser implementado considerandoo enfoque em Moulton e Halsey (1995) que especifica para a resposta obser-vada que

g(yi) = [qi + (1− qi)F (T )](1− Ii) + (1− qi)f(yi)Ii,

onde

15

16 O MODELO TOBIT COM EXCESSO DE ZEROS 3.2

Ii =

{1, se yi > T0, se yi ≤ T.

A situação onde o ponto de truncamento é T = 0 é imediata. Covariadaspodem ser associadas com qi através de uma função de distribuição (ligação)H, ou seja,

qi = H(x′iβ).

Para o caso em que qi = q, isto é, a probabilidade de excesso de zeros éconstante para as unidades amostrais, a função log-verossimilhança parauma amostra y = (y1, . . . , yn)

′ e vetor de parâmetros θ, pode ser escritacomo

l(θ|y) ∝n∑j=1

{(Ij − 1) log(q + (1− q)F (yj)) + Ij log(1− q) + log(f(yj))}.

Uma alternativa ao modelo de Moulton and Halsey (1995) é a alternativaproposta por Cragg (1971) onde é especificado que

g(yi) = qi(1− Ii) + (1− qi)f(yi)Ii,

i = 1, . . . , n, ou seja, os zeros observados são oriundos da massa pontual.Note que o caso particular do modelo tobit padrão segue como um casoparticular dos modelos acima tomando qi = 0, i = 1, . . . , n.

3.2 A distribuição log-α-potência

Conforme visto no caṕıtulo anterior, o modelo tobit (potência) assimétrico(TPA) pode ser definido considerando F como sendo a fda da distribuiçãoα-potência com fdp dada por

fF (z) = αf(x)F (z)α−1.

No caso particular em que F ′ = f = φ, temos, como visto anteriormente,

fN (z) = αφ(z)Φα−1(z).

A distribuição do tempo de vida de um equipamento e a concentração deum elemento qúımico em amostras de solo (água ou sangue) é tipicamentedistribúıda de acordo com a distribuição log-normal. Em muitas dessas

3.2 A DISTRIBUIÇÃO LOG-α-POTÊNCIA 17

situações, contudo, a assimetria da distribuição pode estar acima do esper-ado com a distribuição log-normal.

O modelo log-”skew-normal”é estudado em Gomez et al. (2011) do qual omodelo log-normal é um caso especial. Uma extensão do modelo log-normalpara o modelo log-”skew-normal”é considerado em Chai e Bailey (2008).Contudo uma das dificuldades do modelo log− skew− normal (log-normalassimétrico) é o fato de sua matriz de informação de Fisher ser singular. Adistribuição da estat́ıstica da razão de verossimilhança para testar normali-dade, por exemplo, não segue distribuição χ2.

Como uma alternativa a estas situações, consideramos o modelo log-potência-normal (”log-power-normal”) (LPN), que contém como caso par-ticular a distribuição log-normal. Uma vantagem deste modelo é que elecontém um parâmetro de forma adicional, que o faz mais flex́ıvel em termosde assimetria e curtose para ajustar dados experimentais como os consider-ados nestas notas.

Dizemos que uma variável y, com suporte em R+, segue uma distribuiçãolog-α-potência univariada com parâmetro α, que denotamos por Y ∼ LAP (α),se a variável transformada X = log(Y ) ∼ AP (α).

A fdp de uma variavel Y ∼ LAP (α) pode ser escrita como

g(y;α) =α

yf(log(y)) {F (log(y))}α−1 , y ∈ R+, (3.1)

onde F é uma função de distribuição absolutamente cont́ınua com funçãode densidade f = F ′. Nos referimos a esta distribuição como log-α-potênciapadrão.

No caso especial em que f = φ(·) e F = Φ(·), as funções de densidade ede distribuição da normal padrão (N(0, 1)), respectivamente, a distribuiçãolog-potência-normal segue, com fdp dada por

g(y;α) =α

yφ(log(y)) {Φ(log(y))}α−1 , y ∈ R+, (3.2)

que denotamos por Y ∼ LPN(α). Sua função de distribuição pode serescrita como

FY (y;α) = {Φ(log(y))}α , y ∈ R+. (3.3)

O método de inversão pode ser usado para gerar valores aleatórios davariável com distribuição LPN(α). Isto é, se U ∼ U(0, 1), a distribuição


da variável aleatória Y = eΦ−1(U1/α) é (log-potência-normal) LPN com

parâmetro α.Seja X ∼ PN(µ, σ, α), onde µ ∈ R é um parâmetro de localização e

σ ∈ R+ é um parâmetro de escala. Então, a transformação X = log(Y )leva ao modelo localização-escala log-potência-normal. Usamos a notaçãoY ∼ LPN(µ, σ, α).

No caso particular em que α = 1, isto é,

Z =log(Y )− µ

σ∼ N(0, 1),

pode-se mostrar que (depois de algumas manipulações algébricas que) amatriz de informação de Fisher para θ = (µ, σ, α)′ é dada por

I(θ) =

1/σ2 0 a01/σ0 2/σ2 a11/σa01/σ a11/σ 1

,onde akj = E{zk(φ(z)/Φ(z))j} for k = 0, 1, 2, 3 e j = 1, 2, que coincidecom a matriz de informação de Fisher para a distribuição potência-normal(Pewsey et al., 2012).

Assim, usando procedimentos numéricos, pode-se mostrar que

|IF (θ)| = [2− (a211 + 2a201)]/σ4 6= 0,

de modo que a matriz de informação de Fisher é não singular para α = 1.0. Amatriz de informação completa também foi derivada. Então, para n grande,

θ̂A−→ N3(θ, IF (θ)−1),

implicando na consistência e normalidade assintótica do EMV de θ, cujavariância assintótica d́ada por IF (θ)

−1.Como consequência desta propriedade importante, podemos testar (com

o modelo LPN) log-normalidade (isto é, H0 : LPN = LN), usando pro-priedades para grandes amostras da estat́ıstica da RV que segue distribuiçãoχ2. Este não é o caso, por exemplo da distribuição LSN, para a qual a MIFé singular. A escolha de um modelo conveniente pode ser feito através dosvalores de assimetria e curtose.

3.3 O modelo bernoulli/log-α-potência

Uma extensão importante do modelo log-potência-normal para a situaçãode excesso de zeros é a extensão proposta em Cragg (1971), usualmente

3.3 O MODELO BERNOULLI/LOG-α-POTÊNCIA 19

chamado modelo de duas partes (two-part model), que estabelece uma maneirade relaxar a restrição do truncamento no modelo tobit. Sob o modelo Cragg(1971) a fdp de yi pode ser formalmente escrita como

g(yi) = piIi + (1− pi)f(yi)(1− Ii), (3.4)

onde pi é a probabilidade determinando a contribuição relativa da massapontual na distribuição da mixtura, f é uma fdp com suporte positivo e,

Ii =

{0, se yi > 01, se yi ≤ 0.

Neste modelo os dois componentes são determinados por processos es-tocásticos diferentes de modo que os componentes positivos vem da fdp f .Por outro lado um zero vem da massa pontual. Este modelo não consid-era contudo um limite de deteção mı́nimo e que parte das observações estáabaixo deste limite.

Moulton e Halsey (1995) generalizam o modelo em duas partes per-mitindo que parte das resposta limites resultam de censura intervalar de f .Isto significa que um zero pode vir da massa pontual ou pode ser um valorde f não definido precisamente em (0, T ), com T constante. Formalmente,

g(yi) = [pi + (1− pi)F (T )]Ii + (1− pi)f(yi)(1− Ii), (3.5)

onde F é a fda de f.

Então, uma grande quantidade de modelos são produzidos variando adensidade básica f e a função de ligação correspondente a pi. Diversosmodelos h́ıbridos podem ser considerados como os modelos probit/potência-normal, logit/log-normal, logit/log-gamma e probit/log-skew-normal. Estesmodelos foram considerados em aplicações práticas em biologia, economia,agricultura e muitas outras àreas (Chai and Bailey, 2008). Note que sepi = 0, i = 1, . . . , n, o modelo de Moulton e Halsey (1995) reduz-se ao mod-elo tobit usual (Tobin, 1958).

No caso da medição de concentração de anticorpos por diferentes lab-oratórios, e considerando yi a resposta para a unidade i, é tipicamente deinteresse a situação onde a distribuição de log(yi) é função dos parâmetrosβ0, ..., βp que estão relacionados através do modelo linear

log(yi) = β0 + β1x1i + ...+ βpxpi + �i,


onde �i ∼ PN(0, σ, α) e x1, ..., xp são constantes fixas e conhecidas.

Sob o modelo PN,

E[�i] = ασ

∫ 10

Φ−1(z)zα−1dz 6= 0,

de modo que o valor esperado do termo do erro não é nulo como é o casosob normalidade.

Consequentemente, E[yi] 6= x′iβ e teremos que corrigir o parâmetro inter-cepto, isto é, β∗0 = β0 + µ�, onde µ� = E[�i]. Então,

E[yi] = x′iβ∗, onde β∗ = (β∗0 , β1, ..., βp)

′.

Consideramos agora extensões do modelo Bernoulli/LN para as situaçõesdos modelos logito/LPN e probito/LPN, juntamente com covariadas emcada passo do modelo. Este desenvolvimento está apresentado em Martinezet al. (2012a).

Inicialmente, suponhamos que todas as observações vem do modelo LPNcom parâmetros de localização e escala µ e σ, respectivamente, mas semcovariadas. A contribuição para a verossimilhança de observações não cen-suradas, isto é, para y > T, pode ser representada como

α

σyφ [(log(y)− µ)/σ] {Φ[(log(y)− µ)/σ]}α−1 .

Covariadas são introduzidas para ambas as partes do modelo, ou seja,para as variáveis D e Y , de modo que considerando a ligação logito para avariável D temos que

logit{P [D = 1|x(1)]} = x′(1)β(1),

onde x(1) é o vetor de covariáveis de dimensão p, associados com o vetor deparâmetros β(1). Então, temos que

τi = 1− pi =exp(x′(1)iβ(1))

1 + exp(x′(1)iβ(1)), i = 1, . . . , n.

3.3 O MODELO BERNOULLI/LOG-α-POTÊNCIA 21

Correspondendo a parte LPN temos o vetor de covariáveis x(2) de di-mensão q, possivelmente diferente de x(1), onde temos o vetor de parâmetrosβ(2), para os quais

log(yi) ∼ PN(x′(2)iβ(2), σ, α), yi > 0.

Chamamos atenção para o fato que diferentes distribuições podem levar amodelos de regressão mais informativos (Chai and Bailey, 2008).

O logaritimo da função de verossimilhança para θ = (β′(1)β′(2), σ, α)

′

dados X = (x1, x2) e Y = (y1, . . . , yn), desprezando constantes não infor-mativas, pode ser escrita como

`(θ;X,Y ) =∑i

Ii{log[1 + exp(x′(1)iβ(1)){Φ(zT i)}α]

− log[1 + exp(x′(1)iβ(1))]}

+∑i

(1− Ii){log(α)− log(σyi)

+x′(1)iβ(1) − log[1 + exp(x′(1)iβ(1))

]− 1

2z2i + (α− 1) log(Φ(zi))},

onde zT i =log(T )−x′

(2)iβ(2)

σ e zi =log(yi)−x′(2)iβ(2)

σ .

Usando as equações acima, estimadores de máxima verossimilhança paraos parâmetros do modelo podem ser calculados. Como a MIF para o mod-elo LPN é não singular, inferência em grandes amostras para o modeloBernoulli/LPN podem ser implementadas para os EMV sob condições deregularidades usuais onde o EMV é assintoticamente normal com médiaθ e matriz de covariâncias igual a inversa da MIF, indicando otimalidadeassintótica. Pode-se considerar extensões do modelo acima como a presençade interações.

Considerando agora o modelo probit para a variável de Bernoulli D,temos que

pi = P [yi = 0] = Φ(−x′(1)iβ(1)) = 1− Φ(x′(1)iβ(1))

e

log(yi) ∼ APN(x′(2)iβ(2), σ, α), yi > 0.


O logaritimo da função de verossimilhança (função log-verossimilhança),a menos de constantes, pode ser escrito como

`(θ;X,Y ) =∑i

Ii

{log[1 + Φ(x′(1)iβ(1)){{Φ(zT i)}

α − 1}]}

,

+∑i

(1−Ii){

log(α)− log(η) + log(

Φ(x′(1)iβ(1)))− 1

2z2i + (α− 1) log(Φ(zi))

},

onde

zT i =log(T )− x′(2)iβ(2)

σe zi =

log(yi)− x′(2)iβ(2)σ

.

A função escore é obtida derivando-se a função de log-verossimilhança.

A função log-verossimilhança do modelo tobit (com T=c) considerandoque a distribuição do erro segue distribuição α-potência pode ser escritacomo

`(θ;Y ) = α∑i

(1− Ii) log[F

(c− µσ

)]+

∑i

Ii

{log(α)− log(σ) + log

(f

(yi − µσ

))+ (α− 1) ln

(F

(yi − ξσ

))}onde

Ii =

{1, se yi > c0, se yi ≤ c.

Casos particulares importantes seguem tomando f = φ e f = tν(µ, σ2).

3.4 Aplicação: dados do Haiti

Consideramos a ligação logito e a distribuição log-normal para parte pos-itiva (incluindo respostas limitadas). Os dados são descritos em Moultonand Halsey (1995).

Tabela 1 sintetiza resultados de estimação para os dados de vacinaçãono Haiti sob diferentes modelos considerando ou não mistura e censura.

3.4 APLICAÇÃO: DADOS DO HAITI 23

Variáveis:

EZ (Tipo de vacina, 0: Schwarz, 1: Edmonston-Zagreb);

HI (dose, 0: médio, 1: alto);

FEM (sexo; 0: masculino, 1: feminino);

INT: Termo constante.

A tabela a seguir apresenta análises classica (EMV) e Bayesiana para osdados acima, considerado o modelo Bernoulli/log-normal. As estimativasdas variâncias para o enfoque clássico são apresentadas em Moulton andHalsey (1995) de onde se conclui que das variáveis consideradas no estudo,TIPO e SEXO são significantes.


Tabela 3.1: Estimativas clássicas e BayesianasModelo Método Componente Bernoulli Componente log-normal Component

INT EZ HI FEM INT EZ HI FEMA Clas -0.979

Bay -0.981B Clas -1.287 0.340 0.182 0.115

Bay -0.932 0.203 0.097 0.114C Clas 1.198 -0.273

Bay 1.227 -0.285D Clas 1.178 -0.327 -0.109 -0.037 0.290

Bay 1.226 -0.361 -0.083 -0.025 0.277E Clas 0.732 0.843 0.431 -0.166 -0.274

Bay 0.813 0.950 0.445 -0.244 -0.305F Clas 0.765 0.932 0.433 -0.281 -0.304 -0.192 -0.063 0.329

Bay 0.910 1.112 0.439 -0.425 -0.353 -0.199 -0.055 0.339G Clas 0.648 0.830 0.426 -0.404 0.279

Bay 0.678 0.893 0.440 -0.421 0.266

Tabela 3.2: Ajustes MV e BayesianosModelo −2× loglik DIC pD

A 1115.830 136.600 1.89B 1113.180 120.560 5.17C 1079.320 101.800 2.7D 1075.620 104.500 5.79E 1068.720 95.560 5.08F 1063.360 94.470 9.07G 1065.810 93.840 5.42

Estimadores dos parâmetros para ajustes da mistura logito/LN com ume dois componentes considerando inferência clássica e Bayesiana para osdados do Haiti.

Comparações para dados do Haiti considerando inferência clássica eBayesiana. Note que existe discordância entre os resultados clássicos eBayesianos quanto ao ajuste do modelo. Para o enfoque Bayesiano, o melhormodelo é o modelo G (mais completo), enquanto que para o enfoque clássicoo modelo que melhor se ajusta é o modelo F.

A tabela a seguir apresenta resultados do ajuste Bayesiano dos mode-los log-normal e log-potência-normal incluindo as estimativas dos desviospadrões. Note que o enfoque Bayesiano tanto para os modelos log-normalcomo log-potência-normal indicam significância das variáveis TIPO e SEXO.

Para o modelo completo, Moulton e Halsey (1995) obtiveram os seguintesestimativas (Estimativa/DP):

3.5 APLICAÇÃO: DADOS DE MROZ 25

Tabela 3.3: Ajustes log-potência-normal e log-normalModel Log-Normal Log-Potência-Normal

parameters mean MC error P5 P95 mean MC error P5 P95α 16.69 0.600 3.55 38.38

β(1)0 0.91 0.009 0.42 1.45 0.72 0.009 0.30 1.18β(1)1 1.15 0.058 0.50 1.92 0.86 0.011 0.39 1.35β(1)2 0.44 0.009 -0.06 0.99 0.38 0.008 -0.07 0.85β(1)3 -0.42 0.009 -1.02 0.13 -0.26 0.009 -0.74 0.20β(2)0 -0.35 0.004 -0.66 -0.07 -3.43 0.047 -4.94 -1.70β(2)1 -0.20 0.005 -0.48 0.08 -0.14 0.005 -0.37 0.10β(2)2 -0.06 0.003 -0.34 0.22 0.01 0.005 -0.21 0.26β(2)3 0.35 0.003 0.07 0.63 0.25 0.006 0.01 0.50σ 1.18 0.003 1.06 1.32 1.87 0.009 1.48 2.21τ 0.73 0.003 0.57 0.89 0.30 0.004 0.21 0.46

Dbar 7687.00 7681DIC 7693.00 7687EAIC 7705.00 7701EBIC 7739.19 7739.0

Componente Bernoulli: β̂(1)0 = .77(2.77), β̂(1)1 = .93(2.82), β̂(1)2 =

.43(1.48), β̂(1)3 = −.28(2.82)

Componente log-normal: β̂(2)0 = −.31(−1.89), β̂(2)1 = −.19(−1.20),β̂(2)2 = −.06(−.40), β̂(2)3 = −.33(2.06).

Temos, portanto que os resultados clássicos e Bayesianos concordamquanto a significância dos parâmetros, havendo contudo diferença no melhormodelo ajustado. O enfoque Bayesiano recomenda o modelo G.

3.5 Aplicação: dados de Mroz

Consideramos os dados de Mroz (1987), que analisa as informações de 753mulheres casadas com idade entre 30 e 60 anos, com interesse na relaçãoentre a oferta de trabalho e outras covariáveis, no ano de 1975. Para obteros dados, basta entrar no R com

> library(sampleSelection)

> data(Mroz87)


As variáveis utilizadas no artigo são: Horas de trabalho (variável re-sposta), salário que não é devido ao trabalho da mulher, anos de educação,anos de experiência de trabalho, idade da mulher, número de crianças menoresque 6 anos, nḿero de crianças entre 6 e 18 anos.

Tabela 3.4: Estimadores Bayesianos para parâmetros do componenteBernoulli

Parâmetro Média D.P. Q2.5% Q97.5%

β1(1) -0.05 9.761 -19.31 19.27

β1(2) -0.54 9.68 -19.74 18.46

β1(3) 5.10 7.753 -12.41 19.73

β1(4) -3.80 6.283 -9.98 16.94

β1(5) 6.50 5.866 -8.25 14.48

β1(6) 11.90 5.417 0.023 18.94

β1(7) 1.54 11.62 -17.59 22.21

β1(8) 9.3 6.069 -0.61 20.6

Note que H0 : β1(6) 6= 0 é significante, de modo que existe indicação deque existe excesso de zeros nos dados de Mroz (1976).

Tabela 3.5: Estimadores Bayesianos para parâmetros do componentecont́ınuo

Parâmetro Média D.P. Q2.5% Q97.5%

β2(1) 0.8324 9.921 -19.0 20.02

β2(2) -5.715 9.885 -25.39 13.8

β2(3) 3.111 9.462 -15.05 21.65

β2(4) -8.74 3.444 -15.58 -1.978

β2(5) 23.23 8.355 6.486 39.15

β2(6) -6.308 4.128 -14.13 1.763

β2(7) 38.18 7.592 21.87 52.86

β2(8) 0.7323 0.3108 0.168 1.389

Temos também que σ̂ = 1223, 0. Note que váriáveis significantes para aparte cont́ınua são 1, 4, 5 e 6. Para a parte discreta (pontual), temos quea variável X5 é significativa ao ńıvel de 5%, indicando que existe excesso dezeros nos dados de Mroz.

Caṕıtulo 4

Modelo α-potênciainflacionado de zeros e/ouuns

Neste caṕıtulo consideramos distribuições potência para modelar proporçõesou taxas com inflação de zeros e/ou uns como uma alternativa ao mod-elo de regressão beta. Os modelos considerados são misturas de processosde Bernoulli para explicar o excesso de zeros e/ou uns e uma distribuiçãopotência-normal limitada para explicar a resposta cont́ınua. Consideramosos enfoques de máxima verossimilhança e Bayesiano para a estimação dosparâmetros. Matrizes de informação observadas (MIO) e esperadas (MIF)são derivadas, ilustrando aspectos interessantes destes modelos.

Dada a flexibilidade da distribuição potência-normal, pode-se mostrarem um cenário prático que o modelo tobit modificado pode ser mais precisoque o modelo de regressão beta.

4.1 Modelos duplamente censurados

Modelos estat́ısticos usados para explicar variáveis respostas no intervalo(0, 1) tem recebido considerável atenção na literatura estat́ıstica recente.Entre outros, mencionamos, Ferrari e Cribari-Neto (2004), Brascum et al.(2007) e Bayes et al. (2012). Extensões deste modelos para situações comrespostas no intervalos [0, 1], [0, 1) e (0, 1] são estudadas em Ospina e Ferrari

27

28 MODELO α-POTÊNCIA INFLACIONADO DE ZEROS E/OU UNS 4.2

(2010). Variáveis deste tipo incluem, por exemplo, a proporção de mortescausadas pelo cigarro, a proporção de impostos gastos na educação, a pro-porção de renda familiar gasta em alimetação, etc.

A situação da variável resposta com inflação de zeros e uns é relatado emum conjunto de dados sobre a porcentagem de mortes não explicadas nosmunićıpios brasileiros durante o ano 2000 entre crianças com menos de umaano de idade. Das 5561 observações coletadas, tem-se um total de 3367 zerose 174 uns, que certamente deve ser incorporado no estudo. Para tratar destecenário mais complexo uma extensão do modelo de regressão beta usual foiconsiderado in Ospina (2008) e Ospina e Ferrari (2010), levando a resultadosbastante satisfatórios.

Neste caṕıtulo, propomos um enfoque alternativo ao descrito acima.Ele é uma extensão do modelo tobit censurado (Tobin, 1956) no inter-valo [0, 1], para incorporar inflação de zeros e/ou uns. É considerado queparte dos zeros e/ou uns vem de uma variável Bernoulli ligando posśıveisexcessos de zero e/ou uns com um grupo de covariáveis que podem influ-enciar na probabilidade de de ocorrência de tais valores. Por outro lado,as resposta cont́ınuas podem ser modeladas usando a distribuição potência-normal (Gupta e Gupta, 2008, Pewsey et al., 2012), que são mais flex́ıveisque a distribuição normal em termos de assimetria e curtose com EMVs bemcomportados para os quais as condições de regularidade estão satisfeitas.

Além disso, a extensão do modelo tobit que propomos consiste em sub-stituir a fda da distribuição normal pela fda da distribuição potência-normalque é quase tão simples de se trabalhar quanto o modelo normal usual. Umaalternativa é usar a distribuição normal assimétrica que apresenta as dificul-dades já mencionadas anteriormente e além disso tem fda não tão simplesde ser trabalhada.

Definimos inicialmente o modelo tobit-potência-normal (TPN) dupla-mente censurado no intervalo (0, 1), extendendo o modelo tobit usual parasituações duplamente censuradas. A seguir o modelo é extendido parasituações com excesso de zeros e/ou uns. Situações com dados reais são anal-isadas. Introduzimos o modelo Bernoulli/tobit-potência-normal (Bernoulli/TPN),onde se trata o problema de estimação do ponto de vista Bayesiano.

4.3 DISTRIBUIÇÕES PN PARA DADOS CENSURADOS 29

4.2 Distribuições PN para dados censurados

Em uma situação duplamente censurada, a variável resposta é restrita atomar valores em um intervalo, e eventualmente pode tomar os valores lim-ites para parte significante dos dados. Os valores limites são usualmentechamados de limites de deteção mı́nimo (LDm) e máximo (LDM), respecti-vamente. Temos então o modelo tobit duplamente censurado.

O modelo tobit usual pode não ser adequado em situações onde os valoresobservados para a parte cont́ınua dos dados apresentam assimetria e curtosemaior do que é esperado para o modelo normal. Em tais situações, o modelopotência-normal pode ser uma alternativa viável.

4.3 Modelo potência-normal duplamente censurado

Suponhamos que y∗ ∼ PN(ξ, η;α). Considere uma amostra de tamanho n,(y∗1, y

∗2, ..., y

∗n) e que somente parte dos valores de y

∗ está entre constantesc0 e c2. Para valores de y

∗ ≤ c0 somente o valor c0 é relatado enquantoque para valores de y∗ ≥ c2, somente o valor c2 é relatado. Podemos entãoescrever os dados observados como

yi =

c0, se y

∗i ≤ c0,

y∗i , se c0 < y∗i < c2,

c2, se y∗i ≥ c2,

i = 1, 2, ..., n.

A amostra resultante é dita ser uma amostra PN duplamente censurada.Para observações yi = c0, temos que

P [yi = c0] = P [y∗i ≤ c0] = {Φ (z0)}

α ,

onde z0 = (c0 − µ)/σ; com y∗i = c2 temos

P [yi = c2] = P [y∗i ≥ c2] = 1− {Φ (z2)}

α ,

onde z2 = (c2 − µ)/σ. Para respostas cont́ınuas, isto é, c0 < y∗i < c2, temosque yi ∼ PN(µ, σ, α). Denotamos esta variável por PNDC(µ, σ, α).

Particularmente, para α = 1, o modelo se reduz ao modelo tobit dupla-mente censurado.


Denotando por∑

0,∑

1 and∑

2, as somas correspondendo a y∗ ≤ c0,

c0 < y∗i < c2 e y

∗ ≥ c2 respectivamente, então, o logaritimo da função deverossimilhança correspondente a uma amostra de tamanho n para estimarθ = (µ, σ, α)′ pode ser escrita como

`(θ; Y) = α∑

0

log [Φ (z0)] +∑

2

log [1− {Φ (z2)}α]

+∑

1

{log(α)− log(σ) + log (φ (z1i)) + (α− 1) log (Φ (z1i))} ,

onde zi = (yi − µ)/σ, i = 1, . . . , n.

Portanto, os elementos da função escore são dados por

U(ξ) = − 1σ

∑0

r(z0) +1

σ

∑1

{z1i − (α− 1)w1i}+1

σ

∑2

h(z2),

U(η) = − 1σ

∑0

r(z0)z0 +1

σ

∑1

{−1 + z21i − (α− 1)z1iw1i

}+

1

σ

∑2

z2h(z2),

U(α) =∑

0

log [Φ (z0)]+∑

1

{1

α+ log (Φ (z1i))

}−α−1

∑2

log(Φ(z2))w−12 h(z2),

onde

z0 =c0 − µσ

, z2 =c2 − µσ

, z1i =yi − µσ

, w2 =φ(z2)

Φ(z2), w1i =

φ(z1i)

Φ(z1i),

e h e r são as funções de risco, r(t) = φ(t)/Φ(t), e risco inverso h(t) =φ/(1− Φ(t)).

Pode-se mostrar que as elementos da matriz de informação observadasão dados por

jµµ =1

η2

∑0

r(z0){z0 + α−1r(z0)}

+1

σ2

∑1

{1 + (α− 1)[z1iw1i + w21i]}

+1

σ2

∑2

{h(z2)[−z2 + (α− 1)w2 + h(z2)]},

4.4 MODELO POTÊNCIA-NORMAL DUPLAMENTE CENSURADO 31

jσµ =1

η2

∑0

r(z0){−1 + z20 + α−1z0r(z0)}

+1

σ2

∑1

{2z1i + (α− 1)[−w1i + z21iw1i + z1iw21i]}

+1

σ2

∑2

{h(z2)[1− z22 + (α− 1)z2w2 + z2h(z2)]},

jσσ =1

σ2

∑0

r(z0){−2z0 + α−1z20r(z0) + z30r(z0)}

+1

σ2

∑2

{z2h(z2)[2− z22 + (α− 1)z2w2 + z2h(z2)]}

1

σ2

∑1

{−1 + 3z21i + (α− 1)[−2z1iw1i + z21iw21i + z31iw1i]},

jαµ =1

ασ

∑0

r(z0) +1

σ

∑1

w1i −1

σ

∑2

{h(z2)[α−1

+ log(Φ(z2))[1 + w2]]},

jασ =1

ασ

∑0

z0r(z0)

+1

σ

∑1

z1iw1i −1

σ

∑2

{z2h(z2)[α−1 + log(Φ(z2))[1 + w2]]},

jαα =1

α2

∑1

1 + α−2∑

2

{w−22 log(Φ(z2))h(z2)[αw2 + h(z)]}.

Baseado na função escore, os elementos da matriz de informação observadados parâmetros do modelo podem ser estimados usando algoritmos itera-tivos.

A MIF segue tomando-se esperanças dos componentes acima (multiplicadospor n−1), é importantante no sentido de que a distribuição assintótica doestimador de máxima verossimilança é normal com variância assintótica queé o o inverso da MIF. Temos também que a MIF é não singular.


4.4 A transformação logaŕıtmica

No caso de variáveis respostas tomando somente valores positivos, podemosconsiderar a transformação Z = log(Y ), onde Z ∼ N(µ, σ2).

Considerando agora que Z ∼ PN(µ, σ, α), nos obtemos o modelo log-potência-normal com parâmetros µ, σ e α, denotado por Y ∼ LPN(µ, σ, α).A fdp para este modelo pode ser escrita como: ϕLPN (y;µ, σ, α) = ϕΦ(log(y);µ, σ, α)/y,y > 0. A fda correspondente é dada por FY (y;α) = {Φ((log(y) − µ)/σ)}α.Se os dados censurados em [0,∞), com alta assimetria positiva podemossubstituir y por y + 1 dado que o logaritmo de c0 = 0 não existe.

Para dados duplamente censurados usamos a notação LPNDC(µ, σ, α).A função log-verossimilhança para o modelo LPNDC com c0 = 0 é dado

por

`LPN (θ; Y) = −∑

1

log(y + 1) + `(θ; log(Y + 1)),

onde `(.) é a log-verossimilhança para o modelo PNDC, com z0 = −µ/σ,z1i = (log(yi + 1) − µ)/σ e z2 = (log(c2 + 1) − µ)/σ. A função escore ea matriz de informação observadas podem ser obtidas das correspondentespara o modelo PNDC, substituindo h(z2) por hLPN (z2) = h(log(c2 +1))/y er(z0) por rLPN (z0) = r(z0)/y onde h(.) e r(.) são as funções de risco e riscoinverso do modelo PN.

4.5 O modelo Bernoulli duplamente censurado commistura potência-normal

Para as variáveis resposta distribúıdas no intervalo [0, 1] (c0 = 0 e c2 = 1)o modelo tobit duplamente censurado pode não ser ótimo porque o excessode zeros e uns pode requerer modelos assimétricos capazes de captar taiscaracteristicas especiais.

Introduzimos então o modelo de mistura entre as variáveis resposta dis-creta e cont́ınuas que segue o modelo potência-normal.

Consideramos que a massa pontual no zero pode ser modelada por umavariável de Bernoulli com parâmetro γ, isto é, Ber(y; γ), e que a respostano intervalo (0, 1) pode ser modelada por uma distribuição α-potência (oulog-α-potência) com parâmetro θ = (µ, σ, α)′. A fdp correspondente paraeste modelo pode ser escrita como

4.6 ESTIMAÇÃO 33

g(yi) =

p(1− γ), se yi = 0,(1− p) ϕF (yi,µ,σ,α){F (z2)}α−{F (z0)}α , se 0 < yi < 1,pγ, se yi = 1,

onde 0 < p, γ < 1, σ, α > 0 e µ ∈ R.

Temos também que se ϕF (yi, µ, σ, α) denota a fdp da distribuição potência-normal. Como consequência da construção acima pode-se notar que P [y =0] = p(1− γ) e P [y = 1] = pγ. A fda de yi pode ser escrita como

FY (yi;µ, σ, α) =

p(1− γ), se yi ≤ 0,p(1− γ) + (1− p) {F (zi)}

α−{F (z0)}α{F (z2)}α−{F (z0)}α , se 0 < yi < 1,

1, se yi ≥ 1.

4.6 Estimação

Consideramos inicialmente que F = Φ, a fda da distribuição normal, demodo que temos uma mistura entre a variável aleatória de Bernoulli comparâmetro γ e a distribuição PN(µ, σ, α)). Denotamos este modelo porMBPN(p, γ, µ, σ, α). Logo, para uma amostra de tamanho n, y = (y1, . . . , yn)

T

da distribução MBPN(p, γ, µ, σ, α), denotamos por n0 =∑n

i=1 I0(y), n1 =∑ni=1 I1(y) e n01 =

∑ni=1 I0,1(y), onde IA(y) é a função indicadora do con-

junto A.

Assim, a função log-verossimilhança para θ = (p, γ, µ, σ, α) dado Y podeser escrita como:

`(θ; Y) = n01 log(p) + (n− n01) log(1− p) + n1 log(γ) + n0 log(1− γ)∑1

{log(α)− log(σ) + log(φ(zi)) + (α− 1) log(Φ(zi))

− log({Φ(z2)}α − {Φ(z0)}α)},

onde, zi = (yi − µ)/σ, i = 1, . . . , n.

Portanto, usando um enfoque similar ao de Pewsey et al. (2012), aprimeira derivada com respeito a p, γ, µ, σ e α pode ser escrita como


U(p) =n01p− n− n01

1− p,

U(γ) =n1γ− n0

1− γ,

U(ξ) = (n− n01){z − (α− 1)w

η+ϕΦ(c2, θ)− ϕΦ(c0, θ){Φ(z2)}α − {Φ(z0)}α

},

U(η) = −(n− n01)

{1− z2 + (α− 1)zw

η− z2ϕΦ(c2, θ)− z0ϕΦ(c0, θ){Φ(z2)}α − {Φ(z0)}α

},

U(α) = (n− n01){u+

1

α− {Φ(z2)}

α log(Φ(z2))− {Φ(z0)}α log(Φ(z0)){Φ(z2)}α − {Φ(z0)}α

},

onde wi = φ(zi)/Φ(zi) e ui = log{Φ(zi)}, i = 1, . . . , n.Então, o EMV para o parâmetro θ = (µ, σ, α)′, é obtido resolvendo o sistemade equações que seguem de igualar os escores acima a zero.

Então, obtemos as soluções para p̂ = n01n , γ̂ =n1n01, correspondendo,

respectiveamente, a proporções de zeros e uns na subamostra de zeros e uns.Segue que p̂ é um estimador não viciado para p. Para θ1 = (µ, σ, α)

′, osistema de equações não tem solução anaĺıtica, sendo portanto resolvida pormétodos numéricos.

4.7 Matriz de informação observada

Calculando a derivada segunda da log-verossimilhança obtemos os elementosjpp, jγp, jγγ , jξξ, jξη, . . . , jαα, dados em Martinez et al. (2012b).

Pode-se mostrar que a matriz de informação esperada (MF) para θ =(p, γ, µ, σ, α)′ é dada por

I(θ) = (1− p)

1

p(1−p)2 0 0 0 0

0 pγ(1−γ)(1−p) 0 0 0

0 0 iµµ iµσ iµα0 0 iµσ iσσ iσα0 0 iµα iσα iαα,

4.7 MATRIZ DE INFORMAÇÃO OBSERVADA 35

onde os seus elementos são dados em Martinez et al. (2012b).

Deste resultado segue que os parâmetros (p, γ)′ e (µ, σ, α)′ são ortogo-nais, de modo que a MIF é ortogonal em blocos, e pode ser escrita como

I(θ) = Diag{Ip,γ , Iµ,σ,α}, onde Ip,γ = Diag{

1p(1−p) ,

pγ(1−γ)

}.

Portanto, para n grande,

θ̂A→ N5(θ,Σθθ),

implicando que θ̂ é consistente e assintoticamenete normal com matriz de co-variâncias assintóticas Σθθ = I(θ)

−1 = Diag{I−1p,γ , I−1µ,σ,α} = Diag{Σp,γ ,Σµ,σ,α}.Note que parâmetros nos blocos podem ser estimados separadamente.

A aproximação normal N5(θ,Σ(θ)) pode ser usada para construir inter-valos de cofiança para θr, com coeficiente de confiança γ = 1 − α que sãodados por θ̂r ∓ z1−α/2

√σ̂(θ̂r), com os EMV e quantis da normal correspon-

dentes.

Considerando a reparametrização δ1 = pγ e δ0 = pδ1 podemos escrevero modelo como

g(yi) =

δ0, se yi = 0,

(1− δ0 − δ1) ϕΦ(yi,ξ,η,α){Φ(z2)}α−{Φ(z0)}α , se 0 < yi < 1,δ1, se yi = 1,

onde 0 < δ0 = P [yi = 0], δ1 = prob[yi = 1] < 1 e 0 < δ0 + δ1 < 1.

A função log-verossimilhança para θ = (δ0, δ1, µ, σ, α)′ dado y é dada

por

`(θ; Y) = n0 log(δ0) + n1 log(δ1) + (n− n01) log(1− δ0 − δ1)+∑

1

{log(α)− log(σ) + log (φ (zi))

+(α− 1) log (Φ (zi))− log({Φ(z2)}α − {Φ(z0)}α)},

os elementos do escore são:


U(δ0) =n0δ0− n− n01

1− δ0 − δ1,

U(δ1) =n1δ1− n− n01

1− δ0 − δ1,

U(µ) = (n− n01){z − (α− 1)w

η+ϕΦ(c2, θ)− ϕΦ(c0, θ){Φ(z2)}α − {Φ(z0)}α

},

U(σ) = −(n− n01)

{1− z2 + (α− 1)zw

σ− z2ϕΦ(c2, θ)− z0ϕΦ(c0, θ){Φ(z2)}α − {Φ(z0)}α

},

U(α) = (n− n01){u+

1

α− {Φ(z2)}


}.

Das primeiras duas equações, obtem-se δ̂0 = n0/n, proporção de zeros eδ̂1 = n1/n, a proporções de uns na amostra. Parâmetros restantes devemser estimados numericamente.

A MIF pode ser escrita como

I(θ) = Diag{Iδ0,δ1 , Iµ,σ,α},

onde os elementos de Iδ0,δ1 são dados por

iδ0δ0 =1− δ1

δ0(1− δ0 − δ1), iδ1δ0 =

1

1− δ0 − δ1e

iδ1δ1 =1− δ0

δ1(1− δ0 − δ1),

com Iµ,σ,α computado para o modelo MBPN(p, γ, µ, σ, α). Também temosortogonalidade.

4.8 MODELOS CENSURADOS PARA INFLAÇÃO DE ZEROS E UNS 37

Para n grande,

θ̂A→ N5(θ,Σθθ),

com θ̂ consistente e assintoticamente normal, com

Σθθ = I(θ)−1 = Diag{I−1δ0,δ1 , I

−1µ,σ,α} = Diag{Σδ0,δ1 ,Σµ,σ,α}

a varıância do EMV em grandes amostras.

4.8 Modelos censurados para inflação de zeros euns

Casos particulares são inflação de uns e zeros separadamente. Para o casode inflação de zeros, temos

g(yi) =

{δ0, se yi = 0,

(1− δ0) ϕΦ(yi,µ,σ,α){Φ(z2)}α−{Φ(z0)}α , se 0 < yi ≤ 1.

onde 0 < δ0 = P [yi = 0] e 0 < δ0 < 1.

A função log-verossimilhança para θ = (δ0, µ, σ, α)′ dado y é dada por

`(θ; Y) = n0 log(δ0) + (n− n0) log(1− δ0)+∑

1

{log(α)− log(σ) + log (φ (zi))

+(α− 1) log (Φ (zi))− log({Φ(z2)}α − {Φ(z0)}α)},

de modo que os elementos da função escore são dados por

U(δ0) =n0δ0− n− n0

1− δ0,

U(µ) = (n− n0){z − (α− 1)w

σ+ϕΦ(c2, θ)− ϕΦ(c0, θ){Φ(z2)}α − {Φ(z0)}α

},

U(σ) = −(n− n0)

{1− z2 + (α− 1)zw


},

U(α) = (n− n0){u+

1

α− {Φ(z2)}


}.


Da primeira equação, obtemos o estimator δ̂0 = n0/n, a proporção de ze-ros na amostra. Os parâmetros remanecentes requerem metódos numéricos.

Para o caso de inflação de uns, temos

g(yi) =

{δ1, se yi = 1,

(1− δ1) ϕΦ(yi,µ,σ,α){Φ(z2)}α−{Φ(z0)}α , se 0 ≤ yi < 1,

onde 0 < δ1 = P [yi = 1] e 0 < δ1 < 1, levando a log-verossimilhança paraθ = (δ1, µ, σ, α)

′ dado y pode ser escrita como:

`(θ; Y) = n1 log(δ1) + (n− n1) log(1− δ1)+∑

1

{log(α)− log(σ) + log(φ(zi))

+(α− 1) log(Φ(zi))− log({Φ(z2)}α − {Φ(z0)}α)},

de modo que os elementos da função escore são dados por

U(δ1) =n1δ1− n− n1

1− δ1,

U(µ) = (n− n1){z − (α− 1)w

σ+ϕΦ(c2, θ)− ϕΦ(c0, θ){Φ(z2)}α − {Φ(z0)}α

},

U(σ) = −(n− n1)

{1− z2 + (α− 1)zw


},

U(α) = (n− n1){u+

1

α− {Φ(z2)}


}.

Da primeira equação, obtemos o estimador δ̂1 = n1/n, a proporção deuns na amostra. Os outros parâmetros são estimados numericamente.

4.10 MISTURA BERNOULLI/LPN 39

4.9 Mistura Bernoulli/LPN

Considerando agora ϕF (yi, µ, σ, α)′ como a fdp do modelo LPN, o modelo

Bernoulli/LPN é obtido, que denotamos por MBLPN(p, γ, µ, σ, α). O mod-elo é importante na modelagem de dados com mais assimetria e curtose queos correspondentes da distribuição normal.

A função de log-verossimilhança do modelo reparametrizado pode serescrita como

`MBLPN (θ; Y) = −∑

1

log(yi) + `(θ; log(Y )),

onde `(.) é a função de log-verossimilhança do modelo MBPN e log(Y ) =(log(y1), ..., log(yn))

′. A função escore são como dadas para o modelo MBPNmodel, onde zi = (log(yi)− µ)/σ, i = 1, . . . , n.

4.10 Ilustração com dados reais

Nesta seção illustramos a utilidade das distribuições LPNDC e MBLPN parao ajuste de dados reais. O conjunto de dados que analizamos correspondea proporção de mortes de crianças de menos de um ano por causa não es-clarecidas nos 5561 munićıpios Brasileiros. Dados estão dispońıveis para”download”no site http:www.datasus.gov.br. O conjunto de dados contém3367 zeros (mortes esclarecidas) e 174 uns (mortes não esclarecidas).

Ospina (2008), desenvolve um modelo baseado na regressão beta paramodelar este tipo de dados com inflação de zeros e/ou uns. Como em Os-pina (2008) assumimos a mistura de uma variável de Bernoulli para modelara parte discreta com a regressão beta para a parte cont́ınua (entre zero eum), que é denotada por BIZU(δ0, δ1, ξ, η). Para estimar os parâmetros domodelo BIZU, a rotina GAMLSS no programa R pode ser usado. Nós de-senvolvemos programas no R para ajustar modelos LPNDC e para o modeloreparameterizado MBLPN.

Dada presença de ortogonalidade entre os subconjuntos dos parâmetrospara os modelos mistos, estimadores de máxima verossimilhaça para osparâmetros δ0 e δ1 para os modelos BIZU e MBLPN coincidem e são da-dos por δ̂0 = 0.6055(0.0066) e δ̂1 = 0.0313(0.0023). Para a parte cont́ınua,


os EMV sob o modelo BIZU são dados por µ̂ = 0.2974(0.0043) e σ̂ =0.4562(0.0050). Por outro lado, para o modelo MBLPN os EMVs são dadospor µ̂ = −0.6779(0.0419), σ̂ = 0.4289(0.00001) e α̂ = 29.8227(1.1484).

Para o caso do modelo LPNDC, temos os seguintes EMVs µ̂ = −0.8137(0.1065),η̂ = 0.5834(0.0259) e α̂ = 5.8809(1.4062). A porcentagem de zeros e uns naamostra são 0.6055 e 0.0313, respectivamente, e da função de distribuiçãoacumulada obtem-se 0.6063 e 0.0284, respectivamente, revelando bom ajustedo modelo.

EMVs para os parâmetros no modelo NDC são dados por µ̂ = −0.1556(0.0104)e σ̂ = 0.5420(0.0099), enquanto que para o modelo LNDC as EMVs são da-dos por µ̂ = −0.1375(0.0068) e σ̂ = 0.3239(0.0057). Por outro lado, parao modelo PNDC são dados por ξ̂ = −0.9895(0.1447), η̂ = 0.7394(0.0335) eα̂ = 5.2200(1.3687).

4.11 Testando modelos disjuntos

Para comparar os modelos MBLPN e LPNDC contra o modelo BIZU, umenfoque para modelos disjuntos deve ser utilizado. Sendo Fθ e Gγ doismodelos disjuntos, e f(yi|xi, θ) e g(yi|xi, β) as densidades correspondentes,a estat́ıstica da razão de verossimilhanças pode ser escrita como

LR(θ̂, β̂) ≡ `f (θ̂)− `g(β̂) =n∑i=1

logf(yi|xi, θ̂)g(yi|xi, β̂)

,

que não segue distribuição quiquadrado em grandes amostras.

Consideramos a proposta de Vuong (1989) baseada na divergência deKullback-Leibler (Kullback e Leibler, 1951). Baseando-se na distância entrecada modelo e o verdadeiro processo gerando os dados, ou seja, h0(yi, Xi),temos a estat́ıstica

TLR,NN =1√n

LR(θ̂, β̂)

ω̂2,

onde

ω̂2 =1

n

n∑i=1

(log

f(yi|xi, θ̂)g(yi|xi, β̂)

)2−

(1

n

n∑i=1

(log

f(yi|xi, θ̂)g(yi|xi, β̂)

))2

4.12 CONCLUSÕES 41

é um estimator para a variância de 1√nLR(θ̂, β̂).

Mostra-se que, quando n→∞,

TLR,NNd→ N(0, 1)

sob

H0 : E

[log

f(yi|xi, θ)g(yi|xi, β)

]= 0,

isto é, os modelos são equivalentes. Ao ńıvel de 5%, sendo z0.025 o valorcŕıtico, rejeitamos a equivalência se TLR,NN > z0.025, (ou se TLR,NN <−z0.025).

Para os dados em estudo, sendo Fθ a fda do modelo LPNDC e Gβ, domodelo BIZU, o enfoque de Vuong leva ao valor observado TLR,NN = 21.8608que é maior que o valor cŕıtico z0.025 = 1.96 de modo que BIZU é o melhordos dois modelos.

De maneira similar, comparando os modelos MBLPN e BIZU, temosque TLR,NN = −19.4777, favorecendo o modelo MBLPN levando então aconclusão de que o modelo MBPLN produz melhor ajuste para os dados emquestão.

4.12 Conclusões

Discutimos uma alternativa para a regressão beta para a situação infla-cionada de zeros e uns. O enfoque é baseado em uma extensão do mod-elo tobit com excesso de zeros que está desenvolvida em Moulton e Halsey(1995). Parâmetros são estimados por MV e a matriz de informação ob-servada (Hessiana) é usada para estimar variâncias assintóticas. Aplicaçãoa dados reais indica melhor desempenho do modelo proposto MBPLN, su-perando o mo-delo BIZU.

Caṕıtulo 5

Modelos bimodaiscensurados

Em estudos antiretrovirais de HIV, a concentração viral tem limite de deteção(mı́nimo) podendo ser 20 ou de 50 copias por miliĺıtro. O HIV-1 RNAtem tipicamente dois valores modais correspondendo as concentrações viraisótimas e subotimas, respectivamente. Os modelos podem ser vistos como ex-tensões diretas do modelo tobit censurado adequados para o ajuste de dadosunimodais e bimodais simétricos e assimétricos. Assim, os modelos esten-dem o modelo tobit usual para situações bimodais simétricas e assimétricas.EMV é implementada e MIF é derivada para tais modelos. Applicações a da-dos reais são implementadas ilustrando a performance bastante satisfatóriados modelos considerados.

O problema da concentração de HIV RNA em amostras de sangue (escalalog10) de pacientes com HIV apresenta limite de deteção mı́nimo como noproblema da vacinação no Haiti; para o teste Roche Amplicor este limite éda 50 copias/ml.

Este caṕıtulo está direcionado para uma extensão do modelo tobit paramodelos simétricos e assimétricos bimodais. No estudo de Li et al. (2006),conclui-se que a distribuição do HIV RNA (log10) é bimodal, a qual con-sideram ser uma mistura de duas distribuições normais, refletindo respostasdiferentes para terapias antiretrovirais (HAART). Como trabalhar com mis-turas de distribuições apresenta dificuldades (falta de identificabilidade, porexemplo) (Marin et al., 2005), consideramos um caminho alternativo quesegue da extensão dos modelos normais-assimétricos e potência-normal. Faze-

43

44 MODELOS BIMODAIS CENSURADOS 5.2

mos uso de MV para estimação dos parâmetros. Julgamos ser fact́ıvel o usode inferência Bayesiana.

Seção 6.2 apresenta revisão básica de modelos bimodais simétricos eassimétricos. A Seção 6.3 é direcionada a uma extensão do modelo nor-mal usual para dados censurados (modelos tipo tobit) podendo incorporarsituações uni e bimodais. Estimação é considerada por MV e por métodosBayesianos. Seção 6.6 trata de uma aplicação a um conjunto de dados deuma cĺınica na Colômbia.

5.1 Modelos assimétricos bimodais

Como visto anteriormente, Azzalini (1985) considera a seguinte representaçãogeral para uma distribuição assimétrica:

ϕ(z;λ) = 2f(z)G(λz), z, λ ∈ R,

onde f é uma fdp simétrica em torno de zero e G é fda simétrica e absolu-tamente continua e λ é o parâmetro de assimetria. Mais resultados podemser vistos em Azzalini (1986), Henze (1986), Chiogna (1997) e Pewsey (2000).

Em particular, se f = φ e G = Φ, a fdp e fda da N(0,1), obtemos

ϕ(z;λ) = 2φ(z){Φ(λz)}, z ∈ R,

que denotamos por Z ∼ SN(λ).

5.2 Extensões bimodais para modelos simétricos

Uma modificação para tornar o modelo normal assimétrico bimodal, apareceem Kim (2005),

f(z;λ) = cλφ(z)Φ(λ|z|), z ∈ R,

onde cλ é a constante de normalização, que não é simples de ser obtida.Kim (2005) mostra que este modelo produz densidades simétricas. Umaversão assimétrica do modelo de Kim aparece em Gomez et al. (2009), queconsidera

f(z;λ) = cλφ(z)Φ(λ|z|)Φ(βz), z ∈ R,

5.2 EXTENSÕES BIMODAIS PARA MODELOS SIMÉTRICOS 45

onde cλ é a constante de normalização. Dada a dificuldade de se trabalharcom o modelo acima devido a dificuldade de ser tarbalhar com a constantede normalização, Martinez et al. (2012b) propõe uma modificação bimodal(simétrica) no modelo potência-normal (PN) (Pewsey et al., 2012), con-siderando

f(z|α) = αcαφ(z){Φ(|z|)}α−1,

α > 0, com

cα =2α−1

2α − 1.

Extensão para o caso locação-escala segue fazendo X = ψ + ηZ.Note que neste caso a constante de normalização é bastante simples.

A matriz de informação de Fisher para localização-escala é dada por

IF =

1/η2 0 a01/η2/η2 a11/η(1 + 2(log2)2.

Pode-se mostrar que

|IF | = 2.808/η4.

Para tornar o modelo bimodal assimétrico usamos o enfoque em Gomezet al. (2009), que leva a fdp (Martinez et al., 2012b)

f(z|α, β) = 2αcαφ(z){Φ(|z|)}α−1Φ(βz),

α > 0, z ∈ R, with

cα =2α−1

2α − 1.

A extensão locação-escala segue tomando X = ψ+ηZ. Maximização daverossimilhança deve ser feita numericamente.

A matriz de informação de Fisher para o modelo de locação-escala é dadapor

IF =

1/η2 0

√2π/η a01/η

2/η2 0 a11/η2/π 00 (1 + 2(log2)2)

Pode-se mostra que

|IF | = −0.2999/η4 6= 0.


Pode-se testar normalidade, i.e., H0 : α = 1.0, β = 0, usando a estatisticada razão de verossimilhnaça.

5.2.1 Aplicação: Dados de poluição.

Apresentamos a seguir o ajuste dos modelos acima a um conjunto de dadosreais relacionados com (Y :) poluição nos EUA. O conjunto de dados éapresentado a seguir.

67,54.7,7.0,48.5,14,17.2,20.7,13,43.4,40.2,38.9,54.5,59.8,48.3,

22.9,11.5,34.4,35.1,38.7,30.8,30.6,43.1,56.8,40.8,41.8,42.5,31.0,31.7,

30.2,25.9,49.2,37,35.9,15,30.2,7.2,36.2,45.5,7.8,33.4,36.1,40.2,

42.7,42.5,16.2,39,35,37,31.4,37.6,39.9,36.2,42.8,46.4,24.7,49.1,

46,35.9,7.8,48.2,15.2,32.5,44.7,42.2,38.8,17.4,40.8,29.1,14.6,59.2

Pode-se mostrar que ȳ = 34.9 e s2y = 187.8. Ajustando a normalN(34.9; 187.80), nota-se que não é bom o ajuste deste modelos aos dados.Nota-se também a partir do histograma que os dados apresentam bimodal-idade, de modo que um modelo assimétrico apresentaria uma juste melhoraos dados acima.

Ajustamos então no WinBugs o modelo

f(x|µ, σ, α, beta) ∝ 2αcαφ(z){Φ(|z|)}α−1Φ(βz),

α > 0, z ∈ R, com z = (x− µ)/sigma).

Temos então a notação (µ, σ, α, β) = (mu, sig, lb, beta), com o código

z[i] < −(y[i]−mu)/sig

logLike[i] < −(−log(sig)) + log(lb) + (lb− 1) ∗ log(2)

−log(pow(2, lb)− 1)− (pow(z[i], 2)/2) + (lb− 1) ∗ log(phi(abs(y[i])))

+log(phi(beta ∗ z[i]))

que apresenta as estimativas

µ̂ = 22, σ̂ = 14, α̂ = 4.5, β̂ = 1.0.

5.3 MODELO FLEXÍVEL NORMAL CENSURADO 47

Figura 5.1: Densidade estimada e histograma dos dados.

Veja os gráficos da fda acima para os valores estimados sobre o his-tograma dos dados. Existe indicação de melhor ajuste do modelo bimodal.

5.3 Modelo flex́ıvel normal censurado

Nesta seção estendemos o modelo tobit usual para a situação normal bi-modal. Tomando λ = 0 em Gomez et al. (2009), obtemos a fdp

f(y;λ) = cδφ(|y|+ δ),

onde δ é um número real e cδ = (2(1 − Φ(δ)))−1 é a constante de normal-ização. De maneira similar ao modelo acima, este model é bimodal para δmenor que zero. Denominamos este modelo normal flex́ıvel e denotamos porFN(δ).

Considere agora que y∗ denota a distribuição FN(δ) e que (y∗1, y∗2, ..., y

∗n) é

uma amostra de uma variável aleatória onde somente valores y∗ maiores quea constante c são observados. Para valores y∗ ≤ c somente o valor c é reg-istrado. Deste modo, os valores observados são dados por

yi =

{y∗i , se y

∗i > c

c, se y∗i ≤ c,


i = 1, 2, ..., n.

A amostra resultante é censurada à esquerda. Neste caso dizemos que avariável aleatória Y tem distribuição censurada normal flex́ıvel e denotamospor CNF (δ). A distribuição desta variável aleatória é bimodal para valoresde δ menores que zero e unimodal para valores de δ maiores que zero. Paraδ = 0 temos o modelo normal usual.

5.3.1 Momentos

Os momentos de Z ∼ CFN(δ) são funções dos momentos da distribuiçãonormal, e são dados por

µr(a) =

∫ ∞a

zrφ(z)dz.

O r-ésimo momento da variável aleatória Z ∼ CFN(δ) são dados por

E(Zr) = µr = cδ

r∑k=0

(r

k

)(−δ)r−kµk(c+ δ).

Para c = 0, segue que a esperança e variância da variável aleatória Z sãodadas por

µ = cδ[φ(δ)− δ(1−Φ(δ))] e σ2 = µ2−µ2 = c2δ [2(1−φ(δ))2−φ2(δ)].

5.3.2 Extensão para localização-escala

Para o modelo normal com média µ e variância σ2, dizemos que a variável Xsegue a distribuição flex́ıvel normal de localização-escala se sua fda é dadapor

f(x;λ) =cδσφ

(∣∣∣∣x− µσ∣∣∣∣+ δ) , x ∈ R,

com µ > 0 e σ parâmetros de localização e escala. Assim, definindo

yi =

{xi, se xi > cc, se xi ≤ c,


obtemos a distribuição normal flex́ıvel, que denotamos por NCF (µ, σ, δ).

Também, o r-ésimo momento da variável Y ∼ CNF (µ, σ, δ) é dado por:

E[Y r] = µr = cδ

r∑k=0

(r

k

)δr−k

[µk

(−µ+ σδ

σ,−δ

)+ (−1)r−kµk(δ)

],

onde µr(a, b) =∫ ba z

rφ(z)dz.

5.3.3 Estimação

Denotamos por∑

0 a soma para as observações censuradas e∑

1 a soma paraas observações não censuradas. Assim, para observações com yi = 0 temosque

P [yi = 0] = P [xi ≤ 0] = cδ{

1− Φ(µ+ σδ

σ

)}e para yi > 0, a distribuição de yi é igual a distribuição de xi, isto éyi ∼ NF (µ, σ, δ).

Para uma amostra de n unidades, y1, y2, ..., yn, a função de log-verossimilhançapara θ = (µ, σ, δ)′ é dada por

`(θ; X) =∑

0

log

[cδ

(1− Φ

(µ+ σδ

σ

))]+∑

1

[log(cδ)− log(σ) + log(φ(|zi|+ δ))] ,

onde zi =yi−µσ , i = 1, ..., n.

Temos então o escore

U(µ) = −n0σ

φ(µ+σδσ

)1− Φ

(µ+σδσ

) + 1σ

∑1

yi − µσ− δσ

∑1

sgn(yi − µ),

U(σ) =n0µ

σ

φ(µ+σδσ

)1− Φ

(µ+σδσ

) − n1σ

+1

σ

∑1

(yi − µσ

)2+δ

σ

∑1

∣∣∣∣yi − µσ∣∣∣∣,


U(δ) = −n0φ(µ+σδσ

)1− Φ

(µ+σδσ

) + nφ(δ)1− Φ(δ)

−∑

1

∣∣∣∣yi − µσ∣∣∣∣− n1δ,

onde n0 e n1 como acima denotam o número de observações censuradas enão censuradas, respectivamente. Igualando escore a zero obtem-se sistemade equaç oes (com solução iterativa) que leva aos EMV. A função ”optim”doR pode ser empregada.

5.3.4 Matriz de informação

Nesta subseção apresentamos as matrizes de informação esperadas e ob-servadas para o modelo NFC(µ, σ, δ). Iniciamos com a matriz Hessiana,a saber, a segunda derivada da função log-verossimilhança com respeitoaos parâmetros do modelo (multiplicada por (-1)), para as quais usamos anotação jµµ, jηµ, jδµ jηη, jδσ e jδδ, levando as seguintes expressões:

jµµ =n1σ2

+n0σ2

φ(µ+σδσ

)1− Φ

(µ+σδσ

)−µ+ σδ

σ+

φ(µ+σδσ

)1− Φ

(µ+σδσ

) ,

jηµ =n0σ2

(µ

σ(µ+ σδ

σ)− 1)

φ(µ+σδσ )

1− Φ(µ+σδσ )− n0σ2µ

σ[

φ(µ+σδσ )

1− Φ(µ+σδσ )]2

+2

σ2

∑1

yi − µσ

− δ+σ2∑

1

sgn(yi − µ),

jηη =n0µ

σ2

(1− µ

σ

(µ+ σδ

σ

)) φ(µ+σδσ )1− Φ

(µ+σδσ

) + n0σ2µ

σ

φ(µ+σδσ

)1− Φ

(µ+σδσ

)2

−n1σ2

+3

σ2

∑1

(yi − µσ

)2+

2δ

σ2

∑1

∣∣∣∣yi − µσ∣∣∣∣ ,


jδµ = −n0σ

µ+ σδσ

−φ(µ+σδσ

)1− Φ

(µ+σδσ

) φ

(µ+σδσ

)1− Φ

(µ+σδσ

) + 1σ

∑1

sgn(yi − µ),

jδσ =n0µ

σ

µ+ σδσ

−φ(µ+σδσ

)1− Φ

(µ+σδσ

) φ

(µ+σδσ

)1− Φ

(µ+σδσ

) − 1σ

∑1

∣∣∣∣yi − µσ∣∣∣∣ ,

jδδ = −n0

µ+ σδσ

−φ(µ+σδσ

)1− Φ

(µ+σδσ

) φ

(µ+σδσ

)1− Φ

(µ+σδσ

)+n

(δ − φ(δ)

1− Φ(δ)

)φ(δ)

1− Φ(δ)+ n1.

Para obter a matriz de informação observada avaliamos os elementos daHessiana acima nos EMVs. Para obter MIF calculamos os valores esperadosdos elementos da Hessiana acima, usando a notação iµµ, iηµ, iδµ iηη, iδσ eiδδ, conforme pode ser visto em Martinez et al. (2012b).

iθrθp = n−1E

{−∂

2`(θ; x)

∂θr∂θp

}, r, p = 1, 2, 3,

com θ1 = µ, θ2 = σ e θ3 = δ com:

iµµ =1

σ2

[1− cδ

(1− Φ

(µ+ σδ

σ

))]+cδσ2φ

(µ+ σδ

σ

)−µ+ σδσ

+φ(µ+σδσ

)1− Φ

(µ+σδσ

) ,

iηµ =cδσ2φ

(µ+ σδ

σ

)µσ

µ+ σδσ

−φ(µ+σδσ

)1− Φ

(µ+σδσ

)− 1

− δcδσ2

(1− Φ

(µ+ σδ

σ

))

+2cδσ2

[φ

(µ+ σδ

σ

)+ φ(δ) + δ

(Φ

(µ+ σδ

σ

)+ Φ(δ)− 3

2

)− 1√

2π

],


iηη =µcδσ2

φ

(µ+ σδ

σ

)1 + µσ

−µ+ σδσ

+φ(µ+σδσ

)1− Φ

(µ+σδσ

)− 1

σ2+

cδσ2

[−2δφ(δ) + (1 + 2δ2)

(1− Φ

(µ+ σδ

σ

))− 4δ2(1− Φ(δ))

]+

cδσ2

[(3

(µ− σδσ

)+ 2δ

)φ

(µ+ σδ

σ

)+ 3(1 + δ2)

(1− 2Φ(δ) + Φ

(µ+ σδ

σ

))],

iδµ =cδσφ

(µ+ σδ

σ

)−µ+ σδσ

+φ(µ+σδσ

)1− Φ

(µ+σδσ

)+ cδ

σ

(1− Φ

(µ+ σδ

σ

)),

iδσ =cδµ

σφ

(µ+ σδ

σ

)µ+ σδσ

−φ(µ+σδσ

)1− Φ

(µ+σδσ

)− δcδ

σ

(1− Φ

(µ+ σδ

σ

))

+cδσ

[2δ (1− Φ(δ))− 2φ(δ) + φ

(µ+ σδ

σ

)],

iδδ = cδφ

(µ+ σδ

σ

)−µ+ σδσ

+φ(µ+σδσ

)1− Φ

(µ+σδσ

)+ φ(δ)

1− Φ(δ)

[δ − φ(δ)

1− Φ(δ)

]

+ 1− cδ(

1− Φ(µ+ σδ

σ

)).

Mostramos que a MIF acima não é singular, de modo que o resultadoseguinte segue das condições de regularidade usuais.

Teorema 6.1. Se θ̂ é o EMV de θ, então

θ̂A→ N3(θ, IF (θ)−1),

de modo que a matriz de covariâncias assintóticas do EMV θ̂ é a matrizinversa da MIF I(θ) a qual denotamos por Σθ = I(θ)

−1.

Segue do teorema que podemos testar normalidade (H0 : δ = 0) usandoa estat́ıstica da razão de verossimilhanças. Tal resultado não vale, por ex-emplo, para o modelo em Arnold et al. (2009) para o qual a MIF é singular.

5.4 O MODELO BIMODAL SIMÉTRICO NORMAL CENSURADO 53

5.4 O modelo bimodal simétrico normal censurado

O modelo proposto por Kim (2005),

f(z;λ) = cλφ(z)Φ(λ|z|),

onde λ é um número real,

cλ = 2π/(π + 2arctan(λ))

é a constante de normalização, é uma alternativa viável para o ajuste dedados bimodais simétricos, com λ > 0. Usamos a notação TN(λ).

Pode-se estender o modelo para a situação onde parte das observaçõessão censuradas, considerando Z ∼ TN(λ), onde

yi =

{zi, se zi > cc, se zi ≤ c,

que denotamos por CTN(λ). Assim, para λ > 0 temos o modelo bimodalsimétrico.

A fdp para a variavel Y, truncada a direita, é dada por

f(y|y > c) = 2cλφ(y)Φ(λ|y|)1 + cλ[Φ(c)− 0.5 + π−1 arctan(λ)− 2T (c, λ)]

,

onde T (., λ) é a função de Owen (1956).Os momentos da variável aleatória Y podem ser obtidos a partir dos

momentos da variável aleatória com densidade acima, levando aos seguintesmomentos marginais:

E[Y ] = µ =cλ

2√

2π

[λ√

1 + λ2+ 1

],

E[Y 2] = cλ

[1

4+

1

2πarctanλ+

1

2π

λ√1 + λ2

]e

E[Y 3] =cλ

2√

2π

[2 +

3λ+ 2λ3

(1 + λ2)3/2

].


Temos também

E[Y 4] = cλ

[3

4+

3

2πarctanλ+

1

2π

λ(2λ2 + 5)

(1 + λ2)2

].

Temos então que a variância da variável Y é dada por

σ2 =cλ

4π(π + 2 arctanλ)((π + 2 arctanλ)2

+4λ√

1 + λ2(π + arctanλ)− π

(2λ2 + 1

1 + λ2

)).

5.4.1 Estimação por máxima verossimilhança

A extensão localização-escala para Kim (2005) pode ser escrita como

f(x;µ, σ, λ) =cλσφ

(x− µσ

)Φ

(λ

∣∣∣∣x− µσ∣∣∣∣)

onde cλ = 2π/(π+ 2 arctan(λ)) é a constante de normalização. Sendo∑

0 e∑1 como nas seções anteriores, a função de log-verosssimilhança é dada por

`(θ; Y) =∑

0

log

[1

2

(1− cλ

[Φ(µσ

)− 0.5 + π−1 arctan(λ)− 2T

(µσ, λ)])]

+

∑1

[log(cλ)− log(σ) + log(φ(zi)) + log(Φ(λ|zi|))] ,

onde zi =yi−µσ . Assim, os elementos da função escore são dados por

U(µ) = −2n0cλσ∆

φ(µσ

)Φ

(λµ

σ

)+

1

σ

∑1

yi − µσ

+λ

σ

∑1

sgn(yi − µ)φ(yi−µ

σ

)Φ(∣∣yi−µ

σ

∣∣) ,U(σ) =

2n0µcλσ2∆

φ(µσ

)Φ

(λµ

σ

)− n1

σ

+1

σ

∑1

(yi − µσ

)2− λσ

∑1

yi − µσ

φ(yi−µ

σ

)Φ(∣∣yi−µ

σ

∣∣) ,

5.5 MODELO BIMODAL NORMAL-ASSIMÉTRICO 55

U(λ) = − ncλπ(1 + λ2)

+2n0cλ

(1 + λ2)∆φ(µσ

)φ

(λµ

σ

)

+∑

1

∣∣∣∣yi − µσ∣∣∣∣ φ

(yi−µσ

)Φ(∣∣yi−µ

σ

∣∣) ,onde

∆ = 1− cλ[Φ(µσ

)− 0.5 + π−1 arctan(λ)− 2T

(µσ, λ)],

onde n0 e n1 são como acima. Soluções para as equações obtidas igualandoos escores acima a zero devem ser resolvidas numericamente.

Os elementos da matriz Hessiana são dados em Martinez et al. (2012b).Esta matriz também pode ser obtida diretamente do R quando se usa arotina ”optim”.

5.4.2 Matriz de informação esperada

A matriz de informação esperada (MIF) pode ser calculada a partir damatriz de informação observada tomando esperança para cada um de seuselementos, a saber

Iθrθp = E

{−∂

2`(θ; x)

∂θr∂θp

}, r, p = 1, 2, 3,

con θ1 = µ, θ2 = σ e θ3 = λ. Esta matriz é apresentada em Martinez et al.(2012b).

5.5 Modelo bimodal normal-assimétrico

Como mencionado na seção anterior, o modelo bimodal lá apresentado ajustamodelos simétricos. Não é, portanto, adequado para situações onde os dadossão assimétricos. Para tais situações, propomos usar o modelo propostoem Arnold et al. (2009), que denotamos ETN(λ, β), de modo que para asituação localização-escala, temos que X ∼ ETN(µ, σ, λ, β). Considerandoa situação censurada, onde


yi =

{xi, se xi > cc, se xi ≤ c,

Usamos a notação CETN(µ, σ, λ, β). Então, para c = 0, a contribuiçãopara a verossimilhança de observações menores ou iguais a zero é dada por

Ψ(0) = P [y = 0

Documents

MODELOS DE REGRESSAO COM RESPOSTAS~ PARCIAIS Heleno Bolfarine Jorge Bazan - SBM · 2014. 10. 20. · esta ultima sendo conhecida como raz~ao de Mills. Os resultados para modelos de