講義のーと : データ解析のための統計モデリング. 第３回

Size: px

Start display at page:

Download "講義のーと : データ解析のための統計モデリング. 第３回"

やすもりあると
5 years ago
Views:

1 Title 講義のーと : データ解析のための統計モデリング Author(s) 久保, 拓弥 Issue Date 2008 Doc URL Type learningobject Note この講義資料は, 著者のホームページードできます Note(URL) Additional There Information are other files related to this item in HUSCAP File Information kubostat2008c.pdf ( 第 3 回 ) Instructions for use Hokkaido University Collection of Scholarly and Aca

2 ( :11 ) 1 ( ) 5 (+2) 3 ( ) (GLM) 1 kubo@ees.hokudai.ac.jp (GLM)? 1 2. GLM : link 4 3. GLM : 5 4. : glm(..., family = poisson) Deviance? λ λ λ = ( β 1 + β 2 ( 1) + β 3 ( 2) + ) {β j } = {β 1, β 2, }

3 ( :11 ) 2 1. (GLM)? 1 1 Bayes Bayes MCMC random effects (GLMM) (GLM) +, fixed effects / 8 1: ( ) (general linear model) ANOVA ANCOVA = 2 2.? general linear model

4 ( :11 ) 3 ( ) ? 0, 1, 2, 3,... (1) ( )? (2)? 2: ( ) ( 3) (generalized linear model; GLM) 4 : 5 link ( : ) GLM 6 6.

5 ( :11 ) 4 3: ( ) () 7 logistic GLM : link GLM (probablistic distribution) link (link function) (linear predictor) GLM R glm() logistic 9.? ( ) rbinom() glm(family = binomial) rbinom() glm(family = binomial) rpois() glm(family = poisson) rnbinom() glm.nb() ( ) rgamma() glm(family = gamma) rnorm() glm(family = gaussian)

6 ( :11 ) 5 link (link function) ( ) µ µ = β 1 + β 2 ( 1) + β 3 ( 2) + {β j } (explanatory variable) β j (1 ) j z = j β j (1 ) z GLM link (response variable) µ = z 12 log link log µ = z µ = exp ( β 1 + β 2 ( 1) + β 3 ( 2) + ) ( ) 12. R glm() identity link link GLM link canonical link function (R glm() family ) canonical ( ) 14 link glm() family link (m ) 14. ( ) binomial logit µ(1 µ) ( ): Np(1 p) N binomial logit µ(1 µ) ( ) p poisson log µ ( ) gamma log? µ 2 gaussian identity

7 ( :11 ) 6 3. GLM : 4 GLM f i i C: T: y i x i 4: i (i = 1, 2,, 100) ( ) i (i = 1, 2, 3,, 100, 100 ) y i ; 15 x i ; 50 (i = 1, 2,, 50) (C ) 50 (i = 51, 52,, 100) (T ) spread sheet 16 CSV (comma separated values) 17 data3a.csv 18 R 19 > d <- read.csv("data3a.csv") data.frame d 20 data.frame table R? d ENTER ( ) 17. CSV 18. web page 19. R... R setwd() data3a.csv 20.

8 y x f C C C... ( ) T T ( :11 ) head(d) data.frame head(d, 10) y x f 10 R d x y > d$x [1] [9] ( )... [97] > d$y [1] [17] ( )... [97] 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 (factor) R read.csv() CSV table C T 23 factor f C T (levels) () factor C, T read.csv() 24 R class() ( 24.

9 ( :11 ) 8 ) 25 > 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 data.frame > summary(d) y x f Min. : 2.00 Min. : C:50 1st Qu.: st Qu.: T:50 Median : 8.00 Median : Mean : 7.83 Mean : rd Qu.: rd Qu.: Max. :15.00 Max. : R # data.frame summary f C 50 T 50 summary plot() factor vector summary() x y! (scatter plot) > plot(d$x, d$y, pch = c(21, 19)[d$f]) > legend("topleft", legend = c("c", "T"), pch = c(21, 19)) d$y C T d$x f

