廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

シミュレーションで信頼区間をつくる

Data Analysis Using Regression and Multilevel/HierarchicalModels (Analytical Methods for Social Research)

Data Analysis Using Regression and Multilevel/HierarchicalModels (Analytical Methods for Social Research)

↑を少し読み直してる

頭が熱を帯びたみたいにぼーっとして、脳の働きの一部に鍵がかかったような感じがする。うまく考えられない。つらいなー
とか思いながら一方で、他人の思考を全然信用していなくて、自分の能力を信頼してるような気持ちもあって、なんなの?どうしたいの?おれは何がしたいの?念仏唱えて他力で救われようと思ったけど、歎異抄が見つからない。ふざけんな


さて、本書(歎異抄じゃなくてData Analysis Using Regression and Multilevel/Hierarchical Modelsのほう)の19〜20ページにある例:
男性500人
女性500人
合計1000人に聞きました。死刑賛成ですか?
調査の結果、
賛成700人、反対300人(「わからない」や「無回答」はここではないものとする)
そして、男性は75%が死刑肯定派、女性は65%が死刑肯定派だったとする

0.75/0.65で男性のほうが、女性より約1.15倍、死刑を肯定している、と言えそうだ
信頼区間をつくってみたいが、まともにやると大変なので、シミュレーションを使う。こんなふうだ

 n.men <- 500
 p.hat.men <- 0.75
 se.men <- sqrt (p.hat.men*(1-p.hat.men)/n.men)
#標準誤差
 
 n.women <- 500
 p.hat.women <- 0.65
 se.women <- sqrt(p.hat.women*(1-p.hat.women)/n.women)
 
 n.sims <- 10000
 p.women <-rnorm(n.sims, p.hat.women,se.women)
 p.men <- rnorm(n.sims,p.hat.men, se.men)
 ratio <- p.men/p.women
 int.95  <- quantile (ratio, c(.025,.975))
#ratioの中身を小さい順にならべて、上から2.5%番目と97.5%番目をとる(10000個だから250番目と9750番目)

テキストでは、rnorm(正規乱数)を使っている
しかし、賛成反対の二値的な回答の個数なんだから二項分布(rbinom)でいいんじゃないの?
というより、二項分布が正規分布で近似できるから正規分布でやっているのだと思うけど、どうせRつかうんだったら二項分布でやればいいんじゃないか。こんなふうに

 pbwomen<-rbinom(10000,500,0.65)
 pbmen<-rbinom(10000,500,0.75)
 ratio2 <- pbmen/pbwomen
 newint.95  <- quantile (ratio2, c(.025,.975))

ここでやっていることは、僕なりに例えると

  • 75%の確率で「死刑支持」と書いてあるおみくじが出てくる(残り25%は「死刑反対」とみなす)『男性のつぼ』と
  • 65%の確率で「死刑支持」と書いてあるおみくじが出てくる『女性のつぼ』

を用意してそこから500回くじを引き、「支持」の回数をカウント
(男性の支持数)/(女性の支持数)を計算
これを10000回やる


である。
乱数rbinomが「くじの入ったつぼ」にあたる


テキストの方は、うーん、同じたとえはつかえない
平均が\hat p_{men}=0.75
標準偏差\sqrt{\frac{\hat p_{men}(1-\hat p_{men})}{n_{men}}}
の男性のつぼ
と、
平均が\hat p_{women}=0.65
標準偏差\sqrt{\frac{\hat p_{women}(1-\hat p_{women})}{n_{women}}}
の女性のつぼ

だけど正規分布のつぼ(rnorm)からでてくるのは、「支持」「不支持」ではなく「支持率」「不支持率」が書いてあるくじである


この標準偏差
\sqrt{ \hat p_{men}(1-\hat p_{men})/n_{men}}
(コードだとse.men <- sqrt (p.hat.men*(1-p.hat.men)/n.men))
のところは、母比率pの推定量の標準誤差からこうなっている


うん、だから、直接二項分布をつかうやり方のほうが自然、っていうか普通だと思う



結果をくらべてみると

> int.95
    2.5%    97.5% 
1.064751 1.254016 

> newint.95
    2.5%    97.5% 
1.063218 1.254194 

ほとんどかわらない
(男性の母死刑支持率)/(女性の母死刑支持率)の95%信頼区間は[1.06,1.25]
一応考察っぽいことを書いておく:区間の下の方が1より大きいから、男性の方が死刑を支持する傾向が強いと考えてもよさそうだ


テキストで正規分布を使っている理由は不明だが、たぶんそのほうがどこかほかの部分を説明するときに都合がよくなるからかな?