问题描述
我的问题是双重的:1.如下所示,我试图基于两个变量对子集进行嵌套循环,然后执行t.test,然后使用这些结果填充数据框。就目前而言,我的代码仅迭代一个变量,而不是两个变量。我想念的是不允许这样做吗?
- 我了解向量化在这里会有所帮助,但我对此并不熟悉,并且希望获得一些有关如何实现的反馈。
背景:我一直在研究一个小问题,但是我被困住了。我正在尝试通过使用两个变量进行子集分析一些数据。如果我只是想完成它,我将基于第一个变量将其子集化为数据帧,然后使用新的数据帧和第二个变量继续进行分析以进行进一步的子集设置。有一些循环经验,我想我会尝试使用嵌套循环为我完成此任务。我已经能够使我的循环对于单个变量的子集可以很好地工作,并建立一个单独的日期框架,然后将其用于其他目的。但是,当我尝试使用第二个变量时,它不起作用。现在,循环仅创建4个唯一的子集,而理想情况下应产生12个。我认为我显然缺少了一些东西,我曾尝试搜索该论坛和其他几个论坛,但无济于事。
这是我要开始的代码:
set.seed(10)
graphdata1 <-data.frame("RC" = sample(1:500,1000,replace = T),"Gl" = sample(letters[1:3],"CS" = sample(1:4,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应该先运行内部循环直到完成,然后再增加外部循环。由于我想首先在G1的第一个变量上设置子集,然后对CS的每个变量进行分析,因此我认为在内部循环中包含我的相关G1线是谨慎的。当然,这是行不通的,只会生成此数据帧,该数据帧必须完成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)上递增,而我没有得到想要的结果。
我想要此输出,其中将基于对子集的唯一组合的总数,用每个唯一子集的平均值和CI填充所有12行(下表是一个示例,理想情况下,将n填充为数字,表示前4行的上下CI):
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.我了解向量化在这里会有所帮助,但我对此并不熟悉,希望能获得一些有关如何实现的反馈。
谢谢
达斯汀
解决方法
对您的代码的评论
首先,关于循环,由于调用了错误的索引,因此无法填充数据框。例如:
for(j in 1:3){
for(i in 1:4){
results[j] <- something[j]
}
}
在这种情况下,j
只会在1到3之间循环,每次出现内部循环时都会重写以前的结果(换句话说,您将在results[1]
中写3次,3 results[2]
中的时间,...)。您要做的就是遵循这些原则:
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()
都没有用,也没有任何作用,并且在循环不是一个好主意。
对于第二个,我们的想法是数据帧通常以连续范围存储在内存中(与向量或矩阵相同)。当添加一行时,您需要在已经存储数据帧的地方添加一些内容,如果没有空间,整个数据帧将被复制,这很慢且效率低下。使用for
循环时,您总是想用正确的长度来初始化结果变量:
N <- 12 #the length you want
results <- data.frame(n = rep(NA,N),ameans = rep(NA,CIameanslower = rep(NA,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中,这很少引起关注,因为我们通常可以将操作向量化。
如何矢量化
您可以使用base R进行所有操作,但是为什么不使用可用的最佳工具:使用tidyverse(尤其是软件包dplyr
)会容易得多。
library(tidyverse)
现在我们可以变换原始数据帧了。
graphdata1 %>%
group_by(Gl,CS) %>%
summarize(mean_RC = mean(RC),sd_RC = sd(RC),n = n())
因此,我们很容易得出观测值的平均值,标准差和数量;您可以在此处添加任何摘要统计信息。 但是你想做一个t检验。如果我理解正确,则需要进行一次样本检验,将样本中的均值与0进行比较。您可以尝试将其简单地添加到摘要中:
graphdata1 %>%
group_by(Gl,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,res = list(t.test(RC)),.groups="drop")
我还添加了.groups="drop"
,以避免以后进行分组,这可能会影响后续操作。
剩下要做的就是从存储的测试结果中提取感兴趣的值。还有一个技巧:我们需要指定要使用rowwise()
逐行而不是逐列进行计算。
graphdata1 %>%
group_by(Gl,.groups="drop") %>%
rowwise() %>%
mutate(lower.ci = res$conf.int[1],upper.ci = res$conf.int[2])
我们完成了!我们可以使用select()
删除不再有用的列,并对其重命名和排序,并使用arrange()
用1个或多个变量对行进行排序。
graphdata1 %>%
group_by(Gl,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.vectors,并从循环结构中删除了添加行功能。
-
为了最佳实践,我将数据框放在循环之外。
set.seed(10) graphdata1 <-data.frame("RC" = sample(1:500,1000,replace = T),"Gl" = sample(letters[1:3],"CS" = sample(1:4,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,ncol=6)) names(results) <- c("n","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
再次感谢!