从引导数据的函数生成图

Question

考虑这个数据框：

set.seed(123)
dat1 <- data.frame(Loc = rep(c("a","b","c","d","e","f","g","h"),each = 5),
                   ID = rep(c(1:10), each = 2),
                   var1 = rnorm(200),
                   var2 = rnorm(200),
                   var3 = rnorm(200),
                   var4 = rnorm(200),
                   var5 = rnorm(200),
                   var6 = rnorm(200))
dat1$ID <- factor(dat1$ID)

位置Loc 是每个ID 上测量var1:6 的分组变量。有几对 Locs 彼此非常接近（在地理上），以至于它们可能应该被视为一个单独的组而不是两个独立的组。因此，我编写了一个函数来引导每个变量，以查看这些组是否来自同一分布：

library(tidyverse)
BootT <- function(dat, var, gv1, gv2){
  set.seed(123)
  a<- dplyr::filter(dat, Loc == gv1)
  a2 <- dplyr::select(a, var)
  b <- dplyr::filter(dat, Loc == gv2)
  b2 <- dplyr::select(b, var)
  pooled <- rbind(a2, b2)
  boot.t <- c(1:999)
  for(i in 1:999){
    sample.index <- sample(c(1:length(pooled[,1])), replace = TRUE)
    sample.x <- pooled[sample.index,][1:length(a2[,1])]
    sample.y <- pooled[sample.index,][-c(1:length(b2[,1]))]
    boot.t[i] <- t.test(sample.x, sample.y)$statistic
  }
  p.pooled <-  data.frame(p.pooled = 1 + sum(abs(boot.t) > abs(t.test(a[,var],b[,var])$statistic))) / (999+1) 
 return(p.pooled)
  ids <- data.frame(Group1 = paste0(gv1), Group2 = paste0(gv2), Variable = paste0(var))
  p.pooled <- p.pooled%>%
    dplyr::mutate(Group1 = ids[,1], Group2 = ids[,2], Variable = ids[,3])
  p.pooled <- p.pooled[,c(2,3,4,1)]
 return(p.pooled)
}
#compare 2 locs of interest with a single variable
BootT(dat = dat1, var = "var2", gv1 = "a", gv2 = "g") 
#compare all 6 variables 
vars <- names(dat1[,3:8])
results <- list()
for(i in vars){
  res <- BootT(dat = dat1, var = i, gv1 = "a", gv2 = "b")
  results <- rbind(results, res)
} 

我想修改这个函数，使它输出一个经典的直方图，显示每个变量与观察值的自举分布，并包含图表上的汇总统计信息。我怎样才能修改这个功能来完成这个？ 编辑： 最初，我打算使用引导包来执行此操作，这会更容易，但我不确定我是否理解不同的参数将如何改变采样过程。在两个Locs 具有相等方差的情况下（通过 F 检验评估），我想对合并样本进行抽样，如我在上面演示的那样。但是，当样本是异质的时，我想在创建要比较的合并样本之前减去每个组的平均值（这会强制原假设为真，并且不对同质方差做任何假设）。有关更多信息，请参阅此帖子：https://stats.stackexchange.com/questions/136661/using-bootstrap-under-h0-to-perform-a-test-for-the-difference-of-two-means-repl

我实际上已经做了一个与上面的函数非常相似的函数（另一个非常原始的名称）来处理存在异质方差问题的情况：

BootT2  abs(t.test(a[,var],b[,var])$statistic)) / 999+ 1)-2)
  #p.h0  abs(t.test(a[,var],b[,var])$statistic)) / 999)
  ids %
    变异（Group1 = ids[,1]，Group2 = ids[,2]，变量 = ids[,3]）
  p.h0 
如果有人想解释我如何执行这些程序并使用 boot() 包生成绘图，那就太好了。
    

Answer 1

如果我理解正确，以下将运行数据集 dat1 中变量 var 的 2 个 Loc 的自举 t 检验。它在函数bootTstat 中使用accepted answer 到此CrossValidated post 引导程序，但这是从函数funBoot 调用的。函数funBoot 负责对组gv1 和gv2 行和列var 进行子集化。这样形成的数据集被传递给bootTstat。

bootTstat <- function(x, y, R){
  pool <- c(x, y)
  xt <- x - mean(x) + mean(pool)
  yt <- y - mean(y) + mean(pool)
  boot.t <- numeric(R)
  for (i in seq_len(R)){
    sample.x <- sample(xt, replace = TRUE)
    sample.y <- sample(yt, replace = TRUE)
    boot.t[i] <- t.test(sample.x, sample.y)$statistic
  }
  p.h0 <- (1 + sum(abs(boot.t) > abs(t.test(x, y)$statistic))) / (R + 1)  
  list(
    statistic = boot.t,
    p.value = p.h0
  )
}

funBoot <- function(data, R, var, gv1, gv2){
  i <- data[["Loc"]] == gv1
  j <- data[["Loc"]] == gv2
  x <- data[i, var]
  y <- data[j, var]
  bootTstat(x, y, R)
}

对于"var2" 和组"a" 和"g" 使用整个组数据和R = 1000 测试运行t 检验。

首先是 t 检验。

a <- subset(dat1, Loc == 'a', select = 'var2')
g <- subset(dat1, Loc == 'g', select = 'var2')
t.test(a, g)
#
#        Welch Two Sample t-test
#
#data:  a and g
#t = 1.1002, df = 47, p-value = 0.2769
#alternative hypothesis: true difference in means is not equal to 0
#95 percent confidence interval:
# -0.2585899  0.8828038
#sample estimates:
# mean of x  mean of y 
# 0.1755209 -0.1365860 

还有 bootsrtapped t 检验。 R

b_ag <- funBoot(dat1, R, var = "var2", gv1 = "a", gv2 = "g")
b_ag$p.value
#[1] 0.2737263

这个 p 值类似于之前获得的p.value = 0.2769。
并且可以轻松绘制直方图。

hist(b_ag$statistic, main = "Bootstrapped t-test")

现在对所有变量和组 "a" 和 "b" 运行测试。使用包ggplot2 绘图。

ttest_list <- lapply(names(dat1)[3:8], function(v) {
  b <- funBoot(data = dat1, R = R, var = v, gv1 = "a", gv2 = "b")
  list(
    p.value = b$p.value,
    test = data.frame(var = v, stat = b$statistic)
  )
})

ttest_df <- lapply(ttest_list, '[[', 'test')
ttest_df <- do.call(rbind, ttest_df)

library(ggplot2)

ggplot(ttest_df, aes(stat)) +
  geom_histogram(bins = 25) +
  facet_wrap(~ var)