r - 多个测试

Question

我想对以下格式的数据执行多次测试

第一列是“id”，值（例如）1,1,1,2,2,2

第二列是“比率”，值为 0.2、0.18、0.3、1.5、1.4、1.6

对于“id”的每个实例，我想针对数据框中的所有比率值测试所有比率值

现在我有这个

data <- read.delim("clipboard", stringsAsFactors=FALSE) ##data to test
dist <- as.numeric(readClipboard()) ##distribution to test against


data$Ratio.Mean.H.L <- NA
data$p.value <- NA


for (i in 1:nrow(data))
     if (nrow(data) > 1)
 {
 #welsh t-test
 t.test.result <- t.test(data$ratio[i],dist,
                         alternative = "two.sided",
                         mu = 0, 
                         paired = FALSE, 
                         var.equal = FALSE,
                         conf.level = 0.95)     
 #writes data into the data.frame
 data$p.value[i] <- t.test.result$p.value
 }

write.table(data, file="C:/R_Temp/t-test.txt", sep = "\t")

我知道这不起作用，因为我不确定我只测试共享相同“id”的行。我还手动输入要测试的分布，这是“比率”列中的所有条目。

我该怎么做？并添加多个测试校正（bonferroni）？

score 0 · Accepted Answer

好的，抱歉问题定义不明确。我在其他地方得到了帮助，并将发布对那些感兴趣的人有用的脚本。我想计算蛋白质组学实验中比率变化的 p 值。为此，我对任何给定蛋白质或 PTM 位点的所有比率测量值进行单独的 t 检验。这些测量值与所有测量值的中位数（t.test 函数中的 mu）或测量值的整个分布进行比较。在一个列中，我有每个条目唯一的“id”，在另一列中，我有“值”（比率）。我将进行 t 检验，比较任何给定唯一“id”出现的所有“值”。为了便于使用，我将表格粘贴到脚本中，而不是从文件中调用它（它为我节省了一步）。

data <- read.delim("clipboard", stringsAsFactors=FALSE) ##data to test(two columns "id" and "value") Log-transfrom ratios!!
summary(data)

med <- median(data$value)


# function for the id-grouped t-test 
calc_id_ttest <- function(d) #col1: id, col2:values
  {  
colnames(d) <- c("id", "value")  # reassign the column names

# calculate the number of values for each id
res_N <- as.data.frame(tapply(d$value, d$id, length)) 
colnames(res_N) <- "N"              
res_N$id <- row.names(res_N)

# calculate the number of values for each id
res_med <- as.data.frame(tapply(d$value, d$id, median)) 
colnames(res_med) <- "med"              
res_med$id <- row.names(res_med)

# calculate the pvalues 
res_pval <- as.data.frame(tapply(d$value, d$id, function(x)
{
  if(length(x) < 3) 
    {   # t test requires at least 3 samples
    NA
    } 
  else
    {
    t.test(x, mu=med)$p.value #t.test (Pearson)d$value with other distribution? alternative=less or greater
    }                         #d$value to compare with entire distribution
                              #mu=med for median of values for 1-sided test
}))

colnames(res_pval) <- "pval"  # nominal p value 
res_pval$id <- row.names(res_pval)
res_pval$adj.pval <- p.adjust(res_pval$pval, method = "BH")  #multiple testing correction also "bonferroni"

res <- Reduce(function(x,y)
{
merge(x,y, by = "id", all = TRUE)
}, 
list(res_N, res_med, res_pval))
return (res)
}

data_result <- calc_id_ttest(d = data)
write.table(data_result, file="C:/R_Temp/t-test.txt", quote = FALSE, row.names = FALSE, col.names = TRUE, sep = "\t")

score 0 · Accepted Answer

我怀疑 MattParker 的评论将是这里最重要的事情：您正在将单个数字与向量进行比较，并且t.test会抱怨这一点。由于您建议您要对每个分组变量 ( id) 执行测试，因此在基础 R 中您可能希望使用类似by(或split) 的函数。dplyr（内部也有很好的方法data.table。）

使用mtcars作为样本数据，我将尝试模仿您的数据：

dat <- mtcars[c("cyl", "mpg")]
colnames(dat) <- c("id", "ratio")

不清楚你的意思是用于dist，所以我会使用天真的

dist <- 1:10

现在你可以这样做：

by(dat$ratio, dat$id, function(x) t.test(x, dist, paired = FALSE)$p.value)
# dat$id: 4
# [1] 2.660716e-10
# ------------------------------------------------------------ 
# dat$id: 6
# [1] 4.826322e-09
# ------------------------------------------------------------ 
# dat$id: 8
# [1] 2.367184e-07

如果您想/需要处理的不仅仅是ratio一次，您也可以这样做：

by(dat, dat$id, function(x) t.test(x$ratio, dist, paired = FALSE)$p.value)
# dat$id: 4
# [1] 2.660716e-10
# ------------------------------------------------------------ 
# dat$id: 6
# [1] 4.826322e-09
# ------------------------------------------------------------ 
# dat$id: 8
# [1] 2.367184e-07

调用的结果by是一个 class "by"，它实际上只是一个list带有一些额外属性的重新包装：

res <- by(dat, dat$id, function(x) t.test(x$ratio, dist, paired = FALSE)$p.value)
class(res)
# [1] "by"
str(attributes(res))
# List of 4
#  $ dim     : int 3
#  $ dimnames:List of 1
#   ..$ dat$id: chr [1:3] "4" "6" "8"
#  $ call    : language by.data.frame(data = dat, INDICES = dat$id, FUN = function(x) t.test(x$ratio,      dist, paired = FALSE)$p.value)
#  $ class   : chr "by"

因此，您可以扩展/访问它，但是您可以list：

res[[1]]
# [1] 2.660716e-10
as.numeric(res)
# [1] 2.660716e-10 4.826322e-09 2.367184e-07
names(res)
# [1] "4" "6" "8"

（意识到的不同级别dat$id是整数 4、6 和 8，所以names应该对应于您的$id。）

编辑：

如果您想要 data.frame 中的结果，可以想到两个选项：

对每一行重复 p 值，导致大量重复。我不鼓励这种方法有几个原因；如果您在某些时候需要它，我建议使用选项 2，然后merge.

生成一个具有与 unique 一样多的行的 data.frame id。就像是：

do.call(rbind.data.frame,
        by(dat, dat$id, function(x) list(id=x$id[1], pv=t.test(x, dist, paired=F)$p.value)))
#   id           pv
# 4  4 1.319941e-03
# 6  6 2.877065e-03
# 8  8 6.670216e-05

r - 多个测试

2 回答 2

Related

Reference