17 Intervalos de confiança e função de verossimilhança

Nesta sessão vamos examinar um pouco mais a teoria de intervalos de confiança. São ilustrados os conceitos de:

Serão utilizados conceitos do método da quantidade pivotal, a propriedade de normalidade assintótica dos estimadores de máxima verossimilhança e a distribuição limite da função deviance.

17.1 Inferência para a distribuição Bernoulli

Os dados abaixo são uma amostra aleatória da distribuição Bernoulli(p).

0 0 0 1 1 0 1 1 1 1 0 1 1 0 1 1 1 1 0 1 1 1 1 1 1

Desejamos obter:

(a)
o gráfico da função de verossimilhança para p com base nestes dados
(b)
o estimador de máxima verossimilhança de p, a informação observada e a informação de Fisher
(c)
um intervalo de confiança de 95% para p baseado na normalidade assintótica de ˆp
(d)
compare o intervalo obtido em (b) com um intervalo de confiança de 95% obtido com base na distribuição limite da função deviance
(e)
a probabilidade de cobertura dos intervalos obtidos em (c) e (d). (O verdadeiro valor de p é 0.8)

Primeiramente vamos entrar com os dados na forma de um vetor.

  > y <- c(0, 0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 0,
  +     1, 1, 1, 1, 1, 1)

(a)
Vamos escrever uma "função em R¨p  ara obter a função de verossimilhança usando a função de densidade da distribuição binomial com argumento log=TRUE pois obter a log-verossimilhança.

  > vero.binom <- function(p, dados) {
  +     n <- length(dados)
  +     x <- sum(dados)
  +     return(dbinom(x, size = n, prob = p, log = TRUE))
  + }

Esta função exige dados do tipo 0 ou 1 da distribuição Bernoulli. Entretanto às vezes temos dados binomiais do tipo n e x (número x de sucessos em n observações). Por exemplo, para os dados acima teríamos n = 25 e x = 18. Vamos então escrever a função acima de forma mais geral de forma a poder utilizar dados disponíveis tanto em um formato quanto em outro.

  > vero.binom <- function(p, dados, n = length(dados), x = sum(dados)) {
  +     return(dbinom(x, size = n, prob = p, log = TRUE))
  + }

Para obter o gráfico da função de verossimilhança de um conjunto de dados cria-se uma sequência de valores para o parâmetro p e calcula-se o respectivo valor da (log)verossimilhança. O gráfico da função é obtido com os valores fixados dos parâmetros no eixo-x e o respectivos valores da função no eixo-y e unindo-se os pontos assim obtidos. No R isto pode ser feito com os comandos abaixo que produzem o gráfico mostrado na Figura 38. Evitamos os valores nos extremos do espaço paramétrico (p = 0 ou p = 1) pois nestes casos a verossimilhaça é zero e portanto a log-verossimilhança retornada por dbinom() é -Inf.

  > p.vals <- seq(0.01, 0.99, l = 99)
  > logvero <- sapply(p.vals, vero.binom, dados = y)
  > plot(p.vals, logvero, type = "l", xlab = "p", ylab = "l(p)")

Note que os três comandos acima podem ser substituídos por um único que produz o mesmo resultado:

  > curve(vero.binom(x, dados = y), from = 0, to = 1)


PIC


Figura 38: Função de verossimilhança para o parâmetro p da distribuição Bernoulli.


(b)
Dos resultados para distribuição Bernoulli sabemos que o estimador de máxima verossimilhança é dado por

    ∑n
ˆp = ---i=1-yi
       n
e que a informação esperada coincide com a informação observada e sendo iguais a:
           n
I(ˆp) =  --------
        ˆp(1 - pˆ)
. Para indicar o estimador de MV o gráfico poderíamos usar arrows() w para obter os valores numéricos para a amostra dada utilizamos os comandos a seguir.
  > p.est <- mean(y)
  > arrows(p.est, vero.binom(p.est, dados = y), p.est, min(logvero),
  +     len = 0.1)
  > text(p.est, min(logvero), p.est, cex = 0.8, pos = 1, offset = 0.3)
  > io <- ie <- length(y)/(p.est * (1 - p.est))
  > io
  [1] 124.0079
  > ie
  [1] 124.0079

(c)
O intervalo de confiança baseado na normalidade assintótica do estimador de máxima verossimilhança é dado por:

(        ∘ -----        ∘ ----)
 pˆ- zα∕2  I (pˆ),pˆ+ zα∕2  I (pˆ)
