lme_design

Penalized linear mixed model optimal design

Description

An R6 class for managing penalized optimal experiment designs for random intercept linear mixed models. This class encapsulates a random intercept linear mixed model, the discrete design grid, optimization of an approximate continuous design, augmentation of an existing design, and optimal rounding to exact integer designs including split-plot designs.

Optimal design problem

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

w^* \ = \ \arg \max_{w} \Phi(M(w | X, Z, \eta, v) + \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, Z, \eta, v)$ is the information matrix depending on the weights $w$ , the fixed-effects design matrix $X$ , the random-effects design matrix $Z$ , the variance component ratio $\eta = \sigma_u^2 / \sigma_e^2$ , and in the case of weighted regression the heterogeneous weights $v$ .
$\Phi$ is the chosen optimality criterion, generally a function of the information matrix.

Given a random intercept linear mixed model ( $X\beta + Zu + \epsilon$ ) defined by formula, the information matrix of the fixed-effects model parameters, conditional on the random effects covariance matrix and the error variance, is given by:

M(w | X, Z, \eta, v) \ = \ X^T\Omega^{1/2}V^{-1}\Omega^{1/2}X

where $\Omega^{1/2}$ is the diagonal matrix with $(\sqrt{w_1v_1}, \ldots, \sqrt{w_nv_n})$ on the diagonal, and $V^{-1}$ is the inverse of the block-diagonal covariance matrix $V$ . Partitioned into random effects levels $1, \ldots, S$ , where group $s$ has $n_s$ observations, $V$ consists of the level-specific blocks:

V_s = \sigma_e^2 \left( I_{n_s} + \eta \boldsymbol{1}_{n_s} \boldsymbol{1}_{n_s}^T \right)

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

using the augmented information matrix,

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

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, Z, \eta, 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}

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.

Exact split-plot design

The $round_split_plot() method rounds an approximate design to an exact (integer) split-plot design. The rounding procedure enforces fixed levels of hard-to-change variables within each individual whole plot, where the whole plots are the blocks in the random effects factors (i.e. whole plot variables). Mathematically, for an exact design $m^* = \{m_1, \ldots, m_n \}$ , the split-plot rounding procedure solves:

\begin{array}{rl} m^* \ = \ & \arg \max_{m} \Phi(M(m | X, Z, \eta, v) + \lambda I) - \alpha \cdot c^T (m / \sum m_i), \\ & \mathrm{subject\ to} \ H(m) = Z(m) H_b \mathrm{\ \ and\ \ } \boldsymbol{1}_n^T Z(m) = M \boldsymbol{1}_b^T \end{array}

$H(m)$ and $Z(m)$ are the design matrices of the hard-to-change and random effects factors subject to $m$ . $H_b$ is a matrix with $b$ rows selected from the levels of the hard-to-change variables, where $b$ is the number of whole plots. The second constraint ensures that the number of design points at each whole plot level is exactly $M$ , the whole plot size. The split-plot designs are only available in the form of exact designs, and not approximate designs (as returned by e.g. $optimize()).

Optimality criteria

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

Super class

lm_design -> lme_design

Public fields

formula: a two-part formula of the form ~ fix.eff ~ ran.eff representing a random intercept linear mixed model, where fix.eff is a linear model formula encoding the fixed effects design matrix and ran.eff is an expression encoding the levels of the random effects factor.
data: a data frame or matrix containing all design points (rows) at which formula can be evaluated.
eta: numeric variance component ratio $\sigma_u^2 / \sigma_e^2$ , where $\sigma_u^2$ is the (single) variance component of the random effects and $\sigma_e^2$ is the error variance.
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() or $augment(). 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 fixed effects design matrix obtained from all design points in data and the model formula.
Z: the random effects design matrix obtained from all design points in data and the model formula.

Methods

`lme_design$new()`

Constructor. Initializes the experimental design to optimize with a linear mixed model, discrete design grid, and fixed variance component ratio.

Usage

lme_design$new(
  formula,
  data,
  eta,
  weights = NULL,
  cost = NULL,
  contrasts = NULL,
  alpha = 1,
  lambda = 0,
  upper = 1,
  ...
)

Arguments

