GLMM Reference

library(lme4)

# Binary — random intercept
fit <- glmer(y ~ time + x + (1 | id), data = df, family = binomial)

# Binary — random intercept + slope
fit <- glmer(y ~ time + x + (time | id), data = df, family = binomial)

# Binary — uncorrelated random effects
fit <- glmer(y ~ time + x + (1 | id) + (0 + time | id),
             data = df, family = binomial)

# Poisson count
fit <- glmer(y ~ time + x + (1 | id), data = df, family = poisson)

# AGQ estimation (binary, single random effect only)
fit <- glmer(y ~ time + x + (1 | id), data = df,
             family = binomial, nAGQ = 7)

library(glmmTMB)

# Binary
fit <- glmmTMB(y ~ time + x + (1 | id), data = df, family = binomial)

# Poisson
fit <- glmmTMB(y ~ time + x + (1 | id), data = df, family = poisson)

# Negative binomial (quadratic parameterization)
fit <- glmmTMB(y ~ time + x + (1 | id), data = df, family = nbinom2)

# Negative binomial (linear parameterization)
fit <- glmmTMB(y ~ time + x + (1 | id), data = df, family = nbinom1)

# Zero-inflated negative binomial
fit <- glmmTMB(y ~ time + x + (1 | id),
               zi = ~ 1,
               data = df, family = nbinom2)

# Zero-inflated with predictors of zero-inflation
fit <- glmmTMB(y ~ time + x + (1 | id),
               zi = ~ time + x,
               data = df, family = nbinom2)

# lme4
fixef(fit)
summary(fit)$coefficients

# glmmTMB
fixef(fit)$cond          # Conditional model
fixef(fit)$zi            # Zero-inflation model (if present)

# lme4
exp(fixef(fit))
exp(confint(fit, method = "Wald", parm = "beta_"))

# glmmTMB
exp(fixef(fit)$cond)
exp(confint(fit, parm = "beta_"))

# Variance components
VarCorr(fit)

# Individual random effects (BLUPs / conditional modes)
ranef(fit)

# lme4: random effect SDs
as.data.frame(VarCorr(fit))

# glmmTMB: random effect SDs
sigma(fit)              # Dispersion
VarCorr(fit)$cond$id    # Random effect covariance matrix

# Conditional predictions (includes random effects)
predict(fit, type = "response")

# Population-level predictions (random effects = 0)
predict(fit, newdata = nd, type = "response", re.form = NA)

# Link-scale predictions
predict(fit, type = "link", re.form = NA)

# lme4 — Wald CIs
confint(fit, method = "Wald")

# lme4 — Profile CIs (slower, more accurate)
confint(fit, method = "profile")

# glmmTMB — Wald CIs
confint(fit)

# lme4
anova(fit_simple, fit_complex)

# glmmTMB
anova(fit_simple, fit_complex)

AIC(fit1, fit2)
BIC(fit1, fit2)

# Compare random intercept vs. random intercept + slope
fit1 <- glmer(y ~ time + (1 | id), data = df, family = binomial)
fit2 <- glmer(y ~ time + (time | id), data = df, family = binomial)
anova(fit1, fit2)

# Pearson residual ratio
overdisp <- sum(residuals(fit, type = "pearson")^2) / df.residual(fit)

# performance package
check_overdispersion(fit)

# Compare observed vs expected zeros
observed_zeros <- mean(df$y == 0)
predicted_zeros <- mean(predict(fit, type = "response") < 0.5)  # Rough

# Formal test (performance)
check_zeroinflation(fit)

# QQ-plot of random intercepts
re <- ranef(fit)$id[[1]]   # lme4
qqnorm(re); qqline(re)

# Shapiro-Wilk test (small samples only)
shapiro.test(re)

# Pearson residuals vs fitted
plot(fitted(fit), residuals(fit, type = "pearson"),
     xlab = "Fitted", ylab = "Pearson Residuals")
abline(h = 0, lty = 2)

# lme4: check for warnings
fit@optinfo$conv$lme4

# glmmTMB: check convergence
fit$sdr$pdHess  # Should be TRUE

# lme4: try different optimizers
fit <- glmer(y ~ time + (1 | id), data = df, family = binomial,
             control = glmerControl(optimizer = "bobyqa",
                                    optCtrl = list(maxfun = 100000)))

# glmmTMB: try different optimizer
fit <- glmmTMB(y ~ time + (1 | id), data = df, family = binomial,
               control = glmmTMBControl(optimizer = optim,
                                        optArgs = list(method = "BFGS")))

# Check which component is singular
VarCorr(fit)
# If slope variance ≈ 0, drop it:
fit_simple <- glmer(y ~ time + (1 | id), data = df, family = binomial)

p <- plogis(log_odds)       # = 1 / (1 + exp(-log_odds))

log_odds <- qlogis(p)       # = log(p / (1 - p))

count <- exp(log_count)

OR <- exp(beta)

IRR <- exp(beta)

pnorm(0, mean = fixef(fit)["time"],
      sd = sqrt(VarCorr(fit)$id["time", "time"]),
      lower.tail = FALSE)

fit <- glmer(y ~ time + treatment + sex + (1 | id),
             data = df, family = binomial)

fit <- glmer(y ~ time + stress_current + (1 | id),
             data = df, family = binomial)

# Time × treatment interaction
fit <- glmer(y ~ time * treatment + (1 | id),
             data = df, family = binomial)

# Count model with varying exposure
fit <- glmmTMB(events ~ time + treatment + (1 | id),
               offset = log(exposure_time),
               data = df, family = poisson)

library(ordinal)
fit <- clmm(ordered_y ~ time + treatment + (1 | id), data = df)

# Observations nested in people nested in sites
fit <- glmer(y ~ time + (1 | site/id), data = df, family = binomial)