glm_design

Penalized generalized linear model locally optimal design

Description

An R6 class for managing penalized locally optimal experiment designs for generalized linear models. This class encapsulates the generalized linear 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 on the weights $w$ and $v$ , the design matrix $X$ , 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.

In standard canonical form, the probability density (or mass) function of the exponential family $F$ is:

f(y\ |\ \zeta, \phi) = \exp\left( \frac{y\zeta - b(\zeta)}{a(\phi)} + c(y, \phi) \right)

The cumulant function $b(\zeta)$ admits $\mu = E[Y] = b'(\zeta)$ and $V(\mu) = b''(\zeta)$ , the variance function expressed as a function of the mean. Given the link function $g$ , the systematic component is $g(\mu) = \eta = x^T \theta$ . 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}$ causing the information matrix to depend on the parameter values $\theta$ , unlike the information matrix $M(w | X, v)$ for a standard linear model.

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.

References

Lewis, S. L., Montgomery D. C., and Myers R. H. 2001. Examples of Designed Experiments with Nonnormal Responses. Journal of Quality Technology 3 (33): 265–78.

Super class

lm_design -> glm_design

Public fields

formula: a linear 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.
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 from all design points in data and the model formula.

Methods

`glm_design$new()`

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

Usage

glm_design$new(
  formula,
  data,
  theta,
  family = gaussian(),
  weights = NULL,
  cost = 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.
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.
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 model.matrix.

`glm_design$update()`

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

Usage

glm_design$update(...)

Arguments

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

`glm_design$stratify()`

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

Usage

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

`glm_design$augment()`

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

Usage

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

`glm_design$optimize()`

Optimizes the penalized approximate design based on the selected criterion.

Usage

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

`glm_design$sparsify()`

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

Usage

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

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

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

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

glm_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

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

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

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

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

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

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

`glm_design$clone()`

The objects of this class are cloneable with this method.

Usage

glm_design$clone(deep = FALSE)

Arguments

deep: Whether to make a deep clone.

Examples

## Non-conforming tiles experiment
## Example 3 from [Lewis et al. (2001)]

## design grid
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),
  x7 = c(-1, 1)
)
## initialize model, design grid, family and nominal parameters
design <- glm_design$new(
  formula = ~x1 + x2 + x3 + x4 + x5 + x6 + x7,
  data = data,
  family = binomial(link = "logit"),
  theta = c(-1.419, -0.578, 0.109, 0.266, -0.262, 0.434, -0.633, 0.425)
)
## approximate D-optimal design
design$optimize(criterion = "D")
## integer design w/ 8 design points
(exact_design <- design$round(m = 8, method = "optimal", n_repeats = 25, seed = 1))
design$subset(replicates = exact_design)

## Defect semiconductor wafers
## Example 5 from [Lewis et al. (2001)]

## design grid
data <- rbind(
  expand.grid(x1 = -1:1, x2 = -1:1, x3 = -1:1),
  data.frame(
    x1 = c(-1.732, 1.732, 0, 0, 0, 0),
    x2 = c(0, 0, -1.732, 1.732, 0, 0),
    x3 = c(0, 0, 0, 0, -1.732, 1.732)
  )
)
## initialize model, design grid, family and nominal parameters
design <- glm_design$new(
  formula = ~(x1 + x2 + x3)^2 + I(x1^2) + I(x2^2) + I(x3^2),
  data = data,
  family = poisson(link = "log"),
  theta = c(1.9293, 0.6103, 0.7028, -0.3531, 0, 0.2471, 0, -0.5209, 0, 0)
)
## approximate D-optimal design
design$optimize(criterion = "D")
## integer design w/ 18 design points
(exact_design <- design$round(m = 18, method = "optimal", seed = 1))
design$subset(replicates = exact_design)

glm_design

Description

Optimal design problem

Augmented optimal design

Stratified optimal design

Sparse optimal design

Optimality criteria

References

Seealso

Super class

Public fields

Active bindings

Methods

Public methods

glm_design$new()

Usage

Arguments

glm_design$update()

Usage

Arguments

glm_design$stratify()

Usage

Arguments

glm_design$augment()

Usage

Arguments

glm_design$optimize()

Usage

Arguments

glm_design$sparsify()

Usage

Arguments

Details

glm_design$round()

Usage

Arguments

glm_design$crit()

Usage

Arguments

Returns

glm_design$vcov()

Usage

Arguments

glm_design$sep()

Usage

Arguments

glm_design$corr()

Usage

Arguments

glm_design$clone()

Usage

Arguments

Examples

`glm_design$new()`

`glm_design$update()`

`glm_design$stratify()`

`glm_design$augment()`

`glm_design$optimize()`

`glm_design$sparsify()`

`glm_design$round()`

`glm_design$crit()`

`glm_design$vcov()`

`glm_design$sep()`

`glm_design$corr()`

`glm_design$clone()`