kubostat2015e p.2 how to specify Poisson regression model, a GLM GLM how to specify model, a GLM GLM logistic probability distribution Poisson distrib

Size: px

Start display at page:

Download "kubostat2015e p.2 how to specify Poisson regression model, a GLM GLM how to specify model, a GLM GLM logistic probability distribution Poisson distrib"

みさえしまむね
5 years ago
Views:

kubostat2015e p.1 I 2015 (e) GLM kubo@ees.hokudai.ac.jp http://goo.gl/76c4i 2015 07 22 2015 07 21 16:26 kubostat2015e (http://goo.

1 kubostat2015e p.1 I 2015 (e) GLM kubo@ees.hokudai.ac.jp :26 kubostat2015e ( (e) / 42 1 N k 2 binomial distribution logit link function 3 4 kubostat2015e ( (e) / 42 statistaical models appeared in the class 6 GLM : The development of linear models Hierarchical Bayesian Model Be more flexible Generalized Linear Mixed Model (GLMM) Incoporating random effects such as individuality parameter estimation MCMC MLE Generalized Linear Model (GLM) Always normal distribution? That's non-sense! MSE Linear model Kubo Doctrine: Learn the evolution of linear-model family, firstly! kubostat2015e ( (e) / 42 kubostat2015e ( (e) / 42? how to specify GLM Generalized Linear Model (GLM) (Poisson regression) (logistic regression) (linear regression) Generalized Linear Model (GLM) probability distribution?? link function? kubostat2015e ( (e) / 42 kubostat2015e ( (e) / 42

kubostat2015e p.2 how to specify Poisson regression model, a GLM GLM how to specify model, a GLM GLM logistic probability distribution Poisson distribution : e.g., β 1 + β 2 x i link function log link function -2 0 2 4 6 0.

2 kubostat2015e p.2 how to specify Poisson regression model, a GLM GLM how to specify model, a GLM GLM logistic probability distribution Poisson distribution : e.g., β 1 + β 2 x i link function log link function probability distribution binomial distribution : e.g., β 1 + β 2 x i link function logit yi x i kubostat2015e ( (e) / 42 kubostat2015e ( (e) / 42 N k N k? 1. N k seeds alive 8 y! y i {0, 1, 2,, 8} f i C: T: i N i = 8 y i = 3 (alive) (dead) x i kubostat2015e ( (e) / 42 kubostat2015e ( (e) / 42 N k Reading data file N k data frame d data4a.csv CSV (comma separated value) format file R > d <- read.csv("data4a.csv") or > d <- read.csv( + " d data frame ( ) > summary(d) N y x f Min. :8 Min. :0.00 Min. : C:50 1st Qu.:8 1st Qu.:3.00 1st Qu.: T:50 Median :8 Median :6.00 Median : Mean :8 Mean :5.08 Mean : rd Qu.:8 3rd Qu.:8.00 3rd Qu.: Max. :8 Max. :8.00 Max. : kubostat2015e ( (e) / 42 kubostat2015e ( (e) / 42

3 kubostat2015e p.3 N k binomial distribution logit link function > plot(d$x,, pch = c(21, 19)[d$f]) > legend("topleft", legend = c("c", "T"), pch = c(21, 19)) yi C T 2. binomial distribution logit link function x i fertilization effective? kubostat2015e ( (e) / 42 kubostat2015e ( (e) / 42 binomial distribution logit link function binomial distribution : N y p(y N, q) = ( ) N q y (1 q) N y y ( N ) y N y p(y i 8, q) q = 0.1 q = 0.3 q = 0.8 y i kubostat2015e ( (e) / 42 binomial distribution logit link function logistic curve (z i: q i = logistic(z i ) = e.g. z i = β 1 + β 2x i) exp( z i ) > logistic <- function(z) 1 / (1 + exp(-z)) # > z <- seq(-6, 6, 0.1) > plot(z, logistic(z), type = "l") q q = 1+exp( z) z kubostat2015e ( (e) / 42 binomial distribution logit link function β 1 and β 2 change logistic curve logit link function binomial distribution logit link function {β 1, β 2 } = {0, 2} (A) β 2 = 2 β 1 (B) β 1 = 0 β 2 q (A) β 2 = 2 β 1 = 2 β 1 = x β 1 = (B) β 1 = 0 β 2 = 4 β 2 = x β 2 = 1 {β 1, β 2 } x q 0 q 1 kubostat2015e ( (e) / 42 logistic 1 q = 1 + exp( (β 1 + β 2 x)) = logistic(β 1 + β 2 x) logit q logit(q) = log 1 q = β 1 + β 2 x logit logistic logistic logit logit is the inverse function of logistic function, vice versa kubostat2015e ( (e) / 42

4 kubostat2015e p.4 binomial distribution logit link function MLE for β 1 and β 2 R β 1 β 2 binomial distribution logit link function (A) f i =C (B) y 7 x 7 > glm(cbind(y, N - y) ~ x + f, data = d, family = binomial)... Coefficients: (Intercept) x ft x yi (A) f i =C x i (B) f i =T x i kubostat2015e ( (e) / 42 kubostat2015e ( (e) / 42? 3. q logit(q) = log 1 q = β 1 + β 2 x + β 3 f + β 4 xf... in case that β 4 < 0, sometimes it predicts... y T C x kubostat2015e ( (e) / 42 kubostat2015e ( (e) / 42 in today s example no interaction effect glm(y ~ x + f,...) glm(y ~ x + f + x:f,...) (A) (B) 4. y T C T C x x little difference kubostat2015e ( (e) / 42 kubostat2015e ( (e) / 42

5 kubostat2015e p.5? How to avoid data/data? / : ? ( ) avoidable data/data values probability : N k indices such as densities use statistical model with binomial distribution : specific leaf area (SLA) use offset term! described later offset! kubostat2015e ( (e) / 42 kubostat2015e ( (e) / 42 unfortunately, sometimes fractions appear... example population densities in research plots offset : hard to avoid... outputs from some measuring machines light intensity index x light index {0.1, 0.2,, 1.0} 10 sometimes we have no choice but plot data/data values... glm(..., family = poisson) kubostat2015e ( (e) / 42 kubostat2015e ( (e) / 42 What? Differences in plot size?!?! R data.frame: Area, light index number of plants x, y x A = / glm() offset > load("d2.rdata") > head(d, 8) # 8 Area x y kubostat2015e ( (e) / 42 kubostat2015e ( (e) / 42

6 kubostat2015e p.6 vs A vs y > plot(d$x, / ) > plot(, ) / d$x kubostat2015e ( (e) / 42 A y kubostat2015e ( (e) / 42 x () x > plot(,, cex = d$x * 2)? kubostat2015e ( (e) / 42 y x kubostat2015e ( (e) / 42 = GLM! 1. i y i λ i : y i Pois(λ i ) 2. λ i A i x i λ i = A i exp(β 1 + β 2 x i ) λ i = exp(β 1 + β 2 x i + log(a i )) log(λ i ) = β 1 + β 2 x i + log(a i ) log(a i ) offset ( β ) family: poisson, link "log" : y ~ x offset : log(area) z = β 1 + β 2 x + log(area) a, b λ log(λ) = z λ = exp(z) = exp(β 1 + β 2 x + log(area)) λ : kubostat2015e ( (e) / 42 kubostat2015e ( (e) / 42

7 kubostat2015e p.7 glm() R glm() > fit <- glm(y ~ x, family = poisson(link = "log"), data = d, offset = log(area)) > print(summary(fit)) Call: glm(formula = y ~ x, family = poisson(link = "log"), data = d, offset = log(area)) (......) Coefficients: Estimate Std. Error z value Pr(> z ) (Intercept) x e-06 kubostat2015e ( (e) / 42 kubostat2015e ( (e) / 42 Plotting the model prediction based on estimation : glm() offset x = 0.9 light environment x = 0.1 dark environment offset = exp( ) solid lines prediction glm() dotted lines true model kubostat2015e ( (e) / 42 kubostat2015e ( (e) / 42 Improve your statisitcal model and remove data/data values! avoidable data/data values probability : N k indices such as densities use statistical model with binomial distribution : specific leaf area (SLA) use offset term! Improve your statistical model! offset kubostat2015e ( (e) / 42 yi (A) The next topic (B) x i = x i y i Generalized Linear Mixed Model (GLMM) kubostat2015e ( (e) / 42

kubostat2017e p.1 I 2017 (e) GLM logistic regression : : :02 1 N y count data or

kubostat2017e p.1 I 2017 (e) GLM logistic regression : : :02 1 N y count data or kubostat207e p. I 207 (e) GLM kubo@ees.hokudai.ac.jp https://goo.gl/z9ycjy 207 4 207 6:02 N y 2 binomial distribution logit link function 3 4! offset kubostat207e (https://goo.gl/z9ycjy) 207 (e) 207 4