Type II functional responses via GLM

Re-do Vonesh functional response analysis pp. 310-311. The only real disadvantage of inverse GLM (i.e., GLM with a link function $f(x)=1/x$ is a slightly less direct parameter interpretation: intercept = inverse attack rate $1/a$, slope = handling time $h$. (It's also harder to extend to a different model, e.g. a multi-predator case with multiplicative predator effects, or a Rogers random predator equation — see below.) The advantages are (1) more stable estimation; (2) don't need starting values (although sometimes you still do — inverse-link estimation is less stable than using the canonical (logit) link); (3) you can plug an inverse-link GLM in as a component of a mixed-model (GLMM) analysis, whereas doing this from scratch with a maximum likelihood model is a royal pain (see Chapter 10).

library(emdbook)
data(ReedfrogFuncresp)
with(ReedfrogFuncresp,plot(Initial,Killed/Initial))
## fit GLMs with three different links: logistic (default), log, inverse
m_logist <- glm(cbind(Killed,Initial-Killed)~ Initial, family="binomial",
    data=ReedfrogFuncresp)
m_log <- glm(cbind(Killed,Initial-Killed)~ Initial,
          family=binomial(link="log"),
          data=ReedfrogFuncresp)
m_inv <- glm(cbind(Killed,Initial-Killed)~ Initial,
          family=binomial(link=make.link("inverse")),
          data=ReedfrogFuncresp)
## inverse link is equivalent to 
library(bbmle)
mle2(Killed~dbinom(size=Initial,
                   prob=a/(1+a*h*Initial)),
                   start=list(a=0.25,h=0.02),
                   data=ReedfrogFuncresp)
## but faster, more stable, no need for starting conditions

n.b. starting parameters must be specified as a list, not a vector … (this may have changed since a previous version of bbmle?)

Approx. starting conditions come from asymptote (of number killed, not proportion killed) around 50 (hence h=1/50), initial slope about 50/200 (a=0.25).

Rogers random-predator equation: (could only do this with glm() if we could specify the inverse of this function, which would be ugly. The rogers.pred function gives the functional response curve — i.e. number, not proportion, killed per predator

rogers.pred <- function(N0,a,h,T) {
  N0 - lambertW(a*h*N0*exp(-a*(T-h*N0)))/(a*h)
}

This runs into trouble with default method="BFGS" (rogers.pred is a little bit fragile), so use Nelder-Mead instead:

m_lw <- mle2(Killed~dbinom(size=Initial,
                           prob=rogers.pred(Initial,
                             a,h,T=1)/Initial),
             start=list(a=0.25,h=0.02),
             method="Nelder-Mead",
          data=ReedfrogFuncresp)

Inverse-GLM (i.e. type II func response) is best, LW equivalent, but all fits are pretty much equivalent:

AICtab(m_logist,m_log,m_inv,m_lw,
        delta=TRUE,weights=TRUE,sort=TRUE)

It would normally make sense to consider corrected AIC for small sample sizes (AICc), but irrelevant since all models have equal complexity and AICc only affects the penalty term.

All different functions are more or less indistinguishable within the range of the data:

with(ReedfrogFuncresp,plot(Initial,Killed/Initial,xlim=c(0,1000)))
ivec <- 1:1000
newdat <- data.frame(Initial=ivec)
lines(ivec,predict(m_logist,newdata=newdat,type="response"))
lines(ivec,predict(m_log,newdata=newdat,type="response"),col=2)
lines(ivec,predict(m_inv,newdata=newdat,type="response"),col=4)
lines(ivec,predict(m_lw,newdata=newdat)/ivec,
      col=5)

Differences are somewhat clearer if we extrapolate:

with(ReedfrogFuncresp,plot(Initial,Killed/Initial,xlim=c(0,1000),
                           ylim=c(0,0.7)))
ivec <- 1:1000
newdat <- data.frame(Initial=ivec)
lines(ivec,predict(m_logist,newdata=newdat,type="response"))
lines(ivec,predict(m_log,newdata=newdat,type="response"),col=2)
lines(ivec,predict(m_inv,newdata=newdat,type="response"),col=4)
lines(ivec,predict(m_lw,newdata=newdat)/ivec,
      col=5)

The qualitative difference becomes even clearer if we look at number killed (initial number times mortality prob) rather than the mortality prob:

with(ReedfrogFuncresp,plot(Initial,Killed,xlim=range(ivec),
     ylim=c(0,60)))
lines(ivec,predict(m_logist,newdata=newdat,type="response")*ivec)
lines(ivec,predict(m_log,newdata=newdat,type="response")*ivec,col=2)
lines(ivec,predict(m_inv,newdata=newdat,type="response")*ivec,col=4)
lines(ivec,predict(m_lw,newdata=newdat),col=5)
lw1 <- coef(m_lw)
(inv1 <- with(as.list(coef(m_inv)),
     c(1/`(Intercept)`,Initial)))
clw <- confint(m_lw)
cinv <- confint(m_inv)
cinv[1,] <- rev(1/cinv[1,])
library(plotrix)
par(mfrow=c(1,2),bty="l",las=1)
plotCI(1:2,c(lw1[1],inv1[1]),
       li=c(clw[1,1],cinv[1,1]),
       ui=c(clw[1,2],cinv[1,2]),
       axes=FALSE,xlab="",ylab="attack rate",
       pch=16,xlim=c(0.75,2.25))
## ignore warnings
axis(side=2)
axis(side=1,at=1:2,label=c("RRP","H2"))
box()
## put axis label slightly farther from axis
par(mgp=c(3.5,1,0))
plotCI(1:2,c(lw1[2],inv1[2]),
       li=c(clw[2,1],cinv[2,1]),
       ui=c(clw[2,2],cinv[2,2]),
       axes=FALSE,xlab="",ylab="handling time",
       pch=16,xlim=c(0.75,2.25))
axis(side=2)
axis(side=1,at=1:2,label=c("RRP","H2"))
box()

Attack rate is somewhat higher and more uncertain under RRP than H2.

Or perhaps this is a better plot? Plot from zero, emphasize relative magnitudes:

## back to default mgp
par(mfrow=c(1,2),bty="l",las=1,mgp=c(3,1,0))
plotCI(1:2,c(lw1[1],inv1[1]),
       li=c(clw[1,1],cinv[1,1]),
       ui=c(clw[1,2],cinv[1,2]),
       axes=FALSE,xlab="",ylab="attack rate",
       pch=16,xlim=c(0.75,2.25),
       ylim=c(0,max(c(clw[1,2],cinv[1,2]))))
axis(side=2)
axis(side=1,at=1:2,label=c("RRP","H2"))
box()

par(mgp=c(3.5,1,0))
plotCI(1:2,c(lw1[2],inv1[2]),
       li=c(clw[2,1],cinv[2,1]),
       ui=c(clw[2,2],cinv[2,2]),
       axes=FALSE,xlab="",ylab="handling time",
       pch=16,xlim=c(0.75,2.25),
       ylim=c(0,max(c(clw[2,2],cinv[2,2]))))
axis(side=2)
axis(side=1,at=1:2,label=c("RRP","H2"))
box()
Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-ShareAlike 3.0 License