1.1 Analisando uma resposta binária


Sim ou não. Sucesso ou fracasso. Morte ou sobrevivência. A favor ou contra. As respostas binárias podem ser o tipo mais prevalente de dados categóricos. O objetivo aqui é mostrar como estimar e fazer inferências sobre uma probabilidade de resposta binária e quantidades relacionadas. Começamos examinando uma população homogênea onde há uma probabilidade geral de ser estimada. Generalizamos depois essa situação para o cenário em que os itens amostrados vêm de um dos dois grupos.


1.1.1 As distribuições Bernoulli e binomial de probabilidade


Quase toda análise estatística começa com algum tipo de modelo estatístico. Um modelo estatístico geralmente assume a forma de uma distribuição de probabilidade que tenta quantificar a incerteza que vem com a observação de uma nova resposta. O modelo pretende representar o fenômeno desconhecido que rege o processo de observação. Ao mesmo tempo, o modelo precisa ser conveniente para trabalhar matematicamente, para que procedimentos de inferência como intervalos de confiança e testes de hipóteses possam ser desenvolvidos.

A seleção de um modelo geralmente é um compromisso entre dois objetivos concorrentes: fornecer uma aproximação mais detalhada do processo que gera os dados e fornecer procedimentos de inferência fáceis de usar. No caso de respostas binárias, o modelo natural é a distribuição de Bernoulli. Seja \(Y\) uma variável aleatória de Bernoulli com resultados de 0 e 1. Normalmente, diremos que \(Y = 1\) é um sucesso e \(Y = 0\) é um fracasso. Por exemplo, um sucesso seria uma tentativa de lance livre de basquete que é boa ou um indivíduo que é curado de uma doença por uma nova droga; uma falha seria uma tentativa de lance livre que é perdida ou um indivíduo que não é curado.

Denotamos a probabilidade de sucesso como \(P(Y = 1) = \pi\) e a probabilidade de falha correspondente como \(P(Y = 0) = 1 - \pi\). A função de probabilidade Bernoulli para \(Y\) combina essas duas expressões em uma fórmula:

\[ P(Y=y) = \pi^y (1-\pi)^{1-y} \]

para \(y = 0\) ou 1, onde usamos a convenção padrão de que uma letra maiúscula \(Y\) denota a variável aleatória e a letra minúscula \(y\) denota um possível valor de \(Y\). O valor esperado de \(Y\) é \(\mbox{E}(Y) = \pi\) e a variância de \(Y\) é \(\mbox{Var}(Y) =\pi (1 - \pi)\).

Muitas vezes, observam-se múltiplas respostas de variáveis aleatórias de Bernoulli por meio de amostragem repetida ou tentativas em ambientes idênticos. Isso leva à definição de variáveis aleatórias separadas para cada tentativa, \(Y_1, \cdots, Y_n\), onde \(n\) é o número de tentativas. Se todas as tentativas forem idênticas e independentes, podemos tratar \(W = \sum_{i=1}^n Y_i\) como uma variável aleatória binomial com função de probabilidade \[ P(W = w) = \binom{n}{w} \pi^w (1-\pi)^{1-w} \] para \(w = 0,\cdots,n\).

A função de combinação \(\binom{n}{w} = n!/\big(w!(n - w)!\big)\) conta o número de maneiras que \(w\) sucessos e \(n -w\) falhas podem ser ordenados. O valor esperado de \(W\) é \(\mbox{E}(W) = n\pi\) e a variância de \(W\) é \(\mbox{Var}(W) = n\pi (1 -\pi)\). Observe que a distribuição de Bernoulli é um caso especial da distribuição binomial quando \(n = 1\).

A seguir, mostramos como R pode ser usado para examinar propriedades da distribuição binomial.


Exemplo 1.1: distribuição Binomial no R.


O objetivo deste exemplo é calcular probabilidades simples usando uma distribuição binomial e mostrar como esses cálculos são realizados em R. Seremos muito básicos com o uso de R neste exemplo.

Considere uma variável aleatória binomial contando o número de sucessos de um experimento que é repetido \(n = 5\) vezes e suponha que há uma probabilidade de sucesso de \(\pi= 0.6\). Por exemplo, suponha que um indivíduo tenha essa taxa de sucesso em um determinado jogo de cartas ou arremessando uma bola de basquete em um gol de um local específico. Podemos calcular a probabilidade de cada número de sucessos, \(w\) = 0, 1, 2, 3, 4, 5.

Por exemplo, a probabilidade de 1 sucesso em 5 tentativas é \[ P(W=1) = \binom{5}{1} 0.6^1 (1-0.6)^{5-1} = 0.0768. \]

Este cálculo é realizado em R usando a função dbinom():

dbinom(x = 1, size = 5, prob = 0.6)
## [1] 0.0768

Dentro da função, o argumento \(x\) denota o valor observado da variável aleatória binomial, o que estamos chamando de \(w\), o argumento size é o número de tentativas \(n\) e o argumento prob é \(\pi\). Poderíamos ter usado dbinom(1, 5, 0.6) para obter a mesma probabilidade desde que os valores numéricos estejam na mesma ordem que os argumentos dentro da função, uma lista completa de argumentos e sua ordem está disponível na ajuda do dbinom(). Geralmente, sempre especificaremos os nomes dos argumentos em nosso código, exceto com as funções mais básicas.

Encontramos todas as probabilidades para \(w = 0,\cdots,5\) alterando o argumento \(x\):

dbinom (x = 0:5 , size = 5, prob = 0.6)
## [1] 0.01024 0.07680 0.23040 0.34560 0.25920 0.07776

onde 0:5 significa os inteiros de 0 a 5 por 1. Para exibir essas probabilidades em um formato mais descritivo, nós as salvamos em um objeto e imprimimos de um quadro de dados usando a função data.frame():

pmf <- dbinom(x = 0:5 , size = 5, prob = 0.6)
save <- data.frame(w = 0:5 , prob = round (x = pmf , digits = 4))
save
##   w   prob
## 1 0 0.0102
## 2 1 0.0768
## 3 2 0.2304
## 4 3 0.3456
## 5 4 0.2592
## 6 5 0.0778

A função round() arredonda os valores no objeto pmf para 4 casas decimais.

Traçamos o pmf acima usando as funções plot() e abline():

plot(x = save$w , y = save$prob , type = "h", xlab = "w", ylab = "P(W=w)", 
     main = expression(paste("Gráfico da função de probabilidade binomial para n=5, ",pi ,"=0.6")), 
     panel.first = grid(col = "gray", lty = "dotted"), lwd = 3)
abline (h = 0)

A figura acima apresenta o gráfico resultante. O valor do argumento type = “h” na função plot especifica que as barras verticais devem ser plotadas de 0 aos valores fornecidos no argumento y. O argumento main contém o título do gráfico. A função abline() traça uma linha horizontal em 0 para enfatizar a parte inferior de cada linha vertical.


Suposições


A distribuição binomial é um modelo razoável para a distribuição de sucessos em um determinado número de tentativas, desde que o processo de observação de tentativas repetidas satisfaça certas suposições. Essas suposições são:

Os próximos dois exemplos detalham essas suposições em relação às aplicações e demonstram como pode ser difícil garantir que essas suposições sejam satisfeitas.


Exemplo 1.2: Chute ao gol.


No futebol americano e canadense, os pontos podem ser marcados chutando uma bola através de uma área-alvo (gol) em cada extremidade do campo. Suponha que um experimento seja conduzido onde um placekicker tente sucessivamente cinco chutes de field goal durante o treino. Um sucesso ocorre em uma tentativa quando a bola é chutada por cima do travessão e entre as duas traves dos postes da meta. Uma falha ocorre quando a bola não atinge esse resultado (a bola é chutada para a esquerda ou direita de ambas as traves ou fica aquém do travessão). Queremos usar esses resultados para estimar a verdadeira probabilidade de sucesso do placekicker, então registramos quantos chutes são bem-sucedidos.

Para usar a distribuição binomial aqui, o experimento precisa satisfazer as seguintes condições:


Exemplo 1.3: Prevalência da doença.


A prevalência de uma doença é a proporção de uma população que sofre dessa doença. Isso é equivalente à probabilidade de que um membro da população selecionado aleatoriamente tenha a doença. Muitos estudos de saúde pública são realizados para entender a prevalência da doença, pois conhecer a prevalência é o primeiro passo para resolver os problemas sociais causados pela doença.

Por exemplo, suponha que haja uma preocupação de que uma nova doença infecciosa possa ser transmitida a indivíduos por meio de doações de sangue. Uma maneira de examinar a prevalência da doença seria coletar uma amostra de 1.000 doações de sangue e testar cada uma para a doença.

Para usar a distribuição binomial aqui, essa configuração precisa satisfazer as seguintes condições:


Os exemplos anteriores mostram que pode ser difícil satisfazer todas as suposições para um modelo binomial. No entanto, o modelo binomial ainda pode ser usado como uma aproximação do modelo verdadeiro em um determinado problema, caso em que as suposições violadas precisariam ser identificadas em quaisquer resultados declarados. Alternativamente, se as suposições não forem satisfeitas, existem outros modelos e procedimentos que podem ser usados para analisar respostas binárias. Em particular, se a probabilidade de sucesso não permanecer constante para cada tentativa - por exemplo, se a probabilidade de doença estiver relacionada a certos comportamentos de risco - podemos identificar e medir os fatores que causam as variações e, em seguida, usar um modelo de regressão para o probabilidade de sucesso.


Simulando uma amostra binomial


Como é uma amostra de uma distribuição binomial? Claro, os valores observados só podem ser 0, 1 até \(n\). A proporção de valores observados que são 0, 1 até \(n\) são regidos pela função de probabilidade e o parâmetro \(\pi\) dentro dela. A média e a variância desses valores observados também são controladas pela função de probabilidade e \(\pi\). Essas propriedades são facilmente derivadas matematicamente usando definições básicas de média e variância. Em problemas mais complexos, no entanto, as propriedades estatísticas são muito mais difíceis de derivar matematicamente.

Mostramos no próximo exemplo como simular uma amostra usando R para que possamos verificar se a teoria corresponde ao que realmente acontece. Este exemplo também apresentará a simulação computacional de Monte Carlo como uma ferramenta valiosa para avaliar um procedimento estatístico. Todos os procedimentos estatísticos têm pressupostos subjacentes ao seu quadro matemático. A simulação de Monte Carlo é especialmente útil para avaliar o desempenho dos procedimentos quando essas suposições são violadas. Por exemplo, podemos querer saber se um intervalo de confiança projetado para funcionar em amostras grandes mantém seu nível de confiança declarado quando o tamanho da amostra é pequeno.

Uma simulação de Monte Carlo funciona criando uma versão em computador da população que estamos estudando, amostrando dessa população virtual de maneira prescrita, realizando a análise estatística que estamos estudando e medindo como ela se comporta. Os detalhes de cada etapa variam de um problema para outro. Em todos os casos, extraímos um número “grande” de amostras da população virtual. Ao fazê-lo, a lei dos grandes números nos garante que o desempenho médio medido nas amostras será próximo do desempenho real do procedimento neste contexto.


Exemplo 1.4: Simulação com a distribuição binomial em R.


Abaixo está o código que simula 1.000 observações aleatórias de \(W\) a partir de uma distribuição binomial com \(\pi= 0.6\) e \(n = 5\):

set.seed(4848)
bin5 <- rbinom(n = 1000 , size = 5, prob = 0.6)
bin5[1:10]
##  [1] 3 2 4 1 3 1 3 3 3 4

A função set.seed() define um número de semente para as observações simuladas. Sem entrar nos detalhes por trás da geração de números aleatórios, a especificação do número semente nos permite obter observações simuladas idênticas cada vez que o mesmo código é executado.

A função rbinom() simula as observações, onde o argumento \(n\) fornece o número de observações e não o número de tentativas. O objeto bin5 contém 1.000 valores, onde os 10 primeiros são impressos.

A esperança populacional e a variância para \(W\) são

\[ \mbox{E}(W)=n\times \pi = 5\times 0.6 = 3 \] e \[ \mbox{Var}(W)=n\pi(1-\pi)= 5\times 0.6\times (1-0.6)=1.2. \]

Calculamos a média e a variância amostrais das observações simuladas usando as funções mean() e var():

mean(bin5)
## [1] 2.991
var(bin5)
## [1] 1.236155

A média amostral e a variância são muito semelhantes a \(\mbox{E}(W)\) e \(\mbox{Var}(W)\), como esperado. Se um número maior de observações fosse simulado, digamos 10.000, geralmente esperaríamos que essas medidas amostrais fossem ainda mais próximas de suas quantidades populacionais devido à lei dos grandes números.

Para examinar quão bem a distribuição de frequência observada segue o PMF, usamos table() para encontrar as frequências de cada resposta possível e, em seguida, usamos hist() para traçar um histograma das frequências relativas:

table(x = bin5)
## x
##   0   1   2   3   4   5 
##  12  84 215 362 244  83
hist (x = bin5 , main = expression(paste("Binomial com n=5, ",pi ," =0.6 , 1000 observações")), 
      probability = TRUE , breaks = c ( -0.5:5.5) , ylab = "Frequência relativa")
grid()
box()

Por exemplo, vemos que \(w = 3\) é observado com uma frequência relativa de \(362/1000 = 0.362\). Tínhamos encontrado anteriormente que \(P(W = 3) = 0.3456\), o que é muito semelhante à proporção observada. O histograma está na figura acima e sua forma é muito semelhante ao gráfico da função de probabilidade.

Observe que o argumento probability = TRUE fornece as frequências relativas, probability = FALSE fornece as frequências, que é o padrão e o argumento breaks especifica as classes para as barras serem entre -0.5 a 5.5 por 1, as barras não serão desenhadas corretamente aqui sem especificar quebras.


1.1.2 Inferência para a probabilidade de sucesso


O objetivo desta seção é estimar e fazer inferências sobre o parâmetro de probabilidade de sucesso da distribuição de Bernoulli. Começamos estimando o parâmetro usando sua estimativa de máxima verossimilhança, porque é relativamente fácil de calcular e possui propriedades que o tornam atraente em grandes amostras. A seguir, os intervalos de confiança para a verdadeira probabilidade de sucesso são apresentados e comparados. Muitos intervalos de confiança diferentes têm sido propostos na literatura estatística. Apresentaremos primeiro o procedimento mais simples, e depois apresentaremos várias alternativas melhores. Concluímos esta seção com testes de hipóteses para \(\pi\).


Estimador de máxima verossimilhança e inferência


A função de verossimilhança é uma função de um ou mais parâmetros condicionais aos dados observados. A função de verossimilhança para \(\pi\) quando \(y_1,\cdots,y_n\) são observações de uma distribuição de Bernoulli é \[ L(\pi| y_1,\cdots,y_n) = P(Y_1=y_1)\times \cdots \times P(Y_n=y_n) = \pi^w (1-\pi)^{n-w}\cdot \]

Alternativamente, quando registramos apenas o número de sucessos de um número de tentativas, a função de verossimilhança para \(\pi\) é simplesmente \[ L(\pi|w) = P(W=w) = \binom{n}{w}\pi^w (1-\pi)^{n-w}\cdot \]

O valor de \(\pi\) que maximiza a função de verossimilhança é considerado o valor mais plausível para o parâmetro e é chamado de estimativa de máxima verossimilhança (MLE). Pode-se mostrar que o MLE de \(\pi\) é \(\widehat{\pi}= w/n\), que é simplesmente a proporção observada de sucessos. Isso é verdade tanto para \(L(\pi|y_1,\cdots,y_n)\) quanto para \(L(\pi|w)\) porque o termo \(\binom{n}{w}\) não contém informações sobre \(\pi\).

Como \(\widehat{\pi}\) varia de amostra para amostra é uma estatística e tem uma distribuição de probabilidade correspondente. Como em todos os MLEs, \(\widehat{\pi}\) tem uma distribuição normal aproximada em amostras grandes. A média da distribuição normal é \(\pi\), e a variância é encontrada a partir de \[ \begin{array}{rcl} \widehat{\mbox{Var}}(\widehat{\pi}) & = & -\left. \mbox{E}\left( \dfrac{\partial^2 \log L(\pi|W)}{\partial \pi^2}\right)^{-1}\right|_{\pi=\widehat{\pi}} \\ & = & -\left. \mbox{E}\left( -\dfrac{W}{\pi^2}+\dfrac{n-W}{(1-\pi)^2}\right)^{-1}\right|_{\pi=\widehat{\pi}} \, = \, -\left. \mbox{E}\left( -\dfrac{n}{\pi}-\dfrac{n}{1-\pi}\right)^{-1}\right|_{\pi=\widehat{\pi}} \, = \, \dfrac{\widehat{\pi}(1-\widehat{\pi})}{n}, \end{array} \]

onde \(\log(\cdot )\) é a função log natural. Podemos escrever a distribuição como \[ \widehat{\pi} \sim N(\pi,\widehat{\mbox{Var}}(\widehat{\pi}))\cdot \]

A aproximação tende a ser melhor à medida que o tamanho da amostra aumenta.


Intervalos confidencias de Wald


Utilizando a distribuição normal, podemos considerar \((\widehat{\pi} -\pi)/\widehat{\mbox{Var}}(\widehat{\pi})^{1/2}\) como uma variável normal padrão. Então, para qualquer \(0<\alpha<1\), temos \[ P\left( Z_{\alpha/2}< \dfrac{\widehat{\pi} -\pi}{\sqrt{\widehat{\mbox{Var}}(\widehat{\pi})}} <Z_{1-\alpha/2}\right) \approx 1-\alpha, \]

onde \(Z_{\alpha}\) é o \(\alpha\)-ésimo quantil da distribuição normal padrão, ou seja, \(Z_{0.975}=1.96\). Depois de reorganizar os termos e reconhecer que \(-Z_{\alpha/2} = Z_{1-\alpha/2}\), obtemos \[ P\left( \widehat{\pi} - Z_{1-\alpha/2} \sqrt{\widehat{\mbox{Var}}(\widehat{\pi})} < \pi < \widehat{\pi} + Z_{1-\alpha/2} \sqrt{\widehat{\mbox{Var}}(\widehat{\pi})}\right) \approx 1-\alpha\cdot \]

Agora, temos uma probabilidade aproximada que tem o parâmetro \(\pi\) centrado entre duas estatísticas. Quando substituímos \(\widehat{\pi}\) e \(\widehat{\mbox{Var}}(\widehat{\pi})\) pelos valores observados da amostra, obtemos o \((1-\alpha)100\)% intervalo de confiança para \(\pi\) \[ \widehat{\pi} - Z_{1-\alpha/2} \sqrt{\widehat{\pi}(1-\widehat{\pi})/n} < \pi < \widehat{\pi} + Z_{1-\alpha/2} \sqrt{\widehat{\pi}(1-\widehat{\pi})/n}\cdot \]

Este é o intervalo usual para uma probabilidade de sucesso que é dado na maioria dos livros de estatística introdutórios. Os intervalos de confiança baseados na normalidade aproximada dos MLEs são chamados de “intervalos de confiança de Wald” porque Wald (1943) foi o primeiro a mostrar essa propriedade dos MLEs em grandes amostras.

Quando \(w\) está próximo de 0 ou \(n\), ocorrem dois problemas com este intervalo:

Discutiremos problemas adicionais com o intervalo Wald em breve.


Exemplo 1.5: Intervalos de Wald.


Suponha que \(w = 4\) sucessos sejam observados em \(n = 10\) tentativas. O intervalo de confiança de Wald de 95% para \(\pi\) é \[ 0.4 \pm 1.96 \sqrt{ 0.4(1 - 0.4)/10} = (0.0964; 0.7036), \]

