60 (W30)? 1. ( ) 2. ( ) web site URL ( :41 ) 1/ 77

60 (W30)? 1. ( ) kubo@ees.hokudai.ac.jp 2. ( ) web site URL http://goo.gl/e1cja!! 2013 03 07 (2013 03 07 17 :41 ) 1/ 77

! : :? 2013 03 07 (2013 03 07 17 :41 ) 2/ 77

2013 03 07 (2013 03 07 17 :41 ) 3/ 77!!

2013 03 07 (2013 03 07 17 :41 ) 4/ 77 ( )!

2013 03 07 (2013 03 07 17 :41 ) 5/ 77 : offset : :???

2013 03 07 (2013 03 07 17 :41 ) 6/ 77 GLM

2013 03 07 (2013 03 07 17 :41 ) 7/ 77 : R http://www.r-project.org/ OS free software S http://goo.gl/ufq2 R http://goo.gl/cki2h

2013 03 07 (2013 03 07 17 :41 ) 8/ 77

2013 03 07 (2013 03 07 17 :41 ) 9/ 77

2013 03 07 (2013 03 07 17 :41 ) 10/ 77

2013 03 07 (2013 03 07 17 :41 ) 11/ 77

2013 03 07 (2013 03 07 17 :41 ) 12/ 77

2013 03 07 (2013 03 07 17 :41 ) 13/ 77

2013 03 07 (2013 03 07 17 :41 ) 14/ 77

: : ( ) 2013 03 07 (2013 03 07 17 :41 ) 15/ 77

