lm_design

Penalized linear optimal design

Description

An R6 class for managing penalized optimal experiment designs for linear models. This class encapsulates the linear model formula, 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Φ(M(w∣X,v)+λI)−α⋅cTww^* \ = \ \arg \max_{w} \Phi(M(w | X, v) + \lambda I) - \alpha \cdot c^T w

where,

  • w={w1,…,wn}w = \{w_1, \ldots, w_n \} is the vector of design weights.
  • c={c1,…,cn}c = \{c_1, \ldots, c_n \} is the cost assigned to each design weight with scaling/penalty parameter α\alpha.
  • v={v1,…,vn}v = \{v_1, \ldots, v_n \} is a vector of fixed user weights.
  • M(w∣X,v)M(w | X, v) is the information matrix depending on the weights ww, the design matrix XX created by model.matrix, and in the case of weighted regression the heterogeneous weights vv.
  • λI\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 linear model formula, the information matrix (conditional on the error variance) is:

M(w∣X) ∝ ∑i=1nwivixixiTM(w | X) \ \propto \ \sum_{i = 1}^n w_i v_i x_i x_i^T

Augmented optimal design

Given a vector of initial design weights w0w_0 and associated contribution weight γ∈[0,1]\gamma \in [0, 1], the following augmented penalized optimization problem is solved by $augment():

w∗ = arg⁡max⁡wΦ(M(w∣γ,w0,X,v)+λI)−α⋅cTww^* \ = \ \arg \max_{w} \Phi(M(w | \gamma, w_0, X, v) + \lambda I) - \alpha \cdot c^T w

using the augmented information matrix,

M(w∣γ,w0,X,v) =(1−γ)M(w0∣X,v)+γM(w∣X,v)M(w | \gamma, w_0, X, v) \ = (1 - \gamma) M(w_0 | X, v) + \gamma M(w | X, v)

Stratified optimal design

The $strata() method stratifies the design by requiring that individual strata X1,…,XSX_1, \ldots, X_Sare represented proportionally with given weights p1,…,pSp_1, \ldots, p_S, (∑ipi=1\sum_i p_i = 1). If stratification is applied, the following optimization problem is solved by $optimize():

