BM_project_exp_2_bootstrap

Libraries

library(tidyverse)
library(ggplot2)
library(plotrix)
library(car)
library(rstatix)
library(SimComp)
library(readxl)
library(dplyr)

dataset of the control data

d_BM_exp_2_3_bm_imprint <- read_excel(path = "/mnt/c/Users/Toshiya Matsushima/OneDrive/R projects/BM_project/data analysis using Rstudio/BM_project_dataset.xlsx", sheet = "exp_2_3_bm_imprint")
d_BM_exp_2_3_bm_imprint_00_ctrl <- d_BM_exp_2_3_bm_imprint %>% filter(drug=="00_ctrl") %>%
  arrange(batch) %>%
  mutate(id_num=seq(1:nrow(.)))
d_BM_exp_2_3_bm_imprint_00_ctrl

mean(d_BM_exp_2_3_bm_imprint_00_ctrl$bm)

[1] 136.3089

std.error(d_BM_exp_2_3_bm_imprint_00_ctrl$bm)

[1] 16.13828

mean(d_BM_exp_2_3_bm_imprint_00_ctrl$imprint)

[1] 419.5203

std.error(d_BM_exp_2_3_bm_imprint_00_ctrl$imprint)

[1] 20.25447

BM score

Bootstrap computations

d_BM_exp_2_3_bm_imprint_00_ctrl_BM_bootstrap <- NULL
for (i in 1:10000){
    set.seed(i)
  
    d_BM_exp_2_3_bm_imprint_00_ctrl_test <- d_BM_exp_2_3_bm_imprint_00_ctrl %>% sample_n(size = 10, replace = F) %>% arrange(id_num)
    bm_mean_test_data <- apply(d_BM_exp_2_3_bm_imprint_00_ctrl_test[,4], 2, mean)
    
    d_BM_exp_2_3_bm_imprint_00_ctrl_train <- d_BM_exp_2_3_bm_imprint_00_ctrl %>% setdiff(d_BM_exp_2_3_bm_imprint_00_ctrl_test)
    bm_mean_train_data <- apply(d_BM_exp_2_3_bm_imprint_00_ctrl_train[,4], 2, mean)
    
    temp_diff <- bm_mean_test_data - bm_mean_train_data
    d_BM_exp_2_3_bm_imprint_00_ctrl_BM_bootstrap_temp <- tibble(n_bootstrap = i, mean_diff = temp_diff)
    d_BM_exp_2_3_bm_imprint_00_ctrl_BM_bootstrap <- d_BM_exp_2_3_bm_imprint_00_ctrl_BM_bootstrap %>% bind_rows(d_BM_exp_2_3_bm_imprint_00_ctrl_BM_bootstrap_temp)
}

d_BM_exp_2_3_bm_imprint_00_ctrl_BM_bootstrap

90%, 95%, and 99% confidence intervals (one-sided testing is assumed)

quantile(d_BM_exp_2_3_bm_imprint_00_ctrl_BM_bootstrap$mean_diff, probs = c(0.001, 0.01, 0.05, 0.10))

      0.1%         1%         5%        10% 
-181.13551 -142.16726  -97.65323  -77.08611 

density distribution curve of BM calculated from the control dataset

fig <- d_BM_exp_2_3_bm_imprint_00_ctrl_BM_bootstrap %>%
    ggplot(mapping = aes(x = mean_diff), size =3)+
      geom_vline(xintercept=0, linetype = "solid", color ="blue", size = 1)+
      geom_vline(xintercept=-97.65323, linetype = "solid", color ="dark red", size = 1)+
      geom_vline(xintercept=-142.1673, linetype = "solid", color ="red", size = 1)+
    geom_hline(yintercept=0, linetype = "solid", color = "grey")+
    geom_density()+
    theme_classic(base_size=15)
fig
ggsave(plot = fig, filename = "exp_2_BM_distribution.png", dpi = 300, height = 10, width = 15, units = "cm")

imprint score

Bootstrap computations

d_IMPRINT_exp_2_3_bm_imprint_00_ctrl_BM_bootstrap <- NULL
for (i in 1:10000){
    set.seed(i)
  
    d_BM_exp_2_3_bm_imprint_00_ctrl_test <- d_BM_exp_2_3_bm_imprint_00_ctrl %>% sample_n(size = 10, replace = F) %>% arrange(id_num)
    imprint_mean_test_data <- apply(d_BM_exp_2_3_bm_imprint_00_ctrl_test [,5], 2, mean)
    
    d_BM_exp_2_3_bm_imprint_00_ctrl_train <- d_BM_exp_2_3_bm_imprint_00_ctrl %>% setdiff(d_BM_exp_2_3_bm_imprint_00_ctrl_test)
    imprint_mean_train_data <- apply(d_BM_exp_2_3_bm_imprint_00_ctrl_train [,5], 2, mean)
    
    temp_diff <- imprint_mean_test_data - imprint_mean_train_data
    d_IMPRINT_exp_2_3_bm_imprint_00_ctrl_BM_bootstrap_temp <- tibble(n_bootstrap = i, mean_diff = temp_diff)
    d_IMPRINT_exp_2_3_bm_imprint_00_ctrl_BM_bootstrap <- d_IMPRINT_exp_2_3_bm_imprint_00_ctrl_BM_bootstrap %>% bind_rows(d_IMPRINT_exp_2_3_bm_imprint_00_ctrl_BM_bootstrap_temp)
}

