Package {adlaplace}

Stream per-unit observation gradients (grad-only tapes)

Description

For each observation unit (ELGM stratum column, or row of y when there is no ELGM map), records a temporary one-unit observation tape with hessian_sparsity = FALSE, evaluates the sparse gradient at x, then destroys the tape. Does not run inner_opt.

Usage

.grad_obs_units_cpp(model, x, config, units, batch_size, inner, negative)

Arguments

Value

A dgCMatrix with one column per unit.

Valid model type factor levels

Description

Character vector of allowed levels for the model_role slot on model_term objects.

Usage

.model_role_levels

Value

A character vector of length 3, with elements "fixed", "random", and "response" in that order. These are the allowed levels of the model_role factor on model_term objects, indicating a fixed effect, a random effect, or the response.

Build AD function with Hessian templates

Description

Build AD function with Hessian templates

Usage

ad_pack(
  x,
  config = NULL,
  num_threads = 1L,
  reorder_shards = c("none", "gradient", "hessian", "third"),
  templates = NULL,
  ...
)

## S4 method for signature 'ad_pack_ptr'
ad_pack(
  x,
  config = NULL,
  num_threads = 1L,
  reorder_shards = c("none", "gradient", "hessian", "third"),
  templates = NULL,
  ...
)

## S4 method for signature 'density_data'
ad_pack(
  x,
  config = NULL,
  num_threads = 1L,
  reorder_shards = c("none", "gradient", "hessian", "third"),
  templates = NULL,
  ...
)

## S4 method for signature 'list'
ad_pack(
  x,
  config = NULL,
  num_threads = 1L,
  reorder_shards = c("none", "gradient", "hessian", "third"),
  templates = NULL,
  ...
)

Arguments

Value

Object of class ad_pack.

AD function with Hessian templates attached

Description

AD function with Hessian templates attached

Slots

ptr

Raw combined handle (ad_pack_ptr).

group_sparsity

Per-group inner gradient sparsity patterns.

outer

Outer Hessian template (dgCMatrix).

inner

Inner-gamma Hessian template (dgCMatrix).

map_outer

Outer Hessian shard map.

map_inner

Inner Hessian shard map.

parallel_map

Thread-affinity map (ngCMatrix); rows are shards, columns are threads, with one nonzero per shard row.

chol_inner

Symbolic LDL factor or empty sparse matrix.

chol_inner_list

Numeric LDL list for C++ (L1, Linv, perm, perm_inv, half_H_inv, H_inv), plus trace_columns (per-shard column indices for trace_hinv_t) and hinv_q_index (per-shard slots of Q into H_inv, when a shard defines hinv_trace_index).

sizes

Named numeric vector beta/gamma/theta.

info

List of parameter metadata (beta, gamma, theta, parameters); populated from model_data()$term_data$info when the handle is built from a model-data bundle, otherwise empty.

Retape an AD pack after dropping observation units

Description

Builds a new ad_pack on fit$model_data with observation shards filtered to the remaining units (same coarse sharding as fit$config$obs_groups when present). Random and parameter shards are retaped; reuse of the merged fit pack is not attempted.

Usage

ad_pack_drop(fit, drop, num_threads = NULL, verbose = FALSE)

Arguments

Value

An ad_pack usable with inner_opt at the original outer MLE.

See Also

obs_groups_units, influence_scores, amip_set

Build raw AD handle for one density shard

Description

Constructs CppAD tapes for a single shard. The density kind and name come from data@ad_kind and data@density. For ad_kind = "random", data@precision is passed straight to the backend. For random_diagonal it must be a numeric vector of diagonal precision weights with length == ncol(gamma_map). For random_mult it must be list(Q, log_det, rank) where Q is a square dsCMatrix with uplo = "U" (upper triangle only); use adlaplace:::as_dsC_upper(Q) to coerce a dgCMatrix. Missing precision means the shard is not built (e.g. diffuse rpoly with sd = Inf via model_data()); calling ad_pack_ptr() without it is an error.

Usage

ad_pack_ptr(data, config = NULL, info = NULL)

Arguments

Details

Merge handles for multiple shards with c() before calling ad_pack.

Value

External pointer of class ad_pack_ptr.

Raw AD handle (external pointer)

Description

S4 class for a C++ AD backend handle. Objects are created by ad_pack_ptr, combined with c.ad_pack_ptr, and typically passed to ad_pack to attach Hessian templates.

See Also

ad_pack_ptr, ad_pack, clone_ad_pack_ptr

Rebuild an ad_pack for a new model-data design

Description

Covariate values in observation XTp are baked into CppAD tapes at build time, so they cannot be mutated in place. This helper retapes a new model_data bundle and reuses Hessian templates from ad_pack whenever sizes and shard counts match (identical sparsity across exposure samples). Prefer building a list of packs once for modular inference rather than calling this inside every likelihood evaluation.

Usage

ad_pack_set_X(
  ad_pack,
  model_data,
  config = list(),
  num_threads = 1L,
  reorder_shards = c("none", "gradient", "hessian", "third")
)

Arguments

Value

A new ad_pack.

See Also

ad_pack, replace_vars

Fit a hierarchical model by Laplace-approximate maximum likelihood

Description

One-call interface to the adlaplace pipeline: parse the formula with model_data(), partition observations into AD shards with obs_groups(), build the AD handle with ad_pack(), maximize the Laplace-approximate marginal likelihood with optim() using the exact profile gradient (outer_fn / outer_gr), and return a fitted-model object with standard accessor methods.

Usage

adlaplace(
  formula,
  data,
  config = list(),
  control = list(),
  control_inner = list(maxit = 100L, report.level = 0, report.freq = 0),
  method = "L-BFGS-B",
  hessian = NULL,
  start = NULL,
  verbose = FALSE,
  na_omit = TRUE
)

Arguments

Value

An object of class "adlaplace_fit": a list with components

call, formula

The matched call and model formula.

par_info

Parameter table with mle, se, mle_internal, se_internal, bounds, and log.

gamma

Named random-effect modes at the optimum.

nobs

Number of observations.

details

Laplace evaluation at the optimum, including outer_opt, inner_opt, vcov, log_lik, gradient, and hessian.

model_data, ad_pack, config, cache

Objects needed by the accessor methods and for warm restarts.

If optim throws (for example L-BFGS-B when fn is not finite), the same class is returned with par_info estimates set to NA, details$outer_opt = NULL, and details$error set to the condition message. ad_pack and cache are still present so the Laplace objective can be inspected or re-run without rebuilding tapes:

stats::optim(
  fit$cache$last_par_fn,
  fn = adlaplace::outer_fn, gr = adlaplace::outer_gr,
  method = fit$method, control = fit$control,
  lower = fit$par_info$lower, upper = fit$par_info$upper,
  config = fit$config, ad_pack = fit$ad_pack,
  cache = fit$cache, control_inner = fit$control_inner
)

See Also

summary.adlaplace_fit, predict.adlaplace_fit, log_lik_laplace, model_data, fit_outer

Examples

set.seed(1)
n <- 80
g <- factor(sample(4, n, replace = TRUE))
x <- rbinom(n, 1, 0.5)
eta <- 0.5 + 0.4 * x + rnorm(4, sd = 0.3)[g]
nb_sd <- 0.3
z <- rgamma(n, nb_sd^(-2), nb_sd^(-2))
dat <- data.frame(y = rpois(n, exp(eta) * z), x = x, g = g)
md <- model_data(
  nbinom(y, lower = 1e-9, init = 0.15) ~ x + iid(g, init = 0.3),
  data = dat
)
config <- list(
  transform_theta = TRUE,
  obs_groups = obs_groups(md$term_data$A, num_shards = 5),
  verbose = FALSE
)
ap <- ad_pack(md, config)
cache <- new.env(parent = emptyenv())
x0 <- md$term_data$info$parameters$init
outer_fn(x = x0, config = config, cache = cache, ad_pack = ap)
outer_gr(x = x0, config = config, cache = cache, ad_pack = ap)

fit <- adlaplace(md, config = list(num_shards = 5L, verbose = FALSE))
summary(fit)
confint(fit)

Attach hessian_map() result to an ad_pack handle

Description

Copies outer/inner templates and maps into the C++ ad_pack handle.

Usage

