nlm_design

Penalized nonlinear locally optimal design

Description

An R6 class for managing penalized locally optimal experiment designs for nonlinear models. This class encapsulates the nonlinear 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) + \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)$ is the information matrix depending on the weights $w$ and $v$ , the design matrix $X$ and the nominal parameter values $\theta$ .
$\lambda I$ is a diagonal matrix to regularize the information matrix with ridge penalty parameter $\lambda$ .
$\Phi$ is the chosen optimality criterion, generally a function of the information matrix.

Given the nonlinear model f and nominal parameter values theta, the information matrix (conditional on the error variance) is given by:

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

where $F(x_i, \theta)$ is the gradient vector at $f(x_i, \theta)$ , i.e. the $i$ -th row of the Jacobian matrix.

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) + \lambda I) - \alpha \cdot c^T w

using the augmented information matrix,

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

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) + \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 class

lm_design -> nlm_design

Public fields

formula: a nonlinear model defined as a one-sided formula.
theta: a named vector or list of nominal parameter values at which the information matrix $M(w, \theta)$ is evaluated.
data: a data frame or matrix containing all design points (rows) at which formula can be evaluated.
cost: an optional numeric vector containing costs associated to each design point (row) in data.
weights: an optional numeric vector containing fixed user weights 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 for the alias optimality criterion.
upper: the design weight upper limit. The default is 1, corresponding to no upper limit constraint on the design weights.

Active bindings

X: the Jacobian matrix obtained from all design points in data and the model formula.

Methods

`nlm_design$new()`

Constructor. Initializes the experimental design to optimize with a nonlinear model formula and discrete design grid.

Usage

nlm_design$new(
  formula,
  data,
  theta,
  weights = NULL,
  cost = NULL,
  alpha = 1,
  lambda = 0,
  upper = 1
)

Arguments

formula: a nonlinear model defined as a one-sided formula.
data: a data frame or matrix containing all design points (rows) at which formula can be evaluated.
theta: a named vector or list of nominal parameter values at which the information matrix $M(w, \theta)$ is evaluated.
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.
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.

`nlm_design$update()`

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

Usage

nlm_design$update(...)

Arguments

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

`nlm_design$stratify()`

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

Usage

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

`nlm_design$augment()`

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

Usage

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

`nlm_design$optimize()`

Optimizes the penalized approximate design based on the selected criterion.

Usage

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

`nlm_design$sparsify()`

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

Usage

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

`nlm_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

nlm_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 "A", "I", "D", "G" or "alias".

`nlm_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

nlm_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

`nlm_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 \ = \ (F(X, \theta)^T \Omega F(X, \theta))^{-1}

where $\Omega$ is $\mathrm{diag}(\omega_1, \ldots, \omega_n)$ with $\omega_i = w_i v_i$ , and $F(X, \theta)$ is the Jacobian matrix 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

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

`nlm_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{F(x_i, \theta)^T \Sigma F(x_i, \theta)}, \quad \mathrm{for } i = 1,\ldots,n

where $\Sigma$ is the covariance matrix as returned by $vcov(). 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

nlm_design$sep(m, newdata, 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.
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.
sigma2: numeric error variance used to scale the standard errors. Defaults to 1, returning the standardized standard error of prediction.

`nlm_design$corr()`

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

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

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

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

Furthermore, $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 $J_c^T \Omega J_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)$ . 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

nlm_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 Jacobian columns when evaluating the correlation matrix, defaults to TRUE.

`nlm_design$clone()`

The objects of this class are cloneable with this method.

Usage

nlm_design$clone(deep = FALSE)

Arguments

deep: Whether to make a deep clone.

Examples

## Lubricant example from [Bates & Watts (1998), Appendix 1, A1.8]

## design grid
data <- expand.grid(
  x1 = seq(from = 0, to = 100, by = 10), ## temp. (C)
  x2 = seq(from = 0, to = 7, by = 0.5)   ## pressure (atm)
)
## initialize model, nominal parameters and design grid
design <- nlm_design$new(
  formula = ~b1/(b2 + x1) + b3 * x2 + b4 * x2^2 + b5 * x2^3 +
              (b6 * x2 + b7 * x2^3) * exp(-x1/(b8 + b9 * x2^2)),
  theta = c(b1 = 1054.541, b2 = 206.546, b3 = 1.46, b4 = -0.26,
            b5 = 0.023, b6 = 0.401, b7 = 0.035, b8 = 57.405, b9 = -0.477),
  data = data
)
## approximate I-optimal design
design$optimize(criterion = "I")
## integer design w/ 11 design points
(exact_design <- design$round(m = 11, method = "optimal", seed = 1))
design$subset(replicates = exact_design)

## Penalized optimal design

design$update(
  cost = ~abs(x1 - 25) / 100 + abs(x2 - 1) / 7,                  ## condition weights
  alpha = 0.1,
  upper = 0.1                                                    ## max weight per point
)
## re-optimize
design$optimize(criterion = "I")
## optimal rounding w/ penalization
design$round(m = 14, method = "optimal", seed = 1)

## Stratified optimal design

design$update(
  cost = NULL,   ## drop cost penalties
  data = subset(data, x1 %in% c(0, 20, 40, 100))
)
## partition design space by temperature conditions
## by default, equal weight proportions are assigned to each condition
design$stratify(
  strata = ~x1
)
## stratified optimal design
design$optimize(criterion = "I")
## verify stratification
tapply(design$design_weights, design$data$x1, sum)
## optimal rounding w/ stratification
design$round(m = 12, method = "optimal", seed = 1)