我正在使用重采样计算多个线性模型的系数。之前我使用的是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
最佳答案
你的代码对我来说没有错误。您可能尝试了不同次数的重复,但忘记更新您在循环之前定义的空对象的大小,因为当我更改为 nboots <- 1000
时,我可能会重现您的错误。 .
无需 for
即可做到这一点循环并使用replicate()
相反(它基本上是 lapply()
的包装)。优点是它可能更简单一些,并且您不需要事先定义空对象。
该过程首先定义一个引导函数,例如 bootFun()
,第二,使用 replicate()
的实际 Bootstrap 它产生一个数组,最后第三,将数组中的数据处理成所需的矩阵。
# 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 统计量和 R2 也是矩阵,并且在每个系数的每一行中都有重复。由于您只需要一个,因此只需提取例如第一列:
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
关于r - 获取通过重采样计算的多重回归系数值,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/58830576/