如何创建一个循环来对数据集进行子采样、运行特定方程并为您提供每个子样本的 p 值列表？

Question

之前有一个问题提出了类似的问题（有没有办法创建一个循环，我提供一个函数和数据帧并对它进行子采样，并使用子样本重复该函数 N 次？），但是出现了解决方案提供如何获取数据集两部分之间的简单相关性的 p 值的答案，而不是如我上面所述的如何对特定函数进行二次采样并获取 p 值？

具体来说，我有一个如下所示的函数：

AvianGEE <- compar.gee(Dichromatism ~Temp*Precip, data = BirdData, phy = TrimAvianTree)

其输出提供标题：

Estimate         S.E.       t      Pr(T > |t|)

我希望能够对随机 100 种鸟类进行 999 次子采样，并获得一个列表 广义估计方程每次运行的 Pr(T > |t|) 值。

compar.gee 是 ape 包的一个函数，并生成“compar.gee”类的对象

Answer 1

涉及几个步骤：

1. 随机对数据进行子集化
2. 运行

   compar.gee

   并从结果中提取所需的信息
3. 重复上面

1 有多种方法可以随机对数据帧进行子集化，例如。 G。

sample_n

来自 {dplyr}。让我们用基本 R 来推出一个穷人的版本，如下所示：

sample_N <- \(d, N) d[sample(1:nrow(d), N),]

2 虽然

coef

是提取模型系数的选择函数，但

compar.gee

打印一些额外信息作为副作用，但不会将其作为结果的一部分返回。通过检查

ape:::print.compar.gee

，我们可以将相应的计算拉入自定义函数中，该函数返回而不是打印额外的信息：

coef_summary <- function (x){ ## x being your `compar.gee` output
    nas <- is.na(x$coef)
    coef <- x$coef[!nas]
    cnames <- names(coef)
    coef <- matrix(rep(coef, 4), ncol = 4)
    dimnames(coef) <- list(cnames, c("Estimate", "S.E.", "t", "Pr(T > |t|)"))
    df <- x$dfP - dim(coef)[1]
    coef[, 2] <- sqrt(diag(x$W))
    coef[, 3] <- coef[, 1]/coef[, 2]
    if (df < 0) {
        warning("not enough degrees of freedom to compute P-values.")
        coef[, 4] <- NA
    }
    else coef[, 4] <- 2 * (1 - pt(abs(coef[, 3]), df))
    coef
}

...现在您可以提取所需的系数描述符，例如，e。克：

coef_summary(your_compare.gee_here)[, "Pr(T > |t|)"]

3 我们现在可以使用，e。 G。

replicate

获得引导估计：

replicate(n = 999, {res <- compar.gee(Dichromatism ~Temp*Precip,
                                    data = sample_N(BirdData, 100), 
                                    phy = TrimAvianTree
                                    )
                  coef_summary(res)[, "Pr(T > |t|)"]                 
                  }
)

请注意，由于缺乏示例数据，以上是一次演练。可能需要重新调整

replicate

d 代码的输出。