adlaplace_attach_hessian(handle, hessian_pack)

Arguments

C++ backend entry points

Description

Low-level entry points exposed to R. Accept an ad_pack S4 object (via @ptr) or a raw ad_pack_ptr.

Usage

joint_log_dens(ad_pack, x, ad_shards = NULL, negative = TRUE)

grad(ad_pack, x, ad_shards = NULL, inner = FALSE, negative = TRUE)

hessian(
  ad_pack,
  x,
  ad_shards = NULL,
  inner = FALSE,
  verbose = FALSE,
  negative = TRUE
)

trace_hinv_t(ad_pack, x, LinvPt, LinvPtColumns, verbose = FALSE, H_inv = NULL)

fun_obj_fdfh(ad_pack, parameters, gamma, inner = TRUE, verbose = FALSE)

Arguments

Value

See individual functions.

Sign convention

With default negative = TRUE, joint_log_dens(), grad(), and hessian() match inner_opt() (negative log-density). Set negative = FALSE for the joint log density and its derivatives at the same x.

Thread affinity

These are **serial debug** APIs. They error if the handle was built with ad_pack(..., num_threads > 1). Use a plain ad_pack_ptr (or ad_pack(..., num_threads = 1)), or clone_ad_pack_ptr() before assigning multi-thread affinity. Parallel evaluation uses fun_obj_fdfh / inner_opt / log_lik_laplace on the affined handle.

Methods for fitted adlaplace models

Description

Standard accessor methods for objects of class "adlaplace_fit" returned by adlaplace(). Estimates are stored on the optimization scale (log scale for parameters with par_info$log = TRUE, e.g. standard deviations); methods with a transform argument report the natural scale by exponentiating and, for standard errors, applying the delta method.

Usage

## S3 method for class 'adlaplace_fit'
coef(object, transform = TRUE, ...)

## S3 method for class 'adlaplace_fit'
vcov(object, transform = FALSE, ...)

## S3 method for class 'adlaplace_fit'
logLik(object, ...)

## S3 method for class 'adlaplace_fit'
nobs(object, ...)

## S3 method for class 'adlaplace_fit'
confint(object, parm, level = 0.95, transform = TRUE, ...)

## S3 method for class 'adlaplace_fit'
fitted(object, ...)

## S3 method for class 'adlaplace_fit'
predict(object, newdata, n = 500L, ...)

## S3 method for class 'adlaplace_fit'
print(x, digits = max(3L, getOption("digits") - 3L), ...)

## S3 method for class 'adlaplace_fit'
summary(object, level = 0.95, ...)

## S3 method for class 'summary.adlaplace_fit'
print(x, digits = max(3L, getOption("digits") - 3L), ...)

## S3 method for class 'adlaplace_fit'
plot(x, k = 200L, draws = 100L, ...)

Arguments

Details

predict draws random effects from the Laplace (Gaussian) approximation at the mode via sim_fit and returns the linear predictor on newdata: a nrow(newdata) by n matrix of simulated values (plus fixed-effect contributions from terms whose variables appear in newdata). Without newdata it returns fitted(object).

plot draws each smooth term (terms with more than one knot) on a grid spanning its knots, showing simulated curves from the Laplace approximation with the pointwise mean in red. With no smooth terms it plots the random-effect modes.

Value

coef

Named numeric vector of outer parameter estimates. With transform = TRUE (the default), parameters stored on the log scale (for example standard deviations) are exponentiated to the natural scale.

vcov

Variance-covariance matrix of the outer parameters on the optimization scale. With transform = TRUE, rows and columns for log-scale parameters are mapped to the natural scale by the delta method.

logLik

An object of class "logLik": the Laplace-approximate marginal log-likelihood, with attributes df (number of outer parameters) and nobs.

nobs

Integer number of observations.

confint

Numeric matrix of Wald intervals, one row per outer parameter. With transform = TRUE, endpoints for log-scale parameters are exponentiated to the natural scale.

fitted

Numeric vector of length nobs: the linear predictor at the joint mode of the fixed effects and random effects.

predict

Without newdata, the same vector as fitted. With newdata, a matrix with nrow(newdata) rows and n columns of linear predictors drawn from the Laplace approximation to the random effects.

summary

A list of class "summary.adlaplace_fit" with the coefficient table (natural-scale estimates, standard errors, and interval endpoints), the log-likelihood, degrees of freedom, the number of observations and random effects, the outer optimizer convergence code, and the maximum absolute profile gradient.

print, print.summary, plot

No return value used by the caller. They print or draw and return the input object invisibly. plot draws each smooth term as curves simulated from the Laplace approximation, with the pointwise mean highlighted, or the random-effect modes when the model has no smooth term.

Build aghq objective list from an adlaplace fit

Description

Returns list(fn, gr, he) suitable for aghq, closing over the fitted AD handle and cache. The wrappers return the negated Laplace-approximate log posterior (same convention as outer_fn / outer_gr), so call aghq::default_control(negate = TRUE).

Usage

aghq_ff(fit, control_inner = NULL)

Arguments

Value

A list with fn, gr, and he. The Hessian is a numerical Jacobian of gr (no AD outer Hessian is available).

See Also

outer_fn, outer_gr, adlaplace

Approximate most influential removal set (AMIP)

Description

Greedy rule: drop units with the largest scores of the same sign as phi until \sum_{u \in S} \psi_u matches the sign of \phi and |\sum| \ge |\phi|. Leaving out S changes \phi by approximately -\sum \psi_S.

Usage

amip_set(psi, phi)

Arguments

Value

A list with idx (1-based indices to drop), alpha (fraction of units), and approx_change. If the linear approximation cannot flip the sign, idx is empty and alpha is NA.

See Also

influence_scores, ad_pack_drop

Apply log transform to selected theta columns

Description

When active is TRUE, logs values in cols for rows where theta_info$log is TRUE. Missing or NA log entries default to TRUE.

Usage

apply_theta_log(theta_info, cols = "init", active = TRUE)

Arguments

Value

A copy of theta_info with selected entries logged.

Coerce a sparse matrix to dgCMatrix.

Description

Coerce a sparse matrix to dgCMatrix.

Usage

as_dgC(m)

Coerce a sparse matrix to upper-triangle dsCMatrix (uplo = "U").

Description

Used for random_mult precision Q. The C++ backend rejects dgCMatrix; coerce here before building the tape.

Usage

as_dsC_upper(m)

Coerce a sparse matrix to ngCMatrix (Matrix >= 1.6).

Description

Coerce a sparse matrix to ngCMatrix (Matrix >= 1.6).

Usage

as_ngC(m)

Assign OpenMP owner threads to all shards

Description

By default sets owner_thread = shard_index % num_threads. When owners is supplied (length n_shards, 0-based thread ids in [0, num_threads)), those values are used instead (e.g. LPT balancing from ad_pack(..., reorder_shards=)). Called from ad_pack() after templates are attached. Required before parallel inner_opt / trace_hinv_t.

Usage

assign_owner_threads(handle, num_threads, owners = NULL)

Arguments

Binomial observation term

Description

Model term for the observation-level binomial log density with a logit link, registered as binomial_obs. For each observation with linear predictor eta = X beta + A gamma, trial count N, and success count y, the log density contribution is y * eta - N * log(1 + exp(eta)) (up to a data-only constant). When size is omitted, N = 1 (Bernoulli). There are no observation-level hyperparameters.

Creates a response term wired to the binomial_obs AD shard. A Bernoulli likelihood (size omitted) has no observation-level hyperparameters. With size, trial counts are read from that column of data and stored as observation weights.

Usage

binom(x, size = NULL)

Arguments

Details

Binomial observation model term

Value

A binomial object.

Slots (inherited from model)

ad_pack

Character scalar "binomial_obs".

ad_kind

Character scalar "observations".

Slots

size

Optional column name for per-observation trial counts.

By-Group Classes and Functions

Description

Functions for creating and managing by_group objects for hierarchical model grouping information.

Usage

by_group(term, data = NULL, levels = NULL)

Arguments

Value

A by_group object

Examples

# Create with explicit levels
bg <- by_group(term = "group", levels = c("A", "B", "C"))

# Create from data
dat <- data.frame(group = c("A", "B", "A", "C"))
bg <- by_group(term = "group", data = dat)

By-Group Class

Description

