knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )

This vignette gives a high level overview on how to use the R package LHD to generate different types of efficient Latin hypercube designs (LHDs) with flexible sizes. Popular optimality criteria along with other useful functions will be introduced as well.

\

Load the library first.

devtools::load_all() library(LHD)

The overview starts with basic functions. The first function to be introduced is the "rLHD", which generates a random LHD with dimension $n$ by $k$. Suppose we want a random LHD with 6 rows and 3 columns, and denote it as X. We can run the following code.

set.seed(1) X=rLHD(n=6,k=3);X

If we want to know the inter-site distance between row $i$ and row $j$ of X ($i \neq j$), we could use function "dij". For example, the inter-site distance of the 2nd and the 4th row of X can be got via the following code.

dij(X=X,i=2,j=4) #The default setting is the rectangular distance. #If Euclidean distance is desired, run the following dij(X=X,i=2,j=4,q=2)

The formula of "dij" is given by the following

\begin{equation} d(\boldsymbol{x_i}, \boldsymbol{x_j})= (\sum_{l=1}^{k}|x_{il}-x_{jl}|^q)^{1/q}, \end{equation}

where $\boldsymbol{x_i}$ and $\boldsymbol{x_j}$ denote two design points from LHD matrix X. Having all the inter-site distances calculated, we could further calculate the $\phi_p$ criterion of X according to the formula from Jin, Chen and Sudjianto (2005).

\begin{equation} \phi_{p}= \bigg{\sum_{i=1}^{n-1}\sum_{j=i+1}^{n}d(\boldsymbol{x_i}, \boldsymbol{x_j})^{-p} \bigg} ^{1/p}. \end{equation}

The function "phi_p" provides calculation for above criterion. We can run the following code.

phi_p(X=X,p=15) #If Euclidean distance is desired, run the following phi_p(X=X,p=15,q=2)

$\phi_p$ criterion is popular for maximin distance LHDs. When orthogonal or nearly orthogonal LHDs are desired, average absolute correlation and maximum absolute correlation are widely used (Georgiou (2009)). Their formula are given by the following

\begin{equation} ave(|q|) = \frac{2 \sum_{i=1}^{k-1} \sum_{j=i+1}^{k}|q_{ij}|}{k(k-1)}, \, \, max(|q|)= max_{ij} |q_{ij}|, \end{equation}

where $q_{ij}$ is the columnwise correlation. The function "AvgAbsCor" and "MaxAbsCor" provide calculations for above criteria. We can run the following code.

AvgAbsCor(X=X) #The average absolute correlation of X MaxAbsCor(X=X) #The maximum absolute correlation of X

Another optimality criterion to be introduced is the maximum projection criterion (Joseph, Gul and Ba (2015)), which focuses on the projection property of the design matrix, and its formula is given by the following

\begin{equation} \psi = \Bigg{ \frac{1}{{n \choose 2}} \sum_{i=1}^{n-1} \sum_{j=i+1}^{n} \frac{1}{\Pi_{l=1}^{k}(x_{il}-x_{jl})^2} \Bigg}^{1/k}. \end{equation}

The function "MaxProCriterion" provides calculation for above criterion. We can run the following code.

MaxProCriterion(X=X) #The maximum projection criterion of X

Usually, a random LHD does not guarantee good space-filling property, orthogonality, nor full projection property. Morris and Mitchell (1995) proposed a version of the simulated annealing algorithm (SA) for generating maximin distance LHDs. Our function "SA" implements their algorithm, and we also added OC (optimality criterion) input argument in the "SA" function. Therefore, not only does "SA" generate maximin distance LHDs, but also it generates orthogonal and maximum projection LHDs. This is also true for the other algorithms in our package.

The function "exchange" makes "SA" workable since it is capable of switching two random elements from a given column of X. The following code shows examples of both "exchange" and "SA".

#Suppose we want to exchange two random elements from the 1st column of X. Xnew=exchange(X=X,j=1) #Look and compare X;Xnew #The SA function with default setting set.seed(1) trySA=SA(n=6,k=3) trySA #The phi_p of trySA is much smaller than the phi_p of X, which was 0.2517586 phi_p(X=trySA,p=15) #Now we try SA function with a different optimality criterion set.seed(1) trySA2=SA(n=6,k=3,OC="AvgAbsCor") AvgAbsCor(trySA2) #The average absolute correlation is about 0.05

