## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
has_lme4 <- requireNamespace("lme4", quietly = TRUE)
has_metafor <- requireNamespace("metafor", quietly = TRUE)
has_metadat <- requireNamespace("metadat", quietly = TRUE)
has_ordinal <- requireNamespace("ordinal", quietly = TRUE)
has_glmmTMB <- requireNamespace("glmmTMB", quietly = TRUE)
has_penguins <- requireNamespace("palmerpenguins", quietly = TRUE)

## -----------------------------------------------------------------------------
library(drmTMB)

## ----comparison1-drmtmb, eval = has_lme4--------------------------------------
data(cbpp, package = "lme4")
fit1 <- drmTMB(
  bf(mu = cbind(incidence, size - incidence) ~ period + (1 | herd)),
  data = cbpp, family = binomial()
)
coef(fit1, "mu")

## ----comparison1-comparator, eval = has_lme4----------------------------------
cmp1 <- lme4::glmer(
  cbind(incidence, size - incidence) ~ period + (1 | herd),
  family = binomial, data = cbpp
)
lme4::fixef(cmp1)

## ----comparison1-re, eval = has_lme4------------------------------------------
summary(fit1)$parameters[, c("term", "estimate")]
attr(lme4::VarCorr(cmp1)$herd, "stddev")

## ----comparison1-diff, eval = has_lme4----------------------------------------
coef(fit1, "mu") - lme4::fixef(cmp1)
c(drmTMB = as.numeric(logLik(fit1)), glmer = as.numeric(logLik(cmp1)))

## ----comparison2-data, eval = has_metafor && has_metadat----------------------
data(dat.bcg, package = "metadat")
dat_bcg <- metafor::escalc(
  measure = "RR", ai = tpos, bi = tneg, ci = cpos, di = cneg, data = dat.bcg
)

## ----comparison2-drmtmb, eval = has_metafor && has_metadat--------------------
fit2 <- drmTMB(
  bf(mu = yi ~ 1 + meta_V(V = vi), sigma = ~ 1),
  data = dat_bcg, family = gaussian()
)
coef(fit2, "mu")

## ----comparison2-reml, eval = has_metafor && has_metadat----------------------
cmp2_reml <- metafor::rma.uni(yi, vi, data = dat_bcg)
c(coef(cmp2_reml), tau2 = cmp2_reml$tau2)

## ----comparison2-ml, eval = has_metafor && has_metadat------------------------
cmp2 <- metafor::rma.uni(yi, vi, data = dat_bcg, method = "ML")
c(coef(cmp2), tau2 = cmp2$tau2)

## ----comparison2-convert, eval = has_metafor && has_metadat-------------------
c(drmTMB_tau2 = exp(2 * coef(fit2, "sigma")), metafor_tau2 = cmp2$tau2)
c(drmTMB = as.numeric(logLik(fit2)), metafor = as.numeric(logLik(cmp2)))

## ----comparison3-drmtmb, eval = has_metafor && has_metadat--------------------
fit3 <- drmTMB(
  bf(mu = yi ~ 1 + meta_V(V = vi), sigma = ~ alloc),
  data = dat_bcg, family = gaussian()
)
coef(fit3, "sigma")

## ----comparison3-comparator, eval = has_metafor && has_metadat----------------
dat_bcg$id <- seq_len(nrow(dat_bcg))
cmp3 <- metafor::rma.mv(
  yi, vi, random = ~ alloc | id, struct = "DIAG", data = dat_bcg, method = "ML"
)
cmp3$tau2

## ----comparison3-convert, eval = has_metafor && has_metadat-------------------
sigma3 <- coef(fit3, "sigma")
tau_by_level <- exp(c(
  alternate = sigma3[1],
  random = sigma3[1] + sigma3[2],
  systematic = sigma3[1] + sigma3[3]
))
tau_by_level
sqrt(cmp3$tau2)
c(drmTMB = as.numeric(logLik(fit3)), metafor = as.numeric(logLik(cmp3)))

## ----comparison4-drmtmb, eval = has_ordinal-----------------------------------
data(wine, package = "ordinal")
fit4 <- drmTMB(
  bf(mu = rating ~ temp + contact), data = wine, family = cumulative_logit()
)
summary(fit4)$coefficients
summary(fit4)$ordinal$cutpoints

