A Fixed, Deterministic Single-Allocation Optimal Design

A fixed-sample-size DesignFixed that computes exactly one allocation $w^*$ – the minimizer of a chosen covariate-imbalance or information objective over all allocations with $n_T = \mathrm{round}(n \cdot \mathrm{prob}_T)$ treated subjects – by numerical optimization, rather than drawing from a restricted-randomization distribution the way DesignFixedGreedy/ DesignFixedGreedyDOptimal do.

The objective family and its argument mapping. Write $Z_0 = [1\ X]$, $P = Z_0 (Z_0^\top Z_0)^{-1} Z_0^\top$, and $s(w) = n_T - w^\top P w$. The objectives map to constructor arguments and solved forms as follows:

"D" (determinant / $D_M$, $D_s$) and "A" with interest = "treatment": maximize $s(w)$, solved as the binary quadratic program $\min_w w^\top P w$ – identical criteria and interest/prior_precision semantics to DesignFixedGreedyDOptimal (the same shared construction machinery is used, so the two classes optimize literally the same matrices).
"A" with interest = "all" or a covariate subset ($A$, $A_s$): minimize $(w^\top H w + 1)/s(w)$ with the sibling class's $H$/$H_S$, solved exactly via Dinkelbach's algorithm (Dinkelbach 1967) over product-linearized MILP subproblems.
Bayesian $D_B$/$A_B$: prior_precision = a scalar $\tau$ (covariates only; intercept and treatment unpenalized) or a full matrix $R_0$, replacing $(Z_0^\top Z_0)^{-1}$ with the ridge-regularized inverse – identical to the sibling class.
"mahal_dist" / "abs_sum_diff": DesignFixedGreedy's covariate-imbalance criteria, definitionally identical (column-centered $X$, the same standardization and singular-covariance fallback), translated to exactly solvable forms: the Mahalanobis criterion is the pure quadratic $w^\top Q w$ with $Q = 4 X \Sigma^{-1} X^\top / n^2$, and the absolute-sum criterion is an l1 objective solved by the standard linear MILP.
"custom": a user-compiled black box under the user_compiled_fns.h calling convention (double f(const Eigen::MatrixXd& X, const Eigen::VectorXd& w), minimized), supplied via custom_objective; always solved by the annealing path (no structure to linearize).

Solvers and certificates. solver = "auto" (default) uses the exact "ompr" MILP path (optimum_certificate = "global", a certified global optimum) wherever tractable – always for "abs_sum_diff" (pure linear MILP); up to solver_args$linearization_max_n (default 20, set by a GLPK benchmark: the product linearization adds $n(n-1)/2$ auxiliaries and branch-and-bound cost climbs steeply past $n \approx 20$) for the quadratic and Dinkelbach criteria – and the native simulated-annealing solver beyond it, or always for "custom". The annealing solver is a formal method, not a heuristic: Metropolis acceptance over treated/control swaps with a configurable cooling schedule, for which Hajek (1988) proves convergence in probability to the global optimum under a slow-enough (logarithmic) schedule; the practical geometric schedule used by default is asymptotically motivated only, so its certificate is always "annealing_converged", never "global". solver = "ompr"/"annealing" force a path. Commercial backends extend the exact range via solver_args$roi_solver; see that parameter's wiring guides.

Inference. There is no usable randomization distribution conditional on the observed data (given $X$ there is exactly one $w^*$ up to the mirror coin), so permutation-style randomization tests/CIs are unavailable (supports_randomization_draw() is FALSE); the bootstrap randomization test IS available (the mechanism – "optimize this dataset" – is replayed on each resampled covariate matrix), as is all model-based and plain-resampling inference.

The mirror coin. At prob_T = 0.5, whenever the mirror $1 - w^*$ is a verified co-optimum (checked numerically by evaluating the objective, never by a symmetry table), a fair seeded coin picks between $w^*$ and its mirror (mirror_coin = TRUE, the default). This restores exact treated/control label symmetry – and estimator unbiasedness – at zero cost to balance. A mirror that evaluates strictly better than the solver's answer raises an error (it would be a solver bug).

