Interpreting Generalized Linear Models

Response residuals

For type = "response", the conventional residual on the response level is computed, that is, \[r_i = y_i – \hat{f}(x_i)\,.\] This means that the fitted residuals are transformed by taking the inverse of the link function:

# type = "response"
res.response1 <- residuals(model.pois, type = "response")
res.response2 <- expected - estimates
all.equal(res.response1, res.response2)
## [1] TRUE

Working residuals

For type = "working", the residuals are normalized by the estimates \(\hat{f}(x_i)\):

\[r_i = \frac{y_i – \hat{f}(x_i)}{\hat{f}(x_i)}\,.\]

# type = "working"
res.working1 <- residuals(model.pois, type="working")
res.working2 <- (expected - estimates) / estimates
all.equal(res.working1, res.working2)
## [1] TRUE

Pearson residuals

For type = "pearson", the Pearson residuals are computed. They are obtained by normalizing the residuals by the square root of the estimate:

\[r_i = \frac{y_i – \hat{f}(x_i)}{\sqrt{\hat{f}(x_i)}}\,.\]

# type = "pearson"
res.pearson1 <- residuals(model.pois, type="pearson")
res.pearson2 <- (expected - estimates) / sqrt(estimates)
all.equal(res.pearson1, res.pearson2)
## [1] TRUE

Deviance residuals

Deviance residuals are defined by the deviance. The deviance of a model is given by

\[{D(y,{\hat {\mu }})=2{\Big (}\log {\big (}p(y\mid {\hat {\theta }}_{s}){\big )}-\log {\big (}p(y\mid {\hat {\theta }}_{0}){\big )}{\Big )}.\,}\]

where

• \(y\) is the outcome
• \(\hat{\mu}\) is the estimate of the model
• \(\hat{\theta}_s\) and \(\hat{\theta}_0\) are the parameters of the fitted saturated and proposed models, respectively. A saturated model has as many parameters as it has training points, that is, \(p = n\). Thus, it has a perfect fit. The proposed model can be the any other model.
• \(p(y | \theta)\) is the likelihood of data given the model

The deviance indicates the extent to which the likelihood of the saturated model exceeds the likelihood of the proposed model. If the proposed model has a good fit, the deviance will be small. If the proposed model has a bad fit, the deviance will be high. For example, for the Poisson model, the deviance is

\[D = 2 \cdot \sum_{i = 1}^n y_i \cdot \log \left(\frac{y_i}{\hat{\mu}_i}\right) − (y_i − \hat{\mu}_i)\,.\]

In R, the deviance residuals represent the contributions of individual samples to the deviance \(D\). More specifically, they are defined as the signed square roots of the unit deviances. Thus, the deviance residuals are analogous to the conventional residuals: when they are squared, we obtain the sum of squares that we use for assessing the fit of the model. However, while the sum of squares is the residual sum of squares for linear models, for GLMs, this is the deviance.

How does such a deviance look like in practice? For example, for the Poisson distribution, the deviance residuals are defined as:

\[r_i = \text{sgn}(y – \hat{\mu}_i) \cdot \sqrt{2 \cdot y_i \cdot \log \left(\frac{y_i}{\hat{\mu}_i}\right) − (y_i − \hat{\mu}_i)}\,.\]

Let us verify this in R:

# type = "deviance"
res.dev1 <- residuals(model.pois, type = "deviance")
res.dev2 <- residuals(model.pois)
poisson.dev <- function (y, mu) 
    # unit deviance
    2 * (y * log(ifelse(y == 0, 1, y/mu)) - (y - mu))
res.dev3 <- sqrt(poisson.dev(expected, estimates)) * 
        ifelse(expected > estimates, 1, -1)
all.equal(res.dev1, res.dev2, res.dev3)
## [1] TRUE

Note that, for ordinary least-squares models, the deviance residual is identical to the conventional residual.

Deviance residuals in practice

We can obtain the deviance residuals of our model using the residuals function:

summary(residuals(model.pois))
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -5.3137 -2.0913 -0.6163 -0.1642  1.7141  6.0361

Since the median deviance residual is close to zero, this means that our model is not biased in one direction (i.e. the out come is neither over- nor underestimated).

Null deviance and residual deviance in practice

Let us investigate the null and residual deviance of our model:

paste0(c("Null deviance: ", "Residual deviance: "),
       round(c(model.pois$null.deviance, deviance(model.pois)), 2))
## [1] "Null deviance: 1464.87"    "Residual deviance: 453.25"

These results are somehow reassuring. First, the null deviance is high, which means it makes sense to use more than a single parameter for fitting the model. Second, the residual deviance is relatively low, which indicates that the log likelihood of our model is close to the log likelihood of the saturated model.

However, for a well-fitting model, the residual deviance should be close to the degrees of freedom (74), which is not the case here. For example, this could be a result of overdispersion where the variation is greater than predicted by the model. This can happen for a Poisson model when the actual variance exceeds the assumed mean of \(\mu = Var(Y)\).

