Análise de Dados Categóricos

Apêndice- Métodos de verossimilhança

• 1 Introdução
• 1.1 Estimando a probabilidade de sucesso
• 1.2 O papel da verossimilhança
• 2 Verossimilhança
• 2.1 Definição
• 2.2 Exemplos
• 3 Estimativas de máxima verossimilhança
• 3.1 Maximização matemática da função de log-verossimilhança
• 3.2 Maximização computacional da função de log-verossimilhança
• 3.3 Propriedades do estimador de máxima verossimilhança (EMV) para grandes amostras
• 3.4 Variância do EMV
• 4 Funções de parâmetros
• 4.1 Propriedade de invariância dos MLEs
• 4.2 Método Delta para variâncias de funções
• 5 Inferência com estimadores de máxima verossimilhança (MLEs)
• 5.1 Testes para parâmetros
• 5.2 Intervalos de confiança para parâmetros
• 5.3 Testes para modelos
• 6 Referências

1 Introdução

Métodos de verossimilhança são comumente utilizados em análises estatísticas. Em particular, eles constituem a base para a maioria das técnicas de análise de dados categóricos abordadas neste livro. O conteúdo principal deste livro pressupõe que os leitores já estejam familiarizados com métodos baseados em verossimilhança, pelo menos o suficiente para aplicá-los com certa segurança. No entanto, reconhecemos que os leitores podem vir de diversas áreas do conhecimento e talvez não tenham tido contato prévio com esses métodos. Por isso, apresentamos uma breve introdução a essa importante classe de procedimentos.

Sugerimos que os leitores sem experiência prévia com métodos de verossimilhança consultem, ao menos, este apêndice antes de ler o corpo principal do livro. Estabelecemos referências cruzadas entre este apêndice e as seções em que os métodos de verossimilhança são empregados, permitindo que leitores e instrutores consultem o material conforme a necessidade.

O escopo deste apêndice é deliberadamente limitado. A compreensão da teoria e das propriedades assintóticas dos métodos de verossimilhança não é necessária para a leitura deste livro, embora uma noção intuitiva sobre eles seja útil. Assim, as explicações contidas neste apêndice evitam o excesso de formalismo matemático. Para um tratamento mais completo dos métodos de verossimilhança, recomenda-se consultar obras como Casella and Berger (2002) e Severini (2000).

1.1 Modelo e parâmetros

O objetivo de uma análise estatística é aprender algo sobre uma população, o que poderíamos chamar de “verdade”. Para alcançar esse objetivo, coletam-se dados; no entanto, os dados contêm variabilidade (ou “ruído”) que nos impede de enxergar a verdade com clareza. Para extrair a verdade dos dados, é útil começar com alguma noção de como essa verdade se apresenta e saber algo sobre a origem do ruído.

Por exemplo, em vez de apenas observar um gráfico de uma variável resposta em função de uma variável explicativa, pode ser útil supor que a relação deva ser uma linha reta e que os desvios em torno dessa reta sejam independentes e aproximadamente distribuídos segundo uma distribuição normal. Um modelo estatístico é uma estrutura assumida para a verdade e para o ruído, ou seja, é uma suposição fundamentada. As características do modelo são combinadas em uma distribuição de probabilidade — como a normal, a de Bernoulli ou a de Poisson — que visa servir como uma aproximação útil da realidade. Os modelos geralmente contêm parâmetros.

Os parâmetros do modelo relacionam-se à estrutura da verdade e/ou do ruído e representam características desconhecidas ou flexíveis do modelo. Eles representam quantidades populacionais que frequentemente são de interesse direto do pesquisador — como a média de uma distribuição normal ou a probabilidade de sucesso em uma distribuição binomial —, embora nem sempre ou muitas vezes, não nos importamos realmente com o parâmetro de intercepto em diversos problemas de regressão linear.

O objetivo de uma análise estatística é obter informações sobre os parâmetros do modelo ou sobre algumas funções desses parâmetros, como as previsões em uma regressão, que são uma função dos parâmetros de inclinação e intercepto. Funções de parâmetros do modelo também são parâmetros, uma vez que são igualmente quantidades populacionais desconhecidas; portanto, nossa discussão não fará distinção entre parâmetros do modelo e outros parâmetros.

1.2 O papel da verossimilhança

Um modelo estatístico serve para relacionar os dados aos parâmetros. Ainda precisamos encontrar valores para os parâmetros e utilizá-los para obter informações sobre a população. A primeira etapa — encontrar valores para os parâmetros — é chamada de estimação. A segunda etapa — utilizá-los para tirar conclusões — é chamada de inferência.

Existem muitas maneiras de estimar parâmetros a partir de um modelo estatístico, mas aquela adotada quase universalmente, devido à sua qualidade e flexibilidade, é a estimação por máxima verossimilhança (ML). Isso ocorre porque o procedimento é adaptável a praticamente qualquer modelo estatístico, resultando em um processo automático de estimação. Além disso, ele conta com uma variedade de ferramentas associadas que podem ser utilizadas para inferência.

O procedimento e suas ferramentas apresentam pontos fortes e fracos. Discutimos esses aspectos à medida que surgem em diferentes contextos ao longo do texto.

2 Verossimilhança

2.1 Definição

Modelos estatísticos são geralmente descritos em termos de uma função de probabilidade ou de uma função de densidade de probabilidade. Uma função de probabilidade para uma variável aleatória discreta fornece as probabilidades de cada resultado possível. Para uma variável aleatória contínua, a função densidade de probabilidade correspondente é um pouco mais complexa de interpretar, pois pressupõe que as medições sejam feitas com um número infinito de casas decimais.

De modo geral, uma função de densidade de probabilidade descreve as chances relativas de observar valores provenientes de diferentes áreas da distribuição de probabilidade em questão. A conhecida “curva normal” é um exemplo de função de densidade. Tanto nos casos discretos quanto nos contínuos, valores mais elevados de cada função correspondem a intervalos de valores com maior probabilidade de ocorrência. Ambas funções, as densidades e probabilidades conjuntas para uma amostra de observações pode ser interpretada como a probabilidade de observarmos a amostra completa, dada a distribuição e seus parâmetros.

Na prática, não conhecemos os valores dos parâmetros, mas conhecemos os dados. Portanto, não podemos determinar com exatidão a probabilidade ou densidade conjunta da nossa amostra. Podemos calcular essa quantidade se assumirmos determinados valores para os parâmetros. Se alterarmos os valores dos parâmetros, obteremos um valor diferente para a densidade e ou probabilidade da amostra, pois a probabilidade de ocorrência dos dados varia conforme os parâmetros.

Por exemplo, é muito improvável observarmos 5 sucessos em 10 tentativas de Bernoulli quando a probabilidade real de sucesso é 0.01. Esse resultado seria um pouco mais provável se a probabilidade de sucesso fosse 0.30, e ainda mais provável se fosse 0.50.

Os símbolos reais usados para os parâmetros em um determinado problema variam dependendo do contexto do problema; veja exemplos mais tarde. Essa é a natureza da função de verossimilhança: consideramos a função de probabilidade (PMF) ou a função densidade de probabilidade (PDF) de forma invertida. Observamos como a função varia para diferentes valores dos parâmetros, mantendo os dados fixos.

Podemos, então, utilizar isso para avaliar quais valores dos parâmetros resultam em maiores probabilidades relativas de ocorrência da amostra. Formalmente, se definirmos a PMF ou PDF conjunta de uma amostra como \(f(\pmb{y}|\theta)\) — em que \(\pmb{y} = (y_1, \cdots, y_n)\) representa um vetor contendo os \(n\) valores amostrados e \(\pmb{\theta} = (\theta_1, \cdots, \theta_p)\) representa um vetor de \(p\) parâmetros distintos, sendo que os símbolos específicos utilizados para os parâmetros variam conforme o contexto do problema; veja exemplos mais adiante —, então a função de verossimilhança é \[ L(\pmb{\theta}|\pmb{y}) = f(\pmb{y}|\pmb{\theta})\cdot \]

Valores maiores da função de verossimilhança correspondem a valores dos parâmetros que são relativamente mais bem suportados pelos dados.

Uma função de verossimilhança não é uma probabilidade, porque a única parte aleatória, y, é considerada fixa em sua construção. Em particular, não se espera que a soma da verossimilhança seja igual a 1 para todos os valores de θ. Os valores numéricos reais de uma função de verossimilhança são irrelevantes. O que importa são as magnitudes relativas das funções de verossimilhança para diferentes valores de parâmetros.

Quando as observações são extraídas independentemente, a função de verossimilhança é simplesmente o produto das PDFs ou PMFs. \[ L(\pmb{\theta}|\pmb{y}) = \prod_{i=1}^n f(y_i | \pmb{\theta}), \] onde utilizamos o símbolo \(\displaystyle\prod\) para denotar a multiplicação de termos indexados.

Assim, é muito fácil construir funções de verossimilhança para amostras aleatórias simples, cenário que abrange a maioria dos problemas deste livro. Além disso, observe que o valor da verossimilhança depende da amostra; portanto, as verossimilhanças — e quaisquer características calculadas a partir delas — são estatísticas. Isso significa que tais características são aleatórias e possuem distribuições de probabilidade.

2.2 Exemplos

A seguir, alguns exemplos simples de verossimilhanças comumente usadas na análise de dados categóricos.

Exemplo: Bernoulli.

Suponha que a variável aleatória \(Y\) assuma apenas dois valores possíveis. A função de probabilidade Bernoulli para \(Y\) é \[ f(y | \theta) = \pi^y (1-\pi)^{1-y}, \] com parâmetro de probabilidade de sucesso \(\pi\) (\(0< \pi< 1\)) e \(y=1\) ou 0, denotando “sucesso” ou uma “falha”, respectivamente.

Sejam \(y_1,\cdots,y_n\) observações de variáveis aleatórias de Bernoulli independentes com essa função de probabilidade. A verossimilhança para o parâmetro \(\pi\) é \[ \tag{A.1} L(\pi|\pmb{y})=\prod_{i=1}^n \pi^{y_i}(1-\pi)^{1-y_i}=\pi^\omega(1-\pi)^{n-\omega}, \] onde \(\pmb{y}=(y_1,\cdots,y_n)\) e \(\omega=\sum_{i=1}^n y_i\). Essa função de verossimilhança é utilizada na Seção 1.1.1.

Exemplo: Binomial.

Uma forma alternativa do caso de Bernoulli ocorre quando o número total de sucessos, \(\omega\), é observado em vez dos resultados individuais das tentativas. Nesse caso, a função de probabilidade conjunta de \(y_1,\cdots, y_n\) não pode ser determinada, pois não sabemos quais \(y_i\) devem ser iguais a 1 e quais devem ser iguais a 0.

No entanto, sabemos que existem \(\binom{n}{\omega} = \frac{n!}{\omega!(n-\omega)!}\) maneiras de observar os \(\omega\) sucessos entre as \(n\) observações; portanto, a função de probabilidade de \(\omega\) dado \(\pi\) é \[ f(\omega|\pi)=\dfrac{n!}{\omega! (n-\omega)!} \pi^\omega (1-\pi)^{n-\omega}\cdot \]

