Para os exercícios abaixo, formule as hipóteses a serem testadas (em termos do problema e, depois, em termos estatísticos), defina a variável aleatória a ser medida, construa a região crítica, determine o valor-p e conclua sobre o problema. Fique atento se as suposições de cada teste realizado estão satisfeitas ou não.

library(tidyverse)
## Warning: pacote 'tidyverse' foi compilado no R versão 4.4.3
## Warning: pacote 'ggplot2' foi compilado no R versão 4.4.3
## Warning: pacote 'tidyr' foi compilado no R versão 4.4.3
## Warning: pacote 'readr' foi compilado no R versão 4.4.3
## Warning: pacote 'purrr' foi compilado no R versão 4.4.3
## Warning: pacote 'forcats' foi compilado no R versão 4.4.3
## Warning: pacote 'lubridate' foi compilado no R versão 4.4.3
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.1     ✔ stringr   1.5.2
## ✔ ggplot2   3.5.2     ✔ tibble    3.2.1
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.0.4     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(knitr)
## Warning: pacote 'knitr' foi compilado no R versão 4.4.3
library(boot)
library(bootstrap)

Exercício 1

Calcule as estimativas Jackknife para o erro padrão e viés associados ao coeficiente de correlação para os dados das universidades de direito (dados explorados em sala e disponíveis no código R no AVA). Compare com as estimativas Bootstrap das mesmas quantidades.

data(law) # amostra
plot(law[,1],law[,2])

cor(law[,1],law[,2]) #coeficiente de correlação de Pearson
## [1] 0.7763745
# Bootstrap
B = 3200
cor_boot = NULL
set.seed(100)
for (i in 1:B){
  amostra = sample(seq(1,nrow(law),1),nrow(law),replace=TRUE)
  cor_boot[i] = cor(law[amostra,1],law[amostra,2])
}
hist(cor_boot,breaks=20)

epb<-sqrt(var(cor_boot))
epb
## [1] 0.1362092
vicio_boot = mean(cor_boot) - cor(law[,1],law[,2])
vicio_boot
## [1] -0.004667839
# Jacknife
cor_jack = NULL
for (i in 1:nrow(law)) {
  amostra = law[-i, ]  # remove a i-ésima observação
  cor_jack[i] = cor(amostra[,1], amostra[,2])
}

n = nrow(law)

epj = sqrt(((n-1)^2/(n))*var(cor_jack))
epj
## [1] 0.1425186
vicio_jack = (n-1)*(mean(cor_jack) - cor(law[,1],law[,2]))
vicio_jack
## [1] -0.006473623

Comparando as estimativas do Jacknife com as do Bootstrap, podemos observar que tanto os erros padrão quanto os vícios deram próximos um do outro, sem diferenças muito significativas.

Exercício 2

Já é sabido desde os primeiros cursos de Estatística Descritiva que a variância amostral usando n (tamanho amostral) no denominador é um estimador viesado para a variância populacional, e a correção analítica desse viés é substituir n por n−1 no denominador da variância amostral. Simule amostras de diferentes tamanhos (amostras bem pequenas e maiores) da distribuição Normal(0, \(\sigma^2 = 10\)). Para cada amostra, calcule a variância amostral viesada (com n no denominador), a não viesada (com n−1 no denominador), a não viesada (com viés corrigido pelo método Bootstrap) e a não viesada (com viés corrigido pelo método Jackknife). Faça uma tabela com todas essas estimativas e compare os resultados. O que você observa?

tamanhos = c(5, 10, 20, 30, 50, 100, 500, 1000, 5000, 20000)

variancias = data.frame(
  Var_viesada = numeric(),
  Var_nao_viesada = numeric(),
  Var_nao_viesada_boot = numeric(),
  Var_nao_viesada_jack = numeric()
)

B = 3000
set.seed(400)

for (i in tamanhos) {
  n = i
  amostra = rnorm(n, 0, sqrt(10))
  var_n = ((n-1)/n) * var(amostra)     # viesada
  var_n_1 = var(amostra)               # não viesada
  
  boot_vars = numeric(B)
  for (m in 1:B) {
    amostra_boot = sample(seq_len(n), n, replace = TRUE)
    boot_vars[m] = var(amostra[amostra_boot])
  }
  var_vicio_corrigido_boot = 2 * var_n - mean(boot_vars, na.rm = TRUE)
  
  jack_vars = numeric(n)
  for (j in 1:n) {
    amostra_jack = amostra[-j]
    jack_vars[j] = ((n-2)/(n-1)) * var(amostra_jack)
  }
  var_vicio_corrigido_jack = var_n - ((n-1) * (mean(jack_vars, na.rm = TRUE) - var_n))
  
  variancias = rbind(
    variancias,
    data.frame(
      Var_viesada = var_n,
      Var_nao_viesada = var_n_1,
      Var_nao_viesada_boot = var_vicio_corrigido_boot,
      Var_nao_viesada_jack = var_vicio_corrigido_jack
    )
  )
}

rownames(variancias) = paste0("Amostra n = ", tamanhos)

knitr::kable(variancias, digits = 3, caption = "Variâncias Estimadas")
Variâncias Estimadas
Var_viesada Var_nao_viesada Var_nao_viesada_boot Var_nao_viesada_jack
Amostra n = 5 8.880 11.100 8.844 11.100
Amostra n = 10 3.266 3.629 3.291 3.629
Amostra n = 20 5.883 6.192 5.858 6.192
Amostra n = 30 9.652 9.985 9.644 9.985
Amostra n = 50 14.965 15.270 14.926 15.270
Amostra n = 100 10.137 10.239 10.101 10.239
Amostra n = 500 9.882 9.902 9.865 9.902
Amostra n = 1000 10.509 10.519 10.508 10.519
Amostra n = 5000 9.929 9.931 9.928 9.931
Amostra n = 20000 9.901 9.902 9.900 9.902

Podemos observar que conforme o tamanho da amostra aumenta, as estimativas das 4 variâncias se tornam muito próximas (praticamente iguais). Podemos ver também que a variância não viesada (com viés corrigido pelo método Jackknife) coincide com a variância não viesada (com n−1 no denominador).

Exercício 3

A pressão sanguínea de 11 pacientes foi medida antes e depois de uma nova droga ser administrada. As diferenças na pressão sanguínea sistólica (pressão antes – pressão depois) para cada paciente são: 7; 5; 12; -3; -5; 2; 14; 18; 19; 21; -1.

  1. Use testes não paramétricos (permutação e bootstrap) adequados para verificar se a amostra (assumida aleatória) contradiz a hipótese de não mudança sistemática na pressão sanguínea.

Hipóteses:

\(H_0: \mu_D = 0\) (Não há mudança sistemática na pressão sanguínea)

\(H_1: \mu_D \neq 0\) (Há mudança sistemática na pressão sanguínea)

Variável aleatória a ser medida: D = diferença na pressão sistólica (antes − depois) por paciente.

Método de Permutação:

diff = c(7, 5, 12, -3, -5, 2, 14, 18, 19, 21, -1) # pressão antes - pressão depois

boxplot(diff)

est_obs = mean(diff)
est_obs
## [1] 8.090909
media_rearranjos = NULL

# temos 2^11=2048 combinações diferentes com dif negativas e positivas.
for (i1 in c(1,-1)){
  for (i2 in c(1,-1)){
    for (i3 in c(1,-1)){
      for (i4 in c(1,-1)){
        for (i5 in c(1,-1)){    
          for (i6 in c(1,-1)){
            for (i7 in c(1,-1)){
              for (i8 in c(1,-1)){  
                for (i9 in c(1,-1)){
                  for (i10 in c(1,-1)) {
                    for (i11 in c(1,-1)) {
                      media_rearranjos = c(media_rearranjos,
                                      mean(diff*c(i1,i2,i3,i4,i5,i6,i7,i8,i9,i10, i11)))}}}}}}}}}
                    }
                    
}

hist(media_rearranjos)

quantile(media_rearranjos, 0.025) # quantil da região crítica
## 2.5% 
##   -7
quantile(media_rearranjos, 0.975) # simétrico ao primeiro quantil
## 97.5% 
##     7
valor_p<-2*(sum(media_rearranjos>=est_obs)/length(media_rearranjos)) # valor-p
valor_p
## [1] 0.02246094
# Forma alternativa sem usar os fors

# gera todas as combinações de +1/-1 (2048 linhas, 11 colunas)
combis <- expand.grid(rep(list(c(-1,1)), length(diff)))