onde usamos a notação abreviada entre parênteses para significar \(0.0964 < \pi < 0.7036\). O código R abaixo mostra como esses cálculos são realizados:

w <- 4
n <- 10
alpha <- 0.05
pi.hat <- w/n
var.wald <- pi.hat*(1 - pi.hat)/n
lower <- pi.hat - qnorm (p = 1- alpha /2) * sqrt ( var.wald )
upper <- pi.hat + qnorm (p = 1- alpha /2) * sqrt ( var.wald )
round ( data.frame (lower , upper ), 4)
##    lower  upper
## 1 0.0964 0.7036

No código, usamos a função qnorm() para encontrar o quantil \(1-\alpha/2\) de uma distribuição normal padrão. Podemos calcular o intervalo mais rapidamente aproveitando como R realiza cálculos vetoriais:

round(pi.hat + qnorm (p = c( alpha /2, 1- alpha /2) )*sqrt(var.wald ), 4)
## [1] 0.0964 0.7036

O intervalo de confiança é bastante amplo e pode não ser significativo para algumas aplicações. No entanto, ele fornece informações sobre um intervalo que pode ser útil em situações de teste de hipóteses. Por exemplo, um teste de \(H_0 : \pi = 0.5\) vs. \(H_1 : \pi \neq 0.5\) não rejeitaria \(H_0\) porque 0.5 está dentro desse intervalo. Se, em vez disso, o teste foi \(H_0 : \pi = 0.8\) vs. \(H_1 : \pi \neq 0.8\), há evidências para rejeitar a hipótese nula. Outras maneiras de realizar esses testes com uma estatística de teste e um \(p\)-valor serão discutidas em breve.


As inferências para \(\pi\) do intervalo de confiança de Wald dependem da aproximação da distribuição normal subjacente para o estimador de máxima verossimilhança. Para que essa aproximação funcione bem, precisamos de uma amostra grande e, infelizmente, o tamanho da amostra no último exemplo era bem pequeno. Além disso, observe que \(\widehat{\pi}\) pode assumir apenas 11 valores possíveis diferentes no último exemplo: 0/10, 1/10, \(\cdots\), 10/10, mas uma distribuição normal é uma função contínua.

Esses problemas levam o intervalo de conferência de Wald a ser aproximado, no sentido de que a probabilidade de o intervalo cobrir o parâmetro, sua cobertura ou nível de confiança verdadeiro não é necessariamente igual ao nível declarado \(1-\alpha\). A qualidade da aproximação varia com \(n\) e \(\pi\), e como veremos mais adiante, o intervalo de Wald geralmente tem cobertura < \(1-\alpha\).

Tal intervalo é chamado de intervalo liberal. Por outro lado, um intervalo com cobertura superior ao nível declarado é chamado de conservador. Embora esta última propriedade possa parecer de boa qualidade, ela pode levar a intervalos bastante amplos em comparação com outras. Queremos intervalos de confiança que coloquem o parâmetro dentro de um intervalo tão estreito quanto possível, mantendo pelo menos o nível de confiança declarado. Se quiséssemos intervalos que tivessem maior cobertura, teríamos declarado um nível de confiança maior.

Tem havido muitas pesquisas para encontrar um intervalo como esse para \(\pi\), incluindo Agresti and Caffo (2000), Agresti and Min (2001), Borkowf (2006), Brown et al. (2002), Henderson and Meyer (2001), Newcombe (2001), Suess et al. (2006) e Vos and Hudson (2005). Brown et ai. (2001) apresentam uma revisão completa da maioria dos intervalos concorrentes. Apresentamos suas recomendações a seguir, juntamente com nossos próprios pensamentos sobre os melhores intervalos.


Intervalos confidencias de Wilson


Quando \(n < 40\), Brown et al. (2001) recomendam o intervalo de Wilson ou o intervalo de Jeffreys porque mantêm níveis de confiança verdadeiros mais próximos do nível declarado do que outros intervalos. A fórmula do intervalo de Wilson é encontrada examinando a estatística de teste \[ Z_0 = \dfrac{\widehat{\pi}-\pi_0}{\sqrt{\pi_0(1-\pi_0)/n}}, \]

que é a estatística de teste escore frequentemente usada para um teste de \(H_0 :\pi = \pi_0\) vs. \(H_1 : \pi\neq \pi_0\), onde \(0 < \pi_0 < 1\). A variância no denominador de \(Z_0\) é calculada assumindo que a hipótese nula é verdadeira, em vez de usar a estimativa irrestrita baseada nos dados. Isso leva à vantagem de que o denominador não é 0 sempre que \(w = 0\) ou n.

Podemos aproximar a distribuição de \(Z_0\) com uma normal padrão para obter \[ P(-Z_{1-\alpha/2} < Z_0 < Z_{1-\alpha/2}) \approx 1-\alpha\cdot \]

Tratando a aproximação como uma igualdade, o intervalo de Wilson contém o conjunto de todos os valores possíveis de \(\pi_0\) que satisfazem a equação. Por outro lado, o conjunto de todos os valores possíveis para \(pi_0\) que levam à rejeição da hipótese nula estão fora do intervalo de confiança. O processo de formação de um intervalo a partir de um procedimento de teste de hipótese como esse é frequentemente chamado de “inversão do teste”. Como o intervalo de Wilson é baseado em um teste de pontuação, também é frequentemente chamado de intervalo de pontuação.

Os pontos finais do intervalo são encontrados definindo \(Z_0\) igual a \(\pm Z_{1-\alpha/2}\) e aplicando a fórmula quadrática para resolver \(\pi_0\). Assim, o intervalo \((1-\alpha)100\)% de Wilson é \[ \widetilde{\pi}\pm \dfrac{Z_{1-\alpha/2}\sqrt{n}}{n+Z^2_{1-\alpha/2}}\sqrt{\widehat{\pi}(1-\widehat{\pi})+\dfrac{Z^2_{1-\alpha/2}}{4n}}, \]

onde \[ \widetilde{\pi}=\dfrac{w+Z^2_{1-\alpha/2}/2}{n+Z^2_{1-\alpha/2}} \]

pode ser pensado como uma estimativa ajustada de \(\pi\).

Este intervalo recebeu o nome de Wilson (1927), que primeiro propôs encontrar um intervalo para \(\pi\) dessa maneira. Observe que o intervalo de Wilson sempre tem limites entre 0 e 1.

Os intervalos de confiança de Wald e Wilson discutidos até agora são procedimentos de inferência frequentista. A “confiança” associada a esses tipos de procedimentos de inferência ocorre pela repetição do processo de coleta de uma amostra e cálculo de um intervalo de confiança a cada vez.


Intervalos confidencias de Agresti-Coull


O intervalo de Wilson é nossa escolha preferida para um intervalo de confiança para \(\pi\). No entanto, Brown et al. (2001) recomendam o intervalo de Agresti-Coull (Agresti and Coull 1998) para \(n\geq 40\), principalmente porque é um pouco mais fácil de calcular à mão e se assemelha mais ao popular intervalo de Wald.

O intervalo \((1-\alpha)100\)% Agresti-Coull é \[ \widetilde{\pi}-Z_{1-\alpha/2}\sqrt{\dfrac{\widetilde{\pi}(1-\widetilde{\pi})}{n+Z^2_{1-\alpha/2}}} < \pi< \widetilde{\pi}+Z_{1-\alpha/2}\sqrt{\dfrac{\widetilde{\pi}(1-\widetilde{\pi})}{n+Z^2_{1-\alpha/2}}}\cdot \]

Este intervalo é essencialmente o intervalo Wald onde \(Z^2_{1-\alpha/2}/2\) sucessos e \(Z^2_{1-\alpha/2}/2\) falhas são adicionados aos dados observados. Especificamente, para \(\alpha= 0.05\), isso significa que cerca de dois sucessos e duas falhas são adicionados porque \(Z_{1-0.05/2} = 1.96 \approx 2\). Semelhante ao intervalo Wald, esse intervalo tem a propriedade indesejável que pode ter limites menores que 0 ou maiores que 1.


Exemplo 1.6: Intervalos de Wilson e Agresti-Coull.


Suponha novamente que \(w = 4\) sucessos sejam observados em \(n = 10\) tentativas. Para um intervalo de confiança de 95%, a estimativa ajustada de \(\pi\) é \[ \widetilde{\pi}=\dfrac{w+Z^2_{1-\alpha/2}/2}{n+Z^2_{1-\alpha/2}}=\dfrac{4+1.96^2/2}{10+1.96^2}=0.4278\cdot \]

Os limites dointervalo de 95% de Wilson são \[ \widetilde{\pi}\pm \dfrac{Z_{1-\alpha/2}\sqrt{n}}{n+Z^2_{1-\alpha/2}}\sqrt{\widetilde{\pi}(1-\widetilde{\pi})+\dfrac{Z^2_{1-\alpha/2}}{4n}} = 0.4278\pm \dfrac{1.96\sqrt{10}}{10+1.96^2}\sqrt{0.4278(1-0.4278)+\dfrac{1.96^2}{4\times 10}} \] levando a um intervalo de \(0.1682 < \pi < 0.6873\). Os limites do intervalo de 95% Agresti-Coull são \[ \widetilde{\pi} \pm Z_{1-\alpha/2}\sqrt{\dfrac{\widetilde{\pi}(1-\widetilde{\pi})}{n+Z^2_{1-\alpha/2}}}=0.4278\pm \sqrt{\dfrac{0.4278(1-0.4278)}{10+1.96^2}} \] levando a um intervalo de \(0.1671 < \pi < 0.6884\). Ambos os intervalos de confiança têm limites bastante semelhantes neste caso, mas são bastante diferentes dos limites do intervalo de Wald de (0.0964; 0.7036) que calculamos mais cedo.

Continuando a partir do último exemplo, segue abaixo como os cálculos são realizados em R:

p.tilde <- (w + qnorm(p = 1- alpha/2)^2/2) / (n + qnorm(p = 1- alpha/2)^2)
p.tilde
## [1] 0.4277533

I.C. Wilson

round(p.tilde + qnorm(p = c( alpha/2, 1- alpha/2) ) * sqrt(n)/(n + qnorm(p = 1- alpha/2)^2) 
      * sqrt(p.tilde*(1 - p.tilde) + qnorm (p = 1- alpha/2)^2/(4*n)), 4)
## [1] 0.1663 0.6892

I.C. Agresti - Coull

var.ac <- p.tilde*(1 -p.tilde ) / (n + qnorm(p = 1- alpha/2)^2)
round(p.tilde + qnorm(p = c( alpha/2, 1- alpha/2) ) * sqrt(var.ac), 4)
## [1] 0.1671 0.6884

Após calcular \(\widehat{\pi}\), calculamos os intervalos de Wilson e Agresti-Coull através de uma linha de código para cada um. Observe que a execução de parte de uma linha de código pode ajudar a destacar como ela funciona. Por exemplo, pode-se executar qnorm(p = c(alpha/2, 1-alpha/2)) para ver que ele calcula -1.96 e 1.96.

A função binom.confint() do pacote binom pode ser usada para simplificar os cálculos. Observe que este pacote não está na instalação padrão do R, portanto, ele precisa ser instalado antes de seu uso. Abaixo está o nosso uso da função:

library( package = binom )
binom.confint(x = w, n = n, conf.level = 1-alpha , methods = "all")
##           method x  n      mean      lower     upper
## 1  agresti-coull 4 10 0.4000000 0.16711063 0.6883959
## 2     asymptotic 4 10 0.4000000 0.09636369 0.7036363
## 3          bayes 4 10 0.4090909 0.14256735 0.6838697
## 4        cloglog 4 10 0.4000000 0.12269317 0.6702046
## 5          exact 4 10 0.4000000 0.12155226 0.7376219
## 6          logit 4 10 0.4000000 0.15834201 0.7025951
## 7         probit 4 10 0.4000000 0.14933907 0.7028372
## 8        profile 4 10 0.4000000 0.14570633 0.6999845
## 9            lrt 4 10 0.4000000 0.14564246 0.7000216
## 10     prop.test 4 10 0.4000000 0.13693056 0.7263303
## 11        wilson 4 10 0.4000000 0.16818033 0.6873262

A função calcula 11 intervalos diferentes para quando o argumento methods = “all” é usado. O primeiro, segundo e décimo primeiro intervalos são os intervalos de Agresti-Coull, Wald e Wilson, respectivamente. Consulte a ajuda da função para obter mais informações sobre os outros intervalos.


Intervalos confidencias de Clopper-Pearson


O intervalo de Clopper-Pearson (Clopper and Pearson, 1934) é o último intervalo de confiança para o qual discutiremos nesta seção. Enquanto Brown and cols. (2001) observam que o intervalo é “extremamente conservador e não é uma boa escolha para uso prático”, esse intervalo tem uma propriedade única que os outros intervalos não têm: o nível de confiança verdadeiro é sempre igual ou maior que o declarado nível. Para atingir esse nível de confiança verdadeiro, o intervalo geralmente é mais amplo do que a maioria dos outros intervalos para \(\pi\).

O intervalo usa a relação entre a distribuição binomial e a distribuição beta para atingir seu nível de confiança conservador. De fato, como a distribuição real ou exata de \(W\) é usada, o intervalo é chamado de procedimento de inferência exata. Existem muitos outros procedimentos de inferência exata disponíveis para problemas estatísticos.

Para revisar a distribuição beta, seja \(V\) uma variável aleatória beta. A função de densidade de probabilidade (PDF) para \(V\) é \[ f(v;a,b)=\dfrac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)}v^{a-1}(1-v)^{b-1}, \qquad 0<v<1, \]

onde \(a > 0\) e \(b > 0\) são parâmetros e \(\Gamma(\cdot)\) é a função gama, \[ \Gamma(c) = \int_0^\infty x^{c-1}e^{-x}\mbox{d}x, \]

para \(c > 0\).

Observe que \(\Gamma(c) = (c-1)!\) para um inteiro \(c\). Os parâmetros \(a\) e \(b\) controlam a forma da distribuição. A distribuição é assimétrica à direita para \(a > b\), e a distribuição é assimétrica à esquerda para \(a < b\). Quando \(a = b\), a distribuição é simétrica em torno de \(v = 0.5\).

O quantil \(\alpha\) de uma distribuição beta, denotado por \(v_\alpha\) ou \(\mbox{Beta}(\alpha; a, b)\), é encontrado resolvendo \[ \alpha = \int_0^{v_\alpha} \dfrac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)}v^{a-1}(1-v)^{b-1}\mbox{d}v \] para \(v_\alpha\).

O intervalo \((1-\alpha)100\)% Clopper-Pearson é simplesmente quantis de duas distribuições beta: \[ \mbox{Beta}(\alpha/2;w,n-w+1) < \pi < \mbox{Beta}(1-\alpha/2;w+1,n-w)\cdot \]

Devido à restrição \(a > 0\), o ponto final inferior não pode ser calculado se \(w = 0\). Nesse caso, o limite inferior é considerado 0. Da mesma forma, o limite superior é considerado 1 sempre que \(w = n\) devido a \(b > 0\). Como os quantis beta restantes estão estritamente entre 0 e 1, o intervalo de Clopper-Pearson respeita os limites naturais para probabilidades.

Muitas vezes, o intervalo de Clopper-Pearson é dado em termos de quantis de uma distribuição \(F\) em vez de uma distribuição beta. Isso ocorre por meio de uma relação entre as duas distribuições. Ambas as fórmulas produzem os mesmos limites e os quantis beta são fáceis de calcular, então omitimos a fórmula baseada em \(F\) aqui.

Existem algumas variações no intervalo Clopper-Pearson. O intervalo de Blaker proposto em Blaker (2000, 2001) também garante que o verdadeiro nível de confiança seja sempre igual ou superior ao nível declarado. Um benefício adicional é que o intervalo não é mais amplo que o intervalo de Clopper-Pearson e geralmente é mais estreito. Uma desvantagem é que o intervalo é mais difícil de calcular e requer um procedimento numérico iterativo para encontrar seus limites.

O programa abaixo mostra como calcular o intervalo usando a função binom.blaker.limits() do pacote BlakerCI. Outra variação do intervalo Clopper-Pearson é o intervalo mid-p. Este intervalo não garante mais que o nível de confiança verdadeiro seja maior que o nível declarado, mas será menor que o intervalo de Clopper-Pearson enquanto apresenta um desempenho relativamente bom em relação ao nível de confiança declarado (Brown et al., 2001). O programa mostra como calcular este intervalo usando a função midPci() do pacote PropCIs.


Exemplo 1.7: Intervalos de Wilson e Agresti-Coull.


Suponha que \(w = 4\) sucessos sejam observados em \(n = 10\) tentativas novamente. O intervalo de Clopper-Pearson de 95% é \(\mbox{Beta}(0.025; 4, 7) < \pi < \mbox{Beta}(0.975; 5, 6)\). A função qbeta() em R calcula esses quantis resultando em um intervalo \(0.1216 < \pi < 0.7376\). Observe que este é o mais amplo dos intervalos calculados até agora.

Abaixo está o código R usado para calcular o intervalo:

round( qbeta(p = c( alpha/2, 1- alpha/2) , shape1 = c(w, w+1), shape2 = c(n-w+1, n-w)) ,4)
## [1] 0.1216 0.7376
binom.confint(x = w, n = n, conf.level = 1-alpha, methods = "exact")
##   method x  n mean     lower     upper
## 1  exact 4 10  0.4 0.1215523 0.7376219

Dentro da chamada de função qbeta(), o argumento shape1 é \(a\) e o argumento shape2 é \(b\). Usamos vetores dentro de qbeta() para encontrar os quantis. R corresponde a cada valor de vetor para produzir o equivalente da função qbeta() separada para os limites inferior e superior. A função binom.confint() é usada como uma forma alternativa de encontrar o intervalo onde method = “exact” fornece o intervalo Clopper-Pearson. Observe que esse intervalo também foi dado anteriormente quando usamos method = “all”.


Exemplo 1.8: Prevalência de hepatite C entre doadores de sangue.


As doações de sangue são rastreadas para doenças para evitar a transmissão do doador para o receptor. Para examinar a prevalência da hepatite C entre os doadores de sangue, Liu et al. (1997) concentraram-se em 1.875 doações de sangue na cidade de Xuzhou, China. Eles observaram que 42 doações apresentaram resultado positivo para o antígeno produzido pelo organismo quando infectado pelo vírus. O intervalo de Wilson de 95% é \(0.0166 < \pi < 0.0301\), onde usamos o mesmo tipo de código R dos exemplos anteriores. Assim, com 95% de confiança, a prevalência do antígeno da hepatite C na população de doadores de sangue da cidade de Xuzhou está entre 0.0166 e 0.0301.

Na prática, calcularíamos apenas um dos intervalos discutidos nesta seção. Para fins de demonstração, a Tabela a seguir exibe intervalos de confiança adicionais de 95%. Devido ao grande tamanho da amostra, vemos que os intervalos são semelhantes, sendo o intervalo de Wald o mais diferente dos demais. Os comprimentos dos intervalos também são semelhantes, com o intervalo de Clopper-Pearson sendo um pouco mais longo que os outros.

TABELA: Intervalos de confiança para prevalência de hepatite C.
Método Intervalo Comprimento
Wald (0.0157; 0.0291) 0.0134
Agresti-Coull (0.0165; 0.0302) 0.0137
Wilson (0.0166; 0.0301) 0.0135
Clopper-Pearson (0.0162; 0.0302) 0.0140