Se for realizado apenas um conjunto de \(n\) tentativas, de modo que se observe apenas um número total de sucessos \(\omega\), então \(L(\pi|\omega) = f(\omega|\pi)\). Note que isso é muito semelhante à verossimilhança de Bernoulli.

Exemplo: Poisson.

A função de probabilidade Poisson para uma variável aleatória \(Y\) é \(f(y|\mu) = e^\mu \mu^y / y!\), com parâmetro \(\mu > 0\) e \(y\) assumindo valores inteiros \(0, 1, 2,\cdots\), como na contagem de algo.

Sejam \(y_1,\cdots,y_n\) observações de variáveis aleatórias de Poisson independentes. A função de verossimilhança para o parâmetro \(\mu\) é \[ \tag{A.2} L(\mu|\pmb{y})=\prod_{i=1}^n \dfrac{e^\mu \mu^{y_i}}{y_i!}, \] onde \(\pmb{y}=(y_1,\cdots,y_n)\). Essa função de verossimilhança é utilizada na Seção 4.1.2.

Exemplo: Multinomial.

Considere a variável aleatória \(Y\) com respostas que consistem em uma de \(c\) categorias, rotuladas como \(1,2,\cdots,c\), com respectivas probabilidades de sucesso \(\pi_1,\pi_2,\cdots,\pi_c\), tais que, \[ \sum_{k=1}^c \pi_k =1\cdot \]

Sejam \(y_1,\cdots,y_c\) observações de \(Y\) obtidas em ensaios independentes desse tipo; significa que os valores possíveis para cada \(y_k\) são as categorias \(1,2,\cdots,c\). Funções de verossimilhança podem ser construídas de forma semelhante aos casos de Bernoulli e binomial apresentados anteriormente, dependendo de se são observados os resultados individuais dos ensaios \((y_1,\cdots,y_n)\) ou apenas as contagens resumidas para cada categoria \((\omega_1,\cdots,\omega_c)\).

A distribuição multinomial baseia-se nas contagens resumidas, para as quais a função de probabilidade é \[ \tag{3} f(\omega_1,\cdots,\omega_c|\pi_1,\cdots,\pi_c)=\dfrac{n!}{\omega_1! \times \cdots\times \omega_c!}\pi^{\omega_1}\times \cdots\pi_c^{\omega_c} \]

Se for realizado apenas um conjunto de \(n\) tentativas, de modo que se observe apenas um conjunto de contagens de categorias \(\omega_1,\cdots,\omega_c\), então temos \[ L(\pi_1,\cdots,\pi_c | \omega_1,\cdots,\omega_c) = f(\omega_1,\cdots,\omega_c | \pi_1,\cdots,\pi_c)\cdot \] Essa verossimilhança é utilizada na Seção 3.1.

3 Estimativas de máxima verossimilhança

Definimos a estimativa de máxima verossimilhança (EMV), \(\widehat{\theta}\), de um parâmetro como o valor do parâmetro para o qual a função de verossimilhança amostral dada é maximizada: \[ L(\widehat{\pmb{\theta}}|\pmb{y})\geq L(\widetilde{\pmb{\theta}}|\pmb{y}) \] para qualquer valor possível \(\widetilde{\pmb{\theta}}\) do parâmetro. O exemplo a seguir ilustra como encontrar essa EMV por meio da simples avaliação de uma função de verossimilhança.

Exemplo: Estimativa de Máxima Verossimilhança para uma amostra de Bernoulli

Suponha que sejam observados \(\sum_{i=1}^4 y_i = 4\) sucessos em \(n = 10\) tentativas independentes Bernoulli. Com base nessa informação, queremos determinar o valor mais plausível para \(\pi\).

A função de verossimilhança é \(L(\pi|\pmb{y}) = \pi^4 (1-\pi)^6\). A Tabela 1 apresenta a função de verossimilhança avaliada em alguns valores distintos de \(\pi\) e a Figura 1 exibe o gráfico da função.

# Instalar se necessário: install.packages("ggplot2")
library(ggplot2)

# 1. Definição dos parâmetros básicos
w <- 4
n <- 10

# 2. Dados para a curva contínua (0 a 1)
df_curva <- data.frame(pi = seq(0, 1, by = 0.005))
df_curva$Lik <- df_curva$pi^w * (1 - df_curva$pi)^(n - w)

# 3. Dados para os pontos específicos da amostra
df_pontos <- data.frame(pi = c(0.2, 0.3, 0.35, 0.39, 0.4, 0.41, 0.5))
df_pontos$Lik <- df_pontos$pi^w * (1 - df_pontos$pi)^(n - w)

# Instalar se necessário: install.packages("reactable")
library(reactable)

reactable(
  df_pontos,
  columns = list(
    pi = colDef(name = "Parâmetro (π)", align = "center"),
    Lik = colDef(name = "Verossimilhança L(π | y)", align = "center", format = colFormat(digits = 6))
  ),
  bordered = TRUE,
  highlight = TRUE,
  striped = TRUE
)

Tabela 1: \(L(\pi|\pmb{y})\) avaliado em diferentes valores de \(\pi\).

Observamos que a função de verossimilhança atinge seu valor máximo quando \(\pi = 0.4\). Portanto, o valor mais plausível de \(\pi\), dados os dados observados, é 0.4 — o que faz sentido, visto que essa é a proporção observada de sucessos. Formalmente, dizemos que \(0.4\) é a estimativa de máxima verossimilhança (MLE) de \(\pi\) e a denotamos como \(\widehat{\pi} = 0.4\).

# 4. Criação do gráfico com ggplot2
ggplot() +
  # Linha contínua da verossimilhança
  geom_line(data = df_curva, aes(x = pi, y = Lik), 
            color = "#2c3e50", size = 1.2) +
  
  # Linha vertical tracejada indicando o ponto de máximo (MLE)
  geom_vline(xintercept = w/n, linetype = "dashed", 
             color = "#e74c3c", size = 0.8) +
  
  # Pontos específicos destacados em azul
  geom_point(data = df_pontos, aes(x = pi, y = Lik), 
             color = "#3498db", size = 3) +
  
  # Títulos e Expressões Matemáticas nas legendas dos eixos
  labs(
    title = expression(paste("Função de Verossimilhança para ", pi)),
    subtitle = paste("Estimador de Máxima Verossimilhança (MLE) no pico =", w/n),
    x = expression(pi),
    y = expression(paste("L(", pi, " | y)"))
  ) +
  
  # Tema gráfico minimalista e elegante
  theme_minimal(base_size = 14) +
  theme(
    plot.title = element_text(face = "bold", hjust = 0.5),
    plot.subtitle = element_text(color = "dimgray", hjust = 0.5),
    panel.grid.minor = element_blank() # Remove linhas secundárias de grade
  ) +
  
  # Ajuste exato dos limites do eixo X
  scale_x_continuous(limits = c(0, 1), breaks = seq(0, 1, 0.2))

Figura 1: Função de verossimilhança de Bernoulli avaliada em \(\sum_{i=1}^n y_i = 4\) e \(n = 10\).

Por diversas razões matemáticas, mostra-se mais fácil trabalhar com o logaritmo natural da verossimilhança, \(\log\big(L(\widehat{\pmb{\theta}}|\pmb{y})\big)\), ao buscarmos estimadores de máxima verossimilhança (MLEs). Isso não acarreta problemas, pois a transformação logarítmica não altera a ordenação dos valores de verossimilhança para diferentes valores de \(\pmb{\theta}\); portanto, o MLE também maximiza \(\log\big(L({\pmb{\theta}}|\pmb{y})\big)\).

3.1 Maximização matemática da função de log-verossimilhança

Para modelos simples com um único parâmetro, encontrar o valor de \(\theta\) que maximiza a função de log-verossimilhança é algo que se faz facilmente utilizando cálculo. A técnica padrão consiste em derivar a função de log-verossimilhança em relação ao parâmetro, igualar o resultado a zero e resolver a equação para o parâmetro. Esse processo é demonstrado em alguns dos exemplos apresentados anteriormente.

Exemplo: Bernoulli

A partir da Equação (1), obtemos \[ \log\big(L(\pi|y)\big) = \omega\log(\pi)+(n-\omega)\log(1-\pi)\cdot \]

Ao derivar, encontramos \[ \dfrac{\mbox{d}}{\mbox{d}\pi} \log\big(L(\pi|y)\big) = \dfrac{\omega}{\pi}-\dfrac{n-\omega}{1-\pi}\cdot \]

Igualar essa expressão a zero e resolver para \(\pi\) resulta em \[ \widehat{\pi}= \dfrac{\omega}{n} = \dfrac{\displaystyle \sum_{i=1}^n y_i}{n}\cdot \] Assim, o estimador de máxima verossimilhança (EMV) para \(\pi\) é a proporção amostral de sucessos, o que foi demonstrado no exemplo anterior com \(\sum_{i=1}^n y_i = 4\) e \(n = 10\).

Exemplo: Poisson

A partir da Equação (2), obtemos \[ \log\big(L(\mu|\pmb{y}) \big)=-n\mu +\sum_{i=1}^n y_i \log(\mu)-\sum_{i=1}^n \log(y_i!)\cdot \]

Derivando, encontramos \[ \dfrac{\mbox{d}}{\mbox{d}\pi} \log\big(L(\pi|y)\big) = -n+\dfrac{1}{\mu}\sum_{i=1}^n y_i\cdot \]

Igualando isso a 0 e resolvendo para \(\mu\), obtemos \[ \widehat{\mu}=\frac{1}{n}\sum_{i=1}^n y_i\cdot \] Assim, o estimador de máxima verossimilhança para \(\mu\) é a média amostral.

Muitos modelos dependem de mais de um parâmetro. O modelo multinomial da Equação (3) é um exemplo, pois existe um parâmetro de probabilidade diferente para cada categoria. Outros modelos utilizam funções de regressão para descrever a relação entre médias populacionais ou probabilidades e uma ou mais variáveis explicativas.

Nesses modelos, a maximização é realizada exatamente como descrito anteriormente, calculando-se uma derivada separada em relação a cada parâmetro. Igualar cada derivada a zero resulta em um sistema de equações que deve ser resolvido simultaneamente para encontrar a estimativa de máxima verossimilhança (MLE). Teoricamente, isso não apresenta dificuldades.

Na prática, porém, o processo pode ser desafiador quando realizado manualmente, e as equações podem não admitir soluções de forma fechada. Veja, por exemplo, a função de verossimilhança da regressão logística apresentada na Seção 2.2.

3.2 Maximização computacional da função de log-verossimilhança