2013 03 07 (2013 03 07 17 :41 ) 16/ 77 ( http://goo.gl/76c4i)? : offset

何でも割算! 1. 2. 3. 4. 1, 2, 3! x /x 2013 03 07 (2013 03 07 17 :41 ) 17/ 77

? /? : 10 3 200 60 3? ( ) 2013 03 07 (2013 03 07 17 :41 ) 18/ 77

2013 03 07 (2013 03 07 17 :41 ) 19/ 77 : specific leaf area (SLA) : offset : N k :

2013 03 07 (2013 03 07 17 :41 ) 20/ 77 (1) offset

2013 03 07 (2013 03 07 17 :41 ) 21/ 77 :? x {0.1, 0.2,, 1.0} 10 glm(..., family = poisson)

?!! x A = /! glm() offset 2013 03 07 (2013 03 07 17 :41 ) 22/ 77

2013 03 07 (2013 03 07 17 :41 ) 23/ 77 vs plot(d$x, d$y / d$area) d$y/d$area 0 5 10 15 0.2 0.4 0.6 0.8 1.0 d$x

GLM! family: poisson, link : "log" : y ~ x offset : log(area) z = a + b x + log(area) a, b λ log(λ) = z λ = exp(z) = exp(a + b x + log(area)) λ : 2013 03 07 (2013 03 07 17 :41 ) 24/ 77

2013 03 07 (2013 03 07 17 :41 ) 25/ 77 (2)

2013 03 07 (2013 03 07 17 :41 ) 26/ 77 (1 = ) 1.0 0.8 0.6 0.4 0.2 0.0 0 2 4 6 8 10 seed size [ ] ( )

2013 03 07 (2013 03 07 17 :41 ) 27/ 77 ( ) 1.0 0.8 0.6 0.4 0.2 1.0 0.5 0.0 0.0 0 2 4 6 8 10 seed size 0 2 4 6 8 10 seed size 1. ( 4 ) 2. ; {0, 1} 3. ( or or & )

2013 03 07 (2013 03 07 17 :41 ) 28/ 77? 1 / 2 100 / 200! 1.0 0.5 0.0 0 2 4 6 8 10 seed size? 1? (? )

2013 03 07 (2013 03 07 17 :41 ) 29/ 77 R glm() : 1.0 0.8 0.6 0.4 0.2 0.0 0 2 4 6 8 10 seed size (x) or q q = 1 1 + exp( (a + bx)) (logistic ) a b R glm() ( )

2013 03 07 (2013 03 07 17 :41 ) 30/ 77 (binomial distribution)? y i {0, 1, 2,, N} (paramter: q, N) ( ) N q y (1 q) N y y Nq Nq(1 q) probability 0.30 0.25 0.20 0.15 0.10 0.05 0.00 0.30 prob = 0.2 prob = 0.5 prob = 0.8 0.30 0.25 0.20 0.15 0.10 0.05 0.00 : N y 0.25 0.20 0.15 0.10 0.05 0.00 0 2 4 6 8 10 y

2013 03 07 (2013 03 07 17 :41 ) 31/ 77 R glm() :? ( z): x link : logit family: binomial,

2013 03 07 (2013 03 07 17 :41 ) 32/ 77 ( ) Ending 1.0 germination prob. 0.8 0.6 0.4 0.2 0.0 0 2 4 6 8 10 seed size!!? :!

2013 03 07 (2013 03 07 17 :41 ) 33/ 77 (ESJ59 http://goo.gl/qq1ok) : glm()

: OK http://en.wikipedia.org/wiki/poisson_distribution 2013 03 07 (2013 03 07 17 :41 ) 34/ 77

2013 03 07 (2013 03 07 17 :41 ) 35/ 77 (Poisson distribution)? lambda = 1.4 y {0, 1, 2,, } (paramter: λ) 0.3 0.2 0.1 λ y exp( λ) 0.0 lambda = 3.2 0.3 y! probability 0.2 0.1 λ λ lambda = 6.5 0.0 0.3 : 0.2 0.1 0.0 0 2 4 6 8 10 12 y

2013 03 07 (2013 03 07 17 :41 ) 36/ 77 2 2 glm() log-linear model

: > d2 <- read.csv("d2.csv") # > (ct2 <- xtabs(y ~ x + Spc, data = d2)) Spc x A B 0 286 85 1 378 148 > plot(ct2, col = c("orange", "blue")) ct2 0 1 B Spc A 2013 03 07 (2013 03 07 17 :41 ) 37/ 77 x

GLM ( ) A ( ) y A,x Pois(λ A,x) log(λ A,x) = α A + β A x > summary(glm(y ~ x, data = d2[d2$spc == "A",], family = poisson) > # SpcA (......) Estimate Std. Error z value Pr(> z ) (Intercept) 5.6560 0.0591 95.65 < 2e-16 x 0.2789 0.0784 3.56 0.00037 2013 03 07 (2013 03 07 17 :41 ) 38/ 77

2013 03 07 (2013 03 07 17 :41 ) 39/ 77 GLM ( ) B y B,x Pois(λ B,x) log(λ B,x) = α B + β B x > summary(glm(y ~ x, data = d2[d2$spc == "B",], family = poisson) > # SpcB (......) Estimate Std. Error z value Pr(> z ) (Intercept) 4.443 0.108 40.96 < 2e-16 x 0.555 0.136 4.07 4.6e-05?

2013 03 07 (2013 03 07 17 :41 ) 40/ 77 log(λ A,x ) = 5.66 + 0.279x log(λ B,x ) = 4.44 + 0.555x lambda_a 0 100 200 300 400 500 lambda_b 0 100 200 300 400 500 0 1 0 1 x x AIC AIC λ A,x = α A + β A x 19.3 λ B,x = α B + β B x 17.1 λ A,x = α A 30.1 λ B,x = α B 32.4!

2013 03 07 (2013 03 07 17 :41 ) 41/ 77 GLM GLM GLM: logit(q A,x) = a A + b A x 1 q A,x = 1 + exp[ (a A + b A x)] : log(λ A,x) = α A + β A x λ A,x = exp(α A + β A x) λ B,x = exp(α B + β B x) A? λ A,x λ A,x + λ B,x = = exp(α A + β A x) exp(α A + β A x) + exp(α B + β B x) 1 1 + exp[α B α A + (β B β A )x]

: GLM GLM GLM q A,x = 1 1 + exp[ (a A + b A x)] GLM ( ) λ A,x λ A,x + λ B,x = 1 1 + exp[α B α A + (β B β A )x] GLM GLM a A = α A α B b A = β A β B 2013 03 07 (2013 03 07 17 :41 ) 42/ 77

: GLM GLM GLM a A = 1.213 = α A α B b A = -0.276 = β A β B > GLM (A ) > glm(ct2 ~ c(0, 1), data = d2, family = binomial) (Intercept) c(0, 1) 1.213-0.276 > GLM (A ) > glm(y ~ x, data = d2[d2$spc == "A",], family = poisson) (Intercept) x 5.656 0.279 > GLM (B ) > glm(y ~ x, data = d2[d2$spc == "B",], family = poisson) (Intercept) x 4.443 0.555 2013 03 07 (2013 03 07 17 :41 ) 43/ 77

: GLM GLM GLM (A ) GLM (B ) lambda_a 0 100 200 300 400 500 0 1 x lambda_* 0 100 200 300 400 500 + = lambda_b 0 100 200 300 400 500 0 1 x 0 1 x = q_a 0 1 GLM (A + B ) 0 1 x 2013 03 07 (2013 03 07 17 :41 ) 44/ 77

2013 03 07 (2013 03 07 17 :41 ) 45/ 77 2 3 glm() GLM R package

: 2 3 Spc x A B C 0 286 85 7 1 378 148 17 > plot(ct3, col = c("orange", "blue", "green")) ct3 0 1 C B Spc A 2013 03 07 (2013 03 07 17 :41 ) 46/ 77 x

2013 03 07 (2013 03 07 17 :41 ) 47/ 77 GLM ( ) > glm(y ~ x * Spc, data = d3, family = poisson) (......) Coefficients: (Intercept) x SpcB SpcC x:spcb x:spcc 5.656 0.279-1.213-3.710 0.276 0.608 GLM y A,x Pois(λ A,x ) log(λ A,x ) = α A + β A x α 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

GLM GLM lambda_a 0 100 200 300 400 500 (A ) (B ) (C ) 0 1 x + + lambda_* 0 100 200 300 400 500 lambda_b 0 100 200 300 400 500 0 1 x = 0 1 x lambda_b 0 100 200 300 400 500 = 0 1 x lambda_* 0 1 GLM 0 1 x 2013 03 07 (2013 03 07 17 :41 ) 48/ 77

2013 03 07 (2013 03 07 17 :41 ) 49/ 77 2 9 GLM GLMM!

: 3 9! Spc x A B C D E F G H I 0 62 21 14 11 10 10 2 0 2 1 48 34 22 17 16 7 2 1 1 > plot(ct9, col = c( )) cts 0 1 IHG F E D C Spc B A 2013 03 07 (2013 03 07 17 :41 ) 50/ 77 x

GLM ( ) > ct9 Spc x A B C D E F G H I 0 62 21 14 11 10 10 2 0 2 1 48 34 22 17 16 7 2 1 1 > summary(glm(y ~ x * Spc, data = d9, family = poisson)) (Intercept) x SpcB SpcC 4.127-0.256-1.083-1.488 SpcD SpcE SpcF SpcG -1.729-1.825-1.825-3.434 SpcH SpcI x:spcb x:spcc -26.430-3.434 0.738 0.708 x:spcd x:spce x:spcf x:spcg 0.691 0.726-0.101 0.256 x:spch x:spci 22.559-0.437 H! 2013 03 07 (2013 03 07 17 :41 ) 51/ 77

2013 03 07 (2013 03 07 17 :41 ) 52/ 77 glm() GLMM!??! Spc?: 9 Spc ( )

2013 03 07 (2013 03 07 17 :41 ) 53/ 77 σ best GLMM > (fit.glmm <- glmmml(y ~ x, data = d9, cluster = Spc, family = poisson)) ( ) > fit.glmm$posterior.modes [1] 1.926862 1.230935 0.807098 0.557225 0.483854 [6] 0.067305-1.218522-2.005846-1.426968 individual 1 2 3 small posterior prior σ large σ

( ) データ ( 応答変数 ) 各個体の種子数 Y[i] 個 ポアソン分布平均 lambda[i] 全個体共通 a0 無情報事前分布 b0 データ ( 説明変数 ) X[i] 傾きの差 b[spc[i]] 切片の差 a[spc[i]] 階層事前分布 s[1] 切片差のばらつき 階層事前分布 s[2] 傾き差のばらつき 無情報事前分布 ( 超事前分布 ) WinBUGS (MCMC ) BUGS WinBUGS R 2013 03 07 (2013 03 07 17 :41 ) 54/ 77

2013 03 07 (2013 03 07 17 :41 ) 55/ 77 ( )??

2013 03 07 (2013 03 07 17 :41 ) 56/ 77 : λ = 30 λ = 20 30 30+20 = 0.60 5! p = {0.60, 0.40}

2013 03 07 (2013 03 07 17 :41 ) 57/ 77 : λ = 30 λ = 20 λ = 10 30 30+20+10 = 0.50 20 30+20+10 = 0.33 6! p = {0.50, 0.33, 0.17}

2013 03 07 (2013 03 07 17 :41 ) 58/ 77 (Gamma distribution)? y 0 (paramter: r, s) p(y s, r) = rs Γ(s) ys 1 exp( ry) s/r s/r 2

: 1 : (0.6, 0.4) (0.325, 0.675) : (0.5, 0.3, 0.1) (0.28, 0.54, 0.18) 0.6 a = 12, b = 8 a = 6, b = 4 a = 1.2, b = 0.8 a = 0.6, b = 0.4 0.0 0.5 1.0 y 0.0 0.5 1.0 y 0.0 0.5 1.0 y 0.0 0.5 1.0 y 2013 03 07 (2013 03 07 17 :41 ) 59/ 77

OK? 3 2 3 3+2 = 0.60 3.4 1.8 3.4 3.4+1.8 = 0.654 0.654? p = {0.60, 0.40} 2013 03 07 (2013 03 07 17 :41 ) 60/ 77

OK? 3 2 1 3.4 1.8 1.1 3 3+2+1 = 0.50 2 3+2+1 = 0.33 3.4 3.4+1.8+1.1 = 0.540 0.540, 0.286? p = {0.50, 0.33, 0.17} 2013 03 07 (2013 03 07 17 :41 ) 61/ 77

2013 03 07 (2013 03 07 17 :41 ) 62/ 77?? OK http://en.wikipedia.org/wiki/dirichlet_distribution : R glm(y ~ x, family = Gamma) ( ): θ = 1/r =

> ALake <- ArcticLake > ALake$Y <- DR_data(ALake[,1:3]) > DirichReg(Y ~ depth + I(depth^2), ALake) ---------------------------------------- Coefficients for variable no. 1: sand (Intercept) depth I(depth^2) 1.436197-0.007238 0.000132 ---------------------------------------- Coefficients for variable no. 2: silt (Intercept) depth I(depth^2) -0.025970 0.071745-0.000268 ---------------------------------------- Coefficients for variable no. 3: clay (Intercept) depth I(depth^2) -1.793149 0.110791-0.000487 ---------------------------------------- 2013 03 07 (2013 03 07 17 :41 ) 63/ 77 R library(dirichletreg) package GLM

2013 03 07 (2013 03 07 17 :41 ) 64/ 77 library(dirichletreg) 3 9 100 random effects?

2013 03 07 (2013 03 07 17 :41 ) 65/ 77 (ESJ57 http://goo.gl/pekdu)!

2013 03 07 (2013 03 07 17 :41 ) 66/ 77 : above ground mass Y 0 5 10 15 重量 y x 0 5 10 15 below ground mass X

2013 03 07 (2013 03 07 17 :41 ) 67/ 77 : q w 観測できない世界観測できる世界 parameter process observation 無情報事前分布 w q 1- q y x Y X 地上部地下部 階層的事前分布 個体差 測定時の誤差 重量 y x http://tech-jam.com

BUGS code (process ) for (i in 1:N) { Y[i] ~ dnorm(y[i], Tau.err) # X[i] ~ dnorm(x[i], Tau.err) # y[i] <- q[i] * w[i] x[i] <- (1 - q[i]) * w[i] # logit(q[i]) <- a + b * log.w[i] + re[i] w[i] <- exp(log.w[i]) # log.w[i] + log.w[i] ~ dnorm(0, Tau.noninformative) #!! ( ) 2013 03 07 (2013 03 07 17 :41 ) 68/ 77

2013 03 07 (2013 03 07 17 :41 ) 69/ 77 : MCMC 1. BUGS code (model1.txt) 2. R (runbus1.r) 3. R runbus1.r (source("runbugs1.r")) 4. R library(r2winbugs) WinBUGS 5. WinBUGS Markov chain Monte Carlo (MCMC) 6. R

2013 03 07 (2013 03 07 17 :41 ) 70/ 77 : w b MCMC

2013 03 07 (2013 03 07 17 :41 ) 71/ 77 above ground mass Y 0 5 10 15 重量 y x 0 5 10 15 below ground mass X (median) (95% CI)

2013 03 07 (2013 03 07 17 :41 ) 72/ 77

2013 03 07 (2013 03 07 17 :41 ) 73/ 77 C-N (coverage) GIS

2013 03 07 (2013 03 07 17 :41 ) 74/ 77 : offset : :??

2013 03 07 (2013 03 07 17 :41 ) 75/ 77

60 (W30)? 1. ( ) kubo@ees.hokudai.ac.jp 2. ( ) web site URL http://goo.gl/e1cja!! 2013 03 07 (2013 03 07 17 :41 ) 76/ 77

2013 03 07 (2013 03 07 17 :41 ) 77/ 77