Exemplo 1.9: Intervalos confidenciasi para \(\pi\). 


w<-4  # Sum(y_i)
n<-10
alpha<-0.05
pi.hat<-w/n
# Wald C.I.
var.wald<-pi.hat*(1-pi.hat)/n
lower<-pi.hat - qnorm(p = 1-alpha/2) * sqrt(var.wald)
upper<-pi.hat + qnorm(p = 1-alpha/2) * sqrt(var.wald)
round(data.frame(lower, upper), 4)
##    lower  upper
## 1 0.0964 0.7036
# More simply:
round(pi.hat + qnorm(p = c(alpha/2, 1-alpha/2)) * sqrt(var.wald),4)
## [1] 0.0964 0.7036
# Adjusted estimate of pi
p.tilde<-(w + qnorm(p = 1-alpha/2)^2 /2) / (n + qnorm(p = 1-alpha/2)^2)
p.tilde
## [1] 0.4277533
# Wilson C.I.
round(p.tilde + qnorm(p = c(alpha/2, 1-alpha/2)) * sqrt(n) / 
        (n+qnorm(p = 1-alpha/2)^2) * sqrt(pi.hat*(1-pi.hat) + qnorm(p = 1-alpha/2)^2/(4*n)),4)
## [1] 0.1682 0.6873
# Agresti-Coull C.I.
var.ac<-p.tilde*(1-p.tilde) / (n+qnorm(p = 1-alpha/2)^2)
round(p.tilde + qnorm(p = c(alpha/2, 1-alpha/2)) * sqrt(var.ac),4)
## [1] 0.1671 0.6884
# Binom package
library(package = binom)
binom.confint(x = w, n = n, conf.level = 1-alpha, methods = "all")
##           method x  n      mean      lower     upper
## 1  agresti-coull 4 10 0.4000000 0.16711063 0.6883959
## 2     asymptotic 4 10 0.4000000 0.09636369 0.7036363
## 3          bayes 4 10 0.4090909 0.14256735 0.6838697
## 4        cloglog 4 10 0.4000000 0.12269317 0.6702046
## 5          exact 4 10 0.4000000 0.12155226 0.7376219
## 6          logit 4 10 0.4000000 0.15834201 0.7025951
## 7         probit 4 10 0.4000000 0.14933907 0.7028372
## 8        profile 4 10 0.4000000 0.14570633 0.6999845
## 9            lrt 4 10 0.4000000 0.14564246 0.7000216
## 10     prop.test 4 10 0.4000000 0.13693056 0.7263303
## 11        wilson 4 10 0.4000000 0.16818033 0.6873262
# Could just obtain one interval at a time and save into an object
# Agresti-Coull: 
save.ci<-binom.confint(x = w, n = n, conf.level = 1-alpha, methods = "ac")
save.ci
##          method x  n mean     lower     upper
## 1 agresti-coull 4 10  0.4 0.1671106 0.6883959
# Clopper-Pearson
round(qbeta(p = c(alpha/2, 1-alpha/2), shape1 = c(w, w+1), shape2 = c(n-w+1, n-w)),4)
## [1] 0.1216 0.7376
# lower<-qbeta(p = alpha/2, shape1 = w, shape2 = n-w+1)  #This is a check
# upper<-qbeta(p = 1-alpha/2, shape1 = w+1, shape2 = n-w)   #This is a check
# cat("The C-P C.I. is:", round(lower,4) , "< pi <", round(upper,4), "\n \n")
# Clopper-Pearson
binom.confint(x = w, n = n, conf.level = 1-alpha, methods = "exact")
##   method x  n mean     lower     upper
## 1  exact 4 10  0.4 0.1215523 0.7376219
########################################
# Other intervals
###################
# LR interval
    LRT.int<-binom.confint(x = w, n = n, conf.level = 1-alpha, methods = "lrt")
    LRT.int
##   method x  n mean     lower     upper
## 1    lrt 4 10  0.4 0.1456425 0.7000216
    # -2log(Lambda)
    LRT<-function(pi.0, w, n) {
      pi.hat<-w/n
      -2*(w*log(pi.0/pi.hat) + (n-w)*log((1-pi.0)/(1-pi.hat))) 
    }
    # Plot of -2log(Lambda)
    # the values of pi below the chi^2_1,1-alpha horizontal line corresponds to the interval
    curve(expr = LRT(pi.0 = x, w = w, n = n), from = 0, to = 1, col = "red", xlab = expression(pi),
      ylab = expression(-2*log(Lambda)))
    abline(h = qchisq(p = 1-alpha, df = 1), lty = "dotted")

    # Intersection points for -2log(Lambda) and the chi^2_1,1-alpha horizontal line
    LRT(pi.0 = LRT.int$lower, w = w, n = n)  # Check - should be chi^2_1,1-alpha
## [1] 3.841437
    LRT(pi.0 = LRT.int$upper, w = w, n = n)  # Check - should be chi^2_1,1-alpha
## [1] 3.841458
    # Set -2log(Lambda)
    # chi^2_1,1-alpha equal to 0 and solve for root to find lower and upper bounds for interval
    LRT2<-function(pi.0, w, n, alpha) {
      pi.hat<-w/n
      -2*(w*log(pi.0/pi.hat) + (n-w)*log((1-pi.0)/(1-pi.hat))) - qchisq(p = 1-alpha, df = 1)
    }
    # Confidence interval, differences between here and binom.confint() are very small
    uniroot(f = LRT2, lower = 0, upper = w/n, w = w, n = n, alpha = alpha)  # Lower bound
## $root
## [1] 0.1456271
## 
## $f.root
## [1] 0.0006063786
## 
## $iter
## [1] 6
## 
## $init.it
## [1] NA
## 
## $estim.prec
## [1] 6.103516e-05
    uniroot(f = LRT2, lower = w/n, upper = 1, w = w, n = n, alpha = alpha)  # Upper bound
## $root
## [1] 0.6999695
## 
## $f.root
## [1] -0.001490819
## 
## $iter
## [1] 3
## 
## $init.it
## [1] NA
## 
## $estim.prec
## [1] 0.000116483
    # Logit interval
  binom.confint(x = w, n = n, conf.level = 1-alpha, methods = "logit")
##   method x  n mean    lower     upper
## 1  logit 4 10  0.4 0.158342 0.7025951
  # Blaker interval
  library(package = BlakerCI)
  binom.blaker.limits(x = w, n = n, level = 1-alpha)
## [1] 0.1500282 0.7170653
  # Mid-p
  library(package = PropCIs)
  midPci(x = w, n = n, conf.level = 1 - alpha)
## 
## 
## 
## data:  
## 
## 95 percent confidence interval:
##  0.1426 0.7086

Testes


Quando apenas um parâmetro simples é de interesse, como aqui, geralmente preferimos intervalos de confiança em vez de testes de hipóteses, porque o intervalo fornece uma faixa de valores de parâmetros possíveis. Normalmente, podemos inferir que um valor hipotético para um parâmetro pode ser rejeitado se não estiver dentro do intervalo de confiança para o parâmetro. No entanto, existem situações em que um valor fixo conhecido de \(\pi\), digamos \(\pi_0\), é de interesse especial, levando a uma hipótese formal de teste de \(H_0 : \pi = \pi_0\) vs. \(H_1 : \pi \neq \pi_0\).

Com relação ao intervalo de Wilson, notou-se que a estatística do teste de pontuação \[ Z_0 = \dfrac{\widehat{\pi}-\pi_0}{\sqrt{\pi_0(1-\pi_0)/n}}, \]

é frequentemente utilizado nestas situações. Quando a hipótese nula é verdadeira, \(Z_0\) deve ter aproximadamente uma distribuição normal padrão, onde a aproximação geralmente é melhor para amostras maiores.

A hipótese nula é rejeitada quando um valor incomum de \(Z_0\) é observado em relação a esta distribuição, ou seja, algo menor que \(-Z_{1-\alpha/2}\) ou maior que \(Z_{1-\alpha/2}\). O \(p\)-valor é uma medida de quão extremo o valor da estatística de teste é em relação ao que é esperado quando \(H_0\) é verdadeiro. Este \(p\)-valor é calculado como \(2P(Z > |Z_0|)\) onde \(Z\) tem uma distribuição normal padrão. Observe que este teste é equivalente a rejeitar a hipótese nula quando \(\pi_0\) está fora do intervalo de Wilson. Se desejado, a função prop.test() pode ser usada para calcular \(Z_0\), \(Z^2_0\) é realmente dado e um \(p\)-valor correspondente.

Recomendamos usar o teste de pontuação ao realizar um teste para \(\pi\). No entanto, existem procedimentos de teste alternativos. Em particular, o teste de razão de verossimilhança (LRT) é uma maneira geral de realizar testes de hipóteses e pode ser usado aqui para testar. Informalmente, a estatística LRT é \[ \Lambda = \dfrac{\mbox{Máximo da função de verossimilhança sob } H_0}{\mbox{Máximo da função de verossimilhança sob } H_0 \mbox{ ou } H_1}\cdot \]

Para o teste específico de \(H_0 : \pi = \pi_0\) vs. \(H_1 : \pi \neq \pi_0\), o denominador é \(\widehat{\pi}^w (1-\widehat{\pi} )^{n-w}\), porque o valor máximo possível da função de verossimilhança ocorre quando é avaliada no MLE. O numerador é \(\pi_0^w (1-\pi_0)^{n-w}\) porque existe apenas um valor possível da função de verossimilhança se a hipótese nula for verdadeira.

A estatística transformada \(-2\log(\Lambda)\) acaba por ter uma distribuição aproximada de \(\chi^2_1\) em grandes amostras se a hipótese nula for verdadeira. Para este teste, a estatística transformada pode ser reexpressa como \[ -2\log(\Lambda)=-2\left(w\log\Big(\pi_0/\widehat{\pi}\Big)+(n-w)\log\Big((1-\pi_0)/(1-\widehat{\pi})\Big) \right)\cdot \]

Rejeitamos a hipótese nula se \(-2\log(\Lambda ) > \chi^2_{1,1-\alpha/2}2\), onde \(\chi^2_{1,1-\alpha/2}\) é o quantil \(1-\alpha/2\) de uma distribuição qui-quadrado com 1 grau de liberdade, por exemplo, \(\chi^2_{1,0.095} = 3.84\) quando \(\alpha= 0.05\). O \(p\)-valor é \(P(A > -2\log(\lambda))\) onde \(A\) tem uma distribuição \(\chi^2_1\).

Uma maneira alternativa de calcular um intervalo de confiança para \(\pi\) é inverter o LRT de maneira semelhante ao que foi feito para o intervalo de Wilson. Esse intervalo de razão de verossimilhança (LR) é calculado automaticamente pela função binom.confint() do pacote binom, onde o valor do argumento méthods = “lrt” é usado. Fornecemos código adicional que mostra como encontrar o intervalo sem esta função. O intervalo é geralmente mais difícil de calcular do que os intervalos que recomendamos nesta seção e não tem vantagens sobre eles. Os intervalos de confiança LR geralmente são usados em alguns contextos mais complicados, onde melhores intervalos não estão disponíveis. O intervalo é melhor que o intervalo de Wald na maioria dos problemas.


1.1.3 Níveis de confiança verdadeiros para intervalos de confiança


Conforme discutido, um método de intervalo de confiança pode não atingir seu nível de confiança declarado. As razões para isso são explicadas brevemente. A Figura abaixo fornece uma comparação dos níveis de confiança verdadeiros para os intervalos de Wald, Wilson, Agresti-Coull e Clopper-Pearson. Para cada gráfico, \(n\) é 40 e o nível de confiança declarado é 0.95, \(\alpha = 0.05\). O verdadeiro nível de confiança (cobertura) para cada método de intervalo é graficado em função de \(\pi\). Por exemplo, o nível de confiança verdadeiro em \(\pi= 0.157\) é 0.8760 para o Wald, 0.9507 para o Wilson, 0.9507 para o Agresti-Coull e 0.9740 para os intervalos de Clopper-Pearson, respectivamente. Obviamente, nenhum desses intervalos atinge exatamente o nível de confiança declarado de forma consistente.

alpha <- 0.05
n <- 40
w <- 0:n
pi.hat <- w/n
p.tilde <- (w + qnorm(p = 1-alpha/2)^2/2)/(n+qnorm(1-alpha/2)^2)
# Wald
var.wald<-pi.hat*(1-pi.hat)/n
lower.wald <- pi.hat - qnorm(p = 1-alpha/2)*sqrt(var.wald)
upper.wald <- pi.hat + qnorm(p = 1-alpha/2)*sqrt(var.wald)
# Agresti-Coull
lower.AC <- p.tilde - qnorm(p = 1-alpha/2) * sqrt(p.tilde*(1-p.tilde) / (n+qnorm(1-alpha/2)^2))
upper.AC <- p.tilde + qnorm(p = 1-alpha/2) * sqrt(p.tilde*(1-p.tilde) / (n+qnorm(1-alpha/2)^2))
# Wilson
lower.wilson <- p.tilde - qnorm(p = 1-alpha/2) * sqrt(n) / (n+qnorm(1-alpha/2)^2) * 
  sqrt(pi.hat*(1-pi.hat) + qnorm(1-alpha/2)^2/(4*n))
upper.wilson <- p.tilde + qnorm(p = 1-alpha/2) * sqrt(n) / (n+qnorm(1-alpha/2)^2) * 
  sqrt(pi.hat*(1-pi.hat) + qnorm(1-alpha/2)^2/(4*n))
# Clopper-Pearson - Isso é um pouco mais complicado devido aos casos y = 0 e n
lower.CP <- numeric(n+1)  # Isso inicializa um vetor para salvar os limites inferiores
upper.CP <- numeric(n+1)  # Isso inicializa um vetor para salvar os limites superiores
# y = 0
w0 <- 0  # 
lower.CP[1] <- 0
upper.CP[1] <- qbeta(p = 1-alpha/2, shape1 = w0+1, shape2 = n-w0)
# y = n
wn <- n  # 
lower.CP[n+1] <- qbeta(p = alpha/2, shape1 = wn, shape2 = n-wn+1)
upper.CP[n+1] <- 1
# y = 1, ..., n-1
w.new <- 1:(n-1)
lower.CP[2:n] <- qbeta(p = alpha/2, shape1 = w.new, shape2 = n-w.new+1)
upper.CP[2:n] <- qbeta(p = 1-alpha/2, shape1 = w.new+1, shape2 = n-w.new)
# Para todo pi
pi.seq <- seq(from = 0.001, to = 0.999, by = 0.0005)
# pi.seq <- 0.16 #Teste
# pi.seq <- seq(from = 0.1, to = 0.9, by = 0.1) #Teste
# Salvando níveis de confiança verdadeiros em uma matriz
save.true.conf <- matrix(data = NA, nrow = length(pi.seq), ncol = 5)
# 
counter<-1
# Iterando sobre cada pi em que o verdadeiro nível de confiança é calculado
for(pi in pi.seq) {
  pmf<-dbinom(x = w, size = n, prob = pi)
  # Wald
  save.wald<-ifelse(test = pi>lower.wald, yes = ifelse(test = pi<upper.wald, yes = 1, no = 0), no = 0)
  wald<-sum(save.wald*pmf)
  # Agresti-Coull
  save.AC<-ifelse(test = pi>lower.AC, yes = ifelse(test = pi<upper.AC, yes = 1, no = 0), no = 0)
  AC<-sum(save.AC*pmf)
  # Wilson
  save.wilson = 
    ifelse(test = pi>lower.wilson, yes = ifelse(test = pi<upper.wilson, yes = 1, no = 0), no = 0)
  wilson<-sum(save.wilson*pmf)
  # Clopper-Pearson
  save.CP<-ifelse(test = pi>lower.CP, yes = ifelse(test = pi<upper.CP, yes = 1, no = 0), no = 0)
  CP<-sum(save.CP*pmf)
  save.true.conf[counter,]<-c(pi, wald, AC, wilson, CP)
  counter<-counter+1
}
# Gráficos
par(mar=c(3,4,2,1), pch = 19, mfrow = c(2,2)) 
plot(x = save.true.conf[,1], y = save.true.conf[,2], main = "Wald", xlab = expression(pi),
  ylab = "Nível de confiança verdadeiro", type = "l", ylim = c(0.85,1))
abline(h = 1-alpha, lty = "dotted")
grid()
plot(x = save.true.conf[,1], y = save.true.conf[,3], main = "Agresti-Coull", xlab = expression(pi),
  ylab = "Nível de confiança verdadeiro", type = "l", ylim = c(0.85,1))
abline(h = 1-alpha, lty = "dotted")
grid()
plot(x = save.true.conf[,1], y = save.true.conf[,4], main = "Wilson", xlab = expression(pi),
  ylab = "Nível de confiança verdadeiro", type = "l", ylim = c(0.85,1))
abline(h = 1-alpha, lty = "dotted")
grid()
plot(x = save.true.conf[,1], y = save.true.conf[,5], main = "Clopper-Pearson", xlab = expression(pi),
  ylab = "Nível de confiança verdadeiro", type = "l", ylim = c(0.85,1))
abline(h = 1-alpha, lty = "dotted")
grid()

# Pi = 0.157
save.true.conf[save.true.conf[,1]==0.157,]
## [1] 0.1570000 0.8759050 0.9506532 0.9506532 0.9739795
# Embora AC e Wilson tenham os mesmos níveis de confiança verdadeiros em pi = 0.157, 
# isso nem sempre será o caso
sum(save.true.conf[,3] != save.true.conf[,4])  # Número de diferenças
## [1] 472
length(pi.seq)  # Número de níveis de confiança verdadeiros
## [1] 1997


Abaixo estão algumas conclusões gerais do exame destes gráficos:

Achados semelhantes podem ser mostrados para outros valores de \(n\) e \(\alpha\). Por que esses gráficos na figura acima têm padrões tão estranhos? É tudo por causa da discrição de uma variável aleatória binomial. Para um dado \(n\), existem apenas \(n + 1\) intervalos possíveis que podem ser formados, um para cada valor de \(w = 0, 1,\cdots, n\). Para um valor específico de \(\pi\), alguns desses intervalos contêm \(\pi\) e outros não. Assim, o verdadeiro nível de confiança em \(\pi\), digamos \(C(\pi)\), é a soma das probabilidades binomiais para todos os intervalos que contêm \(\pi\): \[ C(\pi)=\sum_{w=0}^n \pmb{I}(w)\binom{n}{w}\pi^w (1-\pi)^{n-w}, \]

onde \(\pmb{I}(w)=1\) se o intervalo formado com \(w\) contém \(\pi\) e \(\pmb{I}(w) = 0\) se não. Cada uma dessas probabilidades binomiais muda lentamente à medida que muda \(\pi\). Contanto que não ultrapassemos \(\pi\) em nenhum limite de intervalo, o verdadeiro nível de confiança também muda lentamente. No entanto, assim que \(\pi\) ultrapassa um limite de intervalo, uma probabilidade é subitamente adicionada ou subtraída do nível de confiança verdadeiro, resultando nos picos que aparecem em todas as partes da figura acima. Ilustramos como encontrar o verdadeiro nível de confiança e quando esses picos ocorrem no próximo exemplo.


Exemplo 1.10:
Nível de confiança verdadeiro para o intervalo de Wald.


Mostramos neste exemplo como calcular um nível de confiança verdadeiro para o intervalo de Wald com \(n = 40\), \(\pi= 0.157\) e \(\alpha= 0.05\).