# transforma em matriz
combis <- as.matrix(combis)

# calcula as médias rearranjadas
media_rearranjos_2 <- rowMeans(combis * rep(diff, each = nrow(combis)))

hist(media_rearranjos_2)

quantile(media_rearranjos_2,0.025)
## 2.5% 
##   -7
quantile(media_rearranjos_2,0.975)
## 97.5% 
##     7
valor_p<-2*(sum(media_rearranjos_2>=est_obs)/length(media_rearranjos_2)) # valor-p
valor_p
## [1] 0.02246094

O valor \(\bar{D}_{obs} = 8.09\) e, considerando \(\alpha = 0.05\), encontramos através da distribuição de permutação que RC = {\(\bar{D}:\) \(\bar{D} \leq -7\) ou \(\bar{D} \geq 7\)} e valor-p = 0.023. Ou seja, com nível de significância 0.05, rejeitamos \(H_0\) (de não mudança sistemática na pressão sanguínea) e temos evidências de que há mudança sistemática na pressão sanguínea, pois a média das diferenças entre a pressão antes a depois é diferente de zero.

Método Bootstrap:

diff = c(7, 5, 12, -3, -5, 2, 14, 18, 19, 21, -1) # pressão antes - pressão depois
est_obs = mean(diff)

B = 50000
boot_means = NULL

set.seed(1000)
for (i in 1:B){
    amostra_boot = sample(seq(1,length(diff),1),length(diff),replace=TRUE)
    boot_means[i] = mean(diff[amostra_boot])
}

LIP = quantile(boot_means, 0.025)
LSP = quantile(boot_means, 0.975)

IC_percentil = data.frame(
    Limite_Inferior = LIP,
    Limite_Superior = LSP
)

kable(IC_percentil, 
      digits = 3, 
      caption = "IC bootstrap para média das diff pelo percentil - 95% de confiança"
      )
IC bootstrap para média das diff pelo percentil - 95% de confiança
Limite_Inferior Limite_Superior
2.5% 2.909 13.273 
# Pelo valor-p:

diff_ajust = diff - est_obs

boot_means_ajust = NULL

for (i in 1:B){
    amostra = sample(seq(1,length(diff_ajust),1),length(diff_ajust),replace=TRUE)
    boot_means_ajust[i] = mean(diff_ajust[amostra])/(sd(diff_ajust[amostra])/sqrt(11))
}

hist(boot_means_ajust, breaks=50)

valor_p = mean(abs(boot_means_ajust) >= abs(est_obs))
valor_p
## [1] 0.00016

Como o IC não contém zero e o valor-p deu muito menor que 0.05, rejeitamos \(H_0\) ao nível de significância de 5%. Ou seja, temos evidências pela amostra obtida que há mudança sistemática na pressão sanguínea sistólica dos pacientes.

  1. Use um teste paramétrico para verificar a mesma hipótese do item acima. Quais suposições são assumidas nesses testes?
diff = c(7, 5, 12, -3, -5, 2, 14, 18, 19, 21, -1)

shapiro.test(diff)
## 
##  Shapiro-Wilk normality test
## 
## data:  diff
## W = 0.93612, p-value = 0.476
# Teste t pareado. Como já temos as diferenças, paired = FALSE.
t.test(diff, mu = 0, alternative = "two.sided", paired = FALSE)
## 
##  One Sample t-test
## 
## data:  diff
## t = 2.8955, df = 10, p-value = 0.01596
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##   1.86476 14.31706
## sample estimates:
## mean of x 
##  8.090909

Portanto, como p-valor deu menor que 0.05, rejeitamos \(H_0\). Ou seja, ao nível de 5% de significância, temos evidências de que há mudança sistemática na pressão sanguínea sistólica dos pacientes.

As suposições assumidas nesses testes é que as diferenças \(d_1, d_2, ..., d_n\) tem distribuição normal com média \(\mu_D\) e variância \(\sigma^2\). Pelo teste de Shapiro-Wilk, não rejeitamos a hipótese nula e temos a suposição de normalidade atendida. Assume-se também que os pares entre pacientes são independentes.

Exercício 4

Um pesquisador está interessado em verificar se uma nova droga é melhor que o placebo para tratar uma certa doença. Devido à natureza da doença, um número limitado de pacientes foi testado. Entre eles, 5 foram aleatoriamente escolhidos para tomar o placebo e 5 para tomar a nova droga. Suponha que a concentração de uma substância química no sangue é medida e valores menores dessa substância indicam melhor condição. Os dados são:

Droga (X): 3.2, 2.1, 2.3, 1.2, 1.5

Placebo (Y): 3.4, 3.5, 4.1, 1.7, 2.1

  1. Conduza o teste de permutação e o bootstrap para checar a eficiência da nova droga e compare os resultados. Quais suposições foram feitas nesse teste?

Hipóteses:

\(H_0: \mu_{droga} = \mu_{placebo}\) (Não há diferença na concentração média da substância química no sangue entre pacientes sob droga e sob placebo)

\(H_1: \mu_{droga} < \mu_{placebo}\) (A droga leva a valores menores da concentração média da substância química no sangue)

Variável aleatória a ser medida: D = \(\bar{X}\) - \(\bar{Y}\) (diferença entre médias)

Método de Permutação:

droga = c(3.2, 2.1, 2.3, 1.2, 1.5)
placebo = c(3.4, 3.5, 4.1, 1.7, 2.1)

amostra_comb = c(droga, placebo)

boxplot(amostra_comb~c(rep(1,5),rep(2,5)))

N = length(amostra_comb)
n = length(droga)
m = length(placebo)

D_obs = mean(droga) - mean(placebo)
D_obs
## [1] -0.9
D = NULL

for (i in 1:(N-n+1)){
   for (j in (i+1):(N-n+2)){
      for (l in (j+1):(N-n+3)){
        for (k in (l+1):(N-n+4)){
          for (v in (k+1):10){
              grupo_1 = amostra_comb[c(i,j,l,k,v)]
              grupo_2 = amostra_comb[-c(i,j,l,k,v)]
              D = c(D, mean(grupo_1) - mean(grupo_2))
        }
      }
    }
  }
}

