#-------------------------------------------------------------------
# 9.1  CHOOSING THE FAMILY
#-------------------------------------------------------------------

> snow <- c(4.1, 12.2, 0, 14.6, 2.5, 2.1, 4.1, 11.5, 0.7, 0, 23.3, 6.3, 
        0, 0, 0, 3.1, 2.6, 1, 0.3, 0, 0, 0, 2.3, 11.3, 6.5, 0, 4.8, 
        1.1, 0, 17.4, 0, 10.9, 1.3, 0, 0.8, 0, 1, 18.4, 0, 0, 0, 3.8, 
        0, 7, 31, 0, 6.6, 0, 16.9, 0, 0.3, 10.1, 0, 5, 1)
> year <- 6:60
> par(mfrow = c(2, 2))
> plot(year, snow, yaxs = "s", main = "Figure 9.1 Snow")

> poissonexp <- list(family = c(name = "poissonexp", link = "log", 
        variance = "mu^(3/2)"), names = "log", link = function(mu)
log(mu), inverse = function(eta)
care.exp(eta), deriv = function(mu)
1/mu, initialize = expression({
        y <- as.numeric(y)
        mu <- y + 0.167 * (y == 0)
}
), variance = function(mu)
mu^(3/2), deviance = function(mu, y, w, residuals = F)
if(residuals) (sqrt(w) * (sqrt(y) - sqrt(mu)) * 2)/sqrt(sqrt(mu))
         else sum((w * (sqrt(y) - sqrt(mu))^2 * 4)/sqrt(mu)), weight
         = expression(w/((sqrt(family$variance(mu)) * family$deriv(mu)
        )^2)))
> class(poissonexp) <- "family" 

> summary(glm(snow ~ 1, family = poissonexp))

Coefficients:
               Value Std. Error  t value 
(Intercept) 1.497592   0.204564 7.320897

(Dispersion Parameter for poissonexp family taken to be 4.866524 )

Residual Deviance: 289.171 on 54 degrees of freedom

Number of Fisher Scoring Iterations: 0 
> summary(glm(snow ~ year, family = poissonexp))

Coefficients:
                     Value Std. Error      t value 
(Intercept)  1.50028567814 0.47559937  3.154515697
       year -0.00008162735 0.01299116 -0.006283302

(Dispersion Parameter for poissonexp family taken to be 4.949669 )

Residual Deviance: 289.1708 on 53 degrees of freedom

Number of Fisher Scoring Iterations: 3 
> theta <- (2 * pi * year)/11
> sine <- sin(theta)
> cosine <- cos(theta)
> glm <- glm(snow ~ sine + cosine, family = poissonexp)
> summary(glm)

Coefficients:
                 Value Std. Error    t value 
(Intercept)  1.4769288  0.2040109  7.2394583
       sine -0.1615921  0.2876340 -0.5617977
     cosine -0.2376885  0.2879023 -0.8255876

(Dispersion Parameter for poissonexp family taken to be 4.769283 )

Residual Deviance: 284.3678 on 52 degrees of freedom

Number of Fisher Scoring Iterations: 3 
> lines(year, fitted(glm))

#-------------------------------------------------------------------
# 9.1  CHOOSING THE LINK
#-------------------------------------------------------------------

> m <- matrix(scan("/glim/datasets/trees"), ncol = 3, byrow = T)
> diam <- m[, 1]
> height <- m[, 2]
> vol <- m[, 3]
> glm(vol ~ diam + height)

Coefficients:
 (Intercept)     diam    height 
   -57.98766 4.708161 0.3392512

Degrees of Freedom: 31 Total; 28 Residual
Residual Deviance: 421.9214 
> 4.708161/0.3392512
[1] 13.8781

> glm(vol ~ diam + height, family = quasi(variance = constant, link = 
        log))

Coefficients:
 (Intercept)      diam    height 
   0.6794189 0.1341617 0.0111431

Degrees of Freedom: 31 Total; 28 Residual
Residual Deviance: 272.5712 
> 0.1341617/0.0111431
[1] 12.03989

> glm(vol ~ diam + height, family = quasi(variance = constant, link = 
        sqrt))

Coefficients:
 (Intercept)      diam     height 
    -3.10926 0.4106364 0.03913296

Degrees of Freedom: 31 Total; 28 Residual
Residual Deviance: 185.729 
> 0.4106364/0.03913296
[1] 10.49336

> vol
 [1] 10.3 10.3 10.2 16.4 18.8 19.7 15.6 18.2 22.6 19.9 24.2 21.0 21.4
[14] 21.3 19.1 22.2 33.8 27.4 25.7 24.9 34.5 31.7 36.3 38.3 42.6 55.4
[27] 55.7 58.3 51.5 51.0 77.0
> fvol <- floor(vol/5) - 1
> fvol
 [1]  1  1  1  2  2  2  2  2  3  2  3  3  3  3  2  3  5  4  4  3  5  5
[23]  6  6  7 10 10 10  9  9 14
> plot(height, diam, type = "n", main = "Figure 9.2 Trees")
> for(i in 1:14)
        text(height[fvol == i], diam[fvol == i], i)
> Fvol <- factor(fvol)
> glm <- glm(vol ~ diam + height + Fvol)
> glm

Coefficients:
 (Intercept)     diam    height    Fvol1    Fvol2     Fvol3    Fvol4 
    1.625503 1.496181 0.1813517 1.304378 1.480346 0.9626389 1.783727
    Fvol5    Fvol6    Fvol7    Fvol8    Fvol9 
 1.705551 1.647337 1.930808 2.133257 3.197286

Degrees of Freedom: 31 Total; 19 Residual
Residual Deviance: 24.48716 
> cov <- summary(glm)$cov.unscaled[2:3, 2:3]
> cov
             diam      height 
  diam 0.14294736 0.010633641
height 0.01063364 0.002602849
> evec <- eigen(cov)
> evec
$values:
[1] 0.143748475 0.001801731

$vectors:
            [,1]        [,2] 
[1,] -0.99717411  0.07512518
[2,] -0.07512518 -0.99717411

> -0.99717411/-0.07512518
[1] 13.2735
> eta <- 13.2735 * diam + height
> plot(vol, eta, main = "Figure 9.3 Estimated link function")