Segue abaixo a descrição do processo: 
pi <- 0.157
alpha <- 0.05
n <- 40
w <- 0:n
pi.hat <- w/n
pmf <- dbinom (x = w, size = n, prob = pi)
var.wald <- pi.hat *(1 - pi.hat)/n
lower <- pi.hat - qnorm (p = 1- alpha /2) * sqrt(var.wald )
upper <- pi.hat + qnorm (p = 1- alpha /2) * sqrt(var.wald )
save <- ifelse ( test = pi >lower , yes = ifelse ( test = pi <upper, yes = 1, no = 0) , no = 0)
data.frame (w, pi.hat , round ( data.frame(pmf , lower , upper ) ,4), save ) [1:13 ,]
##     w pi.hat    pmf   lower  upper save
## 1   0  0.000 0.0011  0.0000 0.0000    0
## 2   1  0.025 0.0080 -0.0234 0.0734    0
## 3   2  0.050 0.0292 -0.0175 0.1175    0
## 4   3  0.075 0.0689 -0.0066 0.1566    0
## 5   4  0.100 0.1187  0.0070 0.1930    1
## 6   5  0.125 0.1591  0.0225 0.2275    1
## 7   6  0.150 0.1729  0.0393 0.2607    1
## 8   7  0.175 0.1564  0.0572 0.2928    1
## 9   8  0.200 0.1201  0.0760 0.3240    1
## 10  9  0.225 0.0795  0.0956 0.3544    1
## 11 10  0.250 0.0459  0.1158 0.3842    1
## 12 11  0.275 0.0233  0.1366 0.4134    1
## 13 12  0.300 0.0105  0.1580 0.4420    0
sum( save*pmf )
## [1] 0.875905
sum( dbinom (x = 4:11 , size = n, prob = pi))
## [1] 0.875905


O código contém muitos dos mesmos componentes que vimos antes. A principal diferença agora é que estamos calculando um intervalo para cada valor possível de \(w\) em vez de um intervalo para apenas um. Uma nova parte dentro do código é a função ifelse().

Esta função faz uma verificação lógica para saber se está ou não dentro de cada um dos 41 intervalos. Por exemplo, o segundo intervalo é (-0.0234; 0.0734), que não contém \(\pi= 0.157\), portanto, o objeto de save tem um valor de 0 para seu segundo elemento.

O quadro de dados (data.frame) criado no final coloca todos os componentes calculados juntos em uma tabela. Por exemplo, vemos que se \(w = 3\), o intervalo correspondente não contém \(\pi\), mas em \(w = 4\) o intervalo correspondente contém. Examinando outras partes do quadro de dados, vemos que todos os intervalos para \(w = 4\) a 11 contêm \(\pi= 0.157\). A probabilidade de uma variável aleatória binomial estar entre 4 e 11 com \(n = 40\) e \(\pi= 0.157\) é 0.8759, que é o nível de confiança verdadeiro. Obviamente, o intervalo de confiança de Wald não atinge seu nível declarado de 95%.

Observe que o limite superior do intervalo em \(w = 3\) mal contém \(\pi= 0.157\) e \(P(W = 3) = 0.0689\). Usando o mesmo código com a mudança de pi <- 0.156, o limite superior em \(w = 3\) agora contém \(\pi= 0.156\), de modo que \(P(W = 3) = 0.0706\) é incluído ao somar probabilidades para o verdadeiro nível de confiança. No geral, \(w = 3\) a 11 agora têm intervalos de confiança que contêm \(\pi= 0.156\), levando a um nível de confiança verdadeiro de 0.9442. Isso demonstra o que mencionamos anteriormente como a causa dos picos na figura acima.


Em problemas simples como este, podemos determinar exatamente as probabilidades de cada intervalo que contém um dado \(\pi\), de modo que os gráficos como na Figura 1.3 possam ser feitos exatamente. Em outros casos, podemos ter que confiar na simulação de Monte Carlo. Exploramos a abordagem de simulação a seguir para nos permitir comparar um nível de confiança real exato com um estimado por simulação. Isso será útil mais adiante no texto, quando o método de simulação for o único método de avaliação disponível.


Exemplo 1.11:
Nível de confiança estimado para o intervalo


Suponha novamente que \(n = 40\), \(\pi= 0.157\) e \(\alpha= 0.05\). Abaixo está uma descrição do processo para estimar o verdadeiro nível de confiança por meio de simulação: 
numb.bin.samples <- 1000 # Binomial samples of size n
set.seed (4516)
w <- rbinom (n = numb.bin.samples, size = n, prob = pi)
pi.hat <- w/n
var.wald <- pi.hat *(1 - pi.hat)/n
lower <- pi.hat - qnorm (p = 1- alpha /2) * sqrt (var.wald )
upper <- pi.hat + qnorm (p = 1- alpha /2) * sqrt (var.wald )
data.frame(lower , upper ) [1:10 ,]
##           lower     upper
## 1   0.039344453 0.2606555
## 2   0.039344453 0.2606555
## 3   0.057249138 0.2927509
## 4   0.076040994 0.3239590
## 5   0.076040994 0.3239590
## 6   0.039344453 0.2606555
## 7   0.076040994 0.3239590
## 8  -0.006624323 0.1566243
## 9   0.022511030 0.2274890
## 10  0.007030745 0.1929693
save <- ifelse ( test = pi >lower, yes = ifelse ( test = pi < upper, yes = 1, no = 0) , no = 0)
save[1:10]
##  [1] 1 1 1 1 1 1 1 0 1 1
mean( save )
## [1] 0.878


Novamente, estamos usando muito do mesmo código do passado. A função ifelse() é usada para verificar se \(\pi\)= 0.157 está dentro de cada um dos intervalos. Por exemplo, vemos que a amostra #8 resulta em \(\widehat{\pi}\) = 0.075 e um intervalo de (-0.0066; 0.1566), que não contém 0.157, então o valor correspondente de save é 0. A média de todos os 0’s e 1 no save é 0.878. Este é o nosso nível de confiança verdadeiro estimado para o intervalo de Wald em \(n\) = 40 e \(\pi\)= 0.157.

Neste problema de simulação relativamente simples, já sabemos que os intervalos para \(w= 4,\cdots, 11\) contém \(\pi\)= 0.157 enquanto os outros não. Para ver que a simulação está, de fato, estimando \(P(4\leq W\leq 11)\), a função table() é usada para calcular o número de vezes que cada \(w\) ocorre:

counts <- table (w)
counts
## w
##   1   2   3   4   5   6   7   8   9  10  11  12  13 
##   8  35  64 123 147 165 172 123  76  46  26  11   4
sum( counts[4:11]) / numb.bin.samples
## [1] 0.878


Por exemplo, havia 64 das 1.000 observações que resultaram em um \(w = 3\). Isso é muito semelhante ao \(P(W = 3) = 0.0689\) que obtivemos para o exemplo anterior. Somando as contagens para \(w = 4,\cdots,11\) e dividindo por 1000, obtemos a mesma estimativa de 0.878 para o nível de confiança verdadeiro.

A estimativa do nível de confiança verdadeiro aqui é quase igual ao nível de confiança verdadeiro encontrado no exemplo anterior. Devido ao uso de um grande número de amostras, a lei dos grandes números garante que isso aconteça. Poderíamos até chegar a encontrar um intervalo de confiança para o verdadeiro nível de confiança! Para este caso, temos 878 “sucessos” em 1.000 “ensaios”. Um intervalo de Wilson de 95% para o próprio nível de confiança verdadeiro é (0.8563; 0.8969), que por acaso contém 0.8759, o nível de confiança verdadeiro conhecido.


A figura abaixo fornece os verdadeiros níveis de confiança para \(\pi= 0.001,\cdots, 0.999\) por 0.0005, onde usamos linhas retas para preencher os níveis de confiança ausentes entre pontos de plotagem em níveis de não usados.

Para calcular todos esses 1.997 níveis de confiança diferentes para um intervalo específico, repetimos o mesmo código de antes, mas agora para cada um \(\pi\) usando um “loop for” dentro de R. O próximo exemplo ilustra esse processo.


Exemplo 1.12: Gráfico de nível de confiança verdadeiro.


Com \(n = 40\) e \(\pi= 0.05\), calculamos os níveis de confiança verdadeiros para o intervalo de Wald usando o seguinte código:

alpha <- 0.05
n <- 40
w <- 0:n
pi.hat <- w/n
pi.seq <- seq( from = 0.001 , to = 0.999 , by = 0.0005)
# Wald
var.wald <- pi.hat*(1 - pi.hat)/n
lower.wald <- pi.hat - qnorm (p = 1- alpha /2) * sqrt (var.wald )
upper.wald <- pi.hat + qnorm (p = 1- alpha /2) * sqrt (var.wald )
# Save true confidence levels in a matrix
save.true.conf <- matrix ( data = NA , nrow = length (pi.seq), ncol = 2)
# Create counter for the loop
counter <- 1
# Loop over each pi
for(pi in pi.seq) {
  pmf <- dbinom (x = w, size = n, prob = pi)
  save.wald <- 
    ifelse ( test = pi > lower.wald, yes = ifelse ( test = pi < upper.wald, yes = 1, no = 0), no = 0)
  wald <- sum( save.wald * pmf)
  save.true.conf[ counter ,] <- c(pi , wald )
  # print ( save . true . conf [ counter ,])
  counter <- counter +1
}
plot (x = save.true.conf[,1] , y = save.true.conf[,2], main = "Wald", xlab = expression(pi), 
      ylab = "True confidence level", type = "l", ylim = c (0.85 ,1) )
abline (h = 1-alpha , lty = "dotted")
grid()


Criamos um vetor pi.seq que é uma sequência de números de 0.001 a 0.999 por 0.0005. O código de função for(pi in pi.seq), frequentemente chamado de “loop for”, instrui R a tirar um valor de pi.seq por vez. O código entre chaves, então, encontra o verdadeiro nível de confiança para \(\pi\). O objeto save.true.conf é uma matriz criada para ter 1.997 linhas e 2 colunas.

A princípio, todos os seus valores são inicializados como “NA” dentro de R. Seus valores são atualizados uma linha por vez inserindo o valor de e o nível de confiança verdadeiro. Finalmente, o objeto counter nos permite alterar o número da linha de save.true.conf dentro do loop.

Após o loop for, usamos a função plot() para plotar o valor de \(\pi\) no eixo \(x\) e o nível de confiança verdadeiro no eixo \(y\) usando as colunas apropriadas de save.true.conf. O argumento type = “l” instrui R a construir um gráfico de linha onde cada par de nível de confiança verdadeiro é conectado por uma linha. A função abline() desenha uma linha pontilhada horizontal em 0.95, que é o nível de confiança declarado.

O pacote binom em R também pode ser usado para calcular os verdadeiros níveis de confiança. A função binom.coverage() calcula o nível de confiança verdadeiro para um de cada vez, e a função binom.plot() plota os níveis de confiança verdadeiros em um conjunto de valores diferentes de \(\pi\).


1.2 Duas variáveis binárias


Consideramos agora a situação em que o mesmo ensaio de Bernoulli é medido em unidades que podem ser classificadas em grupos. O caso mais simples é quando uma população consiste em dois grupos, como fêmeas e machos, peixes de água doce e salgada, ou empresas americanas e estrangeiras. Abaixo estão dois exemplos com uma resposta binária em ensaios que formam dois grupos.


Exemplo 1.13: Arremesso de lance livre de Larry Bird.


Um lance livre é um arremesso no basquete onde o arremessador pode arremessar livremente (sem oposição de outro jogador) de um local específico na quadra. O arremesso é feito (um sucesso) ou perdido (uma falha). Na maioria das vezes, durante um jogo da National Basketball Association (NBA), os lances livres são feitos em pares. Isso significa que um arremessador de lance livre tem uma tentativa e, posteriormente, tem uma segunda tentativa, não importa o que aconteça na primeira.

O ex-jogador da NBA Larry Bird foi um dos mais bem sucedidos em fazer lances livres durante sua carreira com uma taxa de sucesso de 88,6%. Em comparação, a média da NBA durante esse período foi de cerca de 75%. A excelente taxa de sucesso de Bird está entre suas muitas conquistas, pelas quais ele foi reconhecido como um dos 50 maiores jogadores da história da NBA.

Durante as temporadas da NBA de 1980-81 e 1981-82, foram registrados os resultados das tentativas de lance livre de Bird arremessadas em duplas, e um resumo é mostrado na Tabela abaixo. Por exemplo, Bird fez sua primeira e segunda tentativa 251 vezes. Além disso, Bird fez a primeira tentativa, mas depois posteriormente perdeu a segunda tentativa 34 vezes. No geral, ele fez a primeira tentativa 285 vezes sem levar em conta o que aconteceu na segunda tentativa. No total, Bird acertou 338 pares de lances livres durante a temporada.

TABELA: os resultados dos lances livres de Larry Bird; fonte de dados Wardrop (1995). \[ \begin{array}{ccccc} & & \mbox{Segundo} & \mbox{lance} & \\\hline & & \mbox{Acerto} & \mbox{Erro} & \mbox{Total} \\\hline \mbox{Primeiro lance} & \mbox{Acerto} & 251 & 34 & 285 \\ & \mbox{Erro} & 48 & 5 & 53 \\ & \mbox{Total} & 299 & 39 & 338 \\ \hline \end{array} \]


Fãs e comentaristas de basquete muitas vezes especulam que os resultados de um segundo lance livre podem depender dos resultados do primeiro. Por exemplo, se um atirador erra a primeira tentativa, o desapontamento ou a determinação levarão a alterar sua abordagem para a segunda tentativa?

Se sim, então devemos ver que a probabilidade de sucesso na segunda tentativa é diferente dependendo se a primeira tentativa foi feita ou perdida. Assim, os dois grupos neste problema são formados pelos resultados da primeira tentativa, e os ensaios de Bernoulli que observamos são os resultados da segunda tentativa. Dados os dados da Tabela acima, investigaremos se o resultado da segunda tentativa depende do que acontece na primeira tentativa.


Exemplo 1.14: Ensaio clínico da vacina Salk.


Ensaios clínicos são realizados para determinar a segurança e eficácia de novos medicamentos. Freqüentemente, as respostas de segurança e eficácia são de natureza categórica; por exemplo, a resposta de eficácia pode ser simplesmente se um medicamento cura ou não um paciente de uma doença. Para garantir que um novo medicamento é realmente melhor do que não fazer nada (os pacientes às vezes melhoram sem intervenção), é essencial ter um grupo de controle no estudo. Isso é alcançado em ensaios clínicos randomizando os pacientes em dois grupos: novo medicamento ou controle. O grupo de controle geralmente é um placebo, que é administrado exatamente como o novo medicamento, mas não contém medicação.

Um dos mais famosos e maiores ensaios clínicos já realizados foi em 1954. Mais de 1,8 milhão de crianças participaram do ensaio clínico para determinar a eficácia da vacina contra a poliomielite desenvolvida por Jonas Salk (Francis et al., 1955). Embora o projeto real do estudo tenha suscitado debate (Brownlee, 1955; Dawson, 2004), renunciamos a essa discussão e focamos nos dados obtidos da parte randomizada e controlada por placebo do estudo. Os dados, apresentados na Tabela a seguir, mostram que 57 das 200.745 crianças do grupo vacinado desenvolveram poliomielite durante o período do estudo, em oposição a 142 das 201.229 crianças do grupo placebo. A questão de interesse para o ensaio clínico foi “A vacina ajuda a prevenir a poliomielite?” Desenvolveremos medidas de comparação nesta seção para responder a essa pergunta.

TABELA: Resultados do ensaio clínico da vacina Salk; fonte de dados Francis et al. (1955, pág. 25). \[ \begin{array}{cccc} & \mbox{Polio} & \mbox{Livre de Polio} & \mbox{Total} \\\hline \mbox{Vacina} & 57 & 200,688 & 200,745 \\ \mbox{Placeo} & 142 & 201,087 & 201,229 \\ \mbox{Total} & 199 & 401,775 & 401,974 \\ \hline \end{array} \]


1.2.1 Notação e modelo


O modelo e a notação seguem aqueles usados para uma única variável aleatória binomial. Começamos considerando duas variáveis aleatórias de Bernoulli separadas, \(Y_1\) e \(Y_2\), uma para cada grupo. As probabilidades de sucesso para os dois grupos são denotadas por \(\pi_1\) e \(\pi_2\), respectivamente. Observamos \(n_j\) tentativas de \(Y_j\) levando a \(w_j\) sucessos observados, \(j = 1,2\). Substituímos um subscrito por “+” para denotar uma soma nesse subscrito. Assim, \(n_+ = n_1 + n_2\) é o número total de tentativas e \(w_+ = w_1 + w_2\) é o número total de sucessos observados.

Essa notação está representada na Tabela abaixo. A tabela do lado direito é chamada de tabela de contingência bidirecional, porque fornece uma lista de todas as tabulações cruzadas possíveis de duas variáveis categóricas.

TABELA: Probabilidade e estruturas de contagem observadas para duas variáveis aleatórias binomiais independentes. \[ \begin{array}{ccccccccccc} & & \mbox{Resposta} & & & & \mbox{Resposta} & \\\hline & & 1 & 2 & & & & & 1 & 2 & \\\hline \mbox{Grupo} & 1 & \pi_1 & 1-\pi_1 & 1 & & \mbox{Grupo} & 1 & w_1 & n_1-w_1 & n_1 \\ \mbox{Grupo} & 2 & \pi_2 & 1-\pi_2 & 1 & & \mbox{Grupo} & 2 & w_2 & n_2-w_2 & n_2 \\ \hline & & & & & & & & w_+ & n_+-w_+ & n_+ \end{array} \]


Denotamos a variável aleatória que representa o número de sucessos no grupo \(j\) por \(W_j\) e escrevendo sua função de probabilidade como \[ P(W_j=w_j) = \binom{n_j}{w_j}\pi_j^{w_j}(1-\pi_j)^{n_j-w_j}, \qquad w_j=0,1,\cdots,n_j, \quad j=1,2\cdot \]

Assumimos que as duas variáveis aleatórias \(Y_1\) e \(Y_2\) são independentes, de modo que o resultado de uma não pode afetar o resultado da outra. No ensaio clínico da vacina Salk, por exemplo, isso significa que as crianças designadas para receber a vacina não podem transmitir imunidade ou doença às crianças do grupo placebo e vice-versa. Essa suposição é crítica no que se segue e, portanto, esse modelo é chamado de modelo binomial independente. Quando a independência não é satisfeita, outros modelos precisam ser usados para levar em conta a dependência entre as variáveis aleatórias, por exemplo, dados pareados.

Quando quisermos simular dados deste modelo, usaremos o código R como mostrado abaixo. Esse procedimento de amostragem será importante em breve quando usarmos essas contagens simuladas para avaliar procedimentos de inferência estatística, como intervalos de confiança, para determinar se eles funcionam conforme o esperado.


Exemplo 1.15:
Simulando contagens em uma tabela de contingência.


Considere o caso de \(\pi_1 = 0.2\), \(\pi_2 = 0.4\), \(n_1 = 10\) e \(n_2 = 10\). O código R abaixo mostra como simular um conjunto de contagens para uma tabela de contingência.

pi1 <- 0.2
pi2 <- 0.4
n1 <- 10
n2 <- 10
set.seed (8191)
w1 <- rbinom (n = 1, size = n1 , prob = pi1)
w2 <- rbinom (n = 1, size = n2 , prob = pi2)
c.table <- array ( data = c(w1 , w2 , n1 -w1 , n2 -w2), dim = c(2 ,2) ,
dimnames = list ( Group = c(1 ,2) , Response = c(1, 2)))
c.table
##      Response
## Group 1 2
##     1 1 9
##     2 3 7
c.table[1,1] # w1
## [1] 1
c.table[1,2] # n1 -w1
## [1] 9
c.table[1,] # w1 and n1 -w1
## 1 2 
## 1 9
sum(c.table[1,]) # n1
## [1] 10