References

Dinkelbach, W. (1967). On nonlinear fractional programming. Management Science 13(7):492-498, for the exact A-optimality reduction. Hajek, B. (1988). Cooling schedules for optimal annealing. Mathematics of Operations Research 13(2):311-329, for the annealing solver's formal convergence property. Atkinson, A. C., Donev, A. N., and Tobias, R. D. (2007). Optimum Experimental Designs, with SAS. Oxford University Press, for the D-/A-optimality criteria.

Super classes

Design -> DesignFixed -> DesignFixedOptimal

Methods

Public methods

DesignFixedOptimal$new()
DesignFixedOptimal$get_objective()
DesignFixedOptimal$get_interest()
DesignFixedOptimal$get_prior_precision()
DesignFixedOptimal$get_solver()
DesignFixedOptimal$get_mirror_coin()
DesignFixedOptimal$get_optimization_diagnostics()
DesignFixedOptimal$supports_randomization_draw()
DesignFixedOptimal$prepare_for_resampling_replay()
DesignFixedOptimal$clone()

+ inherited public methods from DesignFixed

+ inherited public methods from Design

`DesignFixedOptimal$new()`

Initialize a deterministic single-allocation optimal fixed experimental design. The optimization itself does not run until assign_w_to_all_subjects() (or draw_ws_according_to_design(r = 1)) is called.

Usage

DesignFixedOptimal$new(
  response_type,
  prob_T = 0.5,
  objective = "D",
  interest = "treatment",
  prior_precision = NULL,
  standardize_covariates = TRUE,
  custom_objective = NULL,
  solver = "auto",
  solver_args = list(),
  mirror_coin = TRUE,
  include_is_missing_as_a_new_feature = TRUE,
  n = NULL,
  verbose = FALSE,
  missingness_method = "impute",
  design_formula = ~.,
  seed = NULL
)

Arguments

response_type

The data type of response values.

prob_T

Probability of treatment assignment, in $(0, 1)$; the solve fixes the treated count at $\mathrm{round}(n \cdot prob_T)$.

objective

"D" (default), "A", "mahal_dist", "abs_sum_diff", or "custom"; see the class documentation.

interest

For objective = "D"/"A" only: "treatment" (default), "all", a one-sided formula, a formula string, or model-matrix column names – identical semantics to DesignFixedGreedyDOptimal.

prior_precision

For objective = "D"/"A" only: NULL (default), a positive scalar $\tau$, or a symmetric prior-precision matrix $R_0$ – identical semantics to DesignFixedGreedyDOptimal.

standardize_covariates

If TRUE (default) and prior_precision is a scalar, covariates are standardized before the penalized criterion matrices are built. ("mahal_dist"/ "abs_sum_diff" standardize internally per their definitions regardless.)

custom_objective

