\(\newcommand{\prob}{\operatorname{P}}\) \(\newcommand{\e}{\operatorname{e}}\)

Bootstrap Paramétrico

Introdução

Um interesse frequente em Estatística é conhecer a distribuição de probabilidades teórica de alguma estatística amostral, ou de características dessa distribuição, como, por exemplo, sua média. desvio-padrão, percentis. O uso de modelos matemáticos dessas distribuições requer muitas vezes considerável manipulação matemática e podem impor suposições bastante severas e muitas vezes intratáveis em aplicações práticas (por exemplo: normalidade, independência, suposições sobre variabilidade, etc.). As simulações computacionais oferecem mais flexibilidade na especificação de distribuições, a qual, por sua vez, permite uma maior flexibilidade na especificação de modelos. Uma técnica de simulação bastante desenvolvida nos últimos anos e muito utilizada é a simulação Markov chain Monte Carlo - MCMC, que permite a especificação arbitrária de distribuições multivariadas complexas, podendo até combinar famílias teóricas de variáveis distintas. Outra técnica de simulação bastante utilizada é o bootstrap que permite construir empiricamente distribuições amostrais de estatísticas convencionais sem suposições que usualmente restringem seus modelos matemáticos teóricos. Ele pode também ser usado construir empiricamente as estatísticas amostrais para novas estatísticas de teste, mesmo para aquelas criadas para um experimento específico.

Assim, o bootstrap é um método intensivo computacionalmente e muito geral, o qual, permite obter estimativas de função amostral de interesse, além de permitir estabelecer sua incerteza. O bootstrap foi proposto por Efron (1979), tornando-se um dos mais utilizados métodos estatísticos. Como tem acontecido freqüentemente com os métodos de Monte Carlo, o bootstrap facilitou a solução de uma série de problemas de estimação, as quais, por outros meios, ou não seriam possíveis ou certamente muito mais difíceis de serem obtidas. Podemos citar dois métodos. Há três formas de bootstrapping, cuja diferença principalestá na maneira de como a população é estimada:

 Não paramétrico (reamostragem)
 Semiparamétrico (adição de ruído)
 Paramétrico (simulação) 

Uma importante diferença entre os procedimentos de bootstrap para métrico e não-paramétrico é que nesse os elementos da amostra bootstrap tomam apenas valores da amostra original, enquanto qe naquele a amostra bootstrap pode tomar qualquer valor do suporte de \(F\).

Maiores referências sobre o tema podem ser encontradas em Beasley and Rodgers (2012), cujo trabalho está disponível em: Bootstrapping and Monte Carlo methods.

Boostrap paramétrico:

Suponha que, dada uma amostra aleatória \(x_1, x_2, \dots x_n\), estima-se \(F\) a partir de uma estimativa paramétrica da função de distribuição distribuição acumulada \(F_X(x| \hat{\boldsymbol{\theta}})\), usando alguma abordagem pelo método dos momentos ou da máxima verossimilhança. O procedimento de bootstrap paramétrico refere-se a selecionar amostras aleatórias de \(\hat{F} = F_X(x| \hat{\boldsymbol{\theta}})\) para calcular a incerteza de uma estatística \(\hat{\theta}\) qualquer. Esse procedimento está formalizado no algorítmo descrito abaixo:

1. Defina \(B\)

2. Selecione \({\mathbf{x}^*}^{1}, {\mathbf{x}^*}^{2}, \dots, {\mathbf{x}^*}^{B}\) amostras bootstrap independentes de \(\hat{F} = F(\mathbf{x}|\hat{\xi})\), em que \({\mathbf{x}^*}^{b} = ({x_1^*}^{b}, {x_2^*}^{b}, \dots, {x_n^*}^{b}), ~ b = 1, 2, \dots, B\) é um vetor de n observações da b-ésima amostra bootstrap.

3. Calcule \(\hat{\theta}^*(b) = T(\hat{F}({\mathbf{x}^*}^b)), ~ b = 1, 2, \dots, B\), as estimativas da estatística de interesse para cada amostra bootstrap.

