如何在 R 中进行参数引导？

Question

我正在使用 R 中 faraway 包中的 orings 数据集。我编写了以下分组二项式模型：

orings_model <- glm(cbind(damage, 6-damage) ~ temp, family = binomial, data = orings)
summary(orings_model)

然后我构建了卡方检验统计量并计算了 p 值：

pchisq(orings_model$null.deviance, orings_model$df.null,lower=FALSE)

首先，我想使用 rbinom 和损坏 O 形圈的平均比例（即变量“损坏”）生成此测试统计的空分布下的数据。其次，我想用这个新数据重新计算上述测试统计数据。我不知道该怎么做。

其次，我想处理1000次以上，保存测试统计每一次。我也不知道该怎么做。我倾向于使用 for 循环，但我不确定如何设置它。任何帮助将不胜感激！

Answer 1

目前还不完全清楚您要在这里做什么，但我们至少可以展示一些关于如何实现这一目标的快速原则，然后希望您能实现您的目标。

1) 模拟空模型

不完全清楚您想在此处模拟空模型。看起来您更像是对模拟实际模型拟合感兴趣。请注意，零模型是具有 cbind(damage, 6-damage) ~ 1 形式的模型，零偏差和 df 来自此模型。无论哪种方式，我们都可以使用基础 R 中的 simulate 函数模拟来自模型的数据。

sims <- simulate(orings_model, 1000)

如果您想采用手动方式估计模型的平均向量并将其用于调用 rbinom

的概率

nsim <- 1000 * nrow(orings)
probs <- predict(orings_model, type = 'response')
sims_man <- matrix(rbinom(nsim, 6, probs), 
                   ncol = 1000)
# Check they are equal:
# rowMeans(sims_man) - probs

在第一个版本中，我们得到一个 data.frame，其中有 1000 个 columns，每个都有一个 n 乘以 2 的矩阵（伤害与非伤害）。在后者中，我们只是召唤 damage 结果。

2) 执行 bootstrapping

您可以使用上面的数据手动执行此操作。

# Data from simulate
statfun <- function(x){
  data <- orings_model$data
  data$damage <- if(length(dim(x)) > 1) 
    x[, 1] 
  else 
    x
  newmod <- update(orings_model, data = data)
  pchisq(newmod$null.deviance, newmod$df.null, lower=FALSE)
}
sapply(sims, statfun)

# data from manual method
apply(sims_man, 2, statfun)

或者可以花一些时间使用 boot 函数，允许以标准化的方式执行 bootstrap:

library(boot)
# See help("boot")
ran_gen <- function(data, mle){
  data$damage <- simulate(orings_model)[[1]][,1]
  data
}
boot_metric <- function(data, w){
  model <- glm(cbind(damage = damage, not_damage = 6 - damage) ~ temp, 
               family = binomial, data = data)
  pchisq(model$null.deviance, 
         model$df.null,
         lower=FALSE)
}
boots <- boot(orings, boot_metric, 
     R = 1000, 
     sim = 'parametric', 
     ran.gen = ran_gen, 
     mle = pchisq(orings_model$null.deviance, 
                  orings_model$df.null,
                  lower=FALSE))

此时我们有 boots$t 中的统计数据和 boots$t0 中的空统计数据，因此可以使用 sum(boots$t > boots$t0) / boots$R 估计一个简单的统计数据（R 是复制数） .

如何在 R 中进行参数引导？

How to conduct parametric bootstrapping in R?

for-loop

r

glm

1) 模拟空模型

2) 执行 bootstrapping