# Ordenando os valores
sort(D)
##   [1] -1.58 -1.50 -1.50 -1.34 -1.26 -1.14 -1.14 -1.14 -1.06 -1.06 -1.06 -1.02
##  [13] -1.02 -0.98 -0.98 -0.94 -0.90 -0.90 -0.90 -0.90 -0.86 -0.82 -0.82 -0.82
##  [25] -0.82 -0.82 -0.78 -0.78 -0.78 -0.78 -0.78 -0.78 -0.74 -0.74 -0.70 -0.70
##  [37] -0.70 -0.70 -0.70 -0.70 -0.66 -0.66 -0.62 -0.62 -0.62 -0.62 -0.58 -0.58
##  [49] -0.58 -0.58 -0.54 -0.54 -0.54 -0.54 -0.54 -0.50 -0.46 -0.46 -0.46 -0.46
##  [61] -0.46 -0.46 -0.42 -0.42 -0.42 -0.42 -0.38 -0.38 -0.38 -0.38 -0.34 -0.34
##  [73] -0.34 -0.34 -0.34 -0.34 -0.34 -0.34 -0.34 -0.30 -0.30 -0.26 -0.26 -0.26
##  [85] -0.26 -0.26 -0.26 -0.26 -0.22 -0.22 -0.22 -0.22 -0.18 -0.18 -0.18 -0.18
##  [97] -0.18 -0.18 -0.14 -0.14 -0.14 -0.14 -0.14 -0.10 -0.10 -0.10 -0.10 -0.10
## [109] -0.10 -0.10 -0.10 -0.10 -0.10 -0.06 -0.06 -0.06 -0.06 -0.02 -0.02 -0.02
## [121] -0.02 -0.02 -0.02 -0.02 -0.02 -0.02  0.02  0.02  0.02  0.02  0.02  0.02
## [133]  0.02  0.02  0.02  0.06  0.06  0.06  0.06  0.10  0.10  0.10  0.10  0.10
## [145]  0.10  0.10  0.10  0.10  0.10  0.14  0.14  0.14  0.14  0.14  0.18  0.18
## [157]  0.18  0.18  0.18  0.18  0.22  0.22  0.22  0.22  0.26  0.26  0.26  0.26
## [169]  0.26  0.26  0.26  0.30  0.30  0.34  0.34  0.34  0.34  0.34  0.34  0.34
## [181]  0.34  0.34  0.38  0.38  0.38  0.38  0.42  0.42  0.42  0.42  0.46  0.46
## [193]  0.46  0.46  0.46  0.46  0.50  0.54  0.54  0.54  0.54  0.54  0.58  0.58
## [205]  0.58  0.58  0.62  0.62  0.62  0.62  0.66  0.66  0.70  0.70  0.70  0.70
## [217]  0.70  0.70  0.74  0.74  0.78  0.78  0.78  0.78  0.78  0.78  0.82  0.82
## [229]  0.82  0.82  0.82  0.86  0.90  0.90  0.90  0.90  0.94  0.98  0.98  1.02
## [241]  1.02  1.06  1.06  1.06  1.14  1.14  1.14  1.26  1.34  1.50  1.50  1.58
dist_D_H0 = round(table(D)/length(D),3)
dist_D_H0
## D
##               -1.58                -1.5               -1.34               -1.26 
##               0.004               0.008               0.004               0.004 
##               -1.14               -1.06               -1.02               -0.98 
##               0.012               0.012               0.008               0.008 
##               -0.94                -0.9               -0.86               -0.82 
##               0.004               0.016               0.004               0.020 
##               -0.78               -0.74                -0.7               -0.66 
##               0.024               0.008               0.024               0.008 
##               -0.62               -0.58               -0.54                -0.5 
##               0.016               0.016               0.020               0.004 
##               -0.46               -0.42               -0.38               -0.34 
##               0.024               0.016               0.016               0.036 
##                -0.3               -0.26               -0.22               -0.18 
##               0.008               0.028               0.016               0.024 
##               -0.14                -0.1 -0.0600000000000001               -0.02 
##               0.020               0.040               0.016               0.036 
##                0.02  0.0600000000000001                 0.1                0.14 
##               0.036               0.016               0.040               0.020 
##                0.18                0.22                0.26                 0.3 
##               0.024               0.016               0.028               0.008 
##                0.34                0.38                0.42                0.46 
##               0.036               0.016               0.016               0.024 
##                 0.5                0.54                0.58                0.62 
##               0.004               0.020               0.016               0.016 
##                0.66                 0.7                0.74                0.78 
##               0.008               0.024               0.008               0.024 
##                0.82                0.86                 0.9                0.94 
##               0.020               0.004               0.016               0.004 
##                0.98                1.02                1.06                1.14 
##               0.008               0.008               0.012               0.012 
##                1.26                1.34                 1.5                1.58 
##               0.004               0.004               0.008               0.004
barplot(dist_D_H0)

quantil = quantile(D, 0.05) # quantil da região crítica
quantil
##     5% 
## -0.998
val_p = sum(D<=D_obs)/length(D) # valor-p
val_p
## [1] 0.07936508

Método Bootstrap:

droga = c(3.2, 2.1, 2.3, 1.2, 1.5)
placebo = c(3.4, 3.5, 4.1, 1.7, 2.1)

var_comb = sqrt((var(droga)/length(droga)) + (var(placebo)/length(placebo)))

T_obs = (mean(droga) - mean(placebo))/var_comb
T_obs
## [1] -1.575798
amostra_comb = c(droga, placebo)

media_comb = mean(amostra_comb)

droga_transf = droga - mean(droga) + media_comb
placebo_transf = placebo - mean(placebo) + media_comb

B = 5000

T_boot = NULL
set.seed(100)

for (i in 1:B){
    amostra1 = sample(seq(1,length(droga_transf),1),length(droga_transf),replace=TRUE)
    amostra2 = sample(seq(1,length(placebo_transf),1),length(placebo_transf),replace=TRUE)
    
    media_d_transf_boot = mean(droga_transf[amostra1])
    media_p_transf_boot = mean(placebo_transf[amostra2])
    
    var_d_transf_boot = var(droga_transf[amostra1])
    var_p_transf_boot = var(placebo_transf[amostra2])
    
    # Calcular o erro padrão da diferença para a amostra bootstrap
    ep_diff_boot = sqrt(var_d_transf_boot/length(droga) + var_p_transf_boot/length(placebo))
    
    if (ep_diff_boot > 0) {
        T_boot[i] = (media_d_transf_boot - media_p_transf_boot) / ep_diff_boot
    } else {
        T_boot[i] = 0 # Se não houver variabilidade, a diferença padronizada é 0
    }
}

hist(T_boot)

IC = quantile(T_boot,c(0.025,0.975))
IC
##      2.5%     97.5% 
## -2.858482  2.233162
valor_critico = quantile(T_boot, 0.05)

p_value = sum(T_boot <= T_obs) / B
p_value
## [1] 0.0938
  1. Conduza o teste paramétrico nessa situação e compare os resultados. O teste paramétrico é adequado nessa situação? Justifique.
shapiro.test(droga)
## 
##  Shapiro-Wilk normality test
## 
## data:  droga
## W = 0.96086, p-value = 0.814
shapiro.test(placebo)
## 
##  Shapiro-Wilk normality test
## 
## data:  placebo
## W = 0.90872, p-value = 0.4599
var.test(droga, placebo)
## 
##  F test to compare two variances
## 
## data:  droga and placebo
## F = 0.58658, num df = 4, denom df = 4, p-value = 0.618
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.06107284 5.63378553
## sample estimates:
## ratio of variances 
##          0.5865759
t.test(droga, placebo, alternative = "less", var.equal = TRUE, paired = FALSE)
## 
##  Two Sample t-test
## 
## data:  droga and placebo
## t = -1.5758, df = 8, p-value = 0.07686
## alternative hypothesis: true difference in means is less than 0
## 95 percent confidence interval:
##       -Inf 0.1620608
## sample estimates:
## mean of x mean of y 
##      2.06      2.96

Exercício 5

Considere o conjunto de dados a seguir:

X.i: -8.7, -8.3, -8.2, -7.2, -6.1, -6.0, -4.1, -2.0, -1.9, -1.6, -1.3, -0.2, 0.7, 1.3, 1.6, 2.1, 2.2, 4.0, 5.6, 5.9, 6.2, 6.6, 6.7, 8.1

Y.i: -0.6, -0.8, -1.3, -1.9, -2.0, -2.1, -4.0, -4.6, -4.7, -5.5, -5.6, -6.0, 4.6, 4.4, 4.2, 3.9, 3.8, 3.5, 3.1, 2.6, 2.0, 1.2, 0.6, 0.4

  1. Faça um gráfico de dispersão para checar o comportamento dos dados. Por ele, você diria que existe associação entre X e Y?
X.i = c(-8.7, -8.3, -8.2, -7.2, -6.1, -6.0, -4.1, -2.0, -1.9, -1.6, -1.3, -0.2, 0.7, 1.3, 1.6,          2.1, 2.2, 4.0, 5.6, 5.9, 6.2, 6.6, 6.7, 8.1)

Y.i = c(-0.6, -0.8, -1.3, -1.9, -2.0, -2.1, -4.0, -4.6, -4.7, -5.5, -5.6, -6.0, 4.6, 4.4, 4.2, 3.9, 3.8, 3.5, 3.1, 2.6, 2.0, 1.2, 0.6, 0.4)

# Gráfico de Dispersão
plot(X.i, Y.i)

Pelo gráfico, podemos dizer que não existe associação entre X e Y, pois a nuvem de pontos está uma parte deslocada para o canto inferior do gráfico e a outra pro canto superior, não parecendo ter uma relação linear entre elas.

  1. Construa um intervalo de confiança bootstrap para o coeficiente de correlação de Pearson entre X e Y (assumindo B = 100,500,1000,5000) e compare os resultados obtidos. O coeficiente de correlação de Pearson é adequado para medir o tipo de associação existente entre essas duas variáveis? Justifique.
dados = data.frame(
    X = X.i,
    Y = Y.i
)

IC_normal = data.frame(
    Limite_Inferior = numeric(),
    Limite_Superior = numeric()
)

IC_percentil = data.frame(
    Limite_Inferior = numeric(),
    Limite_Superior = numeric()
)

ro_chapeu = cor(dados[,1], dados[,2])
ro_chapeu
## [1] 0.4991181
tam = c(100, 500, 1000, 5000)

