Chap 4 - Exercise 4.2

> setwd("C:/Cauldron/garage/R/soulcraft/Volatility/Learn/Dobson-GLM")
> data <- read.csv("test9.csv", header = T, stringsAsFactors = F)
> par(mfrow = c(1, 1))
> plot(data$xi, data$yi, pch = 19, col = "blue", cex = 1)

> model.pois <- glm(data$yi ~ data$xi, family = Gamma(link = log))
> summary(model.pois)
Call:
glm(formula = data$yi ~ data$xi, family = Gamma(link = log))

Deviance Residuals:
     Min        1Q    Median        3Q       Max
-0.18204  -0.07676   0.02860   0.08475   0.21735

Coefficients:
              Estimate Std. Error t value Pr(>|t|)
(Intercept)  1.5278401  0.0436363  35.013 8.45e-16 ***
data$xi     -0.0019717  0.0005338  -3.694  0.00217 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for Gamma family taken to be 0.01346761)

    Null deviance: 0.38385  on 16  degrees of freedom
Residual deviance: 0.20315  on 15  degrees of freedom
AIC: 26.574

Number of Fisher Scoring iterations: 4
> qqnorm(resid(model.pois))

> beta0 <- c(3, 1)
> n <- dim(data)[1]
> poisreg <- function(beta0) {
+     beta <- matrix(data = beta0, ncol = 1)
+     X <- cbind(rep(1, n), data$xi)
+     W <- diag(n)
+     LHS <- t(X) %*% W %*% X
+     z <- matrix(data = NA, nrow = n, ncol = 1)
+     z <- X %*% beta + (data$yi - exp(X %*% beta))/exp(X %*% beta)
+     RHS <- t(X) %*% W %*% z
+     beta.res <- solve(LHS, RHS)
+     return(beta.res)
+ }
> iterations <- matrix(data = NA, nrow = 2, ncol = 150)
> iterations[, 1] <- poisreg(beta0)
> for (i in 2:150) {
+     iterations[, (i)] <- poisreg(iterations[, (i - 1)])
+ }
> print(iterations[, 150])
[1]  1.527840097 -0.001971691

> W <- diag(n)
> X <- cbind(rep(1, n), data$xi)
> J <- t(X) %*% W %*% X
> Jinv <- solve(J)

> sqrt(Jinv[1, 1])
[1] 0.3760131
> sqrt(Jinv[2, 2])
[1] 0.00459955
> summary(model.pois)
Call:
glm(formula = data$yi ~ data$xi, family = Gamma(link = log))

Deviance Residuals:
     Min        1Q    Median        3Q       Max
-0.18204  -0.07676   0.02860   0.08475   0.21735

Coefficients:
              Estimate Std. Error t value Pr(>|t|)
(Intercept)  1.5278401  0.0436363  35.013 8.45e-16 ***
data$xi     -0.0019717  0.0005338  -3.694  0.00217 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for Gamma family taken to be 0.01346761)

    Null deviance: 0.38385  on 16  degrees of freedom
Residual deviance: 0.20315  on 15  degrees of freedom
AIC: 26.574

Number of Fisher Scoring iterations: 4