Quando as equações não possuem solução de forma fechada, elas são resolvidas utilizando-se uma versão aprimorada do método de tentativa e erro. A ideia geral consiste em começar com uma estimativa inicial para o valor do parâmetro, calcular a log-verossimilhança para esse valor e, em seguida, encontrar iterativamente valores de parâmetros que resultem em log-verossimilhanças cada vez maiores, até que não seja possível obter mais nenhuma melhoria.

O aprimoramento da log-verossimilhança pode ser realizado, por exemplo, calculando-se a inclinação da função de log-verossimilhança na estimativa atual e deslocando a estimativa seguinte por uma certa distância na direção que leva a valores maiores de log-verossimilhança.

Alternativamente, pode-se trabalhar com a primeira derivada da função de log-verossimilhança e buscar valores dos parâmetros que a anulem. Existe uma variedade de algoritmos computacionais rápidos e confiáveis para a execução desses procedimentos. Um dos mais amplamente implementados é o algoritmo de Newton-Raphson.

Exemplo: Estimativa de MLE para uma amostra de variáveis aleatórias Bernoulli

Demonstramos o algoritmo de Newton-Raphson em um cenário simples para encontrar a estimativa de máxima verossimilhança (MLE) de \(\pi\) para o exemplo anterior de Bernoulli, no qual \(\sum_{i=1}^{n} y_i = 4\) em \(n = 10\) tentativas.

Vimos anteriormente que \(\widehat{\pi} = \sum_{i=1}^{n} y_i / n = 0.4\). Nesse contexto, o algoritmo mostrado a continuação utiliza a seguinte equação para obter uma estimativa \(\pi^{(i+1)}\) a partir de uma estimativa anterior \(\pi^{(i)}\): \[ \tag{4} \begin{array}{rcl} \pi^{(i+1)} & = & \pi^{(i)}-\dfrac{\displaystyle \left.\dfrac{\mbox{d}}{\mbox{d}\pi} \log\big(L(\pi|y)\big)\right|_{\pi=\pi^{(i)}}}{\displaystyle \left.\dfrac{\mbox{d}^2}{\mbox{d}\pi^2} \log\big(L(\pi|y)\big)\right|_{\pi=\pi^{(i)}}} \\[0.8em] & = & \pi^{(i)}-\dfrac{\displaystyle \dfrac{\omega}{\pi^{(i)}}-\dfrac{n-\omega}{1-\pi^{(i)}}}{\displaystyle -\dfrac{\omega}{\big(\pi^{(i)}\big)^2}-\dfrac{n-\omega}{\big(1-\pi^{(i)} \big)^2}}\cdot \end{array} \]

A Equação (4) deriva de uma expansão em série de Taylor de primeira ordem em torno de \(\pi^{(i)}\)). Suponhmos que queiramos aproximar uma função \(f(x)\) em um ponto \(x_0\). A expansão em série de Taylor de primeira ordem aproxima \(f(x)\) por \(f(x_0) + (x-x_0)f'(x_0)\), em que \(f'(\cdot)\) é a primeira derivada de \(f(\cdot)\) em relação a \(x\).

Para começar a utilizar o algoritmo, escolhemos um valor inicial \(\pi^{(0)}\) que acreditamos estar próximo de \(\widehat{\pi}\) e o substituímos por \(\pi^{(i)}\) na Equação (4) para obter \(\pi^{(1)}\). Se \(\pi^{(1)}\) estiver “suficientemente próximo” de \(\pi^{(0)}\), paramos e utilizamos \(\pi^{(1)}\) como \(\widehat{\pi}\); caso contrário, substituímos \(\pi^{(1)}\) por \(\pi^{(i)}\) na Equação (4) para encontrar \(\pi^{(2)}\). Esse processo continua até que \[ |\pi^{(i+1)}-\pi^{(i)}| < \epsilon \] para algum número pequeno \(\epsilon> 0\), momento em que dizemos que alcançamos a convergência na iteração \(i + 1\).

Com base nos dados observados, suponha que estimemos inicialmente \(\pi^{(0)} = 0.3\) e consideremos que \(\epsilon = 0.0001\) representa um nível de proximidade suficiente. A Tabela 2 apresenta o histórico de iterações, mostrando que a convergência é alcançada após 5 iterações.

De modo geral, os leitores não precisarão implementar diretamente um procedimento de Newton-Raphson como esse. Em vez disso, utilizaremos funções do R que lidam com esses detalhes. Além disso, vale ressaltar que existem outros algoritmos, além do Newton-Raphson, utilizados para obter estimativas de máxima verossimilhança. Muitos deles estão implementados na função optim() do R.

# Purpose: Simple example for Newton-Raphson algorithm               #

  # Data
  sum.y<-4
  n<-10
  w<-sum.y

  # Initialize some values
  epsilon<-0.0001      # Convergence criteria
  pi.hat<-0.3          # Start value
  save.pi.hat<-pi.hat  # Save the results for each iteration here and put pi.hat as first value
  change<-1            # Initialize change for first time through loop

  # Loop to find the MLE (uses Newton-Raphson algorithm)
  while (abs(change) > epsilon) {
    num<-(w-n*pi.hat) / (pi.hat*(1-pi.hat))
    den<- -w/pi.hat^2 - (n-w)/(1-pi.hat)^2
    pi.hat.new<-pi.hat - num/den  # pi^(i+1) = pi^(i) - (1st der log(L))/(2nd der log(L))
    change<-pi.hat.new-pi.hat  # -num/den
    pi.hat<-pi.hat.new  # Same for next time through loop
    save.pi.hat<-c(save.pi.hat, pi.hat)  # Keeps iteration history
  }

  # Print iteration history
  data.frame(iteration = 1:length(save.pi.hat), save.pi.hat)

##   iteration save.pi.hat
## 1         1   0.3000000
## 2         2   0.3840000
## 3         3   0.3997528
## 4         4   0.3999999
## 5         5   0.4000000

Tabela 2: Iterações para o algoritmo de Newton-Raphson.

O código em R acima gerou essa tabela. Mostramos a seguir o código para criar gráficos que ilustram o algoritmo de Newton-Raphson.

  # Log-likelihood function plot
  curve(expr = sum.y*log(x) + (n-sum.y)*log(1-x), from = 0, to = 1,
      xlab = expression(pi), ylab = "Log-likelihood function")
  points(x = save.pi.hat, y = sum.y*log(save.pi.hat) + (n-sum.y)*log(1-save.pi.hat), pch = 1)  # Each pi^(i)

  # Zoomed in Log-likelihood function plot - use when first pi.hat = 0.3
  curve(expr = sum.y*log(x) + (n-sum.y)*log(1-x), from = 0.25, to = 0.5,
        xlab = expression(pi), ylab = "Log-likelihood function")
  points(x = save.pi.hat, y = sum.y*log(save.pi.hat) + (n-sum.y)*log(1-save.pi.hat), pch = 1)  # Each pi^(i)
  segments(x0 = save.pi.hat, y0 = -10, x1 = save.pi.hat, y1 = sum.y*log(save.pi.hat) + (n-sum.y)*log(1-save.pi.hat),
    lty = "dotted")  # Could choose a different y0 to make this more general

  # first derivation of log-likelihood function
  curve(expr = (w-n*x) / (x*(1-x)), from = 0.2, to = 0.6,
      xlab = expression(pi), ylab = "First derivative of log-likelihood function")
  abline(h=0, col = "blue", lty = "dotted")  # We are trying to find the value of pi that makes the first der of L(pi) equal to 0
  points(x = save.pi.hat, y = (w-n*save.pi.hat) / (save.pi.hat*(1-save.pi.hat)), pch = 1)  # Each pi^(i)

  # Linear approximation using 1st order Taylor's series expansion: f(x) = f(x0) + (x-x0)*fp(x0) where x0 is what is being expanded
  #  about and fp() is first derivative of f()
  x0<-0.3 #Start value
  fp.x0<- -( (1-x0)^2*w + x0^2*(n-w) ) / (x0^2*(1-x0)^2)   # Slope of tangent line at x0
  f.x0<-(w-n*x0) / (x0*(1-x0))
  yint<-f.x0 - x0*fp.x0
  abline(a = yint, b = fp.x0, lty = "dashed", col = "red")
  legend(x =0.35, y = 14, bty = "n", legend = c(expression(paste(frac(d, paste(d,pi)), " ", paste("log[L(", pi, "|", bold(y), ")]"))),
    expression(paste("Tangent line at ", pi==0.3)), expression(paste(pi, " estimates"))), lty = c("solid", "dashed", NA), col = c("black", "red"),
    pch = c(NA, NA, 1))

  # Where the tangent line intersects x-axis = 0 corresponds to pi^(i+1)

3.3 Propriedades do estimador de máxima verossimilhança (EMV) para grandes amostras

Para utilizar um EMV na construção de intervalos de confiança e na realização de testes, precisamos conhecer sua distribuição de probabilidade, a distribuição de probabilidade de uma estatística também é conhecida como sua “distribuição amostral”.

É possível demonstrar que todos os EMVs que utilizaremos compartilham certas propriedades relacionadas às suas distribuições amostrais, tornando-os bases muito atraentes para procedimentos de inferência. Essas propriedades geralmente são válidas para grandes amostras; em outras palavras, elas valem assintoticamente, ou seja, à medida que o tamanho da amostra tende ao infinito (\(\infty\)).

A seguir, apresenta-se uma lista dessas propriedades:

1. Os estimadores de máxima verossimilhança têm distribuição normal assintótica – Esse resultado é análogo ao Teorema do Limite Central para médias amostrais. O fato de a normalidade ocorrer assintoticamente significa que, em qualquer amostra específica, a distribuição normal é tipicamente uma aproximação da distribuição amostral correta de \(\widehat{\theta}\), e essa aproximação melhora à medida que o tamanho da amostra aumenta.

2. Os MLEs são consistentes – Isso significa, essencialmente, que se você amostrar toda a população (ou realizar uma amostragem infinita), o MLE será exatamente igual ao parâmetro populacional. Em particular, qualquer viés na estimativa (a diferença entre o valor esperado do MLE e o valor verdadeiro do parâmetro) desaparece à medida que o tamanho da amostra aumenta, e a variância tende a zero.

3. Os MLEs são assintoticamente eficientes – Isso significa que, à medida que o tamanho da amostra cresce em direção ao infinito, eles alcançam a menor variância possível para estimativas desse tipo (por exemplo, entre todas as estimativas assintoticamente normais). Uma implicação importante desse resultado é que intervalos de confiança baseados em MLEs têm o potencial de ser mais estreitos, e os testes, mais poderosos, do que aqueles baseados em outras formas de estimação.

Não se garante que essas propriedades sejam válidas em amostras que não sejam “grandes”. A aproximação pela distribuição normal é, em geral, muito boa em amostras “grandes”, mas pode ser muito ruim em amostras “pequenas”.

Outros estimadores podem apresentar variância menor do que o estimador de máxima verossimilhança (MLE) em amostras finitas. Infelizmente, não existe uma maneira uniforme de definir o que é “grande” ou “pequeno”. No entanto, muitas vezes não é difícil simular dados a partir do modelo escolhido e verificar se o MLE apresenta uma distribuição aparentemente normal e um viés aceitavelmente pequeno.

