Estudo da prevalência da tuberculose: uso de métodos bayesianos

Study of the prevalence of tuberculosis using Bayesian methods

Estudo da prevalência da tuberculose: uso de métodos bayesianos

Study of the prevalence of tuberculosis using Bayesian methods

Jorge Alberto Achcar; Antonio Ruffino Netto

Faculdade de Medicina de Ribeirão Preto, Universidade de São Paulo. Ribeirão Preto, SP CEP 14049-900

RESUMO

Neste artigo, apresentamos estimadores bayesianos para a prevalência de tuberculose usando métodos computacionais de simulação de amostras da distribuição a posteriori de interesse. Em especial, consideramos o uso do amostrador de Gibbs para simular amostras da distribuição a posteriori, e daí encontramos, em uma forma simples, inferências precisas para a prevalência de tuberculose. Em uma aplicação, analisamos os resultados do exame de Rx do tórax no diagnóstico da tuberculose. Com essa aplicação, verificamos que os estimadores bayesianos são simples de se obter e apresentam grande precisão. O uso de métodos computacionais para simulação de amostras como o caso do amostrador de Gibbs tem sido recentemente muito utilizado para análise bayesiana de modelos em bioestatística. Essas técnicas de simulação usando o amostrador de Gibbs são facilmente implementadas e não exigem muito conhecimento computacional, podendo ser programadas em qualquer software disponível. Além disso, essas técnicas podem ser consideradas para o estudo da prevalência de outras doenças.

Palavras-chave: Prevalência de tuberculose. Análise bayesiana. Amostrador de Gibbs.

ABSTRACT

In this paper we present Bayesian estimators of the prevalence of tuberculosis using computational methods for simulation of samples of posterior distribution of interest. We especially considered the Gibbs sampling algorithm to generate samples of posterior distribution, and from these samples we obtained accurate inferences for the prevalence of tuberculosis. In an application, we analyzed the results of lung X-ray tests in the diagnosis of tuberculosis. With this application, we verified that Bayesian estimators are more accurate than some existing estimators usually considered by health researchers. The use of computational methods for simulation of samples as the case of the Gibbs sampling algorithm is becoming very popular for Bayesian analysis in biostatistics. These simulation techniques using the Gibbs sampling algorithm are easily implemented and do not require great computational expertise and usually can be performed using available existing software. We could also consider these techniques for studying the prevalence of other diseases.

Key Words: Prevalence of tuberculosis. Bayesian analysis. Gibbs sampler.

Introdução

Na epidemiologia de qualquer doença infecto-contagiosa é sempre importante a localização e o diagnóstico das fontes de infecção para submetê-la ao tratamento adequado. No caso da tuberculose pulmonar existem vários testes para diagnosticar a doença. No estudo presente, a título de exemplo, vamos analisar os resultados do exame de Rx do tórax no diagnóstico da enfermidade citada.

Para calcular a prevalência da tuberculose em N pessoas submetidas a 3 exames radiológicos do tórax, o que seria equivalente, submetidos a um único exame radiológico lido por 3 observadores diferentes, segundo tenham apresentado 0, 1, 2 ou 3 exames (leituras) do Rx classificadas como positivas. Assim, teríamos a seguinte distribuição de freqüências (ver a Tabela 1):

Definindo-se:

D: o evento representando uma pessoa realmente com a doença (tuberculose);

p: probabilidade de uma pessoa ser realmente doente (p=P(D));

p₁: probabilidade de uma pessoa doente (D) apresentar Rx positivo (+) no primeiro exame;

p₂: probabilidade de uma pessoa não doente apresentar Rx positivo (+) no primeiro exame.

Também vamos assumir p₁ e p₂ constantes para os 3 exames (suposição de independência para o primeiro, segundo e terceiro exames).

Seja X uma variável aleatória representando o número de Rx positivos entre os 3 testes aplicados a cada paciente.

Assim, temos duas situações com as seguintes distribuições de probabilidade condicionais:

(1) Seja X o número de Rx positivos entre os doentes (D); a distribuição de probabilidade da variável aleatória X é dada por uma distribuição binomial b(3, p₁) com função de probabilidade dada por

onde x=0, 1, 2 e 3;

(2) Seja X o número de Rx positivos entre os não-doentes ; a distribuição de probabilidade da variável aleatória X é dada por uma distribuição binomial b(3, p₂) com função de probabilidade dada por

onde x=0, 1, 2 e 3;