10 ( :11 ) 9 > plot(d$f, d$y) C T R plot() 27 (box-whisker %, 50%, 25% plot) ( ) 2.5% ( ) 95% R plot() summary() : glm(..., family = poisson) d y x f ( ) i y i 0, 1, 2 ( ) i λ i = exp(β 1 + β 2 x i ) x i f i f i

11 ( :11 ) 10 log λ i = β 1 + β 2 x i link log link z = β 1 + β 2 x i log link? λ i λ i > 0 λ i = exp(z i ) z (= β 1 + β 2 x i ) exp() z λ i λ i β 1, β 2, x i > xi <- seq(-4, 5, 0.1) λ i 0 > plot(xi, exp( * xi), type = "l", ylab = "lambda", lwd = 2) > lines(xi, exp( * xi), lwd = 2, lty = 2) # > abline(v = 0, lwd = 2, lty = 3) lambda xi λ i = exp(β 1 + β 2 x i ) i y i p(y i {β j }, {x i }) = λy i i exp( λ i ) y i! L({β j } {y i }, {x i }) = 100 i=1 λ y i i exp( λ i ) y i! 100 log L({β j } {y i }, {x i }) = i=1 { λ y i log i exp( λ i ) y i! }

12 ( :11 ) 11 ˆβ 1 ˆβ 2 (GLM) R GLM > fit <- glm(y ~ x, data = d, family = poisson) (fitting) 31 ( 31. ) fit names(fit) 5: glm() 5 family = poisson log link fit > fit # print(fit) Call: glm(formula = y ~ x, family = poisson, data = d) Coefficients: (Intercept) x family = poisson(link = "log") link poisson family default link function 32 "log" 32. canonical link function ( link ) Degrees of Freedom: 99 Total (i.e. Null); 98 Residual Null Deviance: Residual Deviance: AIC: 474.8

13 ( :11 ) 12 summary() > summary(fit) Call: glm(formula = y ~ x, family = poisson, data = d) Deviance Residuals: Min 1Q Median 3Q Max Coefficients: Estimate Std. Error z value Pr(> z ) (Intercept) *** x * --- Signif. codes: 0 *** ** 0.01 * (Dispersion parameter for poisson family taken to be 1) Null deviance: on 99 degrees of freedom Residual deviance: on 98 degrees of freedom AIC: Number of Fisher Scoring iterations: 4 Deviance summary() (coefficient) Coefficients: Estimate Std. Error z value Pr(> z ) (Intercept) *** x * Coefficients ( ) (Intercept) ( ) ( ) β 1 x β 2 Estimate ˆβ ˆβ Std. Errror?

14 ( :11 ) 13 ˆβ 1 ˆβ 2 ˆβ ? (SE) z value 34 ( ) / SE (Wald ) 6 ˆβ ˆβ 1 SE ˆβ 2 SE ˆβ 2 ˆβ :? 34. z Wald 35 Pr(> z ) z value 35. (Neyman- ˆβ Pearson ) 36 ( ) Wald R summary(glm(...)) 37 (prediction) R

15 ( :11 ) 14 > plot(d$x, d$y, pch = c(21, 19)[d$f]) > xx <- seq(min(d$x), max(d$x), length = 100) > lines(xx, exp( * xx), lwd = 2) d$y d$x lines(xx, exp( * xx), lwd = 2) λ > yy <- predict(fit, newdata = data.frame(x = xx), type = "response") > lines(xx, yy, lwd = 2) predict() log L({β j } ) { ˆβ i } { ˆβ i } log L({β j } ) loglik(fit) > loglik(fit) # fit <- glm(y ~ x,...) log Lik (df=2) (df=2) β 1 β x i f i

16 ( :11 ) 15 f (factor) > 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 R glm() numeric, integer, factor > fit.f <- glm(y ~ f, data = d, family = poisson) > fit.f Call: glm(formula = y ~ f, family = poisson, data = d) Coefficients: (Intercept) ft Degrees of Freedom: 99 Total (i.e. Null); 98 Residual Null Deviance: Residual Deviance: AIC: glm(y ~ f,...) Coefficients (Intercept) ( ) f ft f C ( ) T ( ) 2 ( ) 1 C 2 39 i f C T λ i exp( ) = exp(2.05) λ i exp( ) = exp(2.0628) 39. 2