d_IMPRINT_exp_2_3_bm_imprint_00_ctrl_BM_bootstrap

99.9%, 99%, 95%, and 90% confidence intervals (one-sided testing is assumed)

quantile(d_IMPRINT_exp_2_3_bm_imprint_00_ctrl_BM_bootstrap$mean_diff, probs = c(0.001, 0.01, 0.05, 0.10))

      0.1%         1%         5%        10% 
-289.45321 -201.61150 -135.32757  -99.51062 

density distribution curve of BM calculated from the control dataset

fig <- d_IMPRINT_exp_2_3_bm_imprint_00_ctrl_BM_bootstrap %>%
    ggplot(mapping = aes(x = mean_diff), size =3)+
      geom_vline(xintercept=0, linetype = "solid", color ="blue", size = 1)+
      geom_vline(xintercept=-135.32757, linetype = "solid", color ="dark red", size = 1)+
      geom_vline(xintercept=-201.61150, linetype = "solid", color ="red", size = 1)+
    geom_hline(yintercept=0, linetype = "solid", color = "grey")+
    geom_density()+
    theme_classic(base_size=15)
fig
ggsave(plot = fig, filename = "exp_2_IMPRINTING_distribution.png", dpi = 300, height = 10, width = 15, units = "cm")

SessionInfo

sessionInfo()

R version 3.6.3 (2020-02-29)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.3 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0

locale:
 [1] LC_CTYPE=C.UTF-8       LC_NUMERIC=C           LC_TIME=C.UTF-8        LC_COLLATE=C.UTF-8     LC_MONETARY=C.UTF-8    LC_MESSAGES=C.UTF-8    LC_PAPER=C.UTF-8       LC_NAME=C             
 [9] LC_ADDRESS=C           LC_TELEPHONE=C         LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C   

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] plotly_4.10.0    ggsci_2.9        readxl_1.3.1     SimComp_3.3      rstatix_0.7.0    car_3.0-12       carData_3.0-5    plotrix_3.8-2    ggbeeswarm_0.6.0 forcats_0.5.1    stringr_1.4.0   
[12] dplyr_1.0.7      purrr_0.3.4      readr_2.1.1      tidyr_1.1.4      tibble_3.1.6     ggplot2_3.3.5    tidyverse_1.3.1 

loaded via a namespace (and not attached):
 [1] nlme_3.1-144         fs_1.5.2             lubridate_1.8.0      httr_1.4.2           tools_3.6.3          backports_1.4.1      utf8_1.2.2           R6_2.5.1            
 [9] vipor_0.4.5          lazyeval_0.2.2       DBI_1.1.1            mgcv_1.8-31          colorspace_2.0-2     withr_2.4.3          tidyselect_1.1.1     compiler_3.6.3      
[17] survPresmooth_1.1-11 mratios_1.4.2        cli_3.1.0            rvest_1.0.2          xml2_1.3.3           sandwich_3.0-1       labeling_0.4.2       scales_1.1.1        
[25] mvtnorm_1.1-3        digest_0.6.29        rmarkdown_2.11       pkgconfig_2.0.3      htmltools_0.5.2      dbplyr_2.1.1         fastmap_1.1.0        htmlwidgets_1.5.4   
[33] rlang_0.4.12         rstudioapi_0.13      jquerylib_0.1.4      generics_0.1.1       farver_2.1.0         zoo_1.8-9            jsonlite_1.7.2       gtools_3.9.2        
[41] magrittr_2.0.1       Matrix_1.2-18        Rcpp_1.0.7           munsell_0.5.0        fansi_0.5.0          abind_1.4-5          lifecycle_1.0.1      stringi_1.7.6       
[49] multcomp_1.4-18      yaml_2.2.1           MASS_7.3-51.5        grid_3.6.3           crayon_1.4.2         lattice_0.20-40      haven_2.4.3          splines_3.6.3       
[57] hms_1.1.1            knitr_1.37           pillar_1.6.4         codetools_0.2-16     reprex_2.0.1         glue_1.6.0           drc_3.0-1            evaluate_0.14       
[65] data.table_1.14.2    modelr_0.1.8         vctrs_0.3.8          tzdb_0.2.0           cellranger_1.1.0     gtable_0.3.0         assertthat_0.2.1     xfun_0.29           
[73] broom_0.7.10         viridisLite_0.4.0    survival_3.1-8       beeswarm_0.4.0       TH.data_1.1-0        ellipsis_0.3.2