Confidence interval on the fit

We can estimate a confidence interval around our model fit using the se.fit = TRUE argument in the predict() function to obtain the standard errors on the model coefficients and then estimating the 100(1 - \(\alpha\))% CI. Note, however, that one must be careful when doing so, such that we obtain the errors on the predictions with respect to the linear predictor, and not the response.

Here is an example of how not to do it. By specifying type = "response" we are asking for the standard errors in the response space.

## get the SE of the response
se_wrong <- predict(fit_mod, newdata, type = "response", se.fit = TRUE)$se.fit
## wrong upper 95% CI
wrong_up <- p_hat + 1.96*se_wrong
## wrong lower 95% CI
wrong_lo <- p_hat - 1.96*se_wrong

## set plot area
par(mai = c(0.9, 0.9, 0.1, 0.1),
    omi = c(0, 0, 0, 0),
    cex.lab = 1.5)

## plot age vs length
plot(ll, yy, pch = 16, cex = 1.3, col = clr, ylim = c(-0.1, 1.1),
     yaxt = "n", ylab = "Smolt age", xlab = "Length (mm)")
axis(2, at = c(0,1), labels = c(1,2), las = 1)
## add model fit
lines(sl, p_hat, lwd = 2)
## add wrong CI's
lines(sl, wrong_up, lwd = 1, col = "gray")
lines(sl, wrong_lo, lwd = 1, col = "gray")

Notice that the CI’s overlap 0 and 1, which necessarily violates our constraint that probability must lie between 0 and 1.

Now here’s how to estimate the CI’s properly.

## get the SE of the response
se_right <- predict(fit_mod, newdata, type = "link", se.fit = TRUE)$se.fit
## right upper 95% CI
right_up <- 1 / (1 + exp(-(eta + 1.96*se_right)))
## right lower 95% CI
right_lo <- 1 / (1 + exp(-(eta - 1.96*se_right)))

## set plot area
par(mai = c(0.9, 0.9, 0.1, 0.1),
    omi = c(0, 0, 0, 0),
    cex.lab = 1.5)

## plot age vs length
plot(ll, yy, pch = 16, cex = 1.3, col = clr, ylim = c(-0.1, 1.1),
     yaxt = "n", ylab = "Smolt age", xlab = "Length (mm)")
axis(2, at = c(0,1), labels = c(1,2), las = 1)
## add model fit
lines(sl, p_hat, lwd = 2)
## add right CI's
lines(sl, right_up, lwd = 1, col = "gray")
lines(sl, right_lo, lwd = 1, col = "gray")

Finding the probability of 0.5

Sometimes we would like to know at which level of the predictor variable do we expect an equal probability of obtaining a 0 or 1. In this case we can ask, “What is the length at which the probability of smolting at age-2 is 0.5?”

\[ p_i = \frac{1}{1 + \exp(\text{-} \mathbf{X}_i \boldsymbol{\beta})} \\ \Downarrow \\ 0.5 = \frac{1}{1 + \exp(14 - 0.17 L_{0.5})} \\ \Downarrow \\ L_{0.5} \approx 82 ~ \text{mm} \]

Computing the odds

We have seen that we can express the results from a logistic regression in terms of the odds. The odds \(o\) are an unbounded alternative to probability \(p\), such that if we represent the \(k\)-to-1 odds against something as \(1 / k\), then the following holds

\[ o = \frac{p}{1 - p} ~ \Rightarrow ~ p = \frac{o}{1 + o} \]

So, for example, if \(p\) = 0.8, then \(o = \frac{0.8}{1 - 0.8} = 4\). In a logistic regression model, the logit function gives us the log odds, such that

\[ \text{logit}(p) = \mathbf{X} \boldsymbol{\beta} \\ \Downarrow \\ \log \left( \frac{p}{1 - p} \right) = \mathbf{X} \boldsymbol{\beta} \\ \Downarrow \\ \log ( \text{odds} ) = \mathbf{X} \boldsymbol{\beta} \\ \Downarrow \\ \text{odds} = \exp (\mathbf{X} \boldsymbol{\beta}) \]

