2012 11 01 k2 (2012-10-26 16:35 ) 1 6 2 (2012 11 01 k2) (GLM) kubo@ees.hokudai.ac.jp web http://goo.gl/wijx2 web http://goo.gl/ufq2 1 : 2 2 4 3 7 4 9 5 : 11 5.1................... 13 6 14 6.1...................... 14 7 : 14 8 15 9 17 10 19 10.1........................ 19 10.2........................ 20 10.3....................... 24 11 24 12 + 26 12.1 :............ 27 13 28
2012 11 01 k2 (2012-10-26 16:35 ) 2 i y i 1 i y i :!! probability distribution 1 : 1 ( ) 50 0 1, 2 count data R R 1 50 data R R data 2 1 R web data.rdata R load(data.rdata) data R data.rdata 2 print(data)
> data 2012 11 01 k2 (2012-10-26 16:35 ) 3 [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 data 3 50 length() 4 data 50 5 > length(data) # data [1] 50 summary() > summary(data) # data Min. 1st Qu. Median Mean 3rd Qu. Max. 0.00 2.00 3.00 3.56 4.75 7.00 summary(data) Min. Max. data 1st Qu., Median, 3rd Qu. data 25%, 50%, 75% 6 Median Mean sample mean 3.56 7 5 50 R table() > table(data) 0 1 2 3 4 5 6 7 1 3 11 12 10 5 4 4 5 5 6 4 histogram 8 9 3 [1] [26] 1 26 4 R length() () () argument data length(data) R length() length() 5 R # 6 R 7 8 9 seq(-0.5, 9.5, 1) { 0.5, 0.5, 1.5, 2.5,, 8.5, 9.5} hist() breaks -0.5 0.5 0.5 1.5
2012 11 01 k2 (2012-10-26 16:35 ) 4 Histogram of data Frequency 10 data 2 50 R hist() > hist(data, breaks = seq(-0.5, 9.5, 1)) R 2 variability sample variance > var(data) [1] 2.9861 sample standard deviation R > sd(data) [1] 1.7280 > sqrt(var(data)) [1] 1.7280 2 1 2 1 3.56 probability distribution
2012 11 01 k2 (2012-10-26 16:35 ) 5 Poisson distribution 10 random variable i y i y i y i = 2 y i = 2 1 y 1 = 2 parameter 3.56 3.56 3.56 R 3.56 y R > y <- 0:9 > prob <- dpois(y, lambda = 3.56) R dpois(y, lambda = 3.56) prob y 3 > plot(y, prob, type = "b", lty = 2) plot() 4 y prob 3.56 3 4 0.03 3 0.21 3.56 2 3.56 4 5 10 2
2012 11 01 k2 (2012-10-26 16:35 ) 6 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 3 3.56 y {0, 1, 2,, 9} R y p(y λ) prob cbind(y, prob) prob 0.00 0.10 0.20 y 4 λ = 3.56 y prob 3 R plot() type = "b" lty = 2 Histogram of data Frequency 10 data 5 2 y 3.56 4 50
2012 11 01 k2 (2012-10-26 16:35 ) 7 prob 0.00 0.10 0.20 lambda 3.5 7.7 15.1 0 5 10 15 20 y 6 λ λ {3.5, 7.7, 15.1} 11 3 p(y λ) = λy exp( λ) y! p(y λ) λ y 12 y! y 4! 1 2 3 4 λ 3 4 y {0, 1, 2,, } y 1 λ λ 0 y=0 p(y λ) = 1 11 5 hist() lines(y, 50 * prob) 12 p(y λ)
2012 11 01 k2 (2012-10-26 16:35 ) 8 : λ = = λ 6 1 y i {0, 1, 2, } 2 y i 3 13 error measurement error 50 λ 2 14 15 λ 1 λ 2 2 λ 2 3 λ 4 50 λ 3.56 16 13 2.99 14 15 6 16
2012 11 01 k2 (2012-10-26 16:35 ) 9 0 5 10 15 lambda = 2.0 log L = -121.9 0 5 10 15 lambda = 2.4 log L = -109.4 0 5 10 15 lambda = 2.8 log L = -102.0 0 5 10 15 lambda = 3.2 log L = -98.2 0 5 10 15 lambda = 3.6 log L = -97.3 0 5 10 15 lambda = 4.0 log L = -98.5 0 5 10 15 lambda = 4.4 log L = -101.5 0 5 10 15 lambda = 4.8 log L = -106.0 0 5 10 15 lambda = 5.2 log L = -111.8 7 λ lambda logl 2 4 maximum likelihood estimation 17 18 i y i {y 1, y 2, y 3,, y 49, y 50 } = {2, 2, 4,, 2, 3} 50 {y i } Y = {y i } p(y i λ) λ y i 3 λ = 3.56 p(y 1 = 2 λ = 3.56) 0.180 λ λ 17 18
2012 11 01 k2 (2012-10-26 16:35 ) 10 log likelihood -120-110 -100 * ˆλ 2.0 2.5 3.0 3.5 4.0 4.5 5.0 8 50 λ λ = 3.56 ˆλ 7 λ i p(y i λ) 3 {y 1, y 2, y 3 } = {2, 2, 4} 3 0.180 0.180 0.19 = 0.006156 L(λ) 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( λ), i i y i! 19 y 1 2 y 2 2 50 20 L(λ) log likelihood function log L(λ) = i ( y i log λ λ y i k ) log k λ 7 7 log L(λ) λ 8 R 19 20 i {1, 2,, 50}?? i 50 0.5 50
2012 11 01 k2 (2012-10-26 16:35 ) 11 > logl <- function(m) sum(dpois(data, m, log = TRUE)) > lambda <- seq(2, 5, 0.1) > plot(lambda, sapply(lambda, logl), type = "l") 3.5 3.6 λ log L L λ λ ˆλ 21 8 λ log L(λ) λ ˆλ = 1 50 log L(λ) λ i = i y i = y i { yi λ 1 } = 1 λ i y i 50 = = 3.56 ˆλ 3.56 ˆλ maximum likelihood estimator {y 1, y 2, y 3,, y 49, y 50 } = {2, 2, 4,, 2, 3} y i ˆλ = 3.56 maximum likelihood estimate θ y i p(y θ) L(θ Y) = i p(y i θ), log L(θ Y) = i log p(y i θ), ˆθ p(y θ) 22 5 : 21 ˆλ 22
2012 11 01 k2 (2012-10-26 16:35 ) 12 λ = 3.5 y ˆλ = 3.56 y y 9 λ = 3.5 ˆλ = 3.56 y : y 10 9 validation > 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 9 3.5 23 random number 23
2012 11 01 k2 (2012-10-26 16:35 ) 13 sampling λ estimation fitting 9 ˆλ = 3.56 3.56 prediction 24 prediction interval : 5 missing data 25 goodness of prediction 26 10 9 3.56 27 5.1 R 28 29 24 25 26 27 28 29 R
2012 11 01 k2 (2012-10-26 16:35 ) 14 6 : Poisson distribution : binomial distribution : {0, 1, 2,, N} normal distribution : [, + ] gamma distribution : [0, + ] 6.1 7 :
2012 11 01 k2 (2012-10-26 16:35 ) 15 i f i C: T: y i x i 11 i x i f i y i 11 100 i y i body size x i 30 x i 50 i {1, 2,, 50} C 50 i {51, 52,, 100}, T x i x i f i 8 R CSV 31 web data3a.csv R > d <- read.csv("data3a.csv") d (table) 30 31 CSV comma-separated value CSV
2012 11 01 k2 (2012-10-26 16:35 ) 16 R d print(d) 100 32 > 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 100 100 3 y x f d x y > 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 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 f factor R read.csv() CSV C T factor f C T 2 level R Levels f C 1 T 2 read.csv() 33 R class() 32 print() head(d) 6 head(d, 10) 10 33
2012 11 01 k2 (2012-10-26 16:35 ) 17 > 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" R summary() d > 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 summary() summary (numeric) y x f C 50 T 50 9 summary() 34 plot() > plot(d$x, d$y, pch = c(21, 19)[d$f]) > legend("topleft", legend = c("c", "T"), pch = c(21, 19)) 12 x y scatter plot plot(d$f, d$y) 13 box-whisker plot 35 34 35 R 13 ±
2012 11 01 k2 (2012-10-26 16:35 ) 18 d$y 2 4 6 8 10 12 14 C T 7 8 9 10 11 12 d$x 12 y i x i f i C T 2 4 6 8 10 12 14 C T 13 f i plot(d$f, d$y) 75%, 50%, 25% 95% 95%
2012 11 01 k2 (2012-10-26 16:35 ) 19 12 x y f 10 λ λ i x f i x i 36 x i y i f i i y i p(y i λ i ) p(y i λ i ) = λyi i exp( λ i ) y i! 10.1 λ i x i i λ i λ i = exp(β 1 + β 2 x i ) β 1 β 2 parameter β 1 intercept β 2 slope 37 λ i = exp(β 1 + β 2 x i ) 14 38 linear predictor link function GLM λ i log λ i = β 1 + β 2 x i 36 x i y i x i 37 coefficient x i covariate 38 y i x i log x i exp(β 1 + β 2 log x i ) log( ) = β 1 + β 2 log x i x i 0
2012 11 01 k2 (2012-10-26 16:35 ) 20 i λi 0.0 0.5 1.0 1.5 2.0 2.5 {β 1, β 2} = { 2, 0.8} {β 1, β 2} = { 1, 0.4} -4-2 0 2 4 i x i 14 i λ i x i λ i = exp(β 1 + β 2x i) x i x i 7 13 β 1 + β 2 x i β 1 + β 2 x i + β 3 x 2 i {β 1, β 2, β 3 } log λ i = λ i = log link function GLM GLM canonical link function R glm() family GLM λ i = exp( ) 0 14 R 12 10.2 fitting log L ˆβ 1 ˆβ 2 Y log L(β 1, β 2 ) = i log λyi i exp( λ i ) y i! log λ i = β 1 + β 2 x i λ i β 1 β 2
2012 11 01 k2 (2012-10-26 16:35 ) 21 結果を格納するオブジェクト 関数名 確率分布の指定 fit <- glm( y ~ x, モデル式 family = poisson(link = "log"), data = d ) data.frame の指定 リンク関数の指定 ( 省略可 ) 15 glm() {β 1, β 2 } R GLM > fit <- glm(y ~ x, data = d, family = poisson) 39 β 1 β 2 glm() 15 family = poisson 40 fit 41 fit 42 > fit # print(fit) Call: glm(formula = y ~ x, family = poisson, data = d) Coefficients: (Intercept) x 1.29172 0.07566...... summary(fit) 39 argument data = d d data y ~ x formula = y ~ x help(glm) 40 family = poisson(link = "log") poisson family default link function "log" 41 fit names(fit) str(fit) 42 glm() deviance
2012 11 01 k2 (2012-10-26 16:35 ) 22 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 43 (Intercept) β 1 x β 2 Estimate ˆβ 1 = 1.29 ˆβ 2 = 0.0757 Std. Error standard error, SE ˆβ 1 ˆβ 2 SE 44 SE z value SE Wald 16 z Wald Wald statistics Pr(> z ) glm() 45 z 1 z ˆβ 1 ˆβ 2 16 ˆβ 2 2 ˆβ 2 z Pr(> z ) Pr(> z ) P statistical test 46 confidence interval 47 Wald 43 R.Rprofile options(show.signif.stars = FALSE) 44 Wald 45 glm() family glm() family = gaussian t 46 47 α% α%
2012 11 01 k2 (2012-10-26 16:35 ) 23 β 2 β 1 0.0 0.5 1.0 1.5 16 R glm() β 2 2 Pr(> z ) Pr(> z = 0.03358) β 1 Pr(> z ) Pr(> z = 0.00038) d$y 2 4 6 8 10 12 14 7 8 9 10 11 12 d$x 17 λ 12 λ maximum log likelihood goodness of fit log L(β 1, β 2 ) { ˆβ 1, ˆβ 2 } R > loglik(fit) log Lik. -235.3863 (df=2) -235.4 (df=2) degrees of freedom 2 2 β 1 β 2
2012 11 01 k2 (2012-10-26 16:35 ) 24 10.3 x λ prediction x λ { ˆβ 1, ˆβ 2 } λ = exp(1.29 + 0.0757x) R > plot(d$x, d$y, pch = c(21, 19)[d$f]) > xx <- seq(min(d$x), max(d$x), length = 100) > lines(xx, exp(1.29 + 0.0757 * xx), lwd = 2) 17 λ predict() > yy <- predict(fit, newdata = data.frame(x = xx), type = "response") > lines(xx, yy, lwd = 2) 11 f i (p.16 ) R factor C T C 1 T 2 48 GLM R dummy variable 49 glm() f x i f i λ i = exp(β 1 + β 3 d i ) β 1 β 3 f i d i : 48 49 R contrasts
2012 11 01 k2 (2012-10-26 16:35 ) 25 d i = { 0 (f i = C ) 1 (f i = T ) i f i = C f i = T R λ i = exp(β 1 ) λ i = exp(β 1 + β 3 ) glm() glm() > fit.f <- glm(y ~ f, data = d, family = poisson) fit.f Call: glm(formula = y ~ f, family = poisson, data = d) Coefficients: (Intercept) ft 2.05156 0.01277...... Coefficients f i ft f i T f i C T 2 R f i C T i f i C T λ i = exp(2.05 + 0) = exp(2.05) = 7.77 λ i = exp(2.05 + 0.0128) = exp(2.0628) = 7.87 > loglik(fit.f) log Lik. -237.6273 (df=2) (p.23 ) x i -235.4
2012 11 01 k2 (2012-10-26 16:35 ) 26 3 A B f i f i {C, TA, TB} 3 R 2 λ i = exp(β 1 + β 3 d i,a + β 4 d i,b ) β 3 A β 4 B { 0 (f i TA ) d i,a = 1 (f i TA ) { 0 (f i TB ) d i,b = 1 (f i TB ) 3 R glm(y ~ f,...) β 3 β 4 fta ftb 12 + x i f i 50 GLM log λ i = β 1 + β 2 x i + β 3 d i β 1 β 2 x i β 3 (2 f i d i ) 51 R glm() x + f > fit.all <- glm(y ~ x + f, data = d, family = poisson) fit.all 50 multiple regression 51
2012 11 01 k2 (2012-10-26 16:35 ) 27 Call: glm(formula = y ~ x + f, family = poisson, data = d) Coefficients: (Intercept) x ft 1.2631 0.0801-0.0320 Degrees of Freedom: 99 Total (i.e. Null); Null Deviance: 89.51 Residual Deviance: 84.81 AIC: 476.6 97 Residual ft > loglik(fit.all) log Lik. -235.2937 (df=3) p.23 x i (-235.4) 12.1 : 10 + glm() glm(y ~ x + f,...) x i i f i C λ i = exp(1.26 + 0.08x i ) T λ i = exp(1.26 + 0.08x i 0.032) λ i = exp(1.26) exp(0.08x i ) exp( 0.032) = x i x i 1 λ i exp(0.08 1) = 1.08
2012 11 01 k2 (2012-10-26 16:35 ) 28 (A) (B) λi 5 10 15 5 10 15 5 10 15 20 x i 5 10 15 20 x i 18 (A) (B) 3 3 exp( 0.032) = 0.969 0.969 λ i 18 (A) 52 x i R identity link function 53 λ i λ i = 1.27 + 0.661x i 0.205d i 18 (B) 54 0.1 1000 0.205-0.1 999.8 13 GLM 12.1 generalized linear model, LM 52 3 0.032 3 3 53 glm() family = poisson(link = "identity") λ 54 (A) 3 0.205 3
2012 11 01 k2 (2012-10-26 16:35 ) 29 (A) y -2 0 2 4 6 0.5 1.0 1.5 2.0 x (B) y -2 0 2 4 6 0.5 1.0 1.5 2.0 x 19 GLM x x {0.5, 1.1, 1.7} y (A) y β 1 + β 2x (B) y exp(β 1 + β 2x) general lienar model LM linear regression 19 (x i, y i ) 19 (A) GLM 55 : {x 1, x 2,, x n } {y 1, y 2,, y n } X = {x i } Y = {y i } Y µ i σ i µ i = β 1 + β 2 x i LM GLM 56 55?? 56 LM multiple regression x i ANOVA
2012 11 01 k2 (2012-10-26 16:35 ) 30 19 y 0, 1, 2 x y y 19 A : 19 57 58 19 (B) GLM y?? y log y y 59 y GLM GLM ANCOVA 57 R glm() family = gaussian family = gaussian(link = "log") 58 R 2 P 59 log 0 =