Computes an efficient approximate experimental design under general linear constraints using the approach of second-order cone programming.

1 2 |

`F` |
The |

`b, A` |
The real vector of length |

`w0` |
The non-negative vector of length |

`crit` |
The optimality criterion. Possible values are |

`R` |
The region of summation for the IV-optimality criterion. The argument |

`kappa` |
A small non-negative perturbation parameter. |

`tab` |
A vector determining the regressor components to be printed with the resulting design.
This argument should be a subvector of |

`graph` |
A vector determining the regressor components to be plotted with the resulting design.
This argument should be a subvector of |

`t.max` |
The time limit for the computation. |

The procedure computes an efficient approximate design by converting the optimal design problem to a specific problem of second-order cone programming; see the reference for details. The advantage of this approach is the possibility to construct approximate designs under a general system of linear constraints. In particular, this function provides means of computing informative lower bounds on the efficiency of the linearly constrained exact designs computed using methods such as od.RC, od.IQP, and od.MISOCP.

The model should be non-singular in the sense that there exists an approximate design `w`

satisfying the constraints `0<=w0<=w`

and
`A%*%w<=b`

, with a non-singular information matrix, preferably with the reciprocal condition number of at least `1e-5`

. If this requirement
is not satisfied, the computation may fail, or it may produce a deficient design.

If the criterion of IV-optimality is selected, the region `R`

should be chosen such that the associated matrix `L`

(see the help page of the function od.crit) is non-singular, preferably with a reciprocal condition number of at least `1e-5`

.
If this requirement is not satisfied, the computation may fail, or it may produce a deficient design.

The perturbation parameter `kappa`

can be used to add `n*m`

iid random numbers from the uniform distribution
in `[-kappa,kappa]`

to the elements of `F`

before the optimization is executed. This can be helpful for
increasing the numerical stability of the computation or for generating a random design from the potentially large set of optimal or
nearly-optimal designs.

The performance depends on the problem and on the hardware used, but in most cases the function can compute an optimal or nearly-optimal exact design for a problem with a thousand design points within minutes of computing time. If the only constraint on the design is the standard constraint on the size, the function od.AA should be a preferred choice.

A list with the following components:

`method` |
The method used for computing the design |

`w.best` |
the best permissible design found, or |

`Phi.best` |
The value of the criterion of optimality of the design |

`status` |
The status variable of the gurobi optimization procedure; see the gurobi solver documentation for details. |

`t.act` |
The actual time taken by the computation. |

Radoslav Harman, Lenka Filova

Sagnol G, Harman R (2015): Computing exact D-optimal designs by mixed integer second-order cone programming. The Annals of Statistics, Volume 43, Number 5, pp. 2198-2224.

`od.AA, od.MISOCP, od.IQP, od.RC`

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 | ```
if(require("gurobi")){
# Suppose that we model the mean values of the observations of
# a circadian rhythm by a third-degree trigonometric model on
# a discretization of the interval [0, 2*pi]. We would like to
# construct a D-efficient design.
# However, the distance of successive times of observations should not be
# smaller than the (1/72)-th of the interval (20 minutes, if the interval
# represents one day). Also, we cannot perform more than 8 observations
# in total, and more than 4 observations in the interval [2*pi/3, 2*pi]
# (i.e., during the 16 non-working hours).
# Create the matrix of regressors.
F.trig <- F.cube(~I(cos(x1)) + I(sin(x1)) +
I(cos(2 * x1)) + I(sin(2 * x1)) +
I(cos(3 * x1)) + I(sin(3 * x1)),
1 / 144, 2 * pi, 288)
# Create the constraints.
b.trig <- c(rep(1, 285), 12, 4)
A.trig <- matrix(0, nrow=287, ncol=288)
for(i in 1:285) A.trig[i, i:(i+3)] <- 1
A.trig[286,] <- 1; A.trig[287, 97:288] <- 1
# Compute the D-optimal approximate design under the constraints.
res.trig <- od.SOCP(F.trig, b.trig, A.trig, crit="D")
# Inspect the resulting approximate design.
od.plot(res.trig$w.best)
od.print(round(res.trig$w.best,2))
# It is clear that a very efficient exact design of size 8 satisfying
# the constraints performs the observations in design points
# 1, 34, 63, 96, 134, 173, 212, 251, i.e.
w.exact <- rep(0, 288)
w.exact[c(1, 34, 63, 96, 134, 173, 212, 251)] <- 1
# Indeed, the efficiency of this exact design relative to the optimal
# approximate design is:
od.crit(F.trig, w.exact) / od.crit(F.trig, res.trig$w.best)
# Of course, it is also possible to directly use an exact-design function
# such as od.MISOCP for this problem, or it is possible to use
# the optimal approximate design to decrease the support size of
# the candidate design points set, and then use an exact-design
# procedure.
# See also the examples in the help files of functions od.RC, od.IQP
# and od.MISOCP, where od.SOCP is used to compute informative lower
# bounds on the efficiencies of exact designs.
}
``` |

Questions? Problems? Suggestions? Tweet to @rdrrHQ or email at ian@mutexlabs.com.

All documentation is copyright its authors; we didn't write any of that.