2

非常に初歩的な質問で恐縮なのですが、下記コードによる計算結果をtxt等に出力する方法を教えていただけないでしょうか?

下記のように4列、複数行のデータ(csv)があり、これを各列すべての組み合わせでMann-Whitney-testを行います。

       A1 A2 A3 A4
       31 42 33 54
       56 74 48 69
       11 12 13 14
       16 57 18 19
       41 22 63 24
       26 57 28 22
       31 32 33 31
       36 37 78 31
       41 42 43 41

いろいろなサイトを見て回り、下記のようなコードを使えば、解析可能であることがわかりました。

Data <- read.csv("test.csv", header = TRUE, skip = 0)
library(plyr)
combos <- combn(ncol(Data),2)

adply(combos, 2, function(x) {
test <- wilcox.test(Data[, x[1]], Data[, x[2]])

out <- data.frame("var1" = colnames(Data)[x[1]]
              , "var2" = colnames(Data[x[2]])
              ,  "p.value" = sprintf("%.3f", test$p.value)
)
return(out)
})

しかし、データ数(列数)を増やした場合に、大量の結果が出力されてしまい、[ reached getOption("max.print") -- 9397 行を無視しました ] のように出力されます。
最終的には下記にようなテーブルにしたいのですが、ここまで整っていなくとも現時点で出力される結果を、txtなんかにできればと考えています。

        A1   A2   A4   A4
    A1 0.03 0.06 0.05 0.9
    A2 0.99 0.1. 0.56 0.23
    A3 0.32 0.07 0.xx 0.xx
    A4 0.xx 0.xx 0.xx 0.xx

素人考えで下記のように、write,tableを挿入してみましたが、「 0 列 0 行のデータフレーム」と出力され、空のtextが生成されました。

adply(combos, 2, function(x) {
test <- t.test(Data[, x[1]], Data[, x[2]])
out <- data.frame("var1" = colnames(Data)[x[1]]
, "var2" = colnames(Data[x[2]])
,  "p.value" = sprintf("%.100f", test$p.value))
write.table(out, file = "result3.txt", sep = "\t")
})

ご指導のほど、よろしくお願いいたします。

2
  • 実際には adply() の戻り値を TSV(Tab-separated values) 形式でファイルへ保存したいのではないでしょうか? result <- adply(combos, 2, function(x) { ... return(out) }) として、write.table(result, file = "result3.txt", quote = F, row.names = F, sep = "\t") するとか。
    – user39889
    Commented 2017年7月5日 8:54
  • @metropolis ありがとうございます。素人ゆえにadply()を一旦保存するということを想像できませんでした。今後ともご指導のほどよろしくお願いします。
    – Umauma
    Commented 2017年7月6日 1:39

2 件の回答 2

0

こういう感じでしょうか.
adply 内のout には1ペア分の結果しか入ってないので,ファイルに書き出すのはadplyの結果のdata.frameです.

sprintf(".100f", x) としてしまうと,小数点以下が100桁になってしまうので,そのまま数値で書き出せばよいと思います.

Data <- read.csv("test.csv", header=TRUE)
combos <- combn(colnames(Data), 2)
result <- plyr::adply(combos, 2, function(x) {
  test <- wilcox.test(Data[[x[1]]], Data[[x[2]]])
  out <- data.frame(var1 = x[1], var2 = x[2], p.value = test$p.value, stringsAsFactors = FALSE)
  return(out)
}, .id = NULL)
write.table(result, file = "result3.txt", sep = "\t", row.names = FALSE)

また,結果を表形式にするのであれば,例えばこんな感じになります.

n <- ncol(Data)
ret <- matrix(rep(NA, n^2), ncol = n, nrow = n, dimnames = list(colnames(Data), 
  colnames(Data)))
for (i in 1:nrow(result)) {
  ret[result$var1[i], result$var2[i]] <- sprintf("%.3f", result$p.value[i])
}
ret
#>    A1 A2      A3      A4     
#> A1 NA "0.157" "0.426" "0.929"
#> A2 NA NA      "0.930" "0.269"
#> A3 NA NA      NA      "0.536"
#> A4 NA NA      NA      NA
1
  • what_alnkさん ありがとうございます。データを別のものに変えた時、それぞれの郡のn数が違う場合には、wilcox.test()が使えないようでしたのでexactRankTestsパッケーイに搭載されている ` wilcox.exact()` を使用してみました。また表形式への変換方法も求めている通りのものになります。textからexcelへコピーした場合に、左上(A1カラム)にA1がコピーされるようで、headerが一個ずれますが、ここは手動で戻せば問題ありません。(欲を言えば、ここが改善できればと試行錯誤しています。)
    – Umauma
    Commented 2017年7月6日 1:42
0

utils::stack()tidyr::gather()でデータを縦長の形に変形して、標準の pairwise.wilcox.test() を使ってみてはいかがでしょうか。
今回は、多重検定の補正はしていないようですので、

stackedData <- stack(Data)
head(stackedData)
#   values ind
# 1     31  A1
# 2     56  A1
# 3     11  A1
# 4     16  A1
# 5     41  A1
# 6     26  A1

with(stackedData,
     pairwise.wilcox.test(values, ind, exact = FALSE,
                          p.adjust.method = "none"))

#   Pairwise comparisons using Wilcoxon rank sum test 
#
# data:  values and ind 
# 
#    A1   A2   A3  
# A2 0.16 -    -   
# A3 0.43 0.93 -   
# A4 0.93 0.27 0.54
# 
# P value adjustment method: none 

結果を他所に書き出すなら、上記結果の$p.valueだけ取り出せば良いと思います。

temp <- with(stackedData,
             pairwise.wilcox.test(values, ind, exact = FALSE,
                                  p.adjust.method = "none"))$p.value
temp
#           A1        A2        A3
# A2 0.1568459        NA        NA
# A3 0.4260596 0.9295276        NA
# A4 0.9291240 0.2686996 0.5360784

はじめて回答してみましたので、不備などありましたらご指摘ください。

1
  • mokztkさん ありがとうございます。多重検定の補正を行わないという選択肢もあったのですね。思いつきませんでした。上記のwilcox.exact() の結果と比較してみたところ、わずかに(小数点4-5桁くらい)数字が異なる結果があり、現在その関数の違いを調べています。こちらの回答もベストアンサーだと思いますが、先のコメントにチェックさせていただきました。
    – Umauma
    Commented 2017年7月6日 1:51

この質問に回答するには、ログインする必要があります。

求めていた回答ではありませんか? のタグが付いた他の質問を参照する。