获取通过重采样计算的多个回归系数值

Question

我正在使用重采样计算多个线性模型的系数。之前我使用的是 boot 函数，但以后需要在分析中包含新的统计信息，所以我认为这种方式更好。一个可复制的例子：

iris <- iris[,1:4]

nboots <- 100
ncol = ncol(iris)
boot.r.squared <- numeric(nboots)      
boot.p_model <- numeric(nboots)
boot.coef_p <- numeric(nboots)
boot.coef_estimate <- matrix(nrow= nboots,ncol=ncol)
boot.coef_error <- matrix(nrow= nboots,ncol=ncol)
boot.coef_t <- matrix(nrow= nboots,ncol=ncol)
boot.coef_p <- matrix(nrow= nboots,ncol=ncol)

for(i in 1:nboots){
  boot.samp <- iris[sample(nrow(iris),size = 100, replace = TRUE,), ] 
  model <- lm(boot.samp$Sepal.Length ~ .,boot.samp)
  model.sum <- summary(model)

  boot.r.squared[i] <- model.sum$r.squared
  stat <- model.sum$fstatistic
  boot.p_model[i] <- pf(stat[1], stat[2], stat[3], lower.tail = FALSE)

  boot.coef_estimate[i, 1:length(model$coefficients)] <- model$coefficients[1]
  boot.coef_error[i, 1:length(model$coefficients)] <- model$coefficients[2]
  boot.coef_t[i, 1:length(model$coefficients)] <- model$coefficients[3]
  boot.coef_p[i, 1:length(model$coefficients)] <- model$coefficients[4] 
}

但是我无法以矩阵的形式正确保存系数。我想保存 4 个矩阵，每个矩阵包含：参数、误差、统计 t 和 p 值。

使用此代码，所有列都相同。我试着用 put [1] 来保存第一列，但发生了这个错误。我该如何解决这个错误？

Error in model $ coefficients [, 1]: incorrect number of dimensions

Answer 1

你的代码对我来说没有错误。您可能尝试了不同的重复次数并且忘记更新您在循环之前定义的空对象的大小，因为当我更改为 nboots <- 1000.

时我可以重现您的错误

您可以在没有 for 循环的情况下执行此操作，而是使用 replicate()（这基本上是 lapply() 的包装器）。优点是它可能更直接一些，您不需要事先定义空对象。

过程是，first，定义一个bootstrap函数，比如说bootFun()，second，实际 bootstrap 使用 replicate() 生成一个数组，最后 third，将数组中的数据处理成所需的矩阵。

# Step 1 bootstrap function
bootFun <- function(X) {
  fit <- lm(Sepal.Length ~ ., data=X)
  s <- summary(fit)
  r2 <- s$r.squared
  f <- s$fstatistic
  p.f <- pf(f[1], f[2], f[3], lower.tail = FALSE)
  ms <- s$coefficients
  return(cbind(ms, r2, p.f))
}

# Step 2 bootstrap
nboots <- 1e2
set.seed(42)  # for sake of reproducibility
A <- replicate(nboots, bootFun(iris[sample(1:nrow(iris), size=nrow(iris), replace=TRUE), ]))
# Note: I used here `size=nrow(iris)`, but you also could use size=100 if this is your method

# Step 3 process array ("de-transpose" and transform to list)
Ap <- aperm(A, c(3, 1, 2))  # aperm() is similar to t() for matrices
Aps <- asplit(Ap, 3)  # split the array into a handy list

# or in one step
Aps <- asplit(aperm(A, c(3, 1, 2)), 3)

结果

现在你有一个包含所有矩阵的列表，查看列表对象的名称。

names(Aps)
# [1] "Estimate"   "Std. Error" "t value"    "Pr(>|t|)"   "r2"         "p.f"

所以，如果我们例如想要我们估计的矩阵，我们简单地做：

estimates <- Aps$Estimate

estimates[1:3, ]
#      (Intercept) Sepal.Width Petal.Length Petal.Width
# [1,]    1.353531   0.7655760    0.8322749  -0.7775090
# [2,]    1.777431   0.6653308    0.7353491  -0.6024095
# [3,]    2.029428   0.5825554    0.6941457  -0.4795787

请注意，F-统计量和 R² 也是矩阵，并且为每个系数在每一行中重复。由于您只需要一个，只需提取例如第一栏：

Aps$p.f[1:3, ]
#       (Intercept)  Sepal.Width Petal.Length  Petal.Width
# [1,] 2.759019e-65 2.759019e-65 2.759019e-65 2.759019e-65
# [2,] 5.451912e-66 5.451912e-66 5.451912e-66 5.451912e-66
# [3,] 3.288712e-54 3.288712e-54 3.288712e-54 3.288712e-54

Aps$p.f[1:3, 1]
# [1] 2.759019e-65 5.451912e-66 3.288712e-54

基准

这两种方法的速度几乎相同，replicate() 略有优势。这是比较我使用 nboot=1,000 复制的两种方法的基准。

# Unit: seconds
#      expr      min       lq     mean   median       uq      max neval cld
#   forloop 2.215259 2.235797 2.332234 2.381035 2.401933 2.700622   100   a
# replicate 2.218291 2.240570 2.313526 2.257905 2.400532 2.601958   100   a

Answer 2

关于您的原始代码，您混淆了模型和 model.sum。 slot in coefficients, error等需要得到summary(model)的系数

见下文：

for(i in 1:nboots){
  boot.samp <- iris[sample(nrow(iris),size = 100, replace = TRUE,), ] 
  model <- lm(boot.samp$Sepal.Length ~ .,boot.samp)
  model.sum <- summary(model)

  boot.r.squared[i] <- model.sum$r.squared
  stat <- model.sum$fstatistic
  boot.p_model[i] <- pf(stat[1], stat[2], stat[3], lower.tail = FALSE)

  # this is the part where you need to use model.sum
  boot.coef_estimate[i, 1:length(model$coefficients)] <- model.sum$coefficients[,1]
  boot.coef_error[i, 1:length(model$coefficients)] <- model.sum$coefficients[,2]
  boot.coef_t[i, 1:length(model$coefficients)] <- model.sum$coefficients[,3]
  boot.coef_p[i, 1:length(model$coefficients)] <- model.sum$coefficients[,4] 
}

作为旁注，您想知道使用 boot 的系数估计值有多大变化。对于其他统计数据，如 r 平方、p 值，在 bootstrap 之后解释这些并不是很有用。所以想想你想从 bootstrap.

中了解什么

另请查看@jay.sf 的回答，这是包装计算函数的更简洁的方法..

获取通过重采样计算的多个回归系数值

Get multiple regression coefficient values calculated by resampling

r

regression

resampling

statistics-bootstrap

结果

基准