通过将变量改组 1000 次进行回归 -R 的排列测试

Question

我有一个回归，我想使用置换测试。我想改组我的自变量并重新运行回归 1000 次。在运行回归时，我希望通过将置换系数大于或等于初始回归系数的次数相加来计算回归中每个系数的 ap 值（没有改组变量的回归）。

我为此构建了一些代码，但是 1) 它非常慢（尤其是在添加了更多协变量/交互效果的情况下）以及 2) 需要为每个系数运行 for 循环。

initial<-lm(y~x*x1, data=df)
summary(initial)

N=1000
PermuteFunction<-function(y=df$y, x=df$x, x1=df$x1)){
  model.resample=lm(y~sample(x, replace=F)*x1)
  #permutes predictor, then runs model
  sum=summary(model.resample)$coefficients[2]
  sum1=summary(model.resample)$coefficients[3]
  return(sum)
}

sum=numeric(N)
sum1=numeric(N)
for (i in 1:N){
  sum[i]<-PermuteFunction()
  sum1[i]<-PermuteFunction()}


#calculating p-value
length(sum[sum>=initial$coefficients[2]])/N
length(sum[sum1>=initial$coefficients[2]])/N

有没有更有效的方法来做到这一点？我在这个论坛上看到了与置换测试相关的其他问题，但无法找到一个看起来为每个系数计算 p 值的问题。 lmPerm 似乎也不是用来洗牌的（除非我误解了它的功能）

Answer 1

以下并不慢，它在明显延迟后输出sim（不是sum，基本R函数的名称）。另见dcarlon's comment。
该函数的参数是：

1. data - 数据集；
2. regr - 字符串形式的回归器名称；
3. resp - 字符串形式的响应名称。

使用内置数据集 iris 测试。

PermuteFunction <- function(data, regr, resp){
  i <- sample(nrow(data))
  x <- data[i, regr]
  y <- data[[resp]]
  coef(lm(y ~ x))[2]
}

initial <- coef(lm(Sepal.Length ~ Sepal.Width, iris))[2]

set.seed(2021)
R <- 1e3
sim <- replicate(R, PermuteFunction(iris, "Sepal.Width", "Sepal.Length"))

mean(sim >= initial)
#[1] 0.946

---

编辑。

如果需要所有系数，以下函数将使用rowMeans 将初始模型的系数与回归系数进行比较。

PermuteFunction <- function(data, regr, resp){
  i <- sample(nrow(data))
  x <- data[i, regr]
  y <- data[[resp]]
  d <- cbind.data.frame(y, x)
  names(d) <- c(resp, regr)
  regr <- paste(regr, collapse = "+")
  fmla <- as.formula(paste(resp, regr, sep = "~"))
  coef(lm(fmla, data = d))
}

initial <- coef(lm(Sepal.Length ~ Sepal.Width, iris))

现在测试功能。请注意，第一个测试给出了相同的斜率 0.946。

set.seed(2021)
R <- 1e3
sim <- replicate(R, PermuteFunction(iris, "Sepal.Width", "Sepal.Length"))

rowMeans(sim >= initial)
#(Intercept) Sepal.Width 
#      0.054       0.946 

现在是一个有 2 个回归量的测试。

initial2 <- coef(lm(Sepal.Length ~ Sepal.Width + Petal.Length, iris))

sim2 <- replicate(R, PermuteFunction(iris, c("Sepal.Width", "Petal.Length"), "Sepal.Length"))
rowMeans(sim2 >= initial)
# (Intercept)  Sepal.Width Petal.Length 
#       0.562        0.461        0.500