Semelhante à Seção anterior, usamos a função rbinom() para simular valores para \(w_1\) e \(w_2\). Para formar a tabela de contingência, que chamamos de c.table, usamos a função array(). Seu argumento de dados contém as contagens dentro da tabela de contingência. Essas contagens são concatenadas usando a função c().

Observe que os dados são inseridos por colunas \((w_1,w_2, n_1-w_1,n_2-w_2)\). O argumento dim especifica as dimensões da tabela de contingência como número de linhas, número de colunas, onde a função c() é usada novamente. Finalmente, o argumento dimnames dá nomes para as medidas de linha e coluna. Os nomes são fornecidos em formato de lista, o que permite que vários objetos sejam vinculados. Nesse caso, os objetos são Group e Response que contêm os níveis das linhas e colunas, respectivamente.

Para esta amostra em particular, \(w_1 = 1\), \(n_1-w_1 = 9\), \(w_2 = 3\) e \(n_2-w_2 = 7\). Acessamos esses valores de dentro da tabela de contingência especificando um número de linha e coluna com c.table. Por exemplo, c.table[1,2] é igual a \(n_1-w_1\). Omitimos um número de coluna ou linha dentro de [ ] para que uma linha ou coluna inteira, respectivamente, seja exibida. As contagens somadas são encontradas usando a função sum() com as contagens apropriadas. Se quiséssemos repetir esse processo, digamos, 1.000 vezes, o argumento \(n\) das funções rbinom() seria alterado para 1.000. Cada tabela de contingência precisaria ser formada separadamente usando a função array().


Verossimilhança e estimativas


Os principais interesses neste problema são estimar as probabilidades de sucesso, \(\pi_1\) e \(\pi_2\), para cada grupo e comparar essas probabilidades. A estimativa de máxima verossimilhança novamente fornece uma solução conveniente e poderosa. Como \(Y_1\) e \(Y_2\) são independentes, \(W_1\) e \(W_2\) também são. A função de verossimilhança formada por variáveis aleatórias independentes é apenas o produto de suas respectivas funções de verossimilhança. Portanto, a função de verossimilhança é \[ L(\pi_1,\pi_2|w_1,w_2) = L(\pi_1|w_1)\times L(\pi_2|w_2)\cdot \]

Maximizar essa verossimilhança sobre \(\pi_1\) e \(\pi_2\) resulta nas estimativas “óbvias” \(\widehat{\pi}_1 = w_1/n_1\) e \(\widehat{\pi}_2 = w_2/n_2\), as proporções amostrais nos dois grupos. Em outras palavras, quando as variáveis aleatórias são independentes, cada probabilidade é estimada usando apenas os dados de seu próprio grupo.


Exemplo 1.16: Arremesso de lance livre de Larry Bird.


O objetivo deste exemplo é estimar a probabilidade de sucesso usando uma estrutura de tabela de contingência em R. Se os resultados do teste de Bernoulli já estiverem resumidos em contagens como na tabela do exemplo Larry Bird anterior, então uma tabela de contingência é criada em R usando a função array():

c.table <- array ( data = c(251 , 48, 34, 5) , dim = c(2 ,2) , 
            dimnames = list ( First = c(" made ", " missed ") , Second = c(" made ", " missed ")))
c.table
##           Second
## First       made   missed 
##    made       251       34
##    missed      48        5
list ( First = c(" made ", " missed ") , Second = c(" made ", " missed "))
## $First
## [1] " made "   " missed "
## 
## $Second
## [1] " made "   " missed "


Como os níveis das variáveis de linha e coluna têm nomes, usamos esses nomes no argumento dimnames.

As estimativas da probabilidade de sucessos ou proporções amostrais, são encontradas aproveitando-se de como o R realiza os cálculos:

rowSums(c.table ) # n1 and n2
##    made   missed  
##      285       53
pi.hat.table <- c.table/rowSums(c.table)
pi.hat.table
##           Second
## First          made     missed 
##    made    0.8807018 0.11929825
##    missed  0.9056604 0.09433962


A função rowSums() encontra a soma das contagens em cada linha para obter \(n_1\) e \(n_2\). Tomando a tabela de contingência de contagens divididas por essas somas de linhas, obtemos \(\widehat{\pi}_1\) e \(\widehat{\pi}_2\) na primeira coluna e \(1-\widehat{\pi}_1\) e \(1-\widehat{\pi}_2\) na segunda coluna.

Observe que R faz essa divisão corretamente tomando as contagens na primeira (segunda) linha de c.table dividida pelo primeiro (segundo) elemento de rowSums(c.table).


1.2.2 Intervalos de confiança para a diferença de duas probabilidades


Uma abordagem relativamente fácil para comparar \(\pi_1\) e \(\pi_2\) pode ser desenvolvida tomando sua diferença \(\pi_1-\pi_2\). A estimativa correspondente de \(\pi_1-\pi_2\) é \(\widehat{\pi}_1-\widehat{\pi}_2\). Cada estimativa de probabilidade de sucesso tem uma distribuição de probabilidade que é aproximada por uma distribuição normal em grandes amostras \[ \widehat{\pi}_j \sim N(\pi_j,\widehat{\mbox{Var}}(\widehat{\pi}_j)), \]

onde \(\widehat{\mbox{Var}}(\widehat{\pi}_j)=\widehat{\pi}_j(1-\widehat{\pi}_j)/n_j\), \(j=1,2\).

Como as combinações lineares de variáveis aleatórias normais são elas próprias variáveis aleatórias normais (Casella and Berger, 2002), a distribuição de probabilidade para \(\widehat{\pi}_1-\widehat{\pi}_2\) é aproximada por \[ N(\widehat{\pi}_1-\widehat{\pi}_2, \widehat{\mbox{Var}}(\widehat{\pi}_1-\widehat{\pi}_2)), \]

sendo \[ \widehat{\mbox{Var}}(\widehat{\pi}_1-\widehat{\pi}_2) = \widehat{\pi}_1(1-\widehat{\pi}_1)/n_1+\widehat{\pi}_2(1-\widehat{\pi}_2)/n_2\cdot \]

Essa distribuição forma a base para uma série de procedimentos de inferência.

O intervalo de confiança mais fácil de formar para \(\pi_1-\pi_2\) usa a aproximação normal para \(\widehat{\pi}_1-\widehat{\pi}_2\) diretamente para criar um intervalo Wald: \[ \widehat{\pi}_1-\widehat{\pi}_2 \pm Z_{1-\alpha/2}\sqrt{\dfrac{ \widehat{\pi}_1(1-\widehat{\pi}_1)}{n_1}+\dfrac{ \widehat{\pi}_2(1-\widehat{\pi}_2)}{n_2}}\cdot \]

Infelizmente, o intervalo de Wald para \(\pi_1-\pi_2\) tem problemas semelhantes em alcançar o nível de confiança declarado como o intervalo de Wald para \(\pi\). Investigaremos isso em breve.

Devido a esses problemas, vários outros intervalos de confiança para \(\pi_1-\pi_2\) foram propostos. Inspirados pelo bom desempenho geral do intervalo de Agresti-Coull para uma única probabilidade, Agresti and Caffo (2000) investigaram vários intervalos construídos como intervalos do tipo Wald em dados com diferentes números de sucessos e falhas adicionados. Eles descobriram que adicionar um sucesso e um fracasso para cada grupo resulta em um intervalo que faz um bom trabalho para atingir o nível de confiança declarado.

Especificamente, sejam \[ \widetilde{\pi}_1=\dfrac{w_1+1}{n_1+2} \qquad \mbox{e} \qquad \widetilde{\pi}_2=\dfrac{w_2+1}{n_2+2} \]

as estimativas alteradas de \(\pi_1\) e \(\pi_2\). Observe que, diferentemente de \(\widetilde{\pi}\) para o intervalo de Agresti-Coull, as estimativas \(\widetilde{\pi}\) e \(\widetilde{\pi}_2\) não mudam quando o nível de confiança muda. O intervalo de confiança \((1-\alpha)100\)% Agresti-Caffo para \(\pi_1-\pi_2\) é \[ \widetilde{\pi} - \widetilde{\pi}_2 \pm Z_{1-\alpha/2}\sqrt{\dfrac{ \widetilde{\pi}_1(1-\widetilde{\pi}_1)}{n_1+2}+\dfrac{ \widetilde{\pi}_2(1-\widetilde{\pi}_2)}{n_2+2}} \cdot \]

No geral, recomendamos o método Agresti-Caffo para uso geral.

Outros métodos de intervalo de confiança foram desenvolvidos de forma análoga ao caso de parâmetro único discutido anteriormente. Existe um intervalo baseado na inversão da estatística de teste escore para um teste de \(H_0 : \pi_1-\pi_2=d\), ou seja, determinar para quais valores \(d\) de \(\pi_1-\pi_2\) que a hipótese nula não é rejeitada. Isso acaba sendo um problema computacional consideravelmente mais difícil do que no caso de parâmetro único, porque \(H_0\) na verdade não especifica os valores de \(\pi_1\) e \(\pi_2\).

Da mesma forma, há um intervalo Bayesiano confiável de dois grupos semelhante ao método de Jeffreys, mas isso envolve distribuições que não são tão simples de calcular quanto a normal padrão. Detalhes desse intervalo estão disponíveis em Agresti and Min (2005a).


Exemplo 1.17: Arremesso de lance livre de Larry Bird.


O objetivo deste exemplo é calcular um intervalo de confiança para a diferença nas probabilidades de sucesso do segundo lance livre dados os resultados da primeira tentativa. Continuando o código anterior, obtemos o seguinte:

alpha <- 0.05
pi.hat1 <- pi.hat.table[1 ,1]
pi.hat2 <- pi.hat.table[2 ,1]
# Wald
var.wald <- pi.hat1*(1 - pi.hat1) / sum(c.table[1,]) + pi.hat2*(1 - pi.hat2) / sum(c.table[2 ,])
pi.hat1 - pi.hat2 + qnorm (p = c( alpha/2, 1- alpha/2))*sqrt(var.wald )
## [1] -0.11218742  0.06227017
# Agresti - Caffo
pi.tilde1 <- (c.table[1,1] + 1) / (sum(c.table[1 ,]) + 2)
pi.tilde2 <- (c.table [2,1] + 1) / (sum(c.table [2 ,]) + 2)
var.AC <- pi.tilde1*(1 - pi.tilde1 ) / (sum(c.table [1 ,]) + 2) + 
  pi.tilde2*(1 - pi.tilde2 ) / (sum(c.table [2 ,]) + 2)
pi.tilde1 - pi.tilde2 + qnorm (p = c( alpha /2, 1- alpha /2)) * sqrt (var.AC)
## [1] -0.10353254  0.07781192


O intervalo de Wald de 95% é \[ (-0.1122 < \pi_1-\pi_2 < 0.0623) \] e o intervalo de 95% Agresti-Caffo é \[ (-0.1035 < \pi_1-\pi_2 < 0.0778)\cdot \]

Os intervalos são um pouco semelhantes com o intervalo Wald sendo deslocado para a esquerda do intervalo Agresti-Caffo. Usando o intervalo Agresti-Caffo, com 95% de confiança, a diferença nas probabilidades de sucesso do segundo lance livre dado o resultado do primeiro está entre -0.1035 e 0.0778.

Como esse intervalo contém 0, não podemos detectar uma mudança na probabilidade de Bird de um segundo lance livre bem-sucedido após as primeiras tentativas feitas e perdidas. Isso significa que ou não há diferença, ou há uma diferença, mas não foi detectada nesta amostra.

A última situação pode ser causada por má sorte (uma amostra incomum) ou por um tamanho de amostra muito pequeno.

Os mesmos intervalos de confiança podem ser obtidos de outras maneiras. Primeiro, podemos usar o seguinte código quando os dados ainda não estiverem dentro do R por meio da função array():

w1 <- 251
n1 <- 285
w2 <- 48
n2 <- 53
alpha <- 0.05
pi.hat1 <- w1/n1
pi.hat2 <- w2/n2
var.wald <- pi.hat1 *(1 - pi.hat1 ) / n1 + pi.hat2 *(1 - pi.hat2 ) / n2
pi.hat1 - pi.hat2 + qnorm (p = c( alpha /2, 1- alpha /2)) * sqrt (var.wald )
## [1] -0.11218742  0.06227017


Em segundo lugar, a função prop.test() mostrada posteriormente fornece o intervalo de confiança de Wald como parte de sua saída. Por fim, a função wald2ci() no pacote PropCIs também calcula os intervalos de confiança de Wald e Agresti-Caffo.


Para encontrar um nível de confiança verdadeiro para um desses intervalos de confiança, é necessária a distribuição de probabilidade conjunta para todas as combinações possíveis de \((W_1,W_2)\). Essa distribuição é o produto de duas probabilidades binomiais porque essas variáveis aleatórias são independentes.

Para dados \(n_1\) e \(n_2\), existem \((n_1+1)(n_2+1)\) possíveis combinações observadas de \((w_1,w_2)\), e um intervalo de confiança pode ser calculado para cada uma dessas combinações. Para valores definidos de \(\pi_1\) e \(\pi_2\), alguns desses intervalos contêm \(\pi_1-\pi_2\) e outros não. Assim, o verdadeiro nível de confiança em \(\pi_1\) e \(\pi_2\), \(C(\pi_1,\pi_2)\), é a soma das probabilidades conjuntas para todos os intervalos que contêm \(\pi_1-\pi_2\): \[ C(\pi_1,\pi_2)=\sum_{w_2=0}^{n_2}\sum_{w_1=0}^{n_1} \pmb{I}(w_1,w_2)\binom{n_1}{w_1}\pi_1^{n_1}(1-\pi_1)^{n_1-w_1}\binom{n_2}{w_2}\pi_2^{n_2}(1-\pi_2)^{n_2-w_2}, \]

onde a função indicadora \(\pmb{I}(w_1,w_2)\) é 1 se o intervalo correspondente contiver \(\pi_1-\pi_2\) e \(\pmb{I}(w_1,w_2)\) for 0 caso contrário. Os detalhes do cálculo são fornecidos no próximo exemplo.


Exemplo 1.18: Intervalos Wald e Agresti-Caffo.


O verdadeiro nível de confiança para o intervalo de Wald pode ser encontrado de maneira semelhante à discutida anteriormente.

Considere o caso de \(\alpha= 0.05\), \(\pi_1 = 0.2\), \(\pi_2 = 0.4\), \(n_1 = 10\) e \(n_2 = 10\). Para encontrar todas as combinações possíveis de \((w_1,w_2)\), usamos o expand.grid(), que encontra todas as combinações possíveis dos argumentos, separados por vírgulas, dentro de seus parênteses. Repetimos esse mesmo processo para encontrar todos os valores possíveis de \((\widehat{\pi}_1,\widehat{\pi}_2)\) e \(P(W_1 = w_1, W_2 = w_2)\).

Abaixo segue o código R:

alpha <- 0.05
pi1 <- 0.2
pi2 <- 0.4
n1 <- 10
n2 <- 10
# All possible combinations of w1 and w2
w.all <- expand.grid(w1 = 0:n1 , w2 = 0: n2)
# All possible combinations of pi^_1 and pi^_2
pi.hat1 <- (0: n1)/n1
pi.hat2 <- (0: n2)/n2
pi.hat.all <- expand.grid (pi.hat1 = pi.hat1 , pi.hat2 = pi.hat2 )
# Find joint probability for w1 and w2
prob.w1 <- dbinom (x = 0:n1 , size = n1 , prob = pi1)
prob.w2 <- dbinom (x = 0:n2 , size = n2 , prob = pi2)
prob.all <- expand.grid ( prob.w1 = prob.w1 , prob.w2 = prob.w2)
pmf <- prob.all$prob.w1*prob.all$prob.w2
# P(W1 = w1 , W2 = w2)
head ( data.frame(w.all , pmf = round (pmf ,4)))
##   w1 w2    pmf
## 1  0  0 0.0006
## 2  1  0 0.0016
## 3  2  0 0.0018
## 4  3  0 0.0012
## 5  4  0 0.0005
## 6  5  0 0.0002


Assim temos, por exemplo, a probabilidade de observar \(P(W_1 = 1, W_2 = 0) = 0.0016\). Usando essas probabilidades, calculamos o nível de confiança verdadeiro para o intervalo:

var.wald <- pi.hat.all [ ,1] * (1- pi.hat.all [ ,1]) / n1 + pi.hat.all [ ,2] * (1- pi.hat.all [ ,2]) / n2
lower <- pi.hat.all [ ,1] - pi.hat.all [ ,2] - qnorm (p = 1- alpha /2) * sqrt (var.wald )
upper <- pi.hat.all [ ,1] - pi.hat.all [ ,2] + qnorm (p = 1- alpha /2) * sqrt (var.wald )
save <- 
  ifelse ( test = pi1 -pi2 > lower , yes = ifelse ( test = pi1 - pi2 < upper , yes = 1, no = 0) , no = 0)
sum( save *pmf )
## [1] 0.9281274
data.frame (w.all , round ( data.frame (pmf , lower , upper ) ,4) ,  save ) [1:15 ,]
##    w1 w2    pmf   lower  upper save
## 1   0  0 0.0006  0.0000 0.0000    0
## 2   1  0 0.0016 -0.0859 0.2859    0
## 3   2  0 0.0018 -0.0479 0.4479    0
## 4   3  0 0.0012  0.0160 0.5840    0
## 5   4  0 0.0005  0.0964 0.7036    0
## 6   5  0 0.0002  0.1901 0.8099    0
## 7   6  0 0.0000  0.2964 0.9036    0
## 8   7  0 0.0000  0.4160 0.9840    0
## 9   8  0 0.0000  0.5521 1.0479    0
## 10  9  0 0.0000  0.7141 1.0859    0
## 11 10  0 0.0000  1.0000 1.0000    0
## 12  0  1 0.0043 -0.2859 0.0859    1
## 13  1  1 0.0108 -0.2630 0.2630    1
## 14  2  1 0.0122 -0.2099 0.4099    1
## 15  3  1 0.0081 -0.1395 0.5395    0


Todos os intervalos Wald possíveis são calculados e a função ifelse() é usada para verificar se \(\pi_1-\pi_2 = 0.2 - 0.4 = -0.2\) está dentro de cada intervalo. O último quadro de dados mostra os primeiros 15 intervalos com os resultados desta verificação. As probabilidades correspondentes aos intervalos que contêm -0.2 são somadas para produzir o nível de confiança verdadeiro de 0.9281, que é menor que o nível declarado de 0.95.

Também podemos calcular o nível de confiança verdadeiro mantendo uma das probabilidades constante e permitindo que a outra varie. A figura abaixo mostra um gráfico em que \(\pi_2 = 0.4\) e \(\pi_1\) varia de 0.001 a 0.999 por 0.0005. As linhas Agresti-Caffo e Wald são desenhadas no gráfico simultaneamente usando a função plot() primeiro para os níveis de confiança verdadeiros de Wald e depois usando a função lines() para os níveis de confiança verdadeiros do Agresti-Caffo. A função legend() coloca a legenda no gráfico.