## ----comparison4-comparator, eval = has_ordinal-------------------------------
cmp4 <- ordinal::clm(rating ~ temp + contact, data = wine)
coef(cmp4)

## ----comparison4-loglik, eval = has_ordinal-----------------------------------
c(drmTMB = as.numeric(logLik(fit4)), clm = as.numeric(logLik(cmp4)))

## ----comparison4-scale-comparator, eval = has_ordinal-------------------------
cmp4_scale <- ordinal::clm(rating ~ temp + contact, scale = ~ temp, data = wine)
as.numeric(logLik(cmp4_scale))

## ----comparison4-drmtmb-scale, eval = has_ordinal, error = TRUE---------------
try({
drmTMB(
  bf(mu = rating ~ temp + contact, sigma = ~ temp),
  data = wine, family = cumulative_logit()
)
})

## ----comparison5-drmtmb, eval = has_ordinal-----------------------------------
fit5 <- drmTMB(
  bf(mu = rating ~ temp + contact + (1 | judge)),
  data = wine, family = cumulative_logit()
)
summary(fit5)$coefficients
summary(fit5)$ordinal$cutpoints

## ----comparison5-comparator, eval = has_ordinal-------------------------------
cmp5 <- ordinal::clmm(rating ~ temp + contact + (1 | judge), data = wine)
coef(cmp5)

## ----comparison5-re, eval = has_ordinal---------------------------------------
summary(fit5)$parameters[, c("term", "estimate")]
attr(ordinal::VarCorr(cmp5)$judge, "stddev")
c(drmTMB = as.numeric(logLik(fit5)), clmm = as.numeric(logLik(cmp5)))

## ----comparison6-fe-check, eval = has_lme4------------------------------------
data(sleepstudy, package = "lme4")
fit6_fe <- drmTMB(
  bf(mu = Reaction ~ Days, sigma = ~ Days), data = sleepstudy, family = gaussian()
)
mu_hat <- predict_parameters(fit6_fe, dpar = "mu", type = "response")$estimate
eta_sigma <- predict_parameters(fit6_fe, dpar = "sigma", type = "link")$estimate
c(
  reported = as.numeric(logLik(fit6_fe)),
  `hand, sd = exp(eta)` = sum(dnorm(sleepstudy$Reaction, mu_hat, exp(eta_sigma), log = TRUE)),
  `hand, sd = sqrt(exp(eta))` = sum(dnorm(sleepstudy$Reaction, mu_hat, sqrt(exp(eta_sigma)), log = TRUE))
)

## ----comparison6-drmtmb, eval = has_lme4--------------------------------------
fit6 <- drmTMB(
  bf(mu = Reaction ~ Days + (1 + Days | Subject), sigma = ~ Days),
  data = sleepstudy, family = gaussian()
)
coef(fit6, "mu")
coef(fit6, "sigma")

## ----comparison6-comparator, eval = has_lme4 && has_glmmTMB-------------------
cmp6 <- glmmTMB::glmmTMB(
  Reaction ~ Days + (Days | Subject), dispformula = ~ Days, data = sleepstudy
)
glmmTMB::fixef(cmp6)$cond
glmmTMB::fixef(cmp6)$disp

## ----comparison6-loglik, eval = has_lme4 && has_glmmTMB-----------------------
c(drmTMB = as.numeric(logLik(fit6)), glmmTMB = as.numeric(logLik(cmp6)))
abs(coef(fit6, "sigma") - glmmTMB::fixef(cmp6)$disp)

## ----comparison7-drmtmb, eval = has_penguins && has_glmmTMB-------------------
pen <- palmerpenguins::penguins[
  stats::complete.cases(
    palmerpenguins::penguins[, c("species", "sex", "body_mass_g")]
  ),
]
fit7 <- drmTMB(
  bf(mu = body_mass_g ~ species, sigma = ~ sex),
  data = pen, family = drmTMB::lognormal()
)
coef(fit7, "mu")
coef(fit7, "sigma")

## ----comparison7-comparator, eval = has_penguins && has_glmmTMB---------------
cmp7 <- glmmTMB::glmmTMB(
  log(body_mass_g) ~ species, dispformula = ~ sex, data = pen, family = gaussian()
)
glmmTMB::fixef(cmp7)$cond
glmmTMB::fixef(cmp7)$disp

