Paket Rsample sadrži funkcije koje omogućavaju različite tipove resampling-a (npr. cross-validation, bootstrap, itd). Struktura podataka u kojoj se čuvaju resampling podaci je data frame i veoma je zgodna za dalji rad. Više o paketu Rsample možete pročitati na zvaničnoj stranici paketa: https://github.com/tidymodels/rsample. Paket je deo RStudio “tidymodels” projekta, a na čelu projekta je Max Kuhn.
U prvom postu o paketu Rsample, prikazaćemo primenu paketa za ocenu intrevala poverenja srednje vrednosti. Ovaj primer služi da se prikaže osnovna struktura paketa. Zadatak možemo predstaviti na sledeći način. Želimo da ocenimo 95% interval poverenja srednje vrednosti, a nije nam poznat analitički aparat za ocenu intervala, u ovom slučaju t-test. Ideja je da kreiramo veći broj bootstrap uzoraka (uzorak sa ponavljanjem iste veličine kao originalni uzorak) originalnog skupa i da za svaki izračunamo srednju vrednost. Ocena intervala poverenja bi se zatim izračunala prostim izračunavanjem 0.025 i 0.975 kvanitala iz dobijenih srednjih vrednosti.
Prvo ćemo kreirati uzorak iz normalnog rasporeda sa srednjom vrednošću 5 i standardnom devijacijom 1. Veličina uzorka je 100.
library(rsample) library(ggplot2) library(dplyr) library(purrr) set.seed(54321) data <- data_frame(x = rnorm(100, 5)) data
# A tibble: 100 x 1
x
<dbl>
1 4.82
2 4.07
3 4.22
4 3.35
5 4.59
6 3.90
7 3.31
8 7.52
9 6.40
10 5.18
# ... with 90 more rows
Koristićemo funkciju bootstraps iz Rsample paketa da generišemo 1000 bootstrap uzoraka. Generisani uzorci su smešteni u data frame u koloni splits. Kolona id označava naziv svakog uzorka.
bt_data <- bootstraps(data, times = 1000) bt_data
# Bootstrap sampling # A tibble: 1,000 x 2 splits id <list> <chr> 1 <S3: rsplit> Bootstrap0001 2 <S3: rsplit> Bootstrap0002 3 <S3: rsplit> Bootstrap0003 4 <S3: rsplit> Bootstrap0004 5 <S3: rsplit> Bootstrap0005 6 <S3: rsplit> Bootstrap0006 7 <S3: rsplit> Bootstrap0007 8 <S3: rsplit> Bootstrap0008 9 <S3: rsplit> Bootstrap0009 10 <S3: rsplit> Bootstrap0010 # ... with 990 more rows
Strukturu prvog bootstrap uzorka možemo videti preko sledeće komande:
bt_data$splits[[1]]
<100/39/100>
Vidimo da novi uzorak sadrži 100 opservacija (analysis data u terminologiji Rsample paketa), 39 opservacija je ostalo van skupa (assessment data – često se koristi za procenu performansi modela korišćenog u analizi podataka). Poslednji broj označava veličinu originalnog skupa. Ovi podaci su očekivani s obzirom da se radi o uzorku sa ponavljanjem. Svakom novom uzorku možemo pristupiti pomoću funkcije analysis(sample). Podacima koji su ostali van skupa pristupamo pomoću funkcije assessment(sample).
Nakon što smo formirali bootstrap uzorke, potrebno je izračunati srednje vrednosti za svaki uzorak. Kreiraćemo pomoćnu funkciju get_split_mean za izračunavanje srednje vrednosti uzorka i kasnije je primeniti na svaki uzorak.
get_split_mean <- function(split) {
# access to the sample data
split_data <- analysis(split)
# calculate the sample mean value
split_mean <- mean(split_data$x)
return(split_mean)</code>
}
Pomoću funkcije map_dbl iz purrr paketa, proćićemo kroz sve uzorke i izračunati srednje vrednosti. Novi vektor srednjih vrednosti dodaćemo postojećem bt_data data frame-u. Na ovaj način, dobijamo veoma zgodnu strukturu (tidy data) za dalji rad.
bt_data$bt_means <- map_dbl(bt_data$splits, get_split_mean) bt_data
# Bootstrap sampling # A tibble: 1,000 x 3 splits id bt_means <list> <chr> <dbl> 1 <S3: rsplit> Bootstrap0001 5.02 2 <S3: rsplit> Bootstrap0002 4.93 3 <S3: rsplit> Bootstrap0003 4.99 4 <S3: rsplit> Bootstrap0004 4.86 5 <S3: rsplit> Bootstrap0005 4.90 6 <S3: rsplit> Bootstrap0006 5.03 7 <S3: rsplit> Bootstrap0007 5.10 8 <S3: rsplit> Bootstrap0008 5.00 9 <S3: rsplit> Bootstrap0009 4.78 10 <S3: rsplit> Bootstrap0010 5.07 # ... with 990 more rows
Ocenu 95% intervala poverenja dobijamo pomoću funkcije quantile.
bt_ci <- round(quantile(bt_data$bt_means, c(0.025, 0.975)), 3) bt_ci
2.5% 97.5% 4.704 5.117
Interval poverenja originalnog skupa ćemo izračunati pomoću funkcije t.test.
t_test <- t.test(data$x) tt_ci <- round(t_test$conf.int, 3) tt_ci
[1] 4.697 5.117 attr(,"conf.level") [1] 0.95
Vidimo da su intervali gotovo identični. Na sličan način se mogu dobiti ocene parametara za koje ne postoji analitičko rešenje ili je rešenje suviše kompleksno.
Ocena intervala poverenja pomoću bootstrap uzoraka je prikazana na sledećem grafiku. Kod za dobijeni grafik se nalazi u Dodatku.
Dodatak
# Učitavanje potrebnih R paketa.
library(rsample)
library(ggplot2)
library(dplyr)
library(purrr)
# Generisanje uzorka iz normalnog rasporeda.
set.seed(54321)
data <- data_frame(x = rnorm(100, 5))
data
# Generisanje 1000 bootstrap uzoraka.
bt_data <- bootstraps(data, times = 1000)
bt_data
# Prikazati strukturu prvog bootstrap uzorka.
bt_data$splits[[1]]
# Kreiranje pomoćne funkcije za izračunavanje srednje vrednosti uzorka.
get_split_mean <- function(split) {
# pristupanje podacima uzorka
split_data <- analysis(split)
# izračunavanje srednje vrednosti uzorka
split_mean <- mean(split_data$x)
return(split_mean)
}
# Dodavanje novog vektora postojećoj datoteci.
bt_data$bt_means <- map_dbl(bt_data$splits, get_split_mean)
bt_data
# Izračunavanje 95% interval poverenja.
bt_ci <- round(quantile(bt_data$bt_means, c(0.025, 0.975)), 3)
bt_ci
# Izračunavanje 95% interval poverenja originalnog skupa.
t_test <- t.test(data$x)
tt_ci <- round(t_test$conf.int, 3)
tt_ci
######################
# Generisanje grafika
######################
bt_ci_lower <- quantile(bt_data$bt_means, c(0.025))
bt_ci_upper <- quantile(bt_data$bt_means, c(0.975))
bt_mean <- mean(bt_data$bt_means)
bt_data <- bt_data %>%
dplyr::mutate(Color = ifelse(bt_means < bt_ci_lower | bt_means > bt_ci_upper, "Out of 95% CI", "In 95% CI"))
ggplot(bt_data, aes(x = bt_means, y = 0)) +
geom_jitter(aes(color = Color), alpha = 0.6, size = 3, width = 0) +
geom_vline(xintercept = round(c(bt_ci_lower, bt_mean, bt_ci_upper),2), linetype = c(2,1,2)) +
scale_x_continuous(breaks = round(c(bt_ci_lower, bt_mean, bt_ci_upper),2),
labels = c("Lower CI (95%)", "Mean", "Upper CI (95%)"),
sec.axis = sec_axis(~., breaks = round(c(bt_ci_lower, bt_mean, bt_ci_upper),2),
labels = round(c(bt_ci_lower, bt_mean, bt_ci_upper),2))) +
scale_color_manual(values = c("gray40", "firebrick4")) +
theme_void() +
theme(legend.title = element_blank(),
legend.position = "top",
legend.text = element_text(size = 14),
axis.text.x = element_text(size = 14))
Master Statistike. Razvija R pakete i Shiny aplikacije.