我的问题有两个：1。如下所示，我尝试基于两个变量对子集进行嵌套循环，然后执行t.test，然后用这些结果填充数据帧。目前，我的代码只对一个变量进行迭代，而不是对两个变量都进行迭代。我错过了什么不让它发挥作用？
- 我知道矢量化在这里会有所帮助，但我不熟悉这一点，希望能就如何实现这一点提供一些反馈
背景：我处理一个小问题有一段时间了，我被卡住了。我试图通过使用两个变量进行子集设置来分析一些数据。如果我只是想完成它，我只需要根据第一个变量将其子集划分为数据帧，然后使用新的数据帧和第二个变量继续我的分析，以进行进一步的子集划分。有了一些循环的经验，我想我会尝试使用嵌套循环来为我做这件事。我已经能够让我的循环很好地为一个变量子集工作，并构建一个单独的日期框架，然后我可以将其用于其他目的。然而，当我尝试使用第二个变量时，它不起作用。现在，循环只创建4个唯一的子集，而理想情况下应该生成12个。我想我显然遗漏了一些东西，我试过搜索这个论坛和其他几个论坛，但都没有用。
这是我的启动代码：
set.seed(10)
graphdata1 <-data.frame("RC" = sample(1:500, 1000, replace = T), "Gl" = sample(letters[1:3], 1000, replace = T), "CS" = sample(1:4, 1000, replace = T))
responsesGl <- as.vector(levels(as.factor(graphdata1$Gl)))
results <- data.frame("n"=0, "ameans"=0, "CIameanslower"=0, "CIameansupper"=0)
results$Gl<- NA
results$CS <-NA
responsesCS <- as.vector(levels(as.factor(graphdata1$CS)))
for(j in 1:length(responsesGl)) {
for(i in 1:length(responsesCS)) {
results$Gl[j] <- responsesGl[j] #adds in the first subsetting variable to the dataframe
y <- subset(graphdata1, Gl == responsesGl[j]) #creates a subsetted dataframe of the larger data to analyze
results$CS[i] <- responsesCS[i] #adds in the second subsetting variable
x <- subset(y, CS == responsesCS[i]) #further subsets data to obtain only data that is a based on first and second variables
results$n[i] <-length(x$CS) #determines number of responses in this category
ttest <- t.test(x$RC) #this and the next four lines all analyze the data, while amending the analysis to the results dataframe
confidence_interval <- as.vector(unlist(ttest["conf.int"]))
results$ameans[i] <- mean(x$RC, na.rm = TRUE)
results$CIameanslower[i] <- confidence_interval[1]
results$CIameansupper[i] <- confidence_interval[2]
if (length(results$n) == length(responsesCS)*length(responsesGl)) { #adds a row if the results sheet is not as long as the product of the response vectors (12 in this case)
rm(x)
rm(y)} else {
results[nrow(results)+1,] <- NA #adds a row
rm(x)
rm(y)
}
}
}
从我的搜索中，我想我明白了R应该首先运行内部循环直到完成，然后增加外部循环。由于我想先在Gl的第一个变量上子集，然后分析CS的每个变量，我认为在内环中包括我的相关Gl线是谨慎的。当然，它不起作用，只生成4行已完成但8行为空的数据帧(总共12行)：
n ameans CIameanslower CIameansupper Gl CS
1 95 247.7579 218.2211 277.2947 a 1
2 84 257.3929 224.1692 290.6165 b 2
3 88 257.7500 226.3831 289.1169 c 3
4 68 244.8971 206.5598 283.2343 <NA> 4
5 NA NA NA NA <NA> <NA>
6 NA NA NA NA <NA> <NA>
7 NA NA NA NA <NA> <NA>
8 NA NA NA NA <NA> <NA>
9 NA NA NA NA <NA> <NA>
10 NA NA NA NA <NA> <NA>
11 NA NA NA NA <NA> <NA>
12 NA NA NA NA <NA> <NA>
我意识到内部循环也在第一个变量(Gl)上递增，我没有得到我想要的。
我想要这样的输出，其中所有12行都将根据子集的唯一组合总数，用每个唯一子集的平均值和CI填充(下表是一个示例，理想情况下，将为n、ameans、上下CI填充数字，如前4行所示)：
n ameans CIameanslower CIameansupper Gl CS
1 95 247.7579 218.2211 277.2947 a 1
2 84 257.3929 224.1692 290.6165 a 2
3 88 257.7500 226.3831 289.1169 a 3
4 68 244.8971 206.5598 283.2343 a 4
5 NA NA NA NA b 1
6 NA NA NA NA b 2
7 NA NA NA NA b 3
8 NA NA NA NA b 4
9 NA NA NA NA c 1
10 NA NA NA NA c 2
11 NA NA NA NA c 3
12 NA NA NA NA c 4
只是重申我的问题：1。我错过了什么不让它发挥作用？2.我知道矢量化在这里会有所帮助，但我不熟悉这一点，希望能就如何实现这一点提供一些反馈。
谢谢
Dustin
对代码的评论
首先，关于循环，它无法填充数据帧，因为您调用了错误的索引。例如：
for(j in 1:3){
for(i in 1:4){
results[j] <- something[j]
}
}
在这种情况下，j将只在1和3之间循环，在每次出现内部循环时重写之前的结果(换句话说，你将在results[1]中写3次，在results[2]中写3遍，…)
for(j in 0:2){
for(i in 0:3){
results[j*3 + i + 1] <- something[j]
}
}
所以当i=j=0时，你用result[1]写，当i=1,j=0时，你写results[2]。。。，当你用results[4]写i=0,j=1时。。。，当i=3,j=2时，您用results[12]书写。这可能足以使循环执行您想要的操作。
此外，有两件小事不是最佳实践，但不应该影响结果：我认为所有的as.vector()都没有用处，也没有效果，在循环期间向数据帧添加行也不是一个好主意。
对于第二个，其思想是数据帧通常存储在存储器中的连续范围中(对于向量或矩阵也是如此)。当您添加一行时，您需要在数据帧已经存储的位置附加一些内容，如果没有空间，整个数据帧将被复制，这既慢又低效。当使用for循环时，您总是希望用正确的长度初始化结果变量：
N <- 12 #the length you want
results <- data.frame(n = rep(NA, N),
ameans = rep(NA, N),
CIameanslower = rep(NA, N),
CIameansupper = rep(NA, N))
# or an easier equivalent way:
results <- matrix(NA, nrow=N, ncol=4)
results <- as.data.frame(results)
names(results) <- c("n", "ameans", "CIameanslower", "CIameansupper")
但在R中，这很少是一个问题，因为我们通常可以将运算矢量化。
如何矢量化
您可以使用基本R做任何事情，但为什么不使用可用的最佳工具呢？在这里，使用tidyverse(尤其是包dplyr)会更容易。
library(tidyverse)
现在我们可以转换原始数据帧。
graphdata1 %>%
group_by(Gl, CS) %>%
summarize(mean_RC = mean(RC),
sd_RC = sd(RC),
n = n())
因此，我们很容易得到观测值的平均值、sd和数量；您可以在此处添加任何汇总统计信息。但是你想做一个t测试。如果我理解正确的话，你想要一个单样本测试，将你的样本中的平均值与0进行比较。你可以试着简单地添加它作为总结：
graphdata1 %>%
group_by(Gl, CS) %>%
summarize(mean_RC = mean(RC),
sd_RC = sd(RC),
n = n(),
t_test = t.test(RC))
# Error: Problem with `summarise()` input `t_test`.
# x Input `t_test` must be a vector, not a `htest` object.
# i Input `t_test` is `t.test(RC)`.
# i The error occurred in group 1: Gl = "c", CS = "1".
它不起作用。但请看一下错误消息：测试成功了，但不能只将测试结果放在数据帧中。一个魔术是使用一个"；列表列"：数据框架的其中一列将是一个列表，它可以包含任何内容，甚至可以包含整个测试结果。
graphdata1 %>%
group_by(Gl, CS) %>%
summarize(mean_RC = mean(RC),
sd_RC = sd(RC),
n = n(),
res = list(t.test(RC)),
.groups="drop")
我还添加了.groups="drop"，以避免以后进行可能影响后续操作的分组。
我们所要做的就是从存储的测试结果中提取感兴趣的值。还有一个技巧：我们需要指定要使用rowwise()逐行而不是逐列进行计算。
graphdata1 %>%
group_by(Gl, CS) %>%
summarize(mean_RC = mean(RC),
sd_RC = sd(RC),
n = n(),
res = list(t.test(RC)),
.groups="drop") %>%
rowwise() %>%
mutate(lower.ci = res$conf.int[1],
upper.ci = res$conf.int[2])
我们完了！我们可以使用select()删除不再感兴趣的列，并重命名和排序要保留的列，使用arrange()按1个或多个变量对行进行排序。
graphdata1 %>%
group_by(Gl, CS) %>%
summarize(mean_RC = mean(RC),
sd_RC = sd(RC),
n = n(),
res = list(t.test(RC)),
.groups="drop") %>%
rowwise() %>%
mutate(lower.ci = res$conf.int[1],
upper.ci = res$conf.int[2]) %>%
select(Gl, CS, mean_RC,
conf_low = lower.ci, conf_high = upper.ci) %>%
arrange(rev(Gl), CS)
# Gl CS mean_RC conf_low conf_high
# <fct> <fct> <dbl> <dbl> <dbl>
# 1 a 1 213. 181. 245.
# 2 a 2 225. 190. 260.
# 3 a 3 257. 229. 285.
# 4 a 4 221. 184. 257.
# 5 b 1 242. 214. 270.
# 6 b 2 255. 222. 288.
# 7 b 3 225. 196. 255.
# 8 b 4 236. 207. 264.
# 9 c 1 248. 218. 277.
# 10 c 2 257. 224. 291.
# 11 c 3 258. 226. 289.
# 12 c 4 245. 207. 283.
感谢@Alexlok的帮助。看完答案后，我将使用矢量化，因为它的效率要高得多。为了完成，我想我会根据建议发布我的新嵌套循环代码。改进：
-
我调用了正确的索引，使用：(j-1)*3+I+(j-1我发现我需要添加"+(j-1)"；索引的术语，以防止循环自我书写。
-
我去掉了as.vector，并从循环结构中删除了add-rows函数。
-
为了最佳实践，我在循环之外制作了数据帧。
set.seed(10) graphdata1 <-data.frame("RC" = sample(1:500, 1000, replace = T), "Gl" = sample(letters[1:3], 1000, replace = T), "CS" = sample(1:4, 1000, replace = T)) #got rid of as.vector() responsesGl <- levels(factor(graphdata1$Gl)) responsesCS <- levels(factor(graphdata1$CS)) #Create the data frame outside the loop. N <- length(responsesCS)*length(responsesGl) results <- as.data.frame(matrix(NA, nrow=N, ncol=6)) names(results) <- c("n", "ameans", "CIameanslower", "CIameansupper", "Gl", "CS") #The nested loop function. for(j in 1:length(responsesGl)) { for(i in 1:length(responsesCS)) { results$Gl[(j-1)*3+i+(j-1)] <- responsesGl[j] y <- subset(graphdata1, Gl == responsesGl[j]) results$CS[(j-1)*3+i+(j-1)] <- responsesCS[i] x <- subset(y, CS == responsesCS[i]) results$n[(j-1)*3+i+(j-1)] <-length(x$CS) ttest <- t.test(x$RC) confidence_interval <- as.vector(unlist(ttest["conf.int"])) results$ameans[(j-1)*3+i+(j-1)] <- mean(x$RC, na.rm = TRUE) results$CIameanslower[(j-1)*3+i+(j-1)] <- confidence_interval[1] results$CIameansupper[(j-1)*3+i+(j-1)] <- confidence_interval[2] rm(x) rm(y) }}
这是输出：
n ameans CIameanslower CIameansupper Gl CS
1 89 212.8202 181.0133 244.6271 a 1
2 77 224.8961 190.0473 259.7449 a 2
3 95 256.9895 229.0892 284.8897 a 3
4 68 220.5147 183.9511 257.0783 a 4
5 90 242.1667 214.4563 269.8770 b 1
6 75 254.9467 221.7683 288.1250 b 2
7 90 225.4333 195.6203 255.2463 b 3
8 81 235.7037 207.3833 264.0241 b 4
9 95 247.7579 218.2211 277.2947 c 1
10 84 257.3929 224.1692 290.6165 c 2
11 88 257.7500 226.3831 289.1169 c 3
12 68 244.8971 206.5598 283.2343 c 4
再次感谢！