## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")

## ----gaussian-ls-syntax, eval = FALSE-----------------------------------------
# fit <- drmTMB(
#   drm_formula(
#     y ~ x,
#     sigma ~ x
#   ),
#   family = gaussian(),
#   data = dat
# )

## ----lme4-comparator-syntax, eval = FALSE-------------------------------------
# fit <- drmTMB(
#   drm_formula(y ~ x + (1 | id)),
#   family = gaussian(),
#   data = dat
# )
# 
# fit_lme4 <- lme4::lmer(
#   y ~ x + (1 | id),
#   data = dat,
#   REML = FALSE
# )

## ----dense-meta-syntax, eval = FALSE------------------------------------------
# fit <- drmTMB(
#   drm_formula(yi ~ x + meta_V(V = V)),
#   family = gaussian(),
#   data = dat
# )

## ----metafor-comparator-syntax, eval = FALSE----------------------------------
# fit_metafor <- metafor::rma.mv(
#   yi = yi,
#   V = V,
#   mods = ~ x,
#   random = ~ 1 | obs,
#   data = dat,
#   method = "ML"
# )

## ----metafor-comparator-checks, eval = FALSE----------------------------------
# expect_equal(unname(coef(fit, "mu")),
#              unname(stats::coef(fit_metafor)),
#              tolerance = 1e-4)
# expect_equal(stats::sigma(fit)[[1L]]^2,
#              fit_metafor$sigma2[[1L]],
#              tolerance = 1e-4)
# expect_equal(as.numeric(stats::logLik(fit)),
#              as.numeric(stats::logLik(fit_metafor)),
#              tolerance = 1e-4)

## ----mass-nb-comparator-syntax, eval = FALSE----------------------------------
# fit <- drmTMB(
#   drm_formula(count ~ x, sigma ~ 1),
#   family = nbinom2(),
#   data = dat
# )
# 
# fit_mass <- MASS::glm.nb(count ~ x, data = dat)

## ----mass-nb-comparator-checks, eval = FALSE----------------------------------
# expect_equal(unname(coef(fit, "mu")),
#              unname(stats::coef(fit_mass)),
#              tolerance = 1e-4)
# expect_equal(sigma(fit)[[1L]],
#              1 / sqrt(fit_mass$theta),
#              tolerance = 1e-4)
# expect_equal(as.numeric(stats::logLik(fit)),
#              as.numeric(stats::logLik(fit_mass)),
#              tolerance = 1e-4)

## ----beta-comparator-checks, eval = FALSE-------------------------------------
# eta_mu <- as.vector(fit$model$X$mu %*% coef(fit, "mu"))
# eta_sigma <- as.vector(fit$model$X$sigma %*% coef(fit, "sigma"))
# mu <- plogis(eta_mu)
# sigma <- exp(eta_sigma)
# phi <- 1 / sigma^2
# ll_independent <- sum(stats::dbeta(
#   fit$model$y,
#   shape1 = mu * phi,
#   shape2 = (1 - mu) * phi,
#   log = TRUE
# ))
# expect_equal(as.numeric(logLik(fit)), ll_independent, tolerance = 1e-6)

## ----student-syntax, eval = FALSE---------------------------------------------
# fit <- drmTMB(
#   drm_formula(
#     y ~ x,
#     sigma ~ x,
#     nu ~ x
#   ),
#   family = student(),
#   data = dat
# )

## ----independent-student-likelihood, eval = FALSE-----------------------------
# eta_mu <- X_mu %*% beta_mu
# eta_sigma <- X_sigma %*% beta_sigma
# eta_nu <- X_nu %*% beta_nu
# 
# mu <- as.vector(eta_mu)
# sigma <- exp(as.vector(eta_sigma))
# nu <- 2 + exp(as.vector(eta_nu))
# 
# loglik <- sum(
#   stats::dt((y - mu) / sigma, df = nu, log = TRUE) - log(sigma)
# )

## ----independent-truncated-nbinom2-likelihood, eval = FALSE-------------------
# eta_mu <- X_mu %*% beta_mu
# eta_sigma <- X_sigma %*% beta_sigma
# 
# mu <- exp(as.vector(eta_mu))
# sigma <- exp(as.vector(eta_sigma))
# size <- 1 / sigma^2
# 
# log_p0 <- stats::dnbinom(0, mu = mu, size = size, log = TRUE)
# log_positive <- log1p(-exp(log_p0))
# loglik <- sum(
#   stats::dnbinom(y, mu = mu, size = size, log = TRUE) - log_positive
# )

## ----independent-hurdle-nbinom2-likelihood, eval = FALSE----------------------
# eta_mu <- X_mu %*% beta_mu
# eta_sigma <- X_sigma %*% beta_sigma
# eta_hu <- X_hu %*% beta_hu
# 
# mu <- exp(as.vector(eta_mu))
# sigma <- exp(as.vector(eta_sigma))
# hu <- plogis(as.vector(eta_hu))
# size <- 1 / sigma^2
# 
# log_p0 <- stats::dnbinom(0, mu = mu, size = size, log = TRUE)
# log_positive <- log1p(-exp(log_p0))
# loglik <- sum(ifelse(
#   y == 0,
#   log(hu),
#   log1p(-hu) + stats::dnbinom(y, mu = mu, size = size, log = TRUE) -
#     log_positive
# ))

## ----bivariate-rho12-syntax, eval = FALSE-------------------------------------
# fit <- drmTMB(
#   drm_formula(
#     mu1 = y1 ~ x1 + x2,
#     mu2 = y2 ~ x1,
#     sigma1 = ~ x1 + x2,
#     sigma2 = ~ x1,
#     rho12 = ~ x1 + x2
#   ),
#   family = c(gaussian(), gaussian()),
#   data = dat
# )

## ----rejection-test-examples, eval = FALSE------------------------------------
# expect_snapshot(
#   {
#     drmTMB(
#       drm_formula(y ~ x + meta_V(V = V), nu ~ x),
#       family = student(),
#       data = dat
#     )
#   },
#   error = TRUE
# )
# 
# expect_snapshot(
#   {
#     drmTMB(
#       drm_formula(
#         mu1 = y1 ~ x + (1 | id),
#         mu2 = y2 ~ x,
#         sigma1 = ~ 1,
#         sigma2 = ~ 1,
#         rho12 = ~ 1
#       ),
#       family = c(gaussian(), gaussian()),
#       data = dat
#     )
#   },
#   error = TRUE
# )