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

¡@

    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