A class for grouping information in hierarchical model terms. Contains the grouping variable name and its levels/labels.

Slots

term

Character vector of length 1, the grouping variable name

levels

Character vector of unique levels for the grouping variable

labels

Character vector of formatted labels for the levels

Combine ad_pack_ptr handles

Description

Moves AD shards from one or more ad_pack_ptr objects into a single handle. Inputs are cleared after the merge (ownership moves).

Usage

## S3 method for class 'ad_pack_ptr'
c(x, ...)

Arguments

Value

A combined ad_pack_ptr.

See Also

clone_ad_pack_ptr, ad_pack

Merge partial AD handles into one raw handle

Description

Concatenates shards from ad_pack_ptr (or other partial builders) in list order. Does not attach hessian_map; use ad_pack when templates are needed.

Usage

c_ad_pack_ptr(handles)

Arguments

Value

Combined external pointer of class ad_pack_ptr.

See Also

ad_pack_ptr, ad_pack

Column-bind design blocks; coerce to sparse when any block is sparse.

Description

Drops NULL pieces. When any remaining piece inherits sparseMatrix, every piece is coerced to dgCMatrix before cbind so a leading dense block cannot strip the Matrix class. When none are sparse, returns a base dense matrix via cbind.

Usage

cbind_design_blocks(pieces)

Arguments

Value

A design matrix (base matrix or dgCMatrix).

Deep copy of an ad_pack_ptr handle

Description

Returns a new ad_pack_ptr with independent C++ state (CppAD tapes and sparsity patterns). Hessian templates are not copied; call ad_pack on the clone to attach maps. The source handle is left unchanged, unlike c.ad_pack_ptr which moves shards and clears the inputs. OpenMP thread affinity from a prior ad_pack(..., num_threads > 1) is cleared on the clone so joint_log_dens / grad / hessian remain available as serial debug APIs.

Usage

clone_ad_pack_ptr(x)

Arguments

Value

A new ad_pack_ptr.

Deep copy of an ad_pack_ptr handle

Description

Clones CppAD tapes and sparsity patterns into a new external pointer. The source handle is unchanged (unlike c() on ad_pack_ptr, which moves shards and clears sources).

Usage

clone_ad_pack_ptr_(handle)

Arguments

Value

New ad_pack_ptr with independent C++ state.

Coerce map shorthand: length-1 nrow, length-2 c(row, nrow), or list(indices, nrow).

Description

Coerce map shorthand: length-1 nrow, length-2 c(row, nrow), or list(indices, nrow).

Usage

coerce_map_shorthand(x, what)

Parse Model Terms from Formula

Description

Parses a formula and creates model terms by evaluating constructor calls in the formula environment (so pkg::term(...) and terms from attached packages on the search path work). Falls back to the adlaplace namespace for built-in constructors.

Usage

collect_terms(formula, verbose = FALSE)

Arguments

Value

List of model term objects

Build raw AD handle for a random_diagonal shard

Description

Build raw AD handle for a random_diagonal shard

Usage

create_ad_shard_random_diagonal(model, config)

Arguments

Value

External pointer of class ad_pack_ptr.

Build raw AD handle for a random_mult shard

Description

Build raw AD handle for a random_mult shard

Usage

create_ad_shard_random_mult(model, config)

Arguments

Value

External pointer of class ad_pack_ptr.

Default observation shard map (one column, all observations)

Description

Default observation shard map (one column, all observations)

Usage

default_obs_groups(n_obs)

Arguments

Value

A ngCMatrix with one shard column and a structural nonzero mapping each observation row to that shard (0-based row indices).

Per-shard layout for AD density evaluation

Description

Per-shard layout for AD density evaluation

Usage

density_data(
  y = missing(),
  A = NULL,
  X = NULL,
  beta_map = NULL,
  gamma_map = NULL,
  theta_map = NULL,
  elgm_matrix = NULL,
  density = NA_character_,
  ad_kind = NA_character_,
  package = NA_character_,
  precision = NULL,
  weights = numeric(0)
)

Arguments

Value

An S4 object of class "density_data": one automatic-differentiation shard. Slots hold the response, the transposed fixed- and random-effect design blocks, the maps from local design columns into the global parameter vectors, the registered density name, and the shard kind ("observations", "parameters", or "random").

Functions

• density_data(): Construct a density_data object

Slots

y

Response vector.

ATp

Transpose of random-effects design matrix.

XTp

Transpose of fixed-effects design matrix.

beta_map

Beta parameter map (ngCMatrix). nrow is the global beta block size; column j has one structural nonzero at row r meaning local design column j of X uses global beta r. When ncol(X) > 0, require ncol(beta_map) == ncol(X).

gamma_map

Gamma parameter map (ngCMatrix); at most one structural nonzero per row and column. Same column convention as beta_map for local A columns into the global gamma block.

theta_map

Theta parameter map (ngCMatrix). Column j selects global theta row r for this shard's theta component j.

elgm_matrix

Optional exposure-lag map (ngCMatrix; empty by default).

density

Registered AD density name for this shard.

ad_kind

Shard kind ("observations", "parameters", "random").

package

Package name whose shared library records tapes for this shard (defaults to "adlaplace" when missing).

precision

Optional precision payload (any R object).

weights

Optional per-observation weights (e.g. binomial trial counts). Empty means all ones.

Build density_data from designs and config layout sizes.

Description

Used by examples with a plain data list (y, A, X).

Usage

density_data_from_config_matrices(y, A, X, config, theta_local_row = 0L)

Scan joint log-density AD derivatives vs finite differences

Description

Evaluates the joint log density and its AD gradient/Hessian while varying one coordinate of the full parameter vector (beta, gamma, theta). Returns plain matrices shaped for matplot()/matlines().

Usage

deriv_scan_joint(
  ad_pack,
  x,
  index,
  width = 0.1,
  n = 11L,
  ad_shards = NULL,
  negative = FALSE
)

Arguments

Value

A list with:

Scan Laplace profile AD derivatives vs finite differences

Description

Evaluates log_lik_laplace(deriv = TRUE) while varying one coordinate of the outer parameter vector (beta, theta only; gammas are profiled out). Returns plain matrices shaped for matplot()/matlines().

Usage

deriv_scan_laplace(
  ad_pack,
  x,
  index,
  width = 0.5,
  n = 7L,
  config = list(),
  gamma = NULL,
  control = list(maxit = 100L, report.level = 0, report.freq = 0)
)

Arguments

Value

A list with:

Transpose of a design matrix for ATp / XTp slots.

Description

Transpose of a design matrix for ATp / XTp slots.

Usage

design_Tp(M, n_obs, what = "M")

Dirichlet-multinomial case-crossover response term

Description

Dirichlet-multinomial case-crossover response term

Usage

dirichlet_multinom(
  x,
  by,
  init = 0.001,
  lower = 0,
  upper = Inf,
  parscale = 1,
  log = TRUE
)

Arguments

Value

A dirichlet_multinom object.

Dirichlet-multinomial response term

Description

Model term for the case-crossover Dirichlet-multinomial log density registered as dirichlet_multinomial.

Usage

## S4 method for signature 'dirichlet_multinom'
elgm_matrix(term, data)

## S4 method for signature 'dirichlet_multinom'
design(term, data)

## S4 method for signature 'dirichlet_multinom'
precision(term, data)

## S4 method for signature 'dirichlet_multinom'
theta_info(term)

## S4 method for signature 'dirichlet_multinom'
beta_info(term, data)

## S4 method for signature 'dirichlet_multinom'
random_info(term, data)

Arguments

Details

Dirichlet-multinomial case-crossover response term

Value

Sparse ngCMatrix mapping observations to strata columns.

Functions

• elgm_matrix(dirichlet_multinom): Case-crossover stratum map for AD shards.

• design(dirichlet_multinom): Design method (not used).

• precision(dirichlet_multinom): Precision method (not used).

• theta_info(dirichlet_multinom): Theta info for overdispersion sd.

• beta_info(dirichlet_multinom): Beta info method (not used).

• random_info(dirichlet_multinom): Random info method (not used).

Slots

by

Character vector of stratification variables for the case-crossover design.

Slots (inherited from model)

ad_pack

Character scalar "dirichlet_multinomial".

ad_kind

Character scalar "observations".

package

