让我们通过不同的方式来想象做这种事情。
首先,要知道for 循环有它们的位置,如果做得好,它们可以和*apply 函数一样快。尽管您对循环的使用在语法上是正确的,但有不同的使用方式可能更有意义。您正在尝试对列表的多个元素运行一系列命令或函数。想象一下这个简单的计划:对于列表中的每个元素,取第一个元素,然后将其加倍并平方:
invec <- list(c(21,22),c(23,24),c(25,26))
str(invec)
# List of 3
# $ : num [1:2] 21 22
# $ : num [1:2] 23 24
# $ : num [1:2] 25 26
outvec <- replicate(length(invec), NULL) # preallocate same size
for (i in seq_along(invec)) {
outvec[[i]] <- c(2*invec[[i]][1], invec[[i]][1]^2)
}
str(outvec)
# List of 3
# $ : num [1:2] 42 441
# $ : num [1:2] 46 529
# $ : num [1:2] 50 625
看起来很简单。现在让我们看看如何用*apply 函数做同样的事情:
invec <- list(c(21,22),c(23,24),c(25,26))
outvec <- lapply(invec, function(a) c(2*a[1], a[1]^2))
str(outvec)
# List of 3
# $ : num [1:2] 42 441
# $ : num [1:2] 46 529
# $ : num [1:2] 50 625
apply函数的读取方式是“获取向量invec,并在每个元素上调用这个函数,将结果捕获到列表名称outvec”。该函数可以是“匿名的”(就像这里一样),也可以是“命名”的函数,例如
lapply(invec, max)
# [[1]]
# [1] 22
#
# [[2]]
# [1] 24
#
# [[3]]
# [1] 26
那么这对您的抽样问题有何帮助?让我再分歧一下。
您是否知道可以任意索引向量和列表?例如:
str(invec[c(1,3,2,3,2,3)])
# List of 6
# $ : num [1:2] 21 22
# $ : num [1:2] 25 26
# $ : num [1:2] 23 24
# $ : num [1:2] 25 26
# $ : num [1:2] 23 24
# $ : num [1:2] 25 26
有骗子,好吧。假设我们想从这个非常短的列表中获取 1000 个随机样本:
set.seed(3)
ind <- sample(length(invec), size=1000, replace=TRUE)
str(outvec[1:4])
# List of 4
# $ : num 42
# $ : num 46
# $ : num 50
# $ : num 46
outvec <- lapply(invec[ind], function(a) 2*a[1])
str(outvec[1:4])
# List of 4
# $ : num 42
# $ : num 50
# $ : num 46
# $ : num 42
好的,所以我们已经对原始列表进行了 1000 次采样并对其进行了处理 (2*a[1]),并存储了结果。
所以让我们将它应用到您的场景中。由于您的数据是肉眼看不见的,我会补一些。
set.seed(2)
n <- 20
lst <- lapply(1:185, function(ign) data.frame(x=sample(100,size=n), y=sample(100,size=n)))
str(lst[1:2])
# List of 2
# $ :'data.frame': 20 obs. of 2 variables:
# ..$ x: int [1:20] 19 70 57 17 91 90 13 78 44 51 ...
# ..$ y: int [1:20] 67 39 83 15 34 47 97 96 89 13 ...
# $ :'data.frame': 20 obs. of 2 variables:
# ..$ x: int [1:20] 99 30 12 16 91 76 92 33 47 74 ...
# ..$ y: int [1:20] 78 88 62 26 83 42 37 43 21 7 ...
现在我有一个包含 185 个数据帧的列表,每个数据帧都有相同的两个变量 x 和 y。让我们将您的问题应用于此数据。哦,随机性可能很耗时。 (顺便说一句:获得 1000 个随机数比 1000 次获得 1 个随机数要快得多。)
ind <- sample(185, size=1000, replace=TRUE)
现在,lst[ind] 将是一个列表,长度为 1000 个元素,每个元素都是从原始列表中随机选择的。
lms <- lapply(lst[ind], function(a) lm(y~x, data=a))
(lm 部分可以是任何你需要的,只要它是应用于所有元素的相同回归。函数中的代码可以根据你的需要,所以也许可以这样想:
lms <- lapply(lst[ind], function(a) {
z <- lm(y~x, data=a)
return(z)
})
这有意义吗?)好的,让我们看看一些输出:
summary(lms[[1]])
# Call:
# lm(formula = y ~ x, data = a)
# Residuals:
# Min 1Q Median 3Q Max
# -53.944 -13.463 -1.239 15.473 44.430
# Coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 80.4523 15.3577 5.239 5.56e-05 ***
# x -0.4217 0.2499 -1.687 0.109
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Residual standard error: 25.23 on 18 degrees of freedom
# Multiple R-squared: 0.1366, Adjusted R-squared: 0.0886
# F-statistic: 2.847 on 1 and 18 DF, p-value: 0.1088
summary(lms[[2]])
# Call:
# lm(formula = y ~ x, data = a)
# Residuals:
# Min 1Q Median 3Q Max
# -55.108 -20.653 -0.465 18.827 42.747
# Coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 60.7651 12.2366 4.966 1e-04 ***
# x -0.1898 0.2060 -0.922 0.369
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Residual standard error: 27.37 on 18 degrees of freedom
# Multiple R-squared: 0.04506, Adjusted R-squared: -0.007996
# F-statistic: 0.8493 on 1 and 18 DF, p-value: 0.3689
“但我不需要整个模型,我只需要系数!” 当然,你是对的。当你知道你只需要一件事时,你显然可以直接进入并直接得到它(例如coef(lm(y~x,data=a)))。所以,与其重新运行 1000 个随机样本的回归,我可以再做一次 lapply:
coefs <- lapply(lms[1:3], coef)
str(coefs[1:3])
# List of 3
# $ : Named num [1:2] 80.452 -0.422
# ..- attr(*, "names")= chr [1:2] "(Intercept)" "x"
# $ : Named num [1:2] 60.77 -0.19
# ..- attr(*, "names")= chr [1:2] "(Intercept)" "x"
# $ : Named num [1:2] 53.716 -0.189
# ..- attr(*, "names")= chr [1:2] "(Intercept)" "x"
在这种情况下,我实际上有几个选择。我可以坚持使用它并将它们“rbind”(行绑定)在一起,
head(do.call(rbind, coefs))
# (Intercept) x
# [1,] 80.45230 -0.42173749
# [2,] 60.76507 -0.18979726
# [3,] 53.71643 -0.18883933
# [4,] 49.51803 0.01494021
# [5,] 49.51803 0.01494021
# [6,] 68.25463 -0.25840920
或者我之前可以使用“简单应用”(可选,但默认为是)将结果简化为矩阵或向量。如果任何返回值的大小与其他值不同,它将始终返回一个列表。 (正因为如此,不简化它,做一些健全性检查,然后 rbind 可能在编程上更具防御性。)
coefs2 <- t(sapply(lms, coef))
head(coefs2)
# (Intercept) x
# [1,] 80.45230 -0.42173749
# [2,] 60.76507 -0.18979726
# [3,] 53.71643 -0.18883933
# [4,] 49.51803 0.01494021
# [5,] 49.51803 0.01494021
# [6,] 68.25463 -0.25840920
请注意,我必须 transpose 输出:输出(没有 t(...))将有 2 个行(每个回归系数一个),这有点古怪和违反直觉) 和 1000 个列。所以我们转置它,因为我很自然地认为它是每模型行。如果您可以将其作为每个模型的列处理,则不需要这样做。
所以最重要的是,您的 for 循环本身在语法上并没有错误,但是如果您考虑以这种方式对许多事情的向量/列表做一件事情,您将获得显着的速度提升(在这种情况下)而且,可以说,一旦你理解了它,代码的可读性就会大大提高。