No.12750 Rでの分散分析について  【ぽち】 2010/06/01(Tue) 09:54

よろしくお願いいたします。

SPSSやExcelからRへ乗り換えるために,青木先生の御著者などを参考に必死で勉強をしております。

先日も質問させて頂いたのですが,Rでの対応ありの一要因分散分析および二要因分散分析についてなのですが,青木先生や波音さんからご紹介頂きましたページを参考にコピペではありますが,自身で動かしてみました。

一要因分散分析であれば,分散分析表までは計算してくれるのですが,多重比較は計算されません。

以下のようなものです。
RB <- data.frame(
+ a = factor(c(rep(1,8), rep(2,8), rep(3,8), rep(4,8))),
+ s = factor(rep(1:8, 4)),
+ result = c(9,7,8,8,12,11,8,13, 6,5,6,3,6,7,10,9,
+ 10,13,8,13,12,14,14,16, 9,11,13,14,16,12,15,14))
> print(summary(aov(result ~ a + s + Error(a:s), RB)))

Error: a:s
Df Sum Sq Mean Sq F value Pr(>F)
a 3 217.5 72.500 21.4437 1.347e-06 ***
s 7 77.0 11.000 3.2535 0.01682 *
Residuals 21 71.0 3.381
また二要因分散分析でも同様に分散分析表だけです。
CRFpq <- data.frame(
+ a = factor(c(rep(1,20), rep(2,20))),
+ b = factor(rep(c(rep(1,5), rep(2,5), rep(3,5), rep(4,5)),2)),
+ result = c(3,3,1,3,5, 4,3,4,5,7, 6,6,6,4,8, 5,7,8,7,9,
+ 3,5,2,4,6, 2,6,3,6,4, 3,2,3,6,5, 2,3,3,4,6))
> print(summary(aov(result ~ a * b, CRFpq)))
Df Sum Sq Mean Sq F value Pr(>F)
a 1 16.9 16.9000 7.0417 0.01230 *
b 3 19.7 6.5667 2.7361 0.05972 .
a:b 3 30.5 10.1667 4.2361 0.01250 *
Residuals 32 76.8 2.4000
大変申しわけございませんが,対応あり一要因分散分析の多重比較,対応なしの二要因の主効果,単純主効果,多重比較についてのコマンドをご教授願えませんでしょうか?もしくはそれが示されているサイト情報などありましたら,よろしくお願いいたします。

No.12751 Re: Rでの分散分析について  【波音】 2010/06/01(Tue) 10:54

最も手軽で手っ取り早い方法として,pairwise.t.test()という関数が使えます。例えば,

ans <- lm(Y ~ A) # Aはカテゴリカル型の変数
summary(aov(ans))


などとして得られる分散分析について,各群すべての対比較を行うために,

pairwise.t.test(Y, A, p.adj="bonferroni") # ボンフェローニの補正

とするか,あるいはデフォルトではホルム(Holm)の方法によって補正されたp値が得られます。

pairwise.t.test(Y, A) # デフォルトではホルムの補正

---

例えばCRFpqならば,以下のようになります(最もaについては多重比較をする必要がありませんが)。
> attach(CRFpq)
> pairwise.t.test(result, a)

Pairwise comparisons using t tests with pooled SD

data: result and a

1
2 0.030

P value adjustment method: holm
> pairwise.t.test(result, b)

Pairwise comparisons using t tests with pooled SD

data: result and b

1 2 3
2 0.95 - -
3 0.50 1.00 -
4 0.17 0.95 1.00

P value adjustment method: holm

No.12759 Re: Rでの分散分析について  【ぽち】 2010/06/01(Tue) 20:35

波音さん

お返事が遅くなりました。申し訳ありません。

非常に言いにくいことなの ですが,ご教授いただきました対応ある一要因分散分析で,私が参照してコピペを使って練習しました上記コード(?)なのですが,エクセルに入力されたデー タをコピペして使用するには(かつご教授いただきました多重比較を行う場合)どのように書き換えたらよろしいのでしょうか?いろいろと調べてはいるのです が,うまくいかないのです.....。

対応あるデータですから,ExcelにはA列,B列,C列にそれぞれ試験前,試験直後,1週間後と いうタイトル行を設けてそれぞれの数値データが入っています。例えば,このデータをクリップボードにコピーします。そしてxというデータファイル名でコン ソールで読み込む場合,

> x <-read.delim("clipboard")
> x

でデータが読み込めますよね。

その後,ご教授いただきました式を用いて以下のようにしています。