Character scalar "adlaplace".

Build config$obs_groups from design / ELGM when missing

Description

Used by adlaplace and ad_pack so that num_shards alone is enough for case-crossover (ELGM) models. When elgm_matrix has columns, shards index strata, not raw rows of y.

Usage

ensure_config_obs_groups(
  config,
  A,
  elgm_matrix = NULL,
  num_shards = NULL,
  num_threads = 1L,
  n_domain = NULL
)

Arguments

Details

When config$obs_groups is missing and config$obs_units is set (0-based unit indices), builds an identity shard map over those units via obs_groups_units. Optional config$obs_units_grouping selects "identity" (default) or "together".

Value

config with obs_groups filled when possible.

Midpoint finite differences along a coordinate grid

Description

Midpoint finite differences along a coordinate grid

Usage

fd_midpoint(y, grid)

Arguments

Value

List with mid (midpoints) and fd (slopes).

Fill missing config$beta / config$theta / config$gamma from term_data$info.

Description

Explicit config$beta, config$theta, and config$gamma are left unchanged. transform_theta is resolved to a per-theta logical vector (from info$theta$log unless the caller passed a full vector or FALSE). When config$theta is missing, seeds are apply_theta_log(info$theta, cols = "init")$init. When config$gamma is missing, seeds are zeros of length nrow(info$gamma) (info$gamma has no init column).

Usage

fill_config_from_info(config, info)

Outer optimization scaffolding for Laplace profile likelihoods

Description

Runs optim on a user-supplied outer objective / gradient (typically outer_fn / outer_gr, or a modular average of several packs) and returns the standard par_info / vcov / opt structure used by adlaplace.

Usage

fit_outer(
  fn,
  gr,
  par_meta,
  method = "L-BFGS-B",
  control = list(),
  hessian = FALSE,
  start = NULL,
  ...
)

Arguments

Value

A list with components opt, par_info, vcov, and error (character message when optim fails, otherwise NULL).

See Also

adlaplace, outer_fn

Format Laplace estimates using parameter metadata

Description

Maps flat Laplace parameter vectors onto labeled tables from info (as stored on ad_pack@info or model_data()$term_data$info).

Usage

format_parameters(
  info,
  full_parameters = NULL,
  parameters = NULL,
  gamma = NULL
)

Arguments

Value

A list with data frames parameters and gamma, or list() when info is missing or incomplete.

Fixed Polynomial Model Term

Description

Creates and manages fixed polynomial model terms.

Usage

fpoly(
  x,
  p = 2,
  ref_value = 0,
  init = .my_beta_init,
  lower = .my_beta_lower,
  upper = .my_beta_upper,
  parscale = .my_beta_parscale,
  basis_scale = 1
)

## S4 method for signature 'fpoly'
design(term, data)

## S4 method for signature 'fpoly'
beta_info(term, data)

Arguments

Value

A 'fpoly' term object.

A design matrix for the fixed polynomial term, or NULL if p.order is 0.

A data frame containing beta parameter information for the fixed polynomial term.

Functions

• design(fpoly): Design method for fpoly objects

• beta_info(fpoly): Beta info method for fpoly objects

Methods

The following methods are available for 'fpoly' objects:

design(term, data)

Creates design matrix for fpoly term

precision(term, data)

Creates precision matrix for fpoly term

theta_info(term)

Extracts theta parameter information

beta_info(term, data)

Extracts beta parameter information

random_info(term, data)

Extracts random effects information

Examples

# Create a fixed polynomial term
fpoly_term <- fpoly(x = "temp", p = 3, ref_value = 15)

Example GAMM simulation data

Description

Simulation data for the GAMM vignette example.

Format

A data frame dat with 500 rows and columns from mgcv::gamSim(6, n = 500), plus simulated negative-binomial responses used in the GAMM vignette (y, f2, mu, Z, etc.).

Source

Generated with mgcv::gamSim() and a custom negative-binomial response; see the regeneration comments in vignettes/gamm.Rmd.

Examples

data("gamm", package = "adlaplace")
exists("dat")
nrow(dat)

Gaussian observation term

Description

Model term for the observation-level Gaussian log density registered as gaussian_obs. This is the default observation term: a bare response symbol on the left-hand side of a formula (e.g. y ~ x) is coerced to normal(y) by collect_terms.

Creates a response term wired to the gaussian_obs AD shard, with a single residual standard deviation parameter (log scale during optimization).

Usage

normal(
  x,
  init = 1,
  lower = .my_theta_lower,
  upper = .my_theta_upper,
  parscale = .my_theta_parscale,
  log = TRUE
)

## S4 method for signature 'gaussian'
theta_info(term)

Arguments

Details

Gaussian observation model term

Value

A gaussian object.

Functions

• theta_info(gaussian): Theta info for the residual standard deviation.

Slots (inherited from model)

ad_pack

Character scalar "gaussian_obs".

ad_kind

Character scalar "observations".

Germany oral cavity cancer (Besag-York-Mollie example)

Description

Disease-mapping data distributed with INLA (Besag, York, and Mollie 1991). Counts y_i in 544 districts follow y_i \sim \mathrm{Poisson}(E_i e^{\eta_i}) with a BYM spatial random effect. The bundled ICAR precision Q_scaled is the INLA scale.model version with sum-to-zero constraint; prec is the random_mult precision payload (Q, log_det, rank).

Usage

germany

Format

A list with:

Y

Observed case counts (length 544).

E

Expected counts (length 544).

region

District index (length 544).

adj

District adjacency matrix (dgCMatrix, 544 x 544).

Q_scaled

Scaled ICAR precision (dsCMatrix, uplo = "U").

prec

List passed to density = "random_mult": Q (dsCMatrix, uplo = "U"), log_det, and rank.

Details

District boundaries are not stored in this object (to avoid a hard dependency on terra). They ship as system.file("extdata", "germany_map.gpkg", package = "adlaplace") and can be loaded with terra::vect(...) when terra is installed.

Source

INLA Germany data and germany.graph.

References

Besag, J., York, J., and Mollie, A. (1991). Bayesian image restoration, with two applications in spatial statistics. Annals of the Institute of Statistical Mathematics, 43(1), 1–20.

Rue, H., Martino, S., and Chopin, N. (2009). Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations. Journal of the Royal Statistical Society: Series B, 71(2), 319–392.

Build raw AD handle for observation shards only

Description

Build raw AD handle for observation shards only

Usage

get_ad_pack_raw_obs(model, config, name)

Arguments

Value

External pointer of class ad_pack_ptr.

Build raw AD handle for a parameters shard

Description

Build raw AD handle for a parameters shard

Usage

get_ad_pack_raw_parameters(model, config, name)

Arguments

Value

External pointer of class ad_pack_ptr.

Configured num_threads from ad_pack() thread assignment

Description

Configured num_threads from ad_pack() thread assignment

Usage

get_configured_num_threads(handle)

Arguments

Value

Integer num_threads if ad_pack() assigned threads; otherwise NA.

Whether a shard has an assigned OpenMP owner thread

Description

Whether a shard has an assigned OpenMP owner thread

Usage

get_owner_thread_assigned(handle, group)

Arguments

Value

Logical scalar.

Sparse structure sizes for one AD shard

Description

Layout and sparsity sizes for one AD shard

Usage

get_sizes(handle, group)

Arguments

Value

List with n_inner, n_outer, n_beta, n_theta, and nnz_grad_*, nnz_hes_*.

Sparse index patterns for one AD shard

Description

Sparse index patterns for one AD shard

Usage

get_sparse_pattern(handle, group)

Arguments

Value

List with grad, grad_inner, row_hess, col_hess, etc.

CppAD tape sizes for one AD shard

Description

CppAD tape sizes for one AD shard

Usage

get_tape_sizes(handle, group)

Arguments

Value

List with domain, n_global, size_op, size_var, and pattern nnz counts. domain / size_op / size_var are zero when the shard never recorded a tape (e.g. analytic FEM ssq with fem_analytic = TRUE); n_global and the nnz fields remain valid.

Build-time owner thread for one AD shard

Description

Build-time owner thread for one AD shard

Usage

get_thread_owner(handle, group)

Arguments

Value

Integer owner thread id recorded when the shard was prewarmed.

Per-unit observation gradients (streaming, grad-only tapes)