w∗ = arg⁡max⁡wΦ(M(w∣X,v)+λI)−α⋅cTw,subject to ∑i∈Xswi=ps,for s=1,…,S\begin{array}{rl} w^* \ = \ & \arg \max_{w} \Phi(M(w | X, v) + \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 XiX_i is equal to pip_i for each i=1,…,Si = 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 using a custom penalty function P(w)P(w)to promote sparsity in the approximate design w∗w^*. Because the sparse optimization problem is not convex, a sequential convex approximation is used substituting the (concave) sparseness penalty by linearized (convex) costs on the design weights ww similar to the cost penalties above. The following options are supported:

"entropy"-penalty

Shannon entropy penalty, P(w)=−∑inwilog⁡(wi)P(w) = -\sum_i^n w_i \log(w_i), with linearized costs (iteration kk):

ci(k)=−(1+log⁡(wi(k)+ϵ))c_i^{(k)} = -(1 + \log(w_i^{(k)} + \epsilon))

"log-sum"-penalty

Log-sum penalty, P(w)=∑inlog⁡(wi+ϵ)P(w) = \sum_i^n \log(w_i + \epsilon), with linearized costs (iteration kk):

ci(k)=1wi(k)+ϵc_i^{(k)} = \frac{1}{w_i^{(k)} + \epsilon}

"power"-penalty

Power (bridge) penalty, P(w)=∑in(wi+ϵ)pP(w) = \sum_i^n (w_i + \epsilon)^p, with linearized costs (iteration kk):

ci(k)=p(wi(k)+ϵ)p−1,0<p<1c_i^{(k)} = p (w_i^{(k)} + \epsilon)^{p - 1}, \quad 0 < p < 1

Optimality criteria

The following optimality criteria Φ(⋅)\Phi(\cdot) are supported in $optimize(), $augment() and $sparsify():

"D"-optimality

Maximize the volume of the information matrix MM, or equivalently minimize the volume of the confidence ellipsoid of the parameter estimates:

Φ(M)=det⁡(M)1/p\Phi(M) = \det(M)^{1/p}

where pp is the number of columns in the design matrix XX. With non-zero ridge/cost penalties, the (penalized) optimality criterion is calculated as:

Φc(w,M)=exp⁡(log⁡(Φ(M+λI))−α⋅cTw)\Phi_c(w, M) = \exp\left(\log(\Phi(M + \lambda I)) - \alpha \cdot c^T w \right)

"A"-optimality

Maximize the inverse trace of the inverse information M−1M^{-1}, or equivalently minimize the average variance of the parameter estimates:

Φ(M)=(Tr(M−1))−1\Phi(M) = \left(\mathrm{Tr}(M^{-1})\right)^{-1}

With non-zero ridge/cost penalties, the (penalized) optimality criterion is calculated as:

Φc(w,M)=Φ(M+λI)−α⋅cTw\Phi_c(w, M) = \Phi(M + \lambda I) - \alpha \cdot c^T w

"I"-optimality

Maximize the inverse of the average prediction variance across the design space, (equivalently minimize the average prediction variance):

Φ(M)=(∑i=1nxiTM−1xi)−1\Phi(M) = \left(\sum_{i=1}^n x_i^T M^{-1} x_i\right)^{-1}

With non-zero ridge/cost penalties, the (penalized) optimality criterion is calculated as:

Φc(w,M)=Φ(M+λI)−α⋅cTw\Phi_c(w, M) = \Phi(M + \lambda I) - \alpha \cdot c^T w

"G"-optimality

Maximize the inverse of the maximum prediction variance across the design space, (equivalently minimize the maximum prediction variance):

Φ(M)=(max⁡xixiTM−1xi)−1\Phi(M) = \left(\max_{x_i} x_i^T M^{-1} x_i\right)^{-1}

With non-zero ridge/cost penalties, the (penalized) optimality criterion is calculated as:

Φc(w,M)=Φ(M+λI)−α⋅cTw\Phi_c(w, M) = \Phi(M + \lambda I) - \alpha \cdot c^T w

"alias"-optimality

Maximize one minus the average absolute cosine similarity (also uncentered correlation or normalized dot product) between columns of the design matrix XX:

Φ(w,v,X)=1−2p(p−1)∑i<j∣⟨ωxi,xj⟩∥ωxi∥∥ωxj∥∣\Phi(w, v, X) = 1 - \frac{2}{p(p - 1)} \sum_{i < j} \left| \frac{\langle \omega x_i, x_j \rangle}{\Vert \sqrt{\omega} x_i \Vert \Vert \sqrt{\omega} x_j \Vert} \right|

where ωxi\omega x_i is the elementwise-product of ω=(w1v1,…,wnvn)\omega = (w_1 v_1, \ldots, w_n v_n) and the ii-th design column xi=(x1,i,…,xn,i)x_i = (x_{1,i},\ldots,x_{n,i}). If all columns in the design matrix XX are centered, this is equivalent to minimizing the average correlation between columns in XX. With non-zero cost penalties, the (penalized) optimality criterion is calculated as:

Φc(w,v,X)=Φ(w,v,X)−α⋅cTw\Phi_c(w, v, X) = \Phi(w, v, X) - \alpha \cdot c^T w

References

Pukelsheim, F., and Rieder S. (1992) “Efficient rounding of approximate designs”, Biometrika.

Seealso

https://www.cvxgrp.org/scs/https://ergo-code.github.io/HiGHS/dev/installation/ https://mirai.r-lib.org/

Public fields

  • formula: a linear model defined as a one-sided formula.
  • 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.
  • 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.
  • weights: an optional numeric vector containing fixed user weights associated to each design point (row) in data.
  • 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 design matrix obtained from all design points in data and the model formula.
  • orthogonal: the dot product (orthogonality) constraint matrix. This matrix contains all dot product (upper limit) constraints between design columns parsed from the orthogonal formula(s) in $new(orthogonal = ...) or $update(orthogonal = ...). For design columns with no dot product constraint, the matrix entry is Inf.
  • zeros: the proportion of zeros vector. This vector contains the required proportion of zeros for each column in the design matrix parsed from the zeros formula(s) in $new(zeros = ...) or $update(zeros = ...). For design columns with no constraint, the value is NA.
  • col_means: the imposed column mean vector. This vector contains the imposed column mean constraints for each column in the design matrix parsed from the col_means formula in $new(col_means = ...) or $update(col_means = ...). For design columns with no mean constraint, the value is NA.

Methods

Public methods

lm_design$new()

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

Usage

lm_design$new(
  formula,
  data,
  orthogonal = NULL,
  weights = NULL,
  cost = NULL,
  zeros = NULL,
  col_means = NULL,
  contrasts = NULL,
  alpha = 1,
  lambda = 0,
  upper = 1,
  ...
)

Arguments

  • formula: a linear model defined as a one-sided formula.
  • data: a data frame or matrix containing all design points (rows) at which formula can be evaluated.
  • orthogonal: an optional two- or three-part formula or list of formulas imposing dot product (orthogonality) constraints between the columns of the design matrix. In the three-part formula x1 + x2 ~ y1 + y2 ~ 0, dot products are constrained between all combinations of columns in the first part of the formula (x1, x2) and the second part of the formula (y1, y2). The third part (0) specifies the upper limit on the dot product(s). If the formula only has two parts, e.g. x1 + x2 ~ y1 + y2, zero dot product constraints are imposed. Also accepts a list of formulas, in which case constraints are added in the order in which they appear in the list.
  • weights: an optional numeric vector containing fixed user weights associated to each design point (row) in data. Also accepts a one-sided formula using the 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 using the column names present in data. If NULL, the standard (non-penalized) optimal design problem is solved.
  • zeros: an optional two-part formula or list of formulas imposing constraints on the proportion of zero values in the columns of the design matrix. In the two-part formula x1 + x2 ~ 1/3, the proportion of zeros in the columns specified by the first part of the formula (x1, x2) is constrained to be equal to the (nonnegative) numeric value in the second part of the formula (1/3). Note that imposing constraints on the proportion of zeros in the optimized design is only relevant if the design grid actually contain zero entries. For a list of formulas, zero constraints are added in the order in which they appear in the list, analogous to the orthogonal argument.
  • col_means: an optional two-part formula imposing constraints on the column means of the design matrix. In the two-part formula x1 + x2 ~ 0, the means of the columns specified by the first part of the formula (x1, x2) are constrained to be equal to the numeric value in the second part of the formula (0).
  • contrasts: optional contrasts used when evaluating the linear model formula.
  • alpha: numeric value, cost scale 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 model.matrix.

lm_design$update()

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

Usage

lm_design$update(...)

Arguments

  • ...: named arguments formula, data, weigts, orthogonal, cost, zeros, col_means, alpha, lambda and/or upper.

lm_design$stratify()

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

Usage

lm_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 using 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.

lm_design$augment()

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

Usage

lm_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 column names present 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 can be obtained from the design_weights field.

lm_design$optimize()

Optimizes the penalized approximate design based on the selected criterion.

Usage

lm_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 can be obtained from the design_weights field.

lm_design$sparsify()

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

Usage

lm_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,
  recenter = FALSE,
  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 ‘Details’ for the definitions.
  • eps: positive numeric value to ensure finite weights in reweighting functions, see ‘Details’.
  • 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<10 < p < 1. Only used if reweight_fun = "power".
  • design_weights: (optional) vector of design weights used as initial approximate design.
  • recenter: logical, recenter design columns based on the current design weights after each iteration. Defaults to FALSE.
  • 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 can be obtained from the design_weights field.
  • show_progress: logical, display a progress bar during the iterative reweighting procedure. Defaults to TRUE.

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.

lm_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

lm_design$round(
  m,
  method = c("optimal", "efficient", "dist-sum", "dist-min", "orthogonal"),
  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.
    "orthogonal" minimizes aliasing subject to orthogonality constraints, see ‘Details’.
  • 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. Not used for method = "orthogonal". 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. Also used as threshold to ensure positive diagonal entries of M(w∣X)M(w | X) when method = "orthogonal".
  • seed: (optional) random seed used to generate jitter to break ties in rounding the initial exact design. Not used in the case of (mirai) asynchronous evaluation, which requires the random seed to be set via the seed argument in mirai::daemons.
  • 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 monitoring the repeats of the rounding/exchange algorithm. Defaults to TRUE.
  • ...: additional arguments to tune algorithm/solver 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".
    • highs_control, a named list with control arguments supported by the HiGHS solver, only used if method = "orthogonal". Defaults to highs_control = list(time_limit = 180, log_to_console = TRUE).
    • highs_write_mps, .mps file path. If provided, the HiGHS MILP is exported in MPS (Mathematical Programming System) format to the provided file path. Can also be TRUE, in which case a new .mps file is created in the current working directory. Only used if method = "orthogonal".

Details

Parallel evaluation

The repeated rounding/exchange algorithms can be run in parallel using the mirai asynchronous evaluation framework. The parallel background R processes are configured by the user outside of $round() with mirai::daemons. If daemon connections are available, $round() automatically runs the repeated exchange algorithms in parallel using the configured daemons. If the underlying BLAS/LAPACK libraries are already multi-threaded, this may inhibit the use of additional threads by the mirai asynchronous evaluation. In this case, it is recommended to restrict BLAS/LAPACK to a single thread with the OPENBLAS_NUM_THREADS, MKL_NUM_THREADS or OMP_NUM_THREADSenvironment variables, depending on the BLAS library that is used. In addition, note that there is a small overhead in populating the environments of the background R processes, which may result in increased runtime compared to sequential evaluation in the main R process, especially when n_repeats is small or the rounding/exchange algorithm does not take much time to complete.

HiGHS solver

The "orthogonal" rounding procedure uses HiGHS to find an exact design subject to orthogonality constraints. This requires either a local installation of the HiGHS solver or the highs package to be installed. To point to a local HiGHS installation, the environment variable HIGHS_HOME must be set to the highsroot directory when installing cvxoptdes. If installation fails, verify that highs_c_api.h exists, typically present in $HIGHS_HOME/include/highs/interfaces, which is required to compile the source code.

lm_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

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

Returns

numeric criterion value

lm_design$vcov()

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

Σ = (XTΩX)−1\Sigma \ = \ (X^T \Omega X)^{-1}

where Ω\Omega is diag(ω1,…,ωn)\mathrm{diag}(\omega_1, \ldots, \omega_n) with ωi=wivi\omega_i = w_i v_i. 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.

Usage

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

lm_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 (w1,…,wn)(w_1, \ldots, w_n) according to:

SEP(xi) = σxiTΣxi,fori=1,…,n\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(). 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

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

lm_design$corr()

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

R = D−1/2(XcTΩXc)D−1/2R \ = \ D^{-1/2}(X_c^T \Omega X_c)D^{-1/2}

where Ω\Omega is diag(ω1,…,ωn)\mathrm{diag}(\omega_1, \ldots, \omega_n) with ωi=wivi\omega_i = w_i v_i, and XcX_c is the design matrix XX centered by the (weighted) column means:

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

D−1/2D^{-1/2} is diag(1/D1,…,1/Dp)\mathrm{diag}(1/\sqrt{D_1},\ldots, 1/\sqrt{D_p}), with D1,…,DpD_1,\ldots,D_p the diagonal elements of XcTΩXcX_c^T \Omega 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)M(w | X, v). 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

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

lm_design$subset()

Subsets the design points associated to an exact design vector with integer replicates as returned by $round(). The vector must be the same length as the number of design points in data.

Usage

lm_design$subset(replicates)

Arguments

  • replicates: integer vector of design point replicates of the same length as the candidate grid in data

lm_design$clone()

The objects of this class are cloneable with this method.

Usage

lm_design$clone(deep = FALSE)

Arguments

  • deep: Whether to make a deep clone.

Examples

## Getting started

## 3 factors w/ 5, 5 and 3 levels
data <- expand.grid(
    x1 = -2:2,
    x2 = -2:2,
    x3 = -1:1
)
design <- lm_design$new(
    formula = ~(x1 + x2 + x3)^2 + I(x1^2) + I(x2^2) + I(x3^2),   ## full quadratic model
    data = data    ## design points
)
## approximate D-optimal design
design$optimize(criterion = "D")
## integer D-optimal design w/ 13 design points
(exact_design <- design$round(m = 13, method = "optimal", seed = 1))
design$subset(replicates = exact_design)
## integer design covariance and correlation matrices
design$vcov(design_weights = exact_design, sigma2 = 1)
design$corr(design_weights = exact_design)

if(FALSE) {
 ## parallel evaluation
 mirai::daemons(n = parallel::detectCores() / 2, seed = 1)
 design$round(m = 100, method = "optimal", n_repeats = 1000)
 mirai::daemons(0)
}

## Penalized optimal design

## assign cost c = 1 if x3 = 1, and c = 0 otherwise
design$update(
  cost = ~pmax(x3, 0),
  alpha = 1
)
## penalized approximate D-optimal design
design$optimize(criterion = "D")

## assign v = 1 if x3 = 1, and v = 0.1 otherwise
design$update(
 weights = ~pmax(x3, 0.1),
 cost = NULL
)
## penalized approximate D-optimal design
design$optimize(criterion = "D")
## exact design w/ 15 design points
design$round(m = 15, method = "optimal", seed = 1)

## drop weights term
design$update(weights = NULL)

## Stratified optimal design

## add stratification constraints
design$stratify(
  strata = ~x3,
  proportions = c(1/4, 1/2, 1/4)
)
## approximate D-optimal design w/ stratification
design$optimize(criterion = "D")
## verify stratification
tapply(design$design_weights, design$data$x3, sum)
## exact design w/ stratification
design$round(m = 12, method = "optimal", seed = 1)
## drop stratification
design$stratify(strata = NULL)

## Orthogonality constraints

## orthogonal main effects wrt all other columns
design$update(
  orthogonal = x1 + x2 + x3 ~ . ~ 0,
  alpha = 0.075  ## sparseness penalty
)
design$orthogonal
## approximate D-optimal design w/ orthogonality constraints
design$optimize(criterion = "D")
## sparsify support by iterative entropy reweighting
design$sparsify(max_iter = 5, criterion = "D", reweight_fun = "entropy")

if(FALSE) {
 ## integer design w/ orthogonality constraints
 (ortho_design <- design$round(m = 8, method = "orthogonal", seed = 1))
 ## verify correlations
 design$corr(ortho_design)
}

## Optimal design augmentation

## Example from [Goos & Jones (2011), Chapter 3]
## 6 factors w/ 2 levels each
data <- expand.grid(
  x1 = c(-1, 1),
  x2 = c(-1, 1),
  x3 = c(-1, 1),
  x4 = c(-1, 1),
  x5 = c(-1, 1),
  x6 = c(-1, 1)
)

## linear model w/ main effects and selected interactions
design <- lm_design$new(
  formula = ~1 + x1 + x2 + x3 + x4 + x5 + x6 + x2:x3 + x1:x6 + x2:x6 + x3:x6 + x4:x6 + x1:x4,
  data = data
)
## initial design of size n=12
init_design <- c(0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0,
 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0,
 0, 1, 0, 0, 0, 0, 1)
## approximate optimal augmented design
design$augment(design_weights = init_design, gamma = 0.4, criterion = "D")
## exact augmented design
(augment_design <- design$round(m = 8, method = "optimal", seed = 1, augment_design = init_design))
## proposed design points
design$subset(replicates = augment_design - init_design)

## augmented design w/ orthogonality constraints
design_ortho <- design$clone()
design_ortho$update(
  orthogonal = list(
    x1 + x2 + x3 + x4 + x5 + x6 ~ x1 + x2 + x3 + x4 + x5 + x6 ~ 0  ## uncorrelated main effects
  )
)
## approximate augmented design
design_ortho$augment(
  design_weights = init_design, gamma = 0.4, criterion = "D", return_value = FALSE
)

if(FALSE) {
  ## integer design w/ orthogonality constraints
  ortho_augment_design <- design_ortho$round(
    m = 8, method = "orthogonal", augment_design = init_design
  )
  ## verify correlations
  design_ortho$corr(ortho_augment_design)
}