# Initial settings
alpha<-0.05
pi1<-0.2
pi2<-0.4
n1<-10
n2<-10
numb.bin.samples<-1000  # Number of binomial samples
###########################################################################
# Estimated true confidence level
  # Simulate w1 and w2
  set.seed(2349)
  w1<-rbinom(n = numb.bin.samples, size = n1, prob = pi1)
  w2<-rbinom(n = numb.bin.samples, size = n2, prob = pi2)
  pi.hat1<-w1/n1
  pi.hat2<-w2/n2
  # Wald
  var.wald<-pi.hat1*(1-pi.hat1) / n1 + pi.hat2*(1-pi.hat2) / n2
  lower<-pi.hat1 - pi.hat2 - qnorm(p = 1-alpha/2) * sqrt(var.wald)
  upper<-pi.hat1 - pi.hat2 + qnorm(p = 1-alpha/2) * sqrt(var.wald)
  # Intervals 1-5
  data.frame(w1, w2, lower, upper)[1:5,]
##   w1 w2      lower       upper
## 1  0  6 -0.9036363 -0.29636369
## 2  2  3 -0.4770066  0.27700660
## 3  1  4 -0.6560451  0.05604514
## 4  2  1 -0.2098975  0.40989752
## 5  3  4 -0.5157711  0.31577115
  # Calculate estimated true confidence level
  save<-ifelse(test = pi1-pi2 > lower,
               yes = ifelse(test = pi1-pi2 < upper, yes = 1, no = 0), no = 0)
  save[1:5]
## [1] 0 1 1 1 1
  true.conf<-mean(save)
  round(true.conf,4)
## [1] 0.92
  # Agresti-Caffo
  pi.tilde1<-(w1+1)/(n1+2)
  pi.tilde2<-(w2+1)/(n2+2)
  var.AC<-pi.tilde1*(1-pi.tilde1) / (n1+2) + pi.tilde2*(1-pi.tilde2) / (n2+2)
  lower.AC<-pi.tilde1 - pi.tilde2 - qnorm(p = 1-alpha/2) * sqrt(var.AC)
  upper.AC<-pi.tilde1 - pi.tilde2 + qnorm(p = 1-alpha/2) * sqrt(var.AC)
  save.AC<-ifelse(test = pi1-pi2 > lower.AC,
                  yes = ifelse(test = pi1-pi2 < upper.AC, yes = 1, no = 0), no = 0)
  save.AC[1:10]
##  [1] 1 1 1 1 1 1 1 1 1 1
  true.conf.AC<-mean(save.AC)
  round(true.conf.AC,4)
## [1] 0.962
###########################################################################
# Estimated true confidence level holding pi2 fixed at 0.3
  # Number of binomial samples - changed to reduce simulation variability (makes plot look nicer)  
  numb.bin.samples<-10000  
  pi1seq<-seq(from = 0.001, to = 0.999, by = 0.0005)
  # pi1seq<-0.2  # Testing
  # pi1seq<-seq(from = 0.1, to = 0.9, by = 0.1)  # Testing
  # Save true confidence levels in a matrix
  save.true.conf<-matrix(data = NA, nrow = length(pi1seq), ncol = 3)
  # Create counter for the loop
  counter<-1
  set.seed(2114)
  # Loop over each pi1 that the true confidence level is calculated on
  for(pi1 in pi1seq) {
    w1<-rbinom(n = numb.bin.samples, size = n1, prob = pi1)
    w2<-rbinom(n = numb.bin.samples, size = n2, prob = pi2)
    pi.hat1<-w1/n1
    pi.hat2<-w2/n2
    # Wald
    lower<-pi.hat1 - pi.hat2 - qnorm(p = 1-alpha/2) *
      sqrt(pi.hat1*(1-pi.hat1) / n1 + pi.hat2*(1-pi.hat2) / n2)
    upper<-pi.hat1 - pi.hat2 + qnorm(p = 1-alpha/2) *
      sqrt(pi.hat1*(1-pi.hat1) / n1 + pi.hat2*(1-pi.hat2) / n2)
    save<-ifelse(test = pi1-pi2 > lower,
                 yes = ifelse(test = pi1-pi2 < upper, yes = 1, no = 0), no = 0)
    wald<-mean(save)
    # Agresti-Caffo
    pi.tilde1<-(w1+1)/(n1+2)
    pi.tilde2<-(w2+1)/(n2+2)
    lower.AC<-pi.tilde1 - pi.tilde2 - qnorm(p = 1-alpha/2) *
            sqrt(pi.tilde1*(1-pi.tilde1) / (n1+2) +
              pi.tilde2*(1-pi.tilde2) / (n2+2))
    upper.AC<-pi.tilde1 - pi.tilde2 + qnorm(p = 1-alpha/2) *
            sqrt(pi.tilde1*(1-pi.tilde1) / (n1+2) +
              pi.tilde2*(1-pi.tilde2) / (n2+2))
    save.AC<-ifelse(test = pi1-pi2 > lower.AC,
                    yes = ifelse(test = pi1-pi2 < upper.AC, yes = 1, no = 0), no = 0)
    AC<-mean(save.AC)
    save.true.conf[counter,]<-c(pi1, wald, AC)
    counter<-counter+1
  }
  # Plot
  plot(x = save.true.conf[,1], y = save.true.conf[,2], xlab = expression(pi[1]),
       ylab = "Estimated true confidence level", type = "l", 
       ylim = c(0.85,1), lty = "solid", col = "blue",
       main = expression(paste("Verdadeiros níveis de confiança com ",n[1],"=10, ",n[2],"=10, ",
                               pi[2],"=0.2 e ",alpha,"=0.05")))
  lines(x = save.true.conf[,1], y = save.true.conf[,3], lty = "dashed", col = "red")
  abline(h = 1-alpha, lty = "dotted")
  legend(x = 0.1, y = 0.88, legend = c("Wald", "Agresti-Caffo"), lty = c("solid", "dashed"),
    bty = "n", col = c("blue", "red"))
  grid()


A figura acima mostra que o intervalo de Wald nunca atinge o nível de confiança verdadeiro! O intervalo Agresti-Caffo tem um nível de confiança verdadeiro entre 0.93 e 0.97. Incentivamos os leitores a alterar o valor predefinido para \(\pi_2\) no programa para examinar o que acontece em outras situações.

Por exemplo, o intervalo Wald é extremamente liberal e o intervalo Agresti-Caffo é extremamente conservador para \(\pi_1\) pequeno quando \(\pi_2 = 0.1\) com \(n_1 = 10\), \(n_2 = 10\) e \(\alpha= 0.05\).

######################################################################################

  # Alternatively, these true confidence levels can be found using the binom.coverage()
  #  function in the binom package
  library(binom)
  binom.coverage(p = 0.16, n = n, conf.level = 0.95, method = "asymptotic")
##       method    p  n coverage
## 1 asymptotic 0.16 40 0.894039
  # Plots of the true confidence levels can be found from the binom.plot()
  #  function in the binom package. Note that the method argument uses a function name
  #  that calculates the individual confidence interval (binom.asymp is within the binom 
  #  package. Also, the y-axis is the true confidence level and the x-axis is pi. The 
  #  lattice package is used for plotting so the usual arguments (like xlab = ) for basic 
  #  R plots do not always work.  
  binom.plot(n = 40, method = binom.asymp, np = 500, conf.level = 0.95)

#

1.2.3 Teste para a diferença de duas probabilidades


Um teste formal das hipóteses \(H_0 : \pi_1 -\pi_2 = 0\) vs. \(H_1 : \pi_1-\pi_2\neq 0\) pode ser realizado, novamente de várias maneiras. Um teste Wald usa uma estatística de teste \[ Z_w = \dfrac{\widehat{\pi}_1-\widehat{\pi}_2}{\sqrt{\dfrac{ \widetilde{\pi}_1(1-\widetilde{\pi}_1)}{n_1+2}+\dfrac{ \widetilde{\pi}_2(1-\widetilde{\pi}_2)}{n_2+2}}} \]

e compara essa estatística com a distribuição normal padrão. Observe que o denominador contém a variância estimada de \(\widehat{\pi}_1-\widehat{\pi}_2\) sem considerar a hipótese nula.

As distribuições de probabilidade para estatísticas de teste geralmente são calculadas assumindo que a hipótese nula é verdadeira. No presente contexto, isso significa que as probabilidades dos dois grupos são iguais e, portanto, uma variância estimada melhor do que a que está em \(Z_W\) pode ser calculada assumindo que \(\pi_1 = \pi_2\). Observe que essa condição implica que \(Y_1\) e \(Y_2\) tenham a mesma distribuição.

Assim, \(W_1\) e \(W_2\) são ambas contagens de sucessos da mesma variável aleatória de Bernoulli e, portanto, \(w_1\) e \(w_2\) podem ser combinados para representar \(w_+\) sucessos em \(n_+\) tentativas. Seja \(\widetilde{\pi}= w_+/n_+\) a probabilidade estimada de sucesso quando a hipótese nula for verdadeira. Então pode ser mostrado que \[ \widehat{\mbox{Var}}(\widehat{\pi}_1-\widehat{\pi}_2)=\widetilde{\pi}(1-\widetilde{\pi})(1/n_1+1/n_2)\cdot \]

Isso leva a um teste baseado na comparação da estatística \[ Z_0=\dfrac{\widehat{\pi}_1-\widehat{\pi}_2}{\sqrt{\widetilde{\pi}(1-\widetilde{\pi})(1/n_1+1/n_2)}} \]

para uma distribuição normal padrão. Este é o teste escore.

Um procedimento mais geral que é usado para comparar contagens observadas com contagens esperadas estimadas de qualquer modelo hipotético é o teste qui-quadrado de Pearson. Esse teste específico é frequentemente usado para realizar testes de hipóteses em uma configuração de tabelas de contingência e gastaremos muito mais tempo discutindo-o posteriormente. O modelo que está implícito em nossa hipótese nula é um único binômio com \(n_+\) tentativas e probabilidade de sucesso \(\widetilde{\pi}\).

A estatística de teste é formada pelo cálculo da \[ \mbox{"contagem observada - contagem esperada estimada"}^2\mbox{/"contagem esperada estimada"} \] sobre todas as contagens observadas, significando aqui sucessos e falhas nos dois grupos. O número esperado estimado de sucessos no grupo \(j\) sob o modelo de hipótese nula é \(n_j\widetilde{\pi}\) e, da mesma forma, o número esperado de falhas é \(n_j(1-\widetilde{\pi})\). Assim, a estatística do teste qui-quadrado de Pearson é \[ X^2=\sum_{j=1}^2 \left( \dfrac{(w_j-n_j\widetilde{\pi})^2}{n_j\widetilde{\pi}}+\dfrac{(n_j-w_j-n_j(1-\widetilde{\pi}))^2}{n_j(1-\widetilde{\pi})}\right)\cdot \]

Isso pode ser simplificado para \[ X^2=\sum_{j=1}^2 \dfrac{(w_j-n_j\widetilde{\pi})^2}{n_j\widetilde{\pi}(1-\widetilde{\pi})}\cdot \]

A estatística \(X^2\) tem uma distribuição de aproximadamente \(\chi^2_1\) quando \(n_1\) e \(n_2\) são grandes e quando a hipótese nula é verdadeira. Se a hipótese nula for falsa, as contagens observadas tendem a não se aproximar do esperado quando a hipótese nula for verdadeira; assim, grandes valores de \(X^2\) relativos à distribuição \(\chi^2_1\) levam a uma rejeição da hipótese nula.

Pode ser mostrado que os resultados do teste de qui-quadrado e escore de Pearson são idênticos para essa configuração, porque \(X^2 = Z^2_0\) e a distribuição \(\chi^2_1\) é equivalente à distribuição de uma variável aleatória normal padrão quadrada, por exemplo, \(Z^2_{0.975} = \chi^2_{1,0.95} = 3.84\).

Um teste da razão de verossimilhanças (LRT) também pode ser realizado. A estatística de teste pode ser mostrada como \[ -2\log(\Lambda) = \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad\qquad \qquad \qquad \qquad \qquad \qquad \\ -2\left(w_1\log\left(\dfrac{\widetilde{\pi}}{\widehat{\pi}_1}\right) +(n_1-w_1)\log\left(\dfrac{1-\widetilde{\pi}}{1-\widehat{\pi}_1}\right)+w_2\log\left(\dfrac{\widetilde{\pi}}{\widehat{\pi}_2}\right)+(n_2-w_2)\log\left(\dfrac{1-\widetilde{\pi}}{1-\widehat{\pi}_2}\right)\right), \]

onde assumimos que \(0\times \log(\infty) = 0\) por convenção. A hipótese nula é rejeitada se \[ -2\log(\Lambda ) > \chi^2_{ 1,1-\alpha}\cdot \]

Para todos esses testes, o uso da distribuição normal padrão ou qui-quadrado é baseado em aproximações de grandes amostras. Os testes são assintoticamente equivalentes, o que significa que fornecerão essencialmente os mesmos resultados em amostras muito grandes.

Em amostras pequenas, no entanto, as três estatísticas de teste podem ter distribuições sob a hipótese nula que são bem diferentes de suas aproximações. Larntz (1978) comparou a escore, o LRT e três outros testes em várias configurações de pequenas amostras e descobriu que o teste escore claramente mantém seu tamanho melhor do que os outros. Assim, o teste escore é recomendado aqui.


Exemplo 1.19: Arremesso de lance livre de Larry Bird.


O objetivo deste exemplo é mostrar como realizar o teste escore, teste qui-quadrado de Pearson e LRT em R. Podemos usar a função prop.test() para realizar os testes escore e qui-quadrado de Pearson:

prop.test (x = c.table , conf.level = 0.95 , correct = FALSE )
## 
##  2-sample test for equality of proportions without continuity correction
## 
## data:  c.table
## X-squared = 0.27274, df = 1, p-value = 0.6015
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  -0.11218742  0.06227017
## sample estimates:
##    prop 1    prop 2 
## 0.8807018 0.9056604

O valor do argumento para \(x\) é a tabela de contingência. Alternativamente, poderíamos ter atribuído a \(x\) um vetor com \(w_1\) e \(w_2\) dentro dele e usado um novo argumento \(n\) com um valor vetorial de \(n_1\) e \(n_2\).

O argumento correct = FALSE garante que a estatística de teste seja calculada conforme mostrado por \(Z_0\); caso contrário, uma correção de continuidade é aplicada para ajudar a garantir que o teste mantenha seu tamanho igual ou inferior a \(\alpha\).

A saída acima fornece o valor da estatística de teste como \(Z^2_0 = 0.27274\) e um \(p\)-valor de \(P(A > 0.27274) = 0.6015\), onde \(A\) tem uma distribuição \(\chi^2_1\). A decisão é não rejeitar a hipótese nula. A conclusão é que não há uma mudança significativa na porcentagem de sucesso do segundo lance livre de Bird em relação aos possíveis resultados da primeira tentativa.

Observe que a função chisq.test() e o método summary.table() também fornecem maneiras de realizar o teste qui-quadrado de Pearson. O código para o LRT é mostrado abaixo:

pi.bar <- colSums (c.table ) [1]/ sum(c.table )
log.Lambda <- c.table [1 ,1] * log(pi.bar / pi.hat.table [1 ,1]) + 
  c.table [1 ,2] * log ((1 - pi.bar) / (1- pi.hat.table [1 ,1])) +
  c.table [2 ,1] * log(pi.bar / pi.hat.table [2 ,1]) + 
  c.table [2 ,2] * log ((1 - pi.bar) /(1 - pi.hat.table [2 ,1]))
test.stat <- -2* log.Lambda
crit.val <- qchisq (p = 0.95 , df = 1)
p.val <- 1- pchisq (q = test.stat , df = 1)
round ( data.frame (pi.bar , test.stat , crit.val , p.val , row.names = NULL ), 4)
##   pi.bar test.stat crit.val p.val
## 1 0.8846    0.2858   3.8415 0.593


Sob a hipótese nula, a estimativa do parâmetro de probabilidade de sucesso é encontrada usando a soma das contagens na primeira coluna de c.table dividida pelo tamanho total da amostra e o resultado é colocado em pi.bar. O código para o objeto log.Lambda mostra como converter a maior parte da expressão da razão de verossimilhanças na sintaxe R correta. A estatística de teste transformada é \(-2\log(\Lambda ) = 0.2858\), e o \(p\)-valor é \(P(A > 0.2858) = 0.5930\).

Eles são bem impressos usando a função data.frame(), onde o valor do argumento row.names = NULL evita a impressão de um nome de linha errôneo. A conclusão geral é a mesma do teste escore. Observe que a estatística de teste e o \(p\)-valor poderiam ter sido calculados um pouco mais facilmente usando a função assocstats() do pacote vcd, e essa função também fornece a estatística de teste qui-quadrado de Pearson. Mostramos como usar esta função no programa correspondente a este exemplo.

Concluímos este exemplo com algumas notas adicionais:


1.2.4 Riscos relativos


O problema de basear a inferência em \(\pi_1-\pi_2\) é que ela mede uma quantidade cujo significado muda dependendo dos tamanhos de \(\pi_1\) e \(\pi_2\). Por exemplo, considere dois cenários hipotéticos em que a probabilidade de doença é listada para dois grupos de pessoas, digamos, para fumantes (\(\pi_1\)) e para não fumantes (\(\pi_2\)):

  1. \(\pi_1 = 0.51\) e \(\pi_2 = 0.50\)
  2. \(\pi_1 = 0.011\) e \(\pi_2 = 0.001\).

Em ambos os casos \(\pi_1-\pi_2 = 0.01\). Mas no primeiro cenário, um aumento de 0.01 devido ao tabagismo é bastante pequeno em relação ao risco já considerável de doença na população não fumante. Por outro lado, o cenário 2 apresenta fumantes com 11 vezes mais chance de doença do que os não fumantes. Precisamos ser capazes de transmitir as magnitudes relativas dessas mudanças melhor do que as diferenças permitem.

Neste caso, uma forma preferida de comparar duas probabilidades é através do risco relativo, \(RR = \pi_1/\pi_2\), assumindo \(\pi_2\neq 0\). Para o exemplo acima, \(RR = 0.011=0.001 = 11.0\) para o segundo cenário, o que significa que os fumantes têm 11 vezes mais chances de ter a doença do que os não fumantes. Alternativamente, poderíamos dizer que os fumantes têm 10 vezes mais chances de ter a doença do que os não fumantes. Por outro lado, para o primeiro cenário, \(RR = 0.51/0.50 = 1.02\), indicando que os fumantes têm apenas 2% mais chances ou 1.02 vezes mais de ter a doença. Observe que quando \(\pi_1 = \pi_2\), RR = 1.

Esses valores numéricos são baseados em probabilidades populacionais. Para obter um MLE para \(RR\), podemos fazer uso da propriedade de invariância dos MLEs que nos permite substituir as proporções observadas para as probabilidades, \(\widehat{RR} = \widehat{\pi_1}/\widehat{\pi_2}\), assumindo \(\widehat{\pi_2}\neq 0\). É esta estimativa que é frequentemente dada em reportagens que indicam riscos associados a certos fatores, como tabagismo ou obesidade.

Como o \(\widehat{RR}\) é um MLE, a inferência pode ser realizada usando os procedimentos usuais. Acontece, no entanto, que a aproximação normal é bastante pobre para MLEs que são razões, especialmente quando a estimativa no denominador pode ter variabilidade não desprezível como é o caso aqui. Portanto, a inferência baseada em uma aproximação normal para \(\widehat{RR}\) não é recomendada. No entanto, a aproximação normal é um pouco melhor para \[ \log(\widehat{RR}) = \log(\widehat{\pi_1}) - \log(\widehat{\pi_2}), \]

o MLE para \(\log(\widehat{RR})\), então a inferência é geralmente realizada na escala logarítmica.

