# 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