4. Calcule a incerteza de \(\hat{\theta}\) por meio de:

   1. erro padrão amostral, \(\operatorname{ep}_B(\hat{\theta}^*)\), definido como: \[ \operatorname{ep}_B(\hat{\theta}^*) = \left[ \frac{\sum_{b=1}^B (\hat{\theta}^*(b) - {\bar{\theta}}^*)}{B-1} \right],~ \text{com } {\bar{\theta}}^* = \frac{\sum_{b=1}^B \hat{\theta}^*(b)}{B}. \]
   2. histograma de \(\hat{\theta}^*(b) = T(\hat{F}({\mathbf{x}^*}^b))\).
   3. intervalo de confiança do histograma.

Distribuição amostral de uma mediana:

Fonte: Albert and Rizzo (2012).

Seja \(y_1, y_2, \dots, y_n\) os valores obtidos de uma amostra aleatória de uma população exponencial com função de densidade de probabilidade:

\[ f_Y(y) = \e^{-(y - \theta)}, ~ y \geq \theta. \]

e seja \(\tilde{y}\) a mediana observada de \(\{ y_1, y_2, \dots, y_n \}\). Observa-se apenas \(\tilde{x}\). Usando esse dado disponível é possível construir um intervalo com 90% de confiança para \(\theta\)?

Usando-se do método da quantidade pivotal, toma-se \(W = \tilde{Y} - \theta\). Salienta-se que \(W\) não depende de \(\theta\), a transforamção transladou a variável aleatória para a origem, ou seja, \(W\) tem distribuição exponencial padrão. Encontra-se então \(w_1\) e \(w_2\) tal que \(\prob\{ w_1 < W < w_2\} = \gamma\) ou, de outra maneira, \(\prob\{ (w_1, w_2) \text{ cobrir } W \} = \gamma\). Assim, \(\prob\{ \tilde{Y} - w_1 < \theta < \tilde{Y} - w_2\} = \gamma\).

Para montar esse intervalo é necessário obter a distribuição amostral da mediana \(W\) de amostra de tamanha \(n\), proveniente de população com distribuíção exponencial padrão. A função de densidade de probabilidade de mediana de amostra de tamanho ímpar é conhecida: \[\begin{align} f_W(w) & = \frac{n!}{ \left(\frac{n-1}{2}\right)!\left(\frac{n-1}{2}\right)!} \e^{-w}(1-\e^{-w})^{\frac{n-1}{2}} (\e^{-w})^{\frac{n-1}{2}} \\ & = \frac{1}{\operatorname{Beta}\left( \frac{n+1}{2}, \frac{n+1}{2} \right)} \e^{-w}(1-\e^{-w})^{\frac{n-1}{2}} (\e^{-w})^{\frac{n-1}{2}}, ~ n \text{ ímpar}. \end{align}\]

Como obter os percentis \(w_1\) e \(w_2\)? Uma maneira seria simular a distribuição amostral pelo método Monte Carlo.

Considere uma amostra de 21 elementos escolhidos aleatriamente de uma população exponencialmente distribuída com função de densidade de probabilidade \(f_Y \tag{*}\). Vamos conhecer o gráfico de função de densidade (\(f_W\)) e a esperença e variância de \(W= \tilde{X} - \theta\), que serão aproximadas por meio da função integrate().

      # -------- Mediana amostral de amostra exponencial

options(digits = 7)
n <- 21
taxa <- 1

  # função densidade da mediana (tamanho amostral ímpar)
f.median <- function(x){
  cte <- 1/beta((n + 1)/2, (n + 1)/2)
  cte * pexp(x)^((n-1)/2)*(1 - pexp(x))^((n - 1)/2) * dexp(x)
} 
  # gráfico da densidade
curve(f.median, from = 0.2, to = 1.5)

Vamos calcular a esperança e a variância de \(W\) para \(n = 21\):

  # valor aproximado da esperança da mediana amostral