Dispersion parameter

Dispersion (variability/scatter/spread) simply indicates whether a distribution is wide or narrow. The GLM function can use a dispersion parameter to model the variability.

However, for likelihood-based model, the dispersion parameter is always fixed to 1. It is adjusted only for methods that are based on quasi-likelihood estimation such as when family = "quasipoisson" or family = "quasibinomial". These methods are particularly suited for dealing with overdispersion.

AIC

The Akaike information criterion (AIC) is an information-theoretic measure that describes the quality of a model. It is defined as

\[\text{AIC} = 2p – 2 \ln(\hat{L})\]

where \(p\) is the number of model parameters and \(\hat{L}\) is the maximum of the likelihood function. A model with a low AIC is characterized by low complexity (minimizes \(p\)) and a good fit (maximizes \(\hat{L}\)).

Fisher scoring iterations

The information about Fisher scoring iterations is just verbose output of iterative weighted least squares. A high number of iterations may be a cause for concern indicating that the algorithm is not converging properly.

The type argument

Since models obtained via lm do not use a linker function, the predictions from predict.lm are always on the scale of the outcome (except if you have transformed the outcome earlier). For predict.glm this is not generally true. Here, the type parameter determines the scale on which the estimates are returned. The following two settings are important:

• type = "link": the default setting returns the estimates on the scale of the link function. For example, for Poisson regression, the estimates would represent the logarithms of the outcomes. Given the estimates on the link scale, you can transform them to the estimates on the response scale by taking the inverse link function.
• type = "response": returns estimates on the level of the outcomes. This is the option you need if you want to evaluate predictive performance.

Let us see how the returned estimates differ depending on the type argument:

# prediction on link scale (log)
pred.l <- predict(model.pois, newdata = ozone[testset, ])
summary(pred.l)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.688   3.217   3.744   3.695   4.220   4.864
# prediction on respone scale
pred.r <- predict(model.pois, newdata = ozone[testset, ], type = "response")
summary(pred.r)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    5.41   24.95   42.25   50.12   68.01  129.61

Using the link and inverse link functions, we can transform the estimates into each other:

link <- family(model.pois)$linkfun # link function: log for Poisson
ilink <- family(model.pois)$linkinv # inverse link function: exp for Poisson
all.equal(ilink(pred.l), pred.r)
## [1] TRUE
all.equal(pred.l, link(pred.r))
## [1] TRUE

There is also the type = "terms" setting but this one is rarely used an also available in predict.lm.

Obtaining confidence intervals

The predict function of GLMs does not support the output of confidence intervals via interval = "confidence" as for predict.lm. We can still obtain confidence intervals for predictions by accessing the standard errors of the fit by predicting with se.fit = TRUE:

predict.confidence <- function(object, newdata, level = 0.95, ...) {
    if (!is(object, "glm")) {
        stop("Model should be a glm")
    }
    if (!is(newdata, "data.frame")) {
        stop("Plase input a data frame for newdata")
    }
    if (!is.numeric(level) | level < 0 | level > 1) {
        stop("level should be numeric and between 0 and 1")
    }
    ilink <- family(object)$linkinv
    ci.factor <- qnorm(1 - (1 - level)/2)
    # calculate CIs:
    fit <- predict(object, newdata = newdata, level = level, 
                    type = "link", se.fit = TRUE, ...)
    lwr <- ilink(fit$fit - ci.factor * fit$se.fit)
    upr <- ilink(fit$fit + ci.factor * fit$se.fit)
    df <- data.frame("fit" = ilink(fit$fit), "lwr" = lwr, "upr" = upr)
    return(df)
}

Using this function, we get the following confidence intervals for the Poisson model:

conf.df <- predict.confidence(model.pois, ozone[testset,])
head(conf.df, 2)
##        fit      lwr       upr
## 8 9.073256 8.141969 10.111065
## 9 5.409964 4.619623  6.335518

Using the confidence data, we can create a function for plotting the confidence of the estimates in relation to individual features:

plot.confidence <- function(df, feature) {
    library(ggplot2)
    p <- ggplot(df, aes_string(x = feature, 
                       y = "fit")) + 
      geom_line(colour = "blue") + 
      geom_point() + 
      geom_ribbon(aes(ymin = lwr, ymax = upr), 
                    alpha = 0.5) 
      return(p)
}
plot.confidence.features <- function(data, features) {
    plots <- list()
    for (feature in features) {
        p <- plot.confidence(data, feature)
        plots[[feature]] <- p
    }
    library(gridExtra)
    #grid.arrange(plots[[1]], plots[[2]], plots[[3]])
    do.call(grid.arrange, plots)
}

Using these functions, we can generate the following plot:

data <- cbind(ozone[testset,], conf.df)
plot.confidence.features(data, colnames(ozone))