Nothing
## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(echo = FALSE, comment=NA, width=70)
options(show.signif.stars=FALSE)
## ----prepareData, echo=TRUE---------------------------------------------------
HawCon <- qra::HawCon
## Change name "CommonName" to "CN", for more compact output.
CCnum <- match("CommonName", names(HawCon))
names(HawCon)[CCnum] <- "CN"
## trtGp will identify species & lifestage combination
## trtGpRep will identify species, lifestage, and rep
## cTime is centered version of TrtTime
## scTime is centered and scaled version of TrtTime,
## needed to get some mixed model fits to converge
HawCon <- within(HawCon, {
trtGp <- factor(paste0(CN,LifestageTrt, sep=":"))
trtGpRep <- paste0(CN,LifestageTrt,":",RepNumber)
scTime <- scale(TrtTime)
obs <- factor(1:nrow(HawCon))
})
## ----cap1, echo=FALSE---------------------------------------------------------
cap1 <- " Graphs are designed to give an indication of the pattern,
when mortalities are shown on a complementary log-log scale, of
mortality response with days in coolstorage."
## ----plots, fig.width=7, fig.height=7.5, fig.align='center', out.width="75%", fig.cap=cap1----
qra::graphSum(df=HawCon, link="cloglog", logScale=FALSE,
dead="Dead", tot="Total", dosevar="TrtTime",
Rep="RepNumber", fitRep=NULL, fitPanel=NULL,
byFacet=~trtGp, layout=LifestageTrt~Species,
xlab="Days", maint="Hawaian contemporary data")
## ----echo=FALSE---------------------------------------------------------------
pkg <- "glmmTMB"
pcheck <- suppressWarnings(requireNamespace(pkg, quietly = TRUE))
if(pcheck) pcheck & packageVersion("glmmTMB") >= "1.1.2"
if(!pcheck){
message("`glmmTMB` version >= 1.1.2 is not available")
message("Will not continue execution of vignette")
knitr::knit_exit()
}
## ----glmmTMB-altFits, message=FALSE, warning=FALSE, echo=FALSE----------------
if("VGAM" %in% (.packages()))
detach("package:VGAM", unload=TRUE)
# Load packages that will be used
suppressMessages({library(lme4); library(splines)})
form <- cbind(Dead,Live)~0+trtGp/TrtTime+(1|trtGpRep)
form2s <- cbind(Dead,Live)~0+trtGp/TrtTime+ns(scTime,2)+(1|trtGpRep)
HCbb.cll <- glmmTMB::glmmTMB(form, dispformula=~trtGp+ns(scTime,2),
family=glmmTMB::betabinomial(link="cloglog"),
data=HawCon)
HCbb2s.cll <- update(HCbb.cll, formula=form2s)
HCbb.logit <- glmmTMB::glmmTMB(form, dispformula=~trtGp+ns(TrtTime,2),
family=glmmTMB::betabinomial(link="logit"),
data=HawCon)
HCbb2s.logit <- update(HCbb.logit, formula=form2s)
## ----glmmTMB-altFits, eval=FALSE, echo=-(1:2)---------------------------------
# if("VGAM" %in% (.packages()))
# detach("package:VGAM", unload=TRUE)
# # Load packages that will be used
# suppressMessages({library(lme4); library(splines)})
# form <- cbind(Dead,Live)~0+trtGp/TrtTime+(1|trtGpRep)
# form2s <- cbind(Dead,Live)~0+trtGp/TrtTime+ns(scTime,2)+(1|trtGpRep)
# HCbb.cll <- glmmTMB::glmmTMB(form, dispformula=~trtGp+ns(scTime,2),
# family=glmmTMB::betabinomial(link="cloglog"),
# data=HawCon)
# HCbb2s.cll <- update(HCbb.cll, formula=form2s)
# HCbb.logit <- glmmTMB::glmmTMB(form, dispformula=~trtGp+ns(TrtTime,2),
# family=glmmTMB::betabinomial(link="logit"),
# data=HawCon)
# HCbb2s.logit <- update(HCbb.logit, formula=form2s)
## ----cap7, echo=FALSE---------------------------------------------------------
cap7 <- "Differences are shown, between fitted degree 2 normal spline
curves and fittes lines.
Panel A is for the models that use a complementary log-log
(cloglog) link, while Panel B is for a logit link."
## ----glmmTMB-altFits-gph1, fig.width=9, fig.height=4.0, fig.show='hold', top=2, out.width="100%", fig.align='center', message=FALSE, warning=0, fig.cap=cap7, echo=FALSE----
library(lattice)
clog <- make.link('cloglog')$linkfun
logit <- make.link('logit')$linkfun
cloginv <- make.link('cloglog')$linkinv
logitinv <- make.link('logit')$linkinv
panel2 <- function(x, y, ...){
panel.superpose(x, y,
type='l', ...)
panel.xyplot(x, y, type='l', lwd=2, cex=0.6, ...)
}
parset <- simpleTheme(col=rep(1:4,rep(2,4)), lty=rep(1:2,4), lwd=2)
## c('solid','1141')
dat <- expand.grid(trtGp=factor(levels(HawCon$trtGp), levels=levels(HawCon$trtGp)),
TrtTime=pretty(range(HawCon$TrtTime),15), link=c('cloglog','logit'))
dat$scTime <- scale(dat$TrtTime)
dat$trtGpRep <- rep('new',nrow(dat))
hatClog <- predict(HCbb.cll, newdata=subset(dat, link=='cloglog'), se=TRUE, allow.new.levels=TRUE)
hatClog2 <- predict(HCbb2s.cll, newdata=subset(dat, link=='cloglog'), se=TRUE, allow.new.levels=TRUE)
diffClog <- cloginv(hatClog2$fit)-cloginv(hatClog$fit)
hatLgt <- predict(HCbb.logit, newdata=subset(dat, link=='logit'), se=TRUE, allow.new.levels=TRUE)
hatLgt2 <- predict(HCbb2s.logit, newdata=subset(dat, link=='logit'), se=TRUE, allow.new.levels=TRUE)
diffLgt <- logitinv(hatLgt2$fit)-logitinv(hatLgt$fit)
dat <- within(dat, {diff<-c(diffClog,diffLgt)
})
## dat <- within(dat, {lwr<-fit-2*se.fit; upr<-fit+2*se.fit})
gph <- xyplot(diff~TrtTime|link, outer=TRUE, data=dat, groups=trtGp,
# upper = dat$upr, lower = dat$lwr,
panel = panel2,
xlab="Days in coolstorage", ylab="Difference in fitted value",
auto.key=list(text=levels(HawCon$trtGp), columns=4, points=FALSE, lines=TRUE),
par.settings=parset, layout=c(2,1),
scales=list(x=list(at=c(2,6,10,14)),
y=list(relation='free'), alternating=c(1,1)),
strip=strip.custom(factor.levels=
c("A: Complementary log-log link", "B: Logit link")))
gph
## ----echo=FALSE---------------------------------------------------------------
pcheck <- suppressMessages(requireNamespace('DHARMa', quietly = TRUE))
yesDHARMa <- pcheck
if (!yesDHARMa) {
message("The suggested package `DHARMa` is not available/installed.")
message("Code that requires this package will not be executed.")
}
## ----cap3, echo=FALSE---------------------------------------------------------
cap3 <- "Panel A shows the quantile-quantile plot,
for the linear model with a complementary log-log link.
Panel B plots estimated quantiles against mortality,
while Panel C plots estimated quantiles against total
number, on a logarithmic scale."
## ----DHARMa, fig.width=3.5, fig.asp=0.95, bot=-1, top=1.5, out.width="32%", warning=0, fig.align='center', fig.show="hold", message=FALSE, echo=FALSE, fig.cap=cap3, eval=yesDHARMa----
oldpar <- par(mar=c(3.1,3.1,2.6,1.1), mgp=c(2.1, 0.5,0))
set.seed(29)
simRef <- DHARMa::simulateResiduals(HCbb.cll, n=250, seed=FALSE)
DHARMa::plotQQunif(simRef)
DHARMa::plotResiduals(simRef)
DHARMa::plotResiduals(simRef, form=log(HawCon$Total), xlab="log(Total)")
par(oldpar)
## ----cap4, echo=FALSE---------------------------------------------------------
cap4 <- "Diagnostic plots, for the model with a logit link."
## ----residBYgp, fig.width=3.75, fig.asp=0.95, bot=-1, top=1.5, out.width="32%", warning=0, fig.align='center', fig.show="hold", message=FALSE, echo=FALSE, mar=c(3.1,3.1,2.6,1.1), fig.cap=cap4,eval=yesDHARMa----
oldpar <- par(mar=c(3.1,3.1,2.6,1.1), mgp=c(2.1, 0.5,0))
simResBB <- suppressWarnings(DHARMa::simulateResiduals(HCbb.cll, n=250) )
DHARMa::plotQQunif(simResBB)
DHARMa::plotResiduals(simResBB, xlab= "Predictions, Complementary log-log model")
simResLGT <- suppressWarnings(DHARMa::simulateResiduals(HCbb.logit, n=250) )
DHARMa::plotResiduals(simResLGT, xlab= "Predictions, logit model")
par(oldpar)
## Alternative -- group names are not shown:
## plotResiduals(simRes, form=HawCon$trtGp)
## ----cap4.5, echo=FALSE-------------------------------------------------------
cap4.5 <- "Quantile residuals, by treatment group, for the
betabinomial model"
## ----trtGp, echo=FALSE, fig.width=5, fig.asp=0.7, bot=-2.5, warning=FALSE, fig.align='center', fig.show="hold", message=FALSE, warning=FALSE, out.width="49%", echo=FALSE, fig.cap=cap4.5, eval=yesDHARMa----
dotplot(scaledResiduals~HawCon$trtGp, data=simResBB,
scales=list(x=list(rot=30)), ylab="Quantile Residuals",
main=list(expression(plain("A: Residuals, by treatment group")),
x=0.05, y=-0.2, just="left"))
bwplot(scaledResiduals~HawCon$trtGp, data=simResBB, ylab="",
scales=list(x=list(rot=30)),
main=list(expression(plain("B: Boxplot comparison of residuals")),
x=0.05, y=-0.2, just="left"))
## ----cap4.75------------------------------------------------------------------
cap4.75 <- paste0("Diagnostics for model fitted to strongly overdispersed
binomial type data. Notice that the overdispersion results in an
S-shaped distribution of the residuals around the line $y=x$. The
boxplot is, in this case, uninformative.")
## ----overdispSim, echo=FALSE, fig.width=5, fig.asp=0.85, bot=-2.5, warning=FALSE, fig.align='center', fig.show="hold", out.width="49%", echo=FALSE, fig.cap=cap4.75, eval=yesDHARMa----
yes <- rbinom(n=100, size=50, prob=0.5)
sim <- cbind(yes=yes*4, no=200-yes*4)
sim.TMB <- glmmTMB::glmmTMB(sim~1, family=binomial)
sim250 <- DHARMa::simulateResiduals(sim.TMB)
DHARMa::plotQQunif(sim250)
DHARMa::plotResiduals(sim250)
## ----obslevel, echo=FALSE, eval=TRUE------------------------------------------
ctl <- glmmTMB::glmmTMBControl(optimizer=optim,
optArgs=list(method="BFGS"))
HCbiObs.cll <-
glmmTMB::glmmTMB(cbind(Dead, Live) ~ 0 + trtGp/scTime +
(1|obs) + (scTime|trtGpRep),
family=binomial(link='cloglog'),
control=ctl, data=HawCon)
HCbiObs.logit <-
glmmTMB::glmmTMB(cbind(Dead, Live) ~ 0 + trtGp/scTime +
(1|obs) + (scTime|trtGpRep),
family=binomial(link='logit'),
control=ctl, data=HawCon)
## ----glmmTMB-altFits-AIC, echo=FALSE------------------------------------------
aicStats <-
AIC(HCbb.cll,HCbb2s.cll,HCbb.logit, HCbb2s.logit,
HCbiObs.cll, HCbiObs.logit)
rownames(aicStats) <- substring(rownames(aicStats),3)
round(t(aicStats),2)
## ----cap8, echo=FALSE---------------------------------------------------------
cap8 <- "Panels A and B show intra-class correlation estimates
for, respectively, a complementary log-log link and a logit link.
Both models assume a betabinomial error. Panel C shows, for
the complementary log-log model, the dispersion factors that
result."
## ----glmmTMB-altFits-gph-disp, fig.width=9, fig.height=3.5, top=2, out.width="100%", fig.align='center', fig.show="hold", fig.cap=cap8, echo=FALSE----
oldpar <- par(oma=c(0,0,2,0))
parset <- simpleTheme(col=rep(1:4,rep(2,4)),pch=rep(c(1,2), 4), lwd=2)
HawCon$rho2clog <- qra::getRho(HCbb.cll)
HawCon$dispClog <- with(HawCon, 1+(Total-1)*rho2clog)
titles=c(expression("A: "*rho*", cloglog link"),expression("B: "*rho*", logit link"),
"C: Dispersion, cloglog link")
library(lattice)
HawCon$rho2logit <- qra::getRho(HCbb.logit)
xyplot(rho2clog+rho2logit+dispClog ~ TrtTime, groups=trtGp, data=HawCon,
outer=TRUE, between=list(x=0.25),
par.settings=parset,
scales=list(x=list(alternating=FALSE), y=list(relation='free',tick.number=4)),
auto.key=list(columns=4, between.columns=2, between=1),
xlab="Days in coolstorage", ylab="Parameter Estimate",
strip=strip.custom(factor.levels=titles))
par(oldpar)
## ----obslevel, echo=FALSE, eval=TRUE------------------------------------------
ctl <- glmmTMB::glmmTMBControl(optimizer=optim,
optArgs=list(method="BFGS"))
HCbiObs.cll <-
glmmTMB::glmmTMB(cbind(Dead, Live) ~ 0 + trtGp/scTime +
(1|obs) + (scTime|trtGpRep),
family=binomial(link='cloglog'),
control=ctl, data=HawCon)
HCbiObs.logit <-
glmmTMB::glmmTMB(cbind(Dead, Live) ~ 0 + trtGp/scTime +
(1|obs) + (scTime|trtGpRep),
family=binomial(link='logit'),
control=ctl, data=HawCon)
## ----cap2---------------------------------------------------------------------
cap2 <- "AICs are compared between models with binomial family
errors and observation level random effects. The two sets of
points make the comparison for, respectively, data that have
been simulated from a model with logit link and a model with
complementary log-log link. The large triangle makes the
comparison for the models fitted to the 'HawCon' data."
## ----only-obslevel, fig.width=3.75, fig.asp=0.85, bot=-1, top=1.5, out.width="55%", fig.align='center', fig.show="hold", message=FALSE, echo=FALSE, mar=c(3.1,3.1,2.6,1.1), fig.cap = cap2----
## Notice the use of the 'BFGS' optimization method in place
## of the default.
aicData <- setNames(AIC(HCbiObs.logit,HCbiObs.cll)[,2], c('logit','cll'))
set.seed(17)
sim <- simulate(HCbiObs.logit, nsim=10)
simcll <- simulate(HCbiObs.cll, nsim=10)
aic.cll <- aic2.cll <- aic.logit <- aic2.logit <- numeric(10)
for (i in 1:10){
zlogit <-
glmmTMB::glmmTMB(sim[[i]] ~ 0 + trtGp/scTime +
(1|obs) + (1|trtGpRep),
family=binomial(link='logit'),
data=HawCon)
aic.logit[i] <- AIC(zlogit)
zcll <-
glmmTMB::glmmTMB(sim[[i]] ~ 0 + trtGp/scTime +
(1|obs) + (1|trtGpRep),
family=binomial(link='cloglog'),
data=HawCon)
aic.cll[i] <- AIC(zcll)
zlogit2 <-
glmmTMB::glmmTMB(simcll[[i]] ~ 0 + trtGp/scTime +
(1|obs) + (1|trtGpRep),
family=binomial(link='logit'),
data=HawCon)
aic2.logit[i] <- AIC(zlogit2)
zcll2 <-
glmmTMB::glmmTMB(simcll[[i]] ~ 0 + trtGp/scTime +
(1|obs) + (1|trtGpRep),
family=binomial(link='cloglog'),
data=HawCon)
aic2.cll[i] <- AIC(zcll2)
}
gph1 <- xyplot(aic.cll~aic.logit,
key=list(text=list(c("logit model","cll model", "Data")),
points=list(pch=c(1,1,2),
col=c('blue','red','black')),columns=3))
gph2 <- xyplot(aic2.cll~aic2.logit, col='red')
gph3 <- gph1+latticeExtra::as.layer(gph2)+
latticeExtra::layer(panel.abline(0,1),
panel.points(aicData[['logit']],
aicData[['cll']], pch=2, cex=2, col=1))
update(gph3, xlim=range(c(aic.logit, aic2.logit,
aicData['logit']))*c(.98,1.02),
ylim=range(c(aic.cll, aic2.cll,
aicData['cll']))*c(.98,1.02))
## ----cap5, echo=FALSE---------------------------------------------------------
cap5 <- "Diagnostics for model with binomial errors and observation level
random effects."
## ----biObsCLL, fig.width=3.75, fig.asp=0.95, bot=-1, top=2.5, out.width="49%", fig.align='center', fig.show="hold", message=FALSE, warning=0, echo=FALSE, mar=c(3.1,3.1,2.6,1.1), fig.cap=cap5, eval=yesDHARMa----
set.seed(29)
simRefcll <- suppressWarnings(
DHARMa::simulateResiduals(HCbiObs.cll, n=250, seed=FALSE) )
DHARMa::plotResiduals(simRefcll, xlab='cll: model prediction')
DHARMa::plotResiduals(simRefcll, form=log(HawCon$Total),
xlab="cll: log(Total)")
simReflogit <- suppressWarnings(
DHARMa::simulateResiduals(HCbiObs.logit, n=250, seed=FALSE) )
DHARMa::plotResiduals(simReflogit, xlab='logit: model prediction')
DHARMa::plotResiduals(simReflogit, form=log(HawCon$Total),
xlab="logit: log(Total)")
## ----biObsCLL, eval=FALSE, echo=TRUE------------------------------------------
# set.seed(29)
# simRefcll <- suppressWarnings(
# DHARMa::simulateResiduals(HCbiObs.cll, n=250, seed=FALSE) )
# DHARMa::plotResiduals(simRefcll, xlab='cll: model prediction')
# DHARMa::plotResiduals(simRefcll, form=log(HawCon$Total),
# xlab="cll: log(Total)")
# simReflogit <- suppressWarnings(
# DHARMa::simulateResiduals(HCbiObs.logit, n=250, seed=FALSE) )
# DHARMa::plotResiduals(simReflogit, xlab='logit: model prediction')
# DHARMa::plotResiduals(simReflogit, form=log(HawCon$Total),
# xlab="logit: log(Total)")
## ----shorten, echo=FALSE------------------------------------------------------
shorten <- function(nam, leaveout=c('trtGp','Fly',':')){
for(txt in leaveout){
nam <- gsub(txt,'', nam, fixed=TRUE)
}
nam
}
## ----crude-cll, echo=FALSE, warning=F-----------------------------------------
## Fit two simplistic and unsatisfactory models.
HCbbDisp1.cll <- update(HCbb.cll, dispformula=~1)
HCbin.cll <- update(HCbb.cll, family=binomial(link="cloglog"))
## ----extract-BB-LTcll, echo=FALSE---------------------------------------------
LTbb.cll <- qra::extractLT(p=0.99, obj=HCbb.cll, link="cloglog",
a=1:8, b=9:16, eps=0, df.t=NULL)[,-2]
rownames(LTbb.cll) <- shorten(rownames(LTbb.cll))
## ----extract-BI-obsRE, echo=FALSE---------------------------------------------
## NB: The formula has used the scaled value of time.
## Below, `offset` is used to retrieve the scaling parameters
## `center` ## and `scale` in `(x-center)/scale`.
offset <- qra::getScaleCoef(HawCon$scTime)
LTbiObs.cll <- qra::extractLT(p=0.99, obj=HCbiObs.cll,
a=1:8, b=9:16, eps=0, offset=offset,
df.t=NULL)[,-2]
rownames(LTbiObs.cll) <- shorten(rownames(LTbiObs.cll))
## ----extract-BB-LTlogit, echo=FALSE-------------------------------------------
LTbb.logit <- qra::extractLT(p=0.99, obj=HCbb.logit, link="logit",
a=1:8, b=9:16, eps=0, offset=0,
df.t=NULL)[,-2]
rownames(LTbb.logit) <- shorten(rownames(LTbb.logit))
## ----extract-crude-LTcll, echo=FALSE------------------------------------------
LTbbDisp1.cll <-
qra::extractLT(p=0.99, obj=HCbbDisp1.cll,
a=1:8, b=9:16, eps=0, df.t=NULL)[,-2]
rownames(LTbbDisp1.cll) <- shorten(rownames(LTbbDisp1.cll))
LTbin.cll <-
qra::extractLT(p=0.99, obj=HCbin.cll,
a=1:8, b=9:16, eps=0, df.t=NULL)[,-2]
rownames(LTbin.cll) <- shorten(rownames(LTbin.cll))
## ----OKplotrix-plotCI---------------------------------------------------------
OKplotrix <- suppressWarnings(requireNamespace("plotrix", quietly = TRUE))
if (!OKplotrix) {
message("The `plotix' package is not available/installed.")
message("Code that requires the `plotix' package will not be executed.")
}
## ----cap9, echo=FALSE---------------------------------------------------------
cap9 <- "LT99 $95\\%$ confidence intervals are compared between
five different models."
## ----plotCI, echo=FALSE, fig.width=7.0, bot=1.0, fig.asp=0.625, warning=FALSE, fig.align='center', message=FALSE, out.width="75%", echo=FALSE, fig.cap=cap9, eval=OKplotrix----
gpNam <- rownames(LTbb.cll)
ordEst <- order(LTbb.cll[,1])
col5 <- c("blue","lightslateblue","brown",'gray','gray50')
plotrix::plotCI(1:8-0.34, y=LTbb.cll[ordEst,1], ui=LTbb.cll[ordEst,3],
li=LTbb.cll[ordEst,2], lwd=2, col=col5[1], xaxt="n",
fg="gray", xlab="", ylab="LT99 Estimate (days)",
xlim=c(0.8,8.2), ylim=c(0,29), sfrac=0.008)
plotrix::plotCI(1:8-0.17, y=LTbiObs.cll[ordEst,1], ui=LTbiObs.cll[ordEst,3],
li=LTbiObs.cll[ordEst,2], lwd=2, col=col5[2], xaxt="n",
fg="gray", xlab="", ylab="LT99 Estimate (days)",
xlim=c(0.8,8.2), ylim=c(0,29), add=TRUE, sfrac=0.008)
plotrix::plotCI(1:8, y=LTbb.logit[ordEst,1], ui=LTbb.logit[ordEst,3],
li=LTbb.logit[ordEst,2], lwd=2, col=col5[3], xaxt="n",
add=TRUE, sfrac=0.008)
plotrix::plotCI(1:8+0.17, y=LTbbDisp1.cll[ordEst,1], ui=LTbbDisp1.cll[ordEst,3],
li=LTbbDisp1.cll[ordEst,2], lwd=2, col=col5[4], xaxt="n",
add=TRUE, sfrac=0.008)
plotrix::plotCI(1:8+0.34, y=LTbin.cll[ordEst,1], ui=LTbin.cll[ordEst,3],
li=LTbin.cll[ordEst,2], lwd=2, col=col5[5], xaxt="n",
add=TRUE, sfrac=0.008)
axis(1, at=1:8, labels=gpNam[ordEst], las=2, lwd=0,
lwd.ticks=0.5, col.ticks="gray")
legend("topleft", legend=c("BB-cll (cll=cloglog)", "BB-cll-ObsRE", "BB-logit",
"BB-cll, const Disp factor",
"Binomial-cll"),
inset=c(0.01,0.01), lty=c(1,1), col=col5[1:5],
text.col=col5[1:5], bty="n",y.intersp=0.85)
## ----extract-BB-LTcll, eval=FALSE, echo=TRUE----------------------------------
# LTbb.cll <- qra::extractLT(p=0.99, obj=HCbb.cll, link="cloglog",
# a=1:8, b=9:16, eps=0, df.t=NULL)[,-2]
# rownames(LTbb.cll) <- shorten(rownames(LTbb.cll))
## ----obsLevel1, echo=FALSE, eval=FALSE----------------------------------------
# dMedEgg <- with(HawCon, dummy(trtGp,"MedFlyEgg:"))
# dMedL1 <- with(HawCon, dummy(trtGp,"MedFlyL1:"))
# dMedL2 <- with(HawCon, dummy(trtGp,"MedFlyL2:"))
# dMedL3 <- with(HawCon, dummy(trtGp,"MedFlyL3:"))
# dMelEgg <- with(HawCon, dummy(trtGp,"MelonFlyEgg:"))
# dMelL1 <- with(HawCon, dummy(trtGp,"MelonFlyL1:"))
# dMelL2 <- with(HawCon, dummy(trtGp,"MelonFlyL2:"))
# dMelL3 <- with(HawCon, dummy(trtGp,"MelonFlyL3:"))
# formXre <- cbind(Dead, Live) ~ 0 + trtGp/scTime +
# (1|obs) + (1|trtGpRep) +
# (0 + dMedEgg| trtGpRep) + (0 + dMedL1 | trtGpRep) +
# (0 + dMedL2 | trtGpRep) + (0 + dMedL3 | trtGpRep) +
# (0 + dMelEgg| trtGpRep) + (0 + dMelL1 | trtGpRep) +
# (0 + dMelL2 | trtGpRep) + (0 + dMelL3 | trtGpRep)
# ## NB: The formula has used the scaled value of time.
# ## Below, `offset` is used to record the scaling parameters
# ## `center` ## and `scale` in `(x-center)/scale`.
# offset <- with(attributes(HawCon$scTime),
# c(`scaled:center`, `scaled:scale`))
# HCXre.biObs <- gkmmTMB::glmmTMB(formXre, family=binomial(link='cloglog'),
# control=ctl, data=HawCon)
# ## Could not get the following to converge
# ## formXreS <- update(formXre, ~ .+ trtGp/splines::ns(scTime,2))
## ----glmerFits----------------------------------------------------------------
## Comparisons using `glmer()`
HCglmerBIobs.cll <-
lme4::glmer(cbind(Dead, Live) ~ 0 + trtGp/scTime + (1|obs) + (1|trtGpRep),
family=binomial(link='cloglog'), nAGQ=1, data=HawCon,
control=lme4::glmerControl(optimizer='bobyqa'))
HCglmerBIobs2s.cll <- suppressWarnings(
lme4::glmer(cbind(Dead, Live) ~ 0 + trtGp/scTime + splines::ns(scTime,2) +
(1|obs) + (1|trtGpRep), family=binomial(link='cloglog'),
nAGQ=0, data=HawCon))
HCglmerBIobs.logit <- lme4::glmer(cbind(Dead, Live) ~ 0 + trtGp/scTime +
(1|obs) + (1|trtGpRep),
family=binomial(link='logit'), nAGQ=0, data=HawCon,
control=glmerControl(optimizer='bobyqa'))
HCglmerBIobs2s.logit <-
lme4::glmer(cbind(Dead, Live) ~ 0 + trtGp/scTime + splines::ns(scTime,2) +
(1|obs) + (1|trtGpRep), family=binomial(link='logit'), nAGQ=0,
data=HawCon, control=glmerControl(optimizer='bobyqa'))
## ----cfAICs-------------------------------------------------------------------
cfAIC <-
AIC(HCbiObs.cll,HCglmerBIobs.cll,HCbiObs.logit,HCglmerBIobs.logit,
HCglmerBIobs2s.cll,HCglmerBIobs2s.logit)
rownames(cfAIC) <- c('TMB:cll','mer:cll','TMB:lgt','mer:lgt',
'mer:cllCurve', 'mer:lgtCurve')
(tab <- t(round(cfAIC,2)))
## ----allFit, results='hide', echo=TRUE, eval=FALSE----------------------------
# check <- (requireNamespace('dfoptim',quietly=TRUE)&
# requireNamespace('optimx',quietly=TRUE))
# if(check)
# ss<-suppressWarnings(summary(allFit(HCglmerBIobs.cll)))
## ----names-ss, hold=FALSE, echo=-2, eval=FALSE--------------------------------
# stopifnot(check)
# options(width=70)
# names(ss)
# ss$msgs
# ss$llik
## ----LT99gauss.LTcll----------------------------------------------------------
cloglog <- make.link('cloglog')$linkfun
HCgauss.cll <- glmmTMB::glmmTMB(cloglog((PropDead+0.002)/1.004)~0+
trtGp/TrtTime+(TrtTime|trtGpRep),
family=gaussian(), data=HawCon)
LTgauss.cll <- qra::extractLT(p=0.99, obj=HCgauss.cll, link="cloglog",
a=1:8, b=9:16, eps=0.002, offset=c(0,1), df.t=NULL)[,-2]
rownames(LTgauss.cll) <- shorten(rownames(LTgauss.cll))
## ----cap11, echo=FALSE--------------------------------------------------------
cap11 <- "Residuals versus predicted quantile deviations, for
the linear mixed model,
with \\(log(1-log((p+0.002)/(1+0.004)))\\) as the dependent
variable, complementary log-log link, and gaussian error."
## ----Gauss-simRes, echo=FALSE, w=3.5, fig.asp=0.75, left=-0.5, bot=-1, mgp=c(3,1,0), crop=TRUE, fig.align='center', out.width="50%", fig.cap=cap11, eval=yesDHARMa----
sim <- DHARMa::simulateResiduals(HCgauss.cll)
DHARMa::plotResiduals(sim)
## ----cap12, echo=FALSE--------------------------------------------------------
cap12 <- "Comparison of estimates and upper and lower $95\\%$
confidence limits, between the preferred betabinomial
complementary log-log model (BB) and the gaussian error model."
## ----BBvsGauss, out.width="100%", message=FALSE, echo=FALSE------------------
cfLTs <- cbind(LTbb.cll, LTgauss.cll)
colnames(cfLTs) <- c(rep('BB',3),rep('gauss',3))
tab <- round(cfLTs[,c(c(1,4),c(1,4)+1,c(1,4)+2)],1)
pcheck <- suppressWarnings(requireNamespace("kableExtra", quietly = TRUE))
if(pcheck) pcheck & packageVersion("kableExtra") >= "1.2"
if(pcheck){
kableExtra::kbl(tab, font_size=9, caption=cap12) |>
kableExtra::kable_paper("striped", full_width=FALSE) |>
kableExtra::column_spec(6:7, bold=TRUE) |>
kableExtra::row_spec(5:8, color='blue', bold=T) |>
kableExtra::add_header_above(c(' '=1,'Estimate'=2,'Lower'=2, 'Upper'=2), align='c')
} else print(tab)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.