set.seed(100)
for (i in tam){
    B = i
    
    cor_boot = NULL
    
    for (j in 1:B){
        amostra = sample(seq(1,nrow(dados),1),nrow(dados),replace=TRUE)
        cor_boot[j] = cor(dados[amostra,1],dados[amostra,2])
    }
    
    epb = sqrt(var(cor_boot))
    
    LIN = ro_chapeu-1.96*epb
    LSN = ro_chapeu+1.96*epb
    
    nova_linha_normal = data.frame(
        Limite_Inferior = LIN,
        Limite_Superior = LSN
    )
    
    IC_normal = rbind(IC_normal, nova_linha_normal)
    
    LIP = quantile(cor_boot,0.025)
    LSP = quantile(cor_boot,0.975)
    
    nova_linha_percentil = data.frame(
        Limite_Inferior = LIP,
        Limite_Superior = LSP
    )
    
    IC_percentil = rbind(IC_percentil, nova_linha_percentil)
  
}

rownames(IC_normal) = c("B = 100", "B = 500", "B = 1000",
                         "B = 5000")

rownames(IC_percentil) = c("B = 100", "B = 500", "B = 1000",
                         "B = 5000")

kable(IC_normal, 
      digits = 3, 
      caption = "IC bootstrap para coef corr de Pearson supondo normalidade - 95% de confiança"
      )
IC bootstrap para coef corr de Pearson supondo normalidade - 95% de confiança
Limite_Inferior Limite_Superior
B = 100 0.361 0.637
B = 500 0.347 0.652
B = 1000 0.357 0.642
B = 5000 0.358 0.640 
kable(IC_percentil, 
      digits = 3, 
      caption = "IC bootstrap para coef corr de Pearson pelo percentil - 95% de confiança"
      )
IC bootstrap para coef corr de Pearson pelo percentil - 95% de confiança
Limite_Inferior Limite_Superior
B = 100 0.393 0.643
B = 500 0.350 0.655
B = 1000 0.359 0.654
B = 5000 0.363 0.646

Exercício 6

Simule 100 amostras \(X_1,...,X_{20}\) de uma população \(N(\theta, 1)\) com \(\theta=1\).

  1. Para cada amostra, calcule a estimativa da variância e do viés para \(\hat{\theta} = \bar{X}\) através do bootstrap não paramétrico, paramétrico e Jackknife. Compare os resultados.
estimativas = data.frame(
    Amostra = numeric(),
    Var_Boot_NPar = numeric(),
    Var_Boot_Par = numeric(),
    Var_Jack = numeric(),
    Vies_Boot_NPar = numeric(),
    Vies_Boot_Par = numeric(),
    Vies_Jack = numeric()
)

intervalos = data.frame(
    Amostra = numeric(),
    LI_N = numeric(),
    LS_N = numeric(),
    LI_boot_npar = numeric(),
    LS_boot_npar = numeric(),
    LI_boot_par = numeric(),
    LS_boot_par = numeric(),
    LI_jack = numeric(),
    LS_jack = numeric()
)

quant_ic_normal = 0
quant_ic_boot_npar = 0
quant_ic_boot_par = 0
quant_ic_jack = 0

B = 3000
set.seed(400)

for (i in 1:100) {
    amostra = rnorm(20, 1, 1)
    theta_hat = mean(amostra)
    
    # IC paramétrico
    li_n = theta_hat - 1.96*(1/sqrt(20))
    ls_n = theta_hat + 1.96*(1/sqrt(20))
    
    # O verdadeiro valor de theta é 1
    if (li_n < 1 && ls_n > 1)  {
        quant_ic_normal = quant_ic_normal + 1
    }
    
    theta_boot = NULL
    
    # Bootstrap não paramétrico
    for (j in 1:B){
        amostra_boot = sample(seq(1,length(amostra),1),length(amostra),replace=TRUE)
        theta_boot[j] = mean(amostra[amostra_boot])
    }
    
    var_boot_npar = var(theta_boot)
    vies_boot_npar = mean(theta_boot) - theta_hat
    
    # IC Bootstrap não paramétrico
    li_npar = quantile(theta_boot, 0.025)
    ls_npar = quantile(theta_boot, 0.975)
    
    # O verdadeiro valor de theta é 1
    if (li_npar < 1 && ls_npar > 1)  {
        quant_ic_boot_npar = quant_ic_boot_npar + 1
    }
    
    amostra_theta_hat = rnorm(20, theta_hat, 1)
    
    theta_hat_boot = NULL
    
    # Bootstrap paramétrico
    for (k in 1:B){
        amostra_boot = sample(seq(1,length(amostra_theta_hat),1),length(amostra_theta_hat),replace=TRUE)
        theta_hat_boot[k] = mean(amostra_theta_hat[amostra_boot])
    }
    
    var_boot_par = var(theta_hat_boot)
    vies_boot_par = mean(theta_hat_boot) - theta_hat
    
    # IC Bootstrap paramétrico
    li_par = quantile(theta_hat_boot, 0.025)
    ls_par = quantile(theta_hat_boot, 0.975)
    
    # O verdadeiro valor de theta é 1
    if (li_par < 1 && ls_par > 1)  {
        quant_ic_boot_par = quant_ic_boot_par + 1
    }
    
    theta_jack = NULL
    n = length(amostra)
    
    # Jacknife
    for (l in 1:n) {
        amostra_jack = amostra[-l]  # remove a i-ésima observação
        theta_jack[l] = mean(amostra_jack)
    }
    
    var_jack = ((n-1)/n)*sum((theta_jack - mean(theta_jack))^2)
    vies_jack = (n-1)*(mean(theta_jack) - theta_hat)
    
    # IC Jacknife
    li_jack = theta_hat - 1.96*sqrt(var_jack)
    ls_jack = theta_hat + 1.96*sqrt(var_jack)
    
    # O verdadeiro valor de theta é 1
    if (li_jack < 1 && ls_jack > 1)  {
        quant_ic_jack = quant_ic_jack + 1
    }
    
    nova_linha = data.frame(
        Amostra = i,
        Var_Boot_NPar = var_boot_npar,
        Var_Boot_Par = var_boot_par,
        Var_Jack = var_jack,
        Vies_Boot_NPar = vies_boot_npar,
        Vies_Boot_Par = vies_boot_par,
        Vies_Jack = vies_jack
    )
    
    estimativas = rbind(estimativas, nova_linha)
    
    nova_linha_ic = data.frame(
        Amostra = i,
        LI_N = li_n,
        LS_N = ls_n,
        LI_boot_npar = li_npar,
        LS_boot_npar = ls_npar,
        LI_boot_par = li_par,
        LS_boot_par = ls_par,
        LI_jack = li_jack,
        LS_jack = ls_jack
    )
    
    intervalos = rbind(intervalos, nova_linha_ic)
    
}

kable(estimativas, 
      digits = 4, 
      caption = "Estimativas para theta = X_barra"
      )
