★ Rのソースを教えてください。 ★

9163.　Rのソースを教えてください。　souta　2006/01/23 (月) 17:04
└9164.　Re: Rのソースを教えてください。　青木繁伸　2006/01/23 (月) 17:09
　└9184.　Re^2: Rのソースを教えてください。　souta　2006/01/24 (火) 08:34
　　└9191.　Re^3: Rのソースを教えてください。　青木繁伸　2006/01/24 (火) 13:53
　　　└9212.　Re^4: Rのソースを教えてください。　マスオ　2006/01/25 (水) 00:55
　　　　├9287.　Re^5: Rのソースを教えてください。（パート1）　souta　2006/02/01 (水) 10:32
　　　　│└9288.　Re^6: Rのソースを教えてください。（パート2）　souta　2006/02/01 (水) 10:33
　　　　│　└9304.　Re^7: Rのソースを教えてください。（パート2）　マスオ　2006/02/03 (金) 01:14
　　　　│　　└9306.　Re^8: Rのソースを教えてください。（パート2）　souta　2006/02/03 (金) 09:02
　　　　│　　　└9309.　Re^9: Rのソースを教えてください。（パート2）　souta　2006/02/03 (金) 16:30
　　　　│　　　　└9316.　Re^10: Rのソースを教えてください。（パート2）　マスオ　2006/02/04 (土) 01:09
　　　　│　　　　　└9317.　Re^11: Rのソースを教えてください。（パート2）　souta　2006/02/04 (土) 01:33
　　　　└9221.　Re^5: Rのソースを教えてください。　souta　2006/01/25 (水) 18:39
　　　　　└9223.　Re^6: Rのソースを教えてください。　青木繁伸　2006/01/25 (水) 22:11
　　　　　　└9226.　Re^7: Rのソースを教えてください。　マスオ　2006/01/26 (木) 00:46
　　　　　　　└9227.　Re^8: Rのソースを教えてください。　souta　2006/01/26 (木) 01:17
　　　　　　　　└9240.　Re^9: Rのソースを教えてください。　青木繁伸　2006/01/26 (木) 19:20

9163.　Rのソースを教えてください。　souta　　2006/01/23 (月) 17:04
Rを使って，2×2クロスオーバー法の分散分析をしたいのですが（持ち越し効果を確認し，信頼区間を求めるため），そのソースを誰か教えていただけないでしょうか。n1＝n2のものと，n1≠n2の両方知りたいです。青木先生のhttp://aoki2.si.gunma-u.ac.jp/R/index.htmlにあるソースも見たのですが，どれがあてはまるのか？または存在しないのかよく分かりませんでした・・・。

　　　　　[このページのトップへ]

9164.　Re: Rのソースを教えてください。　青木繁伸　　2006/01/23 (月) 17:09
> 青木先生のhttp://aoki2.si.gunma-u.ac.jp/R/index.htmlにあるソースも見たのですが，どれがあてはまるのか？または存在しないのかよく分かりませんでした・・・。

ないと思います。

その分析法を知らないので，書くことができませんからねぇ。

　　　　　[このページのトップへ]

9184.　Re^2: Rのソースを教えてください。　souta　　2006/01/24 (火) 08:34
> ないと思います。
>
> その分析法を知らないので，書くことができませんからねぇ。

ありがとうございました。
「ある」というのは簡単ですけど，「ない」というのはなかなか判断が難しく困っていました。

ところで，他の人で知ってる人がありましたら是非教えてください。もちろんn1＝n2だけでも大変ありがたいです。

　　　　　[このページのトップへ]

9191.　Re^3: Rのソースを教えてください。　青木繁伸　　2006/01/24 (火) 13:53
みんな SAS とかにたよりきっているようですね。
どっかにアルゴリズム書いてあるところでもあるんですかねぇ。
アルゴリズムもわからない R の書き方もわからないでは，どうしようもない。

　　　　　[このページのトップへ]

9212.　Re^4: Rのソースを教えてください。　マスオ　　2006/01/25 (水) 00:55
例数が揃っていれば，

summary(aov(response ~ treatment * period + Error(subject)))

では？

例数不揃いだと Rではちょっとやっかいですね．分散分析のパッケージが欲しい．

それでみんな SASを使ってるんでしょうか．．．

　　　　　[このページのトップへ]

