嵌套循环 - 通过对另外两个变量进行子集化来分析一个变量

Question

我的问题是双重的：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 行都将填充每个唯一子集的平均值和 CIs，基于要子集的唯一组合的总数（下面的 table 是一个例如，理想情况下会为 n、ameans、upper 和 lower 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. 我知道矢量化在这里会有所帮助，但我对此并不熟悉，希望得到一些关于如何实施的反馈。

谢谢

达斯汀

Answer 1

对您的代码的评论

首先，关于您的循环，它无法填充数据框，因为您调用了错误的索引。例如：

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]，...，当i=0,j=1你写成 results[4]，...，当 i=3,j=2 写成 results[12]。这足以使循环执行您想要的操作。

此外，还有两件小事不是最佳实践但不应该影响结果：我认为你所有的 as.vector() 都没有用并且没有效果，以及向数据框添加行在循环期间不是一个好主意。

对于第二个，想法是数据帧通常存储在内存中的连续范围内（对于向量或矩阵也是如此）。当你添加一行时，你需要在数据框已经存储的地方附加一些东西，如果没有 space 整个数据框将被复制，这是缓慢且低效的。使用 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())

所以我们很容易得到平均数、标准差和观察次数；您可以在此处添加任何摘要统计信息。但是您想进行 t 检验。如果我理解正确的话，你想要一个 one-sample 测试，将样本中的平均值与 0 进行比较。你可以尝试简单地将它添加到 summarize:

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".

没用。但是看看错误信息：测试成功了，但是你不能只把测试的结果放在数据框中。一个魔术是使用“list-column”：我们的数据框的其中一列将是一个列表，它可以包含任何内容，甚至是整个测试结果。

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.

Answer 2

感谢@Alexlok 的帮助。查看答案后，我将使用矢量化，因为它效率更高。为了完成，我想我会根据建议 post 我的新嵌套循环代码。改进：

我调用了正确的索引：(j-1)*3+i+(j-1) 我发现我需要将“+(j-1)”项添加到索引中以防止循环覆盖自身。
我摆脱了 as.vectors 并从循环结构中删除了添加行函数。

为了最佳实践，我在循环外制作了数据框。

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

再次感谢！

嵌套循环 - 通过对另外两个变量进行子集化来分析一个变量

Nested loop - analysis of one variable by subsetting on two others variables

r

vectorization

nested-loops

对您的代码的评论

如何矢量化