対応のあるデータを再現する方法

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