9287.　Re^5: Rのソースを教えてください。（パート1）　souta　　2006/02/01 (水) 10:32
だいぶん前の質問の続きで申し訳ありませんが，またまた質問させてください。
マスオ様，青木先生，「今頃になって・・・」とお思いのところをどうにか抑えていただいて，ヨロシクお願いします。

以前のマスオ様のコメントに

> 例数が揃っていれば，
> summary(aov(response ~ treatment * period + Error(subject)))
> では？
> 例数不揃いだと Rではちょっとやっかいですね．分散分析のパッケージが欲しい．

と仰っていましたが，例数不揃いのケースについてマスオ様のソースを使ってみるとちゃんとできました。なぜできるのでしょうか（詳細は，パート2参照）？

わたしも，不揃いのケースでは，例数を何らかのかたちで平均化（調和平均するということでしょうか？このあたりも今検討中です。）しなければならず，また計算式も異なる（特に「群又は持込効果 平方和」の箇所）ので，マスオ様のソースではできないのではと考えていました。

「できるのだったら，イイジャン！」とは仰らずに何卒よろしくお願いします。

＝＝＝＝＝（パート2）に続きます。＝＝＝＝＝＝＝＝＝＝＝＝

　　　　　[このページのトップへ]

9288.　Re^6: Rのソースを教えてください。（パート2）　souta　　2006/02/01 (水) 10:33
＝＝＝＝＝（パート1）の続き＝＝＝＝＝＝＝＝＝＝
Rの解析例
（データは，医薬品研究 P130 表13，14 （1984）より抜粋）
> table1<-read.delim("clipboard")
> table1
response treatment period subject
1 9.76 1 1 1
2 10.10 1 1 2
3 11.93 1 1 3
4 6.09 1 1 4
5 13.18 1 1 5
6 7.53 2 1 6
7 8.85 2 1 7
8 12.22 2 1 8
9 9.87 2 1 9
10 11.64 2 1 10
11 8.85 2 1 11
12 7.88 2 2 1
13 7.58 2 2 2
14 11.51 2 2 3
15 8.47 2 2 4
16 8.61 2 2 5
17 9.82 1 2 6
18 10.37 1 2 7
19 10.48 1 2 8
20 9.14 1 2 9
21 10.20 1 2 10
22 10.17 1 2 11
> response<-c(table1$response) #データ
> treatment<-factor(c(table1$treatment)) #処理
> period<-factor(c(table1$period)) #時期
> subject<-factor(c(table1$subject)) #被験者
> summary(aov(response ~ treatment * period + Error(subject)))

Error: subject
Df SumSq MeanSq Fvalue Pr(>F)
treatment:period 1 0.950 0.950 0.2411 0.6352
Residuals 9 35.462 3.940

Error: Within
Df SumSq MeanSq Fvalue Pr(>F)
treatment 1 3.0788 3.0788 1.3365 0.2774
period 1 1.9593 1.9593 0.8505 0.3805
Residuals 9 20.7327 2.3036
>

（医薬品研究 P130 表13，14 （1984）より抜粋）
変動要因　 自由度 平方和　平均平方　分散比
群又は持込効果　1 0.9508 0.9508 0.241　
被験者／群 9 35.4623 3.9403 1.875

時期 1 1.5238
時期を調整した薬 1 3.5141 3.5141 1.527
薬 1 3.0788
薬を調整した時期 1 1.9594 1.9594 0.852
残差 9 20.7328 2.304
総変動 21 62.1831

＝＝＝＝＝＝＝＝おわり＝＝＝＝＝＝＝＝＝＝

　　　　　[このページのトップへ]

9304.　Re^7: Rのソースを教えてください。（パート2）　マスオ　　2006/02/03 (金) 01:14
> ＝＝＝＝＝（パート1）の続き＝＝＝＝＝＝＝＝＝＝
> Rの解析例
> ．．．

これが同じに見えてしまっては困りますね．
目を皿にして見直してごらんなさい．
肝心のところが．．．

　　　　　[このページのトップへ]

9306.　Re^8: Rのソースを教えてください。（パート2）　souta　　2006/02/03 (金) 09:02
お忙しいところコメントありがとうございます。

> これが同じに見えてしまっては困りますね．
> 目を皿にして見直してごらんなさい．

んんんんんーーーーー。
ただいま目を皿にしておりますので，少々お待ちください。

　　　　　[このページのトップへ]

9309.　Re^9: Rのソースを教えてください。（パート2）　souta　　2006/02/03 (金) 16:30
> 肝心のところが．．．

