一般化線形 (混合) モデル (2) - ロジスティック回帰と GLMM

Size: px

Start display at page:

Download "一般化線形 (混合) モデル (2) - ロジスティック回帰と GLMM"

たみえこしの
7 years ago
Views:

1 .. ( ) (2) GLMM [email protected] I : :29 kubostat2013ou2 ( ( ) (2) / 74

2 I.1 N k.2 binomial distribution logit link function.3.4! offset.5 GLM!.6.7 kubostat2013ou2 ( ( ) (2) / 74

3 II.8 r i.9 GLMM (A) f i =C (B) f i =T yi x i x i kubostat2013ou2 ( ( ) (2) / 74

4 6 7 : : kubostat2013ou2 ( ( ) (2) / 74

5 kubostat2013ou2 ( ( ) (2) / 74

6 ? (GLM) (Poisson regression) (logistic regression) (linear regression) kubostat2013ou2 ( ( ) (2) / 74

7 (GLM)??? kubostat2013ou2 ( ( ) (2) / 74

8 GLM : : e.g., β 1 + β 2 x i : kubostat2013ou2 ( ( ) (2) / 74

9 GLM logistic : : e.g., β 1 + β 2 x i : logit yi x i kubostat2013ou2 ( ( ) (2) / 74

10 N k 1. N k y i {0, 1, 2,, 8} kubostat2013ou2 ( ( ) (2) / 74

11 N k? 8 y! f i C: T: i N i = 8 y i = 3 (alive) (dead) x i kubostat2013ou2 ( ( ) (2) / 74

N k data4a.csv CSV (comma separated value) format file R : > d <- read.csv("data4a.csv") or > d <- read.csv( + "http://hosho.ees.

12 N k data4a.csv CSV (comma separated value) format file R : > d <- read.csv("data4a.csv") or > d <- read.csv( + " d data frame ( ) kubostat2013ou2 ( ( ) (2) / 74

13 N k data frame d > summary(d) N y x f Min. :8 Min. :0.00 Min. : C:50 1st Qu.:8 1st Qu.:3.00 1st Qu.: T:50 Median :8 Median :6.00 Median : Mean :8 Mean :5.08 Mean : rd Qu.:8 3rd Qu.:8.00 3rd Qu.: Max. :8 Max. :8.00 Max. : kubostat2013ou2 ( ( ) (2) / 74

14 N k > plot(d$x, d$y, pch = c(21, 19)[d$f]) > legend("topleft", legend = c("c", "T"), pch = c(21, 19)) yi C T x i? kubostat2013ou2 ( ( ) (2) / 74

15 binomial distribution logit link function 2. binomial distribution logit link function kubostat2013ou2 ( ( ) (2) / 74

16 binomial distribution logit link function : N y ( N y p(y N, q) = ( ) N q y (1 q) N y y ) N y p(y i 8, q) q = 0.1 q = 0.8 q = y i kubostat2013ou2 ( ( ) (2) / 74

17 binomial distribution logit link function (z i : e.g. z i = β 1 + β 2 x 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 q = z 1 1+exp( z) kubostat2013ou2 ( ( ) (2) / 74

18 binomial distribution logit link function {β 1, β 2 } = {0, 2} (A) β 2 = 2 β 1 (B) β 1 = 0 β 2 q (A) β 2 = 2 β 1 = 2 β 1 = x β 1 = (B) β 1 = 0 β 2 = 4 β 2 = x β 2 = 1 {β 1, β 2 } x q 0 q 1 kubostat2013ou2 ( ( ) (2) / 74

19 binomial distribution logit link function logit link function logistic 1 q = 1 + exp( (β 1 + β 2 x)) = logistic(β 1 + β 2 x) logit logit(q) = log q 1 q = β 1 + β 2 x logit logistic logistic logit logit is the inverse function of logistic function, vice versa kubostat2013ou2 ( ( ) (2) / 74

20 binomial distribution logit link function R β 1 β 2 (A) f i =C (B) y x x > glm(cbind(y, N - y) ~ x + f, data = d, family = binomial)... Coefficients: (Intercept) x ft kubostat2013ou2 ( ( ) (2) / 74

21 binomial distribution logit link function : (A) f i =C (B) f i =T yi x i x i kubostat2013ou2 ( ( ) (2) / 74

22 3. kubostat2013ou2 ( ( ) (2) / 74

23 ? logit(q) = log q 1 q = β 1 + β 2 x + β 3 f + β 4 xf... in case that β 4 < 0, sometimes it predicts... y T C x kubostat2013ou2 ( ( ) (2) / 74

24 glm(y ~ x + f,...) glm(y ~ x + f + x:f,...) (A) (B) y T C T C x x kubostat2013ou2 ( ( ) (2) / 74

25 ! offset 4.! offset kubostat2013ou2 ( ( ) (2) / 74

26 ! offset? / : ? ( ) kubostat2013ou2 ( ( ) (2) / 74

27 ! offset : N k : : specific leaf area (SLA) : offset! kubostat2013ou2 ( ( ) (2) / 74

28 ! offset kubostat2013ou2 ( ( ) (2) / 74

29 ! offset offset : x {0.1, 0.2,, 1.0} 10 glm(..., family = poisson) kubostat2013ou2 ( ( ) (2) / 74

30 ! offset?! x A = /! glm() offset kubostat2013ou2 ( ( ) (2) / 74

31 ! offset R data.frame: Area, x, y > load("d2.rdata") > head(d, 8) # 8 Area x y kubostat2013ou2 ( ( ) (2) / 74

32 ! offset vs > plot(d$x, d$y / d$area) d$y/d$area d$x kubostat2013ou2 ( ( ) (2) / 74

33 ! offset A vs y > plot(d$area, d$y) d$y d$area A y kubostat2013ou2 ( ( ) (2) / 74

34 ! offset x ( ) > plot(d$area, d$y, cex = d$x * 2) d$y d$area? kubostat2013ou2 ( ( ) (2) / 74

35 ! offset x y x kubostat2013ou2 ( ( ) (2) / 74

36 ! offset = 1. i y i λ i : y i Pois(λ i ) 2. λ i A i x i λ i = A i exp(β 1 + β 2 x i ) λ i = exp(β 1 + β 2 x i + log(a i )) log(λ i ) = β 1 + β 2 x i + log(a i ) log(a i ) offset ( β ) kubostat2013ou2 ( ( ) (2) / 74

37 ! offset GLM! family: poisson, link : "log" : y ~ x offset : log(area) z = β 1 + β 2 x + log(area) a, b λ log(λ) = z λ = exp(z) = exp(β 1 + β 2 x + log(area)) λ : kubostat2013ou2 ( ( ) (2) / 74

38 ! offset glm() kubostat2013ou2 ( ( ) (2) / 74

39 ! offset R glm() > fit <- glm(y ~ x, family = poisson(link = "log"), data = d, offset = log(area)) > print(summary(fit)) Call: glm(formula = y ~ x, family = poisson(link = "log"), data = d, offset = log(area)) (......) Coefficients: Estimate Std. Error z value Pr(> z ) (Intercept) x e-06 kubostat2013ou2 ( ( ) (2) / 74

40 ! offset d$y x = 0.9 light environment x = 0.1 dark environment d$area glm() kubostat2013ou2 ( ( ) (2) / 74

41 ! offset : glm() offset offset = exp( ) d$y d$area kubostat2013ou2 ( ( ) (2) / 74

42 ! offset : N k : : specific leaf area (SLA) : offset! kubostat2013ou2 ( ( ) (2) / 74

43 GLM! 5. GLM! (overdispersion)? kubostat2013ou2 ( ( ) (2) / 74

44 GLM! :?! (A) i N i = 8 (B) 100 x i y i y i = 3 x i {2, 3, 4, 5, 6} yi x i kubostat2013ou2 ( ( ) (2) / 74

45 GLM! N y? number of alive seeds y i : : β 1 + β 2 x i : logit number of leaves x i kubostat2013ou2 ( ( ) (2) / 74

46 GLM! GLM! yi (A) β 2 (B)! x i = 4 y i x i y i kubostat2013ou2 ( ( ) (2) / 74

47 6.? kubostat2013ou2 ( ( ) (2) / 74

48 (overdispersion)? (A) Not or less overdispersed (B) Overdispersed!! y i y i kubostat2013ou2 ( ( ) (2) /

49 GLM GLM does not take into account individual differences kubostat2013ou2 ( ( ) (2) / 74

50 Almost all real data are overdispersed! kubostat2013ou2 ( ( ) (2) / 74

51 7. kubostat2013ou2 ( ( ) (2) / 74

52 number of alive seeds y i : : β 1 + β 2 x i + r i : logit number of leaves x i kubostat2013ou2 ( ( ) (2) / 74

53 i r i qi r i > 0 r i = 0 r i < x i kubostat2013ou2 ( ( ) (2) / 74

54 {r i } s = 1.0 s = r i s = 3.0 p(r i s) = 1 2πs 2 exp ( r2 i 2s 2 p(r i s) r i r i r i kubostat2013ou2 ( ( ) (2) / 74 )

55 r i (A) (B) p(r i s) s = {r i } s = 3.0 I III I II I I r i II I II I II I I I III I II I I I I I r i 1 q i = 1+exp( r i ) y i p(y i q i ) y i kubostat2013ou2 ( ( ) (2) / 74

56 > # defining logistic function > logistic <- function(z) { 1 / (1 + exp(-z)) } > # random numbers following binomial distribution > rbinom(100, 8, prob = logistic(0)) > # random numbers following Gausssian distribution > rnorm(100, mu = 0, sd = 0.5) > r <- rnorm(100, mu = 0, sd = 0.5) > # random numbers following...? > rbinom(100, 8, prob = logistic(0 + r)) kubostat2013ou2 ( ( ) (2) / 74

57 Generalized Linear Mixed Model (GLMM) Mixed : β 1 + β 2 x i + r i fixed effects: β 1 + β 2 x i random effects: +r i fixed? random?? kubostat2013ou2 ( ( ) (2) / 74

58 : fixed effects random effects kubostat2013ou2 ( ( ) (2) / 74

59 : (linear mixed model, LMM) random effects : : kubostat2013ou2 ( ( ) (2) / 74

60 global parameter, local parameter? Generalized Linear Mixed Model (GLMM) : β 1 + β 2 x i + r i fixed effects: β 1 + β 2 x i global parameter s global parameter random effects: +r i local parameter i ( ) global/local parameter kubostat2013ou2 ( ( ) (2) / 74

61 全データ個体個体 3 3 のデータのデータ個体 1 のデータ個体個体 3 3 のデータのデータ個体 2 のデータ {r 1, r 2, r 3,..., r 100 } β 1 β 2 local parameter global parameter? kubostat2013ou2 ( ( ) (2) / 74 s

62 r i 8. r i kubostat2013ou2 ( ( ) (2) / 74

63 r i r i local parameters: {r 1, r 2,, r 100 } 100 r i > d <- read.csv("data.csv") > head(d) N y x id kubostat2013ou2 ( ( ) (2) / 74

64 r i r i y i ( ) 8 p(y i β 1, β 2 ) = q yi i (1 q i) 8 y i r i i r i L i = p(r i s) = β 1, β 2, s y i ) 1 ( exp r2 i 2πs 2 2s 2 p(y i β 1, β 2, r i ) p(r i s)dr i L(β 1, β 2, s) = i L i kubostat2013ou2 ( ( ) (2) / 74

65 r i global parameter local parameter Generalized Linear Mixed Model (GLMM) Mixed : β 1 + β 2 x i + r i global parameter fixed effects: β 1, β 2 : s local parameter random effects: {r 1, r 2,, r 100 } kubostat2013ou2 ( ( ) (2) / 74

66 r i r i kubostat2013ou2 ( ( ) (2) / 74

67 r i r r = 2.20 q = 0.10 y r = 0.60 q = 0.35 y r = 1.00 q = 0.73 y r = 2.60 q = 0.93 y r p(r s) p(r) = 0.10 r p(r) = 0.13 r p(r) = 0.13 r p(r) = 0.09 r y kubostat2013ou2 ( ( ) (2) / 74

68 r i r r = 1.10 λ = 0.55 y r = 0.30 λ = 1.22 y r = 0.50 λ = 2.72 y r = 1.30 λ = 6.05 y r p(r s) p(r) = 0.22 r p(r) = 0.38 r p(r) = 0.35 r p(r) = 0.17 r y kubostat2013ou2 ( ( ) (2) / 74

69 r i glmmml package GLMM > install.packages("glmmml") # if you don t have glmmml > library(glmmml) > glmmml(cbind(y, N - y) ~ x, data = d, family = binomial, + cluster = id) > d <- read.csv("data.csv") > head(d) N y x id kubostat2013ou2 ( ( ) (2) / 74

70 r i GLMM : ˆβ1, ˆβ 2, ŝ > glmmml(cbind(y, N - y) ~ x, data = d, family = binomial, + cluster = id)...(snip)... coef se(coef) z Pr(> z ) (Intercept) e-06 x e-06 Scale parameter in mixing distribution: 2.49 gaussian Std. Error: Residual deviance: 264 on 97 degrees of freedom AIC: 270 ˆβ 1 = 4.13, ˆβ 2 = 0.99, ŝ = 2.49 kubostat2013ou2 ( ( ) (2) / 74

71 r i GLMM (A) (B) x = 4 yi x i y kubostat2013ou2 ( ( ) (2) / 74

72 GLMM 9. GLMM kubostat2013ou2 ( ( ) (2) / 74

73 GLMM + GLMM I (A) pot A pot A pot B pot B (B) logitq i = β 1 + β 2 x i (GLM) q i : logitq i = β 1 + β 2 x i + r i (A) (B) kubostat2013ou2 ( ( ) (2) / 74

74 GLMM + GLMM II (C) pot A pot B logitq i = β 1 + β 2 x i + r j (D) pot A pot B logitq i = β 1 + β 2 x i + r i + r j kubostat2013ou2 ( ( ) (2) / 74

75 GLMM GLMM random effects global parameter local parameter GLMM global parameter local parameter local parameter (e.g. + ) kubostat2013ou2 ( ( ) (2) / 74

kubostat2017e p.1 I 2017 (e) GLM logistic regression : : :02 1 N y count data or

kubostat2017e p.1 I 2017 (e) GLM logistic regression : : :02 1 N y count data or kubostat207e p. I 207 (e) GLM [email protected] https://goo.gl/z9ycjy 207 4 207 6:02 N y 2 binomial distribution logit link function 3 4! offset kubostat207e (https://goo.gl/z9ycjy) 207 (e) 207 4