Vignette for Package LA

This vignette gives a high level overview on how to use the R package LA to find different types of optimal designs.

devtools::load_all()
library(LA)

We start with order-of-addition designs (OofAs). Suppose we want to generate a full OofA with 4 factors. We can use the following code:

toy=rOofA(k=4);toy
#>       [,1] [,2] [,3] [,4]
#>  [1,]    1    2    3    4
#>  [2,]    1    2    4    3
#>  [3,]    1    3    4    2
#>  [4,]    1    3    2    4
#>  [5,]    1    4    2    3
#>  [6,]    1    4    3    2
#>  [7,]    2    3    4    1
#>  [8,]    2    3    1    4
#>  [9,]    2    4    1    3
#> [10,]    2    4    3    1
#> [11,]    2    1    3    4
#> [12,]    2    1    4    3
#> [13,]    3    4    1    2
#> [14,]    3    4    2    1
#> [15,]    3    1    2    4
#> [16,]    3    1    4    2
#> [17,]    3    2    4    1
#> [18,]    3    2    1    4
#> [19,]    4    1    2    3
#> [20,]    4    1    3    2
#> [21,]    4    2    3    1
#> [22,]    4    2    1    3
#> [23,]    4    3    1    2
#> [24,]    4    3    2    1

Sometimes, practitioners may not need full OofAs. If they want to generate a 12-run random OofA with 4 factors. The following code can be used:

toy=rOofA(n=12,k=4);toy
#>       [,1] [,2] [,3] [,4]
#>  [1,]    2    1    4    3
#>  [2,]    3    2    4    1
#>  [3,]    1    2    3    4
#>  [4,]    2    4    3    1
#>  [5,]    4    1    2    3
#>  [6,]    1    3    4    2
#>  [7,]    3    4    2    1
#>  [8,]    3    2    1    4
#>  [9,]    1    3    2    4
#> [10,]    3    1    4    2
#> [11,]    2    3    4    1
#> [12,]    4    2    3    1

Function “rOofA” is a fundamental support for main function “LA_OofA”. If we want to generate a D-optimal 12-run OofA with 4 factors. The following code can be used:

try.OofA=LA_OofA(n=12,k=4);try.OofA
#>
|
|                                                                      |   0%
|
|=======                                                               |  10%
|
|==============                                                        |  20%
|
|=====================                                                 |  30%
|
|============================                                          |  40%
|
|===================================                                   |  50%
|
|==========================================                            |  60%
|
|=================================================                     |  70%
|
|========================================================              |  80%
|
|===============================================================       |  90%
|
|======================================================================| 100%
#> [1] "average CPU time per iteration is: 0.02 seconds"
#> [1] "the number of iterations completed is: 10"
#>       [,1] [,2] [,3] [,4]
#>  [1,]    1    2    4    3
#>  [2,]    4    2    3    1
#>  [3,]    3    2    4    1
#>  [4,]    3    1    4    2
#>  [5,]    2    1    3    4
#>  [6,]    4    3    2    1
#>  [7,]    3    1    2    4
#>  [8,]    1    3    4    2
#>  [9,]    4    1    2    3
#> [10,]    2    1    4    3
#> [11,]    1    4    3    2
#> [12,]    4    1    3    2

Next, we introduce “LA_LHD” for identifying optimal Latin hypercube designs (LHDs). To generate a 6 by 3 maximin distance LHD, the following code can be used:

try.LHD1=LA_LHD(n=6,k=3);try.LHD1
#>
|
|                                                                      |   0%
|
|=======                                                               |  10%
|
|==============                                                        |  20%
|
|=====================                                                 |  30%
|
|============================                                          |  40%
|
|===================================                                   |  50%
|
|==========================================                            |  60%
|
|=================================================                     |  70%
|
|========================================================              |  80%
|
|===============================================================       |  90%
|
|======================================================================| 100%
#> [1] "average CPU time per iteration is: 0 seconds"
#> [1] "the number of iterations completed is: 10"
#>      [,1] [,2] [,3]
#> [1,]    2    3    2
#> [2,]    4    6    1
#> [3,]    3    2    6
#> [4,]    6    4    4
#> [5,]    1    5    5
#> [6,]    5    1    3

