12/1 ( ) GLM, R MCMC, WinBUGS 12/2 ( ) WinBUGS WinBUGS 12/2 ( ) : 12/3 ( ) :? ( :51 ) 2/ 71

Size: px

Start display at page:

Download "12/1 ( ) GLM, R MCMC, WinBUGS 12/2 ( ) WinBUGS WinBUGS 12/2 ( ) : 12/3 ( ) :? ( :51 ) 2/ 71"

はなとこたに
5 years ago
Views:

1 ( :51 ) 1/ 71 GCOE WinBUGS kubo@ees.hokudai.ac.jp

2 12/1 ( ) GLM, R MCMC, WinBUGS 12/2 ( ) WinBUGS WinBUGS 12/2 ( ) : 12/3 ( ) :? ( :51 ) 2/ 71

3 ( :51 ) 3/ 71 ( ) random effects BUGS

4 ( :51 ) 4/ 71

5 ( :51 ) 5/ 71 GLM (1)? y = 0, 1, 2, 3, (y ) (family = poisson) y = {0, 1}, y = {0, 1, 2,, N} (family = binomial) (family = Gamma) (family = gaussian)

6 R : glm() ( ) rbinom() glm(family = binomial) rbinom() glm(family = binomial) rpois() glm(family = poisson) rnbinom() glm.nb() ( ) rgamma() glm(family = gamma) rnorm() glm(family = gaussian) glm() glm.nb() MASS library GLM ( :51 ) 6/ 71

7 ( :51 ) 7/ 71 GLM (2)!! GLMM! GLMM/ random effects GLM

8 ( :51 ) 8/ 71 (Poisson distribution) (Binomial distribution) ( ) (Normal distribution, Gaussian ) (Gamma distribution)

9 ( :51 ) 9/ 71

10 ( :51 ) 10/ 71

11 ( :51 ) 11/ 71 : + WinBUGS 1. : GLMM 2.

12 1. : GLMM

13 : 8 ( ) ( ) : N i = 8 i y i = 3 : :? ( :51 ) 13/ 71

14 ( :51 ) 14/ 71 :?! ! y i

15 ( :51 ) 15/ 71 (overdispersion) y i 0.5 : overdispersion :?

16 ? : 1. fixed effects 2. random effects ( :51 ) 16/ 71

17 ( :51 ) 17/ 71 fixed random? fixed/random effects : 1. fixed effects : ( ) fixed effects 2. random effects : fixed effects ( ) random effects

18 ( :51 ) 18/ 71 : i N i y i ( ) Ni p(y i q i ) = q y i i (1 q i) N i y i, y i q i

19 ( :51 ) 19/ 71 q i = q(z i ) (logistic) q(z) = 1/{1 + exp( z)} q(z) z z i = a + b i a: b i : i ( )

20 b i (a {b 1, b 2,, b 100 }) /! ( )?? ( ) ( :51 ) 20/ 71

21 ( :51 ) 21/ 71 : b i s p(b i s) = 1 2πs 2 exp b2 i 2s 2, = 1 = 1.5 = 3 b i {b 1, b 2,, b 100 }

22 ( :51 ) 22/ 71 b i s = 1 = 1.5 = 3 b i s b i s b i y i

23 (C) : b i s s ( :51 ) 23/ 71 b i? (A) (B) (C)!? s (A) : b i (B) : b i ( -5 5 )

24 y i = 2 b i b i p(b i y i = 2, s) s p(y i = 2 b i ) b i p(b i s) s b i ( :51 ) 24/ 71

25 y i {2, 3, 5} b i s s s b i ( :51 ) 25/ 71

26 ( :51 ) 26/ 71? データ N[i] Y[i] 胚珠数中の生存 tau は hyper parameter 二項分布生存確率 q[i] 全体の平均 a 無情報事前分布植物の個体差 b[i] 事前分布 tau 個体差のばらつき無情報事前分布 ( 超事前分布 )

27 ( :51 ) 27/ 71 データ N[i] Y[i] 胚珠数中の生存 tau は hyper parameter 二項分布生存確率 q[i] 全体の平均 a 無情報事前分布植物の個体差 b[i] 事前分布 tau 個体差のばらつき無情報事前分布 ( 超事前分布 )

28 ? b i s = 0.1? s = 0.1 s individual posterior prior small large ( ) s ( :51 ) 28/ 71