3.4 Variância do EMV

A variância de qualquer EMV (Estimador de Máxima Verossimilhança) está relacionada à curvatura da função de log-verossimilhança nas imediações do EMV. Se a log-verossimilhança for muito plana próximo ao máximo, há grande incerteza nos dados quanto à localização do parâmetro, ou seja, muitos valores diferentes de \(\theta\) resultam em valores de verossimilhança igualmente elevados.

Por outro lado, se a log-verossimilhança apresentar um pico acentuado, os dados indicam pouca dúvida sobre a região em que o parâmetro deve estar situado.

Formalmente, a variância do estimador de máxima verossimilhança (MLE) se reduz, em amostras grandes, a uma quantidade que é estimada diretamente a partir da segunda derivada da função de verossimilhança em seu pico: \[ \dfrac{\partial^2}{\partial \pmb{\theta}^2} \log\big( L(\pmb{\theta} \, | \, \pmb{y})\big), \] que é conhecida como matriz Hessiana.

A partir da Hessiana, calcula-se a variância assintótica estimada, que chamaremos simplesmente de “variância estimada”, como \[ \tag{5} \widehat{\mbox{Var}}(\widehat{\pmb{\theta}})=-\left. \mbox{E}\left( \dfrac{\partial^2}{\partial \pmb{\theta}^2} \log\big( L(\pmb{\theta} \, | \, \pmb{Y})\big) \right)^{-1} \right|_{\pmb{\theta}=\widehat{\pmb{\theta}}}\cdot \]

A variância estimada pode, em vez disso, ser aproximada por \[ \tag{6} \widehat{\mbox{Var}}(\widehat{\pmb{\theta}}) \approx -\left. \left( \dfrac{\partial^2}{\partial \pmb{\theta}^2} \log\big( L(\pmb{\theta} \, | \, \pmb{y})\big) \right)^{-1} \right|_{\pmb{\theta}=\widehat{\pmb{\theta}}}\cdot \] o que é frequentemente mais fácil de calcular. A Equação (6) é assintoticamente equivalente à Equação (5), o que significa que essas variâncias estimadas serão essencialmente as mesmas em amostras muito grandes.

Outras variantes assintoticamente equivalentes dessas fórmulas também são utilizadas às vezes. O desvio padrão estimado da estatística, isto é, o erro padrão é calculado extraindo-se a raiz quadrada de \(\widehat{\mbox{Var}}(\widehat{\pmb{\theta}})\).

Qundo \(p>1\), \[ \widehat{\mbox{Var}}(\widehat{\pmb{\theta}})=\begin{pmatrix} \widehat{\mbox{Var}}(\widehat{\theta}_1) & \widehat{\mbox{Cov}}(\widehat{\theta}_1,\widehat{\theta}_2) & \cdots & \widehat{\mbox{Cov}}(\widehat{\theta}_1,\widehat{\theta}_p) \\ \widehat{\mbox{Cov}}(\widehat{\theta}_1,\widehat{\theta}_2) & \widehat{\mbox{Var}}(\widehat{\theta}_2) & \cdots & \widehat{\mbox{Cov}}(\widehat{\theta}_2,\widehat{\theta}_p) \\ \vdots & \vdots & \ddots & \vdots \\ \widehat{\mbox{Cov}}(\widehat{\theta}_1,\widehat{\theta}_p) & \widehat{\mbox{Cov}}(\widehat{\theta}_2,\widehat{\theta}_p) & \cdots & \widehat{\mbox{Var}}(\widehat{\theta}_p) \end{pmatrix}\cdot \]

Essa é conhecida como matriz de variância-covariância estimada, às vezes abreviada como matriz de covariância ou matriz de variância. Os elementos da diagonal da matriz, onde os números da linha e da coluna são iguais, são as variâncias estimadas dos MLEs.

Os elementos fora da diagonal são as covariâncias entre pares de MLEs. Essas covariâncias estimadas medem a dependência entre os MLEs e podem ser úteis para encontrar as variâncias de funções dos MLEs (Seção 4). Porque \(\widehat{\mbox{Cov}}(\widehat{\theta}_i,\widehat{\theta}_j)=\widehat{\mbox{Cov}}(\widehat{\theta}_j,\widehat{\theta}_i)\), a matriz é simétrica, de modo que o elemento na linha \(i\), coluna \(j\) é igual ao elemento na linha \(j\), coluna \(i\), para qualquer \(i\neq j\).

Exemplo: Poisson

O objetivo deste exemplo é encontrar a variância de \(\widehat{\mu}\) na função de probabilidade Poisson. A Figura 2 mostra a curvatura da função de log-verossimilhança para duas amostras diferentes, onde a amostra 1 tem \[ y = (3, 5, 6, 6, 7, 10, 13, 15, 18, 22) \] e a amostra 2 tem \[ y = (9, 12)\cdot \]

Para ambas as amostras, \(\widehat{\mu} = 10.5\). Podemos ver que a função de log-verossimilhança da amostra 2 é relativamente plana devido ao seu pequeno tamanho de amostra, enquanto a função de log-verossimilhança da amostra 1 tem muito mais curvatura devido ao seu maior tamanho de amostra.

Como a variância de \(\widehat{\mu}\) é baseada nessa curvatura, esperaríamos que a variância da amostra 1 fosse muito menor que a variância da amostra 2. Formalmente, calculamos a variância como \[ \widehat{\mbox{Var}}(\widehat{{\mu}}) = -\left. \left( \dfrac{\partial^2}{\partial {\mu}^2} \log\big( L({\mu} \, | \, \pmb{y})\big) \right)^{-1} \right|_{{\mu}=\widehat{{\mu}}} = \dfrac{\widehat{\mu}^2}{\displaystyle \sum_{i=1}^n y_i} = \dfrac{\widehat{\mu}}{n}\cdot \]

#########################
# Amostra #1

  y1<-c(3,5,6,6,7,10,13,15,18,22) 
  n1<-length(y1)  # Sample size
  mle1 <- sum(y1)/n1
  mle1

## [1] 10.5

  maxlik1 <- -n1*mle1 + log(mle1)*sum(y1) - sum(lfactorial(y1))  # log(L) = -n*mu +log(mu)*sum(y_i) - sum(log(y_i !)) evaluated at MLE (mu^)
  maxlik1  # Maximum possible value of log likelihood

## [1] -36.79099

  # curve(expr=-n1*x + log(x)*sumy1 - sumlfac1, from=3, to=25, xlab=expression(mu),ylab="Log Likelihood")
  
  mle1^2/sum(y1)  # Var^(mu^)

## [1] 1.05

#########################
# Amostra #2

  y2 <- c(9,12)
  n2<-length(y2)
  
  # Example of evaluating log(L)
  mu <- c(1:25)
  loglik2 <- -n2*mu + log(mu)*sum(y2) - sum(lfactorial(y2)) # log(L) = -n*mu +log(mu)*sum(y_i) - sum(log(y_i !)) 
  data.frame(mu,loglik2)

##    mu    loglik2
## 1   1 -34.789042
## 2   2 -22.232951
## 3   3 -15.718184
## 4   4 -11.676860
## 5   5  -8.990846
## 6   6  -7.162093
## 7   7  -5.924929
## 8   8  -5.120770
## 9   9  -4.647326
## 10 10  -4.434755
## 11 11  -4.433241
## 12 12  -4.606002
## 13 13  -4.925105
## 14 14  -5.368838
## 15 15  -5.919988
## 16 16  -6.564679
## 17 17  -7.291562
## 18 18  -8.091235
## 19 19  -8.955823
## 20 20  -9.878664
## 21 21 -10.854071
## 22 22 -11.877150
## 23 23 -12.943663
## 24 24 -14.049912
## 25 25 -15.192650

  mle2 <- sum(y2)/n2
  mle2

## [1] 10.5

  maxlik2 <- -n2*mle2 + log(mle2)*sum(y2) - sum(lfactorial(y2)) 
  maxlik2  # Maximum possible value of log likelihood

## [1] -4.410162

  mle2^2/sum(y2)  # Var^(mu^)

## [1] 5.25

#########################
# Plot 

# Instale se não tiver: install.packages("ggplot2")
library(ggplot2)

# Criando um data frame com os valores das funções para o ggplot
x_seq <- seq(3, 25, length.out = 200)
f1 <- -n1*x_seq + log(x_seq)*sum(y1) - sum(lfactorial(y1)) - maxlik1
f2 <- -n2*x_seq + log(x_seq)*sum(y2) - sum(lfactorial(y2)) - maxlik2

df_plot <- data.frame(
  mu = rep(x_seq, 2),
  loglik = c(f1, f2),
  Amostra = rep(c(paste0("Amostra 1 (n=", n1, ")"), paste0("Amostra 2 (n=", n2, ")")), each = 200)
)

# Construindo o gráfico
ggplot(df_plot, aes(x = mu, y = loglik, color = Amostra, linetype = Amostra)) +
  geom_line(size = 1.2) +
  # Linha vertical do MLE
  geom_vline(xintercept = mle1, linetype = "dotted", color = "gray40", size = 0.8) +
  # Ponto no valor máximo
  annotate("point", x = mle1, y = 0, color = "red", size = 3) +
  annotate("text", x = mle1 + 1.8, y = 2, label = paste("MLE =", mle1), color = "gray30", fontface = "italic") +
  # Estilização de cores modernas
  scale_color_manual(values = c("#1f77b4", "#ff7f0e")) +
  scale_linetype_manual(values = c("solid", "dashed")) +
  # Limites e Rótulos
  ylim(-50, 5) +
  labs(
    title = "Log-Verossimilhança Relativa para Amostras de Poisson",
    subtitle = "Amostra 2 (menor) apresenta maior incerteza (curva mais aberta)",
    x = expression(paste("Parâmetro ", mu)),
    y = expression(log[L](mu) - log[L](hat(mu)))
  ) +
  # Tema minimalista e limpo
  theme_minimal(base_size = 12) +
  theme(
    legend.position = c(0.85, 0.2),
    legend.title = element_blank(),
    panel.grid.minor = element_blank(),
    plot.title = element_text(face = "bold")
  )

Figura 2: Log-verossimilhanças de Poisson para amostras de tamanho \(n = 10\) (amostra 1) e \(n = 2\) (amostra 2) com estimador de máxima verossimilhança (EMV) comum \(\widehat{\mu}= 10.5\). Observe que as duas curvas foram deslocadas verticalmente de modo que os valores de log-verossimilhança sejam iguais a 0 no EMV.

Para a amostra 1, \[ \widehat{\mbox{Var}}(\widehat{{\mu}}) = 10.5/10 = 1.05, \] enquanto para a amostra 2, \[ \widehat{\mbox{Var}}(\widehat{{\mu}})= 10.5/2 = 5.25\cdot \] Como esperado, a variância de \(\widehat{\mu}\) é maior para a amostra 2 do que para a amostra 1.