To generate a 9 by 4 nearly orthogonal LHD, the following code can be used:

try.LHD2=LA_LHD(n=9,k=4,OC="AvgAbsCor");try.LHD2
#>
|
|                                                                      |   0%
|
|=======                                                               |  10%
|
|==============                                                        |  20%
|
|=====================                                                 |  30%
|
|============================                                          |  40%
|
|===================================                                   |  50%
|
|==========================================                            |  60%
|
|=================================================                     |  70%
|
|========================================================              |  80%
|
|===============================================================       |  90%
|
|======================================================================| 100%
#> [1] "average CPU time per iteration is: 0 seconds"
#> [1] "the number of iterations completed is: 10"
#>       [,1] [,2] [,3] [,4]
#>  [1,]    6    6    1    5
#>  [2,]    8    1    9    7
#>  [3,]    3    2    3    8
#>  [4,]    7    8    4    3
#>  [5,]    1    5    8    1
#>  [6,]    2    9    7    9
#>  [7,]    9    7    6    6
#>  [8,]    5    4    5    2
#>  [9,]    4    3    2    4

To generate a 12 by 5 maximum projection LHD, the following code can be used:

try.LHD2=LA_LHD(n=12,k=5,OC="MaxProCriterion");try.LHD2
#>
|
|                                                                      |   0%
|
|=======                                                               |  10%
|
|==============                                                        |  20%
|
|=====================                                                 |  30%
|
|============================                                          |  40%
|
|===================================                                   |  50%
|
|==========================================                            |  60%
|
|=================================================                     |  70%
|
|========================================================              |  80%
|
|===============================================================       |  90%
|
|======================================================================| 100%
#> [1] "average CPU time per iteration is: 0 seconds"
#> [1] "the number of iterations completed is: 10"
#>       [,1] [,2] [,3] [,4] [,5]
#>  [1,]   10    2   11   11    5
#>  [2,]    1    6    2    6    8
#>  [3,]    8    9    7   12    2
#>  [4,]    7    5   12    2    7
#>  [5,]    6    1    3    3    3
#>  [6,]    3   10    1    1   11
#>  [7,]   11    4    6    4    1
#>  [8,]   12    8    4    7   12
#>  [9,]    5    3    9    8    9
#> [10,]    4   11    5   10    6
#> [11,]    9   12    8    5   10
#> [12,]    2    7   10    9    4

Last, we introduce “LA_opt” for experimental designs with continuous input and function optimization. For function optimization, we have two examples here. Now define an objective function: Sum of Different Powers.

$f(x)=\sum_{i=1}^{d}|x_{i}|^{i+1}$, where $$d$$ is the dimension.

SDP=function(x){i=1:length(x);y=sum(abs(x)^(i=1));return(y)}

Use LA to find the optimal solution under 20-dimensional setting for SDP function with 10 iterations:

try=LA_opt(of=SDP,lb=rep(-1,20),ub=rep(1,20),N=10,type="mini")
#>
|
|                                                                      |   0%
|
|=======                                                               |  10%
|
|==============                                                        |  20%
|
|=====================                                                 |  30%
|
|============================                                          |  40%
|
|===================================                                   |  50%
|
|==========================================                            |  60%
|
|=================================================                     |  70%
|
|========================================================              |  80%
|
|===============================================================       |  90%
|
|======================================================================| 100%
#> [1] "average CPU time per iteration is: 0.02 seconds"
#> [1] "the number of iterations completed is: 10"
SDP(try)  #Note that the true global optimum is 0, but we only have 10 iterations here
#> [1] 3.815988

Another example is the Cross-in-Tray function:

$f(x)=-0.0001(|sin(x_{1})sin(x_{2})exp\{|100-\frac{\sqrt{x_{1}^2+x_{2}^2}}{\pi} | \}|+1)^{0.1}$.

