パッケージのインストール Rには 複雑な解析を便利に行うためのパッケージが容易されています ( 世界中の研究者達が提供してくれる ) 今回は例として多重比較検定用のmultcomp パッケージをインストールしてみます ( 注意 ) 滋賀県立大学のようにプロキシ経由でインターネットに接続する環境で R

Transcription

1 ソフトウェア R を用いた統計解析 清水顕史 R のインストール R の情報 ( 日本語 ) は RjpWikihttp:// にまとめられています 説明に従って最新版の exe ファイルをダウンロード ( し クリックしてインストールします インストール終了後 デスクトップのアイコンをクリックすると R が起動します 起動すると下のようなウィンドウが開きます 左下の端の > 印の右側は カーソル 1 が点 滅しており この位置からキーボード入力が可能です R の操作は コマンドを入力して行います 例えば以下のコマンドを入力し 各行ごとに Enter キーを押して実行してみましょう (# より右はコメント用で入力の必要はありません ) X <- c(5,4,5,4,2,6) mean(x) var(x) sd(x) plot(x) Y <- c(3,4,2,3,4,2) plot(x,y) # ベクトルXに数列 5,4,5,4,2,6 を入力する # ベクトルXの全要素の平均値を計算する # 分散を計算する # 標準偏差を計算する # プロットする # 散布図を描く 以下 R への入力コマンドはゴシック赤字で表記することにします (# 以降はコメント行 のため 入力する必要なし ) 演習上の操作は 入力コマンド部分のみをコピー & ペースト することで再現できます 1 縦棒 があり

2 パッケージのインストール Rには 複雑な解析を便利に行うためのパッケージが容易されています ( 世界中の研究者達が提供してくれる ) 今回は例として多重比較検定用のmultcomp パッケージをインストールしてみます ( 注意 ) 滋賀県立大学のようにプロキシ経由でインターネットに接続する環境で R のパッケージインストールを行う場合 以下のようにアイコンのプロパティ ( アイコンを右クリックして表示 ) で リンク先の Rgui.exe の後ろに 半角空白 --internet2 と入力しておき アイコンをクリックして R を起動し直します ファイルメニューから [ パッケージ ]->[ パッケージのインストール ] を選び CRAN mirror のリストから適当なミラーサイトを選びます ( 例えばJapan(Hyogo)) インストール可能なパッケージリストがアルファベット順で表示されるので multcomp を選びインストールします インストールしたパッケージは R 起動後に library(multcomp) というコマンドを実行すると利用できるようになりま す

3 一要因データの解析 1( 総当たりの要因間の差のチューキー検定 ) 例として 或る栽培条件下で育てたイネ 4 品種 ( 系統番号 WRC1, 23, 57, 60) の或る酵素の活性量を調べたデータ "vartest.txt" を利用します データはエクセルなどを利用して作成されたタブ区切りのテキストファイルで 列 1 は品種名 列 2 が活性量で 一行目は見出しです 各品種 4 個体反復のデータが記録されています 3 種類以上の群の平均値の大小を比較する場合 チューキー検定を用います R を用いた検定手順は以下の通りです まず R のファイルメニューで [ ファイル ] [ ディレクトリの変更 ] で 入力用データの置いてあるフォルダを指定します ( データがデスクトップにある場合はデスクトップを指定 ) その後 以下のコマンドを入力します d <- read.table("vartest1.txt", header=t) Pd <- data.frame(variety=factor(d$val),activity=d$data) plot(activity~variety,data=pd) # データの読み込み # フレームワーク作成 # 箱ヒゲ図の描出

4 箱ヒゲ図 ( 左 ) とその見方 ( 右 ) 各品種の中央値を横太線で 箱の上端はデータを降順にならべた 25% の点 下端は 75% 点を示す 髭はデータの最大値または最小値を示すが ただし上位四分位数 +1.5 四分位範囲および下位四分位数 -1.5 四分位数の範囲を越えたデータは外れ値 ( ) として示す 例データの場合は各品種のデータ数は 4 個なので 降順にならべたデータの一番と二番の平均が上位四分位数 二番と三番の平均が中央値 三番と四番の平均が下位四分位数となる res1 <- aov(activity~variety,data=pd) library(multcomp) # 一元分散分析 # パッケージをロード res2 <- glht(res1, linfct = mcp(variety = "Tukey")) # チューキー クレーマー検定 # 検定結果の表示 検定結果は 全 6 組合せ毎に系統平均間に差がないと仮定した場合の有意確率として表示されます 一般に 有意確率は 0.05 以下の場合に統計学的に有意差があるとみなします ( 考え方については 統計学の入門書などを参考にしましょう ) WRC1 とそれ以外の品種に危険率 5% 以下で有意差がみられました (WRC1 だけ酵素活性が高い )

