(Acest articol a fost publicat pentru prima dată pe R pe Jason Bryerși a contribuit cu drag la R-Bloggers). (Puteți raporta problema despre conținutul de pe această pagină aici)
Doriți să vă împărtășiți conținutul pe R-Bloggers? Faceți clic aici dacă aveți un blog sau aici dacă nu.
În timpul unei clase recente, un student a întrebat dacă intervalele de încredere în bootstrap au fost mai robuste decât intervalele de încredere estimate folosind eroarea standard (adică (Se = frac {s} { sqrt {n}} )) Pentru a răspunde la această întrebare am scris o funcție pentru a simula luarea unei grămadă de eșantioane aleatorii dintr -o populație, calculați intervalul de încredere pentru acel eșantion folosind abordarea standard de eroare ( t Distribuția este utilizată în mod implicit, consultați cv parametru. Pentru a utiliza distribuția normală, de exemplu, setată cv = 1.96.), și apoi, de asemenea, calculând un interval de încredere folosind boostrap.
library(dplyr)
library(ggplot2)
#' Simulate random samples to estimate confidence intervals and bootstrap
#' estimates.
#'
#' @param pop a numeric vector representing the population.
#' @param n sample size for each random sample from the population.
#' @param n_samples the number of random samples.
#' @param n_boot number of bootstrap samples to take for each sample.
#' @param seed a seed to use for the random process.
#' @param cv critical value to use for calculating confidence intervals.
#' @return a data.frame with the sample and bootstrap mean and confidence
#' intervals along with a logical variable indicating whether a Type I
#' error would have occurred with that sample.
bootstrap_clt_simulation <- function(
pop,
n = 30,
n_samples = 500,
n_boot = 500,
cv = abs(qt(0.025, df = n - 1)),
seed,
verbose = interactive()
) {
if(missing(seed)) {
seed <- sample(100000)
}
results <- data.frame(
seed = 1:n_samples,
samp_mean = numeric(n_samples),
samp_se = numeric(n_samples),
samp_ci_low = numeric(n_samples),
samp_ci_high = numeric(n_samples),
samp_type1 = logical(n_samples),
boot_mean = numeric(n_samples),
boot_ci_low = numeric(n_samples),
boot_ci_high = numeric(n_samples),
boot_type1 = logical(n_samples)
)
if(verbose) {
pb <- txtProgressBar(min = 0, max = n_samples, style = 3)
}
for(i in 1:n_samples) {
if(verbose) {
setTxtProgressBar(pb, i)
}
set.seed(seed + i)
samp <- sample(pop, size = n)
boot_samp <- numeric(n_boot)
for(j in 1:n_boot) {
boot_samp(j) <- sample(samp, size = length(samp), replace = TRUE) |>
mean()
}
results(i,)$seed <- seed + i
results(i,)$samp_mean <- mean(samp)
results(i,)$samp_se <- sd(samp) / sqrt(length(samp))
results(i,)$samp_ci_low <- mean(samp) - cv * results(i,)$samp_se
results(i,)$samp_ci_high <- mean(samp) + cv * results(i,)$samp_se
results(i,)$samp_type1 <- results(i,)$samp_ci_low > mean(pop) |
mean(pop) > results(i,)$samp_ci_high
results(i,)$boot_mean <- mean(boot_samp)
results(i,)$boot_ci_low <- mean(boot_samp) - cv * sd(boot_samp)
results(i,)$boot_ci_high <- mean(boot_samp) + cv * sd(boot_samp)
results(i,)$boot_type1 <- results(i,)$boot_ci_low > mean(pop) |
mean(pop) > results(i,)$boot_ci_high
}
if(verbose) {
close(pb)
}
return(results)
}
Distribuție uniformă pentru populație
Să începem cu o distribuție uniformă pentru populația noastră.
pop_unif <- runif(1e5, 0, 1) ggplot(data.frame(x = pop_unif), aes(x = x)) + geom_density()
Media populației este 0,4999484. Acum putem simula probe și estimările lor de bootstrap corespunzătoare.
results_unif <- bootstrap_clt_simulation(pop = pop_unif, seed = 42, verbose = FALSE)
4% din eșantioanele noastre nu au conținut media populației în intervalul de încredere (adică rata de eroare de tip I) în comparație cu rMedia (rezultate_unif $ boot_type1) * 100`% din estimările bootstrap. Următorul tabel compară erorile de tip I pentru fiecare eșantion în comparație cu bootstrap -ul estimat din acel eșantion.
tab <- table(results_unif$samp_type1, results_unif$boot_type1, useNA = 'ifany') tab ## ## FALSE TRUE ## FALSE 477 3 ## TRUE 0 20
În general, comiterea unei erori de tip I este aceeași, indiferent de metodă, deși au existat 3 cazuri în care bootstrap -ul ar fi dus la o rată de eroare de tip I în care abordarea standard de eroare nu ar fi.
Următoarele parcele arată relația dintre media estimată (stânga) și lățimea intervalului de condamnare (dreapta) pentru fiecare eșantion și bootstrap -ul corespunzător.
results_unif |>
ggplot(aes(x = samp_mean, y = boot_mean)) +
geom_vline(xintercept = mean(pop_unif), color = 'blue') +
geom_hline(yintercept = mean(pop_unif), color = 'blue') +
geom_abline() +
geom_point() +
ggtitle("Sample mean vs bootstrap mean")
results_unif |>
dplyr::mutate(samp_ci_width = samp_ci_high - samp_ci_low,
boot_ci_width = boot_ci_high - boot_ci_low) |>
ggplot(aes(x = samp_ci_width, y = boot_ci_width)) +
geom_abline() +
geom_point() +
ggtitle('Sample vs boostrap confidence interval width')
Distribuție înclinată pentru populație
Vom repeta aceeași analiză folosind o distribuție pozitivă înclinată.
pop_skewed <- rnbinom(1e5, 3, .5) ggplot(data.frame(x = pop_skewed), aes(x = x)) + geom_density(bw = 0.75)
Media populației pentru această distribuție este 2.99792
results_skewed <- bootstrap_clt_simulation(pop = pop_skewed, seed = 42, verbose = FALSE)
mean(results_skewed$samp_type1) # Percent of samples with Type I error
## (1) 0.05
mean(results_skewed$boot_type1) # Percent of bootstrap estimates with Type I error
## (1) 0.052
# CLT vs Bootstrap Type I error rate
table(results_skewed$samp_type1, results_skewed$boot_type1, useNA = 'ifany')
##
## FALSE TRUE
## FALSE 473 2
## TRUE 1 24
results_skewed |>
ggplot(aes(x = samp_mean, y = boot_mean)) +
geom_vline(xintercept = mean(pop_skewed), color = 'blue') +
geom_hline(yintercept = mean(pop_skewed), color = 'blue') +
geom_abline() +
geom_point() +
ggtitle("Sample mean vs bootstrap mean")
results_skewed |>
dplyr::mutate(samp_ci_width = samp_ci_high - samp_ci_low,
boot_ci_width = boot_ci_high - boot_ci_low) |>
ggplot(aes(x = samp_ci_width, y = boot_ci_width)) +
geom_abline() +
geom_point() +
ggtitle('Sample vs boostrap confidence interval width')
Putem vedea că rezultatele sunt foarte asemănătoare cu cea a distirubului uniform. Explorarea unui caz în care bootstrap -ul ar fi dus la o eroare de tip I în care abordarea standard de eroare nu ar dezvălui că este foarte aproape, diferența este mai mică de 0,1.
results_differ <- results_skewed |>
dplyr::filter(!samp_type1 & boot_type1)
results_differ
## seed samp_mean samp_se samp_ci_low samp_ci_high samp_type1 boot_mean
## 1 443 3.866667 0.4516466 2.942946 4.790388 FALSE 3.924733
## 2 474 3.933333 0.4816956 2.948155 4.918511 FALSE 3.956800
## boot_ci_low boot_ci_high boot_type1
## 1 3.044802 4.804665 TRUE
## 2 3.018549 4.895051 TRUE
set.seed(results_differ(1,)$seed)
samp <- sample(pop_skewed, size = 30)
boot_samp <- numeric(500)
for(j in 1:500) {
boot_samp(j) <- sample(samp, size = length(samp), replace = TRUE) |>
mean()
}
cv = abs(qt(0.025, df = 30 - 1))
mean(pop_skewed)
## (1) 2.99792
ci <- c(mean(samp) - cv * sd(samp) / sqrt(30), mean(samp) + cv * sd(samp) / sqrt(30))
ci
## (1) 2.942946 4.790388
mean(pop_skewed) < ci(1) | mean(pop_skewed) > ci(2)
## (1) FALSE
ci_boot <- c(mean(boot_samp) - cv * sd(boot_samp), mean(boot_samp) + cv * sd(boot_samp))
ci_boot
## (1) 3.044802 4.804665
mean(pop_skewed) < ci_boot(1) | mean(pop_skewed) > ci_boot(2)
## (1) TRUE
Adăugarea unui mai vechi
Să luăm în considerare un eșantion care obligă cea mai mare valoare din populație să fie în eșantion.
set.seed(2112)
samp_outlier <- c(sample(pop_skewed, size = 29), max(pop_skewed))
boot_samp <- numeric(500)
for(j in 1:500) {
boot_samp(j) <- sample(samp, size = length(samp), replace = TRUE) |>
mean()
}
ci <- c(mean(samp_outlier) - cv * sd(samp_outlier) / sqrt(30), mean(samp_outlier) + cv * sd(samp_outlier) / sqrt(30))
ci
## (1) 1.647006 4.952994
mean(pop_skewed) < ci(1) | mean(pop_skewed) > ci(2)
## (1) FALSE
ci_boot <- c(mean(boot_samp) - cv * sd(boot_samp), mean(boot_samp) + cv * sd(boot_samp))
ci_boot
## (1) 2.905153 4.781381
mean(pop_skewed) < ci_boot(1) | mean(pop_skewed) > ci_boot(2)
## (1) FALSE
În acest exemplu, vedem că președintele celor mai vechi are un impact mai mare asupra intervalului de încredere, cu intervalul de încredere în bootstrap fiind mult mai mic.
Mărimea eșantionului și a bootstrap -ului legat de eroarea standard
Să explorăm și relația dintre nnumărul de probe de bootstrap și eroare standard. Reamintiți -vă formula pentru eroarea standard este:
(Se = frac { sigma} { sqrt {n}} )
Figura de mai jos prezintă eroarea standard împotriva erorii standard, presupunând că Sigma (abaterea standard) este una. După cum puteți vedea, pur și simplu creșterea mărimii eșantionului va scădea eroarea standard (și, prin urmare, intervalul de încredere).
se <- function(n, sigma = 1) {
sigma / sqrt(n)
}
ggplot() + stat_function(fun = se) + xlim(c(0, 100)) +
ylab('Standard Error') + xlab('Sample Size (n)')
Având în vedere din nou o populație cu o distribuție uniformă, următorul cod va atrage mostre aleatorii cu n variind de la 30 până la 50 în pași de 15. Pentru fiecare dintre aceste probe aleatorii, vom estima, de asemenea, erori standard BooStrap cu numărul de probe de bootstrap care variază între 50 și 1.000 în creșteri de 50.
n <- seq(30, 500, by = 15)
n_boots <- seq(50, 1000, by = 50)
results <- expand.grid(n, n_boots)
attributes(results) <- NULL
results <- as.data.frame(results)
names(results) <- c('n', 'n_boots')
results$samp_mean <- NA
results$samp_se <- NA
results$boot_mean <- NA
results$boot_se <- NA
for(i in seq_len(nrow(results))) {
samp <- sample(pop_unif, size = results(i,)$n)
results(i,)$samp_mean <- mean(samp)
results(i,)$samp_se <- sd(samp) / sqrt(length(samp))
boot_samp_dist <- numeric(results(i,)$n_boots)
for(j in seq_len(results(i,)$n_boots)) {
boot_samp_dist(j) <- sample(samp, size = length(samp), replace = TRUE) |> mean()
}
results(i,)$boot_mean <- mean(boot_samp_dist)
results(i,)$boot_se <- sd(boot_samp_dist)
}
Figura din stânga taie dimensiunea eșantionului față de eroarea standard care, ca și mai sus, arată că pe măsură ce dimensiunea eșantionului crește eroarea standard scade. În dreapta este o diagramă a numărului de eșantioane de bootstrap față de eroarea standard în care culorile punctului corespund dimensiunii eșantionului. Aici vedem că eroarea standard este constantă. Adică numărul de probe de bootstrap nu este legat de eroarea standard. Variabilitatea erorii standard este contabilizată de dimensiunea eșantionului.
y_limits <- c(0, 0.075)
p_samp_size_se <- ggplot(results, aes(x = n, y = samp_se)) +
geom_point(fill = '#9ecae1', color = 'grey50', shape = 21) +
geom_smooth(color = 'darkgreen', se = FALSE, method = 'loess', formula = y ~ x) +
ylim(y_limits) +
ylab('Standard Error') +
xlab('Sample size (n)') +
ggtitle(latex2exp::TeX("Standard Error (SE = \frac{\sigma}{\sqrt{n}})")) +
scale_fill_gradient(low = '#deebf7', high = '#3182bd') +
theme(legend.position = 'bottom')
p_boot_size_se <-
ggplot(results, aes(x = n_boots, y = boot_se)) +
geom_point(aes(fill = n), color = 'grey50', shape = 21) +
geom_smooth(color = 'darkgreen', se = FALSE, method = 'loess', formula = y ~ x) +
ylim(y_limits) +
ylab('Standard Error') +
xlab('Number of Bootstrap Samples') +
ggtitle('Bootstrap Standard Error',
subtitle = '(i.e. standard deviation of the bootstrap sample)') +
scale_fill_gradient(low = '#deebf7', high = '#3182bd') #+ theme(legend.position = 'none')
cowplot::plot_grid(p_samp_size_se, p_boot_size_se)
În cele din urmă, putem complota relația dintre cele două estimări de eroare standard; a cărei corelație este extrem de mare cu r = 1.
ggplot(results, aes(x = samp_se, y = boot_se)) +
geom_abline() +
geom_point() +
xlab('Sample Standard Error') +
ylab('Boostrap Standard Error') +
ggtitle(paste0('Correlation between standard errors = ', round(cor(results$samp_se, results$boot_se), digits = 2))) +
coord_equal()
