In jhmaindonald/qra: Quantal Response Analysis for Dose-Mortality Data

library(qra, quietly=TRUE)

pkg <- "glmmTMB"
pcheck <- suppressWarnings(requireNamespace(pkg, quietly = TRUE))
if(pcheck) pcheck & packageVersion("glmmTMB") >= "1.1.2"
noglmmTMB <- !pcheck

Plot data

Plot 1988 and 1989 data separately

qra::graphSum(df=qra::codling1988, link="cloglog", logScale=FALSE,
                     dead="dead", tot="total", dosevar="ct", Rep="rep",
                     fitRep=NULL, fitPanel=NULL,
                     byFacet=~Cultivar,
                     maint="1988: Codling moth, MeBr",
                     xlab=expression(bold("CT ")*"(gm.h."*m^{-3}*")"))

qra::graphSum(df=qra::codling1989, link="cloglog", logScale=FALSE,
                     dead="dead", tot="total", dosevar="ct", Rep="rep",
                     fitRep=NULL, fitPanel=NULL,
                     byFacet=~Cultivar,
                     maint="1989: Codling moth, MeBr",
                     xlab=expression(bold("CT ")*"(gm.h."*m^{-3}*")"))

Correct for control mortality

library(lattice)
cloglog <- make.link("cloglog")$linkfun
cod88 <- subset(qra::codling1988, dose>0)
cm88 <- subset(qra::codling1988, dose==0)
cmMatch <- match(cod88$cultRep,cm88$cultRep)
cod88$cm <- cm88[cmMatch,'PropDead']
cod88$apobs <- with(cod88, (PropDead-cm)/(1-cm))
xyplot(cloglog(apobs)~ct|Cultivar, groups=rep, 
       data=subset(cod88, apobs>0), layout=c(5,1),
       maint="1988: Codling moth, MeBr")

cod89 <- subset(qra::codling1989, dose>0)
cm89 <- subset(qra::codling1989, dose==0)
cmMatch <- match(cod89$cultRep,cm89$cultRep)
cod89$cm <- cm89[cmMatch,'PropDead']
cod89$apobs <- with(cod89, (PropDead-cm)/(1-cm))
xyplot(cloglog(apobs)~ct|Cultivar, groups=rep, data=cod89, layout=c(3,1),
       maint="1989: Codling moth, MeBr")

if (noglmmTMB) {
  message(" This vignette requires a version of `glmmTMB` >= 1.1.2")
  message(" Earlier versions of `glmmTMB` may have issues\n for matching versions of `TMB` and `Matrix`")
knit_exit()
}

We now check the consistency of control mortality estimates across and within cultivars.

ctl <- glmmTMB::glmmTMBControl(optimizer=optim,
                      optArgs=list(method="BFGS"))
cm88.TMB <- glmmTMB::glmmTMB(cbind(dead,total-dead)~Cultivar,
               family=glmmTMB::betabinomial(link='logit'), 
               data=cm88)
cm88.glm <- glm(cbind(dead,total-dead)~Cultivar, 
               family=quasibinomial(link='logit'), 
               data=cm88)
cm89.TMB <- glmmTMB::glmmTMB(cbind(dead,total-dead)~Cultivar, 
                   family=glmmTMB::betabinomial(link='logit'), 
                   data=cm89)
cm89.glm <- glm(cbind(dead,total-dead)~Cultivar, 
                family=quasibinomial(link='logit'), 
                data=cm89)

The following shows the range of estimated variance multipliers for the fit with a betabinomial error, and compares it with the "dispersion" estimate for a glm() fit with quasibinomial error.

mults <- matrix(nrow=2, ncol=6)
dimnames(mults) <- list(c('88','89'),c('phi','nmin','nmax',
                                       'Phimin','Phimax','PhiGLM'))
phi88 <- glmmTMB::sigma(cm88.TMB)
nrange <- range(cm88$total)
Phirange <- 1+phi88^{-1}*(nrange-1)
mults['88',] <- c(phi88, nrange, Phirange, 
                  summary(cm88.glm)$dispersion)
phi89 <- sigma(cm89.TMB)
nrange <- range(cm89$total)
Phirange <- 1+phi89^{-1}*(nrange-1)
mults['89',] <- c(phi89, nrange, Phirange, 
                  summary(cm89.glm)$dispersion)
round(mults,2)

The estimates of $\phi$ (sigma) are 1988: r round(phi88^-1,6); and 1989: r round(phi89^-1,5) . Although very small, multiplication by values of $n$ that are in the thousands lead to a variance multiplier that is noticeably greater than 1.0 on 1988, and substantial in 1989.

Now fit a model for the 1989 control mortality data that allows the betabinomial dispersion parameter to be different for the different cultivars. The dispformula changes from the default dispformula = ~1 to dispformula = ~0+Cultivar. (Use of dispformula = ~Cultivar, equivalent to dispformula = ~1+Cultivar, would give a parameterization that is less convenient for present purposes.)

cm89d.TMB <- update(cm89.TMB, dispformula=~0+Cultivar,    
                   control=ctl)
nranges <- with(cm89,
                sapply(split(total,Cultivar),range))
phi89d <- exp(coef(summary(cm89d.TMB))$disp[,1])
Phiranges <- 1+t(nranges-1)*phi89d^-1
colnames(Phiranges) <- c("Min","Max")
round(Phiranges,2)

Because the model formula allowed different values for different observations, the function sigma() could no longer be used to extract what are now three different dispersion parameters. Instead, it was necessary to specify coef(summary(cm89d.TMB))$disp[,1], and then (because a logarithmic link function is used) take the exponent. (Use of coef(summary(cm89d.TMB))$disp gives, on the logarithmic link scale, standard errors, z values, and $p$-values.)

cults <- unique(cm88$Cultivar)
mat <- matrix(nrow=length(cults),ncol=5)
dimnames(mat) <- list(cults,
                      c("phi","nmin","nmax","PhiMin","PhiMax"))
for(cult in cults){
df <- subset(cm88, cult==Cultivar)
obj <- glmmTMB::glmmTMB(cbind(dead,total-dead)~1, control=ctl,      
               family=glmmTMB::betabinomial(link='logit'), data=df)
phi <- glmmTMB::sigma(obj)
nrange <- range(df$total)
mat[cult,] <- c(phi, nrange, 1+(nrange-1)*phi^-1)
}
round(mat,1)