17 ( :11 ) 16 > loglik(fit.f) # fit.f <- glm(y ~ f,...) log Lik (df=2) 6. + f glm() x > fit.full <- glm(y ~ x + f, data = d, family = poisson) > fit.full Call: glm(formula = y ~ x + f, family = poisson, data = d) Coefficients: (Intercept) x ft Degrees of Freedom: 99 Total (i.e. Null); 97 Residual Null Deviance: Residual Deviance: AIC: ft log link i x i C λ i exp( x i ) T λ i exp( x i 0.032) exp( x i 0.032) = exp( x i ) exp( 0.032) = () ( ) glm(y ~ x + f,...)

18 ( :11 ) 17 exp( 0.032) > loglik(fit.full) # fit.full <- glm(y ~ x + f,...) log Lik (df=3) 7. Deviance?? loglik() deviance 41 log L({β j } ) log L { ˆβ i } log L log L summary(glm(y...)) log L 42 deviance deviance D = 2 log L 41. deviance 42. loglik() log L ? x i λ i = exp(β 1 + β 2 x i )

19 ( :11 ) 18 log L deviance (D = 2 log L ) glm()... ( )... Null Deviance: Residual Deviance: AIC: Null Deviance Residual Deviance AIC Residual Deviance 44 Residual deviance deviance D 44. Residual D ( D) 45 ( D) 45. Deviance residual deviance? full model (!) deviance 46 data.frame d > d$y [1] ( )... x i λ i d$y full model 46. β 1 β 2 2 1, 2, 3 y y y () ?

20 ( :11 ) > sum(dpois(d$y, lambda = d$y, log = TRUE)) [1] log L deviance D = 2 log L ( D) residual deviance D ( D) = 85.0 glm() Residual Deviance: AIC: residual deviance residual deviance residual deviance? Residual deviance full model () λ i = exp(β 1 ) 1 R null model 50 null model > fit.null <- glm(formula = y ~ 1, family = poisson, data = d) > fit.null Call: glm(formula = y ~ 1, family = poisson, data = d) Coefficients: (Intercept) Degrees of Freedom: 99 Total (i.e. Null); 99 Residual Null Deviance: Residual Deviance: AIC: residual deviance 89.5 null model 50. Null (null hypothesis) null model

21 ( :11 ) 20 > loglik(fit.null) # fit.null <- glm(y ~ 1,...) log Lik (df=1) null model residual deviance glm()... ( )... Degrees of Freedom: 99 Total (i.e. Null); 98 Residual Null Deviance: Residual Deviance: AIC: Null Deviance residual deviance 85.0 residual deviance 89.5 x i residual deviance ( ) degree of freedom () 51. Degrees of Freedom: 99 Total (i.e. Null); 98 Residual Null model () 99 ( 100-1) 98 (100-2) (k) (log L ) deviance (D = 2 log L ) residual deviance ( 2 log L D) Model k log L Deviance Residual 2 log L deviance NULL f x x + f FULL

22 ( :11 ) 21 table residual deviance 8. (= deviance ) full model R 2 (model selection) AIC (Akaike s information criterion) AIC = 2( ) + 2() = 2 log L + 2k AIC Model k log L Deviance Residual 2 log L AIC deviance NULL f x x + f FULL (Model x) R stepaic() AIC glm()

23 ( :11 ) 22 > library(mass) # stepaic MASS package > stepaic(fit.full) # Start: AIC= y ~ x + f Df Deviance AIC - f <none> x Step: y ~ x AIC= Df Deviance AIC <none> x Call: glm(formula = y ~ x, family = poisson, data = d) Coefficients: (Intercept) x Degrees of Freedom: 99 Total (i.e. Null); 98 Residual Null Deviance: Residual Deviance: AIC: & 9. (GLM) ( ) ANOVA R glm() GLM : link = (poisson) link = log glm() GLM

24 ( :11 ) 23 (prediction) Deviance -2 AIC = Deviance + 2 ( ) logit link GLM

k2 ( :35 ) ( k2) (GLM) web web 1 :

k2 ( :35 ) ( k2) (GLM) web web 1 : 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......................