Description

Computes sparse gradients of each observation **unit** contribution at a fixed full parameter vector x, without re-running inner_opt. Each unit is taped temporarily with config$hessian_sparsity = FALSE (sparse Jac only), reverse-mode differentiated, then discarded.

Usage

grad_obs_units(
  data,
  x,
  config = list(),
  units = NULL,
  batch_size = 32L,
  inner = TRUE,
  negative = TRUE
)

Arguments

Details

Units are ELGM stratum columns when data@elgm_matrix has columns; otherwise units are rows of the observation vector y.

Value

A dgCMatrix of size n_gamma x length(units) when inner = TRUE, or n_full x length(units) when inner = FALSE. Column j is the gradient for units[j].

Config note

hessian_sparsity = FALSE yields packs usable for grad(..., inner = TRUE) but not for Hessian / inner_opt LDL paths.

See Also

grad, inner_opt, ad_pack

Whether this build was compiled with OpenMP support.

Description

Whether this build was compiled with OpenMP support.

Usage

has_openmp()

Value

TRUE if OpenMP was enabled at compile time, otherwise FALSE.

Build grouped Hessian sparsity structures

Description

Constructs Hessian sparsity templates (global/outer and inner-gamma) and per-group maps from shard sparsity patterns.

Usage

hessian_map(sparsity_list, n_beta, n_gamma, n_theta)

Arguments

Value

A list with outer, inner, optional chol_inner_list (symbolic LDL pattern for C++: L1, Linv, perm, perm (0-based; inner Hessian uses index1 = FALSE), perm_inv (0-based), half_H_inv, H_inv), optional trace_columns (added in ad_pack(), not by this function), chol_inner (dCHMsimpl for R), map_outer, map_inner, and sizes (named beta/gamma/theta; consumed internally by ad_pack()).

Place each shard's hinv_trace_index() result on its AD group.

Description

Observation densities split into several groups and must not define an index. A quadratic-form density is one group; its integer vector is stored at that group and every other group stays NULL.

Usage

hinv_q_index_list(density_shards, group_counts, H_inv)

Positions of a shard's quadratic form inside the joint inverse Hessian

Description

Called from ad_pack() after the symbolic H_inv pattern exists. The default method returns NULL, which leaves that shard on the trace_hinv_t column walk. A backend method may return an integer vector of 0-based slots into H_inv@x, one per stored entry of the shard's precision, with -1L where the pair is absent from the pattern.

Usage

hinv_trace_index(shard, H_inv)

## S4 method for signature 'density_data'
hinv_trace_index(shard, H_inv)

Arguments

Value

An integer vector, or NULL.

IID Random Effects Term

Description

Creates and manages IID (independent and identically distributed) random effects terms.

Usage

iid(
  x,
  init = .my_theta_init,
  lower = .my_theta_lower,
  upper = .my_theta_upper,
  parscale = .my_theta_parscale,
  log = TRUE
)

## S4 method for signature 'iid'
design(term, data)

## S4 method for signature 'iid'
precision(term, data)

## S4 method for signature 'iid'
random_info(term, data)

## S4 method for signature 'iid'
theta_info(term)

Arguments

Value

An 'iid' term object (in a named list).

Functions

• design(iid): Creates design matrix for iid term

• precision(iid): Creates precision matrix for iid term

• random_info(iid): Extracts random effects information for iid term

• theta_info(iid): Extracts theta parameter information for iid term

Per-unit influence scores for an inner linear functional

Description

Freezes outer parameters at the MLE, re-optimizes \gamma, forms v = H^{-1}a from the fit-pack inner Hessian, and streams one-unit observation gradients via grad_obs_units (temporary tapes). Scores are \psi_u = -v^\top g_u.

Usage

influence_scores(
  fit,
  a,
  batch_size = 32L,
  control = list(report.level = 0, report.freq = 0),
  verbose = FALSE
)

Arguments

Value

A list with phi, gamma_hat, H, v, G, psi, x_full, and n_units.

See Also

phi_gamma, amip_set, grad_obs_units, ad_pack_drop

Inner optimization over gamma using trust-region CG (sparse)

Description

Runs the inner optimization over \gamma using the trustOptim sparse trust-region CG solver on a pre-built ad_pack handle.

Usage

inner_opt(
  parameters,
  gamma,
  ad_pack,
  control = NULL,
  deriv = FALSE,
  verbose = FALSE,
  return_hessians = TRUE
)

Arguments

Details

Negates tape log-density values in C++ for minimization via trustOptim.

Value

A list with log_lik, neg_log_lik, parameters, full_parameters, inner_opt, gradient (list inner/outer), and hessian (list inner/outer/chol_inner/half_log_det; with deriv=TRUE and return_hessians=TRUE also half_H_inv, H_inv, trace3). With deriv=TRUE, also deriv (profile pieces including d_neg_log_lik) and dU. Objective and derivatives use the negative log-density convention.

intercept Model Term

Description

Creates and manages intercept model terms for fixed effects. With a finite sd, the intercept is treated as a random effect with a known-SD Gaussian prior (N(0, sd^2)), using random_diagonal without a theta parameter.

Usage

intercept(
  init = .my_beta_init,
  lower = .my_beta_lower,
  upper = .my_beta_upper,
  parscale = .my_beta_parscale,
  sd = Inf
)

## S4 method for signature 'intercept'
design(term, data)

## S4 method for signature 'intercept'
beta_info(term, data)

## S4 method for signature 'intercept'
precision(term, data)

## S4 method for signature 'intercept'
random_info(term, data)

Arguments

Details

intercept Model Term

Value

A intercept term object

Functions

• design(intercept): Creates design matrix for intercept term

• beta_info(intercept): Extracts beta parameter information for intercept term

• precision(intercept): Precision for known-SD intercept prior

• random_info(intercept): Random info for known-SD intercept prior

Methods

The following methods are available for 'intercept' objects:

design(term, data)

Creates design matrix for intercept term

precision(term, data)

Creates precision matrix for intercept term

theta_info(term)

Extracts theta parameter information

beta_info(term, data)

Extracts beta parameter information

random_info(term, data)

Extracts random effects information

Examples

# Create a intercept term
intercept_term <- intercept(init=2)

Integrated Wiener Process Term Constructor

Description

Creates an integrated Wiener process (IWP) model term. The basis is divided by the mean knot spacing to the power p, and the penalty weights are the knot spacings divided by that mean spacing. Design entries and the standard-deviation parameter are then unchanged if x and the knots are rescaled by the same factor.

Usage

iwp(
  x,
  p = 2,
  ref_value = 0,
  knots,
  range = NULL,
  init = .my_theta_init,
  lower = .my_theta_lower,
  upper = .my_theta_upper,
  parscale = .my_theta_parscale,
  log = TRUE,
  boundary_is_random = TRUE,
  include_poly = TRUE
)

Arguments

Value

A list containing the 'iwp' term object and optionally polynomial terms.

Examples

knots <- seq(0, 1, length.out = 5)
iwp_term <- iwp(x = "age", knots = knots)
iwp_term

Integrated Wiener Process Term

Description

Creates and manages integrated Wiener process (IWP) model terms.

Usage

## S4 method for signature 'iwp'
design(term, data)

## S4 method for signature 'iwp'
precision(term, data)

## S4 method for signature 'iwp'
theta_info(term)

## S4 method for signature 'iwp'
random_info(term, data)

Arguments

Value

design

Numeric design matrix for the integrated Wiener process basis, with one column per basis function and one row per observation. Entries are the basis evaluated at data.

precision

Sparse precision matrix for those basis coefficients. The penalty is the integrated Wiener process prior.

theta_info

One-row data frame describing the standard-deviation parameter: label, initial value, bounds, parameter scale, and whether optimization is on the log scale.

random_info

Data frame with one row per random coefficient, including gamma_label, the name of that coefficient in a fitted model.

Functions

• design(iwp): Creates design matrix for IWP term

• precision(iwp): Creates precision matrix for IWP term

• theta_info(iwp): Extracts theta parameter information for IWP term

• random_info(iwp): Extracts random effects information for IWP term

Methods

The following methods are available for 'iwp' objects:

design(term, data)

Creates design matrix for IWP term

precision(term, data)

Creates precision matrix for IWP term

theta_info(term)