integrando <- function(x) x*f.median(x)
E.median <- integrate(integrando, lower = 0, upper = Inf)
str(E.median)

## List of 5
##  $ value       : num 0.716
##  $ abs.error   : num 6.31e-05
##  $ subdivisions: int 3
##  $ message     : chr "OK"
##  $ call        : language integrate(f = integrando, lower = 0, upper = Inf)
##  - attr(*, "class")= chr "integrate"

E.median$value  # exato: 0.70286

## [1] 0.7163905

  # valor aproximado da variância da mediana amostral
integrando2 <- function(x) x^2*f.median(x)
E2.median <- integrate(integrando2, lower = 0, upper = Inf)
  # aproximação do 2o. momeno
E2.median$value 

## [1] 0.5618784

(V.median <- E2.median$value - E.median$value ^ 2)  # exato: 0.01978

## [1] 0.04866309

Serão geradas \(B = 1000\) amostras bootstrap de \(W\) para estimarmos os percentis \(w_1 = 0,05\) e \(w_2 = 0,95\), já que prentendemos contruir um intervalo de confiança para a mediana de amostra de tamanho 21 de distribuição exponencial com distribuição \(f_Y\).

        # ------ distribuição empírica da mediana amostral


# valor estimado da esperança e da variância da mediana amostral
set.seed(666)
  # qte. de amostras bootstrap
B <- 10000
n <- 21

  # medianas amostrais boostrap
medianas <- replicate(B, median(rexp(n = n, rate = taxa)))

  # média e variância das medianas amostrais bootstraps
(media <- mean(medianas))

## [1] 0.713914

(variancia <- var(medianas))

## [1] 0.04746184

  # histograma das B amostras bootstrap
hist(medianas, freq = F, xlab = "Valores observados", ylab = "Densidade",
    main = "Histograma das medianas simuladas")
text(0.2, 1.65, "Valores exatos", cex = 0.8, pos = 3, col = "blue")
text(0.2, 1.50, expression(lambda == 1), cex = 0.8, pos = 3, 
     col = "blue")
text(0.2, 1.30, expression(mu == 0.70286), cex = 0.8, pos = 3, 
     col = "blue")
text(0.2, 1.10, expression(sigma^2 == 0.01978), pos = 3, cex = 0.8, 
     col = "blue")
text(1.2, 1.65, "Valores estimados", pos = 3, cex = 0.8)
text(1.2, 1.5, paste("n =", n), pos = 3, cex = 0.8)
text(1.2, 1.30, bquote(bar(x)[med] == .(round(media, 5))), pos = 3, 
     cex = 0.8)
text(1.2, 1.10, bquote(s[med]^2 == .(round(variancia, 5))), pos = 3, 
     cex = 0.8)
  # densidade exata das medianas amostrais bootstrap
curve(f.median, from = 0.2, to = 1.5, col = "blue", add = T)

As estimativas de \(w_1\) e \(w_2\), respectivamente, percentis \(5\%\) e \(95\%\) são:

      # ---- Percentis estimados de W
percentis <- quantile(medianas, prob = c(0.05, 0.95))
percentis

##        5%       95% 
## 0.3966734 1.1133057