CiT=function(x){x1=x[1];x2=x[2];y=-0.0001*(abs(sin(x1)*sin(x2)*exp(abs(100-sqrt(x1^2+x2^2)/pi)))+1)^0.1;return(y)}

Use LA to find the optimal solution under 2-dimensional setting for CiT function with 10 iterations:

try2=LA_opt(of=CiT,lb=rep(-10,2),ub=rep(10,2),N=10,type="mini")
#>
|
|                                                                      |   0%
|
|=======                                                               |  10%
|
|==============                                                        |  20%
|
|=====================                                                 |  30%
|
|============================                                          |  40%
|
|===================================                                   |  50%
|
|==========================================                            |  60%
|
|=================================================                     |  70%
|
|========================================================              |  80%
|
|===============================================================       |  90%
|
|======================================================================| 100%
#> [1] "average CPU time per iteration is: 0 seconds"
#> [1] "the number of iterations completed is: 10"
CiT(try2) #Note that the true global optimum is -2.06261, but we only have 10 iterations
#> [1] -2.062471

For the optimal design part, we consider the following illustrative example:

A health scientist is designing an experiment to discover the relationship between cooking gas and lung function. The experiment is to be performed on rats. The response variable is maximum lung volume, and the predictor variable is time per day (in hours) exposed to a given level of cooking gas in the atmosphere. There are 20 rats available for the experiment. Assume a simple linear regression model.

Here we want to find a D-optimal 20-run design, where the input variable takes values between 0 and 24. In theory, we know the optimal design is the following:

matrix(c(rep(1,20),rep(0,10),rep(24,10)),ncol=2,nrow=20,byrow=FALSE)
#>       [,1] [,2]
#>  [1,]    1    0
#>  [2,]    1    0
#>  [3,]    1    0
#>  [4,]    1    0
#>  [5,]    1    0
#>  [6,]    1    0
#>  [7,]    1    0
#>  [8,]    1    0
#>  [9,]    1    0
#> [10,]    1    0
#> [11,]    1   24
#> [12,]    1   24
#> [13,]    1   24
#> [14,]    1   24
#> [15,]    1   24
#> [16,]    1   24
#> [17,]    1   24
#> [18,]    1   24
#> [19,]    1   24
#> [20,]    1   24

So Let’s see if LA is able to identify that. Now define the D-optimality criterion in simple linear regression model:

D=function(x){IM=t(x)%*%x;return(det(IM))}

Use LA to find the optimal solution for above problem, and we want to maximize the determinant of the information matrix.

try3=LA_opt(n=20,lb=c(1,0),ub=c(1,24),N=10,OC=D,type="maxi")
#>
|
|                                                                      |   0%
|
|=======                                                               |  10%
|
|==============                                                        |  20%
|
|=====================                                                 |  30%
|
|============================                                          |  40%
|
|===================================                                   |  50%
|
|==========================================                            |  60%
|
|=================================================                     |  70%
|
|========================================================              |  80%
|
|===============================================================       |  90%
|
|======================================================================| 100%
#> [1] "average CPU time per iteration is: 0 seconds"
#> [1] "the number of iterations completed is: 10"
try3 #with more iterations, LA would return even better result.
#>       [,1]        [,2]
#>  [1,]    1 24.00000000
#>  [2,]    1 11.60434053
#>  [3,]    1  3.70118599
#>  [4,]    1 10.50219986
#>  [5,]    1  2.61450112
#>  [6,]    1  0.80862910
#>  [7,]    1 23.76000000
#>  [8,]    1 24.00000000
#>  [9,]    1  1.51419156
#> [10,]    1 24.00000000
#> [11,]    1 24.00000000
#> [12,]    1  0.01700401
#> [13,]    1  9.08338327
#> [14,]    1  2.75103224
#> [15,]    1  9.23860924
#> [16,]    1 11.33949746
#> [17,]    1  0.11429296
#> [18,]    1 14.11705354
#> [19,]    1  3.34816666
#> [20,]    1  1.18134342