k3 ( :07 ) 2 (A) k = 1 (B) k = 7 y x x 1 (k2)?? x y (A) GLM (k

Size: px

Start display at page:

Download "k3 ( :07 ) 2 (A) k = 1 (B) k = 7 y x x 1 (k2)?? x y (A) GLM (k"

かねろうのしろ
5 years ago
Views:

1 k3 ( :07 ) ( k3) kubo@ees.hokudai.ac.jp web web : 4 3 AIC : (1) (2) χ

2 k3 ( :07 ) 2 (A) k = 1 (B) k = 7 y x x 1 (k2)?? x y (A) GLM (k = 1) (B) x 6 GLM (k = 7) GLM 1 (A) log λ = β 1, k = 1 1 (B) x 6 (log λ = β 1 + β 2 x + + β 7 x 6, k = 7) model selection AIC AIC AIC AIC GLM AIC 1 2 R 2 (k2) 100

3 k3 ( :07 ) 3 (A) k = (B) f k = (C) x k = (D) x + f k = x λ k (A) (B) f : (C) x : (D) x + f : y i?? R glm() 2 x i x ; 2 C;?? f i f ; 2 B;?? x + f ; 2 D;?? 3 2 ( ; 2 A 3 ; 2 ) λ exp β 1 β 1 4 maximum log likelihood 3 1 (A)

4 k3 ( :07 ) 4 2 : deviance 4 1 R glm() GLM log L({β j }) log L log L log L D = 2 log L log L -2 5 λ i x i λ i = exp(β 1 + β 2 x i ) x 2 C x log L (D = 2 log L ) log L D 2 log L D D D Null D Null D D glm()... Null Deviance: Residual Deviance: AIC: Null Deviance, Residual Deviance, AIC 3 residual deviance D ( ) R full model deviance -2 χ 2

5 k3 ( :07 ) x Deviance 2 log L () 89.5 (Null Deviance) 85.0 (Residual Deviance) deviance deviance (null deviance) (residual deviance) y i y i = {6, 6, 6, 12, 10, } 100 i {1, 2, 3} y i 6 {λ 1, λ 2, λ 3 } = {6, 6, 6} i = 4 y 4 12 λ 4 = 12 i = 5 y 5 10 λ 5 = log L 6 > sum(log(dpois(d$y, lambda = d$y))) [1] D = 2 log L = x D ( D) = = 85.0 glm() Residual Deviance: AIC: Residual Deviance R dpois() 100 {y i } {λ i } = {y 1, y 2, y 3, } i log L = 0

6 k3 ( :07 ) λ i = exp(β 1 ) β 1 k = 1 R null model 8 log λ i = β 1 2 A R glm() glm(y ~ 1,...) > fit.null <- glm(formula = y ~ 1, family = poisson, data = d) fit.null β Degrees of Freedom: 99 Total (i.e. Null); Null Deviance: Residual Deviance: AIC: Residual 89.5 > loglik(fit.null) log Lik (df=1) D k log L D = 2 log L 2 log L D 2 k 10 3 AIC GLM 8 null hypothesis 9 1 (A) 10 f 2 11??

7 k3 ( :07 ) 7 2 k log L Deviance Residual deviance 2 k log L Deviance Residual 2 log L deviance f x x + f AIC k log L Deviance Residual 2 log L deviance AIC f x x + f model selection model selection criterion AIC (Akaike s information criterion) AIC goodness of fitgoodness of prediction 12 k AIC AIC = 2 { } = 2(log L k) = D + 2k AIC 2 AIC 3 x AIC 13 AIC 4 statistical test 12 13?? GLM R stepaic()

8 k3 ( :07 ) 8 14 AIC likelihood ratio test parametric null hypothesis alternative hypothesis 14 15?? most powerful test AIC

9 k3 ( :07 ) 9 AIC ( ) ( ) 4 AIC 4 test statistic 22 95% 5% significant level 23 Neyman-Pearson 24 6 : 5 λ = exp(β 1 + β 2 x i ) GLM x : λ i x i β 2 = 0; k = 1 x : λ i x i β 2 0; k = Neyman-Pearson

10 k3 ( :07 ) 10 (A) (B) i y i x i yi x = = x i 5 (A) f i (B) 100 x x i ; x x i 4 x AIC 3 k log L Deviance Residual 2 log L deviance AIC x x likelihood ratio : L 1 L 2 = : exp( 237.6) x : exp( 235.4) -2 D 1,2 = 2 (log L 1 log L 2) 27 D 1 = 2 log L 1 D 2 = 2 log L 2 D 1,2 = D 1 D 2 D 1,2 x x D 1,2 = 4.5 x (A) (C) 3-2 D 1,2 χ 2 8.2

11 k3 ( :07 ) 11 5 D 1,2 ( ) ( ) () () 7 5 Neyman-Pearson : k = 1, β 2 = 0 : x k = 2, β Neyman-Pearson 9 5 : D 1,2 = 4.5 x β 2 0 type I error : x D 1,2 = 4.5 x β 2 = 0 type II error Neyman-Pearson 1 2 ˆβ 1 = 2.06 (p.6 ) 3 β 2 = 0(k = 1) β 2 0(k = 2) D 1,2 D 1, x alternative hypothesis

12 k3 ( :07 ) 12 ( ˆβ 1 = 2.06 ) x D 1,2 D 1,2 D 1,2 D 1, x D 1,2 ˆβ 1 = 2.06, p.6 D 1,2 4 x D 1,2 4.5 P D 1,2 = x D 1,2 4.5 P P P value P P : D 1,2 = 4.5 P : D 1,2 = 4.5 x! P Neyman-Pearson α 30 : P α : P < α : α α = P α

13 k3 ( :07 ) (1) P D 1,2 4.5 P P parametric bootstrap 32 χ 2 (PB) 6 R glm() x fit1 fit2 fit1 fit2 > fit2$deviance [1] x x D 1,2 > fit1$deviance - fit2$deviance [1] D 1, rpois() 100 > d$y.rnd <- rpois(100, lambda = mean(d$y)) mean(d$y) 7.85 glm() x glm() ˆβ 1 = 2.06 exp(2.06) = mean(d$y)

14 k3 ( :07 ) 14 > fit1 <- glm(y.rnd ~ 1, data = d, family = poisson) > fit2 <- glm(y.rnd ~ x, data = d, family = poisson) > fit1$deviance - fit2$deviance [1] x i 1.92 x : 1 mean(d$y) d$y.rnd 2 d$y.rnd,x glm() fit1, fit2 3 fit1$deviance - fit2$deviance PB D 1,2 34 PB R pb() 35 get.dd <- function(d) # { n.sample <- nrow(d) # y.mean <- mean(d$y) # d$y.rnd <- rpois(n.sample, lambda = y.mean) fit1 <- glm(y.rnd ~ 1, data = d, family = poisson) fit2 <- glm(y.rnd ~ x, data = d, family = poisson) fit1$deviance - fit2$deviance # } pb <- function(d, n.bootstrap) { sapply(1:n.bootstrap, get.dd, d) } pb.r 36 R R pb.r pb() bootstrap method fit1 fit2$null.deviance - fit2$deviance D 1,2 web ( )

15 k3 ( :07 ) D 1,2 = x D 1,2 7 D 1,2 D 1, x D 1,2 = 4.5 > source("pb.r") # pb.r > dd12 <- pb(d, n.bootstrap = 1000) R D 1, dd12 summary() > summary(dd12) Min. 1st Qu. Median Mean 3rd Qu. Max e e e e e e+01 7 D 1,2 4.5 > hist(dd12, 100) > abline(v = 4.5, lty = 2) 1000 D 1,2 4.5 > sum(dd12 >= 4.5) [1] / 1000 P = P = 0.05 D 1, D 1, n.bootstrap P = α D 1,2 (critical point) D 1,2 (critical region rejection region)

16 > quantile(dd12, 0.95) 95% k3 ( :07 ) 16 5% D 1, P significantly different 40 x 8.2 (2) χ 2 PB 7 41 fit1 fit2 x > fit1 <- glm(y ~ 1, data = d, family = poisson) > fit2 <- glm(y ~ x, data = d, family = poisson) anova() 42 > anova(fit1, fit2, test = "Chisq") Analysis of Deviance Table Model 1: y ~ 1 Model 2: y ~ x Resid. Df Resid. Dev Df Deviance P(> Chi ) D 1, χ 2 χ 2 distribution 39?? 40 P Neyman-Pearson P < α anova() ANOVA analysis of variance analysis of deviance 43 x

17 k3 ( :07 ) 17 "Chisq" χ 2 D 1,2 4.5 P P PB P = χ "Chisq" P PB 44 χ 2 t F 9 α = 0.05 D 1,2 P < α P α fail to reject Neyman-Pearson 45 Neyman-Pearson 7 P < α P α 5 P 2 46 Neyman-Pearson P P 2 P 2 1 P 2 ; power PB β β 1 P 2 47

18 k3 ( :07 ) AIC AIC Neyman-Pearson AIC P < α AIC P effect size

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

講義のーと : データ解析のための統計モデリング. 第５回 Title 講義のーと : データ解析のための統計モデリング Author(s) 久保, 拓弥 Issue Date 2008 Doc URL http://hdl.handle.net/2115/49477 Type learningobject Note この講義資料は, 著者のホームページ http://hosho.ees.hokudai.ac.jp/~kub ードできます Note(URL)http://hosho.ees.hokudai.ac.jp/~kubo/ce/EesLecture20