e para obter o intervalo no R usamos os comandos a seguir.
  > ic1.p <- p.est + qnorm(c(0.025, 0.975)) * sqrt(1/ie)
  > ic1.p
  [1] 0.5439957 0.8960043

(d)
Vamos agora obter e mostrar em um gráfico o intervalo baseado na função deviance. Lembrando que a deviance é definida pela expressão

D (p) = 2{(pˆ) - l(p)},
definimos umaa função dev.binom() para calcular a deviance. Com o comando curve() podemos obter o gráfico de função deviance.
  > dev.binom <- function(p, dados, n = length(dados), x = sum(dados)) {
  +     p.est <- x/n
  +     vero.p.est <- vero.binom(p.est, n = n, x = x)
  +     dev <- 2 * (vero.p.est - vero.binom(p, n = n, x = x))
  +     dev
  + }
  > curve(dev.binom(x, dados = y), 0.35, 0.95, xlab = "p", ylab = "D(p)")
        inf       sup
  0.5283461 0.8686757

A função deviance D(p) tem distribuição assintótica χn-12 e o intervalo de confiança é dado pelos pontos de intersecção entre a função deviance e o valor de quantil da distribuição χ2 para o nível de significância desejado como ilustrado na Figura 39. Nos comandos a seguir primeiro encontramos o ponto de corte para o nível de confiança de 95%. Depois traçamos a linha de corte com abline() e os pontos de corte que definem o intervalo são as raízes de uma função definida como a diferença entre a função deviance e o valor do ponto de corte.

  > L.95 <- qchisq(0.95, df = 1)
  > abline(h = L.95)
  > lim.fc <- function(x) dev.binom(x, dados = y) - L.95
  > ICdev <- c(inf = uniroot(lim.fc, c(0, p.est))$root, sup = uniroot(lim.fc,
  +     c(p.est, 1))$root)
  > ICdev
        inf       sup
  0.5283461 0.8686757
  > arrows(ICdev, L.95, ICdev, 0, len = 0.1)
  > text(ICdev, 0, round(ICdev, dig = 3), cex = 0.8, pos = 1, offset = 0.3)


        inf       sup
  0.5283461 0.8686757

PIC


Figura 39: Função deviance para o parâmetro p da distribuição Bernoulli.


Agora que já vimos as duas formas de obter o IC passo a passo vamos usar os comandos acima para criar uma função geral para encontrar IC para qualquer conjunto de dados e com opções para os dois métodos.

  > ic.binom <- function(dados, n = length(dados), x = sum(dados), nivel = 0.95,
  +     tipo = c("assintotico", "deviance")) {
  +     tipo <- match.arg(tipo)
  +     alfa <- 1 - nivel
  +     p.est <- x/n
  +     if (tipo == "assintotico") {
  +         se.p.est <- sqrt((p.est * (1 - p.est))/n)
  +         ic <- p.est + qnorm(c(alfa/2, 1 - (alfa/2))) * se.p.est
  +     }
  +     if (tipo == "deviance") {
  +         lim.fc <- function(y, ...) dev.binom(y, ...) - qchisq(nivel,
  +             df = 1)
  +         inf <- ifelse(identical(p.est, 0), 0, uniroot(lim.fc, c(0,
  +             p.est), n = n, x = x)$root)
  +         sup <- ifelse(identical(p.est, 1), 1, uniroot(lim.fc, c(p.est,
  +             1), n = n, x = x)$root)
  +         ic <- c(inf, sup)
  +     }
  +     names(ic) <- c("lim.inf", "lim.sup")
  +     ic
  + }

E agora vamos utilizar a função, primeiro com a aproximação assintótica e depois pela deviance. Note que os intervalos são diferentes!

  > ic.binom(dados = y)
    lim.inf   lim.sup
  0.5439957 0.8960043
  > ic.binom(dados = y, tipo = "dev")
    lim.inf   lim.sup
  0.5283461 0.8686757

(e)
O cálculo do intervalo de cobertura consiste em:

1.
simular dados com o valor especificado do parâmetro;
2.
obter o intervalo de confiança;
3.
verificar se o valor está dentro do intervalo
4.
repetir (1) a (3) e verificar a proporção de simulações onde o valor está no intervalo.

Espera-se que a proporção obtida seja o mais próximo possível do nível de confiança definido para o intervalo.

