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

Size: px

Start display at page:

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

ことこつつの
8 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 kubostat2008e.pdf ( 第 5 回 ) Instructions for use Hokkaido University Collection of Scholarly and Aca

2 ( :11 ) 1 ( ) 5 (+2) 5 ( ) [email protected] Deviance (1) parametric bootstrap (2) χ ? : 18 5 ( ) (GLM)

3 ( :11 ) 2 GLM () link ; (deviance) ( ) (AIC )? : : ( ) : discussion : : ( ) : : () / / 7

4 ( :11 ) 3? ( ) 1. 3 GLM () (data3a.csv ) f i 1 i y i 1. 3 f i y i x i 1: i (i = 1, 2,, 100) CSV read.csv() 2 > d <- read.csv("data3a.csv") > head(d) y x f C C C C C C f 3. 3

5 ( :11 ) 4 > (fit2 <- glm(y ~ x, data = d, family = poisson))...( )... Coefficients: (Intercept) x ( )... Residual Deviance: AIC: i λ i = exp(β 1 + β 2 x i ) ( 2 2 ) λ i = exp(β 1 ) ( 1 1 ) > (fit1 <- glm(y ~ 1, data = d, family = poisson))...( )... Coefficients: (Intercept) ( )... Residual Deviance: AIC: > plot(d$x, d$y) # ( ) > xx <- seq(min(d$x), max(d$x), length = 50) # x > abline(h = mean(d$y), lty = 2, lwd = 2) # model 1 > lines(xx, exp( * xx)) # model 2 d$y d$x 1 ( x i ) 2 ( x i ) exp(2.058) = 1 2 mean(d$y) = 7.83

6 ( :11 ) 5 Model k log L Deviance Residual 2 log L AIC deviance Model Model FULL AIC (Akaike s information criterion) AIC = 2 ( ) + 2 ( ) = (Deviance) + 2 ( ) = 2 log L + 2k deviance AIC AIC c AIC c 2. ( (deviance) ( ) ) 7? 6. 2 log D residual deviance () () ()

7 ( :11 ) 6 1 ( 1 x i ) 2 ( 2 x i ) ( 1 2) : 1 2 (AIC ) (!) :? 1 vs x i ( i λ i x i ) 9. x i 1 ( 1) (null hypothesis) 2 ( 2) (alternative hypothesis) 1 ( ) 2 deviance ( )

8 ( :11 ) X deviance 10 deviance AIC deviance 2 (penalty) Deviance? 1 2 (deviance ) deviance 11 deviance (likelihood ratio test) 11. : 1 ( ) vs 2 ( ) 12 (Neyman-Pearson ) 12. Deviance {log( 1 ) log( 2 )} deviance deviance

9 ( :11 ) 8 4. deviance 1 ( ) 2 ( ) deviance 4.5 deviance A: () 2 B: A B 1,2 deviance ( ) () OK OK 1 2 ( B) 1 ( A) : deviance (type I error) 2 : deviance (type II error)

10 ( :11 ) 9? ( ) x i! 5.? 1 vs deviance deviance deviance deviance () 5. deviance 4.5 P P (P value) 17 P P : deviance 4.5 P : deviance 4.5 2! 17. P value?

11 ( :11 ) 10 P α : P < α : 1 x i 2 α α? % α = 0.05 α Neyman-Pearson Neyman-Pearson 1 2 deviance 4.5? 1 deviance 4.5 () P? 6. (1) parametric bootstrap 1 2 deviance P parametric bootstrap 2. (deviance χ 2 ) parametric bootstrap deviance R glm() fit1 fit2? 2 fit2

12 ( :11 ) 11 > (fit2 <- glm(y ~ x, data = d, family = poisson))...( )... Residual Deviance: AIC: (residual) deviance 20 fit2 > sort(names(fit2)) [1] "R" "aic" "boundary" [4] "call" "coefficients" "contrasts" [7] "control" "converged" "data" [10] "deviance" "df.null" "df.residual"... ( )... > fit2$deviance [1] residual deviance 1 2 deviance > fit1$deviance - fit2$deviance [1] deviance deviance rpois() 1 > d$y.rnd <- rpois(100, lambda = mean(d$y)) fit1