Leary, Bhaskar and Keane (2003) modified the work of Morris and Mitchell (1995). They showed that orthogonal-array-based LHD (OALHD) with SA converges faster and generates better designs. Our function "OASA" implements their algorithm.

The function "OA2LHD" makes "OASA" workable since it is capable of transferring an Orthogonal Array (OA) into an LHD with corresponding size, based on the work of Tang (1993). The following code shows examples of both "OA2LHD" and "OASA".

#create an OA(9,2,3,2) OA=matrix(c(rep(1:3,each=3),rep(1:3,times=3)),ncol=2,nrow=9,byrow = F);OA #Transfer the OA above into a 9 by 2 LHD tryOA=OA2LHD(OA=OA) tryOA #The OASA function set.seed(1) tryOASA=OASA(OA=OA) tryOASA #The phi_p of tryOASA is smaller than the phi_p of SA phi_p(X=tryOASA,p=15); phi_p(X=SA(n=9,k=2),p=15) #Now we try OASA function with a different optimality criterion tryOASA2=OASA(OA=OA,OC="AvgAbsCor") AvgAbsCor(tryOASA2) #The average absolute correlation is 0

Joseph and Hung (2008) modified the work of Morris and Mitchell (1995) in another way that is able to reduce both $\phi_{p}$ and column-wise correlations at the same time. Our function "SA2008" implements their algorithm. See the following code and example.

set.seed(1) trySA2008_63=SA2008(n=6,k=3) trySA2008_63 phi_p(X=trySA2008_63,p=15) #Another example with different n and k trySA2008_92=SA2008(n=9,k=2) trySA2008_92 phi_p(X=trySA2008_92,p=15) #Now we try SA2008 function with a different optimality criterion set.seed(1) trySA2008=SA2008(n=6,k=3,OC="AvgAbsCor") AvgAbsCor(trySA2008) #The average absolute correlation is about 0.03

Ba, Myers and Brenneman (2015) extended the idea of Sliced Latin Hypercube Designs (SLHD) from Qian (2012) and had their own modifications of SA from Morris and Mitchell (1995). They called their algorithm "improved two-stage algorithm". Our function "SLHD" implements their algorithm. See the following code and example.

set.seed(1) trySLHD_63F=SLHD(n=6,k=3,t=2) #The default setting (stage2 is FALSE) trySLHD_63F phi_p(X=trySLHD_63F,p=15) #If the second stage is desired, run the following trySLHD_63T=SLHD(n=6,k=3,t=2,stage2=TRUE) trySLHD_63T phi_p(X=trySLHD_63T,p=15) #Result is improved (phi_p is smaller than above) #Now we try SLHD function with a different optimality criterion set.seed(1) trySLHD=SLHD(n=6,k=3,OC="AvgAbsCor",stage2=TRUE) AvgAbsCor(trySLHD) #The average absolute correlation is about 0.07

Chen et al. (2013) followed the structure of Particle Swarm Optimization (PSO) framework and proposed a version of it for constructing maximin distance LHDs, and their algorithm is called LaPSO. Our function "LaPSO" implements their algorithm. See the following code and example.

set.seed(1) tryLaPSO_63=LaPSO(n=6,k=3) tryLaPSO_63 phi_p(X=tryLaPSO_63,p=15) #Another example with different n and k tryLaPSO_92=LaPSO(n=9,k=2) tryLaPSO_92 phi_p(X=tryLaPSO_92,p=15) #Now we try LaPSO function with a different optimality criterion set.seed(1) tryLaPSO=LaPSO(n=6,k=3,OC="AvgAbsCor") AvgAbsCor(tryLaPSO) #The average absolute correlation is about 0.1

Liefvendahl and Stocki (2006) proposed a version of the genetic algorithm (GA) which implemented an elite strategy to focus on the global best directly for constructing maximin distance LHDs. Our function "GA" implements their algorithm. See the following code and example.

set.seed(1) tryGA_63=GA(n=6,k=3) tryGA_63 phi_p(X=tryGA_63,p=15) #Another example with different n and k. tryGA_92=GA(n=9,k=2) tryGA_92 phi_p(X=tryGA_92,p=15) #Now we try GA function with a different optimality criterion set.seed(1) tryGA=GA(n=6,k=3,OC="AvgAbsCor") AvgAbsCor(tryGA) #The average absolute correlation is about 0.07