Extracts theta parameter information

beta_info(term, data)

Extracts beta parameter information

random_info(term, data)

Extracts random effects information

Mean positive gap between finite knots.

Description

Mean positive gap between finite knots.

Usage

iwp_mean_spacing(knots)

Extract inner inverse Cholesky factor from Laplace output

Description

Returns H^{-1/2} for the inner (random-effect) Hessian from the output of log_lik_laplace. Accepts the flat half_H_inv matrix when deriv = TRUE, or reconstructs it from hessian$chol_inner.

Usage

laplace_half_H_inv(laplace)

Arguments

Value

A sparse or dense matrix H^{-1/2} with nrow = Ngamma.

See Also

log_lik_laplace, rmvnldl

Laplace-approximate log-likelihood (positive convention)

Description

Thin wrapper around log_lik_laplace that returns the positive profile log-likelihood and its gradient with respect to the outer parameters c(beta, theta). Prefer this over reading inner_opt(...)$deriv$d_neg_log_lik when building custom outer objectives (e.g. modular / multi-pack likelihoods).

Usage

laplace_loglik(
  ad_pack,
  par,
  gamma,
  config = list(),
  control_inner = list(),
  deriv = TRUE
)

Arguments

Value

A list with components

log_lik

Scalar Laplace-approximate marginal log-likelihood.

gradient

Numeric vector \nabla_\psi \ell, or NULL when deriv = FALSE.

gamma

Inner modes at the optimum.

See Also

log_lik_laplace, inner_opt, outer_fn

Latch (and optionally raise) the process-wide parallel team size.

Description

For requested > 1, updates the process high-water mark. On Windows, returns that mark after raising it so OpenMP/CppAD team setup never shrinks within one R process (a decrease aborts under rtools/libomp). At eval time empty trailing thread groups are padded to this width; shard owners from ad_pack() still follow the requested num_threads. On other platforms the requested count is returned unchanged. Serial requested = 1 is always returned unchanged.

Usage

latch_parallel_threads(requested)

Arguments

Value

Integer effective OpenMP/CppAD team size (may exceed requested on Windows after a larger parallel team was used in this process).

Linear Model Term

Description

Creates and manages linear model terms for fixed effects. With a finite sd, coefficients are random effects with a known-SD Gaussian prior (N(0, sd^2)).

Usage

linear(
  x,
  init = .my_beta_init,
  lower = .my_beta_lower,
  upper = .my_beta_upper,
  parscale = .my_beta_parscale,
  sd = Inf
)

## S4 method for signature 'linear'
design(term, data)

## S4 method for signature 'linear'
beta_info(term, data)

## S4 method for signature 'linear'
precision(term, data)

## S4 method for signature 'linear'
random_info(term, data)

Arguments

Details

Linear Model Term

Value

A linear term object

Functions

• design(linear): Creates design matrix for linear term

• beta_info(linear): Extracts beta parameter information for linear term

• precision(linear): Precision for known-SD linear prior

• random_info(linear): Random info for known-SD linear prior

Methods

The following methods are available for 'linear' objects:

design(term, data)

Creates design matrix for linear term

precision(term, data)

Creates precision matrix for linear term

theta_info(term)

Extracts theta parameter information

beta_info(term, data)

Extracts beta parameter information

random_info(term, data)

Extracts random effects information

Examples

# Create a linear term
linear_term <- linear(x = "temperature")

Log-likelihood with inner Laplace optimization

Description

Evaluates the (profiled) log-likelihood for a hierarchical model by solving an inner optimization problem for the latent vector gamma (e.g. random effects) given outer parameters x (typically beta and theta). The function delegates the inner optimization to inner_opt() from the selected backend package (by default adlaplace).

Usage

log_lik_laplace(
  x,
  config = list(),
  gamma,
  control = list(report.level = 4, report.freq = 1),
  ad_pack,
  data,
  package = c(config[["package"]], "adlaplace")[1L],
  deriv = FALSE,
  return_hessians = TRUE
)

Arguments

Details

The default adlaplace backend uses a single AD handle. This function passes that handle to inner_opt(..., ad_pack = ad_pack). Parameter block sizes are read from ad_pack@sizes; config is not required to contain beta, gamma, or theta when ad_pack is supplied.

When deriv=FALSE, the return value mirrors inner_opt() with inner_opt, top-level gradient, and hessian. When deriv=TRUE, the profile gradient is assembled in C++ inside inner_opt; the same list is augmented with deriv and profile-derivative pieces in gradient$outer.

Value

A list. With deriv=FALSE, same structure as inner_opt (without parameters). With deriv=TRUE, additionally:

deriv

Data frame of profile derivative pieces.

gradient$outer$dU, gradient$outer$trace3

Intermediate profile objects.

Note

The Laplace approximation assumes the inner Hessian is positive definite at the optimum. If inner_opt fails to converge or returns a non-invertible Hessian, log_lik_laplace() issues a warning or error.

See Also

ad_pack, inner_opt

London daily mortality and PM10 (case-crossover example)

Description

Daily all-cause mortality and air-quality covariates for London, 2002–2006, from the case-crossover example of Tobias et al. (2024).

Usage

london

Format

A data frame with 1826 rows and columns:

city

City label ("London").

date

Calendar date.

year, month, day

Integer calendar components.

dow

Day of week (integer encoding from the source file).

all

Daily all-cause death count.

tmean

Mean temperature (degrees Celsius).

rh

Relative humidity (percent).

pm10

PM10 concentration.

Source

https://github.com/aureliotobias/casecrossover

References

Tobias, A., et al. (2024). Case-crossover designs for air pollution epidemiology (example data).

Examples

data("london", package = "adlaplace")
head(london)

Longest-processing-time thread assignment

Description

Longest-processing-time thread assignment

Usage

lpt_assign_owners(costs, num_threads)

Arguments

Value

Integer vector of length n_shards with 0-based owner thread ids.

Posterior mean of random effects via multivariate quadrature

Description

Approximates E[\gamma \mid y] by importance-weighted quadrature over \gamma, using a Gaussian proposal N(\hat\gamma, H^{-1}) (Laplace covariance) and correcting with joint_log_dens.

Usage

mean_mvquad(
  parameters,
  mode,
  cov,
  ad_pack,
  n = 3L,
  type = "nHN",
  ndConstruction = "sparse",
  dec.type = 2L
)

Arguments

Details

Builds a sparse nHN grid with mvQuad, rescales to the Laplace proposal, then forms \sum_i w_i \gamma_i \exp(\ell(\gamma_i) - \log q(\gamma_i)) / \sum_i w_i \exp(\ell(\gamma_i) - \log q(\gamma_i)) with log-sum-exp stabilization.

Value

Numeric vector of length n_gamma.

See Also

log_lik_laplace, inner_opt, joint_log_dens, rmvnldl

Model Term Generics

Description

Generic functions for extracting model matrices and parameter information from model terms.

When non-NULL, model_data emits an extra ad_kind = "parameters" shard (e.g. FEM log-determinant) alongside the random shard. Default is NULL (no companion).

Usage

design(term, data)

precision(term, data)

theta_info(term)

beta_info(term, data)

random_info(term, data)

elgm_matrix(term, data)

## S4 method for signature 'model_term'
design(term, data)

## S4 method for signature 'model_term'
precision(term, data)

## S4 method for signature 'model_term'
theta_info(term)

## S4 method for signature 'model_term'
beta_info(term, data)

## S4 method for signature 'model_term'
random_info(term, data)

extra_density(term)

## S4 method for signature 'model_term'
extra_density(term)

Arguments

Value

Method-specific return values.

Character density name, or NULL.

Assemble model terms and per-shard density_data objects

Description

Parses a formula, builds design matrices, and returns formula terms plus density_data shards for the observation density, parameter densities, and each random-effect term.

Usage

model_data(formula, data, verbose = FALSE, na_omit = TRUE)

Arguments

Value

A list with:

terms

Named list of term objects from collect_terms.

data

Output of term_data_setup, including elgm_matrix when an observation term defines an elgm_matrix method.

observations

Named list of observation density_data shards.

parameters

Named list of parameter density_data shards.

random

Named list of random-effect density_data shards.

Examples

