/ 60 : 1. GLM? 2. A: (pwer functin) x y?

Size: px

Start display at page:

Download "/ 60 : 1. GLM? 2. A: (pwer functin) x y?"

えのみつだ
5 years ago
Views:

1 / 60 ( ) GLM 1. GLM :, link,, deviance (20 ) 2. GLM : (60 ) 3. GLM ( ): ffset (40 ) kub/ce/ecsj2009.html

2 / 60 : 1. GLM? 2. A: (pwer functin) x y?

3 R OS free sftware S R (2007) The R-Tips (2005) Statistics: An Intrductin Using R M. Crawley (2005) / 60

4 1. GLM?

5 : : ( ) / 60

6 / 60 : : (linear predictr): ( ) = ( ) + ( 1) ( 1) + ( 2) ( 2) + ( 3) ( 3) + ( )

7 / 60 (GLMM) GLM + randm effects ( ) (GLM) +

8 GLM / 60

9 / 60 :?! ( )? :? :?

10 / 60! GLM

11 / 60 (Pissn distributin)? lambda = 1.4 y i {0, 1, 2,, } (paramter: λ) λ y exp( λ) 0.0 lambda = y! prbability λ λ lambda = : y

12 (generalized linear mdel; GLM) link : link f () = f ( ) : β 0 + β 1 x 1 + β 2 + x i β i x i (cefficient) ({x i } {y i }) {β i } GLM / 60

13 R : glm() ( ) rbinm() glm(family = binmial) rbinm() glm(family = binmial) rpis() glm(family = pissn) rnbinm() glm.nb() in library(mass) ( ) rgamma() glm(family = gamma) rnrm() glm(family = gaussian) glm() GLM GLM / 60

14 / 60 R glm() :? ( z):? link : z (y)? family:?

15 glm() (1) family: pissn, (0, 1, 2, ) link : "lg" family = pissn link ( z): y ~ x family = pissn(link = "lg")? / 60

16 glm() (2) family: pissn, link : "lg" ( z): y ~ x z = a + bx a, b λ lg(λ) = z λ = exp(z) = exp(a + bx) λ : y Pis(λ) / 60

17 / A: (pwer functin)

18 A:? plant weight (g) number f flwers x y? 100 : i = 1, 2,, 100 x i y i / 60

19 : x y? x (0, 0) :?? / 60

20 1. y i λ i : y i Pis(λ i ) 2. λ i x i : λ i = Ax b i λ i = Ax b i λ i = exp(lg(a) + b lg(x i )) a = lg(a) lg(λ i ) = a + b lg(x i ) / 60

21 GLM! family: pissn, link : "lg" : y ~ lg.x x lg.x z = a + b lg.x a, b λ lg(λ) = z λ = exp(z) = exp(a + b lg.x) λ : y Pis(λ) / 60

22 / 60 R data.frame d = > lad("d.rdata") > head(d) # 6 x y lg.x > d$lg.x <- lg(d$x) > head(d) x y lg.x

23 glm() / 60

24 / 60 R glm() > fit <- glm(y ~ lg.x, data = d, family = pissn) > print(summary(fit)) Call: glm(frmula = y ~ lg.x, family = pissn, data = d) (......) Cefficients: Estimate Std. Errr z value Pr(> z ) (Intercept) lg.x (......) Cefficients

25 GLM plant weight (g) number f flwers plant weight (g) number f flwers / 60

26 : 1. : 2. GLM :, link, 3. R glm() deviance / 60

27 / GLM :, link,, deviance (20 ) 2. GLM : (60 ) 3. GLM ( ): ffset (40 )

28 / 60 : 1. : 2. ffset :

29 / 60 1.

30 / 60 : ( )

31 : ( ) / 60

32 / x x ( )

33 ? ( ) / 60

34 / 60 : : (linear predictr): ( ) = ( ) + ( 1) ( 1) + ( 2) ( 2) + ( 3) ( 3) + ( )

35 / 60

36 / 60! GLM

37 (generalized linear mdel; GLM) link : link f () = f ( ) : β 0 + β 1 x 1 + β 2 + x i β i x i (cefficient) ({x i } {y i }) {β i } GLM / 60

38 / 60 i N i k i i p i = k i /N i j p j = k j /N j i j p

39 ? / : ? ( ) / 60

40 / 60 : specific leaf area (SLA) : ffset : N k :

41 / ffset

42 B:? x {0.1, 0.2,, 1.0} 10 glm(..., family = pissn) / 60

43 ?!! x A = /! glm() ffset / 60

44 / 60 R data.frame: Area, x, y > lad("d2.rdata") > head(d, 8) # 8 Area x y

45 / 60 vs plt(d$x, d$y / d$area) d$y/d$area d$x?

46 / 60 A vs y plt(d$area, d$y) d$y d$area A y

47 / 60 x ( ) plt(d$area, d$y, cex = d$x * 2) d$y d$area?

48 / 60 x y x!

49 = 1. i y i λ i : y i Pis(λ i ) 2. λ i A i x i λ i = A i exp(a + bx i ) λ i = exp(a + bx i + lg(a i )) lg(λ i ) = a + bx i + lg(a i ) lg(a i ) ffset / 60

50 GLM! ffset family: pissn, link : "lg" : y ~ x ffset : lg(area) z = a + b x + lg(area) a, b λ lg(λ) = z λ = exp(z) = exp(a + b x + lg(area)) λ : y Pis(λ) / 60

51 glm() / 60

52 / 60 R glm() > fit <- glm(y ~ x, family = pissn(link = "lg"), data = d, ffset = lg(area)) > print(summary(fit)) Call: glm(frmula = y ~ x, family = pissn(link = "lg"), data = d, ffset = lg(area)) (......) Cefficients: Estimate Std. Errr z value Pr(> z ) (Intercept) x e-06 Cefficients

53 / 60 d$y d$area x = 0.9, x = 0.1 glm()

54 / 60 : glm() ffset ffset = exp( ) d$y d$area

55 / 60 : ffset

56 GLM summary(glm(...)) Wald Pr deviance AIC (verdispersin) GLMM kub/ce/ / 60

57 / 60 GLM (1)? y = 0, 1, 2, 3, (y ) (family = pissn) y = {0, 1}, y = {0, 1, 2,, N} (family = binmial) (family = Gamma) (family = gaussian)

58 R : glm() ( ) rbinm() glm(family = binmial) rbinm() glm(family = binmial) rpis() glm(family = pissn) rnbinm() glm.nb() in library(mass) ( ) rgamma() glm(family = gamma) rnrm() glm(family = gaussian) glm() glm.nb() MASS library GLM / 60

59 GLM (2)!! GLMM! GLMM/ randm effects GLM / 60

60 kub/ce/ / 60

/ 55 2 : : (GLM) 1. 1/23 ( )? GLM? (GLM ) 2.! 1/25 ( ) ffset (GLM )

2012 01 25 1/ 55 ( II) : (2012 1 ) 2 2 (GLM) 2012 01 25! kub@ees.hkudai.ac.jp http://g.gl/76c4i 2012 01 25 2/ 55 2 : : (GLM) 1. 1/23 ( )? GLM? (GLM ) 2.! 1/25 ( ) ffset (GLM ) 2012 01 25 3/ 55 1. : 2.