Estimativas para theta = X_barra
Amostra Var_Boot_NPar Var_Boot_Par Var_Jack Vies_Boot_NPar Vies_Boot_Par Vies_Jack
1 0.0339 0.0370 0.0357 -0.0006 0.1968 0
2 0.0819 0.0486 0.0885 0.0082 0.2971 0
3 0.0386 0.0637 0.0417 0.0029 -0.2136 0
4 0.0481 0.0309 0.0511 -0.0003 -0.3346 0
5 0.0493 0.0562 0.0514 -0.0009 0.0857 0
6 0.0661 0.0335 0.0693 0.0034 0.1368 0
7 0.0639 0.0312 0.0660 0.0015 0.3094 0
8 0.0513 0.0335 0.0537 -0.0010 -0.2429 0
9 0.0248 0.0192 0.0277 0.0014 -0.0868 0
10 0.0389 0.0536 0.0410 0.0027 0.1670 0
11 0.0433 0.0290 0.0468 0.0004 -0.2056 0
12 0.0439 0.0251 0.0477 -0.0016 0.1150 0
13 0.0420 0.0397 0.0452 0.0026 -0.2764 0
14 0.0604 0.0645 0.0637 0.0037 -0.0417 0
15 0.0602 0.0706 0.0642 0.0031 -0.4207 0
16 0.0781 0.0726 0.0847 -0.0088 0.3525 0
17 0.0449 0.0482 0.0466 -0.0026 -0.3507 0
18 0.0529 0.0628 0.0579 0.0027 0.2132 0
19 0.0341 0.0577 0.0361 -0.0001 -0.1195 0
20 0.0634 0.0348 0.0670 0.0026 0.1561 0
21 0.0619 0.0590 0.0662 0.0041 -0.1902 0
22 0.0443 0.0379 0.0461 -0.0083 -0.0112 0
23 0.0275 0.0508 0.0303 -0.0019 -0.2128 0
24 0.0423 0.0554 0.0425 -0.0011 -0.2197 0
25 0.0271 0.0254 0.0291 0.0036 -0.1211 0
26 0.0306 0.0581 0.0332 0.0013 0.1386 0
27 0.0685 0.0646 0.0719 0.0045 0.1004 0
28 0.0613 0.0510 0.0635 -0.0079 0.3447 0
29 0.0335 0.0548 0.0374 0.0006 -0.0152 0
30 0.0408 0.0551 0.0425 0.0050 0.1764 0
31 0.0411 0.0369 0.0416 -0.0005 0.2514 0
32 0.0351 0.0598 0.0372 0.0004 0.0999 0
33 0.0391 0.0673 0.0410 0.0015 0.1178 0
34 0.0187 0.0329 0.0194 0.0007 -0.1711 0
35 0.0352 0.0490 0.0361 0.0069 0.2056 0
36 0.0523 0.0425 0.0548 -0.0028 -0.0885 0
37 0.0646 0.0290 0.0709 0.0002 0.1208 0
38 0.0242 0.0414 0.0244 -0.0016 0.3616 0
39 0.0664 0.0447 0.0698 -0.0026 -0.1401 0
40 0.0454 0.0412 0.0479 0.0049 0.2568 0
41 0.0399 0.0277 0.0415 0.0003 0.1984 0
42 0.0226 0.0342 0.0240 0.0024 -0.1119 0
43 0.0476 0.0631 0.0501 0.0016 -0.2171 0
44 0.0452 0.0793 0.0465 -0.0008 -0.2618 0
45 0.0496 0.0756 0.0532 -0.0029 -0.0901 0
46 0.0717 0.0436 0.0728 0.0007 0.1009 0
47 0.0388 0.0641 0.0406 0.0004 -0.3183 0
48 0.0456 0.0424 0.0487 -0.0003 -0.2306 0
49 0.0488 0.0460 0.0499 -0.0041 0.0400 0
50 0.0488 0.0372 0.0522 0.0007 -0.1278 0
51 0.0719 0.0286 0.0760 -0.0067 -0.1347 0
52 0.0312 0.0363 0.0332 0.0019 0.0572 0
53 0.0573 0.0463 0.0588 0.0018 -0.0551 0
54 0.0579 0.0717 0.0604 0.0068 0.1117 0
55 0.0704 0.0631 0.0729 -0.0020 -0.0780 0
56 0.0578 0.0646 0.0606 -0.0027 0.0789 0
57 0.0670 0.0623 0.0671 0.0004 -0.2867 0
58 0.0438 0.0253 0.0480 -0.0009 -0.0497 0
59 0.0306 0.0677 0.0328 -0.0011 -0.3178 0
60 0.0695 0.0492 0.0726 -0.0011 0.0717 0
61 0.0551 0.0427 0.0576 -0.0039 -0.0413 0
62 0.0628 0.0620 0.0664 -0.0030 0.6110 0
63 0.0507 0.0361 0.0542 -0.0037 -0.0224 0
64 0.0681 0.0633 0.0701 0.0004 0.3694 0
65 0.0511 0.0291 0.0542 0.0060 0.0265 0
66 0.0625 0.0502 0.0642 -0.0105 -0.1018 0
67 0.0502 0.0754 0.0515 0.0038 -0.0430 0
68 0.0407 0.0229 0.0433 -0.0059 -0.0042 0
69 0.0467 0.0481 0.0479 0.0036 0.4003 0
70 0.0258 0.0503 0.0276 -0.0009 -0.3097 0
71 0.0407 0.0626 0.0436 -0.0029 -0.5826 0
72 0.0517 0.0417 0.0526 0.0010 0.1699 0
73 0.0313 0.0362 0.0333 -0.0032 0.1835 0
74 0.0346 0.0450 0.0358 -0.0046 0.1346 0
75 0.0535 0.0313 0.0558 -0.0089 0.0438 0
76 0.0477 0.0431 0.0502 -0.0069 0.0225 0
77 0.0562 0.0484 0.0580 -0.0040 0.0586 0
78 0.0723 0.0438 0.0757 -0.0026 -0.2065 0
79 0.0448 0.0398 0.0473 0.0043 -0.1114 0
80 0.0445 0.0508 0.0484 -0.0025 0.1095 0
81 0.0613 0.0380 0.0665 0.0042 -0.1611 0
82 0.0472 0.0301 0.0492 0.0025 0.0473 0
83 0.1135 0.0283 0.1174 -0.0012 0.1488 0
84 0.0489 0.0325 0.0497 -0.0005 -0.1592 0
85 0.0273 0.0298 0.0278 0.0000 -0.0396 0
86 0.0325 0.0221 0.0356 0.0002 -0.2755 0
87 0.0166 0.1128 0.0183 -0.0023 -0.2543 0
88 0.0507 0.0748 0.0544 0.0034 0.1150 0
89 0.0790 0.0628 0.0854 0.0102 -0.1253 0
90 0.0311 0.0387 0.0315 0.0059 0.5409 0
91 0.0311 0.0719 0.0323 0.0014 -0.0398 0
92 0.0524 0.0239 0.0574 -0.0037 -0.2262 0
93 0.0474 0.0445 0.0492 0.0087 -0.0095 0
94 0.0398 0.0282 0.0442 0.0006 0.6449 0
95 0.0538 0.0785 0.0568 -0.0013 -0.0903 0
96 0.0583 0.0553 0.0630 -0.0029 0.1178 0
97 0.0496 0.0810 0.0514 -0.0001 0.2673 0
98 0.0368 0.0386 0.0390 -0.0038 0.0570 0
99 0.0216 0.0232 0.0223 0.0035 -0.4213 0
100 0.0234 0.0456 0.0245 0.0005 0.0934 0
  1. Para cada amostra, calcule o intervalo de confiança de 95% para \(\theta\), assumindo normalidade dos dados (usando IC paramétrico), bootstrap não paramétrico, paramétrico e Jackknife. Qual a porcentagem de intervalos (para cada método de estimação) que contém o valor verdadeiro do parâmetro? Compare os resultados e apresente justificativas para os resultados obtidos.
kable(intervalos, 
      digits = 4, 
      caption = "Intervalos de Confiança para theta (média populacional)"
      )