ans <- lm(Y ~ A) # Aはカテゴリカル型の変数
summary(aov(ans))
pairwise.t.test(Y, A, p.adj="bonferroni") # ボンフェローニの補正

ただもちろんこれではダメでして,どのように変形すればよろしいのか,ご教授いただけませんでしょうか?

また多重比較においては95%信頼区間は計算されるでしょうか?

本当に申し訳ございませんが,よろしくお願い致します。

No.12760 Re: Rでの分散分析について  【青木繁伸】 2010/06/01(Tue) 21:32

> d
試験前 試験直後 一週間後
1 13 22 23
2 18 16 21
3 14 19 21
4 19 18 24
5 20 13 21
6 10 21 21
7 15 14 19
8 19 12 20
9 16 15 17
10 15 22 15
の ようなデータフレームとして表現されているデータを pairwise.t.test(…, paired=TRUE) で使うためには,stack 関数を使って,以下の形式,すなわち,測定値とその測定値がどの時期のものであるかを示す変数からなるデータフレームにします。
> (d2 <- stack(d))
values ind
1 13 試験前
2 18 試験前
3 14 試験前
4 19 試験前
5 20 試験前
6 10 試験前
: 途中省略
27 19 一週間後
28 20 一週間後
29 17 一週間後
30 15 一週間後

> pairwise.t.test(d2$values, d2$ind, paired=TRUE) # paired=TRUE の指定を忘れずに

Pairwise comparisons using paired t tests

data: d2$values and d2$ind

一週間後 試験前
試験前 0.021 -
試験直後 0.128 0.534

P value adjustment method: holm
のようになるわけです。

No.12761 Re: Rでの分散分析について  【ぽち】 2010/06/01(Tue) 22:16

青木先生

お返事ありがとうございます。いろいろと質問をしてしまいまして,恐縮なのですが,ご教授いただきました多重比較(holmの方法)では,95%信頼区間を求めることはできないのでしょうか?

よろしくお願い致します。

No.12762 Re: Rでの分散分析について  【青木繁伸】 2010/06/01(Tue) 22:32

> 95%信頼区間を求めることはできないのでしょうか

内部的には検定統計量が得られているのですから,自分で計算することは可能でしょう。pooled.sd のときは,
            dif <- xbar[i] - xbar[j]
se.dif <- pooled.sd * sqrt(1/n[i] + 1/n[j])
t.val <- dif/se.dif
if (alternative == "two.sided")
2 * pt(-abs(t.val), total.degf)
else pt(t.val, total.degf, lower.tail = (alternative ==
"less"))
となっているので,例えば,両側検定に対応する信頼区間は dif+c(-1,1)*qt(0.025, total.degf)*se.dif とでもすればよいのかな。他の部分も追加したり直したりしないといけないけど。
paired=TRUE と non pooled-sd のときは,
            t.test(xi, xj, paired = paired, alternative = alternative, 
...)$p.value
としているので,t.test の戻り値から p.value だけでなく,conf.int も取り出せばよい。

まあ,そういうことするよりは,新たにプログラム書いた方がはやいかな。

No.12764 Re: Rでの分散分析について  【ぽち】 2010/06/02(Wed) 11:02

青木先生

お返事が遅くなりましたこと,申しわけございません。

丁寧な解説に感謝いたします。

正直申しまして,先生のようなコードを書ける自信は全くございません.....。

私のようなR素人,統計ソフトエンドユーザーの人間が,こういったコードをコピペして使えるように準備してくれているサイトや著書をご存じでありませんか?

青木先生のご著書は,すでに購入済みでして,最初から拝読させて頂いております。しかし,R初心者ですので,具体的な統計のところに入ってくると,理解するのに相当な時間を要してしまいます。

このままですと,毎回のように,私が書いたコードをここでチェックして頂くような質問になってしまい,きっとご迷惑をおかけするだろうと思います。

私 はExcelでデータを入力することになれていまして,そのデータをそのままコピペして使いたいのです。おそらく様々な関数を覚えることは出来ないと思い まして,例えば,二元配置分散分析について,ExcelデータをコピペしてRコンソールに読み込み,コードを入力すれば,分散分析表,単純主効果,95% 信頼区間を一気に計算してくれるようなものを期待しているのですが。

もちろん私自身がもっと勉強をすれば良いだけの話しなのかもしれませんが.....。

長々と申しわけございません。

● 「統計学関連なんでもあり」の過去ログ--- 043 の目次へジャンプ
● 「統計学関連なんでもあり」の目次へジャンプ
● 直前のページへ戻る