PCI手術は有効か?という奥村先生の記事であるが
テストデータの生成手順が若干甘いのではないか?
データ生成で,y1 と x1,y2 と x2 は対応のあるデータなので,単にそれぞれの平均値と標準偏差だけでは特定できない。
x1 と y1,x2 と y2 の相関係数を定めないといけない。
奥村先生のやり方では,d を通じて両者の相関係数が決定される。
ただし,set.seed の値を変える(2箇所同じにすること)と,d が変わり,相関係数も変わる。
d ではなく,相関係数そのものを指定することで,データが特定できる。
なお,t.test の p 値が 0.2 という制約があるので,r1 か r2 のどちらかを決めれば,もう一方は自動的に決まる。
そこで,以下のようなプログラムを書いた。
library(MASS) # mvrnorm
g = function(r1 = 0.75) { # r1 は PCL での cor(x1, y1)
f = function(r2) {
# set.seed(180210) # 不要
d1 = mvrnorm(n = 104, mu=c(0, 0), Sigma=matrix(c(1,r1,r1,1),2), empirical=TRUE)
x1 = d1[, 1]
y1 = d1[, 2]
x1 = (x1 - mean(x1)) / sd(x1) * 178.7 + 528.0
y1 = (y1 - mean(y1)) / sd(y1) * 178.7 + 556.4
d2 = mvrnorm(n = 90, mu=c(0, 0), Sigma=matrix(c(1,r2,r2,1),2), empirical=TRUE)
x2 = d2[, 1]
y2 = d2[, 2]
x2 = (x2 - mean(x2)) / sd(x2) * 195.0 + 490.0
y2 = (y2 - mean(y2)) / sd(y2) * 190.9 + 501.8
t.test(y1 - x1, y2 - x2)$p.value - 0.2
}
# グローバル変数として r1 が与えられたとき
# t.test(y1 - x1, y2 - x2)$p.value = 0.2 となる r2 = cor(x2, y2) を求める
r2 = uniroot(f, c(0.1, 0.99))$root
#set.seed(180210) # 不要
# 平均値(標準偏差)が 528.0(178.7) と 556.4(178.7) の 2 変数で,
# 両者の相関係数が r1 となるデータ x1, y1 を生成
d1 = mvrnorm(n = 104, mu=c(0, 0), Sigma=matrix(c(1,r1,r1,1),2), empirical=TRUE)
x1 = d1[, 1]
y1 = d1[, 2]
x1 = x1 * 178.7 + 528.0
y1 = y1 * 178.7 + 556.4
# 平均値,標準偏差,相関係数
cat(sprintf("x1 = %.1f(%.1f) y1 = %.1f(%.1f) r(x1, y1) = %.5f\n",
mean(x1), sd(x1), mean(y1), sd(y1), cor(x1, y1)))
# 平均値(標準偏差)が 490.0(195.0) と 501.8(190.9) の 2 変数で,
# 両者の相関係数が r2 となるデータ x2, y2 を生成
d2 = mvrnorm(n = 90, mu=c(0, 0), Sigma=matrix(c(1,r2,r2,1),2), empirical=TRUE)
x2 = d2[, 1]
y2 = d2[, 2]
x2 = x2 * 195.0 + 490.0
y2 = y2 * 190.9 + 501.8
# 平均値,標準偏差,相関係数
cat(sprintf("x2 = %.1f(%.1f) y2 = %.1f(%.1f) r(x2, y2) = %.5f\n",
mean(x2), sd(x2), mean(y2), sd(y2), cor(x2, y2)))
# 念のため,t 検定の p 値が 0.2 になることを確認
cat(sprintf("t.test() p value = %.5f\n", t.test(y1 - x1, y2 - x2)$p.value))
# 重回帰分析を行い,PCL とプラセボの係数の p 値を見る
x = c(x1, x2)
y = c(y1, y2)
z = rep(1:0, c(104, 90))
r = lm(y ~ x + z)
cat(sprintf(" lm() p value = %.5f\n", summary(r)$coefficients[3, 4]))
}
g(0.75)
g(0.80)
g(0.90)
g(0.95)
g(0.97)
g(0.99)
実行結果は以下の通りであるが,まとめると,
・ x1, y1 の相関が大きい場合には,lm の p 値は小さくなる
・ 場合によっては,0.078 程度にまで小さくなる
> g(0.75)
x1 = 528.0(178.7) y1 = 556.4(178.7) r(x1, y1) = 0.75000
x2 = 490.0(195.0) y2 = 501.8(190.9) r(x2, y2) = 0.98524
t.test() p value = 0.19999
lm() p value = 0.10001
> g(0.80)
x1 = 528.0(178.7) y1 = 556.4(178.7) r(x1, y1) = 0.80000
x2 = 490.0(195.0) y2 = 501.8(190.9) r(x2, y2) = 0.94750
t.test() p value = 0.20000
lm() p value = 0.09625
> g(0.90)
x1 = 528.0(178.7) y1 = 556.4(178.7) r(x1, y1) = 0.90000
x2 = 490.0(195.0) y2 = 501.8(190.9) r(x2, y2) = 0.87319
t.test() p value = 0.20000
lm() p value = 0.08751
> g(0.95)
x1 = 528.0(178.7) y1 = 556.4(178.7) r(x1, y1) = 0.95000
x2 = 490.0(195.0) y2 = 501.8(190.9) r(x2, y2) = 0.83660
t.test() p value = 0.20000
lm() p value = 0.08245
> g(0.97)
x1 = 528.0(178.7) y1 = 556.4(178.7) r(x1, y1) = 0.97000
x2 = 490.0(195.0) y2 = 501.8(190.9) r(x2, y2) = 0.82208
t.test() p value = 0.20000
lm() p value = 0.08029
> g(0.99)
x1 = 528.0(178.7) y1 = 556.4(178.7) r(x1, y1) = 0.99000
x2 = 490.0(195.0) y2 = 501.8(190.9) r(x2, y2) = 0.80761
t.test() p value = 0.20000
lm() p value = 0.07804
※コメント投稿者のブログIDはブログ作成者のみに通知されます