RCTで多変量解析#2
http://r-statistics-fan.hatenablog.com/entry/2015/02/19/000248
前回の続き
RCTで多変量解析をするとPower検出力が上がる。
しかし、弱い予測因子を投入するとPowerは下がるはずだというのを
前回実証しようとして失敗した。
今回は、投入する因子を増やして、シミュレーション回数も増やしてみる。
結果、かえってPowerが下がる現象を証明し得た。
#回数増やしたので時間がかかる library(tcltk) ##各変数の相関は今回は考えない set.seed(1) nsim <- 30000 #sim回数 nn <- 200 #各studyのn b0 <- 0 #x0の真の係数 差がない場合 b0d <- 1 ##x0の真の係数 差がある場合 b1 <- 5 #大きな影響因子 b2b <- b2a <- b2 <- 0.01 #小さな影響因子 b3 <- 3 #未測定因子1 b4 <- 2 #未測定因子2 glm_p0_b2 <- glm_p1_b2 <- glm_p1_b1 <- glm_p0_b1 <- p0 <- p1 <- numeric(nsim) pb <- txtProgressBar(min = 1, max = nsim, style = 3) for (i in seq_len(nsim)){ setTxtProgressBar(pb, i) x0 <- sample.int(2, nn, replace = TRUE) - 1 x1 <- rnorm(nn) x2 <- rnorm(nn) x2a <- rnorm(nn) x2b <- rnorm(nn) x3 <- rnorm(nn) x4 <- rnorm(nn) temp <- b1 * x1 + b2 * x2 +b2a * x2a + b2b * x2b + b3 * x3 + b4 + x4 + 0.1 * rnorm(1) y1 <- b0d * x0 + temp y0 <- b0 * x0 + temp p0[i] <- t.test(y0 ~ x0, alternative = "two.sided", var.equal = FALSE)$p.value p1[i] <- t.test(y1 ~ x0, alternative = "two.sided", var.equal = FALSE)$p.value glm_p0_b1[i] <- summary(glm(y0 ~ x0 + x1))$coefficients[2,4] glm_p1_b1[i] <- summary(glm(y1 ~ x0 + x1))$coefficients[2,4] glm_p0_b2[i] <- summary(glm(y0 ~ x0 + x2 + x2a + x2b))$coefficients[2,4] glm_p1_b2[i] <- summary(glm(y1 ~ x0 + x2 + x2a + x2b))$coefficients[2,4] } hist(p0) hist(p1) hist(glm_p0_b1) hist(glm_p1_b1) hist(glm_p0_b2) hist(glm_p1_b2) plot1 <- function(xx){ plot(sort(xx), ylim = c(0,1), pch=20, type = "p") abline(a=0, b=1/nsim, col = "red") } plot1(p0) plot1(p1) plot1(glm_p0_b1) plot1(glm_p1_b1) plot1(glm_p0_b2) plot1(glm_p1_b2) sum(p0 < 0.05)/nsim sum(p1 < 0.05)/nsim sum(glm_p0_b1 < 0.05)/nsim sum(glm_p1_b1 < 0.05)/nsim sum(glm_p0_b2 < 0.05)/nsim sum(glm_p1_b2 < 0.05)/nsim pow2 <- (glm_p1_b2 < 0.05) pow1 <- (p1 < 0.05) table(pow1, pow2) mcnemar.test(pow1, pow2, correct = FALSE) binom.test(367, 819)
> binom.test(367, 819)
Exact binomial test
data: 367 and 819
number of successes = 367, number of trials = 819, p-value = 0.003309
alternative hypothesis: true probability of success is not equal to 0.5
ということで、RCTで弱い予測因子３つを入れて多変量解析をすると
かえって有意に検出力が下がった。
> table(pow1, pow2)
pow2
pow1 FALSE TRUE
FALSE 23118 367
TRUE 452 6063
> binom.test(367, 819)
Exact binomial test
data: 367 and 819
number of successes = 367, number of trials = 819, p-value = 0.003309
alternative hypothesis: true probability of success is not equal to 0.5
> mcnemar.test(pow1, pow2, correct = FALSE)
McNemar's Chi-squared test
data: pow1 and pow2
McNemar's chi-squared = 8.8217, df = 1, p-value = 0.002977
と、いっても、かなり微妙。３万回でたった85回の差でしかない。
RCTで多変量解析
某社のwebカンファを見てきた。
RCTで多変量解析するのは、by-chanceで出た差を調整して
確認するためなどという感じのことを言っていた。
しかし、自分の認識は異なる。
RCTで多変量解析するのは、原則的にあらかじめpowerの増加を
期待して計画された時だけと認識している。
差がでたら困る変数は、層別ランダム化しておけって話。
それでもなお生じる差を更に調整する感じだ。
ここで、多変量解析に投入される因子はアウトカムに多大に影響する
既知の因子である。
しかし、あまり関係ない因子を投入すると、かえってpowerは
下がるのではないか？
ここのところを確認したかったので、シミュレーションを行う
#RCTで多変量解析 library(mvtnorm) library(tcltk) ##各変数の相関は今回は考えない set.seed(1) nsim <- 3000 #sim回数 nn <- 100 #各studyのn b0 <- 0 #x0の真の係数 差がない場合 b0d <- 1 ##x0の真の係数 差がある場合 b1 <- 5 #大きな影響因子 b2 <- 0.02 #小さな影響因子 b3 <- 3 #未測定因子1 b4 <- 2 #未測定因子2 glm_p0_b2 <- glm_p1_b2 <- glm_p1_b1 <- glm_p0_b1 <- p0 <- p1 <- numeric(nsim) pb <- txtProgressBar(min = 1, max = nsim, style = 3) for (i in seq_len(nsim)){ setTxtProgressBar(pb, i) x0 <- sample.int(2, nn, replace = TRUE) - 1 x1 <- rnorm(nn) x2 <- rnorm(nn) x3 <- rnorm(nn) x4 <- rnorm(nn) temp <- b1 * x1 + b2 * x2 + b3 * x3 + b4 + x4 + 0.1 * rnorm(1) y1 <- b0d * x0 + temp y0 <- b0 * x0 + temp p0[i] <- t.test(y0 ~ x0, alternative = "two.sided", var.equal = FALSE)$p.value p1[i] <- t.test(y1 ~ x0, alternative = "two.sided", var.equal = FALSE)$p.value glm_p0_b1[i] <- summary(glm(y0 ~ x0 + x1))$coefficients[2,4] glm_p1_b1[i] <- summary(glm(y1 ~ x0 + x1))$coefficients[2,4] glm_p0_b2[i] <- summary(glm(y0 ~ x0 + x2))$coefficients[2,4] glm_p1_b2[i] <- summary(glm(y1 ~ x0 + x2))$coefficients[2,4] } hist(p0) hist(p1) hist(glm_p0_b1) hist(glm_p1_b1) hist(glm_p0_b2) hist(glm_p1_b2) plot1 <- function(xx){ plot(sort(xx), ylim = c(0,1), pch=20, type = "p") abline(a=0, b=1/nsim, col = "red") } plot1(p0) plot1(p1) plot1(glm_p0_b1) plot1(glm_p1_b1) plot1(glm_p0_b2) plot1(glm_p1_b2) sum(p0 < 0.05)/nsim sum(p1 < 0.05)/nsim sum(glm_p0_b1 < 0.05)/nsim sum(glm_p1_b1 < 0.05)/nsim sum(glm_p0_b2 < 0.05)/nsim sum(glm_p1_b2 < 0.05)/nsim <|| >|| > sum(p0 < 0.05)/nsim [1] 0.05233333 ##RCTのt検定 アルファエラーは0.05に保たれている > sum(p1 < 0.05)/nsim [1] 0.1336667 ##RCTのt検定 Powerは0.13しかない > sum(glm_p0_b1 < 0.05)/nsim [1] 0.05033333 ##RCTの多変量解析 アルファエラーは0.05に保たれている > sum(glm_p1_b1 < 0.05)/nsim [1] 0.342 ##RCTの多変量解析 強い予測因子で調整するとPowerは0.34にUP! > sum(glm_p0_b2 < 0.05)/nsim [1] 0.05233333 ##RCTの多変量解析 アルファエラーは0.05に保たれている > sum(glm_p1_b2 < 0.05)/nsim [1] 0.1333333 ##RCTの多変量解析 弱い予測因子で調整するとPowerは0.13で変わらない。 予測ではかえってPower下がると思ったけど下がらなかったわ。 || ということで、今回のRCTのシミュレーションでは モデルが正しくないとしても、既知の強力な予測因子で多変量解析で 調整すると、単純なRCTよりも、アルファエラーを損なわずにpowerを 高めることが出来る。 弱い予測因子を入れるとかえってPowerが悪化すると思ったけど、今回は 下がらなかった。１個だけじゃ、あまり影響がないのか。多分、沢山弱い因子を 入れたらPowerが下がるんじゃないかな。今日はもう眠いのでやらない。
頭脳王：三山崩しの「逆型」をRで解く。更に高速化一般化
頭脳王：三山崩しの「逆型」
学生時代の旧友よりメールが来た。
数学に関しては私にとっては神様のような方で、学生時代には
よく自作の数学問題を披露してくれたりして、楽しんだものだ。懐かしい。
旧友によると、三山崩しに関しては非「逆型」の必勝法は数理ゲームの理論では
基本に類するそうです。「逆型」でも，三山崩しに関しては，大差はないそうです。
更に面白い問題も教えていただきました
##ゲーム(*)
・黒白2色の碁石がそれぞれいくつかある状態からスタート．
・交互に石を取る．ただし，
「黒だけを任意個数」「白だけを任意個数」「黒白を任意の同数」
のいずれかの取り方だけが可能．
最後の石を取った方が勝ち．(逆型では負け．)
この必勝法はいかに。
##
そして、今回拙ブログを見てR言語に初挑戦したそうです。
素朴に三山崩しの「逆型」を解析するコードを書いてみたそうです。（下記f1()）
R言語ははじめて触るのに、なんかもうすでにこの時点で100倍以上高速。
清々しいほどにぶっちぎりです。
expand.gridの並びの性質を最大限利用して余計な論理判断を省いている
すばらしいコードです。正直、はじめは何でこれで正確な解になるのかが
分かりませんでした。
expand.gridの並びは行ナンバーAの局面を考えると、駒を取れば、
必ずX＜Aとなる局面Xに移行するから、これでも良いわけですね。
しびれる～。本当に楽しい。
そして、evalを使って、expand.gridの行を一般化する手法を例示した後の
コードがf2()です。
更に余計な論理判断を省いてコードが短くなっています。
本当に短い！更に速い！
なお、このようなコーディングは、業務用では落第点らしいです。
それにしてもC言語の素養はあるとはいえ、ちょっとかじっただけで
これほどのコードを作ってしまうとは。いやー、楽しい。
以下コード
#f1 rをはじめて触った友人の初回コード勝手にfunction化したので一部改変 f1 <- function(m=3, n=10){## m = 色の数 # n = 各色の駒数の上限 x <- as.matrix(expand.grid(0:n,0:n,0:n,0)) ## 考察対象の列挙 ## 最後の成分は，0: 未考察 1:必勝形 2:必敗形 x[1,m+1] <- 2 ## (0,0,0) にした人の負け ## 1手で変形できる必勝形があれば必敗，なければ必勝． for (i in 1:nrow(x)) { y <- x[i,]; if (y[m+1]==0) {## 1つの局面を取り出し，まだ未考察なら， x[i,m+1] <- 1+any(x[,m+1]==1 & ## 必勝形で，なおかつ， ## 取り出した局面から1手で移行できるものを探す． ## あれば必敗，なければ必勝 ((x[,1]==y[1] & x[,2]==y[2] & x[,3]<y[3]) | (x[,1]==y[1] & x[,3]==y[3] & x[,2]<y[2]) | (x[,2]==y[2] & x[,3]==y[3] & x[,1]<y[1])))}} ## この記述だと，山が増えると非常にめんどうになりますが， ## 言語をまだよく知らないので，とりあえずこれで．．． return(subset(x,subset=x[,m+1]==1)[,1:m]) ## 得られた必勝形を表示 } #f2 rをはじめて触った友人のコードVer2 勝手にfunction化したので一部改変 f2 <- function(M=3, N=10){ ## m = 色の数 # n = 各色の駒数の上限 x <- as.matrix(eval(parse(text = paste("expand.grid(", paste(rep("0:N", M), sep="", collapse=","), ",0)", sep="")))) x[1,M+1] <- 2 for (i in 1:nrow(x)) {y <- x[i,] if (y[M+1]==0) {x[i,M+1] <- y[M+1] <- 1 x[i,M+1] <- 1+any(rowSums(!!(t(t(x)-y)))==1)}} return(subset(x,subset=x[,M+1]==1)[,1:M]) ## 得られた必勝形を表示 } #自分のコード f3 <- function(n=10){ # n = 駒の最大数 dat <- expand.grid(0:n,0:n,0:n) dat <- dat[-1,] nr <- nrow(dat) row.names(dat) <- seq_len(nr) flag2 <- flag <- numeric(nr) flag3 <- 0 dat <- as.matrix(dat) flag[apply(dat, 1, function(x) sum(x)==1)] <- 1 #flag1=win、1,0,0を相手に渡すと勝利確定 flag=2は負け確定 flag=0は未定 tsugi <- function(x){ #次に考えられる組合せ if (x[1]==0){ temp1 <- NULL }else{ temp1 <- expand.grid(0:(x[1]-1), x[2], x[3]) } if (x[2]==0){ temp2 <- NULL }else{ temp2 <- expand.grid(x[1], 0:(x[2]-1), x[3]) } if (x[3]==0){ temp3 <- NULL }else{ temp3 <- expand.grid(x[1], x[2], 0:(x[3]-1)) } return(rbind(temp1, temp2, temp3)) } mae <- function(x){ #前に考えられる組合せ if (x[1] == n){ temp1 <- NULL }else{ temp1 <- expand.grid((x[1]+1):n, x[2], x[3]) } if (x[2] == n){ temp2 <- NULL }else{ temp2 <- expand.grid(x[1], (x[2]+1):n, x[3]) } if (x[3] == n){ temp3 <- NULL }else{ temp3 <- expand.grid(x[1], x[2], (x[3]+1):n) } return(rbind(temp1, temp2, temp3)) } same <- function(y){ which(apply(dat, 1, function(x) all(x == y))) } while(sum(flag==0) != 0){ cat(sum(flag==0),"st") if (sum(flag == 1 & flag2 == 0) != 0){ make.check <- dat[(flag == 1 & flag2 == 0),] flag2[flag == 1] <- 1 temp1 <- NULL for (i in seq_len(nrow(make.check))){ temp1 <- rbind(temp1, mae(make.check[i,])) } temp1 <- as.matrix(unique(temp1, 1, incomparables = FALSE)) for (i in seq_len(nrow(temp1))){ flag[apply(dat, 1, function(x) all(x == temp1[i,]))] <- 2 } } ##ここまで絶対負けのチェック if (sum(flag==0) != 0){ #未確定が残っていれば temp3 <- dat[flag == 0,] #未確定のもの cat(sum(flag==0),"make", temp3[1,]) for (k in seq_len(nrow(temp3))){ temp4 <- tsugi(temp3[k,]) if( all(flag[apply(temp4, 1, same)] == 2)){ #次がすべて負の場合 flag[same(temp3[k,])] <- 1 #勝ちフラグ達成 } } } } return(dat[flag==1,]) ##勝ち } library(rbenchmark) benchmark(f1(n=5), f2(N=5), f3(n=5), replications = 10) # test replications elapsed relative user.self sys.self user.child sys.child #1 f1(n = 5) 10 0.26 1.857 0.26 0.00 NA NA #2 f2(N = 5) 10 0.14 1.000 0.14 0.00 NA NA #3 f3(n = 5) 10 40.52 289.429 40.17 0.01 NA NA #1 f1,f2=友人 f3=自分
旧友0.14秒と自分45.2秒。圧倒的な差ですね。
しかもnが増えるとますます差が開くという。
C++導入
visual studioの無料版を導入した。
10GBとかでて、気が遠くなった。
とりあえず、ちょっと触ってみた。
#include <stdio.h> #include <limits.h> #include <iostream> int main(){ std::cout << INT_MAX << "INT_MAX\n"; std::cout << LONG_MAX << "LONG_MAX\n"; std::cout << CHAR_BIT << "CHAR_BIT\n"; std::cout << LLONG_MAX << "LLONG_MAX\n"; system("pause"); return(0); }
2147483647INT_MAX 2147483647LONG_MAX 8CHAR_BIT 9223372036854775807LLONG_MAX 続行するには何かキーを押してください . . .
とりあえず、int型もlong型もうちの環境だと同じだった。
RからC＋＋を叩いていた時よりも、補完やエラー表示が親切。
理解しないととんでもないミスに繋がりそうなので、Rにはない概念
例えばポインタなどを勉強していきたい。
でもその前に、Visual studio自体を理解せねば。
それにしてもC言語を触っていると、Rが裏でエラーが出ないように
とっても頑張ってくれているということが身にしみますな。
Rcppを使う：その２
C言語による乱数生成
ここにメルセンヌツイスターのソースコードがあるので利用してみる。
あらかじめMT.hをincludeフォルダにおいておく必要がある。
以下のファイルをpi4.cppとして保存
#include <Rcpp.h> #include <stdio.h> #include <MT.h> using namespace Rcpp; // [[Rcpp::export]] double pi4Cpp(){ long circle_in = 0; init_genrand(10); for (long i = 0; i < 100000000; ++i) { double l = pow(genrand_real3(),2) + pow(genrand_real3(),2); if (l < 1.0){ circle_in++; } } //return(circle_in); return(4.0 * circle_in / 100000000); }
そして、以下のRコードを実行する。
library(Rcpp) library(rbenchmark) sourceCpp("pi4.cpp") cppFunction('double pi3Cpp(){ long circle_in = 0; Rcpp::NumericVector x = runif(200000000); for (long i = 0; i < 100000000; ++i) { double l = x[i]*x[i] + x[i+100000000]*x[i+100000000]; if (l < 1.0){ circle_in++; } } return(4.0 * circle_in / 100000000); }') f0 <- function(){ n = 100000000 x <- runif(n = n) y <- runif(n = n) print(4*sum(x^2 + y^2 <= 1)/n) } benchmark(f0(),pi3Cpp(),pi4Cpp(),replications = 1)
> benchmark(f0(),pi3Cpp(),pi4Cpp(),replications = 1)
[1] 3.141453
[1] 3.141787
test replications elapsed relative user.self sys.self user.child sys.child
1 f0() 1 11.09 8.338 10.17 0.83 NA NA
2 pi3Cpp() 1 2.87 2.158 2.59 0.26 NA NA
3 pi4Cpp() 1 1.33 1.000 1.32 0.00 NA NA
> １１秒だったのが、Rcppを使用して毎ループでメルセンヌツイスター乱数を
呼び出しても、１．３３秒まで短くなった。
Rcppを試す
Rcppを試す
https://rs-fan-jp.shinyapps.io/surv_jp/
H22 survival data in Japan
前に自作したWebアプリで、すでに人生の折り
返し時点を過ぎてしまっていることに気がついた。
親の平均余命もなんか具体的すぎてグッと来る。
時間は無限ではないので、ここらで更に新しい
言語に挑戦するのも良いんじゃないかと思った。
まずはお手軽に、Rccpから挑戦してみる。
http://www.slideshare.net/masakitsuda940/rcpp-36654397
Rcppのすすめ from Masaki Tsuda
Rcppのすすめここを参考にする。
http://rpubs.com/kaz_yos/pi-calcRPubs - pi calculation
このpi計算をやってみる
library(Rcpp) library(rbenchmark) cppFunction('double pi3Cpp(){ long circle_in = 0; Rcpp::NumericVector x = runif(200000000); for (long i = 0; i < 100000000; ++i) { double l = x[i]*x[i] + x[i+100000000]*x[i+100000000]; if (l < 1.0){ circle_in++; } } return(4.0 * circle_in / 100000000); }') f0 <- function(){ n = 100000000 x <- runif(n = n) y <- runif(n = n) print(4*sum(x^2 + y^2 <= 1)/n) } rbenchmark(f0(),pi3Cpp()) > benchmark(f0(),pi3Cpp(),replications = 1) [1] 3.141575 [1] 3.1417 test replications elapsed relative user.self sys.self user.child sys.child 1 f0() 1 9.43 3.532 8.98 0.43 NA NA 2 pi3Cpp() 1 2.67 1.000 2.46 0.17 NA NA >
9.4秒が2.6秒になった。わーい。
ただ使ってみて良くわかったのは、コンパイルエラーは難しい。
こんな単純なことでもエラーで苦労した。
慣れているとはいえ、なんとRは簡単な事か。
本当はC++の一様乱数を使用したかったが、エラーで使えなかった。
#include <random>してもstd::mt19937がエラーで読めない。
結局Rのrunif()を使っている。この呼び出し時間がかかるっぽい。
毎回呼び出すと8.15秒で、ベクトル化したRと変わらん感じだった。
なんだか素のc++を勉強したくなった。
暗殺教室の数学問題をRで描画する
暗殺教室の問題をRで描画する
追記
hoxo_m さんが、ブラウザでグリグリうごかせるファイルをアップしてくれました。
WebGLファイルの作成方法は知っていましたが、なるほどドロップボックスに
出来たファイルを置けばみんなハッピーにうごかせるのか～。
職場にジャンプがころがっていた。
暗殺教室で数学の問題があり、その解の三次元の形状が気になった。
http://maji-w.blog.jp/archives/20064241.html
【画像】暗殺教室の数学問題が難しすぎると話題にｗｗｗｗｗｗｗｗｗｗｗｗｗｗｗｗ : マジワロ速報
Rで描画して３Dでぐるぐる回してみたい。
隣り合う立方体形の各単位格子の中心点同士からの等距離の面は、
立方体形の各単位格子の面そのものである。
したがって、注目する１つの立方体形の単位格子の内部のみ
考えれば良い。（漫画には単位格子の頂点からの距離ばかり考えて
中心点同士の考察が書いてなかったのが気になる）
この中の点のうち、条件をみたすものを力技で判定して描画する。
library(rgl) p0 <- c(0,0,0) temp <- c(0.5, -0.5) p <- expand.grid(temp,temp,temp) dist.p0 <- function(x) sum(x^2) dist.p <- function(x, y){ sum((x - y) ^ 2) } dist.p1 <- function(y){ apply(p, 1, function(x)dist.p(x, y)) } n <- 50 x <- seq(-0.5, 0.5, length = n) x <- expand.grid(x,x,x) p0.range <- function(x){ all(dist.p1(x) > dist.p0(x)) } temp3 <- x[apply(x, 1, p0.range),] nrow(temp3) / n^3 * 100 #何％の点が条件をみたしたか p2 <- matrix(0, nrow=3, ncol=3) diag(p2) <- 1 p2 <- rbind(p2, -p2) dist.p2 <- function(y){ apply(p2, 1, function(x)dist.p(x, y)) } posible.outer <- function(x){ epsiron <- 0.02 any(c(abs(dist.p1(x) - dist.p0(x)) < epsiron, abs(dist.p2(x) - dist.p0(x)) < epsiron)) } temp4 <- temp3[apply(temp3, 1, posible.outer),] col.calc <- function(x){ epsiron <- 0.02 t1 <- abs(dist.p1(x) - dist.p0(x)) < epsiron t2 <- abs(dist.p2(x) - dist.p0(x)) < epsiron sum(t1 * 1:8) + sum(t2 * 9:14) } col0 <- apply(temp4, 1, col.calc) col1 <- rainbow(14) col <- col1[col0] col[is.na(col)] <- "#00000000" plot3d(temp4, col=col)
> nrow(temp3) / n^3 * 100 #何％の点が条件をみたしたか
[1] 48.4992
大体５０％になった。
マウスで、ぐりぐり回せるので楽しい。