4 Funções de parâmetros

4.1 Propriedade de invariância dos MLEs

Frequentemente, nosso interesse não se limita aos parâmetros do modelo estimados diretamente por meio dos procedimentos descritos nos Seções 3.1 ou 3.2 deste Apêndice. Em vez disso, podemos querer estimar parâmetros que sejam funções dos parâmetros do modelo.

Como exemplos, temos:

1. a razão de chances (odds ratio) — ver, por exemplo, as Seções 1.2.5, 2.2.3, 3.3.1 e 4.2.3 —, que pode ser expressa como uma função de probabilidades de duas distribuições binomiais, de probabilidades multinomiais ou de médias de distribuições de Poisson; e

2. os valores preditos em qualquer modelo de regressão, que podem ser expressos como funções dos coeficientes de regressão.

A propriedade de invariância dos estimadores de máxima verossimilhança (EMVs) estabelece que, se \(\widehat{\pmb{\theta}}\) é o EMV de \(\pmb{\theta}\) e \(g(\pmb{\theta})\) é uma função contínua de \(\pmb{\theta}\), então \(g(\widehat{\pmb{\theta}})\) é o EMV de \(g(\pmb{\theta})\). Em outras palavras, o EMV de qualquer função contínua de parâmetros é simplesmente a mesma função aplicada aos EMVs dos parâmetros. Uma implicação importante desse resultado é que as propriedades para grandes amostras listadas acima (Seção 3.3) também se aplicam a funções de EMVs, permitindo que todos os procedimentos de inferência descritos na Seção 5, deste Apêndice, sejam utilizados com essas funções.

Vale ressaltar que essas propriedades são, novamente, de natureza aproximada quando se trata de amostras finitas. Os tamanhos amostrais necessários para que elas sejam válidas de forma satisfatória para uma determinada função de parâmetros podem ser semelhantes ou bastante diferentes daqueles exigidos para os próprios parâmetros. Essa questão pode ser verificada por meio de simulação.

4.2 Método Delta para variâncias de funções

O método delta é um procedimento muito geral e útil para estimar a variância de uma função de uma variável aleatória. Apresentamo-lo aqui para estimar a variância de uma função de um estimador de máxima verossimilhança (MLE), a qual não é simplesmente a mesma função da variância do MLE.

