Case Study

Stat 203 Lecture 28

Dr. Janssen

Setup

Description

An experiment exposed batches of insects to various deposits (in mg) of insecticides. The proportion of insects \(y\) killed after six days of exposure in each batch of size \(m\) is potentially a function of the dose and type of insecticide.

library(GLMsData); data(deposit);
str(deposit)
'data.frame':   18 obs. of  4 variables:
 $ Killed     : int  3 5 19 19 24 35 2 14 20 27 ...
 $ Number     : int  50 49 47 38 29 50 50 49 50 50 ...
 $ Insecticide: Factor w/ 3 levels "A","B","C": 1 1 1 1 1 1 2 2 2 2 ...
 $ Deposit    : num  2 2.64 3.48 4.59 6.06 8 2 2.64 3.48 4.59 ...

Exploring the data

Code 
deposit$Prop <- deposit$Killed / deposit$Number
plot( Prop ~ Deposit, type="n", las=1, ylim=c(0, 1),
     data=deposit, main="Proportion of\ninsects killed",
     xlab="Deposit (in mg)", ylab="Proportion killed")
points( Prop ~ Deposit, pch="A", subset=(Insecticide=="A"), data=deposit)
points( Prop ~ Deposit, pch="B", subset=(Insecticide=="B"), data=deposit)
points( Prop ~ Deposit, pch="C", subset=(Insecticide=="C"), data=deposit)
ins.m1 <- glm(Killed/Number ~ Deposit + Insecticide,
       family = binomial, weights = Number, data = deposit)
newD <- seq( min(deposit$Deposit), max(deposit$Deposit), length=100)
newProp.logA <- predict(ins.m1, type="response",
     newdata=data.frame(Deposit=newD, Insecticide="A") )
newProp.logB <- predict(ins.m1, type="response",
     newdata=data.frame(Deposit=newD, Insecticide="B") )
newProp.logC <- predict(ins.m1, type="response",
     newdata=data.frame(Deposit=newD, Insecticide="C") )
lines( newProp.logA ~ newD, lty=1); lines( newProp.logB ~ newD, lty=2)
lines( newProp.logC ~ newD, lty=3)

Interlude: ED50

Median Effective Dose

Common use of binomial GLMs: examine relationship between does \(d\) of a drug or poison and the proportion \(y\) of insects (or plants or animals) that survive.

The median effective dose is the dose affecting 50% of the individuals; shorthand: “ED50”.

ED50 in R

library(MASS)
cool_glm <- glm(response ~ exp, family="binomial")
dose.p(cool_glm)

Back to our insects

library(MASS)
dose.p(ins.m1, c(1, 2))
             Dose        SE
p = 0.5: 5.099708 0.2468085

Fitting without an intercept

ins.m1A <- update( ins.m1, .~. - 1)   # Do not fit a constant term
coef( ins.m1A )
     Deposit InsecticideA InsecticideB InsecticideC 
   0.6316762   -3.2213638   -2.8518371   -0.5333477

Fitting without an intercept

ins.m1A <- update( ins.m1, .~. - 1)   # Do not fit a constant term
coef( ins.m1A )
     Deposit InsecticideA InsecticideB InsecticideC 
   0.6316762   -3.2213638   -2.8518371   -0.5333477 
ED50s <- cbind( dose.p(ins.m1A, c(2, 1)),    dose.p(ins.m1A, c(3, 1)),
                   dose.p(ins.m1A, c(4, 1)) )
colnames(ED50s) <- c("Insect. A","Insect. B", "Insect. C"); ED50s
         Insect. A Insect. B Insect. C
p = 0.5:  5.099708  4.514714 0.8443372

Exploring Residuals

Recall

Code 
deposit$Prop <- deposit$Killed / deposit$Number
plot( Prop ~ Deposit, type="n", las=1, ylim=c(0, 1),
     data=deposit, main="Proportion of\ninsects killed",
     xlab="Deposit (in mg)", ylab="Proportion killed")
points( Prop ~ Deposit, pch="A", subset=(Insecticide=="A"), data=deposit)
points( Prop ~ Deposit, pch="B", subset=(Insecticide=="B"), data=deposit)
points( Prop ~ Deposit, pch="C", subset=(Insecticide=="C"), data=deposit)
ins.m1 <- glm(Killed/Number ~ Deposit + Insecticide,
       family = binomial, weights = Number, data = deposit)
newD <- seq( min(deposit$Deposit), max(deposit$Deposit), length=100)
newProp.logA <- predict(ins.m1, type="response",
     newdata=data.frame(Deposit=newD, Insecticide="A") )
newProp.logB <- predict(ins.m1, type="response",
     newdata=data.frame(Deposit=newD, Insecticide="B") )
newProp.logC <- predict(ins.m1, type="response",
     newdata=data.frame(Deposit=newD, Insecticide="C") )
lines( newProp.logA ~ newD, lty=1); lines( newProp.logB ~ newD, lty=2)
lines( newProp.logC ~ newD, lty=3)

Residuals

Code 
library(statmod)     # For  qresid()
plot( qresid(ins.m1) ~ fitted(ins.m1), type="n", las=1, ylim=c(-3, 3),
     main="Quantile resids. plotted\nagainst fitted values",
     xlab="Fitted values", ylab="Residuals")