set.seed(1)
dat <- data.frame(
  y = rpois(30, 2),
  x1 = rnorm(30),
  x2 = runif(30)
)
md <- model_data(
  y ~ x1 + iwp(x2, p = 2, knots = seq(0, 1, length.out = 5)),
  data = dat
)
names(md$terms)

Base Model Class

Description

The base S4 class that all specific model term classes inherit from. Subclasses that need hierarchical grouping should include a by slot (either character or by_group type) and optionally by_levels and by_labels slots.

Slots

name

Character vector of length 1, the term name

label

Character scalar term label reused across metadata builders.

formula

Formula object for the term

knots

Numeric vector of knot locations (for spline terms)

ref_value

Numeric reference value for centering

p.order

Integer polynomial order

init

Numeric vector of initial values

lower

Numeric vector of lower bounds

upper

Numeric vector of upper bounds

parscale

Numeric vector of parameter scales for optimization

log

Logical; optimize/tape each theta on the log scale (recycled to the number of thetas for multi-parameter terms).

model_role

Factor indicating term role ("fixed", "random", or "response")

density

Registered AD density name for ad_pack_ptr(), or NA when the term does not map to a density shard.

ad_kind

Shard kind passed to ad_pack_ptr() ("observations", "parameters", "random", or NA).

package

Package whose .so records AD tapes for this term (default "adlaplace" for built-in densities).

Number of AD shards in an ad_pack_ptr handle

Description

Number of AD shards in an ad_pack_ptr handle

Usage

n_groups(handle)

Arguments

Value

Integer count of groups (shards).

Number of removable observation units

Description

Units are ELGM stratum columns when elgm_matrix has columns; otherwise rows of the observation vector y. Same domain as grad_obs_units.

Usage

n_obs_units(data)

Arguments

Value

Integer length-1 count of units.

Negative binomial observation term

Description

Model term for the observation-level negative binomial log density registered as nbinom_obs.

Creates a response term wired to the nbinom_obs AD shard.

Usage

nbinom(
  x,
  init = .my_theta_init,
  lower = .my_theta_lower,
  upper = .my_theta_upper,
  parscale = .my_theta_parscale,
  log = TRUE
)

## S4 method for signature 'nbinom'
theta_info(term)

Arguments

Details

Negative binomial observation model term

Value

A nbinom object.

Functions

• theta_info(nbinom): Theta info for the overdispersion parameter.

Slots (inherited from model)

ad_pack

Character scalar "nbinom_obs".

ad_kind

Character scalar "observations".

Fill missing config$gamma before ad_pack_ptr() tape build.

Description

Fill missing config$gamma before ad_pack_ptr() tape build.

Usage

normalize_config_for_ptr(config, data, kind = NULL)

Normalize random_mult precision list so Q is dsCMatrix.

Description

Normalize random_mult precision list so Q is dsCMatrix.

Usage

normalize_random_mult_precision(precision)

Partition observations into AD shards by sparsity pattern

Description

The function partitions columns into num_shards shards using quantiles of the first right singular vector, optionally using RSpectra for efficiency when available.

Usage

obs_groups(A, elgm_matrix, num_shards, min_shards = 0)

Arguments

Details

The resulting grouping is returned as a pattern sparse matrix whose columns correspond to shards and whose nonzeros mark membership. Loadings are used only to form the partition. Shards are ordered from most heterogeneous to most homogeneous.

If the RSpectra package is available, the leading singular vector is computed using RSpectra::svds; otherwise, a full singular value decomposition via svd is used.

Shard boundaries are defined by empirical quantiles of the loadings. Shards are subsequently reordered so that shards with larger within-shard variability appear first.

Value

A sparse matrix of class "ngCMatrix" (from Matrix), with one column per shard and one row per observation (or ELGM stratum). Nonzero entries mark shard membership; there is no @x slot.

Examples

set.seed(1)
A <- matrix(rnorm(100), 20, 5)
G <- obs_groups(A, num_shards = 3)
G

Build obs_groups for a subset of observation units

Description

Restricts which observation units are recorded on observation tapes. Domain indices are **0-based** (ELGM stratum columns, or rows of y).

Usage

obs_groups_units(
  n_domain,
  units = NULL,
  grouping = c("identity", "together", "filter"),
  obs_groups = NULL
)

Arguments

Value

A ngCMatrix suitable for config$obs_groups.

See Also

n_obs_units, grad_obs_units, ad_pack_drop

Outer objective and gradient wrappers

Description

Convenience wrappers around log_lik_laplace for use with outer optimizers that expect fn / gr callbacks. Both functions solve (or warm-start) the inner problem over gamma via log_lik_laplace and update a mutable cache environment with the latest inner solution.

Usage

outer_fn(x, config, cache, ad_pack, control_inner = list(), ...)

outer_gr(x, config, cache, ad_pack, control_inner = list(), ...)

Arguments

Details

L-BFGS-B requests fn then gr at the same parameter vector. Both wrappers share one inner_opt(..., deriv = TRUE, return_hessians = FALSE) evaluation at that point: the objective is free once the profile gradient has been computed. Cached numerics are stored in cache$fg_x, cache$neg_log_lik, and cache$d_neg_log_lik and reused until x changes.

outer_fn()

Returns the scalar objective negative log likelihood.

outer_gr()

Returns the gradient.

Value

• outer_fn: a numeric scalar neg_log_lik.

• outer_gr: gradient of neg_log_lik (same sign as outer_fn).

See Also

log_lik_laplace

Examples

set.seed(1)
n <- 80
g <- factor(sample(4, n, replace = TRUE))
x <- rbinom(n, 1, 0.5)
eta <- 0.5 + 0.4 * x + rnorm(4, sd = 0.3)[g]
nb_sd <- 0.3
z <- rgamma(n, nb_sd^(-2), nb_sd^(-2))
dat <- data.frame(y = rpois(n, exp(eta) * z), x = x, g = g)
md <- model_data(
  nbinom(y, lower = 1e-9, init = 0.15) ~ x + iid(g, init = 0.3),
  data = dat
)
config <- list(
  transform_theta = TRUE,
  obs_groups = obs_groups(md$term_data$A, num_shards = 5),
  verbose = FALSE
)
ap <- ad_pack(md, config)
cache <- new.env(parent = emptyenv())
x0 <- md$term_data$info$parameters$init
outer_fn(x = x0, config = config, cache = cache, ad_pack = ap)
outer_gr(x = x0, config = config, cache = cache, ad_pack = ap)

Linear functional of the inner field \gamma

Description

Builds a named numeric vector a aligned to fit$model_data$term_data$info$gamma$gamma_label so that \phi(\gamma) = a^\top\gamma.

Usage

phi_gamma(fit, term = NULL, new_x = NULL, coef = NULL)

Arguments

Value

Named numeric vector of length n_gamma.

See Also

influence_scores, amip_set, design, sim_random

Poisson observation term

Description

Model term for the observation-level Poisson log density with a log link, registered as poisson_obs. For each observation with linear predictor eta = X beta + A gamma and optional offset o (via config$offset), the mean is exp(eta + o) and the log density contribution is y * eta - exp(eta + o) + y * o - lgamma(y + 1). There are no observation-level hyperparameters.

Creates a response term wired to the poisson_obs AD shard. There are no observation-level hyperparameters. An exposure or expected-count offset is supplied through config$offset when building AD shards (for example log(E) in disease mapping).

Usage

pois(x)

Arguments

Details

Poisson observation model term

Value

A poisson object.

Slots (inherited from model)

ad_pack

Character scalar "poisson_obs".

ad_kind

Character scalar "observations".

Drop incomplete rows and sort case-crossover strata

Description

Drop incomplete rows and sort case-crossover strata

Usage

prepare_model_rows(data, terms, verbose = FALSE)

Arguments

Value

Subset and reordered data frame.

Hyperparameter prior term

Description

Places a prior on an existing theta (SD) parameter. Does not create a new theta; it emits an ad_kind = "parameters" shard that contributes to the joint log density.

Usage

prior(theta = 0L, dist = "exp", median = 1, rate = NULL)

Arguments

Details

Hyperparameter prior term

Value

A prior term object.

Profile per-shard Reverse3 cost with dummy LinvPt directions

Description

Serial timing helper for ad_pack(..., reorder_shards = "third"). Uses the symbolic half_H_inv CSC pattern (values forced to 1) and trace_columns as LinvPtColumns. Does not require owner threads or a numeric Cholesky.

