Package 'adlaplace'

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.

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.

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

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.

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

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

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.

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)

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.

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.

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

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

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.

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.

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)

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

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

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.

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

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.

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:

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:

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.

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".

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

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.

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)

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)

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".

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.

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

Description

Whether this build was compiled with OpenMP support.

Usage

has_openmp()

Value

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

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.

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

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

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.

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)

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

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

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

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

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")

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

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)

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

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)

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).

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.

Description

Number of AD shards in an ad_pack_ptr handle

Usage

n_groups(handle)

Arguments

Value

Integer count of groups (shards).

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.

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".

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

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

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)

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

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".

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.

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"))

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

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

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).

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

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

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.

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.

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.