GeneralizedLinearModels LinearModel u response variable intercept slope explanatory

Generalized　Linear　Models • Linear　Model u　 　response variable ～　intercept + slope * explanatory variable u　lm(y~ x + f ・・・)，lm(y~x + f -1) (no intercept) require(graphics) ## Annette Dobson (1990) "An Introduction to Generalized Linear Models". ## Page 9: Plant Weight Data. ctl <- c(4. 17, 5. 58, 5. 18, 6. 11, 4. 50, 4. 61, 5. 17, 4. 53, 5. 33, 5. 14) trt <- c(4. 81, 4. 17, 4. 41, 3. 59, 5. 87, 3. 83, 6. 03, 4. 89, 4. 32, 4. 69) group <- gl(2, 10, 20, labels=c("Ctl", "Trt")) weight <- c(ctl, trt) lm. D 9 <- lm(weight ~ group) lm. D 90 <- lm(weight ~ group - 1) 　# omitting intercept anova(lm. D 9) summary(lm. D 90) opar <- par(mfrow = c(2, 2), oma = c(0, 0, 1. 1, 0)) plot(lm. D 9, las = 1) # Residuals, Fitted, . . . Par(opar)

Generalized　Linear　Models • Linear　Model u　 　response variable ～　intercept + slope * explanatory variable u　 lm(y~ x + f ・・・)，lm(y~x + f -1) (no intercept) u. Generalized Linear Model u Model &Link function ～　intercept + slope * explanatory variable u glm(y ~ x, data = d, family = poisson)

Poisson　Model （counting data of occurrence) u. Poisson　Model u λ：mean occurrence in unit time u. Identity link u. Log link(canonical) u正準リンク関数：最も自然なリンク関数：乗法効果） Link function ～　intercept + slope * explanatory variable u glm(y ~ x, data = d, family = poisson(link=“log”)) u. Canonical link function is set as default

Poisson　Model (p 49) （counting data of occurrence) u. Poisson Model for number of seeds of a plant, regressed on plant size and nutrification (p 49) u. Maximize log-likelihood glm(y ~ x + f, data = d, family = poisson) #page 42 plant data d <- read. csv("data 3 a. csv") d$y # number of seeds d$x # plant size (hight) d$f # nutrification (treat-control) plot(d$x, d$y, pch =c(21, 19)[d$f]) # model p 58 fit. all <- glm(y ~ x + f, data=d, family=poisson) print(fit. all) log. Lik(fit. all) 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. 263 + 0. 0801 * xx), lwd=2)

Poisson　Model (p 49) （counting data of occurrence) u. Poisson Model for number of seeds of a plant, regressed on plant size and nutrification (p 49) u. Maximize log-likelihood #page 42 plant data d <- read. csv("data 3 a. csv") d$y # number of seeds d$x # plant size (hight) d$f # nutrification (treat-control) plot(d$x, d$y, pch =c(21, 19)[d$f]) # model p 58 fit. all <- glm(y ~ x + f, data=d, family=poisson) print(fit. all) log. Lik(fit. all) 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. 263 + 0. 0801 * xx), lwd=2)

Other Generalized Linear Models (chap 6 p 114) (discrete) (continuous) Probability Random numbers generation Family in glm() Standard link function Binomial rbinom() binomial logit Poisson rpois() poisson log Negative Binomial rnbinom() (glb. nb() function) log Gamma rgamma() gamma log, inverse Normal rnorm() gaussian identity

Generalized　Linear　Models u. Generalized Linear Model glm(y ~ x, data = d, family = poisson) u. Family （Modelled Probability Distribution) ubinomial(link = “logit“) 2項分布（規定試行中の発生数） ugaussian(link = “identity”)　正規分布 u. Gamma(link = “inverse”)　　ガンマ分布（正のみ） uinverse. gaussian(link = “ 1/mu^2”)　逆ガウス分布 upoisson(link = “log”)　ポアソン分布（一定時間中の発生回数） uquasi(link = “identity”, variance = “constant”)正規分布（不均一） u quasibinomial(link = “logit”)　2項分布（分散不均一） u quasipoisson(link = “log”) 　　ポアソン分布（分散不均一）

Binomial Logistic Model (p 118) （occurrence number in given trials) u. Binomial Model for the number of survived plant in 8 obserbations, regressed on plant size and nutrification (p 118) u. Maximize log-likelihood glm(cbind(y, N-y) ~ x + f, data = d, family = binomial) #page 117 plant data d <- read. csv("data 4 a. csv") d$N # number of trials d$y # number of survived plant d$x # plant size d$f # nutrification (treat-control) # model p 122 fit. all <- glm(cbind(y, N-y) ~ x + f, data=d, family=binomi print(fit. all) log. Lik(fit. all)

Offset Term(p 131) （avoid a division calculation) u. Count data for several zones having different area, or different population u. One way is define a density (occurrence in unit area) and apply Poisson model glm(y ~ x, offset =log(A), data = d, family = poisson) #page 133 plant data d <- read. csv("data 4 b. csv") d$y # number of plants in lot i d$x # brightness at lot I d$A # area of lot i plot(d$A, d$y) # model p 131 fit<- glm(y ~ x, offset = log(A) , data=d, family=poisson) print(fit) log. Lik(fit)

Gamma Distribution Model (p 138) u. Gamma Distribution (continuous positive data) u s: shape parameter, r: rate parameter, theta=1/r: scale parameter u time length before s times occurrence of random events with occurrence rate of r. (average occurrence interval is ) u Average : Variance: u dgamma(y, shape, rate) u. Weight of flower of a plant y (continuous, positive) uaverage weight u. Loglink function of linear estimator uglm(y ~ log(x), data = d, family = gamma(link="log"))

Gamma Distribution Model (p 138) u. Gamma Distribution (continuous positive data) uglm(y ~ log(x), data = d, family = gamma(link="log") # A Gamma example, from Mc. Cullagh & Nelder (1989, pp. 300 -2) clotting <- data. frame( u = c(5, 10, 15, 20, 30, 40, 60, 80, 100), lot 1 = c(118, 58, 42, 35, 27, 25, 21, 19, 18), lot 2 = c(69, 35, 26, 21, 18, 16, 13, 12)) summary(glm(lot 1 ~ log(u), data=clotting, family=Gamma)) summary(glm(lot 2 ~ log(u), data=clotting, family=Gamma)) Call: glm(formula = lot 1 ~ log(u), family = Gamma, data = clotting) Deviance Residuals: Min 1 Q Median 3 Q Max -0. 04008 -0. 03756 -0. 02637 0. 02905 0. 08641 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -0. 0165544 0. 0009275 -17. 85 4. 28 e-07 *** log(u) 0. 0153431 0. 0004150 36. 98 2. 75 e-09 *** ---