Supnha que \(\pmb{\theta}\) consista de \(p\) parâmetros \(\theta_1,\cdots,\theta_p\) e defina \[ \left.g_1'(\widehat{\pmb{\theta}})=\dfrac{\partial}{\partial\theta_1} g(\pmb{\theta})\right|_{\pmb{\theta}=\widehat{\pmb{\theta}}},\left.g_2'(\widehat{\pmb{\theta}})=\dfrac{\partial}{\partial\theta_2} g(\pmb{\theta})\right|_{\pmb{\theta}=\widehat{\pmb{\theta}}},\cdots,\left.g_p'(\widehat{\pmb{\theta}})=\dfrac{\partial}{\partial\theta_p} g(\pmb{\theta})\right|_{\pmb{\theta}=\widehat{\pmb{\theta}}}\cdot \]

Então, a partir de uma aproximação em série de Taylor para \(g(\pmb{\theta})\) centrada em \(\widehat{\pmb{\theta}}\), pode-se mostrar que \[ \tag{7} \widehat{\mbox{Var}}\big(g(\widehat{\pmb{\theta}}) \big)\approx \sum_{i=1}^p \big(g_i'(\widehat{\pmb{\theta}}) \big)^2 \, \widehat{\mbox{Var}}\big(g(\widehat{{\theta}}_i) \big)+\sum_i\sum_{j<i} \big(g_i'(\widehat{\pmb{\theta}}) \big) \big(g_j'(\widehat{\pmb{\theta}}) \big)\widehat{\mbox{Cov}}\big(\widehat{{\theta}}_i,\widehat{{\theta}}_j\big), \] onde $ (g(_i) )$, \(i=1,\cdots,p\) são as variâncias estimadas dos MLEs e \(\widehat{\mbox{Cov}}\big(\widehat{{\theta}}_i,\widehat{{\theta}}_j\big)\), \(i,j=1,\cdots,p\), \(i> j\) as covariâncias estimadas descritas no Apêndice 3.4.

Exemplo: Variância de uma razão de chances

Suponha que existam duas variáveis binomiais independentes: \(W_1\), com \(n_1\) tentativas, probabilidade de sucesso \(\pi_1\) e \(\omega_1\) sucessos observados; e \(W_2\), com \(n_2\) tentativas, probabilidade de sucesso \(\pi_2\) e \(\omega_2\) sucessos observados.

Conforme descrito na Seção 1.2.5, inferência sobre a razão de chances \[ OR=\left(\dfrac{\pi_1}{1-\pi_1}\right)\Big/ \left(\dfrac{\pi_2}{1-\pi_2}\right) \] é feita com base na distribuição amostral de \(\log\big(\widehat{OR}\big)\).

A variância de \(\log\big(\widehat{OR}\big)\) é encontrada utilizando o método delta, e é isso que vamos mostrar a continuação.

Primeiro, observe que temos \(p = 2\) parâmetros aqui, de modo que o vetor \(\pmb{\theta}=(\pi_1,\pi_2)^\top\) e \(g(\pmb{\theta})=\left({\pi_1}/{1-\pi_1}\right)\big/ \left({\pi_2}/{1-\pi_2}\right)\). Então \[ g_1'(\widehat{\pmb{\theta}})=\left.\dfrac{\partial}{\partial \pi_1} \left(\dfrac{\pi_1}{1-\pi_1}\right)\Big/ \left(\dfrac{\pi_2}{1-\pi_2}\right)\right|_{(\pi_1,\pi_2)=(\widehat{\pi}_1,\widehat{\pi}_2)}=\dfrac{1}{\widehat{\pi}_1(1-\widehat{\pi}_1)}, \] e, similarmente, \(g_2'(\widehat{\pmb{\theta}})={1}\big/\big({\widehat{\pi}_2(1-\widehat{\pi}_2)}\big)\), onde \(\widehat{\pi}_i=\omega_i/n_i\).

Além disso, temos que \(\widehat{\mbox{Var}}\big(\widehat{\pi}_i \big)=\widehat{\pi}_i(1-\widehat{\pi}_i)/n_i\) e \(\widehat{\mbox{Cov}}\big(\widehat{{\pi}}_1,\widehat{{\pi}}_2\big)=0\), porque as variáveis aleatórias \(W_1\) e \(W_2\) são independentes. Por isso \[ \begin{array}{rcl} \widehat{\mbox{Var}}\big(g(\widehat{\pmb{\theta}}) \big) & \approx & \left(\dfrac{1}{\widehat{\pi}_1(1-\widehat{\pi}_1)} \right)^2 \dfrac{\widehat{\pi}_1(1-\widehat{\pi}_1)}{n_1}+\left(\dfrac{1}{\widehat{\pi}_2(1-\widehat{\pi}_2)} \right)^2 \dfrac{\widehat{\pi}_2(1-\widehat{\pi}_2)}{n_2}\\[0.8em] & = & \dfrac{1}{\omega_1}+\dfrac{1}{n_1-\omega_1}+\dfrac{1}{\omega_2}+\dfrac{1}{n_2-\omega_2}\cdot \end{array} \]

Esta é a expressão apresentada na Seção 1.2.5.

5 Inferência com estimadores de máxima verossimilhança (MLEs)

Ao longo desta seção, consideramos o problema de realizar inferência sobre um único parâmetro \(\theta\), que pode ser um parâmetro do modelo ou alguma função de parâmetros do modelo. Quando apropriado, mencionamos brevemente extensões para múltiplos parâmetros.

5.1 Testes para parâmetros

Considere as hipóteses \[ \begin{array}{c} H_0 : \theta=\theta_0 \\[0.8em] H_a : \theta\neq \theta_0, \end{array} \] onde \(\theta_0\) é um valor específico de interesse.

Vários procedimentos distintos baseados em princípios de verossimilhança podem ser utilizados para realizar esse teste. Todos eles são procedimentos aproximados, no sentido de que podem não atingir exatamente a taxa de erro do Tipo I (\(\alpha\)) especificada. Essas aproximações são muito precisas em amostras grandes, mas podem apresentar resultados insatisfatórios em amostras pequenas.

Testes de Wald

O teste de Wald (Wald 1943) para um único parâmetro é o procedimento de teste baseado em verossimilhança mais conhecido, pois utiliza as mesmas ideias do teste normal padrão, que faz parte de todos os cursos introdutórios de estatística.

Como o estimador de máxima verossimilhança (EMV) tem distribuição assintoticamente normal com variância estimada conforme a Equação (6) ou a Equação (7), temos que \[ Z_0 = \big(\widehat{\theta}-\theta_0 \big)\big/ \sqrt{\widehat{\mbox{Var}}\big(\widehat{\theta}\big)}\overset{\mathcal{D}}{\sim} N(0,1) \] para uma amostra grande, quando a hipótese nula é verdadeira; onde \(\overset{\mathcal{D}}{\sim}\) significa “aproximadamente distribuído como”.

Portanto, rejeitamos \(H_0\) se \[ \tag{8} |Z_0|= \dfrac{|\widehat{\theta}-\theta_0|}{\sqrt{\widehat{\mbox{Var}}\big(\widehat{\theta}\big)}}>Z_{1-\alpha/2}, \] onde onde o valor crítico \(Z_{1-\alpha/2}\) é um quantil \(1-\alpha/2\) de uma distribuição normal padrão.

Alternativamente, calcula-se um \(p\)-valor como \(2P(Z > |Z_0|)\), em que \(Z\) segue uma distribuição normal padrão. Observe que o valor crítico provém da distribuição normal e não de uma distribuição \(t\)-Student, embora estejamos estimando a variância no denominador de \(Z_0\).

Isso ocorre porque a distribuição \(t\)-Student surge especificamente quando a variância no denominador da estatística de teste baseia-se em um cálculo de soma de quadrados de dados provenientes de uma distribuição normal. A variância em \(Z_0\), mencionada anteriormente, baseia-se na Equação (6), a qual, na maioria das vezes, não envolve um cálculo de soma de quadrados.

O teste de Wald é, geralmente, simples de realizar, mas nem sempre é muito eficaz. Em particular, garante-se que o valor crítico \(Z_{1-\alpha/2}\) esteja próximo do valor crítico correto apenas em amostras muito grandes. Em amostras pequenas, a aproximação pela distribuição normal pode ser precária, e \(Z_{1-\alpha/2}\) pode não ser uma aproximação adequada para o verdadeiro valor crítico do teste.

Assim, esse teste é recomendado apenas quando o tamanho da amostra é grande (no contexto do problema) ou quando não é possível utilizar outro teste.

Uma versão do teste de Wald está disponível para testar hipóteses que envolvem mais de um parâmetro. Isso pode ser útil em problemas de regressão com uma variável explicativa categórica, a qual é representada na regressão por várias variáveis indicadoras, cada uma com um parâmetro distinto.

Uma hipótese nula de ausência de associação entre a variável explicativa e a resposta implica que todos esses parâmetros devem ser simultaneamente iguais a zero. Seja \(\pmb{\theta}_0\) o valor hipotético de um parâmetro \(\pmb{\theta}\) de dimensão \(p\). Seja \(\widehat{\pmb{\theta}}\) a estimativa do parâmetro e \(\widehat{\mbox{Var}}\big(\widehat{\theta}\big)\) a sua variância estimada.

Então, testa-se \(H_0 : \pmb{\theta} = \pmb{\theta}_0\) utilizando a estatística de Wald \[ W = \big(\pmb{\theta} - \pmb{\theta}_0\big)^\top \Big( \widehat{\mbox{Var}}\big(\widehat{\pmb{\theta}}\big)\Big)^{-1}\big(\pmb{\theta} - \pmb{\theta}_0\big), \] que possui uma distribuição aproximadamente \(\chi^2_p\) em grandes amostras.

Testes de razão de verossimilhança

Valores de \(\theta\) que apresentam verossimilhanças próximas ao máximo são estimativas mais plausíveis para o valor verdadeiro do parâmetro do que aqueles cujas verossimilhanças são muito inferiores ao máximo.

Assim, comparar o valor da verossimilhança em seu ponto máximo com o melhor valor de verossimilhança obtido quando os parâmetros estão restritos à hipótese nula constitui uma medida de evidência contra a hipótese nula. Acontece que a melhor maneira de realizar essa comparação é por meio da estatística de razão de verossimilhança (LR) \[ \tag{9} \Lambda = \dfrac{L(\theta_0|\pmb{y})}{L\big(\widehat{\theta}|\pmb{y}\big)}\cdot \]

Note que \(\Lambda\leq 1\), pois \(L\big(\widehat{\theta})|\pmb{y}\big)\) é o valor máximo da verossimilhança para os dados fornecidos. O valor de \(\Lambda\) aproxima-se de 1 quando \(L(\theta_0|\pmb{y})\) está próximo de \(L\big(\widehat{\theta})|\pmb{y}\big)\).

Essa proximidade é avaliada com base no fato de que \[ -2\log(\Lambda) = 2\big(L\big(\widehat{\theta}|\pmb{y}\big)-L(\theta_0|\pmb{y})\big) \] segue uma distribuição \(\chi^2_1\) assintótica quando a hipótese nula é verdadeira.

Assim, o teste da razão de verossimilhanças (LRT) para \(H_0 : \theta=\theta_0\) versus \(H_1 : \theta\neq \theta_0\) rejeita \(H_0\) quando \(-2\log(\Lambda) > \chi^2_{1,1-\alpha/2}\), em que \(\chi^2_{1,1-\alpha/2}\) é o quantil \(1-\alpha\) de uma distribuição \(\chi^2\) com um grau de liberdade.

Exemplo: Poisson

O gráfico à esquerda na Figura 3 apresenta a função de log-verossimilhança para o exemplo anterior da distribuição de Poisson com \(n = 10\) (amostra 1). As evidências contra \(\mu = 9\) são relativamente fracas, pois o valor da log-verossimilhança está bastante próximo do pico em \(\widehat{\mu}=10.5\). Por outro lado, há mais evidências contra \(\mu = 5\), cujo valor de log-verossimilhança está muito mais distante do pico.

# Purpose: Compute and plot the likelihood function                  #

#########################             
# Preliminary calculations

  y1<-c(3,5,6,6,7,10,13,15,18,22) 
  n1<-length(y1) #Sample size
  mle1 <- sum(y1)/n1
  mle1

## [1] 10.5

  # log(L) = -n*mu +log(mu)*sum(y_i) - sum(log(y_i !)) evaluated at MLE (mu^)
  maxlik1 <- -n1*mle1 + log(mle1)*sum(y1) - sum(lfactorial(y1)) 
  maxlik1 #Maximum possible value of log likelihood

## [1] -36.79099

  loglik9 <- -n1*9 + log(9)*sum(y1) - sum(lfactorial(y1)) #Suppose mu = 9
  loglik5 <- -n1*5 + log(5)*sum(y1) - sum(lfactorial(y1)) #Suppose mu = 5

  LRT9<--2*loglik9 + 2*maxlik1
  LRT5<--2*loglik5 + 2*maxlik1
  data.frame(LRT9, LRT5)

##       LRT9     LRT5
## 1 2.371643 45.80684

O testes da razão de verossimilhança de \(H_0 : \mu = 9\) vs. \(H_1 : \mu\neq 9\) leva a \[ \begin{array}{rcl} -2\log(\Lambda) & = & -2\log\left( \dfrac{L(\mu=9|\pmb{y})}{L(\mu=10.5|\pmb{y})}\right)\\[0.8em] & = & -2\log\big(L(\mu=9|\pmb{y})\big)+2\log\big(L(\mu=10.5|\pmb{y})\big) \, = \, 2.37\cdot \end{array} \]

Com \(\chi^2_{1,0.05}=3.84\), não rejeitamos \(H_0 : \mu = 9\). De maneira semelhante para \(H_0 : \mu=5\) vs. \(H_a : \mu\neq 5\), calculamos \(-2\log(\Lambda)=45.81\) levando à rejeição da hipótese nula.

O gráfico à direita na Figura 3 apresenta os resultados dos dois testes de hipóteses de uma maneira diferente. Aqui, representamos a região de rejeição determinando o conjunto de todos os valores possíveis de \(\mu\) para os quais \(-2\log(\Lambda) < \chi^2_{1,1-\alpha}\); resolvendo a equação \[ -2\log\big(L(\mu|\pmb{y})\big) + 2\log\big(L(\mu = 10.5|\pmb{y})) = \chi^2_{1,0.95} \] para encontrar as regiões de rejeição. A figura mostra, mais uma vez, que \(H_0 : \mu = 9\) não é rejeitada, mas \(H_0: \mu = 5\) é rejeitada.

par(mfrow=c(1,2))
#########################
# 1. Configurações iniciais do gráfico (Borda limpa e linha principal mais espessa)
curve(expr = -n1*x + log(x)*sum(y1) - sum(lfactorial(y1)), 
      from = 2, to = 27, ylim = c(-94, -34), 
      lwd = 2.5, col = "gray20", bty = "l", las = 1,
      xlab = expression(paste("Parâmetro ", mu)), 
      ylab = "Log-verossimilhança")

# --- Cenário 1: Máxima Verossimilhança (MLE) ---
# Linhas perpendiculares conectando apenas o ponto da curva aos eixos (Grid limpo)
segments(x0 = 2, y0 = maxlik1, x1 = mle1, y1 = maxlik1, lty = "dashed", col = "firebrick", lwd = 1.2)
segments(x0 = mle1, y0 = -94, x1 = mle1, y1 = maxlik1, lty = "dashed", col = "firebrick", lwd = 1.2)

# Textos posicionados de forma inteligente para evitar sobreposição
text(x = mle1, y = -91, labels = expression(hat(mu) == 10.5), col = "firebrick", font = 2)
text(x = 24, y = maxlik1 + 2.5, labels = expression(log[L](hat(mu) * "|" * bold(y))), col = "firebrick")

# --- Cenário 2: Hipótese mu = 9 ---
segments(x0 = 2, y0 = loglik9, x1 = 9, y1 = loglik9, lty = "dotted", col = "dodgerblue4")
segments(x0 = 9, y0 = -94, x1 = 9, y1 = loglik9, lty = "dotted", col = "dodgerblue4")
text(x = 9, y = -85, labels = expression(mu == 9), col = "dodgerblue4")
text(x = 24, y = loglik9 - 2.5, labels = expression(log[L](9 * "|" * bold(y))), col = "dodgerblue4")

# --- Cenário 3: Hipótese mu = 5 ---
segments(x0 = 2, y0 = loglik5, x1 = 5, y1 = loglik5, lty = "dotted", col = "darkorange2")
segments(x0 = 5, y0 = -94, x1 = 5, y1 = loglik5, lty = "dotted", col = "darkorange2")
text(x = 7, y = -81, labels = expression(mu == 5), col = "darkorange2")
text(x = 24, y = loglik5 + 2.5, labels = expression(log[L](5 * "|" * bold(y))), col = "darkorange2")

# Adiciona pequenos pontos nos locais exatos da curva para dar precisão visual
points(x = c(mle1, 9, 5), y = c(maxlik1, loglik9, loglik5), pch = 19, col = c("firebrick", "dodgerblue4", "darkorange2"))
grid()

#########################
# 1. Base do Gráfico (Borda em 'L', linhas mais espessas e números na horizontal)
curve(expr = -n1*x + log(x)*sum(y1) - sum(lfactorial(y1)), 
      from = 2, to = 27, ylim = c(-94, -34), 
      lwd = 2.5, col = "gray20", bty = "l", las = 1,
      xlab = expression(paste("Parâmetro ", mu)), 
      ylab = "Log-verossimilhança")

# --- Cenário 1: Máxima Verossimilhança (MLE) ---
segments(x0 = 2, y0 = maxlik1, x1 = mle1, y1 = maxlik1, lty = "dashed", col = "firebrick", lwd = 1.2)
segments(x0 = mle1, y0 = -94, x1 = mle1, y1 = maxlik1, lty = "dashed", col = "firebrick", lwd = 1.2)
text(x = mle1, y = -91, labels = expression(hat(mu) == 10.5), col = "firebrick", font = 2)
text(x = 10, y = maxlik1 + 2.5, labels = expression(log[L](hat(mu) * "|" * bold(y))), col = "firebrick")
points(x = mle1, y = maxlik1, pch = 19, col = "firebrick")

# --- Cenário 2: Linha de Corte do LRT (0.5 * chi^2) ---
val_corte <- maxlik1 - qchisq(p = 0.95, df = 1)/2

# Linha horizontal que delimita a aceitação (com cor de destaque roxa/indigo)
segments(x0 = 2, y0 = val_corte, x1 = 27, y1 = val_corte, lty = "dotdash", col = "purple4", lwd = 1.5)
# Texto explicativo da linha posicionado de forma limpa à direita
text(x = 20, y = val_corte + 2.5, labels = expression(log[L](hat(mu)) - frac(1,2) * chi[paste(1, ",", 0.95)]^2), col = "purple4", cex = 0.9)

# --- Cenário 3: Regiões de Rejeição (Limites calculados pelo uniroot) ---
op.func <- function(mu, y, alpha, maxlikMLE) {
  n <- length(y)
  
  # 1. Calcula a log-verossimilhança para o mu atual (sob teste)
  loglik_mu <- -n * mu + log(mu) * sum(y) - sum(lfactorial(y))
  
  # 2. Retorna a estatística do LRT subtraída do valor crítico de corte da Chi-quadrado
  # Queremos encontrar onde essa diferença é ZERO (as raízes)
  return(-2 * loglik_mu + 2 * maxlikMLE - qchisq(p = 1 - alpha, df = 1))
}

lower <- uniroot(f = op.func, interval = c(0, mle1), y = y1, alpha = 0.05, maxlikMLE = maxlik1)  # Limite Inferior
upper <- uniroot(f = op.func, interval = c(mle1, 28), y = y1, alpha = 0.05, maxlikMLE = maxlik1) # Limite Superior

# Limite Inferior (Lower Bound)
segments(x0 = lower$root, y0 = -94, x1 = lower$root, y1 = val_corte, lty = "dashed", col = "darkgreen", lwd = 1.2)
points(x = lower$root, y = val_corte, pch = 19, col = "darkgreen")

# Limite Superior (Upper Bound)
segments(x0 = upper$root, y0 = -94, x1 = upper$root, y1 = val_corte, lty = "dashed", col = "darkgreen", lwd = 1.2)
points(x = upper$root, y = val_corte, pch = 19, col = "darkgreen")

# Indicação Visual da Região de Aceitação/Confiança (Espaço entre as linhas)
text(x = (lower$root + upper$root)/2, y = -80, labels = "Região de\nAceitação", col = "gray30", font = 3, cex = 0.9)

# Setas e Textos para as Regiões de Rejeição (Alinhados perfeitamente)
# Região de Rejeição Esquerda
text(x = lower$root - 2.5, y = -88, labels = "Rejeição", col = "darkgreen", cex = 0.85)
arrows(x0 = lower$root - 0.5, y0 = -91, x1 = 2, y1 = -91, length = 0.08, col = "darkgreen", lwd = 1.2)

# Região de Rejeição Direita
text(x = upper$root + 2.5, y = -88, labels = "Rejeição", col = "darkgreen", cex = 0.85)
arrows(x0 = upper$root + 0.5, y0 = -91, x1 = 27, y1 = -91, length = 0.08, col = "darkgreen", lwd = 1.2)
grid()

Figura 3: À esquerda: Log-verossimilhança de Poisson para a amostra de tamanho 10, com três valores de \(\mu\) indicados. À direita: A região de rejeição para um teste de razão de verossimilhanças utilizando \(\alpha = 0.05\). Para qualquer valor de \(\mu\) situado na região de rejeição, a hipótese nula é rejeitada.

O LRT (Teste da Razão de Verossimilhança) generaliza-se para uma ampla gama de problemas envolvendo múltiplos parâmetros e hipóteses que implicam intervalos ou restrições sobre os parâmetros. A abordagem geral consiste em substituir o numerador da Equação (9) pelo valor máximo da função de verossimilhança considerando todos os parâmetros que satisfazem \(H_0\).

Além disso, o denominador é substituído pelo valor máximo da função de verossimilhança sem restrições. Os graus de liberdade para a distribuição \(\chi^2\) são determinados pelo número de restrições impostas aos parâmetros. Isso é explicado mais detalhadamente ao longo do texto, onde os LRTs são utilizados.

Muitos problemas utilizam modelos que contêm parâmetros adicionais não envolvidos na hipótese nula, por exemplo, hipóteses sobre a média de um modelo normal geralmente não impõem restrições à variância. Nesses casos, os parâmetros adicionais não especificados em \(H_0\) — digamos, \(\phi\) — são fixados em suas estimativas de máxima verossimilhança (MLEs) sob as condições impostas a \(\theta\). Ou seja, no denominador da Equação (9), a verossimilhança é maximizada em relação a \(\theta\) e \(\phi\) simultaneamente. No numerador, a verossimilhança é novamente maximizada para ambos os parâmetros simultaneamente, mas sujeita à restrição de que \(\theta\) satisfaça \(H_0\). Assim, os valores de \(\phi\) que produzem as melhores verossimilhanças no numerador e no denominador podem ser diferentes.

O LRT muitas vezes não é simples de realizar manualmente, mas existem técnicas computacionais muito eficazes capazes de maximizar verossimilhanças sob restrições. Isso torna os LRTs amplamente aplicáveis a uma vasta gama de problemas de teste. Além disso, a precisão do valor crítico do LRT é, em geral, muito superior à do valor crítico do teste de Wald para um dado tamanho de amostra; por isso, o LRT é geralmente preferido em relação aos testes de Wald quando ambos estão disponíveis.

Teste escore

Uma abordagem alternativa para testar hipóteses utilizando verossimilhanças consiste em examinar as propriedades da função de verossimilhança na hipótese nula. Como mostra a Figura 4, a inclinação da log-verossimilhança próxima ao pico deve ter uma magnitude menor do que a inclinação distante do pico. Podemos, portanto, utilizar essa inclinação como uma medida da evidência presente nos dados contra \(H_0\). Essa inclinação, também chamada de score (ou escore), é simplesmente a primeira derivada da log-verossimilhança, avaliada em \(\theta_0\).

# Example: Plot motivating a score test

n1 <- 10 
y1 <- c(3,5,6,6,7,10,13,15,18,22) 
lfac1 <- lfactorial(y1) 
sumlfac1 <- sum(lfac1) 
sumy1 <- sum(y1) 
mu <- c(1:25) 
loglik1 <- -n1*mu + log(mu)*sumy1 - sumlfac1 
mle1 <- sumy1/n1  
loglikmax <- -n1*mle1 + log(mle1)*sumy1 - sumlfac1 
loglik9 <- -n1*9 + log(9)*sumy1 - sumlfac1 
loglik5 <- -n1*5 + log(5)*sumy1 - sumlfac1  

# Ajuste da janela e margens
par(mar = c(5, 5, 2, 2))

# Gráfico principal com cor e linha grossa
curve(expr = -n1*x + log(x)*sumy1 - sumlfac1, xlim = c(3, 25), ylim = c(-94,-34), 
      xlab = expression(mu), ylab = "Log Likelihood", 
      cex.lab = 1.4, cex.axis = 1.3, lwd = 3, col = "#1f78b4", frame.plot = TRUE)

# Grade de fundo leve
grid(col = "lightgray", lty = "dotted")

# Linhas pontilhadas verticais coloridas
lines(x = c(mle1,mle1), y = c(loglikmax,-94), lty = "dotted", lwd = 2, col = "#33a02c") 
lines(x = c(9,9), y = c(loglik9,-94), lty = "dotted", lwd = 2, col = "#e31a1c") 
lines(x = c(5,5), y = c(loglik5,-94), lty = "dotted", lwd = 2, col = "#ff7f00") 

# Linhas tracejadas de intervalo
lines(x = c(2,8), y = c(loglik5-3*(-n1+(sumy1/5)), loglik5+3*(-n1+(sumy1/5))), lty = "dashed", lwd = 2, col = "#ff7f00") 
lines(x = c(6,12), y = c(loglik9-3*(-n1+(sumy1/9)), loglik9+3*(-n1+(sumy1/9))), lty = "dashed", lwd = 2, col = "#e31a1c") 
lines(x = c(7.5,13.5), y = c(loglikmax,loglikmax), lty = "dashed", lwd = 2, col = "#33a02c") 

# Textos com ajustes de posição
text(x = 5, y = -92, expression(mu==5), cex = 1.3, col = "#ff7f00", font = 2) 
text(x = 9, y = -92, expression(mu==9), cex = 1.3, col = "#e31a1c", font = 2) 
text(x = mle1 + 0.8, y = -92, expression(mu==hat(mu)), cex = 1.3, col = "#33a02c", font = 2)

Figura 4: Log-verossimilhança de Poisson para uma amostra artificial de tamanho 10, mostrando a função score (inclinação) em três valores de \(\theta\). O score é igual a 0 no estimador de máxima verossimilhança (MLE). A magnitude é maior para o valor de \(\theta\) mais distante do MLE (\(\theta = 5\)) do que para aquele mais próximo (\(\theta = 9\)).

Em particular, para amostras aleatórias independentes, o escore é \[ U_0=\left.\dfrac{\partial}{\partial \theta} \log\big( L(\theta|\pmb{y})\big)\right|_{\theta=\theta_0}\cdot \]

O Teorema do Limite Central garante que \(U_0\) tenha uma distribuição normal assintótica. Sob a hipótese nula, a inclinação média em todos os conjuntos de dados possíveis é zero: \(\mbox{E}(U_0) = 0\). Pode-se demonstrar que a variância assintótica do score — que mede a variabilidade das inclinações sob \(H_0\), a partir de funções de verossimilhança calculadas com base em diferentes conjuntos de dados — é estimada de forma análoga à Equação (5).

Especificamente \[ \widehat{\mbox{Var}}(U_0)=-\left.\mbox{E}\left( \dfrac{\partial^2}{\partial \theta^2} \log\big( L(\theta|\pmb{y})\big)\right)\right|_{\theta=\theta_0}\cdot \] Em seguida, realiza-se um teste escore semelhante ao teste de Wald, comparando-se os resultados de \(U_0/\sqrt{\widehat{\mbox{Var}}(U_0)}\) e \(Z_{1-\alpha/2}\). Extensões para múltiplos parâmetros são realizadas exatamente como no teste de Wald.

O teste de escore também se baseia em propriedades assintóticas; portanto, o valor crítico é uma aproximação para qualquer amostra finita. Em geral, ele apresenta um desempenho superior ao do teste de Wald, mas não necessariamente melhor do que o do teste da razão de verossimilhança (LRT). Sua principal vantagem é que ele utiliza a função de verossimilhança apenas na hipótese nula.

Em alguns problemas complexos, a hipótese nula representa uma simplificação considerável do modelo geral — por exemplo, ao fixar certos parâmetros em zero —, tornando os cálculos muito mais simples de realizar em \(\theta_0\) do que em qualquer outro ponto.

5.2 Intervalos de confiança para parâmetros

Assim como o teste de Wald, o intervalo de confiança de Wald baseia-se em relações familiares que utilizam a distribuição normal. Podemos escrever \[ \big(\widehat{\theta}-\theta_0 \big)\big/ \sqrt{\widehat{\mbox{Var}}\big(\widehat{\theta}\big)}\overset{\mathcal{D}}{\sim} N(0,1), \] onde \(\widehat{\theta}\) e \(\widehat{\mbox{Var}}\big(\widehat{\theta}\big)\) são os mesmos definidos na Seção 3 deste Apêndice. Assim, \[ P\left( Z_{\alpha/2}\leq \big(\widehat{\theta}-\theta_0 \big)\big/ \sqrt{\widehat{\mbox{Var}}\big(\widehat{\theta}\big)}\leq Z_{1-\alpha/2}\right)\approx 1-\alpha, \] onde \(Z_{1-\alpha/2}\) é um quantil \(1-\alpha/2\) de uma distribuição normal padrão.

Após reorganizar os termos, obtemos \[ P\left( \widehat{\theta}-Z_{1-\alpha/2}\sqrt{\widehat{\mbox{Var}}\big(\widehat{\theta}\big)} \leq \theta\leq \widehat{\theta}-Z_{\alpha/2}\sqrt{\widehat{\mbox{Var}}\big(\widehat{\theta}\big)} \right)\approx 1-\alpha\cdot \]

Reconhecendo que \(-Z_{\alpha/2} = Z_{1-\alpha/2}\), isso leva à forma familiar de um intervalo de confiança de \((1-\alpha)100\%\) para \(\theta\) como \[ \tag{10} \widehat{\theta}\pm Z_{1-\alpha/2}\sqrt{\widehat{\mbox{Var}}\big(\widehat{\theta}\big)} \]

Uma abordagem alternativa para encontrar um intervalo de confiança consiste em “inverter” um teste. Ou seja, busca-se o conjunto de valores de \(\theta_0\) para os quais a hipótese \(H_0 : \theta = \theta_0\) não é rejeitada. Esse procedimento é realizado para o teste de Wald, a partir da Equação (8), e conduz ao mesmo intervalo da Equação (10).

Razão de verossimilhança

Os intervalos de confiança da razão de verossimilhança são encontrados invertendo-se o teste da razão de verossimilhança (LRT). Um intervalo de confiança de \((1-\alpha)100\%\) para \(\theta\) é o conjunto de todos os valores possíveis de \(\theta\) tais que \[ \tag{11} -2\Big(L(\theta|\pmb{y})\big/ L\big(\widehat{\theta}|\pmb{y}\big)\Big) \leq \chi^2_{1,1-\alpha}\cdot \]

Na Figura 3, trata-se do intervalo entre as duas áreas rotuladas como “Região de Rejeição”. Assim como ocorre com os testes, os intervalos de confiança baseados na razão de verossimilhança tendem a ser mais precisos do que os de Wald, no sentido de apresentarem um nível de confiança real mais próximo do nível declarado \(1 - \alpha\) para um dado tamanho de amostra.

No entanto, raramente existe uma solução em forma fechada; portanto, são necessários procedimentos numéricos iterativos para determinar os pontos extremos do intervalo em que a igualdade da Equação (11) é satisfeita. Esses cálculos podem ser difíceis de realizar em problemas mais complexos, razão pela qual nem todos os pacotes de software os executam.

Em situações com múltiplos parâmetros, como em modelos de regressão, frequentemente determinamos intervalos de confiança para parâmetros individuais. Por exemplo, suponha que um modelo contenha dois parâmetros, \(\theta_1\) e \(\theta_2\), e que desejemos encontrar um intervalo de confiança de razão de verossimilhança (LR) de \((1-\alpha)100\%\) para \(\theta_1\).

Como \(L\big(\theta_1,\theta_2 |\pmb{y}\big)\) varia em função de ambos os parâmetros, precisamos levar \(\theta_2\) em consideração de alguma forma ao calcular a Equação (11). Uma abordagem consiste em fixar \(\theta_2\) em um determinado valor, digamos \(d\), e então encontrar os valores de \(\theta_1\) que satisfazem \[ -2\left(L\big(\theta_1,d | \pmb{y} \big)\big/ L\big(\widehat{\theta}_1(d),d | \pmb{y}\big) \right)\leq \chi^2_{1,1-\alpha}, \] em que \(\widehat{\theta}_1(d)\) é o estimador de máxima verossimilhança (MLE) de \(\theta_1\) quando \(\theta_2 = d\). No entanto, isso pode resultar em um intervalo de confiança diferente para cada valor de \(d\). Como alternativa, podemos definir o valor de \(\theta_2\) como seu MLE para cada valor de \(\theta_1\) considerado.

Ou seja, fixamos o denominador da estatística LR como o máximo global, \(L\big(\widehat{\theta}_1,\widehat{\theta}_2 | \pmb{y} \big)\), e, para cada valor diferente de \(\theta_1\) testado no numerador, definimos \(L(\theta_1,\theta_2 |\pmb{y})\) como o valor máximo que essa função atinge em relação a todos os valores de \(\theta_2\). Se definirmos \(\widetilde{\theta}_2(c)\) como o MLE de \(\theta_2\) ao fixar \(\theta_1 = c\), então o intervalo de confiança LR perfilado é o conjunto de valores de \(c\) que satisfazem \[ -2\left( L\big(c,\widetilde{\theta}_2(c) | \pmb{y} \big) \big/ L\big(\widehat{\theta}_1,\widehat{\theta}_2 | \pmb{y} \big)\right)\leq \chi^2_{1,1-\alpha}\cdot \]

Intervalos de confiança baseados no teste de escore

Intervalos de confiança baseados no teste de escore também são obtidos pela inversão desse teste. No entanto, isso não é tão simples de realizar quanto no caso do teste de Wald. Neste último, o erro padrão no denominador da estatística de teste não se altera ao se examinarem diferentes valores de \(\theta_0\), de modo que o rearranjo da Equação (8) é fácil.

Para o teste de escore, o denominador varia com \(\theta_0\); assim, o rearranjo pode resultar em cálculos matemáticos complexos, a menos que a forma de \(\widehat{\mbox{Var}}(U_0)\) seja relativamente simples. Em geral, é necessário determinar o intervalo por meio de procedimentos numéricos iterativos, semelhantes aos descritos na Seção 3.2, neste Apêndice. Um exemplo em que os cálculos podem ser realizados com relativa facilidade é o intervalo de escore de Wilson, apresentado na Seção 1.1.2.

5.3 Testes para modelos

Muitas formas de regressão são utilizadas na análise de dados categóricos. Na prática, sua aplicação frequentemente exige a comparação de diversos modelos que envolvem diferentes subconjuntos de variáveis explicativas, ou a comparação de modelos com e sem grupos de variáveis explicativas, por exemplo, grupos de variáveis dummy que representam uma variável explicativa categórica.

Quando os parâmetros do modelo são estimados por meio do método de máxima verossimilhança, dispõe-se de técnicas padrão de comparação de modelos baseadas em testes de razão de verossimilhança.

Considere dois modelos: um modelo completo, \(M_1\), composto por um conjunto de \(p_1\) variáveis explicativas, e o modelo reduzido, \(M_0\), contendo um subconjunto próprio de \(p_0\) variáveis explicativas. É importante que o modelo reduzido não contenha nenhuma variável que não esteja também presente no modelo completo.

Comparar \(M_0\) e \(M_1\) é equivalente a um teste de hipóteses, especificando \(M_0\) como a hipótese nula e \(M_1\) como a alternativa. Ajuste ambos os modelos \(M_0\) e \(M_1\) aos dados e sejam \(L_{M_1}\) e \(L_{M_0}\) suas respectivas verossimilhanças maximizadas. Se usarmos \(\Lambda(M_0,M_1)\) para denotar a razão de verossimilhança \(L_{M_0}/L_{M_1}\), então o teste da razão de verossimilhança (LRT) para \(H_0: M_0\) vs. \(H_a: M_1\) rejeita a hipótese nula se \[ \tag{12} -2\log\big(\Lambda(M_0,M_1) \big)>\chi^2_{(p_1-p_0),1-\alpha)}\cdot \]