29 ( :51 ) 29/ 71 τ = 1/s 2 s τ (non-informative prior) p(τ ) = τ α 1 e τ β Γ(α)β α, α = β = 10 4

30 ( :51 ) 30/ 71 (1) τ ( 1; 1) ( 1; 100) τ τ

31 (2) a ( 0; 1) ( 0; 100) a (logit) a ( :51 ) 31/ 71

32 ( :51 ) 32/ 71 p(a, {b i }, τ ) = 100 i=1 p(y i q(a + b i )) p(a) p(b i τ ) h(τ ) ( ) db i dτ da p(a, {b i }, τ ) 100 i=1 p(y i q(a + b i )) p(a) p(b i τ ) h(τ ) : : p(a, {b i }, τ ) 100 i=1 p(y i q(a + b i )) : p(a) p(b i τ ) h(τ )

33 ( :51 ) 33/ 71 b i s individual small hyperprior posterior prior (posterior) large hyperparameter s MCMC

34 ? p(a, {b i }, τ ) 100 i=1 p(y i q(a + b i )) p(a) p(b i τ ) h(τ ) p(a, {b i }, τ ) Markov chain Monte Carlo (MCMC)! WinBUGS ( :51 ) 34/ 71

35 Gibbs sampling p(a ) 100 p(y i q(a + b i )) p(a) i=1 p(τ ) 100 i=1 p(b i τ ) h(τ ) p(b 1 ) p(y 1 q(a + b 1 )) p(b 1 τ ) p(b 2 ) p(y 2 q(a + b 2 )) p(b 2 τ ). p(b 100 ) p(y 100 q(a + b 100 )) p(b 100 τ ) ( :51 ) 35/ 71

36 ( :51 ) 36/

37 ( :51 ) 37/ 71 :

38 ( :51 ) 38/ 71 ( ) ( ) ( ) ( ) データ N[i] Y[i] 胚珠数中の生存 tau は hyper parameter 二項分布生存確率 q[i] 全体の平均 a 無情報事前分布植物の個体差 b[i] 事前分布 tau 個体差のばらつき無情報事前分布 ( 超事前分布 ) : Markov Chain Monte Carlo (MCMC)

39 ( :51 ) 39/ 71

40 ( :51 ) 40/ 71 WinBUGS

41 ( :51 ) 41/ 71 WinBUGS BUGS (model.bug.txt) 3. R2WBwrapper R (runbugs.r) 4. R runbugs.r (source(runbugs.r) ) 5. bugs

42 ( :51 ) 42/ 71? データ N[i] Y[i] 胚珠数中の生存 tau は hyper parameter 二項分布生存確率 q[i] 全体の平均無情報事前分布 p(a, {b i }, τ ) 100 i=1 a 植物の個体差 b[i] 事前分布 tau 個体差のばらつき無情報事前分布 ( 超事前分布 ) p( q(a + b i )) p(a) p(b i τ ) h(τ )

