Prediction and Effective Doses (Reading: Faraway (2006, 1st ed.), section 2.10)
¡@
For the Bliss insect data studied in the previous lab, we show below how to predict the response, i.e., the probability of "sucess", using binomial GLM:
> bliss <- read.table("bliss.txt")
> modl <- glm(cbind(dead,alive) ~ conc, family=binomial, data=bliss)
> lmodsum <- summary(modl)
> lmodsum
Call:
glm(formula = cbind(dead, alive) ~ conc, family = binomial, data = bliss)
Deviance Residuals:
1 2 3 4 5
-0.4510 0.3597 0.0000 0.0643 -0.2045
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.3238 0.4179 -5.561 2.69e-08 ***
conc 1.1619 0.1814 6.405 1.51e-10 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 64.76327 on 4 degrees of freedom
Residual deviance: 0.37875 on 3 degrees of freedom
AIC: 20.854
Number of Fisher Scoring iterations: 4
¡@
Let us predict the response at dose of 2.5:
>
x0 <- c(1,2.5) # remember to include the intercept 1
We need to first predict
hx:
>
eta0 <- sum(x0*coef(modl))
> eta0
[1] 0.5809475
So, the predicted value of px is:
> ilogit(eta0)
[1] 0.6412854
A 64.13% predicted chance of death at this dose.
¡@
Now, let us compute a 95% confidence interval (CI) for this predicted probability. We need to first extract the estimated covariance matrix of the coefficients b0 and b1:
> (cm <- lmodsum$cov.unscaled)
(Intercept) conc
(Intercept) 0.17463024 -0.06582336
conc -0.06582336 0.03291168
The standard error of the predicted value of hx is:
>
se <- sqrt( t(x0) %*% cm %*% x0)
> se
[,1]
[1,] 0.2262995
A 95% CI for hx is then:
> c(eta0-1.96*se,eta0+1.96*se)
[1] 0.1374005 1.0244944
So, the 95% CI for the predicted value of px is:
> ilogit(c(0.1374005,1.0244944))
[1] 0.5342962 0.7358471
A more direct way of obtaining the same result is:
> predict(modl,newdata=data.frame(conc=2.5),se=T)
$fit
[1] 0.5809475
$se.fit
[1] 0.2262995
$residual.scale
[1] 1
> ilogit(c(0.5809475-1.96*0.2262995,0.5809475+1.96*0.2262995))
[1] 0.5342962 0.7358472
¡@
In the case of Normal linear model, the CI of an extrapolation is usually much wider than that of an interpolation. Does it still hold in binomial GLM? Let us try predicting px at a very low dose of -5 (note that the range of dose in the dataset is [0, 4]):
>
x0 <- c(1,-5)
>
se <- sqrt( t(x0) %*% cm %*% x0)
>
eta0 <- sum(x0*modl$coef)
A 95% CI for hx is:
> c(eta0-1.96*se,eta0+1.96*se)
[1] -10.655241 -5.611288
Notice that this CI for is much wider in length than that at dose=2.5. The 95% CI for px is:
> ilogit(c(eta0-1.96*se,eta0+1.96*se))
[1] 2.357639e-05 3.643038e-03
Notice that
this is not a wide interval in absolute terms (i.e., the difference of lower and upper limits is small, even smaller than the CI at dose=2.5).
But, in relative terms (i.e., the ratio of lower and upper limits), it certainly is.
The upper limit is about 100 times larger than the lower limit (c.f., for the case of dose=2.5, the ratio is 0.7358472/0.5342962=1.377227).
¡@
We now use the Bliss insect data to predict the lethal dose corresponding to p=1/2, denoted by LD50. Under the logit link, LD50 can be estimated by:
> (ld50 <- -modl$coef[1]/modl$coef[2])
(Intercept)
2
The standard error of LD50 can be approximated by
d-method which, in this example,
works out as:
> dr <- c(-1/modl$coef[2],modl$coef[1]/modl$coef[2]^2)
> dr
conc (Intercept)
-0.860663 -1.721326
> sqrt(dr %*% lmodsum$cov.un %*% dr)[,]
[1] 0.1784367
So the 95% CI of LD50 is given by:
> c(2-1.96*0.1784367,2+1.96*0.1784367)
[1] 1.650264 2.349736
Other levels may be consider --- below is an example for p=0.9:
> ed90 <- (log(0.9/(1-0.9))-modl$coef[1])/modl$coef[2]
> ed90
(Intercept)
3.89107
(exercise: please try calculating its standard error by yourself and check with the correct answer given below.)
More conveniently, we may use the dose.p function in the MASS package:
> library(MASS)
> dose.p(modl,p=c(0.5,0.9))
Dose SE
p = 0.5: 2.00000 0.1784367
p = 0.9: 3.89107 0.3449965