## ----setup, message=FALSE----------------------------------------------------- library(sommer) ## ----basic-call, eval=FALSE--------------------------------------------------- # fit <- mmes(fixed = outcome ~ treatment, # random = ~ subject, # rcov = ~ units, # data = dat, # family = binomial(), # control = mmes_control(update = list(blendWeight = 0))) ## ----direct-pql-call, eval=FALSE---------------------------------------------- # fit_direct <- mmes(outcome ~ treatment, # random = ~ subject, data = dat, family = poisson(), # control = mmes_control(formulation = "direct", # working = list(maxit = 30, tol = 1e-6))) # summary(fit_direct) # fitted(fit_direct, type = "link") ## ----eval=FALSE--------------------------------------------------------------- # predict(fit_direct, D="subject", VarU=FALSE, sed=TRUE, pairwise=TRUE) # fit_direct <- postPEV(fit_direct, mode=1) # diagonal coefficient error covariance # # postPEV(fit_direct, mode=2) returns full joint coefficient error covariance ## ----laplace-call, eval=FALSE------------------------------------------------- # fit_laplace <- mmes( # outcome ~ treatment, # random = ~ subject, # data = dat, # family = binomial(), # method = mmes_method(objective = "ml", integration = "laplace") # ) # summary(fit_laplace) # fit_laplace$covParNative ## ----laplace-backend, eval=FALSE---------------------------------------------- # control <- mmes_control( # likelihood = list(backend = "auto"), # execution = list(denseMemoryMB = 1000)) # fit <- mmes(count ~ treatment, random = ~ plant, # data = dat, family = poisson(), control = control, # method = mmes_method(objective = "ml", integration = "laplace")) # fit$executionPlan$glmm$backend ## ----laplace-acceleration, eval=FALSE----------------------------------------- # control <- mmes_control( # optimizer = list(gradient = "analytic"), # mode = list(warmStart = "auto"), # likelihood = list(backend = "random", # trace = "selected", blocks = "on")) # fit$laplaceDiagnostics$nativeWork ## ----laplace-vcc, eval=FALSE-------------------------------------------------- # laplace_control <- mmes_control( # diagnostics = list(verbose = FALSE, dateWarning = FALSE)) # setup_control <- laplace_control # setup_control$output$setup <- TRUE # setup <- mmes(y ~ x, random = ~ group + block, data = dat, # family = poisson(), control = setup_control, # method = mmes_method(objective = "ml", integration = "laplace")) # setup$vcParams # fit <- mmes(y ~ x, random = ~ group + block, data = dat, # family = poisson(), method=mmes_method(objective="ml", integration="laplace"), control=laplace_control, # vcc = data.frame(parameter = c(1, 2), group = 1, scale = c(2, 1))) ## ----gaussian-residual-blocks, eval=FALSE------------------------------------- # fit <- mmes(y ~ time + x, random = ~ group, # rcov = ~ vsm(ar1m(time), ism(record)), # data = dat, family = gaussian(), method=mmes_method(objective="ml", integration="laplace"), control=laplace_control) # # # Covariance between Gaussian traits paired within a record: # mixed <- mmes(y ~ trait + trait:x - 1, random = ~ group, # rcov = ~ vsm(usm(trait), ism(record)), # data = dat_long, # family = familym("trait", height = gaussian(), weight = gaussian(), # count = poisson()), method=mmes_method(objective="ml", integration="laplace"), control=laplace_control) ## ----zero-inflated-poisson-call, eval=FALSE----------------------------------- # fit_zip <- mmes( # count ~ treatment, # random = ~ subject, # data = dat, # family = zipoisson(pi = 0.1), # method = mmes_method(objective = "ml", integration = "laplace"), # familyParameters = list(zeroInflation = 0.1, estimate = TRUE) # ) # fit_zip$familyParameters ## ----zero-formula-call, eval=FALSE-------------------------------------------- # fit_zip_formula <- mmes( # count ~ treatment, # random = ~ subject, # data = dat, # family = zipoisson(pi=0.1, zero=~treatment), # method = mmes_method(objective = "ml", integration = "laplace") # ) # fit_zip_formula$zeroCoefficients # summary(fit_zip_formula)$zeroCoefficients ## ----native-mixed-call, eval=FALSE-------------------------------------------- # fit_mixed <- mmes( # outcome ~ trait - 1 + trait:covariate, # random = ~ vsm(usm(trait), ism(subject)), # data = long_data, # family = familym( # "trait", binary = binomial("probit"), counts = nbinom2(theta=2), # positive = Gamma("log"), continuous = gaussian() # ), # method = mmes_method(objective="ml", integration="laplace"), # familyParameters=list( # counts=list(theta=2, estimate=TRUE), # positive=list(dispersion=0.3, estimate=TRUE), # continuous=list(dispersion=1, estimate=FALSE) # ) # ) # fit_mixed$familyParameters # fit_mixed$dispersion ## ----ordinal-interface, eval=FALSE-------------------------------------------- # ctrl <- mmes_control() # fit_ord <- mmes(rating ~ treatment + offset(exposure), # random = ~ group, data = dat, family = ordinal(), control = ctrl, # method = mmes_method(objective = "ml", integration = "laplace")) # fitted(fit_ord) # n by K labelled probability matrix # predict(fit_ord, newdata = newdat) # predict(fit_ord, re.form = NA) # zero random effects # predict(fit_ord, re.form = "population") # predict(fit_ord, type = "link") # latent predictor, not category codes # predict(fit_ord, type = "mean") # coding-dependent expected category # summary(fit_ord)$thresholds # cutpoint uncertainty, post-fit only # predict(fit_ord, se.fit = TRUE, cov.fit = TRUE) # simulate(fit_ord, nsim = 2, seed = 12) # ordered factors with original levels ## ----categorical-familym-interface, eval=FALSE-------------------------------- # fit_mix <- mmes(response ~ trait + trait:x - 1, # random = ~ vsm(ism(group), sigma2 = 0.2), # data = longdata, # family = familym("trait", rating = ordinal(), # choice = multinomial(reference = "3"), count = poisson()), # control = ctrl, method=mmes_method(objective="ml", integration="laplace")) # fitted(fit_mix) # named list: probability matrices/vectors # predict(fit_mix, newdata = newlongdata, type = "response") # predict(fit_mix, re.form = NA) # same list, zero random effects # predict(fit_mix, se.fit = TRUE, cov.fit = TRUE, uncertainty = "full") # summary(fit_mix)$thresholds # trait-labelled cutpoint SE, post-fit only # postPEV(fit_mix, mode = 1) # simulate(fit_mix, nsim = 2, seed = 12) # original numeric codes and row order # boot_mmes(fit_mix, nsim = 199, seed = 12) ## ----latent-process-recipes, eval=FALSE--------------------------------------- # laplace_control <- mmes_control( # diagnostics = list(verbose = FALSE, dateWarning = FALSE)) # # # Each plant has its own temporally correlated hidden process. # # Use a factor with explicitly ordered, equally spaced time levels. # dat$week <- factor(dat$week, levels = sort(unique(dat$week))) # temporal <- mmes(count ~ treatment, # random = ~ vsm(ar1m(week), ism(plant)), # data = dat, family = poisson(), method=mmes_method(objective="ml", integration="laplace"), control=laplace_control) # # # Independent fields, with a separable hidden spatial process in each field. # spatial <- mmes(diseased ~ treatment, # random = ~ vsm(ar1m(row), ar1m(column), ism(field)), # data = plots, family = binomial(), method=mmes_method(objective="ml", integration="laplace"), control=laplace_control) # # # In long data, this one process is shared across response types. # shared <- mmes(y ~ trait + trait:treatment - 1, # random = ~ vsm(ar1m(week), ism(plant)), # data = long, # family = familym("trait", count = poisson(), diseased = binomial()), # method=mmes_method(objective="ml", integration="laplace"), control=laplace_control) # # # These are different scientific targets: # predict(temporal, newdata = visits, type = "response") # fitted plant process # predict(temporal, newdata = visits, re.form = NA, type = "response") # process set to zero # predict(temporal, newdata = visits, re.form = "population", # type = "response", se.fit = TRUE, cov.fit = TRUE) # average over fitted process law # postPEV(temporal) # uncertainty about fitted latent coefficients, post-fit only ## ----structured-random, eval=FALSE-------------------------------------------- # fit <- mmes(y ~ environment, # random = ~ vsm(usm(environment), ism(genotype)), # rcov = ~ units, # data = dat, # family = poisson(), # control = mmes_control(update = list(blendWeight = 0))) ## ----offset, eval=FALSE------------------------------------------------------- # fit <- mmes(events ~ treatment + offset(log(exposure)), # random = ~ site, # rcov = ~ units, # data = dat, # family = poisson(), # control = mmes_control(update = list(blendWeight = 0))) ## ----family-data, eval=FALSE-------------------------------------------------- # set.seed(2026) # n_group <- 30 # n_per_group <- 8 # dat <- data.frame( # group = factor(rep(seq_len(n_group), each = n_per_group)), # x = rep(c(0, 1), length.out = n_group * n_per_group), # exposure = runif(n_group * n_per_group, 0.5, 2) # ) ## ----gaussian-example, eval=FALSE--------------------------------------------- # dat$y_gaussian <- 2 + 0.7 * dat$x + rnorm(nrow(dat), sd = 1) # fit_gaussian <- mmes( # y_gaussian ~ x, random = ~ group, rcov = ~ units, data = dat, # family = gaussian() # ) ## ----binomial-example, eval=FALSE--------------------------------------------- # probability <- plogis(-0.7 + 1.1 * dat$x) # dat$y_binomial <- rbinom(nrow(dat), size = 1, prob = probability) # fit_binomial <- mmes(y_binomial ~ x, # random = ~ group, # rcov = ~ units, # data = dat, # family = binomial(), # control = mmes_control(update = list(blendWeight = 0))) ## ----poisson-example, eval=FALSE---------------------------------------------- # rate <- dat$exposure * exp(0.2 + 0.4 * dat$x) # dat$y_poisson <- rpois(nrow(dat), lambda = rate) # fit_poisson <- mmes(y_poisson ~ x + offset(log(exposure)), # random = ~ group, # rcov = ~ units, # data = dat, # family = poisson(), # control = mmes_control(update = list(blendWeight = 0))) ## ----gamma-example, eval=FALSE------------------------------------------------ # mean_gamma <- exp(0.3 + 0.25 * dat$x) # dat$y_gamma <- rgamma(nrow(dat), shape = 4, scale = mean_gamma / 4) # fit_gamma <- mmes(y_gamma ~ x, # random = ~ group, # rcov = ~ units, # data = dat, # family = Gamma(link = "log"), # control = mmes_control(update = list(blendWeight = 0))) ## ----inverse-gaussian-example, eval=FALSE------------------------------------- # dat$y_inverse_gaussian <- exp(0.2 + 0.3 * dat$x + rnorm(nrow(dat), sd = 0.2)) # fit_inverse_gaussian <- mmes(y_inverse_gaussian ~ x, # random = ~ group, # rcov = ~ units, # data = dat, # family = inverse.gaussian(link = "log"), # control = mmes_control(update = list(blendWeight = 0))) ## ----quasi-example, eval=FALSE------------------------------------------------ # dat$y_quasi <- 1 + 0.5 * dat$x + rnorm(nrow(dat), sd = 0.5) # fit_quasi <- mmes(y_quasi ~ x, # random = ~ group, # rcov = ~ units, # data = dat, # family = quasi(link = "identity", variance = "constant"), # control = mmes_control(update = list(blendWeight = 0))) # # fit_quasibinomial <- mmes(y_binomial ~ x, # random = ~ group, # rcov = ~ units, # data = dat, # family = quasibinomial(), # control = mmes_control(update = list(blendWeight = 0))) # # fit_quasipoisson <- mmes(y_poisson ~ x + offset(log(exposure)), # random = ~ group, # rcov = ~ units, # data = dat, # family = quasipoisson(), # control = mmes_control(update = list(blendWeight = 0))) ## ----extract, eval=FALSE------------------------------------------------------ # fitted(fit_poisson) # fitted(fit_poisson, type = "link") # residuals(fit_poisson, type = "deviance") # residuals(fit_poisson, type = "working") # # fit_poisson$pqlMonitor # fit_poisson$pqlConverged # fit_poisson$family ## ----pql-control, eval=FALSE-------------------------------------------------- # fit <- mmes(y_poisson ~ x + offset(log(exposure)), # random = ~ group, # rcov = ~ units, # data = dat, # family = poisson(), # control = mmes_control( # update = list(blendWeight = 0), # working = list(maxit = 30, tol = 1e-6), # convergence = list(maxIters = 20), # solve = list(method = "auto") # ))