## ----comparison7-mu-check, eval = has_penguins && has_glmmTMB-----------------
mu_hat <- predict_parameters(fit7, dpar = "mu", type = "link")$estimate
eta_sigma <- predict_parameters(fit7, dpar = "sigma", type = "link")$estimate
y <- pen$body_mass_g
c(
  reported = as.numeric(logLik(fit7)),
  `hand, meanlog = mu` = sum(dlnorm(y, mu_hat, exp(eta_sigma), log = TRUE)),
  `hand, meanlog = mu - sigma^2/2` = sum(dlnorm(y, mu_hat - exp(eta_sigma)^2 / 2, exp(eta_sigma), log = TRUE))
)

## ----comparison7-loglik, eval = has_penguins && has_glmmTMB-------------------
c(drmTMB = as.numeric(logLik(fit7)), glmmTMB = as.numeric(logLik(cmp7)))
sum(-log(y))

## ----comparison8-fe-check, eval = has_glmmTMB---------------------------------
data(Owls, package = "glmmTMB")
fit8_fe <- drmTMB(
  bf(mu = SiblingNegotiation ~ FoodTreatment, sigma = ~ FoodTreatment),
  data = Owls, family = drmTMB::nbinom2()
)
mu_hat <- predict_parameters(fit8_fe, dpar = "mu", type = "response")$estimate
eta_sigma <- predict_parameters(fit8_fe, dpar = "sigma", type = "link")$estimate
y <- Owls$SiblingNegotiation
c(
  reported = as.numeric(logLik(fit8_fe)),
  `hand, size = 1/sigma^2` = sum(dnbinom(y, mu = mu_hat, size = 1 / exp(eta_sigma)^2, log = TRUE)),
  `hand, size = 1/sigma` = sum(dnbinom(y, mu = mu_hat, size = 1 / exp(eta_sigma), log = TRUE)),
  `hand, size = sigma` = sum(dnbinom(y, mu = mu_hat, size = exp(eta_sigma), log = TRUE))
)

## ----comparison8-clash, eval = has_glmmTMB, error = TRUE----------------------
try({
drmTMB(
  bf(mu = SiblingNegotiation ~ FoodTreatment, sigma = ~ FoodTreatment),
  data = Owls, family = glmmTMB::nbinom2()
)
})

## ----comparison8-drmtmb, eval = has_glmmTMB-----------------------------------
fit8 <- drmTMB(
  bf(
    mu = SiblingNegotiation ~ FoodTreatment * SexParent + (1 | Nest),
    sigma = ~ FoodTreatment
  ),
  data = Owls, family = drmTMB::nbinom2()
)
coef(fit8, "mu")
coef(fit8, "sigma")

## ----comparison8-comparator, eval = has_glmmTMB-------------------------------
cmp8 <- glmmTMB::glmmTMB(
  SiblingNegotiation ~ FoodTreatment * SexParent + (1 | Nest),
  dispformula = ~ FoodTreatment, family = glmmTMB::nbinom2, data = Owls
)
glmmTMB::fixef(cmp8)$cond
glmmTMB::fixef(cmp8)$disp

## ----comparison8-convert, eval = has_glmmTMB----------------------------------
c(drmTMB = as.numeric(logLik(fit8)), glmmTMB = as.numeric(logLik(cmp8)))
sigma8 <- coef(fit8, "sigma")
abs(-2 * sigma8 - glmmTMB::fixef(cmp8)$disp)

## ----reproducibility----------------------------------------------------------
pkgs <- c(
  "drmTMB", "TMB", "Matrix", "lme4", "metafor", "metadat", "ordinal", "glmmTMB",
  "palmerpenguins"
)
data.frame(
  package = pkgs,
  version = vapply(
    pkgs,
    function(p) {
      if (requireNamespace(p, quietly = TRUE)) {
        as.character(utils::packageVersion(p))
      } else {
        NA_character_
      }
    },
    character(1)
  ),
  row.names = NULL
)

## ----reproducibility-platform-------------------------------------------------
c(
  R = R.version.string,
  platform = R.version$platform,
  BLAS = basename(extSoftVersion()[["BLAS"]])
)