Intervalos de Confiança para theta (média populacional)
Amostra LI_N LS_N LI_boot_npar LS_boot_npar LI_boot_par LS_boot_par LI_jack LS_jack
2.5% 1 0.3059 1.1825 0.3758 1.1143 0.5462 1.2953 0.3736 1.1147
2.5%1 2 0.7479 1.6244 0.6174 1.7348 1.0415 1.8827 0.6031 1.7691
2.5%2 3 0.2807 1.1572 0.3555 1.1219 0.0001 0.9908 0.3185 1.1193
2.5%3 4 0.6800 1.5565 0.6853 1.5442 0.4462 1.1459 0.6753 1.5612
2.5%4 5 0.7135 1.5900 0.6993 1.5692 0.7728 1.6956 0.7075 1.5960
2.5%5 6 0.3624 1.2389 0.2920 1.3004 0.5909 1.2914 0.2848 1.3166
2.5%6 7 0.4838 1.3604 0.4187 1.4194 0.9057 1.5924 0.4185 1.4257
2.5%7 8 0.7011 1.5777 0.7038 1.5780 0.5450 1.2490 0.6851 1.5937
2.5%8 9 0.7428 1.6194 0.8837 1.4984 0.8289 1.3693 0.8547 1.5075
2.5%9 10 0.6057 1.4822 0.6824 1.4375 0.7531 1.6627 0.6472 1.4407
2.5%10 11 0.0254 0.9019 0.0469 0.8595 -0.0749 0.5915 0.0397 0.8876
2.5%11 12 0.5976 1.4741 0.6481 1.4461 0.8335 1.4461 0.6078 1.4639
2.5%12 13 0.6339 1.5104 0.6610 1.4619 0.4125 1.1905 0.6553 1.4891
2.5%13 14 0.4343 1.3109 0.4037 1.3667 0.3331 1.3360 0.3777 1.3675
2.5%14 15 0.1622 1.0387 0.1288 1.0949 -0.3343 0.7080 0.1038 1.0970
2.5%15 16 0.5866 1.4632 0.4967 1.5863 0.8715 1.9053 0.4546 1.5952
2.5%16 17 0.5101 1.3867 0.5504 1.3668 0.1576 1.0105 0.5253 1.3715
2.5%17 18 0.4155 1.2920 0.4083 1.3108 0.5854 1.5853 0.3823 1.3252
2.5%18 19 0.5328 1.4094 0.6117 1.3367 0.4003 1.3345 0.5985 1.3436
2.5%19 20 0.6757 1.5523 0.6329 1.6151 0.9209 1.6398 0.6067 1.6213
2.5%20 21 0.3051 1.1816 0.2975 1.2568 0.0758 1.0285 0.2392 1.2475
2.5%21 22 0.6485 1.5250 0.6748 1.4807 0.6799 1.4307 0.6659 1.5076
2.5%22 23 0.4727 1.3492 0.5893 1.2356 0.2840 1.1744 0.5700 1.2520
2.5%23 24 0.0982 0.9747 0.1013 0.9317 -0.1714 0.7793 0.1326 0.9403
2.5%24 25 0.9267 1.8033 1.0401 1.6778 0.9285 1.5642 1.0308 1.6992
2.5%25 26 0.1186 0.9952 0.2141 0.9072 0.2485 1.1934 0.1999 0.9138
2.5%26 27 0.8530 1.7296 0.7572 1.7914 0.8839 1.8693 0.7656 1.8170
2.5%27 28 0.4418 1.3183 0.3812 1.3396 0.7941 1.6750 0.3862 1.3739
2.5%28 29 0.8667 1.7432 0.9525 1.6650 0.8262 1.7444 0.9258 1.6841
2.5%29 30 0.4040 1.2805 0.4500 1.2275 0.5458 1.4647 0.4382 1.2462
2.5%30 31 0.3731 1.2497 0.4199 1.2026 0.6947 1.4322 0.4115 1.2113
2.5%31 32 0.6064 1.4829 0.6936 1.4125 0.6661 1.6158 0.6664 1.4229
2.5%32 33 0.7256 1.6022 0.7702 1.5679 0.7706 1.8068 0.7671 1.5607
2.5%33 34 0.4814 1.3579 0.6418 1.1926 0.4026 1.1205 0.6468 1.1925
2.5%34 35 0.3131 1.1896 0.3831 1.1280 0.5226 1.3853 0.3788 1.1239
2.5%35 36 0.9961 1.8727 1.0005 1.8983 0.9239 1.7297 0.9757 1.8931
2.5%36 37 0.7160 1.5926 0.6513 1.6468 0.9483 1.6205 0.6325 1.6760
2.5%37 38 0.6024 1.4789 0.7436 1.3527 1.0112 1.7891 0.7345 1.3468
2.5%38 39 0.2227 1.0993 0.1401 1.1525 0.1151 0.9379 0.1432 1.1787
2.5%39 40 0.4673 1.3438 0.4763 1.3215 0.7796 1.5737 0.4765 1.3345
2.5%40 41 1.0548 1.9313 1.0840 1.8636 1.3517 2.0192 1.0940 1.8921
2.5%41 42 0.6236 1.5002 0.7834 1.3642 0.5945 1.3231 0.7581 1.3657
2.5%42 43 0.7289 1.6054 0.7265 1.5820 0.4491 1.4119 0.7285 1.6059
2.5%43 44 0.5667 1.4432 0.5800 1.4148 0.2250 1.3268 0.5823 1.4275
2.5%44 45 0.4338 1.3103 0.4388 1.2983 0.2413 1.2968 0.4201 1.3240
2.5%45 46 0.7455 1.6221 0.6717 1.7033 0.8630 1.6927 0.6549 1.7127
2.5%46 47 0.4739 1.3505 0.5398 1.2983 0.0871 1.0693 0.5172 1.3072
2.5%47 48 0.4910 1.3676 0.5041 1.3321 0.2902 1.1077 0.4968 1.3618
2.5%48 49 0.5440 1.4205 0.5217 1.3810 0.6169 1.4413 0.5443 1.4201
2.5%49 50 0.2484 1.1250 0.2401 1.1078 0.1969 0.9496 0.2389 1.1345
2.5%50 51 0.7622 1.6387 0.6788 1.7199 0.7261 1.3855 0.6600 1.7409
2.5%51 52 0.7247 1.6012 0.7941 1.4794 0.8283 1.5768 0.8056 1.5202
2.5%52 53 0.2218 1.0984 0.1820 1.1236 0.1780 1.0088 0.1847 1.1355
2.5%53 54 0.3973 1.2739 0.3553 1.3030 0.4106 1.4533 0.3538 1.3173
2.5%54 55 0.7609 1.6374 0.6904 1.7267 0.6078 1.6145 0.6700 1.7282
2.5%55 56 0.6608 1.5374 0.6128 1.5710 0.6806 1.6607 0.6165 1.5817
2.5%56 57 0.9413 1.8179 0.8762 1.8940 0.6207 1.6086 0.8720 1.8872
2.5%57 58 0.6427 1.5193 0.6468 1.4639 0.7280 1.3455 0.6515 1.5105
2.5%58 59 0.5281 1.4047 0.6132 1.2857 0.1734 1.1512 0.6113 1.3215
2.5%59 60 0.4611 1.3376 0.3644 1.4102 0.5253 1.3901 0.3712 1.4275
2.5%60 61 0.9854 1.8620 0.9677 1.8812 0.9525 1.7556 0.9532 1.8942
2.5%61 62 0.6714 1.5479 0.6056 1.5838 1.2334 2.1936 0.6046 1.6147
2.5%62 63 0.9035 1.7801 0.9072 1.7808 0.9401 1.6822 0.8854 1.7982
2.5%63 64 0.5462 1.4228 0.4611 1.4726 0.8730 1.8462 0.4655 1.5036
2.5%64 65 0.6683 1.5448 0.6529 1.5529 0.7919 1.4580 0.6502 1.5629
2.5%65 66 0.5656 1.4421 0.4945 1.4846 0.4556 1.3547 0.5073 1.5005
2.5%66 67 0.4822 1.3588 0.4699 1.3527 0.3295 1.3908 0.4756 1.3654
2.5%67 68 0.4774 1.3539 0.5439 1.3241 0.6259 1.2132 0.5079 1.3234
2.5%68 69 0.6961 1.5726 0.7103 1.5566 1.1307 1.9969 0.7052 1.5635
2.5%69 70 0.9381 1.8146 1.0550 1.6908 0.6369 1.5109 1.0508 1.7019
2.5%70 71 0.1734 1.0499 0.1856 0.9688 -0.4659 0.5033 0.2025 1.0208
2.5%71 72 0.7333 1.6098 0.7283 1.6032 0.9361 1.7443 0.7221 1.6211
2.5%72 73 0.4748 1.3513 0.5628 1.2607 0.7194 1.4645 0.5552 1.2709
2.5%73 74 0.8068 1.6833 0.8846 1.6132 0.9355 1.7639 0.8744 1.6157
2.5%74 75 0.4799 1.3565 0.4316 1.3288 0.6126 1.3010 0.4554 1.3810
2.5%75 76 0.7289 1.6054 0.7223 1.5876 0.7654 1.5635 0.7279 1.6064
2.5%76 77 0.8110 1.6876 0.7797 1.6831 0.8710 1.7305 0.7772 1.7214
2.5%77 78 0.4349 1.3115 0.3579 1.4140 0.2616 1.0884 0.3339 1.4125
2.5%78 79 0.5366 1.4131 0.5776 1.4069 0.4671 1.2390 0.5485 1.4013
2.5%79 80 0.2985 1.1751 0.3321 1.1705 0.3736 1.2641 0.3054 1.1681
2.5%80 81 0.6641 1.5406 0.6141 1.5726 0.5768 1.3402 0.5971 1.6076
2.5%81 82 0.4904 1.3670 0.4967 1.3613 0.6405 1.3176 0.4937 1.3637
2.5%82 83 0.5196 1.3962 0.2967 1.6227 0.7738 1.4291 0.2864 1.6294
2.5%83 84 0.5843 1.4609 0.5646 1.4418 0.5209 1.2235 0.5857 1.4595
2.5%84 85 0.6624 1.5389 0.7926 1.4292 0.7279 1.3987 0.7739 1.4274
2.5%85 86 0.3698 1.2463 0.4373 1.1592 0.2534 0.8345 0.4383 1.1779
2.5%86 87 0.6828 1.5594 0.8805 1.3768 0.1975 1.5054 0.8557 1.3865
2.5%87 88 0.2517 1.1283 0.2633 1.1499 0.2273 1.3146 0.2327 1.1473
2.5%88 89 0.2602 1.1367 0.1825 1.2706 0.0711 1.0483 0.1255 1.2714
2.5%89 90 0.5012 1.3777 0.6062 1.2890 1.0885 1.8460 0.5915 1.2874
2.5%90 91 0.3590 1.2355 0.4499 1.1417 0.2223 1.2492 0.4452 1.1493
2.5%91 92 0.7648 1.6413 0.7379 1.6386 0.6658 1.2684 0.7334 1.6727
2.5%92 93 0.1365 1.0130 0.1624 1.0111 0.1689 0.9920 0.1399 1.0096
2.5%93 94 0.6533 1.5299 0.6853 1.4628 1.4213 2.0759 0.6797 1.5035
2.5%94 95 0.4724 1.3489 0.4282 1.3348 0.2656 1.3467 0.4437 1.3776
2.5%95 96 0.1591 1.0356 0.1370 1.0791 0.2459 1.1637 0.1053 1.0894
2.5%96 97 0.5163 1.3928 0.5304 1.4129 0.6337 1.7486 0.5103 1.3987
2.5%97 98 0.3522 1.2287 0.4297 1.1735 0.4716 1.2453 0.4034 1.1775
2.5%98 99 -0.0485 0.8281 0.0963 0.6747 -0.3315 0.2615 0.0973 0.6823
2.5%99 100 0.5595 1.4361 0.6965 1.2943 0.6637 1.5068 0.6911 1.3045 
cat("A porcentagem de IC's paramétricos (normal) é de ", quant_ic_normal, "%.\n")
## A porcentagem de IC's paramétricos (normal) é de  95 %.
cat("A porcentagem de IC's Bootstrap não paramétrico é de ", quant_ic_boot_npar, "%.\n")
## A porcentagem de IC's Bootstrap não paramétrico é de  91 %.
cat("A porcentagem de IC's Bootstrap paramétrico é de ", quant_ic_boot_par, "%.\n")
## A porcentagem de IC's Bootstrap paramétrico é de  83 %.
cat("A porcentagem de IC's Jacknife é de ", quant_ic_jack, "%.\n")
## A porcentagem de IC's Jacknife é de  93 %.
  1. Repita o item 1 para \(\hat{\theta} = \bar{X^2}\) e compare os resultados.