Para isto vamos escrever uma função implementando estes passos e que utiliza internamente ic.binom() definida acima.

  > cobertura.binom <- function(n, p, nsim, ...) {
  +     conta <- 0
  +     for (i in 1:nsim) {
  +         ysim <- rbinom(1, size = n, prob = p)
  +         ic <- ic.binom(n = n, x = ysim, ...)
  +         if (p > ic[1] & p < ic[2])
  +             conta <- conta + 1
  +     }
  +     return(conta/nsim)
  + }

E agora vamos utilizar esta função para cada um dos métodos de obtenção dos intervalos.

  > set.seed(3214)
  > cobertura.binom(n = length(y), p = 0.8, nsim = 1000)
  [1] 0.897
  > set.seed(3214)
  > cobertura.binom(n = length(y), p = 0.8, nsim = 1000, tipo = "dev")
  [1] 0.96

Note que a cobertura do método baseado na deviance é muito mais próxima do nível de 95%, o que pode ser explicado pelo tamanho da amostra. O IC assintótico tende a se aproximar do nível nominal de confiança na medida que aumenta o tamanho da amostra.

17.2 Exercícios

1.
Refaça o ítem (e) do exemplo acima com n = 10, n = 50 e n = 200. Discuta os resultados.
2.
Seja X1,X2,⋅⋅⋅,Xn uma amostra aleatória da distribuição U(0). Encontre uma quantidade pivotal e:
(a)
construa um intervalo de confiança de 90% para θ
(b)
construa um intervalo de confiança de 90% para log θ
(c)
gere uma amostra de tamanho n = 10 da distribuição U(0) com θ = 1 e obtenha o intervalo de confiança de 90% para θ. Verifique se o intervalo cobre o verdadeiro valor de θ.
(d)
verifique se a probabilidade de cobertura do intervalo é consistente com o valor declarado de 90%. Para isto gere 1000 amostras de tamanho n = 10. Calcule intervalos de confiança de 90% para cada uma das amostras geradas e finalmente, obtenha a proporção dos intervalos que cobrem o verdadeiro valor de θ. Espera-se que este valor seja próximo do nível de confiança fixado de 90%.
(e)
repita o item (d) para amostras de tamanho n = 100. Houve alguma mudança na probabilidade de cobertura?

Note que se - in log F(x i; θ) ~ Γ(n, 1) então -2 in log F(x i; θ) ~ χ2n2.

3.
Acredita-se que o número de trens atrasados para uma certa estação de trem por dia segue uma distribuição Poisson(θ), além disso acredita-se que o número de trens atrasados em cada dia seja independente do valor de todos os outros dias. Em 10 dias sucessivos, o número de trens atrasados foi registrado em:

5 0 3 2 1 2 1 1 2 1

Obtenha:

(a)
o gráfico da função de verossimilhança para θ com base nestes dados
(b)
o estimador de máxima verossimilhança de θ, a informação observada e a informação de Fisher
(c)
um intervalo de confiança de 95% para o número médio de trens atrasados por dia baseando-se na normalidade assintótica de ˆθ
(d)
compare o intervalo obtido em (c) com um intervalo de confiança obtido com base na distribuição limite da função deviance
(e)
o estimador de máxima verossimilhança de ϕ, onde ϕ é a probabilidade de que não hajam trens atrasados num particular dia. Construa intervalos de confiança de 95% para ϕ como nos itens (c) e (d).
4.
Encontre intervalos de confiança de 95% para a média de uma distribuição Normal com variância 1 dada a amostra

9.5 10.8 9.3 10.7 10.9 10.5 10.7 9.0 11.0 8.4
10.9 9.8 11.4 10.6 9.2 9.7 8.3 10.8 9.8 9.0

baseando-se:

(a)
na distribuição assintótica de μˆ
(b)
na distribuição limite da função deviance
5.
Acredita-se que a produção de trigo, Xi, da área i é normalmente distribuída com média θzi, onde zi é quantidade (conhecida) de fertilizante utilizado na área. Assumindo que as produções em diferentes áreas são independentes, e que a variância é conhecida e igual a 1, ou seja, Xi ~ N(θzi, 1), para i = 1,⋅⋅⋅,n:
(a)
simule dados sob esta distribuição assumindo que θ = 1.5, e z = (1, 2, 3, 4, 5). Visualize os dados simulados através de um gráfico de (z × x)
(b)
encontre o EMV de θ, ˆθ
(c)
mostre que ˆθ é um estimador não viciado para θ (lembre-se que os valores de zi são constantes)
(d)
obtenha um intervalo de aproximadamente 95% de confiança para θ baseado na distribuição assintótica de ˆθ