A estimativa de variância de \(\log(\widehat{RR})\) pode ser derivada pelo método delta \[ \widehat{\mbox{Var}}(\log(\widehat{RR})) = \dfrac{1-\widehat{\pi}_1}{n_1\widehat{\pi}_1}+\dfrac{1-\widehat{\pi}_2}{n_2\widehat{\pi}_2}=\dfrac{1}{w_1}-\dfrac{1}{n_1}+\dfrac{1}{w_2}-\dfrac{1}{n_2}\cdot \]

Um \((1-\alpha)100\)% intervalo de confiança de Wald para o risco relativo da população é encontrado calculando primeiro o intervalo de confiança para \(\log(\pi_1/\pi_2)\), \[ \log\left( \dfrac{\widehat{\pi}_1}{\widehat{\pi}_2}\right) \pm Z_{1-\alpha/2} \sqrt{\dfrac{1}{w_1}-\dfrac{1}{n_1}+\dfrac{1}{w_2}-\dfrac{1}{n_2}}\cdot \]

A transformação exponencial é então usada para encontrar o intervalo de Wald para o próprio risco relativo: \[ \exp\left( \log\left( \dfrac{\widehat{\pi}_1}{\widehat{\pi}_2}\right) \pm Z_{1-\alpha/2} \sqrt{\dfrac{1}{w_1}-\dfrac{1}{n_1}+\dfrac{1}{w_2}-\dfrac{1}{n_2}}\right), \] onde \(\exp(\cdot )\) é o inverso da função logarítmica natural, ou seja, \(b = \exp(a)\) é equivalente a \(a = \log(b)\). Quando \(w_1\) e/ou \(w_2\) são iguais a 0, o intervalo de confiança não pode ser calculado. Um ajuste ad-hoc é adicionar uma pequena constante, como 0.5, à contagem de 0 células e ao total de linhas correspondente. Por exemplo, se \(w_1 = 0\), substitua \(w_1\) por 0.5 e \(n_1\) por \(n_1 + 0.5\).

Observe que o intervalo \(\log(\widehat{RR}) \pm Z_{1-\alpha/2} \sqrt{\widehat{\mbox{Var}}(\log(\widehat{RR}))}\) é simétrico, mas o intervalo de Wald para o risco relativo acima não é; isso pelo uso da transformação exponencial.


Exemplo 1.20: Ensaio clínico da vacina Salk.


O objetivo deste exemplo é calcular o risco relativo estimado e o intervalo de confiança para o risco relativo da população para determinar a eficácia da vacina Salk na prevenção da poliomielite. Abaixo está o código R usado para inserir os dados em uma matriz e realizar os cálculos necessários:

c.table <- array ( data = c(57 , 142 , 200688 , 201087) , dim = c(2 ,2) , 
                   dimnames = list ( Treatment = c(" vaccine ", " placebo ") ,
                   Result = c(" polio ", " polio free ")))
c.table
##            Result
## Treatment    polio   polio free 
##    vaccine       57       200688
##    placebo      142       201087
pi.hat.table <- c.table / rowSums (c.table )
pi.hat.table
##            Result
## Treatment         polio   polio free 
##    vaccine  0.0002839423    0.9997161
##    placebo  0.0007056637    0.9992943
pi.hat1 <- pi.hat.table [1 ,1]
pi.hat2 <- pi.hat.table [2 ,1]
round (pi.hat1 /pi.hat2 , 4)
## [1] 0.4024
round (1/( pi.hat1 /pi.hat2 ), 4) # inverted
## [1] 2.4852
alpha <- 0.05
n1 <- sum (c.table [1 ,])
n2 <- sum (c.table [2 ,])
# Wald confidence interval
var.log.rr <- (1- pi.hat1 )/( n1*pi.hat1 ) + (1- pi.hat2 )/( n2*pi.hat2 )
ci <- exp (log(pi.hat1 /pi.hat2 ) + qnorm (p = c( alpha /2, 1- alpha /2) ) * sqrt (var.log.rr))
round (ci , 4)
## [1] 0.2959 0.5471
rev( round (1/ci , 4)) # inverted
## [1] 1.8278 3.3792


Definindo o índice 1 para representar o grupo vacina e 2 o grupo placebo, encontramos \(\widehat{RR} = 0.40\). A probabilidade estimada de contrair poliomielite é apenas 0.4 vezes ou 40% menor para o grupo da vacina do que para o placebo.

Observe que usamos a palavra estimado com essa interpretação porque estamos usando estimativas de parâmetros. O intervalo de confiança é \(0.30 < RR < 0.55\). Portanto, com 95% de confiança, podemos dizer que a vacina reduz o risco populacional de poliomielite em 45-70%. Além disso, observe que a função exponencial é calculada usando exp() em R, por exemplo, exp(1) é 2.718.


A razão para o risco relativo é muitas vezes organizada de modo que \(\widehat{RR}\geq 1\). Isso permite uma interpretação atraente de que o grupo representado no numerador tem um risco que é, por exemplo, “11 vezes maior” que o grupo denominador. Isso pode ser mais fácil para um público-alvo apreciar do que a alternativa “0.091 vezes maior”. No entanto, o aplicativo pode ditar qual grupo deve ser o numerador, independentemente das proporções da amostra. Foi o caso do exemplo da vacina Salk, onde é natural pensar na vacina em termos de sua redução de risco.

Além disso, o risco relativo é geralmente uma medida mais útil do que a diferença entre probabilidades quando as probabilidades são bastante pequenas. É de uso limitado de outra forma. Por exemplo, se \(\pi_2 = 0.8\); então o risco relativo máximo possível é 1/0.8=1.25. Portanto, é útil ter uma estatística alternativa para comparar probabilidades que seja aplicável independentemente do tamanho das probabilidades. Isso faz parte da motivação para a próxima medida.


1.2.5 Odds ratios ou Razões de chaces


Nós nos concentramos até agora em usar probabilidades para medir a chance de um evento ocorrer. Odds podem ser usadas como uma medida semelhante. Odds são simplesmente a probabilidade de um sucesso dividido pela probabilidade de uma falha, \[ odds = \pi/(1-\pi)\cdot \]

Em algumas áreas de aplicação, como apostas, as probabilidades são usadas quase exclusivamente. Por exemplo, se uma probabilidade for \(\pi= 0.1\), então as chances de sucesso correspondentes serão \(0.1/(1-0.1) = 1/9\). Isso será chamado de “probabilidades de 9 para 1 contra”, porque a probabilidade de falha é 9 vezes a probabilidade de sucesso. Observe que há uma relação de 1 para 1 entre probabilidade e odds: se você conhece uma, pode encontrar a outra. Observe também que as odds não têm limite superior, ao contrário das probabilidades. Assim como os riscos relativos, as odds são estimadas com MLEs, substituindo as probabilidades por suas estimativas correspondentes.

Quando há dois grupos, podemos calcular as chances separadamente em cada grupo: \(odds_1 =\pi_1/(1-\pi_1)\) e \(odds_2 = \pi_2/(1-\pi_2)\). Em seguida, é feita uma comparação das chances usando uma razão de chances. Formalmente, é definido como \[ OR = \dfrac{odds_1}{odds_2}=\dfrac{\pi_1/(1-\pi_1)}{\pi_2(1-\pi_2)}=\dfrac{\pi_1(1-\pi_2)}{\pi_2(1-\pi_1)}\cdot \]

A razão de chances pode ser estimada substituindo os parâmetros com suas estimativas correspondentes para obter o MLE: \[ \widehat{OR} = \dfrac{\widehat{odds}_1}{\widehat{odds}_2}=\dfrac{\widehat{\pi}_1(1-\widehat{\pi}_2)}{\widehat{\pi}_2(1-\widehat{\pi}_1)}=\dfrac{(w_1/n_1)\Big((n_2-w_2)/n_2\Big)}{(w_2/n_2)\Big((n_1-w_1)/n_1\Big)}=\dfrac{w_1(n_2-w_2)}{w_2(n_1-w_1)}\cdot \]

Quando as contagens são escritas na forma de uma tabela de contingência, a razão de chances estimada é um produto das contagens na “diagonal” superior esquerda para inferior direita da tabela, dividida por um produto das contagens fora da diagonal.


Interpretação


Determinar se uma razão de chances é ou não igual a 1, maior que 1 ou menor que 1 costuma ser interessante. Uma razão de chances igual a 1 significa que as chances do grupo 1 são iguais às chances do grupo 2. Assim, podemos dizer que as chances não dependem do grupo; ou seja, as chances de sucesso são independentes da designação do grupo. Uma razão de chances maior que 1 significa que as chances de sucesso são maiores para o grupo 1 do que para o grupo 2. O oposto é verdadeiro para uma razão de chances menor que 1. Observe que \(OR = 1\) implica que \(RR = 1\) e \(\pi_1 =\pi_2\) e vice-versa. Por causa dessa equivalência, os procedimentos de teste de hipóteses da Seção 1.2.3 podem ser usados equivalentemente para testar \(H_0 : OR = 1\) vs. \(H_1 : OR \neq 1\).

Como as razões de chances são combinações de quatro probabilidades, elas são fáceis de interpretar erroneamente. A interpretação padrão de uma razão de chances estimada é a seguinte: \[ \mbox{As chances estimadas de sucesso são } \widehat{OR} \mbox{ vezes maiores no grupo 1 do que no grupo 2}\cdot \]

Obviamente, “sucesso”, “grupo 1” e “grupo 2” seriam substituídos por termos significativos no contexto do problema. Tal como acontece com os riscos relativos, alterar a ordem da razão, para que tenhamos \(\widehat{odds}_2/\widehat{odds}_1\), às vezes pode facilitar a interpretação. Referimo-nos a isso como invertendo a razão de chances. Com essa razão de chances invertida, podemos dizer agora que \[ \mbox{As chances estimadas de sucesso são } 1/\widehat{OR} \mbox{ vezes maior no grupo 2 do que no grupo 1}\cdot \]

A razão das chances de falha para os dois grupos também pode ser de interesse. Assim, \((1-\widehat{\pi}_1)/\widehat{\pi}_1\) é a chance de falha para o grupo 1, e \((1-\widehat{\pi}_2)/\widehat{\pi}_2\) é a chance de falha para o grupo 2. A razão do grupo 1 para o grupo 2 é \[ \dfrac{(1-\widehat{\pi}_1)/\widehat{\pi}_1}{(1-\widehat{\pi}_2)/\widehat{\pi}_2} = \dfrac{\widehat{\pi}_2(1-\widehat{\pi}_1)}{\widehat{\pi}_1(1-\widehat{\pi}_2)}, \]

que é simplesmente \(1/\widehat{OR}\) de antes. A interpretação torna-se então

\[ \mbox{As chances estimadas de falha são } \widehat{OR} \mbox{ vezes maior no grupo 1 do que no grupo 2}\cdot \]

Finalmente, também podemos inverter essa razão de duas chances de falha para obter a interpretação de \[ \mbox{As chances estimadas de falha são } \widehat{OR} \mbox{ vezes maiores no grupo 2 do que no grupo 1}, \]

onde a razão de chances é numericamente a mesma que tínhamos no início.

Essa “simetria” na interpretação das razões de chances pode ser confusa às vezes, mas tem uma vantagem muito distinta. Nos estudos de caso-controle, as medições são feitas retrospectivamente em números fixos de “casos” e “controles”, que são equivalentes a “sucessos” e “fracassos”, respectivamente. Ou seja, a proporção de casos e controles na amostra, por exemplo, pacientes que tiveram ou não AVC, é escolhida pelo pesquisador e normalmente não corresponde à proporção de casos na população. Se eles se enquadram ou não no grupo 1 ou 2, por exemplo, eles tomaram ou não um determinado medicamento no ano passado, é a medida de interesse.

Nesse contexto, as definições usuais de “grupo” e “resposta” são invertidas: queremos saber como tomar medicação pode se relacionar com a probabilidade de alguém ter um AVC, mas medimos como ter um AVC está relacionado com a probabilidade de eles tomaram o remédio. Se calcularmos uma razão de chances neste caso, então é fácil ver que se as chances de uso de medicamentos são \(x\) vezes maiores para aqueles que tiveram um AVC do que para aqueles que não tiveram, então as chances de ter um AVC também são \(x\) vezes maiores para usuários de medicamentos do que para não usuários de medicamentos! Assim, podemos obter uma resposta para a questão de interesse sem realmente medir as chances de ter um acidente vascular cerebral para qualquer grupo de medicamentos.

De fato, dada apenas a razão de chances, não podemos estimar diretamente nenhuma de suas probabilidades constituintes. É incorreto interpretar uma razão de chances como relacionada diretamente à probabilidade de sucesso. Não podemos dizer que “A probabilidade estimada de sucesso é \(\widehat{OR}\) vezes maior no grupo 1 do que no grupo 2”. Essa é uma desvantagem substancial de usar razões de probabilidade em vez de riscos relativos, porque a maioria das pessoas se sente mais à vontade discutindo probabilidades do que odds.

No entanto, é fácil ver que \(OR\approx RR\) se as probabilidades dos dois grupos forem bastante baixas. Nesse caso, ambas as probabilidades de falha são quase 1 e se cancelam na razão de chances. Por esta razão, as razões de chance podem ser usadas como um substituto para os riscos relativos quando os sucessos contados são relativamente incomuns, por exemplo, ambos \(\pi_j \leq 0.1\).


Intervalos confidenciais


Assim como com os riscos relativos, a distribuição de probabilidade de \(\log(\widehat{OR})\) é melhor aproximada por uma distribuição normal do que a própria distribuição de probabilidade de \(\widehat{OR}\). A variância de \(\log(\widehat{OR})\) é novamente encontrada usando o método delta e tem uma forma particularmente fácil: \[ \widehat{\mbox{Var}}(\log(\widehat{OR}))=\dfrac{1}{w_1}+\dfrac{1}{n_1-w_1}+\dfrac{1}{w_2}+\dfrac{1}{n_2-w_2}\cdot \]

O \((1-\alpha)100\)% intervalo de confiança de Wald para \(OR\) torna-se \[ \exp\left( \log(\widehat{OR})\pm Z_{1-\alpha/2}\sqrt{\dfrac{1}{w_1}+\dfrac{1}{n_1-w_1}+\dfrac{1}{w_2}+\dfrac{1}{n_2-w_2}}\right)\cdot \]

Lui and Lin (2003) mostram que o intervalo é conservador em relação ao seu verdadeiro nível de confiança, onde o grau de sua conservatividade é dependente de \(\pi_1\), \(\pi_2\) e \(n_+\). O verdadeiro nível de confiança está um pouco acima do nível declarado na maioria das vezes. No entanto, o intervalo é muito conservador para pequenos e grandes \(\pi_1\) e \(\pi_2\) e pequenos \(n_+\). Esses são casos que provavelmente produzirão contagens de células muito pequenas na tabela de contingência observada correspondente, e isso leva a alguma instabilidade nas estimativas de \(OR\) e sua variância correspondente.

Contagens zero em qualquer lugar na tabela de contingência fazem com que a razão de chances estimada seja 0 ou indefinida e a variância correspondente seja indefinida. Pequenos ajustes nas contagens geralmente são feitos nessas situações. Por exemplo, podemos usar \[ \widetilde{OR} = \dfrac{(w_1+0.5)(n_2-w_2+0.5)}{(w_2+0.5)(n_1-w_1+0.5)} \]

como uma estimativa da razão de chances.

A adição de 0.5 a cada contagem de células leva a razão de chances estimada a ser mais próxima de 1 do que sem o ajuste. A variancia estimada para \(\log(\widehat{OR})\) também pode ser alterada para ter 0.5 adicionado a cada contagem dentro dela. Um intervalo de confiança de Wald ajustado pode então ser formado com \(\widetilde{OR}\) e essa nova variância fazendo as substituições apropriadas. Observe que algumas pessoas defendem fazer esse ajuste o tempo todo, mesmo quando não há contagens de 0.

A principal vantagem dessa abordagem é que o intervalo de confiança é mais curto, devido às observações “adicionais”, do que sem esse ajuste. Lui and Lin (2003) mostram esse intervalo de Wald ajustado para atingir um nível de confiança verdadeiro muito semelhante ao intervalo sem o ajuste.


Exemplo 1.21: Ensaio clínico da vacina Salk.


O objetivo deste exemplo é estimar a razão de chances que descreve a relação entre vacina e poliomielite. Codificamos os dados para que \(w_1\) seja a contagem de crianças contraindo poliomielite no grupo da vacina e \(w_2\) seja a mesma contagem para placebo. Continuando o código anterior, obtemos o seguinte:

# OR
OR.hat<-c.table[1,1]*c.table[2,2] / (c.table[2,1]*c.table[1,2])
round(OR.hat, 2)
## [1] 0.4
round(1/OR.hat, 2)
## [1] 2.49
alpha<-0.05
var.log.or<-1/c.table[1,1] + 1/c.table[1,2] + 1/c.table[2,1] + 1/c.table[2,2]
OR.CI<-exp(log(OR.hat) + qnorm(p = c(alpha/2, 1-alpha/2)) * sqrt(var.log.or))
round(OR.CI, 2)
## [1] 0.30 0.55
rev(round(1/OR.CI, 2))
## [1] 1.83 3.38
# Another way to get the OR
#  Note that the function below automatically adds 0.5 to each cell for variance 
#  but not for the odds ratio itself unless there are 0 counts (see code in function).
library(vcd)  # Visualizing categorical data package
save.OR<-oddsratio(x = c.table, log = TRUE)
attributes(save.OR)  # names( ) does not work
## $names
## [1] "coefficients" "dimnames"     "dim"          "vcov"         "contrasts"   
## [6] "log"         
## 
## $class
## [1] "loddsratio"
summary(save.OR) 
## 
## z test of coefficients:
## 
##                                          Estimate Std. Error z value  Pr(>|z|)
##  vaccine : placebo / polio : polio free  -0.91079    0.15683 -5.8074 6.343e-09
##                                             
##  vaccine : placebo / polio : polio free  ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
confint(save.OR, level = 0.95)
##                                              2.5 %     97.5 %
##  vaccine : placebo / polio : polio free  -1.218173 -0.6034058
OR.tilde<-(c.table[1,1]+0.5)*(c.table[2,2]+0.5)/((c.table[1,2]+0.5)*(c.table[2,1]+0.5))
log(OR.tilde)  # Does not match vcd's log(OR^)
## [1] -0.9055709
sqrt(1/(c.table[1,1]+0.5) + 1/(c.table[2,2]+0.5) + 
       1/(c.table[1,2]+0.5) + 1/(c.table[2,1]+0.5))  # Matches vcd's standard error (sqrt(var^(log(OR^)))
## [1] 0.1562652


Encontramos \(\widehat{OR} = 0.40\), então as chances estimadas de contrair poliomielite são 0.4 vezes ou 60% menores quando a vacina é administrada do que quando um placebo é administrado. O intervalo de confiança de 95% é \(0.30 < OR < 0.55\). Como 1 não está dentro do intervalo, há evidências suficientes para indicar que a vacina diminui as verdadeiras chances e, portanto, a probabilidade de contrair poliomielite. Observe que não incluímos a palavra estimada na última frase porque interpretamos \(OR\) em vez de \(\widehat{OR}\).

Alternativamente, poderíamos interpretar as razões de chance invertidas em termos de proteção devido à vacinação: As chances estimadas de estar livre da pólio são 2.49 vezes maiores quando a vacina é administrada do que quando o placebo é administrado. O intervalo de confiança de 95% para a razão de chances é (1.83; 3.38). Observe que a função rev() inverte a ordem do elementos dentro do objeto entre parênteses. Isso é feito para que o valor mais baixo do intervalo de confiança apareça primeiro no vetor resultante.

