14

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

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

原始样本的形式为：

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)

我有两个相关的问题：

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

这是戴维森（第 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] 戴维森，AC，欣克利，DV（1997 年）。引导方法及其应用。剑桥大学出版社。

4

2 回答 2

7

如果我对您的理解正确，那么您正在尝试将其c.data作为输入：

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

这是一个实现此目的的脚本，您可以将其包装到一个函数中以重复它 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
于 2013-01-03T01:26:03.597 回答
0

我尝试通过以下方式解决问题。虽然它有效，但它可能在速度和“优雅”方面有所改进。此外，如果可能的话，我宁愿找到一种使用该-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 模型）。

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

于 2013-01-03T01:34:09.223 回答