Usage

profile_shard_trace3_times(
  handle,
  x,
  half_H_inv,
  trace_columns,
  n_rep = 4L,
  n_warmup = 1L
)

Arguments

Value

Numeric vector of length n_shards (mean elapsed seconds).

Align Reference Value to Nearest Knot

Description

Aligns a reference value to the nearest knot in a set of knots. This is useful for ensuring that reference values used in polynomial terms are aligned with the actual knot positions.

Usage

ref_align(ref_value, knots)

Arguments

Value

The knot value that is closest to the specified reference value.

Examples

# Align reference value to nearest knot
knots <- c(10, 20, 30, 40, 50)
adlaplace::ref_align(23, knots) # Returns 20

Replace data-column names in model terms

Description

Rewrites the data-column references on model terms (and formulas / term lists) according to a named character map. Used by modular-inference workflows that share one model structure across exposure sample columns (e.g. sqrt_pm -> sqrt_pm_s3).

Usage

replace_vars(x, map)

## S4 method for signature 'model_term'
replace_vars(x, map)

## S4 method for signature 'dirichlet_multinom'
replace_vars(x, map)

## S4 method for signature 'formula'
replace_vars(x, map)

## S4 method for signature 'list'
replace_vars(x, map)

Arguments

Value

An object of the same type as x, with column references updated and labels / formulas regenerated where applicable.

Examples

terms <- collect_terms(nbinom(y) ~ fpoly(sqrt_pm, p = 1))
replace_vars(terms, c(sqrt_pm = "sqrt_pm_s1"))

Resolve inner gamma starting values for optimization.

Description

Resolve inner gamma starting values for optimization.

Usage

resolve_gamma_start(config, cache = NULL, num_gamma, gamma = missing())

Simulate from a multivariate normal using LDL of the precision matrix

Description

Draws n samples from N(\mu, H^{-1}) where H is the precision matrix, given its LDL factorization (unit lower L1, diagonal D, and permutation perm).

Usage

rmvnldl(n, mean = 0, chol_prec, fit)

Arguments

Details

Uses H^{-1/2} = P^\top L^{-1} D^{-1/2} (same as reformat_chol() in log_lik_deriv), and rnorm only: x = \mu + H^{-1/2} z, z \sim N(0, I).

Value

An n-by-p matrix; each row is one draw. When n = 1, a length-p vector. Column / element names come from names(mean) when present and of length p.

See Also

Cholesky, log_lik_laplace, adlaplace

Random Polynomial Model Term

Description

Creates and manages random polynomial model terms.

Usage

rpoly(x, p = 2, ref_value = 0, sd = Inf, basis_scale = 1)

## S4 method for signature 'rpoly'
design(term, data)

## S4 method for signature 'rpoly'
precision(term, data)

## S4 method for signature 'rpoly'
random_info(term, data)

Arguments

Details

Random Polynomial Model Term

Value

A 'rpoly' term object.

A design matrix for the random polynomial term.

A precision matrix for the random polynomial term.

A data frame containing random effects information for the random polynomial term.

Functions

• design(rpoly): Design method for rpoly objects

• precision(rpoly): Precision method for rpoly objects

• random_info(rpoly): Random info method for rpoly objects

Methods

The following methods are available for 'rpoly' objects:

design(term, data)

Creates design matrix for rpoly term

precision(term, data)

Creates precision matrix for rpoly term

theta_info(term)

Extracts theta parameter information

beta_info(term, data)

Extracts beta parameter information

random_info(term, data)

Extracts random effects information

Simulate linear predictors on a new covariate grid

Description

Draws random effects from the Laplace approximation and evaluates the linear predictor on x by summing fixed (\beta) and random (\gamma) contributions from each model term.

Usage

sim_fit(x, data, fit, n = 500L)

Arguments

Value

The combined linear predictor: a nrow(x) by n matrix of simulated values (fixed-effect contributions are constant across columns).

Simulate IID random effects with pointwise and global intervals

Description

Identifies IID components in an adlaplace fit, forms their marginal Laplace covariance (Woodbury identity when mixed with other random effects), draws samples, and summarizes pointwise quantiles. When GET is available, also returns a global central region across levels.

Usage

sim_iid(fit, n = 500L, probs = c(0.1, 0.5, 0.9), coverage = 0.8, ...)

Arguments

Details

IID rows are those with model == "iid" in fit$model_data$term_data$info$gamma. The marginal covariance uses WoodburyMatrix when installed; otherwise the IID block of fit$details$hessian$H_inv (if present) or a dense partitioned solve.

Value

A named list with one element per IID term. Each element is a list with x (basis / factor levels), sim (levels-by-draws matrix), pointwise (quantile matrix with columns named by probs), and global (matrix with columns lo, central, hi, or NULL if GET is unavailable). Attributes var_iid, iid_index, and iid_labels store the joint marginal covariance and index metadata.

See Also

sim_random, rmvnldl, iid

Simulate a random-effect contribution for one covariate

Description

Builds a prediction grid for variable x, forms the design matrix for all random terms on that variable, and multiplies by Laplace draws of the random effects. Intended for plotting smooths (iwp / rpoly) and similar random contributions from an adlaplace fit.

Usage

sim_random(
  x,
  fit,
  new_x,
  gamma_sims,
  num_grid = 101L,
  num_sim = 100L,
  probs = c(0.1, 0.5, 0.9),
  coverage = 0.8,
  transform = identity
)

Arguments

Details

Only terms with model_role == "random" whose @name equals x are included (so companion rpoly bases for an iwp are included, while fixed fpoly bases are not). Column names of the design must match columns of gamma_sims.

IID prediction grids use factor levels stored in the gamma info table; levels not present at fit time are not invented here.

Value

A list with x (numeric or character grid), sim (an length(x) by draws matrix after transform), pointwise (quantile matrix), and global (GET envelope with columns lo, central, hi, or NULL if GET is unavailable).

See Also

rmvnldl, sim_fit, sim_iid, design

Parameter-block row counts for an density_data

Description

Parameter-block row counts for an density_data

Usage

sizes(x, ...)

## S4 method for signature 'density_data'
sizes(x)

Arguments

Value

Named list with n_beta, n_gamma, and n_theta.

Per-shard layout for AD density evaluation

Description

Per-shard layout for AD density evaluation

Usage

strip_term_name(x)

Arguments

Value

Character vector with surrounding quotes removed.

Sparse inverse subset of a symmetric positive-definite matrix

Description

Entries of Q^{-1} on the sparsity pattern of the Cholesky factor, in the original variable order. Positions the factor does not store are absent, not zeros of the inverse.

Usage

takahashi_davis(Q, chol = NULL)

Arguments

Details

C++ callers 'LinkingTo: adlaplace' and include 'adlaplace/takahashi_impl.hpp'. That function is 'adlaplace::chol::takahashi_davis' and takes the unit-lower factor L and the LDL diagonal D. A caller that has only Q factors it with 'chol_update_csc' first.

Value

A 'dgCMatrix' of the selected inverse.

Build design matrices and parameter metadata from a formula

Description

Build design matrices and parameter metadata from a formula

Usage

term_data_setup(formula, data, verbose = FALSE)

Arguments

Value

List with terms, y, A, X, data, info. The info$theta data frame includes a logical log column: TRUE means the parameter is stored on the log scale when config$transform_theta is TRUE (default for all rows unless a model term sets otherwise). info$parameters stacks beta and theta rows with the same log column (FALSE for fixed effects); init, lower, and upper are already on the optimization (log) scale where log is TRUE.

Per-shard column indices into half_H_inv for trace_hinv_t

Description

Built from group_sparsity and the symbolic half_H_inv pattern in chol_inner_list. Added to that list in new_ad_pack_from_ptr().

Usage

trace_columns_from_pattern(group_sparsity, n_beta, n_gamma, half_H_inv_pat)

Arguments

Value

Sparse ngCMatrix (index1 = FALSE) with dims = c(n_gamma, length(group_sparsity)).

Touch OpenMP from the main thread after dyn.load.

Description

Intended for fresh R processes (e.g. R CMD check vignette re-builds) so Homebrew libomp TLS is initialized before the first multi-thread CppAD session.

Usage

warm_openmp_runtime()