kubo2017sep16a p.1 ( 1 ) : : :55 kubo ( ( 1 ) / 10

kubo2017sep16a p.1 ( 1 ) kubo@ees.hokudai.ac.jp 2017 09 16 : http://goo.gl/8je5wh : 2017 09 13 16:55 kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 1 / 106 kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 2 / 106 :! http://goo.gl/8je5wh 階層ベイズモデル もっと自由な統計モデリングを! 線形モデルの発展 (HBM) 一般化線形混合モデル 個体差 場所差といった変量効果をあつかいたい 一般化線形モデル 正規分布以外の確率分布をあつかいたい (GLMM) 推定計算方法 MCMC 最尤推定法 (GLM) 最小二乗法 線形モデル : GLM GLIM General Linear Model kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 3 / 106 kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 4 / 106 () 1 2 :? 3 GLM: GLM R GLM 4 GLM GLM 1. kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 5 / 106 kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 6 / 106

kubo2017sep16a p.2 : :?! ( ) method section p < 0.05 OK kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 7 / 106 kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 8 / 106? (): kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 9 / 106 kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 10 / 106 : 0, 1, 2 (y {0, 1, 2, 3, } ) y -2 0 2 4 6 0.5 1.0 1.5 2.0 x y? kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 11 / 106 x y -2 0 2 4 6 0.5 1.0 1.5 2.0? y 0? x NO! kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 12 / 106

kubo2017sep16a p.3? y -2 0 2 4 6 0.5 1.0 1.5 2.0 x YES! kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 13 / 106 階層ベイズモデル もっと自由な統計モデリングを! 線形モデルの発展 (HBM) 一般化線形混合モデル 個体差 場所差といった変量効果をあつかいたい 一般化線形モデル 正規分布以外の確率分布をあつかいたい (GLMM) 推定計算方法 MCMC 最尤推定法 (GLM) 最小二乗法 線形モデル kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 14 / 106 : : 1. 2. (GLM) 3. MCMC 2. : : 1, 2, 3 4. GLM kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 15 / 106 kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 16 / 106 : : () : : kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 17 / 106 kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 18 / 106

kubo2017sep16a p.4 : () : R i 50 i {1, 2, 3,, 50} y i {y i }! {y i } = {y 1, y 2,, y 50 } {y i } R > data [1] 2 2 4 6 4 5 2 3 1 2 0 4 3 3 3 3 4 2 7 2 4 3 3 3 4 [26] 3 7 5 3 1 7 6 4 6 5 2 4 7 2 2 6 2 4 5 4 5 1 3 2 3! kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 19 / 106 kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 20 / 106 : R : table() > table(data) 0 1 2 3 4 5 6 7 1 3 11 12 10 5 4 4 ( 5 5 6 4 ) > hist(data, breaks = seq(-0.5, 9.5, 1))! kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 21 / 106 kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 22 / 106 : > mean(data) [1] 3.56 > abline(v = mean(data), col = "red") kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 23 / 106 : : > var(data) [1] 2.9861 (SD = variance) > sd(data) [1] 1.7280 > sqrt(var(data)) [1] 1.7280 : : kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 24 / 106

kubo2017sep16a p.5 :? :?? kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 25 / 106 kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 26 / 106 :? :? > data.table <- table(factor(data, levels = 0:10)) > cbind(y = data.table, prob = data.table / 50) y prob 0 1 0.02 1 3 0.06 2 11 0.22 3 12 0.24 4 10 0.20 5 5 0.10 6 4 0.08 7 4 0.08 8 0 0.00 9 0 0.00 10 0 0.00 kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 27 / 106 kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 28 / 106 :? :?? () y p(y λ) = λy exp( λ) y! y! y 4! 1 2 3 4 exp( λ) = e λ (e = 2.718 ) kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 29 / 106 kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 30 / 106

kubo2017sep16a p.6 :? :?? R > y <- 0:9 # () > prob <- dpois(y, lambda = 3.56) # > plot(y, prob, type = "b", lty = 2) > # cbind > cbind(y, prob) (λ) 3.56 Poisson distribution y prob 1 0 0.02843882 2 1 0.10124222 3 2 0.18021114 4 3 0.21385056 5 4 0.19032700 6 5 0.13551282 7 6 0.08040427 8 7 0.04089132 9 8 0.01819664 10 9 0.00719778 kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 31 / 106 > hist(data, seq(-0.5, 8.5, 0.5)) # > lines(y, prob, type = "b", lty = 2) # kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 32 / 106 :? :? λ? λ λ λ 0 : λ = = > # cbind > cbind(y, prob) y prob 1 0 0.02843882 2 1 0.10124222 3 2 0.18021114 4 3 0.21385056 5 4 0.19032700 6 5 0.13551282 7 6 0.08040427 8 7 0.04089132 9 8 0.01819664 10 9 0.00719778 y {0, 1, 2,, } y 1 p(y λ) = 1 y=0 kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 33 / 106 : y i {0, 1, 2, } y i : kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 34 / 106 :? : λ p(y λ) = λy exp( λ) y! λ kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 35 / 106!? kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 36 / 106

kubo2017sep16a p.7 : (likelihood)? : L(λ) λ λ λ 3 {y 1, y 2, y 3 } = {2, 2, 4} 0.180 0.180 0.19 = 0.006156 (the likelihood definition for the example): L(λ) = (y 1 2 ) (y 2 2 ) (y 50 3 ) = p(y 1 λ) p(y 2 λ) p(y 3 λ) p(y 50 λ) = p(y i λ) = λ yi exp( λ), y i i i! kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 37 / 106 kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 38 / 106 : : λ () (!) (log likelihood function) log L(λ) = i ( ) y i y i log λ λ log k k log L(λ) L(λ) λ kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 39 / 106 kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 40 / 106 : : ˆλ log L(λ) = ( i yi log λ λ y i k log k) log likelihood -120-110 -100 * ˆλ = 3.56 2.0 2.5 3.0 3.5 4.0 4.5 5.0 (ML estimator): i y i/50 d log L dλ (ML estimate): ˆλ = 3.56 = 0! λ λ 3.5 50 0 100 300 3000 ˆλ 2.5 3.0 3.5 4.0 4.5 ˆλ λ 50 kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 41 / 106 kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 42 / 106

kubo2017sep16a p.8 : : : λ = 3.5 0 2 4 6 8 y ˆλ = 3.56 0 2 4 6 8 y 0 2 4 6 8 y kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 43 / 106 kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 44 / 106 : : λ = 3.5 0 2 4 6 8 y ˆλ = 3.56 0 2 4 6 8 y : ( 0 2 4 6 8 0 2 4 6 8 0 2 4 6 8 : y {0, 1, 2, 3, } y : y {0, 1, 2,, N} N y : < y < kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 45 / 106 kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 46 / 106 : : GLMM 線形モデルの発展 階層ベイズモデル もっと自由な統計モデリングを! (HBM) 一般化線形混合モデル 個体差 場所差といった変量効果をあつかいたい 一般化線形モデル 正規分布以外の確率分布をあつかいたい (GLMM) 推定計算方法 MCMC 最尤推定法 (GLM) 最小二乗法線形モデル? (GLM) (Poisson regression) (logistic regression) (linear regression) kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 47 / 106 kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 48 / 106

kubo2017sep16a p.9 GLM: GLM: i 3. GLM:? : {y i } : {x i } {f i } f i C: T: y i x i (f i = C): 50 sample (i {1, 2, 50}) (f i = T): 50 sample (i {51, 52, 100}) kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 49 / 106 kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 50 / 106 GLM: GLM: data frame d : data: http://hosho.ees.hokudai.ac.jp/~kubo/ce/eeslecture2017.html#toc4 data3a.csv CSV (comma separated value) format file R : > d <- read.csv("data3a.csv") d data frame () data frame d > d y x f 1 6 8.31 C 2 6 9.44 C 3 6 9.50 C... 99 7 10.86 T 100 9 9.97 T > d$x [1] 8.31 9.44 9.50 9.07 10.16 8.32 10.61 10.06 [9] 9.93 10.43 10.36 10.15 10.92 8.85 9.42 11.11... [97] 8.52 10.24 10.86 9.97 > d$y [1] 6 6 6 12 10 4 9 9 9 11 6 10 6 10 11 8 [17] 3 8 5 5 4 11 5 10 6 6 7 9 3 10 2 9... [97] 6 8 7 9 kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 51 / 106 kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 52 / 106 GLM: data frame d : GLM: R f > d$f [1] C C C C C C C C C C C C C C C C C C C C C C C C C [26] C C C C C C C C C C C C C C C C C C C C C C C C C [51] T T T T T T T T T T T T T T T T T T T T T T T T T [76] T T T T T T T T T T T T T T T T T T T T T T T T T Levels: C T : C T 2 > class(d) # d data.frame [1] "data.frame" > class(d$y) # y integer [1] "integer" > class(d$x) # x numeric [1] "numeric" > class(d$f) # f factor [1] "factor" kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 53 / 106 kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 54 / 106

kubo2017sep16a p.10 GLM: data frame summary() GLM:! Generate Data Plots! Always! > summary(d) y x f Min. : 2.00 Min. : 7.190 C:50 1st Qu.: 6.00 1st Qu.: 9.428 T:50 Median : 8.00 Median :10.155 Mean : 7.83 Mean :10.089 3rd Qu.:10.00 3rd Qu.:10.685 Max. :15.00 Max. :12.400 kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 55 / 106 > plot(d$x, d$y, pch = c(21, 19)[d$f]) > legend("topleft", legend = c("c", "T"), pch = c(21, 19)) kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 56 / 106 GLM: GLM: GLM f (box-whisker plot) > plot(d$f, d$y) # note that d$f is factor type! GLM kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 57 / 106 kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 58 / 106 GLM: GLM GLM: GLM GLM (GLM)??? : : e.g., β 1 + β 2 x i : () + () x i : -2 0 2 4 6 0.5 1.0 1.5 2.0 kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 59 / 106 kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 60 / 106

kubo2017sep16a p.11 GLM: GLM GLM: GLM (?) GLM : (response variable) : (explanatory variable) (linear predictor): () = (, intercept) + ( 1) ( 1) + ( 2) ( 2) : : e.g., β 1 + β 2 x i : -2 0 2 4 6 0.5 1.0 1.5 2.0 + ( 3) ( 3) + kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 61 / 106 kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 62 / 106 GLM: GLM GLM: GLM GLM logistic R (GLM) GLM () rbinom() glm(family = binomial) : : e.g., β 1 + β 2 x i : logit yi x i rbinom() glm(family = binomial) rpois() glm(family = poisson) rnbinom() glm.nb() in library(mass) () rgamma() glm(family = gamma) rnorm() glm(family = gaussian) glm() GLM kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 63 / 106 kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 64 / 106 GLM: GLM Go buck to the seed example y i λ i p(y i λ i ) = λyi i exp( λ i ) y i! i λ i? λ i = exp(β 1 + β 2 x i ) β 1 β 2 () x i i f i kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 65 / 106 GLM:? i λi GLM λ i = exp(β 1 + β 2 x i ) {β 1, β 2} = { 2, 0.8} {β 1, β 2} = { 1, 0.4} i x i kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 66 / 106

kubo2017sep16a p.12 GLM: GLM GLM: GLM GLM () i λ i λ i = exp(β 1 + β 2 x i ) log(λ i ) = β 1 + β 2 x i log() = : : β 1 + β 2 x i : log kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 67 / 106 kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 68 / 106 GLM: R GLM GLM: R GLM glm() R GLM > d y x f 1 6 8.31 C 2 6 9.44 C 3 6 9.50 C... 99 7 10.86 T 100 9 9.97 T! > fit <- glm(y ~ x, data = d, family = poisson) kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 69 / 106 kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 70 / 106 GLM: R GLM GLM: R GLM glm() glm() > fit <- glm(y ~ x, data = d, family = poisson) all: glm(formula = y ~ x, family = poisson, data = d) Coefficients: (Intercept) x 1.2917 0.0757 ( z):? link : z (y)? family:? Degrees of Freedom: 99 Total (i.e. Null); Null Deviance: 89.5 Residual Deviance: 85 AIC: 475 98 Residual kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 71 / 106 kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 72 / 106

kubo2017sep16a p.13 GLM: R GLM GLM: R GLM glm() > summary(fit) Call: glm(formula = y ~ x, family = poisson, data = d) Deviance Residuals: Min 1Q Median 3Q Max -2.368-0.735-0.177 0.699 2.376 Coefficients: Estimate Std. Error z value Pr(> z ) (Intercept) 1.2917 0.3637 3.55 0.00038 x 0.0757 0.0356 2.13 0.03358 () kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 73 / 106 () β 2 (Estimate 0.0757, SE 0.0356) (Estimate 1.29, SE 0.364) β 1 0.0 0.5 1.0 1.5 p p ˆβ p 0.5 ˆβ (: ) kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 74 / 106 GLM: R GLM GLM: R GLM (?) β 2 (Estimate 0.0757, SE 0.0356) (Estimate 1.29, SE 0.364) β 1 0.0 0.5 1.0 1.5 95%? kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 75 / 106 > fit <- glm(y ~ x, data = d, family = poisson)... Coefficients: (Intercept) x 1.2917 0.0757 > plot(d$x, d$y, pch = c(21, 19)[d$f]) # data > xp <- seq(min(d$x), max(d$x), length = 100) > lines(xp, exp(1.2917 + 0.0757 * xp)) kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 76 / 106 GLM: GLM: f i GLM + y i λ i i λ i β 3 f i p(y i λ i ) = λyi i exp( λ i ) y i! λ i = exp(β 1 + β 2 x i + β 3 d i ) d i = { 0 (f i = C ) 1 (f i = T ) kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 77 / 106 kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 78 / 106

kubo2017sep16a p.14 GLM: GLM: glm(y x + f,...) > summary(glm(y ~ x + f, data = d, family = poisson))...()... Coefficients: Estimate Std. Error z value Pr(> z ) (Intercept) 1.2631 0.3696 3.42 0.00063 x 0.0801 0.0370 2.16 0.03062 ft -0.0320 0.0744-0.43 0.66703 () x + f > plot(d$x, d$y, pch = c(21, 19)[d$f]) # data > xp <- seq(min(d$x), max(d$x), length = 100) > lines(xp, exp(1.2631 + 0.0801 * xp), col = "blue", lwd = 3) # C > lines(xp, exp(1.2631 + 0.0801 * xp - 0.032), col = "red", lwd = 3) # T kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 79 / 106 kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 80 / 106 GLM: GLM: λi f i = C: λ i = exp(1.26 + 0.0801x i ) f i = T: λ i = exp(1.26 + 0.0801x i 0.032) = exp(1.26 + 0.0801x i ) exp( 0.032) x i exp( 0.032)! kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 81 / 106 λi (A) (B) λ = exp(β 1 + β 2 x + ) λ = β 1 + β 2 x + x i x i kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 82 / 106 GLM: GLM: GLM: y -2 0 2 4 6 0.5 1.0 1.5 2.0 x log y -2 0 2 4 6 0.5 1.0 1.5 2.0 x 階層ベイズモデル もっと自由な統計モデリングを! 線形モデルの発展 (HBM) 一般化線形混合モデル 個体差 場所差といった変量効果をあつかいたい 一般化線形モデル 正規分布以外の確率分布をあつかいたい (GLMM) 推定計算方法 MCMC 最尤推定法 (GLM) 最小二乗法 線形モデル kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 83 / 106 kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 84 / 106

kubo2017sep16a p.15 GLM GLM : 4. GLM N i = 8 : 8 : i y i = 3 : 20 160 73 kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 85 / 106 kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 86 / 106 GLM example data: dead or alive of seeds 0 1 2 3 4 5 6 7 8 1 2 1 3 6 6 1 0 0 y i (no individual difference!) kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 87 / 106 GLM q q i N i y i ) p(y i q) = In this example... ( Ni y i q y i (1 q) N i y i, no individual difference q kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 88 / 106 GLM : 20 {y i } q 20 q L(q {y i }) = 20 i=1 p(y i q) 0 1 2 3 4 5 6 7 8 1 2 1 3 6 6 1 0 0 kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 89 / 106 GLM L(q ) ˆq 20 log L(q ) = log i=1 ( Ni 20 + {y i log(q) + (N i y i ) log(1 q)} i=1 q kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 90 / 106 y i )

kubo2017sep16a p.16 GLM (MLE) do you remember? GLM 8 y i L(q ) q log L(q ) q 0 ˆq log L(q )/ q = 0 q log likelihood -50-45 -40 ˆq = = 73 160 = 0.456 0.3 0.4 0.5 0.6 * ˆq = 0.46 ( ) 8 y 0.46 y 0.54 8 y y i kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 91 / 106 kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 92 / 106 GLM GLM GLM GLM GLM logistic GLM : : e.g., β 1 + β 2 x i yi : logit x i kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 93 / 106 kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 94 / 106 GLM GLM GLM GLM N y 8 y! f i C: T: i N i = 8 y i = 3 (alive) (dead) x i data4a.csv CSV (comma separated value) format file R : > d <- read.csv("data4a.csv") or > d <- read.csv( + "http://hosho.ees.hokudai.ac.jp/~kubo/stat/2014/fig/binomial/data4a.csv") d data frame () kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 95 / 106 kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 96 / 106

kubo2017sep16a p.17 GLM GLM GLM GLM data frame d > summary(d) N y x f Min. :8 Min. :0.00 Min. : 7.660 C:50 1st Qu.:8 1st Qu.:3.00 1st Qu.: 9.338 T:50 Median :8 Median :6.00 Median : 9.965 Mean :8 Mean :5.08 Mean : 9.967 3rd Qu.:8 3rd Qu.:8.00 3rd Qu.:10.770 Max. :8 Max. :8.00 Max. :12.440 > plot(d$x, d$y, pch = c(21, 19)[d$f]) > legend("topleft", legend = c("c", "T"), pch = c(21, 19)) yi x i? kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 97 / 106 kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 98 / 106 GLM GLM : N y p(y N, q) = ( ) N q y (1 q) N y y ( N ) y N y logit p(y i 8, q) 0.0 0.1 0.2 0.3 0.4 q = 0.1 q = 0.8 q = 0.3 0 2 4 6 8 y i kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 99 / 106 kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 100 / 106 GLM GLM (z i: e.g. z i = β 1 + β 2x i) 1 q i = logistic(z i ) = 1 + exp( z i ) > logistic <- function(z) 1 / (1 + exp(-z)) # > z <- seq(-6, 6, 0.1) > plot(z, logistic(z), type = "l") q 0.0 0.2 0.4 0.6 0.8 1.0 q = 1 1+exp( z) -6-4 -2 0 2 4 6 z kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 101 / 106 {β 1, β 2 } = {0, 2}(A) β 2 = 2 β 1 (B) β 1 = 0 β 2 q 0.0 0.2 0.4 0.6 0.8 1.0 (A) β 2 = 2 β 1 = 2 β 1 = 0-3 -2-1 0 1 2 3 x β 1 = 3 0.0 0.2 0.4 0.6 0.8 1.0 (B) β 1 = 0 β 2 = 4 β 2 = 2-3 -2-1 0 1 2 3 x β 2 = 1 {β 1, β 2 } x q 0 q 1 kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 102 / 106

kubo2017sep16a p.18 GLM logit link function GLM R β 1 β 2 logistic 1 q = 1 + exp( (β 1 + β 2 x)) = logistic(β 1 + β 2 x) logit q logit(q) = log 1 q = β 1 + β 2 x logit logistic logistic logit logit is the inverse function of logistic function, vice versa y (A) f i =C x (B) > glm(cbind(y, N - y) ~ x + f, data = d, family = binomial)... Coefficients: (Intercept) x ft -19.536 1.952 2.022 x kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 103 / 106 kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 104 / 106 GLM GLM : : yi (A) f i =C x i (B) f i =T x i 1. 2. (GLM) 3. MCMC 4. GLM kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 105 / 106 kubo (http://goo.gl/ufq2) ( 1 ) 2017 09 16 106 / 106

kubo2017sep16b p.1 ( 2 ) GLM kubo@ees.hokudai.ac.jp 2017 09 16 : http://goo.gl/8je5wh : 2017 09 13 16:58 kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 1 / 90 階層ベイズモデル もっと自由な統計モデリングを! 線形モデルの発展 (HBM) 一般化線形混合モデル 個体差 場所差といった変量効果をあつかいたい 一般化線形モデル 正規分布以外の確率分布をあつかいたい (GLMM) 推定計算方法 MCMC 最尤推定法 (GLM) 最小二乗法 線形モデル kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 2 / 90 : 1. 2. (GLM) 3. MCMC () 1 MCMC MCMC: 2 MCMC 3 4 4. GLM kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 3 / 90 kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 4 / 90 MCMC MCMC : () 1. MCMC : () 0 1 2 3 4 5 6 7 8 1 2 1 3 6 6 1 0 0 y i kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 5 / 90 kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 6 / 90

kubo2017sep16b p.2 MCMC q () MCMC (MLE) () i N i y i ) p(y i q) = ( Ni y i q y i (1 q) N i y i, q L(q ) q log L(q ) q 0 ˆq log L(q )/ q = 0 q log likelihood -50-45 -40 ˆq = = 73 160 = 0.456 0.3 0.4 0.5 0.6 * kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 7 / 90 kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 8 / 90 MCMC 8 y i () MCMC ˆq = 0.46 ( ) 8 y 0.46 y 0.54 8 y y i kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 9 / 90 kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 10 / 90 MCMC MCMC : p(q Y) = () p(y q) p(q) p(y) () p(y q): q Y () p(q): q () p(q Y): Y q () p(y): Y (q!) 1 p(y q) 2 p(q) 3 p(q Y ) 1 p(y q) 2 p(q) 3 p(q Y ) p(y q) 4 p(q)? kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 11 / 90 kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 12 / 90

kubo2017sep16b p.3 MCMC MCMC log likelihood -50-45 -40 log L(q) * 0.3 0.4 0.5 0.6 q L(q) p(y q) L(q) q? p(q Y) Y p(y q) L(q) p(y q) q p(q) OK? 0.0 0.2 0.4 0.6 0.8 1.0? : q q? kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 13 / 90 kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 14 / 90 MCMC MCMC MCMC: : GLM MCMC:? kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 15 / 90 kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 16 / 90 MCMC MCMC: MCMC MCMC: : MCMC () Markov chain Monte Carlo (MCMC) Metropolis Method (Metropolis method) :?? MCMC Metropolis Method kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 17 / 90 MCMC log likelihood -50-45 -40 log L(q) * 0.3 0.4 0.5 0.6 : q log likelihood -50-45 -40 0.3 0.4 0.5 0.6 q () kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 18 / 90

kubo2017sep16b p.4 MCMC MCMC: MCMC MCMC: q log likelihood -49-48 -47-46 -45-44 -47.62-46.38-45.24 0.28 0.29 0.30 0.31 0.32 q: seed sruvivorship 1 q 2 q 3 (1), (2) kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 19 / 90 Metropolis Method : 1 q ( q 0.3) 2 q ( q q new ) 3 q new L(q new ) L(q) L(q new ) L(q) (): q q new L(q new ) < L(q) (): r = L(q new)/l(q) q qnew 1 r q 4 2. (q = 0.01 q = 0.99 ) kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 20 / 90 MCMC MCMC: MCMC MCMC: Metropolis Method q log likelihood -49-48 -47-46 -45-44 -47.62-46.38-45.24 log likelihood -46-44 -42 0.28 0.29 0.30 0.31 0.32 (MCMC) 0.30 0.31 0.32 0.33 0.34 0.35 Metropolis Method kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 21 / 90 q Metropolis Method ( MCMC) log likelihood -50-45 -40 0.3 0.4 0.5 0.6? q kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 22 / 90 MCMC MCMC: MCMC MCMC: q q? q q? q MCMC MCMC?? kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 23 / 90 kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 24 / 90

kubo2017sep16b p.5 MCMC MCMC: MCMC MCMC: q MCMC q q? kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 25 / 90 MCMC? log L(q) log likelihood -50-45 -40 * 0.3 0.4 0.5 0.6 L(q) MCMC () kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 26 / 90 MCMC MCMC: MCMC MCMC: MCMC q () MCMC p(q) q 95%! kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 27 / 90 q Metropolis Method p(q Y) L(q) p(q) 0.0 0.2 0.4 0.6 0.8 1.0 q kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 28 / 90 MCMC MCMC MCMC? 2. MCMC Gibbs sampling 1 : : MCMC 2 R package : package : 3 BUGS Gibbs sampler : Gibbs sampler? kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 29 / 90 kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 30 / 90

kubo2017sep16b p.6 MCMC のためのソフトウェア MCMC のためのソフトウェア さまざまな MCMC アルゴリズム Gibbs sampling とは何か? MCMC アルゴリズムのひとつ いろいろな MCMC Metropolis Method: 試行錯誤で値を変化させて いく MCMC Metropolis-Hastings: その改良版 複数のパラメーターの MCMC サンプリングに使う 例: パラメーター β1 と β2 の Gibbs sampling 1 2 Gibbs sampling: 条件つき確率分布を使った MCMC β2 に何か適当な値を与える β2 の値はそのままにして その条件のもとでの β1 の MCMC sampling をする (条件つき事後分布) 3 複数の変数 (パラメーター 状態) を効率よくサンプリング an β1 の値はそのままにして その条件のもとでの β2 の MCMC sampling をする (条件つき事後分布) efficient method for sampling of parameter values 4 2. 3. をくりかえす 教科書の第 9 章の例題で説明 kubo (http://goo.gl/ufq2) 統計モデリング入門 (第 2 部) 2017 09 16 31 / 90 kubo (http://goo.gl/ufq2) MCMC のためのソフトウェア 図解: Gibbs sampling β2 のサンプリング BUGS 言語 (+ っぽいもの) でベイズモデルを 8 10 記述できるソフトウェア 6 8 10 4 5 6-0.02 4 6 4 1.6 1.8 2.0 2.2 2.4 β1 3 7 3 4 5 6 0.02 β2 0.06 7 8 10 4 1.6 1.8 2.0 2.2 2.4 β1 3 4 5 6-0.02 4 6 6 8 10 step 2 7 3 4 5 6 0.02 β2 0.06 7 WinBUGS 歴史を変えて さようなら? OpenBUGS 予算が足りなくて停滞? JAGS お手軽で良い どんな OS でも動く Stan いま一番の注目 今日は紹介しませんが 8 10 リンク集: 6 4 1.6 1.8 2.0 2.2 2.4 β1 3 kubo (http://goo.gl/ufq2) 4 5 6 7-0.02 4 6 8 10 step 3 32 / 90 便利な BUGS 汎用 Gibbs sampler たち (統計モデリング入門の第 9 章) step 1 2017 09 16 MCMC のためのソフトウェア β1 のサンプリング MCMC 統計モデリング入門 (第 2 部) 3 4 5 6 0.02 β2 0.06 http://hosho.ees.hokudai.ac.jp/~kubo/ce/bayesianmcmc.html 7 統計モデリング入門 (第 2 部) えーと BUGS 言語って何? 2017 09 16 33 / 90 MCMC のためのソフトウェア kubo (http://goo.gl/ufq2) 統計モデリング入門 (第 2 部) 2017 09 16 34 / 90 MCMC のためのソフトウェア いろいろな OS で使える JAGS 4.2.0 このベイズモデルを BUGS 言語で記述したい データ Y[i] 種子数８個のうちの生存数 二項分布 dbin(q,8) BUGS 言語コード R core team のひとり Martyn Plummer さんが開発 Just Another Gibbs Sampler C++ で実装されている R がインストールされていることが必要 生存確率 q Linux, Windows, Mac OS X バイナリ版もある 無情報事前分布 開発進行中 R からの使う: library(rjags) 矢印は手順ではなく 依存関係をあらわしている BUGS 言語: ベイズモデルを記述する言語 Spiegelhalter et al. 1995. BUGS: Bayesian Using Gibbs Sampling version 0.50. kubo (http://goo.gl/ufq2) 統計モデリング入門 (第 2 部) 2017 09 16 35 / 90 kubo (http://goo.gl/ufq2) 統計モデリング入門 (第 2 部) 2017 09 16 36 / 90

kubo2017sep16b p.7 MCMC JAGS R モデルの構造 データとパラメーターの初期値 サンプリングの詳細 Input BUGS 言語 JAGS MCMC sampling from posterior distributions 151.5 153.0 3 5 7 事後分布からのランダムサンプル Trace of beta[1] 0 500 1500 2500 Iterations Trace of beta[2] 0 500 1500 2500 Iterations 0.0 0.6 1.2 0.0 0.4 0.8 Output Density of beta[1] 151.5 152.5 153.5 154.5 N = 3000 Bandwidth = 0.0479 Density of beta[2] 3 4 5 6 7 8 N = 3000 Bandwidth = 0.08514 kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 37 / 90 MCMC R JAGS (1 / 3) library(rjags) library(r2winbugs) # to use write.model() model.bugs <- function() { for (i in 1:N.data) { Y[i] ~ dbin(q, 8) # } q ~ dunif(0.0, 1.0) # q } file.model <- "model.bug.txt" write.model(model.bugs, file.model) # # kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 38 / 90 MCMC R JAGS (2 / 3) load("mcmc.rdata") # (data.rdata mcmc.rdata!!) list.data <- list(y = data, N.data = length(data)) inits <- list(q = 0.5) n.burnin <- 1000 n.chain <- 3 n.thin <- 1 n.iter <- n.thin * 1000 model <- jags.model( file = file.model, data = list.data, inits = inits, n.chain = n.chain ) # kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 39 / 90 MCMC R JAGS (3 / 3) # burn-in update(model, n.burnin) # burn in # post.mcmc.list post.mcmc.list <- coda.samples( model = model, variable.names = names(inits), n.iter = n.iter, thin = n.thin ) # kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 40 / 90 MCMC burn in? 0.3 0.4 0.5 0.6 0 100 200 300 400 500 MCMC step kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 41 / 90 MCMC MCMC ˆR = 1.019 MCMC ˆR = 2.520 MCMC MCMC step! kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 42 / 90

kubo2017sep16b p.8 MCMC ˆR MCMC Gibbs sampling plot(post.mcmc.list) gelman.diag(post.mcmc.list) R-hat Gelman-Rubin ˆR var ˆ + (ψ y) = W var ˆ + (ψ y) = n 1 n W + 1 n B W : variance B : variance Gelman et al. 2004. Bayesian Data Analysis. Chapman & Hall/CRC 0.35 0.45 0.55 Trace of q 0 2 4 6 8 10 Density of q 1000 1400 1800 Iterations 0.30 0.40 0.50 0.60 N = 1000 Bandwidth = 0.00838 kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 43 / 90 kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 44 / 90! 3. 100 800 403 0.50 GLM 0 5 10 15 20 25! 0 2 4 6 8 y i? ( 10 ) kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 45 / 90 kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 46 / 90 (overdispersion) 0 5 10 15 20 25 0 2 4 6 8 y i 0.5 : overdispersion : Modeling of individual difference i N i y i ) p(y i q i ) = ( Ni y i q y i i (1 q i) N i y i, q i kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 47 / 90 kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 48 / 90

kubo2017sep16b p.9 Survivorship: logistic regression (GLM) q i = q(z i ) q(z) = 1/{1 + exp( z)} 0.0 0.2 0.4 0.6 0.8 1.0 q(z) -6-4 -2 0 2 4 6 z i = a + r i a: r i : i () z r i > 100 101 (a {r 1, r 2,, r 100 }) /! () kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 49 / 90 kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 50 / 90 {r i } s = 1.0 s = 1.5 s = 3.0 r i ) 1 p(r i s) = ( exp r2 i 2πs 2 2s 2 p(r i s) r i r i r i kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 51 / 90 : r i (A) 0 5 10 15 I III I II I I -6-4 -2 0 2 4 6 0 2 4 6 8 y i (B) p(r i s) s = 0.5 50 {r i} s = 3.0 r i 1 q i = 1+exp( r i) 0 5 10 15 II I II I II I I I III I II I I I I I -6-4 -2 0 2 4 6 2.9 9.9 p(y i q i) 0 2 4 6 8 y i kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 52 / 90 r i r i {r i } 100 r i s = 1.0 s = 1.5 r i p(r i s) = s = 3.0 1 2πs 2 exp ( r2 i 2s 2 ) (A) (!) (B) (C) s kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 53 / 90 kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 54 / 90

kubo2017sep16b p.10 r i! s = 1.0 s = 1.5 r i s = 3.0 p(r i s) = 1 2πs 2 exp ( r2 i 2s 2 kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 55 / 90 ) 全データ 個体個体 3 3 のデータのデータ個体 1 のデータ個体個体 3 3 のデータのデータ個体 2 のデータ {r 1, r 2, r 3,..., r 100 } s local parameter random effects a global parameter fixed effects? kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 56 / 90 {r i } s (B) (C) s = 1.0 a, s {r i } s s = 1.5 s = 3.0 s s (non-informative prior) 0 < s < 10 4 kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 57 / 90 kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 58 / 90 a : Hierarchical and non-informative priors 0.0 0.1 0.2 0.3 0.4 ( 0; 1) ( 0; 100) -10-5 0 5 10 (logit) a a 種子 8 個のうち Y[i] が生存 生存 q[i] 全個体共通の 平均 a r[i] s hyper parameter kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 59 / 90 kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 60 / 90

kubo2017sep16b p.11 JAGS R JAGS BUGS model { for (i in 1:N.data) { Y[i] ~ dbin(q[i], 8) logit(q[i]) <- a + r[i] } a ~ dnorm(0, 1.0E-4) for (i in 1:N.data) { r[i] ~ dnorm(0, tau) } tau <- 1 / (s * s) s ~ dunif(0, 1.0E+4) } 種子 8 個のうち Y[i] が生存 生存 q[i] 全個体共通の 平均 a r[i] s hyper parameter kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 61 / 90 kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 62 / 90 JAGS > source("mcmc.list2bugs.r") # > post.bugs <- mcmc.list2bugs(post.mcmc.list) # bugs 80% 10 interval 5 0 for each 5 10 chain 1 R hat 1.5 2+ a * r[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] [38] [39] [40] [41] [42] [43] [44] [45] [46] [47] [48] [49] [50] [51] [52] [53] [54] [55] [56] [57] [58] [59] [60] [61] [62] [63] [64] [65] [66] [67] [68] [69] [70] [71] [72] [73] [74] [75] s [76] 10 5 0 5 10 1 1.5 2+ * array truncated for lack of space 3 chains, each with 4000 iterations (first 2000 discarded) 0.02 0.01 a 0 0.01 0.02 10 5 * r 0 5 10 3.5 s 3 2.5 medians and 80% intervals 1 2 3 4 5 6 7 8 910 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 63 / 90 bugs post.bugs print(post.bugs, digits.summary = 3) 95% 3 chains, each with 4000 iterations (first 2000 discarded), n.thin = 2 n.sims = 3000 iterations saved mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff a 0.020 0.321-0.618-0.190 0.028 0.236 0.651 1.007 380 s 3.015 0.359 2.406 2.757 2.990 3.235 3.749 1.002 1200 r[1] -3.778 1.713-7.619-4.763-3.524-2.568-1.062 1.001 3000 r[2] -1.147 0.885-2.997-1.700-1.118-0.531 0.464 1.001 3000 r[3] 2.014 1.074 0.203 1.282 1.923 2.648 4.410 1.001 3000 r[4] 3.765 1.722 0.998 2.533 3.558 4.840 7.592 1.001 3000 r[5] -2.108 1.111-4.480-2.775-2.047-1.342-0.164 1.001 2300... () r[99] 2.054 1.103 0.184 1.270 1.996 2.716 4.414 1.001 3000 r[100] -3.828 1.766-7.993-4.829-3.544-2.588-1.082 1.002 1100 kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 64 / 90 R Trace of a Density of a post.mcmc <- to.mcmc(post.bugs) 1.0 0.5 2.0 4.0 0 200 600 1000 Iterations Trace of s 0 200 600 1000 Iterations 0.0 0.8 0.0 0.8 1.0 0.0 0.5 1.0 N = 1000 Bandwidth = 0.06795 Density of s 2.0 3.0 4.0 5.0 N = 1000 Bandwidth = 0.07627 matrix 0 5 10 15 20 25 0 2 4 6 8 kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 65 / 90 kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 66 / 90

kubo2017sep16b p.12 GLMM? 線形モデルの発展推定計算方法階層ベイズモデル (HBM) MCMC もっと自由な統計モデリングを! 一般化線形混合モデル (GLMM) 最尤推定法個体差 場所差といった変量効果をあつかいたい一般化線形モデル (GLM) 正規分布以外の確率分布をあつかいたい 最小二乗法線形モデル (Generalized Linear Mixed Model) 4. + GLMM local parameter GLMM kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 67 / 90 kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 68 / 90 :? pot G pot E pot B pot H pot C pot I pot A pot D pot J pot F y i 10 10 ( 100 ) (f j = C) 5 ( 50 ) (f j = T) 5 ( 50 ) kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 69 / 90 pot A pot G pot J pot H pot E pot D pot C pot I pot F pot B y i 10 10 ( 100 ) (f j = C) 5 ( 50 ) (f j = T) 5 ( 50 ) kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 70 / 90!! I > d <- read.csv("d1.csv") > head(d) id pot f y id : {1, 2, 3,, 100} 1 1 A C 6 pot : {A, B, C, 2 2 A C 3, J} 3 3 A C 19 f : : C, 4 4 A C 5 T 5 5 A C 0 y : () 6 6 A C 19 d$y 0 10 20 30 D D J B AA D G I D E I A G II A AA A A B DE E G C D E H B C A B C B C CC D EE D E EF F F F I G H H I II J J JJ J 0 20 40 60 80 100 d$id plot(d$id, d$y, pch = as.character(d$pot),...)? kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 71 / 90 kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 72 / 90

kubo2017sep16b p.13? I d$y 0 10 20 30 D D J B AA D G I D E I A G II A AA A A B DE E G C D E H B C A B C B C CC D EE D E EF F F F I G H H I II J J JJ J 0 20 40 60 80 100 d$id 0 10 20 30 A B C D E F G H I J? () plot(d$pot, d$y, col = rep(c("blue", "red"), each = 5)) random effects kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 73 / 90 kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 74 / 90 () 階層ベイズモデル もっと自由な統計モデリングを! 線形モデルの発展 (HBM) 一般化線形混合モデル 個体差 場所差といった変量効果をあつかいたい 一般化線形モデル 正規分布以外の確率分布をあつかいたい (GLMM) 推定計算方法 MCMC 最尤推定法 (GLM) 最小二乗法 線形モデル GLM: > summary(glm(y ~ f, data = d, family = poisson))...()... Coefficients: Estimate Std. Error z value Pr(> z ) (Intercept) 1.8931 0.0549 34.49 < 2e-16 ft -0.4115 0.0869-4.73 2.2e-06...()... (f)? AIC kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 75 / 90 kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 76 / 90 GLMM:! > library(glmmml) > summary(glmmml(y ~ f, data = d, family = poisson, + cluster = id))...()... coef se(coef) z Pr(> z ) (Intercept) 1.351 0.192 7.05 1.8e-12 ft -0.737 0.280-2.63 8.4e-03...()...?? kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 77 / 90,! kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 78 / 90

kubo2017sep16b p.14 R! Use JAGS! + R GLMM library(glmmml) glmmml() library(lme4) lmer() library(nlme) nlme() () GLMM + () log log(λ i ) = a + bf i + () + () a f i b () ( σ 1, σ 2 ) σ ([0, 10 4 ] ) kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 79 / 90 kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 80 / 90 データ各個体の種子数 Y[i] 個 ポアソン分布平均 lambda[i] 全個体共通 beta[1] 無情報事前分布 データ ( 説明変数 ) 施肥処理 F[i] 植木鉢差 rp[pot[i]] 個体差 r[i] 階層事前分布階層事前分布 s[1] 個体差の s[2] 植木鉢差の beta[2] ばらつきばらつき 無情報事前分布 ( 超事前分布 ) kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 81 / 90 + BUGS code (1) model { for (i in 1:N.sample) { Y[i] ~ dpois(lambda[i]) log(lambda[i]) <- a + b * F[i] + r[i] + rp[pot[i]] } # BUGS coding f i {C, T} F[i] 0, 1 Pot[i] 1, 2,..., 10 rp[...] kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 82 / 90 + BUGS code (2) # a ~ dnorm(0, 1.0E-4) # b ~ dnorm(0, 1.0E-4) # for (i in 1:N.sample) { r[i] ~ dnorm(0, tau[1]) # } for (j in 1:N.pot) { rp[j] ~ dnorm(0, tau[2]) # () } for (k in 1:N.tau) { tau[k] <- 1.0 / (sigma[k] * sigma[k]) # sigma[k] ~ dunif(0, 1.0E+4) } } kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 83 / 90 ( b)? mean sd 2.5% 25% 50% 75% 97.5% Rhat n. a 1.501 0.529 0.482 1.157 1.493 1.852 2.565 1.032 b -1.016 0.706-2.436-1.476-0.993-0.565 0.395 1.019 sigma[1] 1.020 0.114 0.822 0.939 1.014 1.089 1.265 1.004...()... kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 84 / 90

kubo2017sep16b p.15 () rp[j]! random effects random effects fixed effects! kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 85 / 90 kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 86 / 90 : 1. 2. (GLM) 3. MCMC 4. GLM kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 87 / 90 kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 88 / 90 kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 89 / 90 階層ベイズモデル もっと自由な統計モデリングを! 線形モデルの発展 (HBM) 一般化線形混合モデル 個体差 場所差といった変量効果をあつかいたい 一般化線形モデル 正規分布以外の確率分布をあつかいたい (GLMM) 推定計算方法 MCMC 最尤推定法 (GLM) 最小二乗法 線形モデル () kubo (http://goo.gl/ufq2) ( 2 ) 2017 09 16 90 / 90