rstan神がやっている人口ピラミッドのAge Heapingを階層ベイズで補正するという解析について、自己申告に基づくデータにheapingが見られるという別の話があったので、3番煎じだがrstanを使わねば!!
上記のブログにある分布は、ギザギザしすぎているし、データを作成するにしても難しいので竿の直径のデータを使う。
最大で5cm くらいだが、いくつかのポイントで突出している。元々は50万人の自己申告に基づくデータを収集したらしく、元データほしいと思いながらも、分布の尖っているところをlocator 関数でポチポチする。
0.017 0.000 1.062 0.111 1.140 0.590 1.238 0.093 1.365 0.119 1.414 0.198 1.521 0.119 1.863 0.198 1.970 1.707 2.088 0.303 2.215 0.224 2.517 0.320 2.586 0.608 2.723 0.355 2.908 0.425 2.996 3.190 3.152 0.555 3.279 0.634 3.455 0.555 3.504 1.218 3.670 0.651 3.846 0.555 3.963 0.643 3.992 2.195 4.139 2.056 4.197 0.198 4.324 0.276 4.461 0.154 4.529 0.381 4.656 0.119 4.774 0.398 4.852 0.093 4.940 0.041 4.989 0.241 5.116 0.041
その後、線形補間して適当に水増しする。
# n にデータが入っている l <- function(x, v1, v2){ coef <- solve(n[c(i+1, i),], c(1, 1)) -coef[1]/coef[2] + 1/coef[2] } x <- seq(0, 5, by=0.05) # x を水増し y <- rep(0, length(x)) # x に対応する点 for(j in 2:length(x)){ i <- tail(which(n[,1] < x[j]), 1) coef <- solve(n[c(i+1, i),], c(1, 1)) # ax + by = 1 を解く y[j] <- -coef[1]/coef[2]*x[j] + 1/coef[2] } plot(x, y, type="l", lwd=3, col=4, xlab="Penile diameter (cm)", ylab="人数 (万人)") abline(v=0:5, h=0:4);abline(v=0:4+0.5, lty=3, h=0:3+0.5)
さて、ここでデータについて考えよう。
局部についてのデータなので、かなりデリケートである。匿名のアンケート、データ収集といえども、羞恥心や見栄をはりたい、という気持ちはよくわかる。
さらに、局部をものさしなどの測定器具でしっかり挟んで直径を測定しない限り、ミリ単位では計測できないし、つまんだりにぎったりしたときのだいたいの大きさでcm を報告してしまうかもしれない。せめて、なんとなくわかりそうな0.5cm 刻みでは報告しそうである。
ということを考慮すると、rstan神の年齢解析から改変すると、
・下のカテゴリーから上へ過剰申告する
・2.4cm の人は、2.1とか2.2cm の人よりかは2.5cm にしてしまったほうが見栄えがいいので、よく2.5cm と報告する(これはまあ、2.1cmでもより上に報告するとかはありそうだけれども)
さらに、0.5cm 刻みだと楽なんだが、グラフでは0.5cm 刻みより外れていそうなところがあるので
・Ncm + 0.5cm + cm のなんらかの項が存在する。
みたいなモデルを考えたい。
けどとりあえず、0.5cm 刻みで、データとしては0.05cm 刻みで報告されているデータの、とあるN.5cm より小さい5群くらいが過剰報告してくる(つまり、fromが常に下側、上のtoへ流入してくる)とする。
やってみたけど、全然滑らかにならなかった。もっといじりたかったけど時間がないし、こういうのは旬が大事()なので。
95% 信用区間と推定値。
diam_idx <- which(x %in% seq(from=0, 5.5, by=0.5)) x_idx <- cut(x, diam, include.lowest=TRUE) li <- lapply(-(5:1), function(j) data.frame(to=diam_idx, from=diam_idx + j)) diam_heap <- subset(do.call(rbind, li), from > 0) stanmodel <- stan_model("penile_stan.stan") # rstan 神のスクリプトそのまま dat <- list(A=length(x), Y=round(y*10000), J=nrow(diam_heap), From=diam_heap$from, To=diam_heap$to) fit <- sampling(stanmodel, data=dat, chain=3, seed=123) ex0 <- extract(fit) qu <- apply(ex0$q, 2, quantile, c(0.025, 0.5, 0.975))*sum(y) plot(x, y, type="n", xlab="Penile diameter (cm)", ylab="推定") polygon(c(x, rev(x)), c(qu[1,], rev(qu[3,])), border=NA, col=grey(0.9)) lines(x, qu[2,], lty=3, lwd=3, col=4) abline(v=0:5, h=0:4);abline(v=0:4+0.5, lty=3, h=0:3+0.5)
