Practical Exercise: The second Beetles exercise - Solutions

This exercise is useful for the following courses:

Aims of this practical exercise

In this exercise you will:

  • fit a variety of models to the Beetles data involving a variety of link functions.

Get the data

data(Beetles, package="ELMER")
Beetles = Beetles |> mutate(NoLiving = NoBeetles - NoKilled) |> glimpse()
Rows: 8
Columns: 6
$ Dose      <dbl> 1.6907, 1.7242, 1.7552, 1.7842, 1.8113, 1.8369, 1.8610, 1.88…
$ NoBeetles <int> 59, 60, 62, 56, 63, 59, 62, 60
$ NoKilled  <int> 6, 13, 18, 28, 52, 53, 61, 60
$ Prop      <dbl> 0.1016949, 0.2166667, 0.2903226, 0.5000000, 0.8253968, 0.898…
$ TProp     <dbl> 0.3245649, 0.4841708, 0.5690309, 0.7853982, 1.1397128, 1.246…
$ NoLiving  <int> 53, 47, 44, 28, 11, 6, 1, 0

The Exercise

Revisit the Beetle mortality data. Investigate the predictions that arise through use of different link functions.

The solution

First obtain the models using different link functions, and then obtain and compare the predictions.

Beetles.logit = glm(cbind(NoKilled,NoLiving)~Dose,data=Beetles,family=binomial(link="logit"))
Beetles.logit |> predict(type = "response")
         1          2          3          4          5          6          7 
0.05860103 0.16402787 0.36211901 0.60531491 0.79517177 0.90323582 0.95519611 
         8 
0.97904934 
Beetles.probit = glm(cbind(NoKilled,NoLiving)~Dose,data=Beetles,family=binomial(link="probit"))
Beetles.probit |> predict(type = "response")
         1          2          3          4          5          6          7 
0.05691142 0.17869349 0.37874084 0.60384830 0.78754962 0.90370973 0.96233307 
         8 
0.98713279 
Beetles.cloglog = glm(cbind(NoKilled,NoLiving)~Dose,data=Beetles,family=binomial(link="cloglog"))
Beetles.cloglog |> predict(type = "response")
         1          2          3          4          5          6          7 
0.09473644 0.18801129 0.33797124 0.54231139 0.75835580 0.91767335 0.98569859 
         8 
0.99912042 
Beetles.cauchy = glm(cbind(NoKilled,NoLiving)~Dose,data=Beetles,family=binomial(link="cauchit"))
Beetles.cauchy |> predict(type = "response")
         1          2          3          4          5          6          7 
0.08336583 0.13195791 0.26271969 0.60407630 0.81464759 0.88446172 0.91558421 
         8 
0.93297414

We cannot easily obtain interval estimates. The intervals are symmetric on the linear scale, and asymmetric once back-transformed onto the original proportions scale.

You should see though that the pattern for the fitted values is dependent on the link used.