estimativas_2 = data.frame(
    Amostra = numeric(),
    Var_Boot_NPar = numeric(),
    Var_Boot_Par = numeric(),
    Var_Jack = numeric(),
    Vies_Boot_NPar = numeric(),
    Vies_Boot_Par = numeric(),
    Vies_Jack = numeric()
)

B = 3000
set.seed(400)

for (i in 1:100) {
    amostra = rnorm(20, 1, 1)
    theta_hat = mean(amostra)^2
    
    theta_boot = NULL
    
    # Bootstrap não paramétrico
    for (j in 1:B){
        amostra_boot = sample(seq(1,length(amostra),1),length(amostra),replace=TRUE)
        theta_boot[j] = mean(amostra[amostra_boot])^2
    }
    
    var_boot_npar = var(theta_boot)
    vies_boot_npar = mean(theta_boot) - theta_hat
    
    amostra_theta_hat = rnorm(20, theta_hat, 1)
    
    theta_hat_boot = NULL
    
    # Bootstrap paramétrico
    for (k in 1:B){
        amostra_boot = sample(seq(1,length(amostra_theta_hat),1),length(amostra_theta_hat),replace=TRUE)
        theta_hat_boot[k] = mean(amostra_theta_hat[amostra_boot])^2
    }
    
    var_boot_par = var(theta_hat_boot)
    vies_boot_par = mean(theta_hat_boot) - theta_hat
    
    theta_jack = NULL
    n = length(amostra)
    
    # Jacknife
    for (l in 1:n) {
        amostra_jack = amostra[-l]  # remove a i-ésima observação
        theta_jack[l] = mean(amostra_jack)^2
    }
    
    var_jack = ((n-1)/n)*sum((theta_jack - mean(theta_jack))^2)
    vies_jack = (n-1)*(mean(theta_jack) - theta_hat)
    
    nova_linha = data.frame(
        Amostra = i,
        Var_Boot_NPar = var_boot_npar,
        Var_Boot_Par = var_boot_par,
        Var_Jack = var_jack,
        Vies_Boot_NPar = vies_boot_npar,
        Vies_Boot_Par = vies_boot_par,
        Vies_Jack = vies_jack
    )
    
    estimativas_2 = rbind(estimativas_2, nova_linha)
    
}

kable(estimativas_2, 
      digits = 4, 
      caption = "Estimativas para theta = X_barra^2"
      )