A rejeição de \(H_0\) significa que pelo menos um dos parâmetros e, portanto, uma das variáveis que compõem a diferença entre \(M_0\) e \(M_1\) é importante para ser incluído em um modelo que já contém \(M_0\). A não rejeição de \(H_0\) sugere que esse modelo mais simples pode ser suficiente e que as variáveis adicionais não contribuem significativamente para o poder explicativo do modelo.

Naturalmente, nunca podemos concluir que a hipótese nula é verdadeira em qualquer teste de hipóteses; portanto, não é possível afirmar que \(M_0\) é um “modelo significativamente melhor” do que \(M_1\).

Deviance

Na regressão linear com erros normalmente distribuídos e variância constante, a soma dos quadrados dos erros (SSE) de qualquer modelo mede, de forma agregada, a proximidade entre as previsões e as respostas observadas.

Em modelos de regressão mais gerais que envolvem diferentes distribuições, a medida comum de ajuste de um modelo é a deviance. A deviance é simplesmente \(D_M = \Lambda(M, M_\text{SAT})\), conforme a Equação (12), em que \(M\) é o modelo que está sendo ajustado e \(M_\text{SAT}\) é o modelo saturado — aquele que ajusta um parâmetro distinto para cada observação.

O modelo saturado fornece uma previsão perfeita para os dados observados, mas frequentemente não é considerado um modelo viável para prever novos dados, pois é improvável que novos dados apresentem os mesmos “erros” aleatórios dos dados atuais. Veja a Figura 5 para um exemplo simples de modelo saturado.