A probabilidade não condicional da variável aleatória X (número de testes positivos) é dada por,

Como P(D) = p e P = 1-p, temos:

onde x= 0, 1, 2 , 3.

Assim, temos:

Pela fórmula de Bayes, podemos atualizar a probabilidade de um indivíduo ter a doença dado o número de testes positivos:

onde x=1, 2, 3.

Assim, observamos que

onde x = 0, 1, 2, 3.

Podemos notar que p = P(D) é a probabilidade a priori de uma pessoa ter tuberculose e P(D / X=x) é a probabilidade a posteriori de uma pessoa ter tuberculose, dado que foi observado x resultados positivos entre os 3 testes de Rx.

Como ilustração, vamos considerar os dados introduzidos por Ruffino-Neto¹ apresentados na Tabela 2, relativos às freqüências de resultados positivos de N = 10000 pessoas submetidas a 3 testes de Rx de tórax.

Como o objetivo do pesquisador de saúde publica seria a estimação dos parâmetros desconhecidos p, p₁ e p₂, usualmente se resolve o sistema de equações obtidas a partir de (5) com as freqüências relativas apresentadas na Tabela 2:

Após considerável álgebra, encontramos a solução do sistema (8) dada por

p = 0,0290; p₁ = 0,7130 e p₂ = 0,0008.

A partir desses estimadores, é possível prever o número esperado de doentes N_i P(D/X=x), onde P(D/X=x) é dado por (7).

Na tabela 3, temos os números esperados de doentes com tuberculose entre os N_i indivíduos com i testes de Rx positivos. O método das freqüências relativas usado para a obtenção dos estimadores pontuais para p, p₁ e p₂ em geral não apresenta bons estimadores (ver por exemplo, Bickel e Doksun⁸). Além disso, há necessidade da obtenção do erro-padrão associado ao estimador obtido.

É importante salientarmos que outros estimadores para p, p₁ e p₂ podem ser obtidos usando-se outros métodos de inferência estatística, entre os quais destacamos os métodos bayesianos (ver por exemplo, Press²). O uso de métodos Bayesianos via simulação MCMC (Monte Carlo em Cadeias de Markov) apresenta grande precisão e simplicidade para analisar dados médicos (ver por exemplo, Gelfand e Smith³).

Análise bayesiana do modelo

Vamos considerar X como o número de testes de Rx positivos entre os 3 testes aplicados a cada indivíduo. Neste caso, observamos que a variável aleatória X tem uma distribuição de probabilidade dada por uma mistura de duas distribuições binomiais:

onde x=0, 1, 2, 3; l₁ + l₂ = 1; p₁ = P( Rx positivo / D) e p₂ = P( Rx positivo / ).

Considerando-se uma amostra aleatória de N indivíduos, a função de verossimilhança para l₁ , p₁ e p₂ (ver,por exemplo, Mood, Graybill e Boes⁴) é dada por

Para uma análise Bayesiana do modelo, vamos considerar as seguintes distribuições a priori para l₁ , p₁ e p₂:

(i) l₁ ~ Beta(a,b), onde a e b são conhecidos;

(ii) p₁ ~ Beta(c₁,d₁), onde c₁ e d₁ são conhecidos; (11)

(iii) p₂ ~ Beta(c₂,d₂), onde c₂ e d₂ são conhecidos;

onde Beta (a,b) denota uma distribuição beta com média a/(a+b) e variância ab/[(a+b)²(a+b+1)]. Também vamos assumir independência a priori entre l₁, p₁ e p₂ .

Observamos em (11) que se a=b=1, temos uma distribuição uniforme U(0,1). Este caso deve ser considerado quando não temos informação a priori sobre os parâmetros l₁, p₁ e p₂ do modelo, isto é, assumimos uma priori não informativa para l₁, p₁ e p₂. Outras distribuições a priori também poderiam ser consideradas para l₁, p₁ e p₂. A escolha das distribuições a priori (11) foi fundamentada observando-se a variação dos parâmetros (valores entre 0 e 1).

Usando a fórmula de Bayes (ver, por exemplo, Box e Tiao⁵), encontramos a distribuição a posteriori conjunta para l₁, p₁ e p₂ dado os dados x' =(x₁, x₂ , ....,x_N ):

Para encontrarmos amostras de l₁, p₁ e p₂ da distribuição a posteriori conjunta (12), utilizamos o amostrador de Gibbs (ver, por exemplo, Gelfand e Smith³).

