Plot MeBr Apple Fumigation Data

John Maindonald

22 May 2025

library(qra, quietly=TRUE)
## [1] TRUE

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")

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)
##        phi nmin nmax Phimin Phimax PhiGLM
## 88 9697.52 1067 3277   1.11   1.34   2.33
## 89  447.54 1474 4177   4.29  10.33  12.30

The estimates of \(\phi\) (sigma) are 1988: 1.03^{-4}; and 1989: 0.00223 . 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)
##                Min   Max
## Gala          1.03  1.05
## Red Delicious 6.08 10.09
## Splendour     8.07 13.70

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)
}
## Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
## problem; . See vignette('troubleshooting'), help('diagnose')
round(mat,1)
##                    phi nmin nmax PhiMin PhiMax
## ROYAL         577332.6 1597 1676    1.0    1.0
## BRAEBURN      649760.1 1662 2123    1.0    1.0
## FUJI          804367.4 1392 1773    1.0    1.0
## GRANNY          2778.0 2284 3277    1.8    2.2
## Red Delicious    592.6 1067 1147    2.8    2.9