In addition to above algorithms, Wang, L., Xiao, Q., and Xu, H. (2018) proposed a construction method for maximin $L_1$ distance LHDs based on good lattice point designs. One big advantage is that their method gives instant solutions, and this is very helpful for practitioners who seek immediate results. See the following code and example.

#n by n design when 2n+1 is prime try=FastMmLHD(8,8) try phi_p(try) #calculate the phi_p of "try". #n by n design when n+1 is prime try2=FastMmLHD(12,12) try2 phi_p(try2) #calculate the phi_p of "try2". #n by n-1 design when n is prime try3=FastMmLHD(7,6) try3 phi_p(try3) #calculate the phi_p of "try3". #General cases try4=FastMmLHD(24,8) try4 phi_p(try4) #calculate the phi_p of "try4".

Besides maximin distance LHDs, orthogonal Latin hypercube designs (OLHDs) are popular as well. Ye (1998) proposed a construction method for generating OHLDs. Our function "OLHD.Y1998" implements this method. See the following code and example.

#create an orthogonal LHD with m=4. So n=2^m+1=17 and k=2*m-2=6 tryOLHD=OLHD.Y1998(m=4) tryOLHD MaxAbsCor(tryOLHD) #zero columnwise correlations

Cioppa and Lucas (2007) extended Ye's idea, and their method is able to accommodate more factors with the same run size. Our function "OLHD.C2007" implements this method. See the following code and example.

#create an orthogonal LHD with m=4. So n=2^m+1=17 and k=4+3=7 tryOLHD=OLHD.C2007(m=4) tryOLHD MaxAbsCor(tryOLHD) #zero columnwise correlations

Lin et al. (2009) proposed to couple OLHDs with OAs to generate OLHDs with large factor size. Our function "OLHD.L2009" implements this method. See the following code and example.

#create a 5 by 2 OLHD OLHD=OLHD.C2007(m=2) #create an OA(25,6,5,2) OA=matrix(c(2,2,2,2,2,1,2,1,5,4,3,5,3,2,1,5,4,5,1,5,4,3,2,5, 4,1,3,5,2,3,1,2,3,4,5,2,1,3,5,2,4,3,1,1,1,1,1,1,4,3,2,1,5,5, 5,5,5,5,5,1,4,4,4,4,4,1,3,1,4,2,5,4,3,3,3,3,3,1,3,5,2,4,1,3, 3,4,5,1,2,2,5,4,3,2,1,5,2,3,4,5,1,2,2,5,3,1,4,4,1,4,2,5,3,4, 4,2,5,3,1,4,2,4,1,3,5,3,5,3,1,4,2,4,5,2,4,1,3,3,5,1,2,3,4,2, 4,5,1,2,3,2),ncol=6,nrow=25,byrow=T) #Construct a 25 by 12 OLHD tryOLHD=OLHD.L2009(OLHD,OA) tryOLHD MaxAbsCor(tryOLHD) #zero columnwise correlations

Sun et al. (2010) proposed to construct OLHDs with flexible run sizes via matrices manipulations. Our function "OLHD.S2010" implements this method. See the following code and example.

#create an orthogonal LHD with C=3, r=3, type="odd". #So n=3*2^4+1=49 and k=2^3=8 tryOLHD1=OLHD.S2010(C=3,r=3,type="odd") tryOLHD1 #create an orthogonal LHD with C=3, r=3, type="even". #So n=3*2^4=48 and k=2^3=8 tryOLHD2=OLHD.S2010(C=3,r=3,type="even") tryOLHD2 MaxAbsCor(tryOLHD1) #zero columnwise correlations MaxAbsCor(tryOLHD2) #zero columnwise correlations

Butler (2001) proposed to construct OLHDs using Williams transformation. Our function "OLHD.B2001" implements this method. See the following code and example.

#create an orthogonal LHD with n=11 and k=5 OLHD.B2001(n=11,k=5) #create an orthogonal LHD with n=7 and k=6 OLHD.B2001(n=7,k=6)

**Any scripts or data that you put into this service are public.**

Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.