引导分层/多级数据（重采样集群）

Question

我正在生成一个脚本，用于从 cats 数据集（来自 -MASS- 包）创建引导样本。

按照 Davidson 和 Hinkley 教科书 [1]，我运行了一个简单的线性回归，并采用了一种基本的非参数程序来从 iid 观察中引导，即pairs resampling。

原始样本格式为：

Bwt   Hwt

2.0   7.0
2.1   7.2

...

1.9    6.8

通过一个单变量线性模型，我们想通过它们的大脑重量来解释猫的壁炉重量。

代码是：

library(MASS)
library(boot)


##################
#   CATS MODEL   #
##################

cats.lm <- glm(Hwt ~ Bwt, data=cats)
cats.diag <- glm.diag.plots(cats.lm, ret=T)


#######################
#   CASE resampling   #
#######################

cats.fit <- function(data) coef(glm(data$Hwt ~ data$Bwt)) 
statistic.coef <- function(data, i) cats.fit(data[i,]) 

bootl <- boot(data=cats, statistic=statistic.coef, R=999)

现在假设存在一个聚类变量cluster = 1, 2,..., 24（例如，每只猫都属于给定的垃圾）。为简单起见，假设数据是平衡的：每个集群有 6 个观察值。因此，24 窝中的每一窝都由 6 只猫组成（即n_cluster = 6 和n = 144）。

可以通过以下方式创建一个假的cluster变量：

q <- rep(1:24, times=6)
cluster <- sample(q)
c.data <- cbind(cats, cluster)

我有两个相关的问题：

如何根据（集群）数据集结构模拟样本？也就是说，如何在集群级别重新采样？我想对集群进行替换采样，并将每个选定集群内的观察值设置为原始数据集中（即用替换集群进行采样而没有替换每个集群内的观察值）。

这是 Davidson 提出的策略（第 100 页）。 假设我们抽取B = 100 样本。它们中的每一个都应该由 24 个可能经常出现的集群（例如cluster = 3, 3, 1, 4, 12, 11, 12, 5, 6, 8, 17, 19, 10, 9, 7, 7, 16, 18, 24, 23, 11, 15, 20, 1）组成，并且每个集群应该包含与原始数据集相同的 6 个观察值。如何在R 中做到这一点？ （有或没有-boot- 包。）您有其他建议吗？

第二个问题涉及初始回归模型。假设我采用固定效应模型，具有集群级截距。 是否改变了采用的重采样程序？

[1] Davidson, A. C., Hinkley, D. V. (1997)。 引导方法及其应用。剑桥大学出版社。

Answer 1

如果我理解正确，这就是您尝试使用 c.data 作为输入的操作：

• 使用替换重新采样集群
• 维护随机样本中每个聚类与其原始数据集（即 c.data）中的点之间的关联
• 使用采样集群创建引导程序

这是一个实现此目的的脚本，您可以将其包装到一个函数中以重复它 R 次，其中 R 是引导程序复制的数量

q <- rep(1:24, times=6)
cluster <- sample(q)
c.data <- cbind(cats, cluster)

# get a vector with all clusters
c <- sort(unique(c.data$cluster))

# group the data points per cluster
clust.group <- function(c) {
    c.data[c.data$cluster==c,]
}

clust.list <- lapply(c,clust.group)

# resample clusters with replacement
c.sample <- sample(c, replace=T)

clust.sample <- clust.list[c.sample]

clust.size <- 6

# combine the cluster list back to a single data matrix
clust.bind <- function(c) {
    matrix(unlist(c),nrow=clust.size)
}

c.boot <- do.call(rbind,lapply(clust.sample,clust.bind))

# Just to maintain columns name
colnames(c.boot) <- names(c.data)

# the new data set (single bootstrap replicate)
c.boot

Answer 2

我尝试通过以下方式解决问题。 虽然它有效，但它可能在速度和“优雅”方面有所改进。另外，如果可能的话，我宁愿找到一种使用-boot- 包的方法，因为它允许通过boot.ci 自动计算一些自举置信区间...

为简单起见，起始数据集包含嵌套在 6 个实验室（聚类变量）中的 18 只猫（“较低级别”观察结果）。数据集是平衡的（每个集群为n_cluster = 3）。我们有一个回归器x，用于解释y。

假数据集和存储结果的矩阵是：

  # fake sample 
  dat <- expand.grid(cat=factor(1:3), lab=factor(1:6))
  dat <- cbind(dat, x=runif(18), y=runif(18, 2, 5))

  # empty matrix for storing coefficients estimates and standard errors of x
  B <- 50 # number of bootstrap samples
  b.sample <- matrix(nrow=B, ncol=3, dimnames=list(c(), c("sim", "b_x", "se_x")))
  b.sample[,1] <- rep(1:B)

在每次B 迭代中，以下循环对 6 个带替换的集群进行采样，每个集群由 3 只未替换采样的猫组成（即集群的内部组成保持不变）。回归系数及其标准误差的估计值存储在先前创建的矩阵中：

  ####################################
  #   loop through "b.sample" rows   #
  ####################################

  for (i in seq(1:B)) {

  ###   sampling with replacement from the clustering variable   

    # sampling with replacement from "cluster" 
    cls <- sample(unique(dat$lab), replace=TRUE)
    cls.col <- data.frame(lab=cls)

    # reconstructing the overall simulated sample
    cls.resample <- merge(cls.col, dat, by="lab")


  ###   fitting linear model to simulated data    

    # model fit
    mod.fit <- function(data) glm(data$y ~ data$x)

    # estimated coefficients and standard errors
    b_x <- summary(mod.fit(data=cls.resample))$coefficients[2,1]
    se_x <- summary(mod.fit(data=cls.resample))$coefficients[2,2]

    b.sample[i,2] <- b_x
    b.sample[i,3] <- se_x

  }

最终的自举标准错误是：

 boot_se_x <- sum(b.sample[,3])/(B-1) 
 boot_se_x

有什么方法可以采用-boot- 包来做到这一点？

另外，关于在集群级别采用固定效应而不是简单的线性回归，我的主要疑问是，在一些模拟样本中，我们可能没有选择一些初始集群，因此无法识别相关的特定于集群的截距。如果您查看我发布的代码，从“机械”的角度来看，这应该不是问题（在每次迭代中，我们可以仅使用采样集群的截距拟合不同的 FE 模型）。

我想知道这一切是否存在统计问题