gam_design

Penalized generalized additive model optimal design

Description

An R6 class for managing penalized locally optimal experiment designs for generalized additive models. This class encapsulates the generalized additive model at a set of nominal parameter values, the discrete design grid, optimization of an approximate continuous design, augmentation of an existing design, and optimal rounding to exact integer designs.

Optimal design problem

The following penalized optimization problem is solved by $optimize():

w^* \ = \ \arg \max_{w} \Phi(M(w | X, v, \theta, F, g) + \lambda I) - \alpha \cdot c^T w

where,

$w = \{w_1, \ldots, w_n \}$ is the vector of design weights.
$c = \{c_1, \ldots, c_n \}$ is the cost assigned to each design weight with scaling/penalty parameter $\alpha$ .
$v = \{v_1, \ldots, v_n \}$ is a vector of fixed user weights.
$M(w | X, v, \theta, F, g)$ is the information matrix depending only on the weights $w$ and $v$ , and the design matrix $X$ by default. If a non-trivial distribution family $F$ is chosen, the information matrix also depends on the nominal parameter values $\theta$ , the chosen (exponential) distribution family $F$ and the link function $g$ .
$\Phi$ is the chosen optimality criterion, generally a function of the information matrix.

By default, the error distribution family $F$ is Gaussian with identity link function $g$ . This implies that the information matrix (conditional on the error variance) is independent of the model parameters and given by:

M(w | X, v) \ \propto \ \sum_{i = 1}^n w_i v_i x_i x_i^T

If $F$ is not Gaussian or Gaussian with a non-canonical link function $g$ , the information matrix is given by:

M(w | X, v, \theta, F, g) \ = \ \sum_{i = 1}^n w_i v_i u(x_i, \theta) x_i x_i^T

