logo标签列表

在R中进行pairwise.t.test的数据操作

rdata-manipulation
3

我正在尝试制作一个成对比较表，并保留每个成对的p值。需要注意的是，我还是R的初学者。我的数据看起来像这样（尽管更大）：

a <- factor(c("ID1","ID2","ID3","ID4","ID5"))
b <- runif(5)
b1 <- runif(5)
b2 <- runif(5)
b3 <- runif(5)
c1 <- runif(5)
c2 <- runif(5)
c3 <- runif(5)
df <- data.frame(a,b1,b2,b3,c1,c2,c3)

b1、b2、b3应该与每一行中的c1、c2、c3进行比较（对于a列中的每个ID）。

最终结果应该类似于：

a <- cbind(a,Adjusted_P_Values)

head(a,1)的例子如下：

head(a,1)
    a        b1        b2        b3        c1        c2
1 ID1 0.1337694 0.7347543 0.5808391 0.4324976 0.5378458
         c3        Adjusted_P_value
1 0.6368778        0.99

每行都有其相应的P值。 我找到了一个可以胜任的函数 pairwise.t.test。 （目前，我只是为每行运行一个循环，并进行普通的t检验，然后用p.adjust进行校正，但我无法进行汇总的sd——这是我想要的。）

所以我的问题现在是如何构建数据，使得R能够接受它。我可以使用reshape2库中的函数，但它不会给我正确的结构。 我像这样使用它：

Test_Data <- melt(df, "a", c("b1","b2","b3","c1","c2","c3"))

但是我失去了行对称性。 因为现在当我使用pairwise.t.test时，我必须使用由melt创建的“a”列或“variable”列之一，因此我要么比较复制品，要么比较ID。 因此，我的问题很简单： 如何构造数据以便测试每一行，并为每一行获取p值，每个处理（b或c）都有一个基于所有行的标准差（所有b的标准差和所有c的标准差各一个）？ 我已经搜索了很多类似问题（以及关于pairwise.t.test的教程），但没有成功。
尝试使用 lapply(split(Test_Data, Test_Data$a), function(TD) t.test(value ~ grepl("b", variable), TD))。 - Rui Barradas
3个回答
3

我的方法与其他答案略有不同，通过时间度量将数据分成两列 b 和c (1-3)，然后使用t.test(...,paired=TRUE)进行成对t检验。

set.seed(1234)
a <- factor(c("ID1","ID2","ID3","ID4","ID5"))
b <- runif(5)
b1 <- runif(5)
b2 <- runif(5)
b3 <- runif(5)
c1 <- runif(5)
c2 <- runif(5)
c3 <- runif(5)
df <- data.frame(a,b1,b2,b3,c1,c2,c3)
library(tidyr)
library(dplyr)
df %>% 
 gather(.,key="variable",value="value",-a) %>%
     extract(.,variable,into = c("measure", "time"), 
             regex = "([A-Za-z]+)([0-9]+)") %>%
      spread(.,measure,value) -> spreadData
# split by ID to conduct paired t-tests by ID
dataList <- split(spreadData,spreadData$a)
pValues <- unlist(lapply(dataList,function(x){
   t.test(x$b,x$c,paired=TRUE)$p.value
}))
df$p.value <- pValues
df

...并输出：

> df
    a          b1        b2        b3         c1         c2
1 ID1 0.640310605 0.6935913 0.8372956 0.31661245 0.81059855
2 ID2 0.009495756 0.5449748 0.2862233 0.30269337 0.52569755
3 ID3 0.232550506 0.2827336 0.2668208 0.15904600 0.91465817
4 ID4 0.666083758 0.9234335 0.1867228 0.03999592 0.83134505
5 ID5 0.514251141 0.2923158 0.2322259 0.21879954 0.04577026
         c3   p.value
1 0.4560915 0.3391364
2 0.2651867 0.5043753
3 0.3046722 0.4598274
4 0.5073069 0.6764142
5 0.1810962 0.1178471
>

注意：如果将其他答案中的代码修改以包括paired=TRUE参数，则两个解决方案中的p值匹配。

另一种方法：在c和b之间运行t检验的差异

考虑到这篇文章关于成对t检验的评论，我想说明一下成对测试中发生了什么。基本上，对于每个时间段1-3，我们从c值中减去b值，并对差异运行t检验。由于我们将数据缩减为单列，因此不需要paired=参数，但是与通过向t.test()传递2个列和paired=TRUE参数相同的结果。

# alternative 2: subtract b from c and use regular t-test
# to show how pairwise works
spreadData$difference <- spreadData$c - spreadData$b
dataList <- split(spreadData,spreadData$a)
pValues <- unlist(lapply(dataList,function(x){
     t.test(x$difference)$p.value
}))
df$p.value <- pValues
df

...和输出：

> spreadData$difference <- spreadData$c - spreadData$b
> dataList <- split(spreadData,spreadData$a)
> pValues <- unlist(lapply(dataList,function(x){
+      t.test(x$difference)$p.value
+ }))
> df$p.value <- pValues
> df
    a          b1        b2        b3         c1         c2
1 ID1 0.640310605 0.6935913 0.8372956 0.31661245 0.81059855
2 ID2 0.009495756 0.5449748 0.2862233 0.30269337 0.52569755
3 ID3 0.232550506 0.2827336 0.2668208 0.15904600 0.91465817
4 ID4 0.666083758 0.9234335 0.1867228 0.03999592 0.83134505
5 ID5 0.514251141 0.2923158 0.2322259 0.21879954 0.04577026
         c3   p.value