43 BUGS model.bug.txt ( ) model{ for (i in 1:N.sample) { Y[i] ~ dbin(q[i], N[i]) # logit(q[i]) <- a + b[i] # q[i] } a ~ dnorm(0, 1.0E-4) # for (i in 1:N.sample) { b[i] ~ dnorm(0, tau) # } tau ~ dgamma(1.0e-4, 1.0E-4) # sigma <- sqrt(1 / tau) # tau SD } ( :51 ) 43/ 71

44 ( :51 ) 44/ 71 (hierarchical) random effects (non-informative) fixed effects (subjective) ( )

45 BUGS BUGS ( ) node ( ) 1. ~ sthochastic node 2. <- deterministic node ( :51 ) 45/ 71

46 R2WBwrapper R runbugs.r ( ) source("r2wbwrapper.r") # R2WBwrapper d <- read.csv("data.csv") # clear.data.param() # ( ) set.data("n.sample", nrow(d)) # set.data("n", d$n) # set.data("y", d$y) # ( :51 ) 46/ 71

47 R2WBwrapper R runbugs.r ( ) set.param("a", 0) # set.param("sigma", NA) # set.param("b", rep(0, N.sample)) # set.param("tau", 1, save = FALSE) # set.param("p", NA) # post.bugs <- call.bugs( # WinBUGS file = "model.bug.txt", n.iter = 2000, n.burnin = 1000, n.thin = 5 ) ( :51 ) 47/ 71

48 WinBUGS post.bugs <- call.bugs( file = "model.bug.txt", # WinBUGS n.iter = 2000, n.burnin = 1000, n.thin = 5 ) default ( ) 3 (n.chains = 3) MCMC sampling ( ) cf. PC MCMC chain 2000 step (n.iter = 2000) 1000 step (n.burnin = 1000) step 5 step (n.thin = 5) ( :51 ) 48/ 71

49 ? R source("runbugs.r") WinBUGS MCMC sampling (WinBUGS ) WinBUGS WinBUGS R post.bugs ( :51 ) 49/ 71

50 ( :51 ) 50/ 71 R a a?

51 ( :51 ) 51/ 71 bugs post.bugs (1) plot(post.bugs), R-hat Gelman-Rubin var ˆ ˆR + (ψ y) = W var ˆ + (ψ y) = n 1 n W + 1 n B W : chain variance B : chain variance Gelman et al Bayesian Data Analysis. Chapman & Hall/CRC

52 ( :51 ) 52/ 71 ubothinkpad/public_html/stat/2009/ism/winbugs/model.bug.txt", fit using WinBUGS, 3 chains, each with 1300 i 80% interval for each chain R-hat a sigma * b[1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] tau * q[1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] * array truncated for lack of space a sigma * b tau * q deviance medians and 80% intervals

53 mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff a sigma b[1] b[2] b[3] b[4] b[5] b[6] b[7] b[8] b[9] b[10] ( :51 ) 53/ 71 bugs post.bugs (2) print(post.bugs, digits.summary = 3) 95%

54 ( :51 ) 54/ 71 mcmc.list post.list <- to.list(post.bugs) plot(post.list[,1:4,], smooth = F),

55 ( :51 ) 55/ 71

56 ( :51 ) 56/ 71 mcmc post.mcmc <- to.mcmc(post.bugs) matrix :

57 2.

58 ( :51 ) 58/ 71 : abundance ( ) location :

59 ( :51 ) 59/ 71 : abundance location ( 95% )

60 β : β Normal(0, 10 2 ), ( :51 ) 60/ i λ i : y i Poisson(λ i ) λ i ( ) + ( ) : log λ i = β + r i

61 ( :51 ) 61/ 71 ( ) Conditional Autoregressive (CAR) r i (N i i, J i i ): j J r i Normal( i r j, σ ) σ N i N i σ : τ = 1/σ Gamma(1.0 2, ) p(β, {r i }, τ {y i }) = p({y i} β, {r i }, τ ) ( ( )dβdr1 dr 50 dτ

62 τ - τ (σ ) - τ (σ ) - τ ( ) abundance tau = abundance tau = tau = location location ( :51 ) 62/ 71

63 ( :51 ) 63/ 71 BUGS model { # BUGS for (i in 1:N.site) { Y[i] ~ dpois(mean[i]) # log(mean[i]) <- beta + re[i] # ( ) + ( ) } # re[i] CAR model re[1:n.site] ~ car.normal(adj[], Weights[], Num[], tau) beta ~ dnorm(0, 1.0E-2) # tau ~ dgamma(1.0e-2, 1.0E-2) # }

64 ( :51 ) 64/ 71 abundance β beta τ location tau ( 95% )

65 ( :51 ) 65/ 71 abundance location GLMM OK?

66 ( :51 ) 66/ 71 vs abundance abundance location location?

67 ( :51 ) 67/ 71 ( ):?!! abundance location (

68 ( :51 ) 68/ 71 abundance abundance location location!

69 ( :51 ) 69/ 71 abundance abundance location location CAR

70 ( :51 ) 70/ 71 : abundance abundance location location

71 ( :51 ) 71/ 71 : 1. : GLMM 2.

kubo2015ngt6 p.2 ( ( (MLE 8 y i L(q q log L(q q 0 ˆq log L(q / q = 0 q ˆq = = = * ˆq = 0.46 ( 8 y 0.46 y y y i kubo (ht

kubo2015ngt6 p.2 ( ( (MLE 8 y i L(q q log L(q q 0 ˆq log L(q / q = 0 q ˆq = = = * ˆq = 0.46 ( 8 y 0.46 y y y i kubo (ht kubo2015ngt6 p.1 2015 (6 MCMC kubo@ees.hokudai.ac.jp, @KuboBook http://goo.gl/m8hsbm 1 ( 2 3 4 5 JAGS : 2015 05 18 16:48 kubo (http://goo.gl/m8hsbm 2015 (6 1 / 70 kubo (http://goo.gl/m8hsbm 2015 (6 2 /