肝心のところとは，R解析例中の薬変動（treatment）の値（3.0788）が，参考資料中の時期を調整した薬変動の値（3.5141）になっていない，要するに「薬に関する変動」の分散比が異なる（Rは1.3365で，資料は1.527）ことでしょうか？

いままで漠然と，「時期変動＋時期を調整した薬変動」か「薬変動＋薬を調整した時期変動」のどちらかが分かればそれでいいと思っていました。

と いうのは，平方和を見ると「時期＋時期を調整した薬＝1.5238＋3.5141＝5.0379」と「薬＋薬を調整した時期＝3.0788＋1.9594 ＝5.0382」で同じ値となり，「被験者間変動（0.9508＋35.4623）＋5.0379＋残差変動（20.7328）=62.1838＝総変動 （62.1831）」となるので，時期と薬のうちどちらか一方が調整されたら，例数不揃いのケースの分散分析としてこと足りると思っていました。

ただ改めて資料をよく見ると，時期（調整なし）及び薬（調整なし）の行の分散比欄は空欄ですよね・・・。

例数不揃いを調整しているのだとは思いますが，「・・・を調整した・・・」とはどういう意味なのでしょうか？
なぜ，「調整した薬（3.5141）」及び「調整した時期（1.9594）」を用いて，残差を「（5.0379＋20.7328）−3.5141−1.9594＝20.2972」とならないのでしょうか？

お忙しいとお察し致しますが何卒ご教示を・・・・。

　　　　　[このページのトップへ]

9316.　Re^10: Rのソースを教えてください。（パート2）　マスオ　　2006/02/04 (土) 01:09
> 肝心のところとは，R解析例中の薬変動（treatment）の値（3.0788）が，参考資料中の時期を調整した薬変動の値（3.5141）になってい ない，要するに「薬に関する変動」の分散比が異なる（Rは1.3365で，資料は1.527）ことでしょうか？
>
> ただ改めて資料をよく見ると，時期（調整なし）及び薬（調整なし）の行の分散比欄は空欄ですよね・・・。

ピンポ〜ン！　正解です．

> 例数不揃いを調整しているのだとは思いますが，「・・・を調整した・・・」とはどういう意味なのでしょうか？
> なぜ，「調整した薬（3.5141）」及び「調整した時期（1.9594）」を用いて，残差を「（5.0379＋20.7328）−3.5141−1.9594＝20.2972」とならないのでしょうか？

