r - 引导分层/多级数据(重采样集群)

Question

我正在生成一个脚本，用于从 cats 创建 bootstrap 示例数据集(来自 -MASS- 包)。

按照戴维森和欣克利教科书 [1]，我运行了一个简单的线性回归，并采用了一个基本的非参数程序从 iid 观察中进行引导，即 对重采样 .

原始样本的形式为:

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) 提出的策略 (p. 100)。
假设我们画了 B = 100 sample 。它们中的每一个都应该由 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)。 Bootstrap 方法及其应用。剑桥大学出版社。

Answer 1

如果我理解正确的话，这就是你试图用 c.data 做的事情。作为输入:

用替换重新采样集群

保持随机样本中的每个簇与其原始数据集(即 c.data)中的点之间的关联

使用采样集群创建 bootstrap

这是一个实现此目的的脚本，您可以将其包装成一个函数以重复 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