5 一要因データの解析 2( ダネット検定による対照に対する多重比較 ) CSSLs( 染色体断片置換系統 ) は染色体領域の一部が供与親に置換した系統で それ以外の殆どのゲノムは反復親と同じ情報を持ちます CSSLs は置換領域に含まれる遺伝情報が形質にもたらす影響を調べるのに非常に便利な遺伝解析材料で 形質に差があるかどうかは反復親と各置換系統との差を比較します この場合 総当たりの比較をするチューキー検定ではなく 対照群とその他の各群を比較するダネット検定を用います 例として コシヒカリ ( 反復親 ) とその 39 種類の CSSL(SL201~SL239) を用いて 或る酵素活性量を調べたデータ "CSSLtest.txt" を解析してみます データはエクセルなどを利用して作成されたタブ区切りのテキストファイルで 列 1 は系統名 列 2 が活性量で 一行目は見出しです 各品種 8 個体反復のデータが記録されています multcomp パッケージを利用したダネット検定では 群の名前をアルファベット順に並べたときに一番小さいものが自動的に対照群に割り当てられます 例では Koshi が対照群になり その他の SL 群と比較されます ) 検定は以下のコードで行います s <- read.table("cssltest.txt ", header=t) Sd <- data.frame(line=factor(s$val),activity=s$data) plot(activity~line,data=sd) res1 <- aov(activity~line,data=sd) library(multcomp) res2 <- glht(res1, linfct = mcp(line = "Dunnett")) # ファイルの読み込み # フレームワーク作成 # 箱ヒゲ図 # 一元分散分析 # パッケージ起動 # ダネット検定 # 結果を出力

6 Simultaneous Tests for General Linear Hypotheses Multiple Comparisons of Means: Dunnett Contrasts Fit: aov(formula = ACTIVITY ~ LINE, data = Sd) Linear Hypotheses: Estimate Std. Error t value Pr(> t ) SL201 - Koshihikari == SL202 - Koshihikari == SL203 - Koshihikari == ** SL204 - Koshihikari == SL205 - Koshihikari == SL206 - Koshihikari == SL207 - Koshihikari == SL208 - Koshihikari == SL209 - Koshihikari == SL210 - Koshihikari == SL211 - Koshihikari == SL212 - Koshihikari == SL213 - Koshihikari == SL214 - Koshihikari == SL215 - Koshihikari == SL216 - Koshihikari == SL217 - Koshihikari == SL218 - Koshihikari == SL219 - Koshihikari == SL220 - Koshihikari == SL221 - Koshihikari == SL222 - Koshihikari == SL223 - Koshihikari == SL224 - Koshihikari == < *** SL225 - Koshihikari == SL226 - Koshihikari == SL227 - Koshihikari == SL228 - Koshihikari == SL229 - Koshihikari == SL230 - Koshihikari == SL231 - Koshihikari == SL232 - Koshihikari == SL233 - Koshihikari == SL234 - Koshihikari == SL235 - Koshihikari == SL236 - Koshihikari == SL237 - Koshihikari == SL238 - Koshihikari == SL239 - Koshihikari == < *** --- Signif. codes: 0 *** ** 0.01 * (Adjusted p values reported -- single-step method) 検定の結果より SL203 は危険率 1% で SL224 および 239 は危険率 0.1% でコシヒカリ との有意差が検出された つまりそれらの系統の染色体置換領域に形質に影響を及ぼす遺 伝子が存在する可能性が示唆された なお コシヒカリよりも平均値が高い置換系統のみを片側検定する場合 res2 <- glht(res1, linfct = mcp(line = "Dunnett"), alternative = "greater") コシヒカリよりも平均値が低い置換系統のみを片側検定する場合 res2 <- glht(res1, linfct = mcp(line = "Dunnett"), alternative = "less") とします