Julie Zhou (University of Victoria)
Jason Hou-Liu (University of Waterloo)
Chi-Kuang Yeh (Georgia State University)
Julie Zhou (University of Victoria)
Jason Hou-Liu (University of Waterloo)
April 04, 2026
This is a package to compute the optimal regression design under the second-order Least squares estimator
SLSEdesign is now available on CRAN. Hence you may install it by typing
install.packages("SLSEdesign")
or you may download the develop version by typing
# devtools::install_github("chikuang/SLSEdesign")
install.packages("SLSEdesign", repos = c('https://chikuang.r-universe.dev', 'https://cloud.r-project.org')) # or pak::pkg_install("chikuang/SLSEdesign")
library(SLSEdesign)
Consider a general regression model, $$y_i=\eta(\mathbf{x}_i, \mathbf{\theta})+ \epsilon_i, \quad i=1, \ldots, n,$$ where $y_i$ is the $i$-th observation of a response variable $y$ at design point $\mathbf{x}_i \in S \subset \mathbb{R}^p$, $S$ is a design space, $\mathbf{\theta} \in \mathbb{R}^q$ is the unknown regression parameter vector, response function $\eta(\mathbf{x}_i, \mathbf{\theta})$ can be a linear or nonlinear function of $\mathbf{\theta}$, and the errors $\epsilon_i$ are assumed to be uncorrelated with mean zero and finite variance $\sigma^2$.
Let $\hat{\mathbf{\theta}}$ be an estimator of $\mathbf{\theta}$, such as the least squares estimator. Various optimal designs are defined by minimizing $\phi \left( \mathbb{c}ov(\hat{\mathbf{\theta}}) \right)$ over the design points $\mathbf{x}_1, \ldots, \mathbf{x}_n$, where function $\phi(\cdot)$ can be determinant, trace, or other scalar functions. The resulting designs are called optimal exact designs (OEDs), which depend on the response function $\eta(\cdot,\cdot)$, the design space $S$, the estimator $\hat{\mathbf{\theta}}$, the scalar function $\phi(\cdot)$, and the number of points $n$.
Second order least-squares estimator is defined as
$$ (\mathbf{\hat{\theta}}^\top,\hat{\sigma}^2)^\top:=\underset{\mathbf{\theta},\sigma^2}{\mathrm{argmin}}\sum_{i=1}^n \begin{pmatrix} y_i-\eta(\mathbf{x}_i;\mathbf{\theta}) & y_i^2-\eta^2(\mathbf{x}_i;\mathbf{\theta})-\sigma^2 \end{pmatrix} W(\mathbf{x}_i) \begin{pmatrix} y_i-\eta(\mathbf{x}_i;\mathbf{\theta})\ y_i^2-\eta^2(\mathbf{x}_i;\mathbf{\theta})-\sigma^2 \end{pmatrix}. $$
Note that $W(\mathbf{x}_i)$ is a $2\times 2$ non-negative
semi-definite matrix which may or may not depend on $\mathbf{x}_i$ (Wang
and Leblanc, 2008). It is clear that SLSE is a natural extension of the
OLSE which is defined based on the first-order difference function
(i.e. $y_i-\mathbb{E}[y_i]=y_i-\eta(\mathbf{x}_i;\mathbf{\theta})$).
On the other hand, SLSE is defined using not only the first-order
difference function, but also second-order difference function
(i.e. $y_i^2-\mathbb{E}[y_i^2]=y_i^2-(\eta^2(\mathbf{x}_i;\mathbf{\theta})+\sigma^2))$.
One might think about the downsides of SLSE after discussing its
advantages over OLSE. SLSE does have its disadvantages. It is not a
linear estimator, and there is no closed-form solution. It requires more
computational resources compared to OLSE due to its nonlinearity.
However, numerical results can be easily computed for SLSE nowadays. As
a result, SLSE is a powerful alternative estimator to be considered in
research studies and real-life applications.
In particular, if we set the skewness parameter $t$ to be zero, the resulting optimal designs under SLSE and OLSE will be the same!
$$ y_i = \beta_0 + \beta_1 x_i + \beta_2 x_i^2 + \beta_3 x_i^3 +\varepsilon_i $$
A partial derivative of the mean function is required:
poly3 <- function(xi,theta){
matrix(c(1, xi, xi^2, xi^3), ncol = 1)
}
We first calculate the D-optimal design when the skewness parameter t
is set to be zero. The resulting D-optimal design should be the same as
the optimal design under the ordinary least-squares estimator.
my_design <- Dopt(N = 31, u = seq(-1, 1, length.out = 31),
tt = 0, FUN = poly3, theta = rep(1, 4))
my_design$design
# location weight
# 1 -1.0 0.2615264
# 10 -0.4 0.2373288
# 22 0.4 0.2373288
# 31 1.0 0.2615264
my_design$val
# 5.133616
Now we look at the situation when the skewness parameter t is in the
interval (0, 1], for instance, $0.7$.
my_design <- Dopt(N = 31, u = seq(-1, 1, length.out = 31),
tt = 0.7, FUN = poly3, theta = rep(1, 4))
my_design$design
# location weight
# 1 -1.0 0.2714088
# 10 -0.4 0.2287621
# 22 0.4 0.2287621
# 31 1.0 0.2714088
my_design$val
# 6.27293
Add equivalence theorem plot for D-optimal design:
design = data.frame(location = c(-1, -0.447, 0.447, 1),
weight = rep(0.25, 4))
u = seq(-1, 1, length.out = 201)
plot_dispersion(u, design, tt = 0, FUN = poly3,
theta = rep(0, 4), criterion = "D")
$$ y_i = \beta_1 x_i + \beta_2 x_i^2 + \beta_3 x_i^3 +\varepsilon_i $$
In the last example, the support points did not change as t increases.
However, it is not always the case, and. the optimal design may be
depending on t.
poly3_no_intercept <- function(xi, theta){
matrix(c(xi, xi^2, xi^3), ncol = 1)
}
my_design <- Dopt(N = 31, u = seq(-1, 1, length.out = 31),
tt = 0, FUN = poly3_no_intercept, theta = rep(1, 3))
my_design$design
# location weight
# 1 -1.0 0.3275005
# 7 -0.6 0.1565560
# 25 0.6 0.1565560
# 31 1.0 0.3275005
my_design$val
# 3.651524
my_design <- Dopt(N = 31, u = seq(-1, 1, length.out = 31),
tt = 0.9, FUN = poly3_no_intercept, theta = rep(1, 3))
my_design$design
# location weight
# 1 -1.0 0.2888423
# 10 -0.4 0.2096781
# 22 0.4 0.2096781
# 31 1.0 0.2888423
my_design$val
# 4.892601
CVXR::exclude_solvers("MOSEK")
peleg_grad <- function(xi, theta){
th1 <- theta[1]
th2 <- theta[2]
den <- (th1 + th2 * xi)^2
matrix(c(
-xi / den,
-xi^2 / den
), ncol = 1)
}
u <- seq(0, 180, length.out = 1001)
theta0 <- c(0.5, 0.05)
res0 <- SLSEdesign::Dopt(
N = length(u),
u = u,
tt = 0,
FUN = peleg_grad,
theta = theta0
)
res0$design |> round(3)
## location weight
## 51 9 0.5
## 1001 180 0.5
res0$val
## [,1]
## [1,] 14.87739