kubostat2017j p.1 2017 (j) Categorical Data Analsis kubo@ees.hokudai.ac.jp http://goo.gl/76c4i 2017 11 15 : 2017 11 08 17:11 kubostat2017j (http://goo.gl/76c4i) 2017 (j) 2017 11 15 1 / 63 A B C D E F G H I 0 62 21 14 11 10 2 0 2 1 48 34 22 17 16 7 2 1 1!!!! 5?!!!!! 2012 03 19 (2013 03 02 11 :27 ) 5/ 51 kubostat2017j (http://goo.gl/76c4i) 2017 (j) 2017 11 15 2 / 63 GLM ( : GLM ) ( : MCMC ) ( : Neman-Pearson ) 2 2 2012 03 19 (2013 03 02 11 :27 ) 6/ 51 kubostat2017j (http://goo.gl/76c4i) 2017 (j) 2017 11 15 3 / 63 kubostat2017j 2012 03 19 (http://goo.gl/76c4i) (2013 03 02 11 :27 ) 2017 (j) 2017 11 15 7/ 51 4 / 63 : (?) ( = 0) ( = 1) () A B { A,0, B,0, A,1, B,1} : A B?? R data.frame() tabs() 2012 03 19 (2013 03 02 11 :27 ) 8/ 51 kubostat2017j (http://goo.gl/76c4i) 2017 (j) 2017 11 15 5 / 63 kubostat2017j 2012 03 19 (http://goo.gl/76c4i) (2013 03 02 11 :27 ) 2017 (j) 2017 11 15 9/ 51 6 / 63
kubostat2017j p.2 CSV CSV (!) d2.csv d2.csv,, 286,0,A 85,0,B 378,1,A 148,1,B 2012 03 19 (2013 03 02 11 :27 ) 10/ 51 kubostat2017j (http://goo.gl/76c4i) 2017 (j) 2017 11 15 7 / 63 R data.frame > d2 <- read.csv("d2.csv") # data.frame! > d2 # d2 data.frame 1 286 0 A 2 85 0 B 4 378 1 A 5 148 1 B 2012 03 19 (2013 03 02 11 :27 ) 11/ 51 kubostat2017j (http://goo.gl/76c4i) 2017 (j) 2017 11 15 8 / 63 tabs: R 1 286 0 A 2 85 0 B 4 378 1 A 5 148 1 B > (ct2 <- tabs( ~ +, data = d2)) A B 0 286 85 1 378 148 2012 03 19 (2013 03 02 11 :27 ) 12/ 51 kubostat2017j (http://goo.gl/76c4i) 2017 (j) 2017 11 15 9 / 63 tabs: > tabs( ~, data = d2) 371 526 > tabs( ~, data = d2) A B 664 233 > tabs( ~ +, data = d2) A 286 378 B 85 148 kubostat2017j 2012 03 19 (http://goo.gl/76c4i) (2013 03 02 11 :27 ) 2017 (j) 2017 11 15 13/ 51 10 / 63 tabs: A B 0 286 85 1 378 148 > plot(ct2, col = c("orange", "blue")) ct2 librar(lattice) A B 0 286 85 1 378 148 > librar(lattice) > plot(log() ~ factor(), data = d2, groups =, tpe = "b") 6.0 A 5.5 log() 5.0 B 4.5 2012 03 19 (2013 03 02 11 :27 ) 14/ 51 kubostat2017j (http://goo.gl/76c4i) 2017 (j) 2017 11 15 11 / 63 factor() 2012 03 19 (2013 03 02 11 :27 ) 15/ 51 kubostat2017j (http://goo.gl/76c4i) 2017 (j) 2017 11 15 12 / 63
B A kubostat2017j p.3 GLM 2 2 GLM logistic regression A, Binom(q A,, A, + B,) logit(q A,) = a A + b A A B 0 286 85 1 378 148 > summar(glm(ct2 ~ c(0, 1), data = d2, famil = binomial)) Estimate Std. Error z value Pr(> z ) (Intercept) 1.213 0.124 9.82 <2e-16 c(0, 1) -0.276 0.157-1.76 0.079 kubostat2017j 2012 03 19 (http://goo.gl/76c4i) (2013 03 02 11 :27 ) 2017 (j) 2017 11 15 16/ 51 13 / 63 2012 03 19 (2013 03 02 11 :27 ) 17/ 51 kubostat2017j (http://goo.gl/76c4i) 2017 (j) 2017 11 15 14 / 63 logit(q A,) = 1.213 + ( 0.276) q_a ct2 plot(ct2) AIC logit(q A,) = a A + b B 16.5 logit(q A,) = a A 17.6 2 2 GLM Poisson regression log-linear model 2012 03 19 (2013 03 02 11 :27 ) 18/ 51 kubostat2017j (http://goo.gl/76c4i) 2017 (j) 2017 11 15 15 / 63 kubostat2017j 2012 03 19 (http://goo.gl/76c4i) (2013 03 02 11 :27 ) 2017 (j) 2017 11 15 19/ 51 16 / 63 GLM ( ) A GLM ( ) B A, Pois(λ A,) log(λ A,) = α A + β A B, Pois(λ B,) log(λ B,) = α B + β B > # A > summar(glm( ~, data = d2[d2$ == "A",], famil = poisson)) Estimate Std. Error z value Pr(> z ) (Intercept) 5.6560 0.0591 95.65 < 2e-16 0.2789 0.0784 3.56 0.00037 > # B > summar(glm( ~, data = d2[d2$ == "B",], famil = poisson)) Estimate Std. Error z value Pr(> z ) (Intercept) 4.443 0.108 40.96 < 2e-16 0.555 0.136 4.07 4.6e-05 2012 03 19 (2013 03 02 11 :27 ) 20/ 51 kubostat2017j (http://goo.gl/76c4i) 2017 (j) 2017 11 15 17 / 63 2012 03 19 (2013 03 02 11 :27 ) 21/ 51 kubostat2017j (http://goo.gl/76c4i) 2017 (j) 2017 11 15 18 / 63
kubostat2017j p.4 log(λ A,) = 5.66 + 0.279 log(λ B,) = 4.44 + 0.555 lambda_a 00 200 300 400 500 lambda_b 00 200 300 400 500 AIC AIC λ A, = α A + β A 19.3 λ B, = α B + β B 17.1 λ A, = α A 30.1 λ B, = α B 32.4 2012 03 19 (2013 03 02 11 :27 ) 22/ 51 kubostat2017j (http://goo.gl/76c4i) 2017 (j) 2017 11 15 19 / 63 GLM GLM GLM: logit(q A,) = a A + b A 1 q A, = 1 + ep[ (a A + b A)] : log(λ A,) = α A + β A λ A, = ep(α A + β A) λ B, = ep(α B + β B) A? λ A, ep(α A + β A ) = λ A, + λ B, ep(α A + β A ) + ep(α B + β B ) 1 = 1 + ep[α B α A + (β B β A )] 2012 03 19 (2013 03 02 11 :27 ) 23/ 51 kubostat2017j (http://goo.gl/76c4i) 2017 (j) 2017 11 15 20 / 63 : GLM GLM GLM q A, = GLM ( ) 1 1 + ep[ (a A + b A)] λ A, 1 = λ A, + λ B, 1 + ep[α B α A + (β B β A)] GLM GLM a A = α A α B b A = β A β B 2012 03 19 (2013 03 02 11 :27 ) 24/ 51 kubostat2017j (http://goo.gl/76c4i) 2017 (j) 2017 11 15 21 / 63 : GLM GLM GLM a A = 1.213 = α A α B b A = -0.276 = β A β B > GLM (A ) > glm(ct2 ~ c(0, 1), data = d2, famil = binomial) (Intercept) c(0, 1) 1.213-0.276 > GLM (A ) > glm( ~, data = d2[d2$ == "A",], famil = poisson) (Intercept) 5.656 0.279 > GLM (B ) > glm( ~, data = d2[d2$ == "B",], famil = poisson) (Intercept) 4.443 0.555 kubostat2017j 2012 03 19 (http://goo.gl/76c4i) (2013 03 02 11 :27 ) 2017 (j) 2017 11 15 25/ 51 22 / 63 : GLM GLM GLM (A ) GLM (B ) lambda_a 00 200 300 400 500 lambda_* 00 200 300 400 500 + lambda_b 00 200 300 400 500 GLM (A + B ) kubostat2017j 2012 03 19 (http://goo.gl/76c4i) (2013 03 02 11 :27 ) 2017 (j) 2017 11 15 26/ 51 23 / 63 q_a 2 2 GLM? kubostat2017j 2012 03 19 (http://goo.gl/76c4i) (2013 03 02 11 :27 ) 2017 (j) 2017 11 15 27/ 51 24 / 63
kubostat2017j p.5 GLM ( ) > summar(glm( ~ *, data = d2, famil = poisson)) Coefficients: Estimate Std. Error z value Pr(> z ) (Intercept) 5.6560 0.0591 95.65 < 2e-16 0.2789 0.0784 3.56 0.00037 B -1.2133 0.1235-9.82 < 2e-16 :B 0.2757 0.157.76 0.07921 GLM GLM α A = 5.66 α B = 5.66 1.21 β A = 0.279 β B = 0.279 + 0.276 2012 03 19 (2013 03 02 11 :27 ) 28/ 51 kubostat2017j (http://goo.gl/76c4i) 2017 (j) 2017 11 15 25 / 63 GLM ( ) GLM (A ) lambda_a 00 200 300 400 500 lambda_* 00 200 300 400 500 + GLM (B ) lambda_b 00 200 300 400 500 GLM (A + B ) kubostat2017j 2012 03 19 (http://goo.gl/76c4i) (2013 03 02 11 :27 ) 2017 (j) 2017 11 15 29/ 51 26 / 63 q_a tabs: 2 3 GLM? GLM? A B C 0 286 85 7 1 378 148 17 > plot(ct3, col = c("orange", "blue", "green")) ct3 A C B kubostat2017j 2012 03 19 (http://goo.gl/76c4i) (2013 03 02 11 :27 ) 2017 (j) 2017 11 15 30/ 51 27 / 63 2012 03 19 (2013 03 02 11 :27 ) 31/ 51 kubostat2017j (http://goo.gl/76c4i) 2017 (j) 2017 11 15 28 / 63 librar(lattice) A B C 0 286 85 7 1 378 148 17 > librar(lattice) > plot(log() ~ factor(), data = d3, groups =, tpe = "b") 6 GLM ( ) > glm( ~ *, data = d3, famil = poisson) Coefficients: (Intercept) B C :B :C 5.656 0.279-1.213-3.710 0.276 0.608 GLM log() 5 4 3 2 A, Pois(λ A,) log(λ A,) = α A + β A α A = 5.66 α B = 5.66 1.21 α C = 5.66 3.71 β A = 0.279 β B = 0.279 + 0.276 β C = 0.279 + 0.608 factor() 2012 03 19 (2013 03 02 11 :27 ) 32/ 51 kubostat2017j (http://goo.gl/76c4i) 2017 (j) 2017 11 15 29 / 63 2012 03 19 (2013 03 02 11 :27 ) 33/ 51 kubostat2017j (http://goo.gl/76c4i) 2017 (j) 2017 11 15 30 / 63
kubostat2017j p.6 GLM GLM (A ) (B ) (C ) lambda_a 00 200 300 400 500 + + lambda_* 00 200 300 400 500 lambda_b 00 200 300 400 500 lambda_b 00 200 300 400 500 GLM kubostat2017j 2012 03 19 (http://goo.gl/76c4i) (2013 03 02 11 :27 ) 2017 (j) 2017 11 15 34/ 51 31 / 63 lambda_* GLM GLM librar(nnet) multinom() kubostat2017j 2012 03 19 (http://goo.gl/76c4i) (2013 03 02 11 :27 ) 2017 (j) 2017 11 15 35/ 51 32 / 63 GLM > ct3 # A B C 0 286 85 7 1 378 148 17 > librar(nnet) # nnet package > multinom(ct3 ~ c(0, 1)) GLM Coefficients: B, Multinom(q B,, 3 ) (Intercept) c(0, 1) C, Multinom(q C,, 3 ) B -1.2133 0.27552 logit(q B,) = a B + b B C -3.7097 0.60763 logit(q C,) = a C + b C > # GLM! 2012 03 19 (2013 03 02 11 :27 ) 36/ 51 kubostat2017j (http://goo.gl/76c4i) 2017 (j) 2017 11 15 33 / 63? cf. 2004 2012 03 19 (2013 03 02 11 :27 ) 37/ 51 kubostat2017j (http://goo.gl/76c4i) 2017 (j) 2017 11 15 34 / 63 2 9 GLM? kubostat2017j 2012 03 19 (http://goo.gl/76c4i) (2013 03 02 11 :27 ) 2017 (j) 2017 11 15 38/ 51 35 / 63 : 3 9! > d9 1 62 0 A 2 21 0 B 3 14 0 C 4 11 0 D 5 10 0 E 6 10 0 F 7 2 0 G 8 0 0 H 9 2 0 I 10 48 1 A 15 7 1 F 16 2 1 G 17 1 1 H 18 1 1 I > ct9 A B C D E F G H I 0 62 21 14 11 10 2 0 2 1 48 34 22 17 16 7 2 1 1 2012 03 19 (2013 03 02 11 :27 ) 39/ 51 kubostat2017j (http://goo.gl/76c4i) 2017 (j) 2017 11 15 36 / 63
kubostat2017j p.7 tabs: A B C D E F G H I 0 62 21 14 11 10 2 0 2 1 48 34 22 17 16 7 2 1 1 > plot(ct9, col = c( )) cts librar(lattice) A B C D E F G H I 0 62 21 14 11 10 2 0 2 1 48 34 22 17 16 7 2 1 1 > librar(lattice) > plot(sqrt() ~ factor(), data = d9, groups =, tpe = "b") 8 E B D C A 6 sqrt() 4 2 F 0 IHG 2012 03 19 (2013 03 02 11 :27 ) 40/ 51 kubostat2017j (http://goo.gl/76c4i) 2017 (j) 2017 11 15 37 / 63 factor() 2012 03 19 (2013 03 02 11 :27 ) 41/ 51 kubostat2017j (http://goo.gl/76c4i) 2017 (j) 2017 11 15 38 / 63 GLM ( ) > ct9 A B C D E F G H I 0 62 21 14 11 10 2 0 2 1 48 34 22 17 16 7 2 1 1 > summar(glm( ~ *, data = d9, famil = poisson)) (Intercept) B C 4.127-0.256-1.083-1.488 D E F G -1.729-1.825-1.825-3.434 H I :B :C -26.430-3.434 0.738 0.708 :D :E :F :G 0.691 0.726-0.101 0.256 :H :I 22.559-0.437 H! kubostat2017j 2012 03 19 (http://goo.gl/76c4i) (2013 03 02 11 :27 ) 2017 (j) 2017 11 15 42/ 51 39 / 63 glm() 2012 03 19 (2013 03 02 11 :27 ) 43/ 51 kubostat2017j (http://goo.gl/76c4i) 2017 (j) 2017 11 15 40 / 63 GLMM? ( ) > (fit.glmm <- glmmml( ~, data = d9, cluster =, famil = poisson)) coef se(coef) z Pr(> z ) (Intercept) 2.012 0.465 4.323 1.5e-05 0.114 0.120 0.956 3.4e-01 > fit.glmm$posterior.modes [1] 1.926862 1.230935 0.807098 0.557225 0.483854 [6] 0.067305-1.218522-2.005846-1.426968 データ ( 応答変数 ) 各個体の種子数 Y[i] 個 ポアソン分布平均 lambda[i] 全個体共通 a0 無情報事前分布 b0 データ ( 説明変数 ) X[i] 傾きの差 b[[i]] 切片の差 a[[i]] 階層事前分布階層事前分布 s[1] 切片差の s[2] 傾き差のばらつきばらつき 無情報事前分布 ( 超事前分布 ) 2012 03 19 (2013 03 02 11 :27 ) 44/ 51 kubostat2017j (http://goo.gl/76c4i) 2017 (j) 2017 11 15 41 / 63 WinBUGS (MCMC ) BUGS WinBUGS R 2012 03 19 (2013 03 02 11 :27 ) 45/ 51 kubostat2017j (http://goo.gl/76c4i) 2017 (j) 2017 11 15 42 / 63
kubostat2017j p.8 ( WinBUGS ) 1 2012 03 19 (2013 03 02 11 :27 ) 46/ 51 kubostat2017j (http://goo.gl/76c4i) 2017 (j) 2017 11 15 43 / 63 2010 03 15 (2010-03-23 16:29 ) 11/ 33 kubostat2017j (http://goo.gl/76c4i) 2017 (j) 2017 11 15 44 / 63 1: 0 5 15 0 5 15? (isometric)?? 2010 03 15 (2010-03-23 16:29 ) 12/ 33 kubostat2017j (http://goo.gl/76c4i) 2017 (j) 2017 11 15 45 / 63 0 5 15 0 5 15 above ground mass log(y) 0.5 1..5 2.0 2.5 0.5 1..5 2.0 2.5 below ground mass log(x) 2010 03 15 (2010-03-23 16:29 ) 13/ 33 kubostat2017j (http://goo.gl/76c4i) 2017 (j) 2017 11 15 46 / 63 ( ): 観測できない世界観測できる世界 parameter process observation 無情報事前分布 w q 1- q Y X 地上部 地下部 http://tech-jam.com 2010 03 15 (2010-03-23 16:29 ) 14/ 33 kubostat2017j (http://goo.gl/76c4i) 2017 (j) 2017 11 15 47 / 63 ( ): 観測できない世界観測できる世界 parameter process observation 無情報事前分布 階層的事前分布 w 1- q 個体差 q Y X 地上部 地下部 測定時の誤差 http://tech-jam.com 2010 03 15 (2010-03-23 16:29 ) 15/ 33 kubostat2017j (http://goo.gl/76c4i) 2017 (j) 2017 11 15 48 / 63
kubostat2017j p.9 BUGS code (process ) for (i in 1:N) { Y[i] ~ dnorm([i], Tau.err) # X[i] ~ dnorm([i], Tau.err) # [i] <- q[i] * w[i] [i] <- (1 - q[i]) * w[i] logit(q[i]) <- a + re[i] w[i] <- ep(log.w[i]) log.w[i] ~ dnorm(0, Tau.noninformative) #!! }# log.w[i] +! 2010 03 15 (2010-03-23 16:29 ) 16/ 33 kubostat2017j (http://goo.gl/76c4i) 2017 (j) 2017 11 15 49 / 63 : MCMC 1. BUGS code (model1.tt) 2. R (runbus1.r) 3. R runbus1.r (source("runbugs1.r")) 4. R librar(r2winbugs) WinBUGS 5. WinBUGS Markov chain Monte Carlo (MCMC) 6. R 2010 03 15 (2010-03-23 16:29 ) 17/ 33 kubostat2017j (http://goo.gl/76c4i) 2017 (j) 2017 11 15 50 / 63 0 5 15 0 5 15 above ground mass log(y) 0.5 1..5 2.0 2.5 0 5 15 (median) (95% CI) 2010 03 15 (2010-03-23 16:29 ) 18/ 33 kubostat2017j (http://goo.gl/76c4i) 2017 (j) 2017 11 15 51 / 63 0 5 15 0.5 1..5 2.0 2.5 below ground mass log(x) 2010 03 15 (2010-03-23 16:29 ) 19/ 33 kubostat2017j (http://goo.gl/76c4i) 2017 (j) 2017 11 15 52 / 63??? ( ) 2? 2010 03 15 (2010-03-23 16:29 ) 20/ 33 kubostat2017j (http://goo.gl/76c4i) 2017 (j) 2017 11 15 53 / 63 2010 03 15 (2010-03-23 16:29 ) 21/ 33 kubostat2017j (http://goo.gl/76c4i) 2017 (j) 2017 11 15 54 / 63
kubostat2017j p.10 2:? 0 5 15 0 5 15 above ground mass log(y) 0.5 1..5 2.0 2.5 0 5 15 2010 03 15 (2010-03-23 16:29 ) 22/ 33 kubostat2017j (http://goo.gl/76c4i) 2017 (j) 2017 11 15 55 / 63 0 5 15? -2-1 2 3 below ground mass log(x) 2010 03 15 (2010-03-23 16:29 ) 23/ 33 kubostat2017j (http://goo.gl/76c4i) 2017 (j) 2017 11 15 56 / 63 : q w 観測できない世界観測できる世界 parameter process observation 無情報事前分布 階層的事前分布 w q 1- q 個体差 Y X 地上部 地下部 測定時の誤差 http://tech-jam.com 2010 03 15 (2010-03-23 16:29 ) 24/ 33 kubostat2017j (http://goo.gl/76c4i) 2017 (j) 2017 11 15 57 / 63 BUGS code ( ) + ( ) logit(q[i]) <- a + re[i] (w) (b dnorm(0, Tau.noninformative) ) logit(q[i]) <- ( a + b * (log.w[i] - Mean.log.w) + re[i] ) # Mean.log.w WinBUGS R WinBUGS MCMC 2010 03 15 (2010-03-23 16:29 ) 25/ 33 kubostat2017j (http://goo.gl/76c4i) 2017 (j) 2017 11 15 58 / 63 : 0 5 15 w b MCMC 2010 03 15 (2010-03-23 16:29 ) 26/ 33 kubostat2017j (http://goo.gl/76c4i) 2017 (j) 2017 11 15 59 / 63 0 5 15 (median) (95% CI) 2010 03 15 (2010-03-23 16:29 ) 27/ 33 kubostat2017j (http://goo.gl/76c4i) 2017 (j) 2017 11 15 60 / 63
kubostat2017j p.11 (!) 0 5 15 above ground mass log(y) 0.5 1..5 2.0 2.5 0 5 15-2 -1 2 3 below ground mass log(x) ( ) 2010 03 15 (2010-03-23 16:29 ) 28/ 33 kubostat2017j (http://goo.gl/76c4i) 2017 (j) 2017 11 15 61 / 63 2010 03 15 (2010-03-23 16:29 ) 29/ 33 kubostat2017j (http://goo.gl/76c4i) 2017 (j) 2017 11 15 62 / 63 random effects 3 1: ( )? 2: ( logistic ) 2010 03 15 (2010-03-23 16:29 ) 30/ 33 kubostat2017j (http://goo.gl/76c4i) 2017 (j) 2017 11 15 63 / 63