formula: a two-part formula of the form ~ fix.eff ~ ran.eff, where fix.eff is a linear model formula encoding the fixed effects design matrix and ran.eff is an expression encoding the random effects levels. ran.eff is evaluated as a factor, with the factor levels determining the random effects groups.
data: a data frame or matrix containing all design points (rows) at which formula can be evaluated.
eta: numeric variance component ratio $\sigma_u^2 / \sigma_e^2$ , where $\sigma_u^2$ is the (single) variance component of the random effects and $\sigma_e^2$ is the error variance.
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 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 fixed effects formula fix.eff.
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 when evaluating the fixed effects formula fix.eff.

`lme_design$update()`

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

Usage

lme_design$update(...)

Arguments

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

`lme_design$stratify()`

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

Usage

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

`lme_design$optimize()`

Optimizes the penalized approximate design based on the selected criterion by iterative reweighting of a convex approximation of the optimization problem. "alias"-optimality is included for completeness only, as the optimality criterion does not vary with the variance component ratio $\eta$ .

Usage

lme_design$optimize(
  criterion = c("D", "A", "I", "G", "alias"),
  max_iter = 1,
  tol = .Machine$double.eps^(1/4),
  scs_control = scs_control_dflt(),
  verbose = FALSE,
  return_value = TRUE,
  show_progress = 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.
max_iter: integer number of iterations sequential convex approximation. The default is 1, which returns the approximate design corresponding to a variance component ratio $\eta = 0$ .
tol: numeric tolerance used to scale the minimum improvement threshold in the backtracking line search to determine a new set of approximate design weights.
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.
show_progress: logical, display a progress bar during the rounding procedure (only if max_iter > 1). Defaults to TRUE.


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

`lme_design$augment()`

Optimizes the augmented penalized approximate design given a vector of initial design weights by iterative reweighting of a convex approximation of the optimization problem. "alias"-optimality is included for completeness only, as the optimality criterion does not vary with the variance component ratio $eta$ .

Usage

lme_design$augment(
  design_weights,
  max_iter = 1,
  gamma = 0.5,
  criterion = c("D", "A", "I", "G", "alias"),
  tol = .Machine$double.eps^(1/4),
  scs_control = scs_control_dflt(),
  verbose = FALSE,
  return_value = TRUE,
  show_progress = 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.
max_iter: integer number of iterations sequential convex approximation. The default is 1, which returns the approximate design corresponding to a variance component ratio $\eta = 0$ .
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.
tol: numeric tolerance used to scale the minimum improvement threshold in the backtracking line search to determine a new set of approximate design weights.
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.
show_progress: logical, display a progress bar during the rounding procedure (only if max_iter > 1). Defaults to TRUE.


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

`lme_design$round()`

Round penalized approximate design to an exact (integer) design. The approximate design is initialized from the design_weights field, which is the optimal (approximate) design calculated by $optimize(). For method "optimal", cost penalties and upper limit constraints present in $cost and $upper(if specified) are included in the criterion optimized by the point-exchange procedure. 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

lme_design$round(
  m,
  method = c("optimal", "efficient"),
  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).
max_iter: integer number of maximum iterations in point-exchange algorithms. Only used when method = "optimal".
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 algorithms.
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 during the rounding procedure (only if n_repeats > 1). 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".


`"optimal"` optimality criterion point-exchange algorithm (default).
`"efficient"` efficient rounding method of Pukelsheim & Reider (1992).

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.

`lme_design$round_split_plot()`

Round penalized approximate design to an exact (integer) split-plot design. This rounding procedure enforces fixed levels of the defined hard_to_change factors within each whole plot, where the whole plots are determined by the random effects factors in the model $formula. Stratification constraints imposed with $stratify() are not taken into account as these are overridden by the split-plot design constraints. The approximate design is initialized from the design_weights field, which is the optimal (approximate) design calculated by $optimize(). 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.

Usage

lme_design$round_split_plot(
  m,
  hard_to_change,
  max_iter = 100,
  n_repeats = 10,
  replicates = TRUE,
  tol = .Machine$double.eps^(1/3),
  seed = NULL,
  design_weights = NULL,
  augment_design = NULL,
  show_progress = TRUE,
  ...
)

Arguments

m: integer whole plot size(s). If a single integer, the number of design points per level of the whole plot factor(s) is set equal to m, i.e. the total number of design points in the split-plot design is m times the number of whole plot levels. If a vector, the length of m must be the same as the number of whole plot levels, with the elements of m specifying the sizes of the individual whole plots.
hard_to_change: hard-to-change factor variable of the same length as the number of rows in data. Also accepts a one-sided formula made up of column names present in data, in which case the factor levels are defined by the interaction between terms in the model.frame obtained from the formula.
max_iter: integer number of maximum iterations in point-exchange algorithms.
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 a lower threshold to continue iterations in point-exchange algorithms.
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.
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 procedure. 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", or "G".

Details

Parallel evaluation

`lme_design$crit()`

Evaluate optimality criterion of a penalized approximate or rounded (exact) design. If cost penalties are present in the $cost field, the penalized optimality criterion as described in the ‘Optimality criteria’ section of lm_design are returned. 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

lme_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

`lme_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\Omega^{1/2}V^{-1}\Omega^{1/2}X)^{-1}

Here $\Omega^{1/2} = \mathrm{diag}(\sqrt{w_1v_1}, \ldots, \sqrt{w_nv_n})$ and $V^{-1}$ is the inverse of the block-diagonal covariance matrix $V = \mathrm{diag}[V_1,\ldots,V_S]$ , consisting of the level-specific blocks $V_s = \sigma_e^2 \left( I_{n_s} + \eta \boldsymbol{1}_{n_s} \boldsymbol{1}_{n_s}^T \right)$ , with $n_s$ is the number of observations for random effects level $s$ . 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

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

`lme_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(). 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

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

`lme_design$corr()`

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

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

where $\Omega = \mathrm{diag}(\sqrt{w_1v_1}, \ldots, \sqrt{w_nv_n})$ and $V^{-1}$ is the inverse of the block-diagonal covariance matrix $V = \mathrm{diag}[V_1,\ldots,V_S]$ . $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 \Omega^{1/2} V^{-1} \Omega^{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)$ . 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

lme_design$corr(design_weights = NULL, center = TRUE, unscaled = FALSE)

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.
unscaled: logical value, whether to include $V^{-1}$ in the calculation of $R$ . If TRUE, the ordinary correlation matrix (without random effects) $R \ = \ D^{-1/2}(X_c^T \Omega X_c)D^{-1/2}$ is returned.

`lme_design$clone()`

The objects of this class are cloneable with this method.

Usage

lme_design$clone(deep = FALSE)

Arguments

deep: Whether to make a deep clone.

Examples

## Getting started

## unbalanced design grid
data <- data.frame(
  Treatment = c(-1, -1, 0, 0, 1, 1, -1, 0, 0, 1, -1, 1),
  Block = factor(rep(c("B1", "B2", "B3", "B4", "B5", "B6"), each = 2))
)
## fixed effects quadratic model in Treatment ~ random intercepts by Block
design <- lme_design$new(
  formula = ~ Treatment + I(Treatment^2) ~ Block,
  eta = 10,       ## variance component ratio
  data = data
)

## approximate D-optimal design
design$optimize(criterion = "D", max_iter = 10)
## integer design w/ 18 design points
(exact_design <- design$round(m = 6, method = "optimal", seed = 1))
design$subset(replicates = exact_design)

## Split-plot design

## Case study from [Otava & Mylona (2025), Section 4.1]

data <- expand.grid(
  A = -1:1,
  B = -1:1,
  C = -1:1,
  D = -1:1,
  week = factor(1:7, levels = 1:7)
)

design <- lme_design$new(
  formula = ~ (A + B + C + D)^2 + I(A^2) + I(B^2) + I(C^2) + I(D^2) ~ week,
  eta = 1,
  data = data
)
## initial D-optimal approximate design
design$optimize(criterion = "D", max_iter = 5, return_value = FALSE)
if(FALSE) {
  ## parallel evaluation
  mirai::daemons(n = parallel::detectCores() / 2, seed = 1)
  ## D-optimal split-plot design w/ 3 points per whole plot
  exact_design <- design$round_split_plot(
   m = 3,
   hard_to_change = ~A,
   n_repeats = 100
   )
  mirai::daemons(0)
}

## Penalized optimal design

design$update(
  cost = ~A,         ## increase cost of high levels
  alpha = 1          ## penalty scaling term
)
## initial D-optimal approximate design
design$optimize(criterion = "D", max_iter = 5, return_value = FALSE)
if(FALSE) {
  mirai::daemons(n = parallel::detectCores() / 2, seed = 1)
  exact_design <- design$round_split_plot(
   m = 3,
   hard_to_change = ~A,
   n_repeats = 100
   )
  mirai::daemons(0)
}