Required iff objective = "custom" (and forbidden otherwise): an RcppXPtrUtils::cppXPtr() external pointer or a C++ source string under the user_compiled_fns.h calling convention (double f(const Eigen::MatrixXd& X, const Eigen::VectorXd& w); X is the design's model matrix, w a candidate 0/1 allocation; the returned value is minimized). A source string is compiled through the same cppXPtr mechanism and retained so parallel workers can recompile locally. A plain R function is not accepted, and cannot be: the annealing solver evaluates the objective once per candidate swap – typically thousands of times per chain, times n_chains, and again per BRT replicate – and an R-call round-trip on every one of those evaluations is orders of magnitude too slow to be usable, not merely slower. Since "custom" is always solved by annealing (never the MILP path), there is no lower-frequency code path where an R closure would be merely inconvenient; the restriction is a hard performance requirement of this objective's only execution path. See the class examples for a worked cppXPtr() construction. Save/ reload note: a raw cppXPtr() object does not survive saveRDS()/readRDS() (see Design's "Saving and loading" section); pass a C++ source string instead if this design needs to be reloadable.

solver

"auto" (default), "ompr", or "annealing".

solver_args

A named list of solver tuning arguments. Supported: roi_solver ("glpk"/"gurobi"/"cplex" – a closed set; arbitrary ROI plugin names are rejected), linearization_max_n, max_dinkelbach_iter, n_chains, max_iter, initial_temp, cooling_rate, and (consumed by the BRT replicate path) brt_max_iter, brt_n_chains, brt_solver.

Wiring up Gurobi (roi_solver = "gurobi"): (1) obtain a Gurobi license (free academic licenses are available) and install the Gurobi Optimizer itself – this sets up GUROBI_HOME and the license file, entirely outside this package's control; (2) install Gurobi's own R package, which is not on CRAN – it ships inside the Gurobi installation: R CMD INSTALL "$GUROBI_HOME/R/gurobi_<version>_R_<Rmajor.minor>.tar.gz" (exact filename depends on your Gurobi version and platform); (3) install.packages("ROI.plugin.gurobi") from CRAN; (4) verify "gurobi" %in% ROI::ROI_registered_solvers() after loading the plugin; (5) pass solver_args = list(roi_solver = "gurobi").

Wiring up CPLEX (roi_solver = "cplex"): (1) obtain an IBM CPLEX license (free academic licenses are available) and install IBM ILOG CPLEX Optimization Studio; (2) install Rcplex (CRAN) – unlike the Gurobi bridge, it compiles from source against your local CPLEX SDK and must be pointed at your CPLEX version's include/lib directories at install time; follow Rcplex's own INSTALL instructions for your CPLEX version rather than a fixed command, since the flags change across CPLEX releases; (3) install.packages("ROI.plugin.cplex") from CRAN; (4) verify "cplex" %in% ROI::ROI_registered_solvers(); (5) pass solver_args = list(roi_solver = "cplex").

ROI.plugin.gurobi/ROI.plugin.cplex/Rcplex are deliberately never listed in this package's Suggests: declaring them would misrepresent the dependency as something install.packages() could satisfy, when the vendor installation/license underneath cannot be. Availability is checked lazily at solve time; if the plugin loads but the solve fails, the likely cause is a missing vendor installation or license.

mirror_coin

If TRUE (default), flip a fair seeded coin between $w^*$ and a verified co-optimal mirror $1 - w^*$ after every solve (only possible at prob_T = 0.5); see the class documentation.

include_is_missing_as_a_new_feature

Flag for missingness indicators.

n

Sample size (if fixed).

verbose

Flag for verbosity.

missingness_method

How to handle missing values in covariates.

design_formula

A formula object.

seed

DesignFixedOptimal$clone(deep = FALSE)

Arguments

deep: Whether to make a deep clone.

Examples

if (FALSE) { # \dontrun{
des = DesignFixedOptimal$new(n = 14, response_type = 'continuous', objective = "mahal_dist")
des$add_all_subjects_to_experiment(data.frame(x1 = rnorm(14)))
des$assign_w_to_all_subjects()
des$get_optimization_diagnostics()

# A custom compiled objective (the user_compiled_fns.h calling convention),
# built with RcppXPtrUtils::cppXPtr() -- here, squared imbalance of the
# centered covariate sums:
fobj = RcppXPtrUtils::cppXPtr(
  "double f(const Eigen::MatrixXd& X, const Eigen::VectorXd& w) {
    Eigen::RowVectorXd mu = X.colwise().mean();
    Eigen::MatrixXd Xc = X.rowwise() - mu;
    Eigen::VectorXd s = 2.0 * w - Eigen::VectorXd::Ones(X.rows());
    return (Xc.transpose() * s).squaredNorm();
  }", depends = "RcppEigen")
des2 = DesignFixedOptimal$new(n = 14, response_type = 'continuous',
  objective = "custom", custom_objective = fobj)
des2$add_all_subjects_to_experiment(data.frame(x1 = rnorm(14)))
des2$assign_w_to_all_subjects()
} # }