Assim, o intervalo com 90% de confiança para o parâmetro \(\theta\) (de \(f_Y\) é dado por \([ \tilde{x} - 1.113; \tilde{x} - 0.397]\), onde \(\tilde{y}\) é a mediana de amostra de tamanho 21 dessa população. Seja por exemplo a amostra aleatória de tamanho 21 relacionada abaixo (oriunda de uma população exponencial com parâmetro \(\theta = 4,5\)):

##  [1] 6.072468 6.630780 4.862047 5.151346 5.808789 4.833432 4.938401 4.913377
##  [9] 4.533592 5.432880 7.361616 6.466409 7.621179 6.659363 5.061442 4.832503
## [17] 4.918104 4.923814 4.719661 6.298154 7.045329

A mediana amostral é \(\tilde{x} = 5.151\). O intervalo com 90% de confiança para \(\theta\) é \([5.151 - 1.113;\, 5.151 - 0.397]\), ou seja \([4.038;\, 4.754]\).

Repita esse procedimento um grande número de vezes: qual é a procentagem de vezes em que o intervalo de confiança cobre o parâmetro \(\theta = 4.5\).

Importante salientar que esse método de construção de intervalo de confiança depende do tamanho amostral. Seria necessário novo procedimento de simulação Monte Carlo, caso o tamanho amostral fosse diferente de $n = 21 $.

Os extremos do intervalo de confiança construído relacionam-se com as estimativas simuladas dos percentis 5% e 95% de \(W\). Há erros nessas estimativas, pois, elas não são percentis exatos (teóricos) dessa distribuição. É possível obter procedimentos bootstrap para obter os erros padrão desses percentis estimados.

Comparando estimadores

Problema dos táxis

Fonte: Albert and Rizzo (2012).

Num certo dia, uma pessoa anota ao acaso cinco números que identificam os táxis de uma ciadade, digamos, \(\{34, 100, 65, 81, 120\}\). Vamos assumir qye os táxis estão numerado de 1 a \(N\), a quantidade total de táxis dessa cidade. Considere também que são observados ao acaso (e sem qualquer viés) \(n\) deles. De que maneira essa pessoa poderia estimar a quantidade total de táxis na cidade (\(N\)) a partir da amostra observada de tamanho \(n\)?

Formalizando o problema, seja \(Y_1, Y_2, \dots, Y_n\) uma sequência de variáveis aleatórias independentes e com distribuição discreta uniforme, com função de probabilidade: \[ f_Y(y)= \frac{1}{N}, ~ \text{para} ~ y = 1, 2, \dots, N. \]

Considere dois estiamdores possíveis de \(N\) (valor desconhecido da população):

1. \(\hat{N}_1\): o maior valor observado na amostra, ou seja:

\[ \hat{N}_1 = \max\{Y_1, Y_2, \dots, Y_n \}. \] (2) \(\hat{N}_2\): o dobro da média amostral, ou seja:

\[ \hat{N}_2 = 2 \bar{y} \]

Qual deles é o melhor estimador da quantidade total (desconhecida) de táxis (\(N\))?

Os desempenhos destes dois estimadores serão comparados por meio de simulação, Os \(n\) números dos táxis observados são escolhidos aleatoriamente a partir de uma distribuição discreta uniforme com a quantidade total de táxis da cidade (\(N\)), calculando-se as estimativas de \(\hat{N}_1\) e \(\hat{N}_1\). Repetindo esse processo de simulação uma quantidade muito grande de vezes (\(B\)), obtém-se as distribuições amostrais empíricas dos dois estimadores. Os dois estimadores são então comparados, examinando várias características de suas respectivas distribuições amostrais. Para facilitar a montagem da simulação, será implementada a função taxi() para gerar uma amostra e para calcular as duas estimativas dessa amostra. A função tem dois argumenro: Ene, a quantidade real de táxis da cidade (\(N\)) e ene, o tamanho amostral (\(n\)). A saída da função taxi são as duas estimativas amostrais, em forma de vetor:

               # -----  Função para simular os táxis observados

taxi <- function(Ene, ene){
  
          # Entradas:
          #   Ene: quantidade populacional de táxis (N)
          #   ene  quantidade observada de táxis  (n)
          # Saída:
          #   vetor com as estimativas N^1 e N^2
  
  observado <- sample(Ene, size = ene, replace = TRUE)
  est1 <- max(observado)
  est2 <- 2 * mean(observado)
  c(N1 = est1, N2 = est2)
}

Supondo que a quantidade verdadeira de táxis da cidade seja \(N = 100\) e que sejam observados \(n = 5\) táxis, Seja pro exemplo uma realização desta função com esses parâmteros:

             # simulação de uma amostra

  # parâmetros de entrada da função
N <- 100
n <- 5

  # números observados dos táxis
estimativa <- taxi(Ene = N, ene = n)
estimativa

##    N1    N2 
##  74.0 104.8

Os valores simulados das duas estimativas são: \(\hat{N}_1 = 74\) e \(\hat{N}_2 = 104.8\). Sera utilizada a função replicate() para simular esse processo de amostragem uma quantidade grande de vezes (\(B = 1000\)). Os resultados das estimativas das \(B\) simulações serão armazenados na matriz EST:

              # Simulação do processo de amostragem

  # qte de amostras 
B = 1000

  # conjunto de B estimativas de N^1 e N^2
EST <- replicate(B, taxi(N, n))
  # dimensão da matriz EST
dim(EST)

## [1]    2 1000

A matriz EST, saída desse processo de simulação, tem 2 linhas, cada uma contendo as estimativas das 1000 simulações de um dos estimadores (\(\hat{N}_1\), na primeira linha e \(\hat{N}_2\), na segunda). Salienta-se que as saídas múltiplas da função replicate() será uma matriz, em que as linhas correspondem aos resultados e as colunas, as réplicas.

Análise do desempenho dos estimadores:

Uma propriedade desejável de um estimador é ele ser não-viciado, ou seja, quando a média das estimativas obtidas é igual ao parãmetro que se quer estimar. Em nosso caso, o estimador será não viciado se \(\operatorname{E}(\hat{N}) = N\). Quando o estimador é viciado é conveniente seu estimar (ou calcular) seu viés: \[ \operatorname{Viés} = \operatorname{E}(\hat{N}) - N \]

O valor esperado do estimador, \(\operatorname{E}(\hat{N})\) é estimado pela média amostral de suas \(B\) estimativas. O vício é então estimado subatraindo-se o verdadeiro valor do parâmetro (\(N\)) dessa média (como estamos simulando, o parâmetro \(N\) é conhecido): \[ \hat{\operatorname{Viés}} = \frac{1}{B} \sum_{j=1}^B \hat{N}^{(j)} - N, \]

onde \(\hat{N}^{(j)}\) é o valor da estimativa na j-ésima réplica. Por sua vez, o erro padrão estimado dessa simulação é obtido pelo desvio-padrão das \(B\) estimativas(\(\{\hat{N^{(j)}} \}\)), dividido pela raiz quadrada do tamanho das amostras: \[ \operatorname{ep}(\hat{\operatorname{Viés}}) = \frac{s(\hat{N})}{\sqrt{B}} \]

Em nosso exemplo de simulação, as estimativas obtidas são armazenadas na matriz EST e o valor verdadeiro da quantidade de táxis da cidade é N = 100. A estimativa do viés e do erro padrão de cada um dos estimadores são obtidos com os comandos abaixo:

             # Inferência sobre os viés dos estimadores

  # estimativas N^1
vies.1 <- round(c(mean(EST["N1", ]) - N, sd(EST["N1", ]) / sqrt(B)), 2)
vies.1

## [1] -15.69   0.42

  # estimativas N^2
vies.2 <- round(c(mean(EST["N2", ]) - N, sd(EST["N2", ]) / sqrt(B)), 2)
vies.2

## [1] -0.11  0.76

Pode-se avaliar os verdadeiros valores dos vieses, determinando seus intervalos de 95% de confiança. A estimativa intervalar para o vício de \(\hat{N}_1\) é \([-15.69 - 1,96 \times 0.42, -15.69 + 1,96 \times 0.42] = [-16.5132, -14.8668]\). Por sua vez, a estimativa intervalar do viés de \(\hat{N}_2\) é \([-0.11 - 1,96 \times 0.76, -0.11 + 1,96 \times 0.76] = [-1.5996, 1.3796]\). O estimador \(\hat{N}_1\) apresenta viès negativo, indicando qie ele subestima consistentemente a verdadeira quantidade de táxis, \(N\). O intervalo com 95% de confiança do viés de \(\hat{N}_2\) inclue o valor zero e não podemos então concluir se ele tem viés positivo ou negativo.

Embora o estimador \(\hat{N}_2\) seja preferível do ponto de vista do viés, não está claro se ele é o melhor estmador. É desejável que um estimador ofereça estimativas próximas ao verdadeiro valor do parâmetro, assim é melhor compará-los com relação à distância média ao parãmetro (em geral, denominado erro absoluto médio). Essa distância é uma medida que combina a exatidão (viés) e a precisão (erro padrão) do estimador, podendo usá-la como uma medida do desempenho desse dois estimadores. O erro absoluto médio é definido como \(D = \operatorname{E}\left(|\hat{N} - N \right)\) que pode ser estimado por:

\[ \hat{D}= \frac{1}{B} \sum_{j=1}^B|\hat{N}^{(j)} - N| \]

            # erro absoluto médio

  # erro absoluto de cada estimativa
erro.abs <- abs(EST - N)

  # distribuições dos erros absolutos de cada estimativa
boxplot(t(erro.abs))

Do gráfico, percebemos que \(\hat{N}_1 = \max\{y_1, y_2, \dots, y_n \}\) tende a ter erros absolutos médios menores que \(\hat{N}_2 = 2\bar{y}\). A distância média dos dois estimadores ao verdaeiro valor, ou seja, seus erros absoluto médios são:

  # erro absoluto médio
apply(t(erro.abs), 2, mean)

##      N1      N2 
## 15.6890 19.6676

De maneira similar, obtemos as estimativas do erro padrão do erro absoluto médio dos dois estimadores. Calculamos o desvio padrão das estimativas dos erros absolutos, dividido pela raiz quadrada da quantidade de amostras simuladas.

  # erro padrão do erro absoluto médio
apply(t(erro.abs), 2, sd) / sqrt(B)

##        N1        N2 
## 0.4214854 0.4390877

Fica claro que o estimador \(\hat{N}_1\) é preferível a \(\hat{N}_2\), pois, apresenta o menor erro absoluto médio.

Probabilidade de cobertura de intervalo de confiança:

Construção de intervalo de confiança para proporção:

Fonte: Albert and Rizzo (2012).

Deseja-se conhecer o parâmetro populacional \(p\), ou seja, a proporção de sucessos (ocorrência de eventos de interesse) nessa população. Suponha que seja obtida um amostra aleatória de tamanho \(n\) dessa população, calculando-se a proporção amostral, \(\hat{p}\). O intervalo com \(\gamma \times 100 \%\) confiança de Wald para \(p\) é um procedimento assintótico bastante utilizado e baseia-se na normalidade assintótica de \(\hat{p}\): \[ \operatorname{IC}_{\text{W}}= \left( \hat{p}- z\sqrt{\frac{\hat{p}(1 - \hat{p})}{n}}, \hat{p}+ z\sqrt{\frac{\hat{p}(1 - \hat{p})}{n}} \right), \]

em que \(z\) é o percentil\(1 -(1-\gamma)/2\) da variável aleatória normal padrão. Como esse procedimento tem nível de confiança \(\gamma\), deseja-se que esse intervalo contenham o parâmetro desconhecido \(p\) com probabilidade de cobretura \(\gamma\), ou seja, para qualquer valor de proporção \(p \in (0,1)\), espera-se que: \[ \tag{*} \prob \{ p \in \operatorname{IC}_{\text{W}} \} \geq \gamma \]

Esse procedimento pode não apresentar a probabilidade de cobertura nominal \(\gamma\), ou seja, há valores de \(p\) e de \(n\) em que a probabilidade de cobertura real \(\tag{*}\) é menor que o nível de confiança nominal \(\gamma\), pois, esse é um procedimento assintótico. Podemos levantar algumas questões de ordem prática inportantes:

• Há regiões do espaço paramétrico em que a probabilidade de cobertura real \(\tag{*}\) está abaixo do nível de confiança nominal \(\gamma\) do procedimento?

• Há regiões específicas do espaço paramétrico em que o intervalo de confiança de Wald tem probabilidade de cobertura baixa?

• Sob a perspectiva da probabilidade de cobertura de intervalos de confiança para a verdadeira proporção \(p\), há métodos de construção alternativos que sejam superiores ao popular intervalo de confiança de Wald.

Experimento de Monte Carlo para estimação da probabilidade de cobertura

Será construída a função wald() para calcular o intervalo de confiança de Wald, tendo como entrada: a quantidade de sucessos y, o tamanho amostral enee o nível de confiança gama. A saída é a linha de uma matriz com duas colunas, em que cada elemento amostral é um dos extremos do intervalo de confiança:

wald <- function(y, ene, gama = 0.95){

  # Entradas:
  #   ene:  tamanho amostral (n)
  #   y:    quantidade observada de sucessos
  #   gama: nível de confiança do IC de Wald
  # Saída:
  #   linha de matriz com limites inferior (li) e superior (ls)
  
  prop <- y/ene
  z <- qnorm(1 -(1 - gama)/2)
  semi <- z * sqrt(prop*(1-prop)/ene)
  li <- prop - semi
  ls <- prop + semi
  cbind(li, ls)
}

Na função, qnorm() calcula o percentil desejado da normal padrão e "prop` é a proporção de sucessos na amostra. Suponha que observamos 5 sucessos em uma amostra de tamanho 20 e desejamos uma estimativa intervalar da verdadeira proporção p, com 95% de confiança:

  # intervalo de 95%
intervalo <- wald(5, 20, 0.95)
intervalo

##             li        ls
## [1,] 0.0602273 0.4397727

is.matrix(intervalo)

## [1] TRUE

Perceba que a função wald() é vetorizada (boa propriedadde!) e podemos então usá-la para encontrar os intervalos de confiança de um vetor de quantidades de sucessos amostrai y. Seja por exemplo, calcular os intervalos com 95% de confiança para os valores de \(y\) iguais a 2, 4, 6 e 8, em uma amostra de tamanho \(n = 20\). Perceba que a saída será uma matriz em que as linhas serão osextremos dos intervalos de confiança para cada proporção amostral correspondente

  # intervalos de confiança para vetor de qtes
quantidades <- c(2, 4, 6, 8)
wald(quantidades, ene = 20, gama = 0.90)

##               li        ls
## [1,] -0.01034014 0.2103401
## [2,]  0.05287982 0.3471202
## [3,]  0.13145266 0.4685473
## [4,]  0.21981531 0.5801847

O fato de wald() ser vetorizada facilita a montagem da simulação. Para tanto, será construída a função cobertura.mc() para estimar por simulação Monte Carlo a probabilidade de cobertura de intervalos de Wald para uma particular proporção \(p\). Os valores de entrada desta função são: pe, a proporção populacional que se quer simular, ene, o tamanho amostral, gama, o nível de confiança \(\gamma\) do intervalo de confiança de wald e iter, a quantidade de intervalos que serão simulados (\(B\)):

            # ------------------------------------------
            # Probabilidade de cobertura - simulação

cobertura.mc <- function(pe, ene, gama, iter = 1E4){
  #   pe:   proporção verdadeira (populacional)
  #   ene:  tamanho amostral (n)
  #   gama: nível de confiança
  #   iter: iterações MonteCarlo (B)
    # sorteia qte de sucesso
  y<- rbinom(iter, ene, pe)
  ic <- wald(y, ene, gama)
  ic
  mean((ic[ , 1] < pe) & (pe < ic[ , 2]))
}

Os valores de entrada função construída, rbinom() simula iter valores de uma variável aleatória binomial com parâmetros prob = pe e tamanho size = ene. Note que a função utiliza como 1E04 como valor default da quantidade de iterações. Os valores simulados são armzenados no vetor y. Os extremos dos iter intervalos de confiança são calculados pela função wald() e armazenados na matriz ic. Por último, a proporção de intervalos que contém a verdadeira proporção \(p\) é calculada pela função mean() aplicada a vetor lógico criado por verificação de pertinência de pe a cada uma dos iterintervalos de confiança de Wald simulados. Deseja-se uma etimativa por simulação da probabilidade de cobertura de intervalo de Wald com 90% de confiança de uma amostra aleatória de uma população com proporção de sucessos \(p= 0,15\), com tamanho \(n = 20\).

 # simulação para p = 0,15 e n = 20 e B = 1E04 iterações
B <- 1E04
cobre <- cobertura.mc(pe = 0.15, ene = 20, gama = 0.90, iter = B)
cobre

## [1] 0.7962

O erro padrão da estimativa Monte Carlo da real probabilidade de cobertura é estimado por:

  # erro padrão da estimativa da probabilidade de cobertura
cobre.ep <- sqrt(cobre*(1 - cobre)/B)
cobre.ep

## [1] 0.00402822

A estimativa pontual Monte Carlo da probabilidade de cobertura desse intervalo é 0.7962, com um erro padrão de 0.004. É visível que a probabilidade de cobertura é menor que o nível de confiança nominal \(\gamma = 0.9\), fato que pode comprometer as inferências realizadas com esse intervalo.

Pretende-se também utilizar simulação Monte Carlo para estudar a relação entre todos os valores de proporção populacional em todo o espaço paramétroco (intervalo \([0,1]\)) e a probabilidade de cobertura efetiva do intervalo de confiança de Wald. A simulação será efetuada pela função cobertura.muitos() que terá como entrada um vetor de de valores da proporção populacional \(p\) e os demais parâmetros de entrada da função cobertura.mc(), teno como saída vetor com as estimativas das probabilidades de cobertura correspondentes.

 # função para estimação vetorizada de probabilidade de cobertura
cobertura.muitas <- function(pe.vetor, ene, gama) 
  sapply(pe.vetor, cobertura.mc, ene, gama)

Com essa função construída, podemos plotar, com a função curve(). as probabilidades e cobertura do Intervalo de Wald construídos com confiança \(\gamma = 0.90\), tamanho amostral \(n = 20\) e para valores de \(p\) entre 0,001 e 0,999.

  # gráfico dos probabilidades de cobertura
n <- 20
confia <- 0.90
curve(cobertura.muitas(x, ene = n, gama = confia), 
      from = 0.001, to = 0.999,
      xlab = "p", ylab = "Probabilidade de Cobertura",
      ylim = c(0.7, 1))
  # título
title(main = "Cobertura de Procedimento de Wald")
  # 2a, linha de título em tamanho menor
mtext(bquote(n==.(n)~", " ~gamma == .(confia)), line = 0.25, 
      font = 2, cex = 0.90)
  # linha indicando o nível d confiança nominal
abline(h = confia, lty = 2, lwd = 2, col = "gray")

A linha horizontal sinaliza o nível nominal de confiança dos intervalos. Verifica-se na figura qie o intervalo de Wald não é verdadeiramente um intervalo com 90% de confiança já que a a probabilidade de cobertura não é uniformemente maior que 0.90 no espaço paramétrico das proporções populacionais, sendo na realidade menor que 0,90 na maior parte do intervalo \([0, 1]\), com valores próximos de zero nas extremidades.

Referências

Albert, Jim, and Maria Rizzo. 2012. R by Example. Springer Science & Business Media.

Beasley, William Howard, and Joseph Lee Rodgers. 2012. “Bootstrapping and Monte Carlo Methods.” In APA Handbook of Research Methods in Psychology, Vol. 2. Research Designs: Quantitative, Qualitative, Neuropsychological, and Biological, edited by H. Cooper, P. M. Camic, D. L. Long, A. T. Panter, D. Rindskopf, and K. J. Sher, 407–425). New York: American Psychological Association. https://doi.org/10.1037/13620-022.

Efron, B. 1979. “Bootstrap Methods: Another Look at the Jackknife.” Ann. Statist. 7 (1): 1–26. https://doi.org/10.1214/aos/1176344552.