Para simplificação do algoritmo, introduzimos variáveis artificiais ou latentes (ver, por exemplo, Tanner e Wong⁶) para eliminarmos o produto de somatórios na distribuição a posteriori conjunta (12).

Para isso, definimos = (m_1i,m_2i), i=1,...,N onde m_1i + m_2i = 1 e m_1i tem uma distribuição de Bernoulli com probabilidade de sucesso h_1i dada por,

onde i=1,....N.

Assim, obtemos,

onde m_2i = 1 — m_1i , i=1,....,N.

Isto é,

Como m_2i + m_1i = 1, temos,

Combinando-se (12) com (16), temos:

A partir de (17), encontramos as distribuições condicionais a posteriori necessárias para o amostrador de Gibbs:

(i) l₁/p₁, p₂, x, m ~ Beta(a + r₁, b + r₂)

(ii) p₁/l_1,, p₂, x, m ~ Beta(c₁ + s₁ , d₁ + 3r₁ — s₁) (18)

(iii) p₂/l₁ , p₁, x, m ~ Beta(c₂ + s₂, d₂ +3r₂ — s₂ )

A partir das distribuições condicionais a posteriori (18), gerar amostras da distribuição a posteriori conjunta para l₁ , p₁ e p₂ seguindo o seguinte algoritmo:

(i) Dado um valor inicial l, p e p , gerar N observações m_1i da distribuição de Bernoulli com probabilidade de sucesso h_1i dada em (13);

(ii) Gerar l, p e p das distribuições condicionais a posteriori

P(l₁ / p₁ , p₂, x , m ), P(p₁ / l₁ , p₂, x , m ) e P(p₂ / l₁ , p₁, x , m ).

(iii) Repetir (i) e (ii) até conseguir uma distribuição estacionária.

Aplicação com os dados de prevalência de tuberculose

Vamos considerar os dados de prevalência da tuberculose apresentados na Tabela 2. Assumindo ignorância a priori sobre os parâmetros l₁, p₁ e p₂, considerar a=b=c₁=d₁=c₂=d₂=1 nas distribuições a priori (11) para l₁ , p₁ e p₂. Neste caso, estamos considerando distribuições a priori uniformes para os parâmetros. Na tabela 2, temos 9600 valores x_i = 0; 150 valores x_i = 1; 50 valores x_i = 2; e 50 valores x_i = 3, o que totaliza N = 10000 observações. Assim, em cada ciclo do amostrador de Gibbs, geramos 9600 valores m_1i da distribuição de Bernoulli com probabilidade de sucesso h_1i = l₁(1 - p₁)³ / [l₁(1 - p₁)³ + l₂(1 — p₂)³ ], i=1,....,9600; 150 valores m_1i da distribuição de Bernoulli com probabilidade de sucesso h_1i = l₁p₁(1 - p₁)² / [l₁p₁(1 - p₁)² + l₂p₂(1 — p₂)² ], i=9601,....,9750; 50 valores m_1i da distribuição de Bernoulli com probabilidade de sucesso h_1i = l₁p(1 - p₁) / [l₁p(1 - p₁) + l₂p(1 — p₂) ], i=9751,...,9800 e 200 valores m_1i da distribuição de Bernoulli com probabilidade de sucesso h_1i = l₁p/ [l₁p+ l₂p] para i=9801,....,10000.

Como valor inicial para o algoritmo, consideramos l= 0,03; p=0,7 e p=0,001. É importante salientar que outros valores iniciais poderiam ser considerados.

Daí, geramos l, pe pdas distribuições condicionais a posteriori dadas em (18). Na Tabela 4, temos os sumários a posteriori para os parâmetros de interesse, considerando-se 1000 amostras geradas pelo amostrador de Gibbs. Na geração das amostras, simulamos cadeias com 1200 amostras e descartamos as primeiras 200 amostras para eliminar a influência dos valores iniciais para l₁ , p₁ , p₂ . A convergência do algoritmo foi verificada a partir de métodos gráficos e de outros métodos existentes (ver, por exemplo, Gelman e Rubin⁷).

Na Tabela 4, também temos intervalos de credibilidade bayesianos com probabilidade de cobertura igual a 0,95. Outros valores para os hiperparâmetros foram considerados, obtendo-se inferências bayesianas similares para l₁, p₁, p₂. Considerando-se a=b=c₁=d₁=c₂=d₂=10, encontramos as médias a posteriori para l₁, p₁ e p₂ dadas por 0,02637; 0,90837 e 0,00530, respectivamente.

