---

No.21435　Re: シミュレーションによる検出力　　【青木繁伸】　2014/11/05(Wed) 19:12

DF <- (Time-1)*(N-1) # DFは，上記Meanの自由度
は，関数の中で定義すべきものでしょう？
今の状態では，関数の中で使われている DF は関数の外で定義されたもの，関数の引数で定義した N を反映しないものになっている。

なお，その他のプログラム部分が正しいかどうかは見ていません。

---

No.21436　シミュレーションによる検出力　　【赤岳】　2014/11/05(Wed) 22:10

青木先生，いつもご指導ありがとうございます。
確かにご指摘のところはおかしかったので，修正しましたが，結果は依然期待と違っています。

> Time <- 5 #W0, W2,W4, W8, W12 #観察時期
> Mean <- -0.18 # Mixed effect modelで算出した時期と薬剤の交互作用の平均
> SE <- 0.076 # 上記誤差
> Power <- function(loop, N) # loopはシミュレーションの回数，Nは例数
+ {
+ counta <- 0
+ for(i in 1:loop){
+ DF <- (Time-1)*(N-1) # DFは，上記Meanの自由度
+ x <- rnorm(N, mean= Mean, sd= sqrt(N)*SE)
+ newSE <- sd(x)/sqrt(N)
+ T <- mean(x)/newSE
+ P.value <- 2*pt(abs(T), DF, lower.tail=FALSE)
+ if (P.value < 0.05)
+ counta <- counta+1 }
+ return(counta/loop) }
> Power(10000, 20)
[1] 0.6578
> Power(10000, 50)
[1] 0.6467
> Power(10000, 100)
[1] 0.6546

---

No.21437　Re: シミュレーションによる検出力　　【青木繁伸】　2014/11/05(Wed) 22:56

なんの検定をやっているのか？
DF <- (Time-1)*(N-1)
なのに
x <- rnorm(N, mean= Mean, sd= sqrt(N)*SE)
で N 個のデータしか作っていないのか？
power は 1-βということが分かっているのか？
分かっているなら，
Mean <- -0.18
として
rnorm(N, mean= Mean, sd= sqrt(N)*SE)
でデータを作るのが不可解

母平均が0ではない母集団からデータを取って，母平均=0の検定をやって帰無仮説が棄却される確率が，power でしょ？

> Mean <- 0.5
> SE <- 10
> Power <- function(loop, N) {
+ 	DF <- N-1
+ 	counta <- 0
+ 	for (i in 1:loop) {
+ 		x <- rnorm(N, mean=Mean)
+ 		P.value <- t.test(x, mu=0)$p.value
+ 		if (P.value < 0.05) 
+ 			counta <- counta + 1
+ 	}
+ 	return(counta/loop)
+ }
> Power(10000, 20)
[1] 0.5703
> Power(10000, 50)
[1] 0.9354
> Power(10000, 100)
[1] 0.9989
> Power(10000, 1000)
[1] 1

---

No.21439　Re: シミュレーションによる検出力　　【赤岳】　2014/11/06(Thu) 11:07

青木先生
説明不足ですみません。
試験は，4群（プラセボ，5mg，10mg，20mg）を12週間投与し，効果を0W，2W，4W，8W，12Wの5ポイント測定したものでした。このデータを，次のMixed effect modelで解析しました。
Effect(ij)=β0 + β1*Time(ij) + β2*Treatment(i) + β3*Time(ij)*Treatment(i) + b0(i) + eij
Rのプログラムは，
MIX<- lme(Effect ~ Time + Treatment + Time*Treatment, random=~1|Pt, method="REML", data=Dataset)
この結果，興味ある交互作用項（傾き）は，つぎのようでした。
Fixed effects: Effect ~ Time + Treatment + Time * Treatment
Value Std.Error DF t-value p-value
(Intercept，Timeは省略）
Time:Treatment[T.5mg] -1.96750 1.49467 68 -1.3163460 0.1925
Time:Treatment[T.10mg] -1.43599 1.49467 68 -0.9607422 0.3401
Time:Treatment[T.20mg] -3.97589 1.49467 68 -2.6600547 0.0097

このValueとSEを使って，症例数（DF）を変化させて検出力をシミュレーションしたいと考えました。
rnormがseではなくsdを指定しないといけないこと，DFを変化させたいことから，当初のプルグラムを組みました。
meanとseとDF｛（Time-1)*（N-1）｝を既知として，シミュレーションしたいと考えています。
ちなみに，先生の下記の式で，seとDFをどのように指定すればよいのでしょうか。
x <- rnorm(N, mean=Mean)
+ P.value <- t.test(x, mu=0)$p.value

---

No.21440　Re: シミュレーションによる検出力　　【青木繁伸】　2014/11/06(Thu) 22:42

母平均の検定ですから，あなたのとは全然別です。
あなたの検定の場合は，かなり複雑なので，書かれているプログラムは適切なものではないのでしょう。

---

No.21441　Re: シミュレーションによる検出力　　【赤岳】　2014/11/06(Thu) 23:09

青木先生，
少し考えてみます。
ご指摘ありがとうございました。

---

No.21442　Re: シミュレーションによる検出力　　【赤岳】　2014/11/07(Fri) 21:10

解決したと思います。
2点ミスがありました。
DFを求める式とxの式のsdです。
下記のように修正しましたら，期待した検出力となりました。
お騒がせいたしました。
青木先生，今後ともよろしくご指導お願いいたします。

> Time <- 5 #W0, W2,W4, W8, W12 #観察時期
> Mean <- -0.18 # Mixed effect modelで算出した時期と薬剤の交互作用の平均
> SE <- 0.076 # 上記誤差
> SD <- SE*sqrt(36*5) #36例x5時点
> Treatment<-4 #0，5，10，20mg
> Power <- function(loop, N) # loopはシミュレーションの回数，Nは例数
+ {
+ DF <- N*(Time-1)-Treatment #DFは，上記Meanの自由度
+ counta <- 0
+ for(i in 1:loop){
+ x <- rnorm(N*Time, mean= Mean, sd=SD)
+ newSE <- sd(x)/sqrt(N*Time)
+ T <- mean(x)/newSE
+ P.value <- 2*pt(abs(T), DF, lower.tail=FALSE)
+ if (P.value < 0.05)
+ counta <- counta+1 }
+ return(counta/loop) }
> Power(1000,20)
[1] 0.421
> Power(1000,50)
[1] 0.775
> Power(1000,100)
[1] 0.975

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