13 ( :11 ) glm() 1 2 > 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 1 deviance ( ) 2 2 deviance ( ) 21. glm() coefficient exp(2.05) deviance mean(d$y)a d$y.rnd 2. d$y.rnd 1, 2 GLM fit1, fit2 3. deviance fit1$deviance - fit2$deviance 1 deviance parametric bootstrap bootstrap deviance parametric bootstrap pb() pb() pb.r

14 ( :11 ) 13 pb <- function(d, n.bootstrap) { n.sample <- nrow(d) # y.mean <- mean(d$y) # v.d.dev12 <- sapply( # PB deviance 1:n.bootstrap, function(i) { 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 # deviance } ) v.d.dev12 # deviance vector } # : fit1 # fit2$null.deviance - fit2$deviance # deviance ( fit1, fit2 ) web page pb.r pb.r R 24 source("pb.r")24. R R pb() R working directory bootstrap 1000 > source("pb.r") > diff.dev12 <- pb(d, n.bootstrap = 1000) #......( )... GLM deviance 1000 diff.dev12 vector > summary(diff.dev12) Min. 1st Qu. Median Mean 3rd Qu. Max e e e e e e deviance deviance 4.51

15 ( :11 ) 14 > hist(diff.dev12, 100) > abline(v = 4.51, lty = 2) Histogram of diff.dev12 Frequency diff.dev deviance 4.51? > sum(diff.dev12 >= 4.51) [1] deviance / 1000 P = P = 0.05 deviance > quantile(diff.dev12, 0.95) 95% deviance (significantly different) 27 ( 1 1) 2 parametric bootstrap 28 P 25. 5% significant 28. deviance

16 ( :11 ) 15 parametric bootstrap ; R 7. (2) χ 2 parametric bootstrap deviance R : anova() ANOVA > fit1 <- glm(y ~ 1, data = d, family = poisson) # 1 > fit2 <- glm(y ~ x, data = d, family = poisson) # 2 analysis of variance deviance ( > anova(fit1, fit2, test = "Chisq") # ) Analysis of Deviance Table ANOVA analysis of deviance Model 1: y ~ 1 Model 2: y ~ x Resid. Df Resid. Dev Df Deviance P(> Chi ) parametric bootstrap 1 2 deviance ( 1 ) pb() diff.dev12 vector (diff.dev12 ) deviance 1 31 χ 2 32 χ ( "Chisq" ) deviance 4.51 P

17 ( :11 ) 16 parametric bootstrap deviance ? 1 ( ) 2 ( ) 1 P 1,2 deviance ( ) () OK OK = 1 ( ) / 7 α P < α P < α : 1 x i α P α? P α : ! P < α P α P α

18 ( :11 ) 17 β () β Neyman-Pearson P 36 P (?) P < α : deviance ( ) χ 2 χ 2 table P = 0.05 deviance? ? ( )?

19 ( :11 ) 18 : 40 P P 0.05 P 40. ( ) ( ; power) 41 : β (experimental design) x i () x i

20 ( :11 ) : GLM GLM / random effects random effects (generalized linear mixed model; GLMM) GLMM GLMM? ( )? ( ) random effects

21 ( :11 ) 20 Bayes Bayes MCMC random effects (GLMM) (GLM) +, fixed effects / 7

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

k3 ( :07 ) 2 (A) k = 1 (B) k = 7 y x x 1 (k2)?? x y (A) GLM (k 2012 11 01 k3 (2012-10-24 14:07 ) 1 6 3 (2012 11 01 k3) [email protected] web http://goo.gl/wijx2 web http://goo.gl/ufq2 1 3 2 : 4 3 AIC 6 4 7 5 8 6 : 9 7 11 8 12 8.1 (1)........ 13 8.2 (2) χ 2....................