Estimativas para theta = X_barra^2
Amostra Var_Boot_NPar Var_Boot_Par Var_Jack Vies_Boot_NPar Vies_Boot_Par Vies_Jack
1 0.0780 0.0816 0.0778 0.0329 0.0465 0.0357
2 0.4674 0.5568 0.5116 0.1014 1.5454 0.0885
3 0.0857 0.0295 0.0833 0.0427 -0.3612 0.0417
4 0.2433 0.1082 0.2572 0.0474 -0.3808 0.0511
5 0.2597 0.4563 0.2805 0.0471 0.7239 0.0514
6 0.1777 0.0843 0.1849 0.0715 -0.0025 0.0693
7 0.2261 0.1754 0.2234 0.0667 0.5256 0.0660
8 0.2713 0.1528 0.2820 0.0490 -0.1511 0.0537
9 0.1425 0.1336 0.1534 0.0281 0.3355 0.0277
10 0.1784 0.3405 0.1706 0.0445 0.5433 0.0410
11 0.0395 0.0016 0.0415 0.0436 -0.1859 0.0468
12 0.1927 0.1404 0.2059 0.0407 0.3635 0.0477
13 0.1930 0.1261 0.2145 0.0475 -0.3475 0.0452
14 0.1960 0.1392 0.1895 0.0668 -0.1789 0.0637
15 0.0975 0.0101 0.0903 0.0639 -0.2863 0.0642
16 0.3498 0.5903 0.3473 0.0601 0.9905 0.0847
17 0.1675 0.0599 0.1665 0.0400 -0.5501 0.0466
18 0.1589 0.2435 0.1697 0.0576 0.2214 0.0579
19 0.1310 0.1676 0.1355 0.0339 -0.2071 0.0361
20 0.3275 0.2773 0.3250 0.0691 0.7456 0.0670
21 0.1595 0.0385 0.1345 0.0680 -0.3622 0.0662
22 0.2072 0.2039 0.2220 0.0264 0.2252 0.0461
23 0.0931 0.0889 0.1002 0.0240 -0.3983 0.0303
24 0.0493 0.0069 0.0515 0.0411 -0.2278 0.0425
25 0.2031 0.3094 0.2184 0.0370 1.1969 0.0291
26 0.0403 0.0575 0.0407 0.0320 -0.0507 0.0332
27 0.4591 0.7962 0.4829 0.0801 1.5223 0.0719
28 0.1891 0.2644 0.2046 0.0475 0.5290 0.0635
29 0.2303 0.6261 0.2543 0.0352 1.2002 0.0374
30 0.1174 0.1745 0.1237 0.0492 0.1303 0.0425
31 0.1097 0.1274 0.1117 0.0404 0.2062 0.0416
32 0.1582 0.3437 0.1606 0.0360 0.3873 0.0372
33 0.2178 0.5911 0.2204 0.0427 0.8808 0.0410
34 0.0636 0.0633 0.0658 0.0200 -0.3577 0.0194
35 0.0837 0.1178 0.0829 0.0456 0.0776 0.0361
36 0.4427 0.6519 0.4427 0.0442 1.8619 0.0548
37 0.3466 0.2482 0.3906 0.0650 0.8084 0.0709
38 0.1065 0.3495 0.1064 0.0209 1.0450 0.0244
39 0.1204 0.0199 0.1241 0.0629 -0.3041 0.0698
40 0.1518 0.1988 0.1606 0.0544 0.3806 0.0479
41 0.3476 0.6504 0.3833 0.0409 3.6915 0.0415
42 0.1039 0.1449 0.1075 0.0276 -0.0618 0.0240
43 0.2596 0.3318 0.2789 0.0513 0.0123 0.0501
44 0.1842 0.2028 0.1883 0.0435 -0.3709 0.0465
45 0.1555 0.1375 0.1620 0.0445 -0.2355 0.0532
46 0.4103 0.3931 0.4104 0.0733 0.8991 0.0728
47 0.1335 0.0709 0.1341 0.0394 -0.5041 0.0406
48 0.1563 0.0711 0.1718 0.0450 -0.4205 0.0487
49 0.1816 0.1915 0.2025 0.0406 0.0909 0.0499
50 0.0951 0.0215 0.1021 0.0498 -0.3161 0.0522
51 0.4279 0.1943 0.4231 0.0559 0.2940 0.0760
52 0.1651 0.2875 0.1878 0.0356 0.6710 0.0332
53 0.1039 0.0300 0.1076 0.0596 -0.2446 0.0588
54 0.1685 0.1934 0.1702 0.0693 0.0294 0.0604
55 0.4196 0.4649 0.4126 0.0657 0.4746 0.0729
56 0.2839 0.4262 0.2950 0.0518 0.5127 0.0606
57 0.5285 0.6700 0.5038 0.0680 0.7723 0.0671
58 0.1971 0.1292 0.2391 0.0418 0.1085 0.0480
59 0.1115 0.1153 0.1289 0.0285 -0.4867 0.0328
60 0.2292 0.1542 0.2346 0.0675 0.0158 0.0726
61 0.4544 0.6561 0.4671 0.0442 1.9583 0.0576
62 0.3077 0.8438 0.3370 0.0562 2.2246 0.0664
63 0.3697 0.4552 0.3861 0.0407 1.3969 0.0542
64 0.2642 0.4642 0.2804 0.0689 0.8859 0.0701
65 0.2514 0.1804 0.2702 0.0644 0.3696 0.0542
66 0.2552 0.1702 0.2541 0.0416 -0.1368 0.0642
67 0.1723 0.1993 0.1764 0.0572 -0.1251 0.0515
68 0.1444 0.0659 0.1395 0.0300 -0.1197 0.0433
69 0.2420 0.5683 0.2504 0.0549 1.6076 0.0479
70 0.1948 0.5166 0.2097 0.0234 0.6671 0.0276
71 0.0591 0.0208 0.0709 0.0372 -0.2681 0.0436
72 0.2841 0.4001 0.2924 0.0541 1.0482 0.0526
73 0.1061 0.1509 0.1119 0.0254 0.2372 0.0333
74 0.2170 0.4945 0.2213 0.0230 1.3335 0.0358
75 0.1708 0.0985 0.2006 0.0372 -0.0252 0.0558
76 0.2628 0.3268 0.2760 0.0316 0.5984 0.0502
77 0.3449 0.5123 0.3680 0.0463 1.1098 0.0580
78 0.2336 0.0600 0.2255 0.0678 -0.4096 0.0757
79 0.1799 0.1125 0.1735 0.0532 -0.2067 0.0473
80 0.1051 0.0849 0.0985 0.0408 -0.0665 0.0484
81 0.3019 0.1740 0.3263 0.0704 -0.0661 0.0665
82 0.1686 0.1014 0.1694 0.0519 -0.0047 0.0492
83 0.4468 0.1302 0.4289 0.1112 0.2480 0.1174
84 0.2025 0.1074 0.2152 0.0479 -0.2273 0.0497
85 0.1364 0.1669 0.1316 0.0272 0.1916 0.0278
86 0.0857 0.0143 0.0946 0.0328 -0.4884 0.0356
87 0.0847 0.4698 0.0907 0.0114 -0.1390 0.0183
88 0.1059 0.1064 0.1012 0.0554 -0.0520 0.0544
89 0.1770 0.0390 0.1625 0.0934 -0.2937 0.0854
90 0.1132 0.3139 0.1117 0.0421 1.1824 0.0315
91 0.0805 0.1053 0.0829 0.0332 -0.2088 0.0323
92 0.3038 0.1422 0.3392 0.0434 0.0678 0.0574
93 0.0714 0.0242 0.0630 0.0574 -0.1829 0.0492
94 0.1868 0.3854 0.2191 0.0411 2.2093 0.0442
95 0.1757 0.1741 0.1956 0.0514 -0.2047 0.0568
96 0.0915 0.0557 0.0858 0.0549 -0.0763 0.0630
97 0.1878 0.4489 0.1829 0.0494 0.5585 0.0514
98 0.0954 0.0753 0.0961 0.0309 -0.1215 0.0390
99 0.0137 0.0079 0.0141 0.0243 -0.0562 0.0223
100 0.0937 0.2183 0.0970 0.0244 0.2359 0.0245

Exercício 7

Curvas de dose-resposta, como a apresentada abaixo, são muito utilizadas em estudos biológicos e na indústria farmacêutica. Suponha que uma certa droga (medida em mililitros) é administrada em porcos para verificar se uma reação específica ocorre (câncer, diabetes etc.). Cinco porcos recebem cada um de alguns níveis de concentração da droga (X) e a porcentagem (Y) de animais que apresentam a reação (em cada dose) é anotada.

Dose (X): 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5 % de resposta (Y): 0, 0, 20, 0, 40, 60, 40, 80, 100, 100

  1. Faça um gráfico de dispersão. O valor esperado da variável resposta parece ser uma função linear da dosagem da droga?
x = c(0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5)
y = c(0, 0, 20, 0, 40, 60, 40, 80, 100, 100)

plot(x, y,
     main = "Curva Dose-Resposta (Gráfico de Dispersão)",
     xlab = "Dose (ml)",
     ylab = "% de resposta",
     pch = 19,
     col = "blue",
     ylim = c(0, 110))

  1. Estime a reta de regressão de Y em função de X pelo método dos mínimos quadrados. Plote a curva no mesmo gráfico do item 1 e verifique se a estimação está boa.
# Ajuste do modelo linear
modelo = lm(y ~ x)

# Resumo do modelo
summary(modelo)
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -25.4545  -0.4545   3.6364   9.0909  12.7273 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -24.000      9.748  -2.462   0.0392 *  
## x             24.727      3.142   7.869 4.91e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14.27 on 8 degrees of freedom
## Multiple R-squared:  0.8856, Adjusted R-squared:  0.8713 
## F-statistic: 61.93 on 1 and 8 DF,  p-value: 4.915e-05
plot(x, y,
     main = "Curva Dose-Resposta (Gráfico de Dispersão)",
     xlab = "Dose (ml)",
     ylab = "% de resposta",
     pch = 19,
     col = "blue",
     ylim = c(0, 110))
abline(modelo, col = "red", lwd = 2)

  1. Use a reta de regressão dos mínimos quadrados para estimar o valor esperado de Y quando X = 0.5 mililitros. Essa estimativa parece razoável para você?
# Estimativa pontual para X = 0.5
dose = data.frame(x = 0.5)
predict(modelo, dose)
##         1 
## -11.63636
  1. Verifique, usando intervalos de confiança do bootstrap com reamostragem de resíduos (assumindo o estimador de mínimos quadrados de uma função linear) e bootstrap com reamostragem dos dados (assumindo estimadores de mínimos quadrados de uma função linear), se um aumento de 1 mililitro da dose é acompanhado por um aumento médio de 20 pontos percentuais na proporção de porcos que apresentam a reação.
set.seed(123)
B <- 2000
beta1_resid <- numeric(B)

# Valores ajustados e resíduos
yhat <- fitted(modelo)
resid <- residuals(modelo)

for (b in 1:B) {
  e_star <- sample(resid, size = length(resid), replace = TRUE)
  y_star <- yhat + e_star
  modelo_star <- lm(y_star ~ x)
  beta1_resid[b] <- coef(modelo_star)[2]
}

# IC bootstrap percentil
quantile(beta1_resid, c(0.025, 0.975))
##     2.5%    97.5% 
## 19.19967 30.01278
# -------------------------------------------------------------

set.seed(123)
beta1_pairs <- numeric(B)

for (b in 1:B) {
  idx <- sample(1:length(x), size = length(x), replace = TRUE)
  x_star <- x[idx]
  y_star <- y[idx]
  modelo_star <- lm(y_star ~ x_star)
  beta1_pairs[b] <- coef(modelo_star)[2]
}

# IC bootstrap percentil
quantile(beta1_pairs, c(0.025, 0.975))
##     2.5%    97.5% 
## 19.53435 31.62834