1

我想对尺寸（50x752）的数据框进行 chisq.test 。我想为所有列的所有可能的成对比较获取 pvalues（通过多次测试调整）。最后，我想取回一个矩阵 (50x50) 以生成调整后的 chisq pvalues 的热图。这是我目前所做的，但这远非理想。

Step1：进行成对比较

function(data,p.adjust.method="holm")
{
cor.mat <- cor(data)
x<-ncol(data)#nb of column in matrix here 50
y<-nrow(data)#nb of column in matrix here 758
index<-t(combn(x, 2)) #create the matrix position of output for all possible combination
nindex <- nrow(index)
pvals <- numeric(nindex)

for (i in 1:nindex)
{
pvals[i]<-chisq.test(data[, index[i, 1]], data[, index[i,2]])$p.value   
}
pvals<-p.adjust(pvals,method = p.adjust.method)
out <- as.data.frame(cbind(index, pvals))
}

Step2：将输出表转换为矩阵

   dcast(df,V2~V1,fill=1) # thanx to Roland for this function!

但这效果不佳，因为我没有在最终矩阵中镜像 pvalue，并且我必须操纵第一个函数的输出以使对角线填充为 0（将列与自身进行比较时）。对你的帮助表示感谢！

4

1 回答 1

2

像这样？

#some data
set.seed(42)
df <- data.frame(a=rbinom(1000,5,0.3),
                 b=rbinom(1000,5,0.001),
                 c=rbinom(1000,5,0.1),
                 d=rbinom(1000,5,0.9))

#function to calculate the adj. p-value
fun <- function(x,y) {
  p.adjust(chisq.test(df[,x],df[,y])$p.value,method="holm",n=choose(ncol(df),2))
}

p.adj <- outer(names(df),names(df),FUN=Vectorize(fun)) #use outer to get a matrix
diag(p.adj) <- 1  #you should find out why chisq.test returns zero at the diag
rownames(p.adj) <- names(df)
colnames(p.adj) <- names(df)

p.adj
#  a         b c         d
#a 1 1.0000000 1 1.0000000
#b 1 1.0000000 1 0.6152165
#c 1 1.0000000 1 1.0000000
#d 1 0.6152165 1 1.0000000
于 2012-10-26T16:41:22.513 回答