É importante salientar que a programação computacional para o amostrador de Gibbs não requer grande conhecimento de programação e pode ser feita usando-se linguagens usuais de programação. Neste exemplo, usamos o software MINITAB (Ryan e Joiner¹⁰), que não exige muito conhecimento de programação.

Observamos que os resultados Bayesianos da Tabela 4 apresentam grande precisão e não dependem de resultados assintóticos. Além disso, temos na Tabela 4 intervalos de credibilidade para os parâmetros do modelo, o que pode ser de grande interesse prático para pesquisadores de saúde publica. Esses estimadores são obtidos das 1000 amostras geradas pelo amostrador de Gibbs (estimativas de Monte Carlo das médias a posteriori).

Também podemos usar as amostras geradas de Gibbs para achar estimativas de Monte Carlo das probabilidades a posteriori P(D /X=x) dadas por (7). Essas estimativas são apresentadas na Tabela 5. As estimativas bayesianas para P(D /X=x) foram encontradas a partir de estimadores de Monte Carlo baseados nas 1000 amostras simuladas de Gibbs e dadas por:

para x=0, 1, 2 e 3.

Na Tabela 5, também temos as estimativas dos valores esperados de doentes, dado o número de Rx positivos.

Observamos que estimadores precisos foram obtidos usando-se o método bayesiano via amostrador de Gibbs. Além disso, podemos obter, de forma bem simples, intervalos de credibilidade para as probabilidades condicionais de doentes dados os resultados dos 3 testes de Rx. Intervalos de credibilidade com probabilidade de cobertura igual a 0,95 são apresentados na Tabela 6 para P(D / X=x) e para os números esperados de doentes.

Algumas notas conclusivas

O uso de métodos bayesianos pode ser de grande interesse prático para a obtenção de boas inferências na área de epidemiologia. Os métodos computacionais recentes de Monte Carlo em Cadeias de Markov com ênfase no amostrador de Gibbs levam a resultados de grande precisão e simplicidade, e modelos mais complexos, possivelmente incluindo alguns fatores de risco poderiam ser considerados sem muita dificuldade computacional adicional. Além disso, os métodos bayesianos podem incorporar opinião de especialistas da área de saúde publica, o que pode levar a melhores inferências. É importante salientar que um dos métodos clássicos de inferência mais utilizados nas aplicações é apresentado pelo método de máxima verossimilhança onde as inferências obtidas em geral são baseadas em resultados assintóticos. Observamos também que, para o modelo considerado neste artigo, podemos ter dificuldades práticas na maximização da função de verossimilhança (10), pois o modelo consiste de uma mistura de duas distribuições binomiais (ver por exemplo,Titterington et al⁹).

Recebido em: 03/07/2003

Versão final reapresentada em: 15/10/2003

Aprovação em: 12/11/2003

• 1. Ruffino-Netto A . Estudo da Prevalência da Tuberculose, Revista de Ciência e Cultura, 1977; 29: 1429-33.
• 2. Press SJ. Bayesian Statistics: Principles, Models and Applications New York: John Wiley; 1989.
• 3. Gelfand AE, Smith AFM. Sampling Based Approaches to Calculating Marginal Densities, J Am Stat Assoc 1990; 85: 398-409.
• 4. Mood AM, Graybill FA, Boes DC. Introduction to the Theory of Statistics. London: McGraw-Hill; 1974.
• 5. Box GEP, Tiao GC. Bayesian Inference in Statistical Analysis. New York: Addison-Wesley; 1973.
• 6. Tanner MA, Wong WH. The Calculation of Posterior Distributions by Data Augmentation. J Am Stat Assoc 1987; 82: 528-50.
• 7. Gelman A, Rubin DB. Inference from Iterative Simulation Using Multiple Sequences. Stat Sci, 1972; 7: 457-511.
• 8. Bickel PJ, Doksun KA. Mathematical Statistics: Basic Ideas and Selected Topics San Francisco: Holden Day; 1977.
• 9. Titterington DM, Smith AFM, Makov UE. Statistical Analysis of Finite Mixture Distributions. New York: Wiley; 1985.
• 10. Ryan B.F.; Joiner BL. Minitab Handbook, Belmont, California: Duxbury Press; 1994.

Datas de Publicação

• Publicação nesta coleção
  05 Abr 2005
• Data do Fascículo
  Dez 2003

Histórico

• Recebido
  03 Jul 2003
• Revisado
  15 Out 2003
• Aceito
  12 Nov 2003