with weights $u(x_i, \theta) = \frac{1}{V(\mu_i)[g'(\mu_i)]^2}$ depending on the parameter values $\theta$ . Here, $\mu = E[Y] = b'(\zeta)$ and $V(\mu) = b''(\zeta)$ , with $b(\zeta)$ is the cumulant function of the exponential family $F$ . See glm_design for the exact details.

Augmented optimal design

Given a vector of initial design weights $w_0$ and associated contribution weight $\gamma$ , the following augmented penalized optimization problem is solved by $augment():

w^* \ = \ \arg \max_{w} \Phi(M(w | \gamma, w_0, X, v, \theta, F, g) + \lambda I) - \alpha \cdot c^T w

using the augmented information matrix,

M(w | \gamma, w_0, X, v, \theta, F, g) \ = (1 - \gamma) M(w_0 | X, v, \theta, F, g) + \gamma M(w | X, v, \theta, F, g)

Stratified optimal design

The $strata() method stratifies the design by requiring that individual strata $X_1, \ldots, X_S$ are represented proportionally with given weights $p_1, \ldots, p_S$ , ( $\sum_i p_i = 1$ ). If stratification is applied, the following optimization problem is solved by $optimize():

\begin{array}{rl} w^* \ = \ & \arg \max_{w} \Phi(M(w | X, v, \theta, F, g) + \lambda I) - \alpha \cdot c^T w, \\ & \mathrm{subject\ to} \ \sum_{i \in X_s} w_i = p_s, \quad \quad \mathrm{for} \ s = 1,\ldots, S \end{array}

using the same notation as in the ‘Optimal design problem’ section above. This means that the total weight assigned to the design points in stratum $X_i$ is equal to $p_i$ for each $i = 1,\ldots,S$ . Similarly, if stratification is applied, $augment() solves the augmented optimization problem with stratification constraints.

Sparse optimal design

The $sparsify() method solves the penalized optimization problem with an extra sparseness penalty $P(w)$ in the objective to promote sparsity in the approximate design $w^*$ . See lm_design for the available penalty functions.

Optimality criteria

The optimality criteria $\Phi(\cdot)$ supported by $optimize(), $augment() and $sparsify() are: D-, A-, I-, G-, and alias-optimality, see lm_design for the definitions.

Super classes

lm_design -> glm_design -> gam_design

Public fields

formula: a GAM model defined as a one-sided formula.
theta: a named vector or list of nominal parameter values at which the information matrix $M(w | X, \theta, F, g)$ is evaluated. Only used if $family is not gaussian(link = "identity") (the default).
family: a description of the model distribution family and link function. See family for details of family functions.
data: a data frame or matrix containing all design points (rows) at which formula can be evaluated.
weights: an optional numeric vector containing fixed user weights associated to each design point (row) in data.
cost: an optional numeric vector containing costs associated to each design point (row) in data.
design_weights: a numeric vector with optimal design weights populated after optimization with $optimize(), $augment(), or $sparsify(). The attribute criterion lists the used design criterion and the achieved optimal criterion value.
strata: a list with elements strata and proportions imposing a partition of the design points into individual (disjoint) strata with given weight proportions.
alpha: the cost scale/penalty parameter in the penalized optimization problem, only used if $cost is available.
lambda: the ridge penalty added to the diagonal of the information matrix in the optimization problem. The ridge penalty is not used
upper: the design weight upper limit. The default is 1, corresponding to no upper limit constraint on the design weights.

Active bindings

X: the design matrix obtained via mgcv::gam(formula, data = data, fit = FALSE, ...) given the design points in data and the model formula.

Methods

`gam_design$new()`

Constructor. Initializes the experimental design to optimize with a linear model formula, discrete design grid, and optionally a distribution family with link function (and nominal set of parameter values).

Usage

gam_design$new(
  formula,
  data,
  theta = NULL,
  family = NULL,
  weights = NULL,
  cost = NULL,
  contrasts = NULL,
  alpha = 1,
  lambda = 0,
  upper = 1,
  ...
)

Arguments

formula: a GAM model defined as a one-sided formula, see gam.
data: a data frame or matrix containing all design points (rows) at which formula can be evaluated.
theta: a named or unnamed vector or list of nominal parameter values at which the information matrix $M(w | X, \theta, F, g)$ is evaluated. Only needed if $family is not NULL.
family: a family function describing the error distribution and link function in the model. Defaults to the gaussian family with "identity" link function.
weights: an optional numeric vector containing fixed user weights associated to each design point (row) in data. Also accepts a one-sided formula made up of column names present in data. If NULL, the unweighted optimal design problem is solved.
cost: an optional numeric vector containing costs associated to each design point (row) in data. Also accepts a one-sided formula made up of column names present in data. If NULL, the standard (non-penalized) optimal design problem is solved.
contrasts: optional contrasts used when evaluating the linear model formula.
alpha: numeric value, cost scale/penalty parameter in the penalized optimization problem, only used if $cost is available. The default is 1.
lambda: numeric value, ridge penalty added to the diagonal of the information matrix in the penalized optimization problem. The default is 0.
upper: numeric value in (0, 1], design weight upper limit. The default is 1, corresponding to no constraint on the design weights.
...: any additional arguments passed to mgcv::gam(formula, data = data, fit = FALSE, ...) to evaluate the design matrix.

`gam_design$update()`

Updates the fields of an existing object. Inputs should be named formula, data, theta, family, weights, cost, alpha, lambda and/or upper respectively.

Usage

gam_design$update(...)

Arguments

...: named arguments formula, data, theta, family, weights, cost, alpha, lambda and/or upper.

`gam_design$stratify()`

Stratify the design by requiring that individual strata are represented proportionally. Call with strata = NULL to unset previously defined strata.

Usage

gam_design$stratify(strata, proportions = NULL)

Arguments

strata: a factor variable of the same length as the number of rows in data. The factor levels determine the partition of the design points into individual (disjoint) strata. Also accepts a one-sided formula made up of column names present in data, in which case the strata are defined by the interaction between terms in the model.frame obtained from the formula.
proportions: an optional numeric vector containing the weight proportions of the individual strata. Must be of the same length as the number of factor levels in strata. If NULL, each stratum is assigned the same weight proportion.

`gam_design$augment()`

Optimizes the augmented penalized approximate design given a vector of initial design weights

Usage

gam_design$augment(
  design_weights,
  gamma = 0.5,
  criterion = c("D", "A", "I", "G", "alias"),
  scs_control = scs_control_dflt(),
  verbose = FALSE,
  return_value = TRUE
)

Arguments

design_weights: a numeric vector with initial design weights, must be the same length as the number of design points in data. Also accepts a one-sided formula made up of parameters in theta and column names in data.
gamma: numeric contribution weight of the initial design in the augmented design problem, see also ‘Augmented optimal design’ above. Defaults to 0.5.
criterion: character string specifying the optimality criterion to use, see also ‘Details’. The following choices are supported:
"D" D-optimality (default).
"A" A-optimality.
"I" I-optimality.
"G" G-optimality.
"alias" alias-optimality.
scs_control: an optional list of control parameters to tune the SCS optimization. See scs_control_dflt for the available control parameters and their default values.
verbose: display verbose messages printed by the scs solver. Defaults to FALSE.
return_value: return the optimal design weights. If set to FALSE, the design weights should be retrieved from the design_weights field.


`"D"` D-optimality (default).
`"A"` A-optimality.
`"I"` I-optimality.
`"G"` G-optimality.
`"alias"` alias-optimality.

`gam_design$optimize()`

Optimizes the penalized approximate design based on the selected criterion.

Usage

gam_design$optimize(
  criterion = c("D", "A", "I", "G", "alias"),
  scs_control = scs_control_dflt(),
  verbose = FALSE,
  return_value = TRUE
)

Arguments

criterion: character string specifying the optimality criterion to use, see also ‘Details’. The following choices are supported:
"D" D-optimality (default).
"A" A-optimality.
"I" I-optimality.
"G" G-optimality.
"alias" alias-optimality.
scs_control: an optional list of control parameters to tune the SCS optimization. See scs_control_dflt for the available control parameters and their default values.
verbose: display verbose messages printed by the scs solver. Defaults to FALSE.
return_value: return the optimal design weights. If set to FALSE, the design weights should be retrieved from the design_weights field.


`"D"` D-optimality (default).
`"A"` A-optimality.
`"I"` I-optimality.
`"G"` G-optimality.
`"alias"` alias-optimality.

`gam_design$sparsify()`

Sparsify the penalized approximate design by iterative reweighting of the optimization problem.

Usage

gam_design$sparsify(
  max_iter,
  criterion = c("D", "A", "I", "G", "alias"),
  reweight_fun = c("entropy", "log-sum", "power"),
  eps = .Machine$double.eps,
  tol = .Machine$double.eps^(1/4),
  p = 0.5,
  design_weights = NULL,
  scs_control = scs_control_dflt(),
  return_value = TRUE,
  show_progress = TRUE
)

Arguments

max_iter: positive integer, maximum number of reweighted optimizations/iterations to run.
criterion: character string specifying the optimality criterion to use, see also ‘Details’. The following choices are supported:
"D" D-optimality (default).
"A" A-optimality.
"I" I-optimality.
"G" G-optimality.
"alias" alias-optimality.
reweight_fun: character reweighting function, one of "entropy", "log-sum", "power", see lm_design for the definitions.
eps: positive numeric value to ensure finite weights in reweighting functions.
tol: numeric tolerance used as a threshold for mean absolute relative difference in design weights between iterations to stop early.
p: positive numeric power with $0 < p < 1$ . Only used if reweight_fun = "power".
design_weights: (optional) vector of design weights used as initial approximate design.
scs_control: an optional list of control parameters to tune the SCS optimization. See scs_control_dflt for the available control parameters and their default values.
return_value: return the optimal design weights. If set to FALSE, the design weights should be retrieved from the design_weights field.
show_progress: logical, display a progress bar during the iterative reweighting procedure. Defaults to TRUE.


`"D"` D-optimality (default).
`"A"` A-optimality.
`"I"` I-optimality.
`"G"` G-optimality.
`"alias"` alias-optimality.

Details

The sparseness cost penalties are scaled by the same alpha tuning parameter as for a user-specific cost vector, see scs_control_dflt. If a cost vector c_1,…,c_n associated to the design weights is already specified. The sparseness penalty costs are multiplied by the existing cost vector.

`gam_design$round()`

Round penalized approximate design to an exact (integer) design. By default, the approximate design is initialized from the design_weights field, which is the optimal (approximate) design calculated by $optimize(). For method "optimal", ridge and cost penalties and upper limit constraints present in $lambda, $cost and $upper(if specified) are included in the criterion optimized by the point-exchange algorithm. To override the initial approximate (optimal) design, a vector of design weights can be passed using the design_weights argument. The vector must be the same length as the number of design points in data. If design_weights is set to FALSE, the initial design is generated by random sampling.

Usage

gam_design$round(
  m,
  method = c("optimal", "efficient", "dist-sum", "dist-min"),
  max_iter = 100,
  n_repeats = 10,
  replicates = TRUE,
  tol = .Machine$double.eps^(1/3),
  seed = NULL,
  design_weights = TRUE,
  augment_design = NULL,
  show_progress = TRUE,
  ...
)

Arguments

m: integer number of design points included in the exact design.

method: the method used to calculate the exact design. The following choices are supported:


`"optimal"` optimality criterion point-exchange algorithm (default).
`"efficient"` efficient rounding method of Pukelsheim & Reider (1992).
`"dist-sum"` maximize total pairwise (Euclidean) distance between design points.
`"dist-min"` maximize minimum pairwise (Euclidean) distance between design points.

max_iter: integer number of maximum iterations in point-exchange algorithm. Only used when method = "optimal", method = "dist-sum" or method = "dist-min".
n_repeats: integer number of (random) repeats of the rounding/exchange algorithm. The exact design with the highest criterion value is returned. Defaults to 10.
replicates: logical, allow replicate design points in the exact design, defaults to TRUE. Can also be an integer number, limiting the maximum number of allowed replicates per design point. This argument overrides the $upper value constraint if existing.
tol: numeric tolerance used as lower threshold to continue iterations in point-exchange algorithm.
seed: (optional) random seed used to generate jitter to break ties in rounding the initial exact design.
design_weights: (optional) vector of design weights used as initial approximate design, with length equal to the number of available design points. If TRUE (default), the approximate design is initialized from the design_weights field. If FALSE, the initial approximate design is generated by random sampling.
augment_design: (optional) integer design to augment, with length equal to the number of available design points.
show_progress: logical, display a progress bar during the rounding procedure (only if n_repeats > 1). Defaults to TRUE.
...: additional arguments to tune algorithm settings:
- criterion, character string specifying the optimality criterion for the rounding procedure in case design_weights has no "criterion" attribute, or to override the criterion specified in the "criterion" attribute. If specified, must be one of "D", "A", "I", "G" or "alias".

`gam_design$crit()`

Evaluate optimality criterion of a penalized approximate or rounded (exact) design. The approximate design is initialized from the design_weights field, which is the optimal (approximate) design calculated by $optimize(). To override the default approximate (optimal) design, a vector of design weights can be passed using the design_weights argument.

Usage

gam_design$crit(
  design_weights = NULL,
  criterion = c("D", "A", "I", "G", "alias"),
  lambda = TRUE,
  alpha = TRUE
)

Arguments

design_weights: (optional) vector of design weights for which to evaluate optimality criterion.
criterion: character string specifying the optimality criterion to use, see also ‘Details’. The following choices are supported:
"D" D-optimality (default).
"A" A-optimality.
"I" I-optimality.
"G" G-optimality.
"alias" alias-optimality.
lambda: logical, whether to include the ridge penalty in the calculated optimality criterion (see ‘Optimality criteria’). Only relevant if $lambda is positive.
alpha: logical, whether to include cost penalties in the calculated optimality criterion (see ‘Optimality criteria’). Only relevant if cost penalties are present in the $cost field.


`"D"` D-optimality (default).
`"A"` A-optimality.
`"I"` I-optimality.
`"G"` G-optimality.
`"alias"` alias-optimality.

Returns

numeric criterion value

`gam_design$vcov()`

Calculates the variance-covariance matrix of the model parameters based on an approximate or rounded (exact) design $(w_1, \ldots, w_n)$ according to:

\Sigma \ = \ (X^T U^{1/2} \Omega U^{1/2} X)^{-1}

where $\Omega$ is $\mathrm{diag}(\omega_1, \ldots, \omega_n)$ with $\omega_i = w_i v_i$ , and $U^{1/2}$ is $\mathrm{diag}(\sqrt{u(x_1, \theta)}, \ldots, \sqrt{u(x_n, \theta)})$ , with GLM weights $u(x_i, \theta) = \frac{1}{V(\mu_i)[g'(\mu_i)]^2}$ evaluated at the nominal parameter values. The design weights are initialized from the design_weights field, which is the optimal (approximate) design calculated by $optimize(). Use the m argument to evaluate the covariance matrix for a specific number of design points, in which case the sum of the design weights is scaled to m. To override the default approximate (optimal) design weights, a vector of design weights can be passed using the design_weights argument. The vector must be the same length as the number of design points in data. If the vector of design weights is already an integer design, leave m unspecified, in which case the design is not rescaled.

Usage

gam_design$vcov(m, design_weights = NULL, sigma2 = 1)

Arguments

m: (optional) integer number of design points at which to evaluate the estimated covariance matrix. If missing, the design weights remain unscaled, which for the default approximate (optimal) design is equivalent to m = 1.
design_weights: (optional) vector of design weights for which to evaluate the estimated covariance matrix.
sigma2: numeric error variance used to scale the covariance matrix. Defaults to 1, returning the unscaled covariance matrix.

`gam_design$sep()`

Calculates the standard error of prediction or the square root prediction variance across the design grid based on an approximate or rounded (exact) design $(w_1, \ldots, w_n)$ according to:

\mathrm{SEP}(x_i) \ = \ \sigma \sqrt{x_i^T \Sigma x_i}, \quad \mathrm{for } i = 1,\ldots,n

where $\Sigma$ is the covariance matrix as returned by $vcov(). If type = "response", the standard errors of prediction are transformed to the response scale by an application of the Delta method:

\mathrm{SEP}_{y}(x_i) \ = \ \mathrm{SEP}(x_i) \Big| \frac{\partial \mu}{\partial \eta }\Big|

Use the newdata argument to evaluate the standard error of prediction at design points not in the candidate design grid defined by $data. The design weights are initialized from the design_weights field, which is the optimal (approximate) design calculated by $optimize(). To override the default approximate (optimal) design weights, a vector of design weights can be passed using the design_weights argument. The vector must be the same length as the number of design points in data.

Usage

gam_design$sep(
  m,
  newdata,
  design_weights = NULL,
  type = c("link", "response"),
  sigma2 = 1
)

Arguments

m: (optional) integer number of design points at which to evaluate the estimated covariance matrix. If missing, the design weights remain unscaled, which for the default approximate (optimal) design is equivalent to m = 1.
newdata: (optional) data frame or matrix with design points (rows) at which to return the standard error of prediction and which can be evaluated by $formula.
design_weights: (optional) vector of design weights for which to evaluate the estimated covariance matrix.
type: the type of prediction required. The default choice "link" returns the standard error of prediction on the scale of the linear predictors, whereas "response" returns the standard error of prediction on the scale of the response variable.
sigma2: numeric error variance used to scale the standard errors. Defaults to 1, returning the standardized standard error of prediction.

`gam_design$corr()`

Calculates the correlation matrix of the columns of the design matrix based on an approximate or rounded (exact) design $(w_1, \ldots, w_n)$ according to:

R \ = \ D^{-1/2}(X_c^T U^{1/2} \Omega U^{1/2} X_c)D^{-1/2}

where $\Omega$ is $\mathrm{diag}(\omega_1, \ldots, \omega_n)$ with $\omega_i = w_i v_i$ , and $X_c$ is the design matrix $X$ centered by the (weighted) column means:

\frac{\sum_{i=1}^n \omega_i X_{ij}}{\sum_{i=1}^n \omega_i}, \quad j = 1,\ldots, p

$D^{-1/2}$ is $\mathrm{diag}(1/\sqrt{D_1},\ldots, 1/\sqrt{D_p})$ , with $D_1,\ldots,D_p$ the diagonal elements of $X_c^T U^{1/2} \Omega U^{1/2} X_c$ . If center = FALSE, the uncentered correlation matrix is returned, which corresponds to the scale-invariant (unit diagonal) version of the information matrix $M(w | X, v, \theta, F, g)$ . The design weights are initialized from the design_weights field, which is the optimal (approximate) design calculated by $optimize(). To override the default approximate (optimal) design weights, a vector of design weights can be passed using the design_weights argument. The vector must be the same length as the number of design points in data.

Usage

gam_design$corr(design_weights = NULL, center = TRUE)

Arguments

design_weights: (optional) vector of design weights for which to evaluate the information matrix.
center: logical, whether to center the design columns when evaluating the correlation matrix, defaults to TRUE.

`gam_design$clone()`

The objects of this class are cloneable with this method.

Usage

gam_design$clone(deep = FALSE)

Arguments

deep: Whether to make a deep clone.

Examples

## Ethanol dataset

## initialize model
design <- gam_design$new(
  formula = ~s(E) + C,
  data = ethanol
)
## approximate I-optimal design
design$optimize(criterion = "I")
## integer optimal design w/ 11 design points
(exact_design <- design$round(m = 11, method = "optimal", criterion = "I", seed = 1))
design$subset(replicates = exact_design)

if(FALSE) {
  ## GAM fitted to integer design
  gam_fit <- mgcv::gam(NOx ~ s(E) + C, data = ethanol[exact_design > 0, ])
  summary(gam_fit)
}

## Volcano height dataset

## pre-process data
data(volcano)
volcano <- volcano[c(TRUE, FALSE), c(TRUE, FALSE)]
data <- expand.grid(
  x1 = seq_len(dim(volcano)[1]),
  x2 = seq_len(dim(volcano)[2])
)
data[["y"]] <- c(volcano)

## initialize model
design <- gam_design$new(
  formula = ~s(x1, x2),
  data = data
)
## integer I-optimal design w/ 30 design points
exact_design <- design$round(
  m = 30, method = "optimal", criterion = "I", design_weights = FALSE
)
design$crit(exact_design, criterion = "I")

if(FALSE) {
  ## GAM fitted to integer design
  gam_fit <- mgcv::gam(y ~ s(x1, x2), data = data[exact_design > 0, ])
  summary(gam_fit)
}

gam_design

Description

Optimal design problem

Augmented optimal design

Stratified optimal design

Sparse optimal design

Optimality criteria

Seealso

Super classes

Public fields

Active bindings

Methods

Public methods

gam_design$new()

Usage

Arguments

gam_design$update()

Usage

Arguments

gam_design$stratify()

Usage

Arguments

gam_design$augment()

Usage

Arguments

gam_design$optimize()

Usage

Arguments

gam_design$sparsify()

Usage

Arguments

Details

gam_design$round()

Usage

Arguments

gam_design$crit()

Usage

Arguments

Returns

gam_design$vcov()

Usage

Arguments

gam_design$sep()

Usage

Arguments

gam_design$corr()

Usage

Arguments

gam_design$clone()

Usage

Arguments

Examples

`gam_design$new()`

`gam_design$update()`

`gam_design$stratify()`

`gam_design$augment()`

`gam_design$optimize()`

`gam_design$sparsify()`

`gam_design$round()`

`gam_design$crit()`

`gam_design$vcov()`

`gam_design$sep()`

`gam_design$corr()`

`gam_design$clone()`