5 Lab 9 — December 1
5.2 Logistic regression with binomial (\(n > 1\)) data
Note: For your assignment, you should use the code provided in the logitfit1 file on
Brightspace. However, the code and output shown below are equivalent and also fitted via
iteratively reweighted least squares (IRWLS).
5.2.1 Shocking some cows
Consider the data below, that is similar to Table 12.3 in the handout, but with the number of
responses in the first experiment changed from 0 to 2. Note that the code in the logitfit1 file
uses a value of 0.0001 rather than 0 to prevent its resulting model deviance from being NaN.
| Current (mA) | Trials | Responses | Proportion of Responses |
|---|---|---|---|
| 0 | 70 | 2 | 0.029 |
| 1 | 70 | 9 | 0.129 |
| 2 | 70 | 21 | 0.300 |
| 3 | 70 | 47 | 0.671 |
| 4 | 70 | 60 | 0.857 |
| 5 | 70 | 63 | 0.900 |
For the \(i^{\text{th}}\) experiment:
- \(x_{i}\) is the amount of current applied
- \(n_{i}\) is the number of trials
- \(y_{i}\) is the number of responses
- \(n_{i} - y_{i}\) is the number of no-responses
- \(y_{i} / n_{i}\) is the proportion of responses
We wish to fit the model:
\[E\left(\frac{y_{i}}{n_{i}}\right) \,=\, \theta_{i} \,=\, \frac{\exp{(\beta_{0} + \beta_{1}x)}} {1 + \exp{(\beta_{0} + \beta_{1}x)}}\]
Inputting the data into R:
cows <- tibble(
current = 0:5,
trials = rep(70, 6),
responses = c(2, 9, 21, 47, 60, 63),
proportion = responses / trials
)
cows## # A tibble: 6 x 4
## current trials responses proportion
## <int> <dbl> <dbl> <dbl>
## 1 0 70 2 0.0286
## 2 1 70 9 0.129
## 3 2 70 21 0.3
## 4 3 70 47 0.671
## 5 4 70 60 0.857
## 6 5 70 63 0.9
5.2.2 Fitting the model
When fitting the model, we represent the dependent variable as a matrix of successes and failures. As such the representation of our dependent variable will look like:
## responses no_responses
## [1,] 2 68
## [2,] 9 61
## [3,] 21 49
## [4,] 47 23
## [5,] 60 10
## [6,] 63 7Fitting the model now:
Note that in the family argument, we specify binomial but we do not need to wrap it in quotes
as it is a special object.
summary(full_model)##
## Call:
## glm(formula = cbind(responses, trials - responses) ~ current,
## family = binomial, data = cows)
##
## Deviance Residuals:
## 1 2 3 4 5 6
## -0.6334 0.0131 -0.4350 1.0587 0.4748 -1.4322
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.1020 0.3075 -10.09 <2e-16 ***
## current 1.1837 0.1067 11.10 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 232.6660 on 5 degrees of freedom
## Residual deviance: 3.9881 on 4 degrees of freedom
## AIC: 31.313
##
## Number of Fisher Scoring iterations: 4This gives us the fitted equation on the predictor scale (log-odds):
\[\widehat{\eta}_{i} \,=\, \widehat{\beta}_{0} + \widehat{\beta}_{1}X_{i} \,=\, -3.1020 + 1.1837X_{i}\]
To transform it onto the response scale (probabilities), we use the relationship that:
\[\widehat{\eta}_{i} \,=\, g(\widehat{\theta}_{i}) \,=\, \log{\left(\frac{\widehat{\theta}_{i}}{1 - \widehat{\theta}_{i}}\right)}\]
thus,
\[\widehat{\theta}_{i} \,=\, g^{-1}(\widehat{\eta}_{i}) \,=\, \frac{\exp{(\widehat{\eta}_{i})}}{1 + \exp{(\widehat{\eta}_{i})}}\]
5.2.3 Obtaining the predicted probabilities
As with the OLS linear models, we can pass a GLM into the predict() function which will then
call predict.glm().
predict(full_model)## 1 2 3 4 5 6
## -3.1019825 -1.9183311 -0.7346797 0.4489717 1.6326231 2.8162745However, notice that that some of these values are outside of \([0,1]\)! This is because these are
the predicted values on the predictor scale (log-odds). If we want the predicted values on the
response scale (probabilities), from the documentation (?predict.glm), we need to specify that
type="response".
predict(full_model, type="response")## 1 2 3 4 5 6
## 0.04302555 0.12804778 0.32416863 0.61039471 0.83652865 0.943548965.2.4 Visualising the fit
Since our data is quite small, unlike with straight line fits where we can just proceed to connect the dots and still obtain a straight line, here, we need to supplement our data with a sufficient number of grid points to get as close as possible to the true shape of the predicted curve.
fitted_vals <- full_model %>%
predict(newdata = tibble(current=seq(0, 5, 0.01)), type="response") %>%
as_tibble_col(column_name=".fitted") %>%
add_column(current=seq(0, 5, 0.01), .before=".fitted")
fitted_vals## # A tibble: 501 x 2
## current .fitted
## <dbl> <dbl>
## 1 0 0.0430
## 2 0.01 0.0435
## 3 0.02 0.0440
## 4 0.03 0.0445
## 5 0.04 0.0450
## 6 0.05 0.0455
## 7 0.06 0.0460
## 8 0.07 0.0466
## 9 0.08 0.0471
## 10 0.09 0.0476
## # ... with 491 more rows
ggplot(data=NULL, aes(x=current))+
geom_point(data=cows, aes(y=proportion))+
geom_line(data=fitted_vals, aes(y=.fitted), colour="#3366FF", lwd=1.5, alpha=0.6)+
labs(x="Current (mA)", y="Proportion of responses to electric shock")
5.2.5 Analysis of deviance
The analogue of the residual sum of squares for GLMs is the deviance. We can compute the deviance of a GLM by passing our GLM into the deviance function.
deviance(full_model)## [1] 3.988104Similar to the residual sum of squares in OLS, the deviance has degrees of freedom \(n - p - 1\), where \(n\) is the number of observations, and \(p\) is the number of non-intercept parameters in the model.
We can perform a test for model usefulness, i.e.
\[H_{0}: \beta_{1} = 0, \quad H_{A}: \beta_{1} \neq 0\]
by computing the difference in the deviance of the reduced model (with degrees of freedom \(n - q - 1\)) with the deviance of the full model (with degrees of freedom \(n - p - 1\)). Then the test statistic is
\[G^{2} \,=\, \mathcal{D}_{\text{reduced}} \,-\, \mathcal{D}_{\text{full}}\]
which has a chi-square distribution on \(p - q\) degrees of freedom.
Similar to OLS, we can pass our reduced and full GLMs into the anova() function, with the
additional specification of test="Chisq".
## Analysis of Deviance Table
##
## Model 1: cbind(responses, trials - responses) ~ 1
## Model 2: cbind(responses, trials - responses) ~ current
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 5 232.666
## 2 4 3.988 1 228.68 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The test statistic has a value of 228.68, and follows a chi-square distribution on 1 degree of freedom. As the \(p\)-value is less than 0.05, we can reject the null hypothesis and conclude that the full model is an improvement over the reduced model.