Intervalos de confiança pelo método Bootstrap
Neste arquivo, vamos apresentar os métodos utilizados no cálculo dos intervalos de confiança apresentados na tese de Mestrado: ‘Produção normativa, participação e tipos de regras: uma análise das audiências públicas da Comissão de Valores Mobiliários’, de Luiz Cantarelli. Os intervalos foram calculados utilizando a linguagem de programação R.1
Começamos pelo download dos pacotes utilizados. 2 3 4
library(readxl)
library(boot)
library(tidyverse)## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──## ✓ ggplot2 3.3.3 ✓ purrr 0.3.4
## ✓ tibble 3.0.4 ✓ dplyr 1.0.2
## ✓ tidyr 1.1.2 ✓ stringr 1.4.0
## ✓ readr 1.4.0 ✓ forcats 0.5.0## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()Antes de iniciarmos o cálculo dos intervalos de confiança, fixamos uma referência para a randomização dos dados, com set.seed(), de modo que os resultados aqui obtidos possam ser reproduzidos posteriormente.5
set.seed(12345)Em seguida, carregamos os dados da base, que se encontram em uma planilha limpa de Excel, utilizando a função read_xlsx() do pacote readxl, e visualizamos a estrutura da planilha.
base <- read_xlsx('base.xlsx')
glimpse(base)## Rows: 3,313
## Columns: 12
## $ Audiência <chr> "CVM_2018_Audiencia Pública SDM 01/2018_2…
## $ Participante <chr> "Associação Brasileira das Companhias Abe…
## $ `Categoria de Participante` <chr> "Agente Econômico Regulado", "Agente Econ…
## $ Contribuição <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13…
## $ Substância <chr> "Sim", "Sim", "Sim", "Sim", "Sim", "Não",…
## $ Caráter <chr> "Não", "Não", "Não", "Não", "Não", "Não",…
## $ Status <chr> "Não", "Não", "Não", "Não", "Não", "Não",…
## $ Estrutura <chr> "Não", "Sim", "Não", "Sim", "Sim", "Sim",…
## $ `Precisão/Vagueza` <dbl> 0, 1, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ `Simplicidade/Complexidade` <dbl> 0, 0, 0, 1, 1, 0, -1, 0, 0, 1, 0, 0, 0, 0…
## $ `Clareza/Opacidade` <dbl> 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1,…
## $ Impacto <chr> "Sim", "Não", "Sim", "Não", "Não", "Não",…Podemos observar que os dados da planilha são heterogêneos. Em algumas colunas, temos apenas duas respostas possíveis, ‘Sim’ ou ‘Não’, que se encontram em formato de texto. Em outras linhas, temos colunas numéricas com 3 respostas possíveis: 1, 0 ou -1.
Assim sendo, vamos recodificar as variáveis de interesse que estão em formato de texto para um formato lógico (‘Sim’ -> TRUE, ‘Não’ -> FALSE). De forma semelhante, as variáveis classificadas em 1, 0 ou -1 serão desagregadas em 3 colunas distintas, com valores lógicos para cada possibilidade de classificação.
Obtemos uma tabela com os nomes das audiências, e todas as demais variáveis em formato booleano. Essa padronização dos dados facilitará a aplicação do modelo.
base_recode <- tibble(audiencia = base$Audiência,
substância_sim = ifelse(base$Substância == "Sim", TRUE, FALSE),
caráter_sim = ifelse(base$Caráter == "Sim", TRUE, FALSE),
status_sim = ifelse(base$Status == "Sim", TRUE, FALSE),
estrutura_sim = ifelse(base$Estrutura == "Sim", TRUE, FALSE),
impacto_sim = ifelse(base$Impacto == "Sim", TRUE, FALSE),
mais_precisão = ifelse(base$`Precisão/Vagueza` == 1, TRUE, FALSE),
neutro_precisão = ifelse(base$`Precisão/Vagueza` == 0, TRUE, FALSE),
menos_precisão = ifelse(base$`Precisão/Vagueza` == -1, TRUE, FALSE),
mais_complexidade = ifelse(base$`Simplicidade/Complexidade` == 1, TRUE, FALSE),
neutro_complexidade = ifelse(base$`Simplicidade/Complexidade` == 0, TRUE, FALSE),
menos_complexidade = ifelse(base$`Simplicidade/Complexidade` == -1, TRUE, FALSE),
mais_clareza = ifelse(base$`Clareza/Opacidade` == 1, TRUE, FALSE),
neutro_clareza = ifelse(base$`Clareza/Opacidade` == 0, TRUE, FALSE),
menos_clareza = ifelse(base$`Clareza/Opacidade` == -1, TRUE, FALSE))
glimpse(base_recode)## Rows: 3,313
## Columns: 15
## $ audiencia <chr> "CVM_2018_Audiencia Pública SDM 01/2018_258", "CV…
## $ substância_sim <lgl> TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, TRUE, TRUE, …
## $ caráter_sim <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, …
## $ status_sim <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, …
## $ estrutura_sim <lgl> FALSE, TRUE, FALSE, TRUE, TRUE, TRUE, TRUE, FALSE…
## $ impacto_sim <lgl> TRUE, FALSE, TRUE, FALSE, FALSE, FALSE, FALSE, FA…
## $ mais_precisão <lgl> FALSE, TRUE, FALSE, TRUE, FALSE, FALSE, FALSE, FA…
## $ neutro_precisão <lgl> TRUE, FALSE, TRUE, FALSE, FALSE, TRUE, TRUE, TRUE…
## $ menos_precisão <lgl> FALSE, FALSE, FALSE, FALSE, TRUE, FALSE, FALSE, F…
## $ mais_complexidade <lgl> FALSE, FALSE, FALSE, TRUE, TRUE, FALSE, FALSE, FA…
## $ neutro_complexidade <lgl> TRUE, TRUE, TRUE, FALSE, FALSE, TRUE, FALSE, TRUE…
## $ menos_complexidade <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, TRUE, F…
## $ mais_clareza <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, TRUE, FALSE, F…
## $ neutro_clareza <lgl> TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, TRUE, TRUE, …
## $ menos_clareza <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, …Como dito, nosso objetivo é calcular, pelo método Bootstrap,6 os intervalos de confiança para a proporção de diferentes classificações encontradas na amostra. Como nossa amostragem inicial foi construída de forma agregada (clustering) no nível das audiências, a randomização do modelo Bootstrap deve ser feita também no nível das audiências.
Com as variáveis da tabela padronizadas, torna-se possível criar uma única função básica para o cálculo de todas essas proporções (calc_prop()). Nossa função receberá um vetor com as 37 audiências da amostra (audiencias), os indices gerados randomicamente pelo modelo Bootstrap (idx), a tabela recodificada (dados) e o índice da coluna em que se encontra, na tabela recodificada, a variável que queremos analisar (variavel). Com esses dados, o modelo deve calcular a proporção de casos em que, para as audiências selecionadas aleatoriamente, o valor da variável de interesse é TRUE.
calc_prop <- function (audiencias, idx, dados, variavel) {
# Filtra apenas as linhas referentes às audiências randomicamente selecionadas pelo modelo
selecao <- filter(dados, audiencia %in% audiencias[idx])
# Calcula a proporção de valores TRUE na coluna desejada
sum(selecao[[variavel]]) / length(selecao[[variavel]])
}Agora podemos passar nossa função e as colunas desejadas para a o modelo (boot()), indicando quantas amostras aleatórias serão criadas para testar desvios da proporção original. Vamos agregar os resultados em uma lista, aplicando a função boot() a cada variável de interesse. Em todos os casos, utilizaremos 1000 repetições (R = 1000).
apply_boot <- function (input, tabela = base_recode, STAT = calc_prop, BOOT = boot) {
resultado <- vector('list', (length(tabela) - 1))
for (i in seq_along(resultado)) { # <- Aplicando o modelo Bootstrap a cada coluna
resultado[[i]] <- BOOT(data = input,
statistic = STAT,
dados = tabela,
variavel = (i + 1),
R = 1000)
}
names(resultado) <- names(tabela)[2:length(tabela)] # <- Nomeando a lista com as colunas
resultado
}
boot_list <- apply_boot(input = unique(base$Audiência))Obtemos assim uma lista com os resultados do modelo. Podemos acessar dados gerais sobre cada resultado, ou a própria distribuição de repetições calculadas para cada variável e representá-las graficamente. Os exemplos a seguir mostram os resultados do modelo para a coluna substância_sim.
boot_list$substância_sim##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## BOOT(data = input, statistic = STAT, R = 1000, dados = tabela,
## variavel = (i + 1))
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 0.7802596 -0.0004646734 0.01068734plot(boot_list$substância_sim)
Com esses dados, podemos enfim calcular intervalos de confiança por meio do método Bootstrap, utilizando a função boot.ci(). O método será o da aproximação à distribuição normal (intervalo de confiança de Wald), indicado na chamada da função por type = 'norm'. O nível de confiança, indicado pelo parâmetro conf, será de 95%.
O código abaixo salva os intervalos de confiança para todos as variáveis em uma nova lista (ci_list).
ci_list <- boot_list %>% map(~boot.ci(boot.out = ., conf = 0.95, type = 'norm'))
names(ci_list) <- names(boot_list)Por último, vamos calcular os intervalos de confiança para todas as variáveis e apresentar os resultados obtidos em uma tabela final, com as seguintes colunas:
- variável = nome da variável,
- proporção = indica qual proporção foi calculada (…sim = “Sim”, mais… = 1, neutro… = 0, menos… = -1),
- conf = nível de confiança adotado (95% em todos os casos),
- intervalo = intervalo para o nível de confiança adotado (inferior - superior),
- erro_padrão = o erro padrão das repetições geradas pelo modelo,
- viés = o viés das repetições geradas pelo modelo.
res <- unlist(lapply(boot_list, `[[`, 1))
ci_1 <- unlist(lapply(lapply(ci_list, `[[`, 4), `[[`, 2))
ci_2 <- unlist(lapply(lapply(ci_list, `[[`, 4), `[[`, 3))
ep <- unlist(lapply(lapply(boot_list, `[[`, 2), sd))
vs <- unlist(lapply(lapply(boot_list, `[[`, 2), mean)) - res
resultados <- tibble(
variável = names(base)[c(5, 6, 7, 8, 12, 9, 9, 9, 10, 10, 10, 11, 11, 11)],
proporção = names(boot_list),
resultado = round(res, digits = 4),
conf = 0.95,
intervalo = paste(round(ci_1, digits = 4), "-", round(ci_2, digits = 4)),
erro_padrão = round(ep, digits = 4),
viés = round(vs, digits = 4)
)
resultadosApêndice: Resultados completos para todas as variáveis
A seguir, apresentamos resultados mais completos, incluindo o erro, o desvio e a representação gráfica da distribuição aleatória gerada pelo modelo, para todas as variáveis.
for (i in 1:14) { #
writeLines(paste("****************************************\n",
str_to_upper(names(boot_list)[i]),
"\n****************************************"))
print(boot_list[[i]])
writeLines("\n")
print(ci_list[[i]])
if (i < 14) {
plot(boot_list[[i]], title = "names(boot_list)[i]")
}
}## ****************************************
## SUBSTÂNCIA_SIM
## ****************************************
##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## BOOT(data = input, statistic = STAT, R = 1000, dados = tabela,
## variavel = (i + 1))
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 0.7802596 -0.0004646734 0.01068734
##
##
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = ., conf = 0.95, type = "norm")
##
## Intervals :
## Level Normal
## 95% ( 0.7598, 0.8017 )
## Calculations and Intervals on Original Scale
## ****************************************
## CARÁTER_SIM
## ****************************************
##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## BOOT(data = input, statistic = STAT, R = 1000, dados = tabela,
## variavel = (i + 1))
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 0.01448838 0.0001191986 0.002397516
##
##
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = ., conf = 0.95, type = "norm")
##
## Intervals :
## Level Normal
## 95% ( 0.0097, 0.0191 )
## Calculations and Intervals on Original Scale
## ****************************************
## STATUS_SIM
## ****************************************
##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## BOOT(data = input, statistic = STAT, R = 1000, dados = tabela,
## variavel = (i + 1))
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 0.01448838 -9.194127e-05 0.005703981
##
##
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = ., conf = 0.95, type = "norm")
##
## Intervals :
## Level Normal
## 95% ( 0.0034, 0.0258 )
## Calculations and Intervals on Original Scale
## ****************************************
## ESTRUTURA_SIM
## ****************************************
##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## BOOT(data = input, statistic = STAT, R = 1000, dados = tabela,
## variavel = (i + 1))
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 0.5348627 0.0001364379 0.02180179
##
##
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = ., conf = 0.95, type = "norm")
##
## Intervals :
## Level Normal
## 95% ( 0.4920, 0.5775 )
## Calculations and Intervals on Original Scale
## ****************************************
## IMPACTO_SIM
## ****************************************
##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## BOOT(data = input, statistic = STAT, R = 1000, dados = tabela,
## variavel = (i + 1))
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 0.3673408 0.0006365688 0.02329513
##
##
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = ., conf = 0.95, type = "norm")
##
## Intervals :
## Level Normal
## 95% ( 0.3210, 0.4124 )
## Calculations and Intervals on Original Scale
## ****************************************
## MAIS_PRECISÃO
## ****************************************
##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## BOOT(data = input, statistic = STAT, R = 1000, dados = tabela,
## variavel = (i + 1))
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 0.1762753 -0.0003038381 0.01321758
##
##
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = ., conf = 0.95, type = "norm")
##
## Intervals :
## Level Normal
## 95% ( 0.1507, 0.2025 )
## Calculations and Intervals on Original Scale
## ****************************************
## NEUTRO_PRECISÃO
## ****************************************
##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## BOOT(data = input, statistic = STAT, R = 1000, dados = tabela,
## variavel = (i + 1))
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 0.7872019 0.0002786182 0.01445291
##
##
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = ., conf = 0.95, type = "norm")
##
## Intervals :
## Level Normal
## 95% ( 0.7586, 0.8153 )
## Calculations and Intervals on Original Scale
## ****************************************
## MENOS_PRECISÃO
## ****************************************
##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## BOOT(data = input, statistic = STAT, R = 1000, dados = tabela,
## variavel = (i + 1))
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 0.03652279 1.417438e-05 0.004670633
##
##
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = ., conf = 0.95, type = "norm")
##
## Intervals :
## Level Normal
## 95% ( 0.0274, 0.0457 )
## Calculations and Intervals on Original Scale
## ****************************************
## MAIS_COMPLEXIDADE
## ****************************************
##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## BOOT(data = input, statistic = STAT, R = 1000, dados = tabela,
## variavel = (i + 1))
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 0.2390583 -8.712823e-05 0.01347906
##
##
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = ., conf = 0.95, type = "norm")
##
## Intervals :
## Level Normal
## 95% ( 0.2127, 0.2656 )
## Calculations and Intervals on Original Scale
## ****************************************
## NEUTRO_COMPLEXIDADE
## ****************************************
##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## BOOT(data = input, statistic = STAT, R = 1000, dados = tabela,
## variavel = (i + 1))
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 0.6715967 0.0006185935 0.01578974
##
##
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = ., conf = 0.95, type = "norm")
##
## Intervals :
## Level Normal
## 95% ( 0.6400, 0.7019 )
## Calculations and Intervals on Original Scale
## ****************************************
## MENOS_COMPLEXIDADE
## ****************************************
##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## BOOT(data = input, statistic = STAT, R = 1000, dados = tabela,
## variavel = (i + 1))
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 0.089345 -0.0005457691 0.01021075
##
##
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = ., conf = 0.95, type = "norm")
##
## Intervals :
## Level Normal
## 95% ( 0.0699, 0.1099 )
## Calculations and Intervals on Original Scale
## ****************************************
## MAIS_CLAREZA
## ****************************************
##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## BOOT(data = input, statistic = STAT, R = 1000, dados = tabela,
## variavel = (i + 1))
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 0.127377 0.0004991654 0.008456321
##
##
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = ., conf = 0.95, type = "norm")
##
## Intervals :
## Level Normal
## 95% ( 0.1103, 0.1435 )
## Calculations and Intervals on Original Scale
## ****************************************
## NEUTRO_CLAREZA
## ****************************************
##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## BOOT(data = input, statistic = STAT, R = 1000, dados = tabela,
## variavel = (i + 1))
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 0.872623 0.0004026434 0.008275368
##
##
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = ., conf = 0.95, type = "norm")
##
## Intervals :
## Level Normal
## 95% ( 0.8560, 0.8884 )
## Calculations and Intervals on Original Scale
## ****************************************
## MENOS_CLAREZA
## ****************************************
##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## BOOT(data = input, statistic = STAT, R = 1000, dados = tabela,
## variavel = (i + 1))
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 0 0 0
##
##
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = ., conf = 0.95, type = "norm")
##
## Intervals :
## Level Normal
## 95% ( 0, 0 )
## Calculations and Intervals on Original ScaleR Core Team (2020). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/.↩︎
Hadley Wickham and Jennifer Bryan (2019). readxl: Read Excel Files. R package version 1.3.1. https://CRAN.R-project.org/package=readxl↩︎
Angelo Canty and Brian Ripley (2021). boot: Bootstrap R (S-Plus) Functions. R package version 1.3-27.↩︎
Wickham et al., (2019). Welcome to the tidyverse. Journal of Open Source Software, 4(43), 1686, https://doi.org/10.21105/joss.01686.↩︎
Cf. Stodden, Victoria. 2014. “The Reproducible Research Movement in Statistics”. Statistical Journal of the IAOS 30: 91–93.↩︎
Cf. Efron, B. 1979. “Bootstrap Methods: Another Look at the Jackknife”. The Annals of Statistics 7 (1): 1–26. https://doi.org/10.1214/aos/1176344552.↩︎