これが例数不揃いの時には別の方法があり，(一般）線形モデルがもてはやされた（そして SASが売れた？）理由ですね．
こちらで勉強しましょう．
http://aoki2.si.gunma-u.ac.jp/lecture/TwoWayANOVA/TwoWay4.html

今の場合，Rの作法に従えば，

mod.t.p <- aov(response ~ treatment + period + Error(subject))
mod.t <- aov(response ~ treatment + Error(subject))
mod.p <- aov(response ~ period + Error(subject))
anova(mod.p$Within, mod.t.p$Within)
anova(mod.t$Within, mod.t.p$Within)

ということでしょう．

　　　　　[このページのトップへ]

9317.　Re^11: Rのソースを教えてください。（パート2）　souta　　2006/02/04 (土) 01:33
お忙しいところコメントありがとうございました。

> ピンポ〜ン！　正解です．
ホッ。

> こちらで勉強しましょう．
> http://aoki2.si.gunma-u.ac.jp/lecture/TwoWayANOVA/TwoWay4.html
よく読ませて頂きます（青木先生，拝見させて頂きます）。

> 今の場合，Rの作法に従えば，
　・
　・
> ということでしょう．
ありがたや・・・！
Rで試してみます。

　　　　　[このページのトップへ]

9221.　Re^5: Rのソースを教えてください。　souta　　2006/01/25 (水) 18:39
お忙しいところコメントありがとうございます。

> 例数が揃っていれば，
> summary(aov(response ~ treatment * period + Error(subject)))

私はRを使った経験がないのですが（SASも），青木先生のHP等の資料を参考にして使ってみました。その結果を次に示します。

どうも数値が参考にした文献値とあわないのですが（特に残差のところが・・・）。お忙しいと存じますが，ご指導よろしくお願いします。

文献値（医薬品研究 1106〜1119 （1982））
被験者　1期　2期　　
　A 　　137 　170
　B　　 89 61
C 82 59
D 115 97
E 154 92

被験者　1期　2期
F 94 99
G 48 107
H 110 53
I 110 63
J 146 101
〔A〜E（1期）：対照薬，A〜E（2期）：試験薬〕
〔F〜J（1期）：試験薬，F〜J（2期）：対照薬〕

> response<-c(137,89,82,115,154,94,48,110,110,146,170,61,59,97,92,99,107,53,63,101) #データ
> treatment<-c(1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,1,1,1,1,1) #試験薬と対照薬　
> period<-c(1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2) #時期
> subject<-c(1,1,1,1,1,2,2,2,2,2,1,1,1,1,1,2,2,2,2,2) #群　
> summary(aov(response ~ treatment * period + Error(subject)))

Error: subject
Df Sum Sq Mean Sq
treatment:period 1 781.25 781.25

Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
treatment 1 8.4 8.4 0.0070 0.9344
period 1 1674.4 1674.4 1.3851 0.2564
Residuals 16 19342.4 1208.9

マスオ様には以前も大変お世話になりましたが，今回もヨロシクお願いします。

　　　　　[このページのトップへ]

9223.　Re^6: Rのソースを教えてください。　青木繁伸　　2006/01/25 (水) 22:11
> どうも数値が参考にした文献値とあわないのですが（特に残差のところが・・・）。お忙しいと存じますが，ご指導よろしくお願いします。

その文献にある結果の数値も書いてあげた方が親切なんじゃないでしょうか？

　　　　　[このページのトップへ]

9226.　Re^7: Rのソースを教えてください。　マスオ　　2006/01/26 (木) 00:46
> > どうも数値が参考にした文献値とあわないのですが（特に残差のところが・・・）。お忙しいと存じますが，ご指導よろしくお願いします。
>
> その文献にある結果の数値も書いてあげた方が親切なんじゃないでしょうか？

まあ，残差のところということで見当はつきますが，あった方が話が早いですね．

subjectは被験者です．群の効果？は，処理*時期 交互作用にあたります．
それと，Rでは因子は factor()で名義変数に変換したものを使います．
今はたまたま二値変数しかなかったので，数値のままでもちゃんと計算されているようですが．．．

ということで，以下でいかがでしょうか．

> response<-c(137,89,82,115,154,94,48,110,110,146,170,61,59,97,92,99,107,53,63,101) #データ
> treatment<-factor(c(1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,1,1,1,1,1)) #処理
> period<-factor(c(1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2)) #時期
> subject<-factor(c(1,2,3,4,5,6,7,8,9,10,1,2,3,4,5,6,7,8,9,10)) #被験者
> summary(aov(response ~ treatment * period + Error(subject)))

Error: subject
Df Sum Sq Mean Sq F value Pr(>F)
treatment:period 1 781.3 781.3 0.5104 0.4953
Residuals 8 12245.8 1530.7

Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
treatment 1 8.5 8.5 0.0095 0.9247
period 1 1674.4 1674.4 1.8876 0.2067
Residuals 8 7096.6 887.1

　　　　　[このページのトップへ]

9227.　Re^8: Rのソースを教えてください。　souta　　2006/01/26 (木) 01:17
お忙しいところありがとうございました。
結果は完璧に文献値と一致しています。
本当にありがとうございました。
これで安心して眠れます。

>> その文献にある結果の数値も書いてあげた方が親切なんじゃな>> いでしょうか？
ただ実際のデータを正確に伝えることが一番重要だと考えたのと，幾分コメントが長くなったため字数制限のことを考えて省きました。今考えると確かに不親切だったと思います（コメントを分けたらよかったと思います）。

> まあ，残差のところということで見当はつきますが，あった方> が話が早いですね．
申し訳ありませんでした。

> subjectは被験者です．
英語が苦手なのがばれたようでお恥ずかしい。

> それと，Rでは因子は factor()で名義変数に変換したものを使> います．
ご指導ありがとうございます。

> subject<-factor(c (1,2,3,4,5,6,7,8,9,10,1,2,3,4,5,6,7,8,9,10)) #被験者
こうとらえるのですか。なるほどーーー。

　　　　　[このページのトップへ]

9240.　Re^9: Rのソースを教えてください。　青木繁伸　　2006/01/26 (木) 19:20
あちらの方へも，n1=n2 のときの解答があったということは報告しておくのがよいでしょうね。

　　　　　[このページのトップへ]

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