A razão para utilizar um modelo saturado ao calcular deviances é que qualquer outro modelo é, garantidamente, um subconjunto próprio do modelo saturado; portanto, \(\Lambda(M, M_\text{SAT})\) pode sempre ser calculado. Assim, a deviance frequentemente faz parte da saída padrão de muitos procedimentos de modelagem na análise de dados categóricos.

set.seed(2895612)
x = c(1:10)
set1 <- data.frame(x, y = x + rnorm(n = 10))
lin <- lm(formula = y~x, data = set1)

# 1. Plot inicial (ajustando margens e cores mais suaves)
plot(x = set1$x, y = set1$y, type = "n", xlab = "Eixo X", ylab = "Eixo Y", 
     main = "Ajuste Linear", bty = "l", las = 1)

# 2. Grid primeiro para ficar ao fundo do gráfico
grid(lwd = 1, col = "gray90")

# 3. Desenha as linhas e pontos do Modelo Saturado
lines(x = set1$x, y = set1$y, lty = "dotted", lwd = 2, col = "#2B6CB0")
points(x = set1$x, y = set1$y, pch = 16, col = "#2B6CB0", cex = 1.3)

# 4. Desenha a linha de Regressão
abline(lin, lty = "solid", lwd = 2.5, col = "#C53030")

# 5. Legenda posicionada dinamicamente no topo esquerdo (topleft)
legend("topleft", legend = c("Regressão Linear", "Modelo Saturado"), 
       lty = c("solid", "dotted"), col = c("#C53030", "#2B6CB0"), 
       bty = "n", lwd = 2, cex = 0.9)

Figura 5: Comparação de uma regressão linear simples com o modelo saturado.

Se dispusermos das deviances dos dois modelos, \(M_0\) e \(M_1\), que desejamos comparar conforme descrito acima, podemos calcular \(\Lambda(M_0,M_1) = D_{M_0}-D_{M_1}\) para realizar essa comparação e efetuar o teste, como na Equação (12). Por razões descritas no Capítulo 5, a comparação com \(\chi_{(p_1-p_0),1-\alpha}\) nem sempre é feita.

6 Referências

Casella, G., and R. Berger. 2002. Statistical Inference. Duxbury Press.

Severini, T. 2000. Likelihood Methods in Statistics. Oxford University Press.

Wald, A. 1943. “Tests of Statistical Hypotheses Concerning Several Parameters When the Number of Observations Is Large.” Transactions of the American Mathematical Society, no. 54: 426–82.