For our smolt age model, the log(odds) are

\[ \log \left( \frac{p}{1 - p} \right) = \text{-}14 + 0.17 L \\ \Downarrow \\ \log ( \text{odds} ) = \text{-}14 + 0.17 L \]

which mean that a 1 mm increase in length \(L\) increases the log-odds by 0.17. By the same reasoning, the odds themselves as given by

\[ \log \left( \frac{p}{1 - p} \right) = \text{-}14 + 0.17 L \\ \Downarrow \\ \log ( \text{odds} ) = \text{-}14 + 0.17 L \\ \Downarrow \\ \text{odds} = \exp (\text{-}14 + 0.17 L) \]

such that a 1 mm increase in length \(L\) increases the odds by exp(0.17) \(\approx\) 1.19 = 19%.

Significance test for \(\beta_i\)

An alternative to the likelihood ratio test is a \(z\) test, where

\[ z = \frac{\hat{\beta_i}}{\text{SE}(\hat{\beta_i})} \sim z_{\alpha / 2} \]

The summary() output from a fitted GLM includes this the results of this test.

## summary table
faraway::sumary(fit_mod)

##             Estimate Std. Error z value  Pr(>|z|)
## (Intercept) 13.55885    3.18639  4.2552 2.088e-05
## length      -0.16735    0.03870 -4.3243 1.531e-05
## 
## n = 80 p = 2
## Deviance = 42.55364 Null Deviance = 110.85354 (Difference = 68.29991)

As you can see, the \(p\)-value for the \(z\) test is also very small.

Confidence interval for \(\beta_i\)

We can also compute a 100(1 - \(\alpha\))% confidence interval on the regression coefficient as

\[ \hat{\beta_i} \pm z_{\alpha / 2} \text{SE}(\hat{\beta_i}) \]

which we do in R via

## beta
beta_1 <- coef(fit_mod)[2]
## SE of beta
se_beta_1 <- sqrt(diag(vcov(fit_mod)))[2]
## 95% CI
beta_1 + c(-1,1) * 1.96 * se_beta_1

## [1] -0.2432010 -0.0914967

We can also compute CI’s based on the profile likelihood with confint()

## 95% CI via profile likelihood
confint(fit_mod)

##                  2.5 %     97.5 %
## (Intercept)  8.3933573 21.1334846
## length      -0.2593926 -0.1045879

Due to possible bias in \(\text{SE}(\beta)\), however, we can also compute CI’s based on the profile likelihood, which we obtain by evaluating the likelihood along a sequence of possible parameter values.

## number of points to profile
nb <- 200
## possible beta's
beta_hat <- seq(-0.4, 0, length = nb)
## calculate neg-LL of possible beta's
pl <- rep(NA, nb)
for(i in 1:nb) {
  mm <- glm(age ~ 1 + offset(beta_hat[i] * length), data = df,
            family = binomial(link = "logit"))
  pl[i] <- -logLik(mm)
}

## set plot area
par(mai = c(0.9, 0.9, 0.1, 0.1),
    omi = c(0, 0, 0, 0),
    cex.lab = 1.5)

## plot likelihood profile
plot(beta_hat, pl, type = "l", las = 1,
     ylab = "Negative log-likelihood", xlab = expression(beta))
crit <- -(logLik(fit_mod) - qchisq(0.95, 1) / 2)
abline(h = crit, lty = "dashed")
points(confint(fit_mod)[2,], c(crit, crit), pch = 16, col = "blue")

Model selection via AIC

As with other linear models, we can use AIC to help select among possible models.

\[ \begin{align} AIC &= 2 k - 2 \log \mathcal{L} \\ &= 2 k + D \end{align} \]

Here are three different ways to compute the AIC for our smolt age model.

## AIC
AIC(fit_mod) 
## AIC via likelihood
(2 * 2) - 2 * logLik(fit_mod)
## AIC via deviance
(2 * 2) + summary(fit_mod)$deviance