1 0.4560915 0.3391364
2 0.2651867 0.5043753
3 0.3046722 0.4598274
4 0.5073069 0.6764142
5 0.1810962 0.1178471
>
我打赌这也是一个非常好的回答（如果我声望更高，我会点赞），尽管我很难理解正在发生什么；我需要在我的虚拟数据集上稍微试验一下。从两个答案中可以清楚地看出，我需要更多地了解这个tidyr包。一个问题浮现在脑海中，就是如何使用regex函数选择特定的列名？ - Baraliuh
@Baraliuh - 没问题，欢迎来到StackOverflow。输入和输出列名在extract()函数中指定，而regex=参数指定了一个正则表达式，用于将输入列中的数据拆分成into=参数中指定的输出列。有关Tidyverse的更多信息，请访问Tidyverse网站。 - Len Greski
2
一种可能的解决方案是使用软件包。
首先，将数据框的格式调整为以下结构。
library(tidyverse)

df2 <- df %>%
  gather(Column, Value, -a) %>%
  extract(Column, into = c("Group", "Number"), regex = "([A-Za-z]+)([0-9]+)")
df2
#      a Group Number       Value
# 1  ID1     b      1 0.640310605
# 2  ID2     b      1 0.009495756
# 3  ID3     b      1 0.232550506
# 4  ID4     b      1 0.666083758
# 5  ID5     b      1 0.514251141
# 6  ID1     b      2 0.693591292
# 7  ID2     b      2 0.544974836
# 8  ID3     b      2 0.282733584
# 9  ID4     b      2 0.923433484
# 10 ID5     b      2 0.292315840
# 11 ID1     b      3 0.837295628
# 12 ID2     b      3 0.286223285
# 13 ID3     b      3 0.266820780
# 14 ID4     b      3 0.186722790
# 15 ID5     b      3 0.232225911
# 16 ID1     c      1 0.316612455
# 17 ID2     c      1 0.302693371
# 18 ID3     c      1 0.159046003
# 19 ID4     c      1 0.039995918
# 20 ID5     c      1 0.218799541
# 21 ID1     c      2 0.810598552
# 22 ID2     c      2 0.525697547
# 23 ID3     c      2 0.914658166
# 24 ID4     c      2 0.831345047
# 25 ID5     c      2 0.045770263
# 26 ID1     c      3 0.456091482
# 27 ID2     c      3 0.265186672
# 28 ID3     c      3 0.304672203
# 29 ID4     c      3 0.507306870
# 30 ID5     c      3 0.181096208

第二步，将数据框分割并进行pairwise.t.test，然后提取P值。
p_value <- df2 %>%
  split(.$a) %>%
  map(function(x) pairwise.t.test(x$Value, x$Group, paired = TRUE)) %>%
  map_dbl("p.value")
p_value
#       ID1       ID2       ID3       ID4       ID5 
# 0.3391364 0.5043753 0.4598274 0.6764142 0.1178471

最后，将P值作为新列添加到原始数据框中。
df_final <- df %>% mutate(Adjusted_P_value = p_value)
df_final
#     a          b1        b2        b3         c1         c2        c3 Adjusted_P_value
# 1 ID1 0.640310605 0.6935913 0.8372956 0.31661245 0.81059855 0.4560915        0.3391364
# 2 ID2 0.009495756 0.5449748 0.2862233 0.30269337 0.52569755 0.2651867        0.5043753
# 3 ID3 0.232550506 0.2827336 0.2668208 0.15904600 0.91465817 0.3046722        0.4598274
# 4 ID4 0.666083758 0.9234335 0.1867228 0.03999592 0.83134505 0.5073069        0.6764142
# 5 ID5 0.514251141 0.2923158 0.2322259 0.21879954 0.04577026 0.1810962        0.1178471

DATA

set.seed(1234)

a <- factor(c("ID1","ID2","ID3","ID4","ID5"))
b <- runif(5)
b1 <- runif(5)
b2 <- runif(5)
b3 <- runif(5)
c1 <- runif(5)
c2 <- runif(5)
c3 <- runif(5)
df <- data.frame(a,b1,b2,b3,c1,c2,c3)

编辑：

为了正确地将P值映射回数据框，数据框必须按照“a”列进行排序。

1
好的，我认为这正是我想要的。需要阅读和学习一些细节。有一件事我不明白，就是如何选择要测试哪些列？因为我简化了我发布的数据集，需要更加有选择性（我有比b和c更多的组）。我会接受这个答案，因为我可以根据这个解释制作数据框架！:)非常感谢。 - Baraliuh
@LenGreski 感谢您的评论。我已经更新了我的帖子，将 paired=TRUE 设置为真。 - www
1
@www - pairwise.t.test() 的参数有些奇怪，因为人们会期望该函数的默认值为paired=TRUE。我正要发布我的答案，但是我注意到你在网站上发布了帖子，所以不得不弄清楚为什么我们的答案中的p值不匹配。 - Len Greski
1
@LenGreski 我认为 pairwise.t.test() 中的 "pairwise" 意味着进行多重比较时组水平之间的成对比较，而不是每一对之间的比较。这有点令人困惑。这是我第一次学习这个函数，感谢 OP 的问题。再次感谢。 - www
1
我得承认，我有些“老派”，多年前我学习如何进行成对 t 检验时，是通过将数据分成两列，每个观察值都有一个“之前”和“之后”的列，然后通过相减并对结果进行 t 检验来计算的。 - Len Greski
显示剩余3条评论
0

只需添加到Baraliuh的解决方案中：

map_dbl("p.value")不起作用，但是在我的情况下，map_df("p.value")可以使用

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接