function to create SOAs of strength 2+ from regular s-level designs

creates an array in s^k runs with columns in s^2 levels for prime or prime power s

SOAs2plus_regular(
  s,
  k,
  m = NULL,
  orth = TRUE,
  old = FALSE,
  noptim.rounds = 1,
  noptim.repeats = 1,
  optimize = TRUE,
  dmethod = "manhattan",
  p = 50
)

Arguments

s: prime or prime power
k: array will have n=s^k runs; for s=2, k>=4 is needed, for s>2, k>=3 is sufficient
m: optional integer: number of columns requested; if NULL, the maximum possible number of columns is created, which is (s^k-1)/(s-1) - ((s-1)^k-1)/(s-2) for s>2 and s^k-s^k1 - s^(k-k1) + 2, with k1=floor(k/2), for s=2; specifying a smaller m is beneficial not only for run time but also for possibly achieving a column-orthogonal array (see Details section)
orth: logical: if FALSE, suppresses attempts for orthogonal columns and selects the first permissible column for each column of B (see Details section)
old: logical, relevant for orth=TRUE only: if TRUE, limits possible columns for B to the columns not eligible for A (instead of the columns not used in A); should only be used for reproducing designs created by version 1.1 or earlier
noptim.rounds: the number of optimization rounds for each independent restart
noptim.repeats: the number of independent restarts of optimizations with noptim.rounds rounds each
optimize: logical: should optimization be applied? default TRUE
dmethod: method for the distance in phi_p, "manhattan" (default) or "euclidean"
p: p for phi_p (the larger, the closer to maximin distance)

Value

matrix of class SOA with the attributes that are listed below. All attributes can be accessed using function attributes, or individual attributes can be accessed using function attr. These are the attributes:

type: the type of array (SOA or OSOA)
strength: character string that gives the strength
phi_p: the phi_p value (smaller=better)
optimized: logical indicating whether optimization was applied
permpick: matrix that lists the id numbers of the permutations used
perms2pickfrom: optional element, when optimization was conducted: the overall permutation list to which the numbers in permlist refer
call: the call that created the object

Details

The construction is by He, Cheng and Tang (2018), Prop.1 (C2) / Theorem 2 for s=2 and Theorem 4 for s>2.
B is chosen as an OA of strength 2, if possible, which yields orthogonal columns according to Zhou and Tang (2019). This is implemented using a matching algorithm for bipartite graphs from package igraph; the smaller m, the more likely that orthogonality can be achieved. However, strength 2+ SOAs are not usually advisable for m small enough that a strength 3 OA exists.
Optimization according to Weng has been added (separate level permutations in columns of A and B, noptim.rounds times). Limited tests suggest that a single round (noptim.rounds=1) often does a very good job (e.g. for s=2 and k=4), and further rounds do not yield too much improvement; there are also cases (e.g. s=5 with k=3), for which the unoptimized array has a better phi_p than what can be achieved by most optimization attempts from a random start.

The search for orthogonal columns can take a long time for larger arrays, even without optimization. If this is prohibitive (or not considered valuable), orth=FALSE causes the function to create the matrix B for equation D=2A+B with less computational effort.
The subsequent optimization, if not switched off, is of the same complexity, regardless of the value for orth. Its duration heavily depends on the number of optimization steps that are needed before the algorithm stops. This has not been systematically investigated; cases for which the total run time with optimization is shorter for orth=TRUE than for orth=FALSE have been observed.

With package version 1.2, the creation of SOAs has changed: Up to version 1.1, the columns of B were chosen only from those columns that were not eligible for A, whereas the new version chooses them from those columns that are not used for A. This increases the chance to achieve geometrically orthogonal columns.
Users who want to reproduce a design from an earlier version can use argument old.

Note

Strength 2+ SOAs can accommodate a large number of factors with reasonable stratified balance behavior. Note that their use is not usually advisable for m small enough that a strength 3 OA with s^2 level factors exists.

References

For full detail, see SOAs-package.

Groemping (2023a) He, Cheng and Tang (2018)
Weng (2014)
Zhou and Tang (2019)

Author

Ulrike Groemping

Examples

# \donttest{
## unoptimized OSOA with 8 16-level columns in 64 runs
## (maximum possible number of columns)
plan64 <- SOAs2plus_regular(4, 3, optimize=FALSE)
ocheck(plan64)   ## the array has orthogonal columns
#> [1] TRUE

## optimized SOA with 20 9-level columns in 81 runs
## (up to 25 columns are possible)
plan <- SOAs2plus_regular(3, 4, 20)
#> Optimization 1 of 1 started
#> Optimization round 1 of 1 started
## many column pairs have only 27 level pairs covered
count_npairs(plan)
#> $paircounts
#>   [1] 27 27 27 27 27 27 27 81 27 27 27 81 27 27 27 27 27 27 27 27 27 27 27 27 27
#>  [26] 81 81 81 27 81 81 81 81 81 81 81 27 27 27 27 27 27 27 81 81 27 27 81 81 81
#>  [51] 81 81 81 27 27 27 81 81 81 27 81 81 81 27 27 81 81 27 27 81 27 81 81 27 27
#>  [76] 81 81 27 27 27 81 81 27 27 81 81 81 27 81 27 81 27 81 81 27 27 81 81 81 27
#> [101] 81 27 27 27 27 81 81 81 81 81 81 27 81 27 27 27 81 81 81 81 81 81 81 27 81
#> [126] 81 27 27 81 27 81 81 81 27 27 27 81 27 27 27 81 81 27 27 81 81 81 81 81 27
#> [151] 27 81 81 81 81 27 27 27 81 81 81 27 81 81 81 81 81 81 81 27 27 81 27 27 81
#> [176] 27 81 27 27 81 27 81 81 81 27 27 27 27 27 27
#> 
#> $columnpaircounts
#>  [1]  621 1053  945  999  891  999 1053 1107 1107  945 1161  999 1215 1053  999
#> [16] 1161 1161 1053  999  999
#> 
## an OA would exist for 10 9-level factors (DoE.base::L81.9.10)
## it would cover all pairs
## (SOAs are not for situations for which pair coverage
## is of primary interest)
# }