Description
1. The table below reports results from a toxicological experiment, including the number of beetles killed (yi) after 5 hours exposure to gaseous carbon disulphide at various concentrations.
Concentration (log dose, xi) is given on the log10 scale.
Log Dose, xi Number of beetles, mi Number killed, yi
1.6907 59 6
1.7242 60 13
1.7552 62 18
1.7842 56 28
1.8113 63 52
1.8369 59 53
1.8610 62 61
1.8839 60 60
Consider a binomial response distribution, and assume that the yi are independent realizations
from Bin(mi
, πi), i = 1,…,n. The objective is to study the effect of the choice of link function
g(·), where πi = g
−1
(ηi) = g
−1
(β1 + β2xi).
(a) Using R, fit 3 binomial GLMs for these data corresponding to 3 link functions, logit, probit
and complementary log-log. Perform residual analysis for each model, using the deviance residuals. Obtain fitted values, ˆπi
, under each model and compare with observed proportions, yi/mi
.
Obtain the estimated dose-response curve under each model by evaluating ˆπ(x) = g
−1
(βˆ
1 +βˆ
2x)
over a grid of values x for log dose in the interval (1.65, 1.9). Plot these curves and compare
with the scatter plot of the observed xi plotted against the observed proportions. Based on all
the results above, discuss the fit of the different models.
(b) One of the more general (parametric) link functions for binomial GLMs that has been
suggested in the literature is defined through
g
−1
α (ηi) = exp(αηi)
{1 + exp(ηi)}
α
for α > 0. (1.1)
Note that the logit link arises as a special case of (1.1), when α = 1. Discuss the effect of the
additional model parameter α, in particular, for values 0 < α < 1 and α > 1. Provide the
expression for the log-likelihood for β1, β2 and α under the link in (1.1), and discuss the complications that arise for maximum likelihood estimation under this more general model compared
with the logit GLM. (You do not need to fit the model, estimates are given below.)
(c) The MLEs under the model with link given in (1.1) are βˆ
1 = −113.625, βˆ
2 = 62.5 and
αˆ = 0.279. (The MLEs can be obtained using the Newton-Raphson method.) Using these estimates, obtain the fitted values ˆπi and the estimated dose-response curve under the link (1.1).
Compare with the corresponding results under the 3 models in (a). Obtain the deviance residuals from the model with link (1.1) and analyze them graphically.
(d) Compute the AIC and BIC for the 4 models considered above to compare them.
2. This problem involves Bayesian analysis of the beetle mortality data from the previous problem.
(a) Consider a Bayesian binomial GLM with a complementary log-log link, i.e., assume that,
given β1 and β2, the yi are independent from Bin(mi
, π(xi)), i = 1,…,8, where
π(x) ≡ π(x; β1, β2) = 1 − exp{− exp(β1 + β2x)}.
Design and implement an MCMC method to sample from the posterior distribution of (β1, β2).
Study the effect of the prior for (β1, β2), for example, consider a flat prior as well as (independent) normal priors. Under the flat prior, obtain the posterior distribution for the median
lethal dose, LD50, that is, the dose level at which the probability of response is 0.5. Finally, plot
point and interval estimates for the dose-response curve π(x) (over a grid of values x for log dose).
(b) Next, consider a binomial GLM with a logit link, i.e., now the yi are assumed independent, given β1 and β2, from Bin(mi
, π(xi)), i = 1,…,8, where
π(x) ≡ π(x; β1, β2) = exp(β1 + β2x)/{1 + exp(β1 + β2x)}.
Working with a flat prior for (β1, β2), obtain MCMC samples from the posterior distributions
for β1, β2, and for LD50, along with point and interval estimates for the dose-response curve π(x).
(c) As a third model, consider the binomial GLM with the parametric link given in (1.1).
Develop an MCMC method to sample from the posterior distribution of (β1, β2, α), and obtain
the posterior distribution for LD50, and point and interval estimates for π(x).
(d) Use the results from parts (a), (b) and (c) for an empirical comparison of the three Bayesian
binomial GLMs for the beetle mortality data. Moreover, perform a residual analysis for each
model using the Bayesian residuals: (yi/mi)−π(xi
; β1, β2) for the first two models, and (yi/mi)−
π(xi
; β1, β2, α) for the third. Finally, use the quadratic loss L measure for formal comparison of
the three models.