「統計学関連なんでもあり」の過去ログ--- 042

No.09071　Ｒで多数の遺伝子のテューキーがしたい　　【tukey】　2009/01/30(Fri) 16:50

Rでテューキーをしたいと思っています。
40個ほどの遺伝子が，3群（コントロール，処理1，処理2）あります。
データのtxtファイル[data.txt]は，最初の3遺伝子のみコピーしますと，以下です。

Master	cont1	cont2	cont3	cont4	treat11	treat12	treat13	treat14	treat21	treat22	treat23	treat24
179	0.034187121	0.056523724	NA	NA	NA	NA	0.048247532	0.021396057	0.043990666	0.013342904	NA	0.039969027
252	0.00877045	0.1069724	0.019199401	0.119519818	0.034588354	0.180613251	0.05095946	0.033463966	0.102124856	0.159025491	0.081311295	NA
350	0.044382947	NA	0.054842779	0.065355601	0.024198232	0.023705043	0.027186462	0.007107383	0.017283822	0.038461196	0.061791639	0.06699642

一つ目　http://aoki2.si.gunma-u.ac.jp/R/tukey.html

のウェブサイトを参考にさせて頂きながら進めていたのですが，うまくいきません。

data <- read.table("data.txt", header=TRUE, row.names=1)

このあとの，コントロール，処理1，処理2のgroup指定がまずできません。

group <- c(rep(1, 4), rep(2, 4), rep(3, 4))

でよいのでしょうか？？

ためしにTukey以下をコピペして，計算させた後に，

> tukey (data, group)

としても，『引数の中に異なった長さのものがあります』と出てデータがでません。
おそらく，グループの設定だおかしいのだと思いますが，これ以上がわかりません。

大変素人質問だと思うのですが，これ以上一人では解決できないと思いましたので，質問させていただきました。
よろしくお願いいたします。

No.09072　Re: Rで多数の遺伝子のテューキーがしたい　　【青木繁伸】　2009/01/30(Fri) 17:42

たとえば，179という遺伝子について，cont, treat1, treat2 で4検体ずつ測定があって，cont, treat1, treat2 について差があるかどうか検定するということですか。

> このあとの，コントロール，処理1，処理2のgroup指定がまずできません。
group <- c(rep(1, 4), rep(2, 4), rep(3, 4))
でよいのでしょうか？？

できないって，どうできないのですか。
> (group <- c(rep(1, 4), rep(2, 4), rep(3, 4)))
 [1] 1 1 1 1 2 2 2 2 3 3 3 3
でいいんじゃないですか？
tukey(data, group)
はだめでしょう。例題をよく見ないと。分からないで出鱈目に試行錯誤してもだめですよ。
data     観察値ベクトル
group    群変数ベクトル
ですよ。あなたがそうやったときの data はデータフレームでしょう。

あなたがやりたいことは，i行目にある遺伝子について，tukey(data[i,], group) のようにするのでは？
しかし，NA があるので，単に tukey(data[1,], group) のようにしてもだめです。
func <- function(i)
{
	d <- cbind(as.numeric(as.matrix(data[i,])), group)
	d <- na.omit(d) # 欠損値を除く（対応する group 要素も）
	tukey(d[,1], d[,2]) # d の1列目がデータ，2列目がグループ
}
のような関数を定義して，for でぐるぐる回して呼び出す
for (i in 1:3) {
	print(func(i))
}
1行目のデータは，func の中でまず，以下のようにしたものが d に代入される
                 group
 [1,] 0.03418712     1
 [2,] 0.05652372     1
 [3,]         NA     1
 [4,]         NA     1
 [5,]         NA     2
 [6,]         NA     2
 [7,] 0.04824753     2
 [8,] 0.02139606     2
 [9,] 0.04399067     3
[10,] 0.01334290     3
[11,]         NA     3
[12,] 0.03996903     3
次に na.omit で NA のある行が取り除かれ，
                group
[1,] 0.03418712     1
[2,] 0.05652372     1
[3,] 0.04824753     2
[4,] 0.02139606     2
[5,] 0.04399067     3
[6,] 0.01334290     3
[7,] 0.03996903     3
のようになる。この1列目がデータ，2列目が群を表すので，
tukey の引数の説明を良く読めば tukey(d[,1], d[,2]) とすれば良いことがわかる。
そして，以下のような結果が表示される。
$result1
       n       Mean     Variance
Group1 2 0.04535542 0.0002494619
Group2 2 0.03482179 0.0003605009
Group3 3 0.03243420 0.0002774016

$result2
            t         p
1:2 0.6172892 0.8190923
1:3 0.8294781 0.7069682
2:3 0.1532717 0.9871737

$phi
[1] 4

$v
[1] 0.0002911915

No.09079　Re: Rで多数の遺伝子のテューキーがしたい　　【Tukey】　2009/01/31(Sat) 11:18

青木先生，お返事が遅くなりまして，すみません。
大変詳細に教えていただき，ありがとうございました。

はじめは，データフレームの意味やベクトルの意味からわからなかったのですが，一つ一つ調べながら理解していったら，無事に全ての遺伝子について，Tukeyすることができました。

functionで欠損値を定義して，データのi行目を順次処理するforの使い方など，一人だとできなかったと思います。

急にRを使えるような気がしてきましたので，もっと勉強してみようと思います。
本当にありがとうございました。