abline(h = 0, col="grey")
points( qresid(ins.m1) ~ fitted(ins.m1), pch="A", type="b", lty=1,
           subset=(deposit$Insecticide=="A") )
points( qresid(ins.m1) ~ fitted(ins.m1), pch="B", type="b", lty=2,
           subset=(deposit$Insecticide=="B") )
points( qresid(ins.m1) ~ fitted(ins.m1), pch="C", type="b", lty=3,
           subset=(deposit$Insecticide=="C"))

Plotting log-odds

Code 
LogOdds <- with(deposit, log(Prop/(1-Prop)) )
plot( LogOdds ~ Deposit, type="n", xlab="Deposit", data=deposit,
     main="Logits plotted\nagainst Deposit", las=1)
points( LogOdds ~ Deposit, pch="A", type="b", lty=1,
           data=deposit, subset=(Insecticide=="A") )
points( LogOdds ~ Deposit, pch="B", type="b", lty=2,
           data=deposit, subset=(Insecticide=="B") )
points( LogOdds ~ Deposit, pch="C", type="b", lty=3,
           data=deposit, subset=(Insecticide=="C") )

What to do?

Often in dose-response models, the logarithm of the dose is used.

logDep

deposit$logDep <- log( deposit$Deposit )
ins.m2 <- glm(Killed/Number ~ logDep + Insecticide - 1,
       family = binomial, weights = Number, data = deposit)

logDep

library(MASS)
deposit$logDep <- log( deposit$Deposit )
ins.m2 <- glm(Killed/Number ~ logDep + Insecticide - 1,
       family = binomial, weights = Number, data = deposit)
ED50s <- cbind( dose.p(ins.m2, c(2, 1)),    dose.p(ins.m2, c(3, 1)),
                   dose.p(ins.m2, c(4, 1)) )
colnames(ED50s) <- c("Insect. A","Insect. B", "Insect. C"); exp(ED50s)
         Insect. A Insect. B Insect. C
p = 0.5:  4.688232  4.154625  1.753202

Checking ins.m2

Residuals

Code 
library(statmod)     # For  qresid()
plot( qresid(ins.m2) ~ fitted(ins.m2), type="n", las=1, ylim=c(-3, 3),
     main="Quantile resids. plotted\nagainst fitted values",
     xlab="Fitted values", ylab="Residuals")
abline(h = 0, col="grey")
points( qresid(ins.m2) ~ fitted(ins.m2), pch="A", type="b", lty=1,
           subset=(deposit$Insecticide=="A") )
points( qresid(ins.m2) ~ fitted(ins.m2), pch="B", type="b", lty=2,
           subset=(deposit$Insecticide=="B") )
points( qresid(ins.m2) ~ fitted(ins.m2), pch="C", type="b", lty=3,
           subset=(deposit$Insecticide=="C"))

Plotting log-odds against Deposit

Code 
LogOdds <- with(deposit, log(Prop/(1-Prop)) )
plot( LogOdds ~ log(Deposit), type="n", xlab="log(Deposit)", data=deposit,
     main="Logits plotted\nagainst log(Deposit)", las=1)
points( LogOdds ~ log(Deposit), pch="A", type="b", lty=1,
           data=deposit, subset=(Insecticide=="A") )
points( LogOdds ~ log(Deposit), pch="B", type="b", lty=2,
           data=deposit, subset=(Insecticide=="B") )
points( LogOdds ~ log(Deposit), pch="C", type="b", lty=3,
           data=deposit, subset=(Insecticide=="C") )

A third model

ins.m3 <- glm(Killed/Number ~ poly(logDep, 2) + Insecticide,
       family = binomial, weights = Number, data = deposit)
anova( ins.m2, ins.m3, test="Chisq")
Analysis of Deviance Table

Model 1: Killed/Number ~ logDep + Insecticide - 1
Model 2: Killed/Number ~ poly(logDep, 2) + Insecticide
  Resid. Df Resid. Dev Df Deviance Pr(>Chi)   
1        14     23.385                        
2        13     15.090  1   8.2949 0.003976 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Evaluating the model

Code 
newD <- seq( min(deposit$logDep), max(deposit$logDep), length=200)
plot( LogOdds ~ log(Deposit) , type="n", xlab="log(Deposit)^2", data=deposit,
  main="Logits plotted\nagainst log(Deposit)^2", las=1)
newProp4.logA <- predict(ins.m3, type="response",
     newdata=data.frame(logDep=newD, Insecticide="A") )
newProp4.logB <- predict(ins.m3, type="response",
     newdata=data.frame(logDep=newD, Insecticide="B") )
newProp4.logC <- predict(ins.m3, type="response",
     newdata=data.frame(logDep=newD, Insecticide="C") )
lines( newProp4.logA ~ newD, lty=1); lines( newProp4.logB ~ newD, lty=2)
lines( newProp4.logC ~ newD, lty=3)

Overdispersion?

No.

c( deviance( ins.m3 ), df.residual( ins.m3 ) )
[1] 15.09036 13.00000