## [1] 46.55364
## 'log Lik.' 46.55364 (df=2)
## [1] 46.55364

We can compare our model with length as a predictor to a null model with no predictors via the following code.

## fit null model
fit_null <- glm(age ~ 1, data = df,
                family = binomial(link = "logit"))
faraway::sumary(fit_null)

##             Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.05001    0.22368 -0.2236   0.8231
## 
## n = 80 p = 1
## Deviance = 110.85354 Null Deviance = 110.85354 (Difference = 0.00000)

## difference in AIC
AIC(fit_null) - AIC(fit_mod)

## [1] 66.29991

Residuals

As with other models, it’s important to examine diagnostic checks for our fitted models. One of the first things we can do is inspect the model residuals. We usually think about residuals \(e\) as

\[ e = y - \hat{y} \]

but with logistic regression, the response can take 1 of 2 possible values. Let’s plot them to see how they look.

## raw residuals
res_raw <- df$age - predict(fit_mod, type = "response")
## equivalently
res_raw <- residuals(fit_mod, type = "response")

## set plot area
par(mfrow = c(1, 2),
    mai = c(0.9, 0.9, 0.1, 0.1),
    omi = c(0, 0, 0, 0),
    cex.lab = 1.5)

## plot resids vs eta
plot(predict(fit_mod), res_raw, las = 1, pch = 16,
     ylab = "Residuals", xlab = "Linear predictor",
     main = "")
## plot resids vs response
plot(predict(fit_mod, type = "response"), res_raw, las = 1, pch = 16,
     ylab = "Residuals", xlab = "Response",
     main = "")

Deviance residuals

We can instead use the deviance residuals, which we define as

\[ e_i = (2 y_i - 1) D_i \]

where \(2 y - 1 = 1\) if y is 1, and \(2 y - 1 = -1\) if y is 0. These are the default residuals returned by residuals(). We then place the deviance residuals into bins for easier inspection, noting that

• it’s sensitive to the number of bins (~30/bin is good)

• the mean of \(e\) is not constrained to be 0

• ~95% of points should fall within the CI

We can use binnedplot() from the arm package to do this automatically for us.

library(arm)

## set plot area
par(mai = c(0.9, 0.9, 0.1, 0.1),
    omi = c(0, 0, 0, 0),
    cex.lab = 1.5)

## plot resids vs eta
binnedplot(fitted(fit_mod), residuals(fit_mod), las = 1, pch = 16,
     ylab = "Residuals", xlab = "Fitted values",
     main = "")

Leverage

We can also calculate the leverages \(h\) to look for unusual observation in predictor space. Recall we are potentially concerned about any

\[ h_i > 2 \frac{k}{n} \]

We can use hatvalues() to calculate the leverages and faraway::halfnorm() to plot them.

## set plot area
par(mai = c(0.9, 0.9, 0.1, 0.1),
    omi = c(0, 0, 0, 0),
    cex.lab = 1.5)

## leverages
levs <- hatvalues(fit_mod)

## threshhold value
h_crit <- 2 * length(coef(fit_mod)) / nn

## halfnormal plot
faraway::halfnorm(levs, las = 1)
text(0, 0.92*par("usr")[4], substitute(italic(h[crit]) == h_crit, list(h_crit = h_crit)), pos = 4)

Cook’s Distance

Recall that we can also use Cook’s \(D\) to identify potentially influential points. Cook’s \(D\) is defined as

\[ D_{i}=e_{i}^{2} \frac{1}{k}\left(\frac{h_{i}}{1-h_{i}}\right) \]

In general, we should be concerned about an \(D_i > F^{(0.5)}_{n, n - k} \approx 1\). We can compute the \(D_i\) with cooks.distance() and again use faraway::halfnorm() to plot them.

## set plot area
par(mfrow = c(1, 2),
    mai = c(0.9, 0.9, 0.1, 0.1),
    omi = c(0, 0, 0, 0),
    cex.lab = 1.5)