Esses valores calculados são essencialmente os mesmos de quando exploramos o mesmo exemplo usando o risco relativo. Isso ocorre porque (a) a probabilidade de contração da pólio é muito baixa em ambos os grupos (cerca de 1/1.000) e (b) o tamanho da amostra é muito grande aqui. Em amostras menores, as estimativas \(\widehat{RR}\) e \(\widehat{OR}\) ainda seriam bastante semelhantes, mas os intervalos de confiança teriam diferido. A diferença nos intervalos de confiança decorre da diferença nas variâncias: em amostras muito grandes, tanto \(1/(n_j - w_j)\) quanto \(-1/n_j\) são insignificantes em relação a \(1/w_j\), \(j = 1,2\).


Exemplo 1.22: Arremesso de lance livre de Larry Bird.


O objetivo deste exemplo é calcular uma razão de chances comparando as chances de uma segunda tentativa de lance livre bem-sucedida quando a primeira é feita e quando a primeira é perdida. A razão de chances estimada é \(\widehat{OR} = 0.77\), e o intervalo de confiança de 95% para OR é (0.29, 2.07). Como a razão de chances estimada é menor que 1, decidimos invertê-la para ajudar na sua interpretação.

As chances estimadas de uma segunda tentativa de lance livre bem-sucedida são 1.30 vezes maiores do que quando o primeiro lance livre é perdido do que quando o primeiro lance livre é feito. Além disso, com 95% de confiança, as chances de uma tentativa bem-sucedida de lance livre são entre 0.48 e 3.49 vezes maiores do que quando o primeiro lance livre é perdido do que quando o primeiro lance livre é feito. Como 1 está dentro do intervalo, não há evidências suficientes para indicar uma diferença real nas chances de sucesso.


1.2.6 Dados de pares combinados


A Seção 1.2.2 examina um intervalo de confiança para a diferença de duas probabilidades. Esses intervalos foram desenvolvidos usando o fato de que \(W_1\) e \(W_2\) são variáveis aleatórias independentes. Existem outras situações em que as duas probabilidades comparadas são variáveis aleatórias dependentes. Isso ocorre com dados de pares combinados, onde duas observações de resposta binária, digamos \(X\) e \(Y\), são feitas em cada unidade amostrada, e uma comparação desejada entre as probabilidades de sucesso para \(X\) e \(Y\) envolve duas estatísticas correlacionadas. Abaixo está um exemplo para ilustrar este problema.


Exemplo 1.23:
Procedimentos de diagnóstico de câncer de próstata.


Zhou and Qin (2005) discutem um estudo que foi usado para comparar a precisão diagnóstica de ressonância magnética (MRI) e ultra-som em pacientes que foram estabelecidos como tendo câncer de próstata localizado por um teste padrão-ouro. Os dados de um local de estudo são fornecidos na tabela abaixo. Estamos interessados em comparar a probabilidade de um diagnóstico localizado por uma ressonância magnética com a probabilidade de um diagnóstico localizado por um ultrassom. Observe que, como ambos os procedimentos são realizados nos mesmos pacientes, essas estimativas de probabilidade - 10/16 e 7/16, respectivamente - são baseadas nas mesmas contagens e, portanto, suas estatísticas correspondentes não são independentes. Isso viola a suposição de contagens binomiais independentes usadas no modelo da Seção 1.2.2.

TABELA. Diagnósticos estabelecidos usando tecnologia de ressonância magnética e ultra-som entre indivíduos que realmente têm câncer de próstata localizado. A fonte de dados é Zhou and Qin (2005). \[ \begin{array}{ccccc} & & & \mbox{Ultrassom} & \\\hline & & & \mbox{Localizado} & \mbox{Avançado} & & \\\hline \mbox{MRI} & \mbox{Localizado} & 4 & 6 & 10\\ & \mbox{Avançado} & 3 & 3 & 6 \\ \hline & & 7 & 9 & 16 \end{array} \]


Para desenvolver os métodos de análise apropriados, observe que qualquer combinação de duas medições pode ocorrer para cada sujeito. Ou seja, podemos observar \[ (X = 1, Y = 1), \quad (X = 1, Y = 2), \quad (X = 2, Y = 1) \qquad \mbox{ou} \qquad (X = 2, Y = 2), \]

onde os valores de \(X\) e \(Y\) correspondem aos números de linha e coluna, respectivamente.

Cada um desses eventos tem probabilidade \(\pi_{ij} = P(X = i, Y = j)\); \(i = 1,2\), \(j = 1,2\). Estas são chamadas de probabilidades conjuntas. Como o exemplo anterior destaca, estamos interessados em resumos dessas probabilidades conjuntas, e não nas próprias probabilidades conjuntas. Ou seja, estamos interessados nas probabilidades marginais, \(\pi_{1+} = \pi_{11} + \pi_{12}\) e \(\pi_{+1} = \pi_{11} + \pi_{21}\), assim chamadas porque se referem a contagens nas margens da tabela. Com relação ao exemplo de diagnóstico de câncer de próstata, \(\pi_{1+}\) é a probabilidade de um diagnóstico localizado por ressonância magnética (MRI) e \(\pi_{+1}\) é a probabilidade de um diagnóstico localizado por ultrassom.

Suponha que \(n\) unidades sejam classificadas cruzadamente de acordo com \(X\) e \(Y\). Em seguida, observamos as contagens \(n_{11}\), \(n_{12}\), \(n_{21}\) e \(n_{22}\), sendo \(n_{ij}\) o número de ocorrências conjuntas de \((X = i, Y = j)\). Isso é exatamente como a configuração binomial, exceto com quatro resultados possíveis em vez de dois. O modelo que descreve as contagens produzidas dessa maneira é chamado de distribuição multinomial. Por enquanto, basta saber que MLEs e variâncias seguem o mesmo padrão do modelo binomial: \[ \widehat{\pi}_{ij} = n_{ij}/n, \] com \[ \widehat{\mbox{Var}}(\widehat{\pi}_{ij}) = \widehat{\pi}_{ij}(1-\widehat{\pi}_{ij})/n, \]

\[ \widehat{\pi}_{1+} = (n_{11}+n_{12})/n \]

com \[ \widehat{\mbox{Var}}(\widehat{\pi}_{1+}) = \widehat{\pi}_{1+}(1-\widehat{\pi}_{1+})/n \] e \[ \widehat{\pi}_{+1} = (n_{11}+n_{21})/n \]

com \[ \widehat{\mbox{Var}}(\widehat{\pi}_{+1}) = \widehat{\pi}_{+1}(1-\widehat{\pi}_{+1})/n, \] onde \[ n = n_{11} + n_{12} + n_{21} + n_{22}\cdot \]

A comparação de probabilidades marginais em pares combinados geralmente é feita através da diferença \(\pi_{1+}-\pi_{+1}\), porque essa comparação é matematicamente muito mais fácil de trabalhar do que outras comparações possíveis. Em algumas situações, especialmente com probabilidades marginais muito pequenas, a razão das duas probabilidades pode ser mais interessante. Veja Bonett and Price (2006) para saber como formar intervalos para razões de probabilidades marginais.

Se \(\pi_{+1}-\pi_{1+} = 0\), os duas probabilidades marginais são iguais. Alternativamente, se \(\pi_{+1}-\pi{1+} > 0\), isso significa que o sucesso, categoria de primeira resposta, é mais provável de acontecer para \(Y\) do que para \(X\). O oposto é verdadeiro se \(\pi_{+1}-\pi_{1+} < 0\).


Intervalo de confiança para a diferença de probabilidades marginais


Para desenvolver um intervalo de confiança para \(\pi_{+1}-\pi_{1+}\), começamos com o MLE \(\widehat{\pi}_{+1}-\widehat{\pi}_{1+}\) e encontramos sua distribuição de probabilidade. Mais uma vez nos valemos dos fatos que:

No entanto, esse problema difere daquele apresentado na Seção 1.2.2 porque as duas variáveis aleatórias \(\widehat{\pi}_{+1}\) e \(\widehat{\pi}_{1+}\) são dependentes. Portanto temos que \[ \widehat{\mbox{Var}}(\widehat{\pi}_{+1}-\widehat{\pi}_{1+})=\widehat{\mbox{Var}}(\widehat{\pi}_{+1})+\widehat{\mbox{Var}}(\widehat{\pi}_{1+})-2\widehat{\mbox{Cov}}(\widehat{\pi}_{+1},\widehat{\pi}_{1+})\cdot \]

Usando a distribuição multinomial subjacente para as contagens de células, a covariância pode ser mostrada como \[ \widehat{\mbox{Cov}}(\widehat{\pi}_{+1},\widehat{\pi}_{1+}) = \dfrac{\widehat{\pi}_{11}\widehat{\pi}_{22}-\widehat{\pi}_{12}\widehat{\pi}_{21}}{n}\cdot \]

Usando as variâncias para as probabilidades marginais dadas anteriormente, obtemos \[ \widehat{\mbox{Var}}(\widehat{\pi}_{+1}-\widehat{\pi}_{1+})=\dfrac{\widehat{\pi}_{+1}(1-\widehat{\pi}_{+1})+ \widehat{\pi}_{1+}(1-\widehat{\pi}_{1+}) -2(\widehat{\pi}_{11}\widehat{\pi}_{22}-\widehat{\pi}_{12}\widehat{\pi}_{21})}{n} = \dfrac{\widehat{\pi}_{21}+\widehat{\pi}_{12}+(\widehat{\pi}_{21}-\widehat{\pi}_{12})^2}{n}\cdot \]

Isso leva ao \((1-\alpha)100\)% intervalo de confiança de Wald 100% para \({\pi}_{+1}-{\pi}_{1+}\): \[ \widehat{\pi}_{+1}-\widehat{\pi}_{1+} \pm Z_{1-\alpha/2} \sqrt{\widehat{\mbox{Var}}(\widehat{\pi}_{+1}-\widehat{\pi}_{1+})}\cdot \]

Não surpreendentemente, o intervalo de Wald tem problemas para atingir o nível de confiança declarado. Houve uma série de outros intervalos propostos e estes são resumidos e avaliados em artigos como Newcombe (1998), Tango (1998) e Zhou and Qin (2005). Focamos no intervalo proposto por Agresti and Min (2005b) devido à sua simplicidade e desempenho.

Motivados pelo intervalo Agresti-Caffo, Agresti and Min (2005b) propõem adicionar 0.5 a cada contagem na tabela de contingência e calcular um intervalo de confiança de Wald usando essas contagens ajustadas. Especificamente, seja \(\widetilde{\pi}_{ij} = (n_{ij} +0.5)/(n+2)\), de modo que \[ \widetilde{\pi}_{+1} = (n_{11} +n_{21} +1)/(n+2) \] e \[ \widetilde{\pi}_{1+} = (n_{11} + n_{12} + 1)/(n + 2)\cdot \]

O intervalo \((1-\alpha)100\)% Agresti-Min é \[ \widehat{\pi}_{+1}-\widehat{\pi}_{1+} \pm Z_{1-\alpha/2} \sqrt{\dfrac{\widetilde{\pi}_{+1}(1-\widetilde{\pi}_{+1})+ \widetilde{\pi}_{1+}(1-\widetilde{\pi}_{1+}) -2(\widetilde{\pi}_{11}\widetilde{\pi}_{22}-\widetilde{\pi}_{12}\widetilde{\pi}_{21})}{n}}\cdot \]

Agresti and Min (2005b) examinaram a adição de outras constantes a cada célula e concluíram que 0.5 era o melhor em termos de nível de confiança verdadeiro e comprimento esperado. Veja a Seção 3 de seu artigo para mais justificativas com relação a um intervalo escore e um intervalo Bayesiano confiável. A principal desvantagem de um intervalo escore aqui é que ele não possui uma expressão de forma fechada e às vezes é mais longo que o intervalo Agresti-Min. A função scoreci.mp() no pacote PropCIs calcula o intervalo escore.

Mostramos um exemplo de sua utilização no programa correspondente ao exemplo de diagnóstico de câncer de próstata.


Exemplo 1.24:
Procedimentos de diagnóstico de câncer de próstata.


O objetivo deste exemplo é calcular intervalos de confiança para \({\pi}_{+1}-{\pi}_{1+}\). Começamos inserindo os dados em R usando a função array():

c.table <- array ( data = c(4, 3, 6, 3) , dim = c(2 ,2) , 
                    dimnames = list (MRI = c(" Localized ", " Advanced "), 
                                     Ultrasound = c(" Localized ", " Advanced ")))
c.table
##              Ultrasound
## MRI            Localized   Advanced 
##    Localized            4          6
##    Advanced             3          3
n <- sum (c.table )
pi.hat.plus1 <- sum(c.table [ ,1])/n
pi.hat.1plus <- sum(c.table [1 ,])/n
data.frame (pi.hat.plus1 , pi.hat.1plus , diff = pi.hat.plus1 - pi.hat.1plus )
##   pi.hat.plus1 pi.hat.1plus    diff
## 1       0.4375        0.625 -0.1875


A estimativa de \({\pi}_{+1}-{\pi}_{1+}\) é -0.1875. Para encontrar o intervalo de confiança para \({\pi}_{+1}-{\pi}_{1+}\), usamos as funções diffpropci.Wald.mp() e diffpropci.mp() do pacote PropCIs:

library ( package = PropCIs )
diffpropci.Wald.mp(b = c.table [1 ,2] , c = c.table [2 ,1] , n = sum (c.table ), conf.level = 0.95)
## 
## 
## 
## data:  
## 
## 95 percent confidence interval:
##  -0.5433238  0.1683238
## sample estimates:
## [1] -0.1875
diffpropci.mp(b = c.table [1 ,2] , c = c.table [2 ,1] , n = sum (c.table ), conf.level = 0.95)
## 
## 
## 
## data:  
## 
## 95 percent confidence interval:
##  -0.5022786  0.1689453
## sample estimates:
## [1] -0.1666667


Os argumentos \(b\) e \(c\) nas funções correspondem a \(n_{12}\) e \(n_{21}\), respectivamente, o que leva ao cálculo de \(\widehat{\pi}_{+1}-\widehat{\pi}_{1+}=\widehat{\pi}_{21}-\widehat{\pi}_{12}\). A ordenação das estimativas nesta expressão é importante para determinar que o intervalo é para \({\pi}_{+1}-{\pi}_{1+}\) em vez de \({\pi}_{1+}-{\pi}_{+1}\).

Os intervalos 95% Wald e Agresti-Min são \[ -0.5433 < {\pi}_{+1}-{\pi}_{1+} < 0.1683 \] e \[ -0.5023 < {\pi}_{1+}-{\pi}_{1+} < 0.1689, \] respectivamente.

Mostramos no programa correspondente para este exemplo como calcular esses intervalos codificando as fórmulas diretamente. No geral, vemos que o intervalo é bastante amplo, deixando muita incerteza sobre a diferença real - especialmente, observe que 0 está dentro do intervalo.


Teste para a diferença entre probabilidades marginais


Especialmente para o exemplo de diagnóstico de câncer de próstata, saber se \(\pi_{+1}-\pi_{1+}\) é diferente de 0 é um caso especial muito importante do problema de dados de pares correspondentes. Para testar formalmente \(H_0 : \pi_{+1}-\pi_{1+} =0\) vs. \(H_0 : \pi_{+1}-\pi_{1+} \neq 0\), poderíamos usar a estatística do teste de Wald de \[ Z_0 = \dfrac{\widehat{\pi}_{+1}-\widehat{\pi}_{1+}-0}{\sqrt{\widehat{\mbox{Var}}(\widehat{\pi}_{+1}-\widehat{\pi}_{1+})}}=\dfrac{\widehat{\pi}_{21}-\widehat{\pi}_{12}}{\sqrt{\Big(\widehat{\pi}_{21}+\widehat{\pi}_{12}+(\widehat{\pi}_{21}-\widehat{\pi}_{12})^2\Big)/n}} \cdot \]

Em vez disso, uma estatística de teste diferente é frequentemente usada, incorporando simplificações que ocorrem sob a hipótese nula. Como \({\pi}_{+1}-{\pi}_{1+} = 0\) implica que \({\pi}_{21}={\pi}_{12}\), a variância no denominador pode ser simplificada para \((\widehat{\pi}_{21}-\widehat{\pi}_{12})/n\). Se elevarmos ao quadrado a estatística resultante, obteremos o que é conhecido como estatística de teste de McNemar (McNemar, 1947): \[ M=\dfrac{(\widehat{\pi}_{21}-\widehat{\pi}_{12})^2}{(\widehat{\pi}_{21}+\widehat{\pi}_{12})/n}=\dfrac{(n_{21}-n_{12})^2}{n_{21}+n_{12}}\cdot \]

Sob a hipótese nula, a estatística tem um valor aproximado da distribuição \(\chi^2_1\) para grandes amostras. Quando \(\pi_{21}\neq \pi_{12}\), a hipótese nula é falsa, e esperaríamos que isso fosse refletido por grandes valores de \((\widehat{\pi}_{21}-\widehat{\pi}_{12})^2\) em relação ao denominador. Assim, a hipótese nula é rejeitada quando \(M\) tem um valor observado incomumente grande em relação à distribuição \(\chi^2_1\); em outras palavras, rejeite \(H_0\) se \(M>\chi^2_{1,1-\alpha}\).


Exemplo 1.25:
Procedimentos de diagnóstico de câncer de próstata.


Realizamos formalmente um teste de \(H_0 : \pi_{+1}-\pi_{1+} = 0\) vs. \(H_1 : \pi_{+1}-\pi_{1+}\neq 0\) usando a função mcnemar.test() do pacote stats em R:

mcnemar.test (x = c.table , correct = FALSE )
## 
##  McNemar's Chi-squared test
## 
## data:  c.table
## McNemar's chi-squared = 1, df = 1, p-value = 0.3173


O argumento \(x\) é para a tabela de contingência. O valor argumento correct = FALSE impede que o numerador seja modificado para \((|n_{21}-n_{12}|-1)^2\), que às vezes é usado para ajudar na aproximação da distribuição \(\chi^2_1\) em amostras pequenas. Devido ao grande \(p\)-valor, não há evidências suficientes para indicar que as probabilidades marginais são diferentes.


1.2.7 Tabelas de contingência maiores


Embora a Seção 1.2 se concentre em dois grupos com respostas binárias, há muitos casos em que existem mais de dois grupos, ou seja, há mais de duas linhas em uma tabela de contingência. Por exemplo, além de comparar um medicamento recém-desenvolvido com um placebo em um ensaio clínico, medicamentos aprovados anteriormente também podem ser incluídos. O objetivo usual neste caso é mostrar que o novo medicamento tem um desempenho melhor do que o placebo e um desempenho semelhante ou melhor do que qualquer outro medicamento aprovado anteriormente. Depois que mais linhas são adicionadas a uma tabela de contingência, geralmente é mais fácil realizar a análise com um modelo de regressão binária. Normalmente, os mesmos tipos de perguntas da Seção 1.2 podem ser respondidos com esses modelos de regressão. Além disso, cenários mais complexos podem ser considerados, como incluir variáveis adicionais na análise.

Também pode haver mais de duas categorias de resposta. Combinada com mais de duas linhas, uma tabela de contingência com \(I\) linhas e \(J\) colunas pode ser formada. A amostragem binomial independente não pode mais ser usada como a estrutura de distribuição de probabilidade subjacente. Em vez disso, geralmente usamos distribuições multinomiais ou de Poisson.


1.3 Exercícios