7 Mais sobre intervalos de confiança

Nesta aula vamos nos aprofundar um pouco mais na teoria de intervalos de confiança. São ilustrados os conceitos de:

• obtenção de intervalos de confiança pelo método da quantidade pivotal,
• resultados diversos da teoria de verossimilhança,
• intervalos de cobertura.

Voce vai precisar conhecer de 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.

7.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 $\hat{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 probabildade 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 para obter a função de 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 que possamos utilizar dados disponíveis tanto em um quanto em ou outro formato.

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

Agora vamos obter o gráfico da função de verossimilhança para estes dados. Uma forma de fazer isto é criar uma sequência de valores para o parâmetro $p$ e calcular o valor da verossimilhança para cada um deles. Depois fazemos o gráfico dos valores obtidos contra os valores do parâmetro. No R isto pode ser feito com os comandos abaixo que produzem o gráfico mostrado na Figura 14.

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

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)

Figura : 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

$$\hat{p} = \frac{\sum_{i=1}^n y_i}{n}$$

e que a informação esperada coincide com a esperança observada e sendo iguais a:

$$I(\hat{p}) = \frac{n}{\hat{p}(1-\hat{p})}$$

. Para obter os valores numéricos para a amostra dada utilizamos os comandos:

> p.est <- mean(y)
> arrows(p.est, vero.binom(p.est,dados=y), p.est, min(logvero))
> 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:

$$\left(\hat{p} - z_{\alpha/2} \sqrt{I(\hat{p})}\;\; , \;\;\hat{p} + z_{\alpha/2} \sqrt{I(\hat{p})}\right)$$

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 o intervalo baseado na função deviance graficamente. Primeiro vamos escrever uma função para calcular a deviance que vamos chamar de dev.binom, lembrando que a deviance é definida pela expressão:

$$D(p) = 2\{\l (\hat{p}) - l(p)\}.$$

> 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
+ }

E agora vamos fazer o gráfico de forma similar ao que fizemos para função de verossimilhança, definindo uma sequência de valores, calculando as valores da deviance e traçando a curva.

> p.vals <- seq(0.3, 0.95, l=101)
> dev.p <- dev.binom(p.vals, dados=y)
> plot(p.vals, dev.p, typ="l")

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

Agora usando esta função vamos obter o intervalo gráficamente. Para isto definimos o ponto de corte da função usando o fato que a função deviance $D(p)$ tem distribuição assintótica $\chi^2$. 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 o comando abline. Os comandos seguintes consistem em uma forma simples e aproximada para encontrar os pontos onde a linha conta a função, que definem o intervalo de confiança.

> corte <- qchisq(0.95, df=1)
> abline(h=corte)
> dif <- abs(dev.p - corte)
> inf <- ifelse(p.est==0, 0, p.vals[p.vals<p.est][which.min(dif[p.vals<p.est])])
> sup <- ifelse(p.est==1, 1, p.vals[p.vals>p.est][which.min(dif[p.vals>p.est])])
> ic2.p <- c(inf, sup)
> ic2.p
[1] 0.5275 0.8655
> segments(ic2.p, 0, ic2.p, corte)

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"){
+     p.vals <- seq(0,1,l=1001)
+     dev.p <- dev.binom(p.vals, n = n, x = x)
+     corte <- qchisq(nivel, df=1)
+     dif <- abs(dev.p - corte)
+     inf <- ifelse(p.est==0, 0, p.vals[p.vals<p.est][which.min(dif[p.vals<p.est])])
+     sup <- ifelse(p.est==1, 1, p.vals[p.vals>p.est][which.min(dif[p.vals>p.est])])
+     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.528   0.869

(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(123)
> cobertura.binom(n=length(y), p=0.8, nsim=1000)
[1] 0.885
> set.seed(123)
> cobertura.binom(n=length(y), p=0.8, nsim=1000, tipo = "dev")
[1] 0.954

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 a amostra cresce.

7.2 Exercícios

1. Re-faça o ítem (e) do exemplo acima com $n=10$, $n = 50$ e $n=200$. Discuta os resultados.

2. Seja $X_1, X_2, \cdots, X_n$ uma amostra aleatória da distribuição $U(0,\theta)$. Encontre uma quantidade pivotal e:

   (a)

   construa um intervalo de confiança de 90% para $\theta$

   (b)

   construa um intervalo de confiança de 90% para $\log \theta$

   (c)

   gere uma amostra de tamanho $n=10$ da distribuição $U(0,\theta)$ com $\theta=1$ e obtenha o intervalo de confiança de 90% para $\theta$. Verifique se o intervalo cobre o verdadeiro valor de $\theta$.

   (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 $\theta$. 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 $-\sum_i^n \log F(x_i;\theta) \sim \Gamma(n,1)$ então $-2 \sum_i^n \log F(x_i;\theta) \sim \chi^2_{2n}$.

3. Acredita-se que o número de trens atrasados para uma certa estação de trem por dia segue uma distribuição Poisson($\theta$), 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 $\theta$ com base nestes dados

   (b)

   o estimador de máxima verossimilhança de $\theta$, 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 $\hat{\theta}$

   (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 $\phi$, onde $\phi$ é a probabilidade de que não hajam trens atrasados num particular dia. Construa intervalos de confiança de 95% para $\phi$ 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 $\hat{\mu}$

   (b)

   na distribuição limite da função deviance

5. Acredita-se que a produção de trigo, $X_i$, da área $i$ é normalmente distribuída com média $\theta z_i$, onde $z_i$ é 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, $X_i \sim N(\theta z_i, 1)$, para $i=1,\cdots,n$:

   (a)

   simule dados sob esta distribuição assumindo que $\theta=1.5$, e $z=(1,2,3,4,5)$. Visualize os dados simulados através de um gráfico de $(z \times x)$

   (b)

   encontre o EMV de $\theta$, $\hat{\theta}$

   (c)

   mostre que $\hat{\theta}$ é um estimador não viciado para $\theta$ (lembre-se que os valores de $z_i$ são constantes)

   (d)

   obtenha um intervalo de aproximadamente 95% de confiança para $\theta$ baseado na distribuição assintótica de $\hat{\theta}$