## Cook's D
CD <- cooks.distance(fit_mod)

## halfnormal plot
faraway::halfnorm(CD, las = 1)

## plot age v
plot(ll, yy + 1, las = 1, pch = 16, cex = 1.3, col = clr,
     yaxt = "n", ylab = "", xlab = "Length (mm)")
lines(sl, p_hat + 1, lwd = 2)
points(df[43,"length"], df[43, "age"] + 1, pch = 1, cex = 2)
text(df[43,"length"], df[43, "age"] + 1, "43", pos = 3, offset = 0.8)
mtext(side = 2, "Smolt age", line = 1, cex = 1.5)
axis(2, at = c(1,2), las = 1)

The plot on the right shows the data, model fit, and the point with the largest \(D\).

Binned predictions

It’s hard to compare our predictions on the real interval [0,1] to discrete binary outcomes that take a value of 0 or 1. To help, we can compute \(\hat{y}\) for bins of data. Here are the smolt age data split into 8 equal sized bins along the length predictor.

## cut data
l_cut <- cut(df$length, seq(40, 120, 10))
y_bin <- by(df$age, l_cut, mean)
y_cut <- cut(fitted(fit_mod), quantile(fitted(fit_mod), probs = seq(0, 1, 1/8)))
p_bin <- by(fitted(fit_mod), y_cut, mean)
  
## set plot area
par(mai = c(0.9, 0.9, 0.1, 0.1),
    omi = c(0, 0, 0, 0),
    cex.lab = 1.5)
## plot age v
plot(ll, yy + 1, las = 1, pch = 16, cex = 1.3, col = clr,
     yaxt = "n", ylab = "Smolt age", xlab = "Length (mm)")
lines(seq(40, 120), p_hat + 1, lwd = 2, col = "gray")
points(seq(45, 115, 10), y_bin + 1, pch = 16, cex = 1.2)
abline(v = c(seq(40, 120, 10)), lty = "dashed")
axis(2, at = c(1,2), las = 1)

The black points are the estimated proportion of smolts, which seem to agree rather well with the model fit (gray line). We can now compare them to one another.

## set plot area
par(mai = c(0.9, 0.9, 0.1, 0.1),
    omi = c(0, 0, 0, 0),
    cex.lab = 1.5)
## plot age v
plot(y_bin, rev(p_bin), las = 1, pch = 16, cex = 1.3,
     ylab = "Observed", xlab = "Predicted")
abline(a = 0, b = 1, col = "gray")

Hosmer-Lemeshow test

We can formalize this binned comparison with the Hosmer-Lemeshow test, which is given by

\[ HL = \sum_{j = 1}^J \frac{(y_j - m_j \hat{p}_J)^2}{m_j \hat{p}_J(1 - \hat{p}_J)} \sim \chi^2_{(J - 1)} \]

where \(J\) is the number of groups and \(y_j = \sum y_{i = j}\). We can perform the Hosmer-Lemeshow test with generalhoslem::logitgof(). The null hypothesis is the model fits the data well.

## H-L test with 8 groups
generalhoslem::logitgof(obs = df$age, exp = fitted(fit_mod), g = 8)

## 
##  Hosmer and Lemeshow test (binary model)
## 
## data:  df$age, fitted(fit_mod)
## X-squared = 1.7874, df = 6, p-value = 0.9382

The \(p\)-value is large so we conclude an adequate fit.

Proportion of variance explained

Calculating \(R^2\) for logistic models is not the same as linear models. Given the deviance \(D_M\) for our model and a null model \(D_0\), we define

\[ R^2 = \frac{1 - \exp \left( [D_M - D_0]/n \right)}{1 - \exp(\text{-}D_0 / n)} \]

Here is the \(R^2\) for our smolt-at-age model.

## deviances
DM <- fit_mod$deviance
D0 <- fit_mod$null.deviance
# R^2
R2 <- (1 - exp((DM - D0) / nn)) / (1 - exp(-D0 / nn))
round(R2, 2)

## [1] 0.77