Code run with :

R.Version()
## $platform
## [1] "x86_64-pc-linux-gnu"
## 
## $arch
## [1] "x86_64"
## 
## $os
## [1] "linux-gnu"
## 
## $system
## [1] "x86_64, linux-gnu"
## 
## $status
## [1] ""
## 
## $major
## [1] "4"
## 
## $minor
## [1] "3.1"
## 
## $year
## [1] "2023"
## 
## $month
## [1] "06"
## 
## $day
## [1] "16"
## 
## $`svn rev`
## [1] "84548"
## 
## $language
## [1] "R"
## 
## $version.string
## [1] "R version 4.3.1 (2023-06-16)"
## 
## $nickname
## [1] "Beagle Scouts"

Building sampling designs

The sampling in an equilateral triangle embedded within the trigonometric circle (i.e. center at \((0,0)\) and radius of \(1\)). One tip of the triangle is located at coordinates \((0,1)\). The sampling area is visualized below :

drawFrame <- function(){
  vecTheTip <- c(-pi/6,pi/2,7*pi/6) #polar angles of sampling area tips
  vecAbsTip <- cos(vecTheTip) #abscissas of tips
  vecOrdTip <- sin(vecTheTip) #ordinates of tips
  plot(cos(seq(0,2*pi,0.01)),sin(seq(0,2*pi,0.01)),lty="dashed",type="l",
       xlab="Longitude",ylab="Latitude",asp=1)
  lines(vecAbsTip[c(1:3,1)],vecOrdTip[c(1:3,1)])}
drawFrame()

where dashed lines show the trigonometric circle while continuous lines show the sampling area.

Random design

Random sampling of a point within the triangular area depicted above is not completely straightforward. We use a two-step procedure, first sampling the abscissa \(X\) of the point, then the ordinate \(Y\).

Abscissa \(X\) has the following probability distribution function: \[ \begin{array}{l} \mathbb{P}(X<x)=2(x+\sqrt{3}/2)^2/3 \quad \text{if}\quad x\leq 0 \\ \mathbb{P}(X<x)=1-2(\sqrt{3}/2-x)^2/3 \quad \text{if}\quad x> 0 \end{array} \]

Abscissa \(X\) can therefore be simulated applying the inverse of its distribution function to a uniform random variable \(U\) on \(\left[0,1\right]\) : \[ \begin{array}{l} X = \frac{\sqrt{3}}{2}(\sqrt{2U}-1) \quad \text{if}\quad U \leq 1/2 \\ X=\frac{\sqrt{3}}{2}(1-\sqrt{2(1-U)}) \quad \text{if}\quad U > 1/2 \end{array} \]

Ordinate \(Y\) can then be simulated conditionally to \(X\) through a uniform sampling in the height associated to abscissa \(X\) :

\[ \begin{array}{l} Y = V (X+\sqrt{3}/2) \sqrt{3}-1/2 \quad \text{if}\quad X \leq 0 \\ Y=V (\sqrt{3}/2-X) \sqrt{3}-1/2 \quad \text{if}\quad X > 0 \end{array} \]

The function used to generate random sampling design is :

trRa <- function(nTo){
  
  #OBJECTIVE : generate a random sampling design within a triangle within circle 
  #of radius 1
  
  #INPUTS :
  #nTo : number of sites in the design
  
  #OUTPUT :
  #tabZ : coordinates of sampling points; tabZ[1,i] (resp. tabZ[2,i]) is the 
    #longitude (resp. latitude) of the ith points
  
  if(nTo==0) return(NULL)
  vU <- runif(nTo,0,1); vV <- runif(nTo,0,1) 
  tabZ <- sapply(1:nTo,function(re){
    u <- vU[re]; v <- vV[re]
    if(u<1/2){x <- sqrt(3)/2*(sqrt(2*u)-1); y <- v*(x+sqrt(3)/2)*sqrt(3)-1/2} 
    else{x <- sqrt(3)/2*(1-sqrt(2*(1-u))); y <- v*(sqrt(3)/2-x)*sqrt(3)-1/2}
    return(c(x,y))})
  return(tabZ)}

We also define a generic function to plot designs :

plotDes <- function(tabZ){

  #OBJECTIVE :
    #plots a sampling design
  
  #INPUTS :
    #tabZ is the table of sites coordinates (see tabZ in functions below)
  
  drawFrame()
  points(tabZ[1,],tabZ[2,],pch="+")}

We plot an example of random design with \(100\) points:

tabZ_exTr <- trRa(100)
plotDes(tabZ_exTr)

Grid design

The grid design considered in our study has a triangular mesh. The mesh is equilateral. The side length of the sampling area (\(\sqrt{3}\)) is a multiple of grid mesh size, which makes grid design exactly embedded in the sampling area.

For a grid design with mesh size \(\sqrt{3}/n\), \(n\in \mathbb{N}\), sampling points are all the points with position of the type \[\left(\begin{array}{l}-\frac{\sqrt{3}}{2} \\ -\frac{1}{2} \end{array}\right) +k\left(\begin{array}{l}\frac{\sqrt{3}}{n} \\ 0 \end{array}\right) +l\left(\begin{array}{l}\frac{\sqrt{3}}{2n} \\ \frac{3}{2n} \end{array}\right)\] where \(k,l \in \mathbb{N}\) and \(k+l\leq n\). There are \(\frac{(n+1)(n+2)}{2}\) such points.

The function below generates the coordinates of the sampling points of a grid design associated to a target number of sampling points on the side of sampling area (i.e. using notations above \(n+1\)):

trGr <- function(nSi){
  
  #OBJECTIVE : 
    #generate a triangular grid sampling design within a circle 
    #of radius 1
  
  #INPUTS :
    #nSi : number of sites on the side of the sampling area
  
  #OUTPUT :
    #tabZ : coordinates of sampling points; tabZ[1,i] (resp. tabZ[2,i]) is the 
      #longitude (resp. latitude) of the ith points
  
  vecK <- unlist(sapply(1:nSi,function(i) 1:i))-1
  vecL <- unlist(sapply(1:nSi,function(i) rep(nSi-i+1,i)))-1
  nTo <- nSi*(nSi+1)/2
  tabZ <- sapply(1:nTo,function(i){
    k <- vecK[i]; l <- vecL[i]
    x <- -sqrt(3)/2+k*sqrt(3)/(nSi-1)+l*sqrt(3)/(2*(nSi-1))
    y <- -1/2+l*3/(2*(nSi-1))
    return(c(x,y))})
  return(tabZ)}

Here is an example of simulated design with \(14\) points on the side of sampling area (i.e. \(105\) sampling sites in total):

tabZ_exGr <- trGr(14)
plotDes(tabZ_exGr)

Hybrid design

Hybrid designs are obtained through merging sampling sites from a grid design and sampling sites from a random design. A hybrid design is defined by its target number of sampling sites \(n\) and the fraction \(p\) of the sampling sites that comes from a random design (the rest coming from a grid design). Number of sites coming from grid design is rounded if needed, although throughout our study we use values of \(p\) that lead to integer number of sites.

The grid design from which the \((1-p)n\) sampling sites are taken is the smallest grid design that has enough sampling sites to make a hybrid design with \(p=0\) possible. The mesh size of this grid design is \(\sqrt{3}/n_g\), where \(n_g=\lceil x_g \rceil\) and \(x_g\) is the unique positive root of : \[\frac{(x_g+1)(x_g+2)}{2}=n\] which is equal to: \[x_g=\frac{-3+\sqrt{1+8n}}{2}\] The function below generates the coordinates of the sampling points of a grid design associated to a target number of sampling points \(n\) and a fraction of random points \(p\):

trHy <- function(nTo,p){
  
  #OBJECTIVE : 
    #generate a hybrid sampling design within a circle 
    #of radius 1
  
  #INPUTS :
    #nTo : total number of sites
    #p   : fraction of sites coming from a random design
  
  #OUTPUT :
    #tabZ : coordinates of sampling points; tabZ[1,i] (resp. tabZ[2,i]) is the 
      #longitude (resp. latitude) of the ith points
  
  nSi_min <- ceiling((sqrt(1+8*nTo)-1)/2) #minimum number of sampling sites on 
                                          #the side of grid sampling design to 
                                          #make p=0 possible
  desGr0 <- trGr(nSi_min)
  tabZ_1 <- desGr0[,sample(1:ncol(desGr0),round((1-p)*nTo),replace=FALSE)]
  tabZ_2 <- trRa(nTo-round((1-p)*nTo))
  tabZ <- cbind(tabZ_1,tabZ_2)  
  return(tabZ)}

Here is an example of simulated design with \(100\) points and \(20\%\) of random points (i.e. \(p=0.2\)):

tabZ_exHy <- trHy(100,0.2)
plotDes(tabZ_exHy)

Figure 1 of maintext

Using the code to generate hybrid designs above, we can generate the Figure 1 of maintext :

nTo <- 27; vecP <- c(0,5/27,1)
vLab <- paste("p =",c(0,round(5/27,2),1))
LL <- lapply(vecP,function(p) trHy(nTo,p))
par(mfrow=c(1,3),cex.lab=1.4,cex.axis=1,lwd=1.5,cex=1.2)
nul <- sapply(1:3,function(j){
  lis <- LL[[j]]; plotDes(lis)
  mtext(vLab[j],3,adj=0.5,padj=-1,font=2,cex=1.7)})

## png 
##   2

Fractal design

The triangular fractal sampling designs considered in our study are generated by iteratively applying three similarities on the set of sampling sites. It is simpler to present these similarities in the complex plane : \[ \begin{array}{l} s_1 : z\rightarrow \rho z +(1-\rho)e^{-i\frac{\pi}{6}} \\ s_2 : z\rightarrow \rho z +(1-\rho)e^{i\frac{\pi}{2}} \\ s_3 : z\rightarrow \rho z +(1-\rho)e^{-i\frac{5\pi}{6}} \end{array} \]

where \(\rho\) is a positive real parameter that drives the level of contraction at each iteration We implement these similarities using the following generic function coding for similarities in the complex plane :

IS <- function(z,a,b){
  
  #OBJECTIVE : 
    #apply the similarity z -> s = a*z+b on the complex number z
  
  #INPUTS :
    #z : a complex number; z[1] and z[2] are the real and imaginary part of z
    #a : a complex number; a[1] and a[2] are the real and imaginary part of a
    #b : a complex number; b[1] and b[2] are the real and imaginary part of b
    
  #OUTPUT :
    #s : a complex number; s[1] and s[2] are the real and imaginary part of b
  
  s <- c(a[1]*z[1]-a[2]*z[2]+b[1],a[1]*z[2]+a[2]*z[1]+b[2])
  return(s)
}

Applying an iteration on a set of sampling points \(\mathcal{S}\) consists in switching from \(\mathcal{S}\) to the set \(s_1(\mathcal{S}) \cup s_2(\mathcal{S}) \cup s_3(\mathcal{S})\). Here again, we provide a generic function applying an iteration of a set of similarities on a set of sampling points :

nextGen <- function(tabZ,tabA,tabB){
  
  #OBJECTIVE :
    #apply a set of similarities on a set of complex numbers
  
  #INPUTS :
    #cf. IS function documentation for notations
    #tabZ : a set of complex numbers z ; tabZ[1,i] (resp. tabZ[2,i]) is the 
      #real (resp. imaginary) part of the ith complex number z_i
    #tabA : a set of complex numbers a ; tabA[1,j] (resp. tabA[2,j]) is the 
      #real (resp. imaginary) part of the ith complex number a_j
    #tabB : a set of complex numbers b ; tabB[1,j] (resp. tabB[2,j]) is the 
      #real (resp. imaginary) part of the jth complex number b_j
    
  #OUTPUT : 
    #tabS : a set of complex numbers s ; tabS[1,k] (resp. tabS[2,k]) is the 
      #real (resp. imaginary) part of the kth complex number s_k
    
  tabS <- NULL
  for(j in 1:ncol(tabB)){
    a <- tabA[,j]; b <- tabB[,j]
    tabZ2 <- apply(tabZ,2,function(z) IS(z,a,b))
    tabS <- cbind(tabS,tabZ2)}
  return(tabS)
}

The iterative process is initiated on a single sampling point at the center of the area of study. We provide an example of a sequence of iterations below:

par(mfrow=c(2,3))
zIni <- rbind(0,0); rho <- 0.4
#Similarity parameters corresponding to s_1,s_2,s_3
tabA <- rbind(rep(rho,3),rep(0,3))
tabB <- rbind(sqrt(3)/2*(1-rho)*c(-1,0,1),(1-rho)*c(-1/2,1,-1/2))
tabZ <- zIni
drawFrame()
points(tabZ[1,],tabZ[2,],pch="+")
title(main=paste("Iteration",0))
for(i in 1:5){
  tabZ <- nextGen(tabZ,tabA,tabB)
  plotDes(tabZ)
  title(main=paste("Iteration",i))}

After \(k\) iterations, the design harbours \(3^k\) sampling points. We embedded the iteration scheme in a function that applies the iteration scheme a given number of times with a given value of parameter \(\rho\) :

trFr <- function(nIt,rho=sqrt(3)/(sqrt(3)+2)){
  
  #OBJECTIVE :
    #generate a triangular fractal sampling design within a circle 
    #of radius 1
  
  #INPUTS :
    #nIt : number of iteration of the similarities (total number of sites in
      #the design is 3^nIt)
    #rho : contraction used in similarities
  
  #OUTPUT :
    #tabZ : coordinates of sampling points; tabZ[1,i] (resp. tabZ[2,i]) is the 
      #longitude (resp. latitude) of the ith points
  
  if(rho>sqrt(3)/(sqrt(3)+2)) print("Fractal self-overlaps : decrease rho !")
  zIni <- rbind(0,0)
  if(nIt==0) return(zIni)
  tabA <- rbind(rep(rho,3),rep(0,3))
  tabB <- rbind(sqrt(3)/2*(1-rho)*c(-1,0,1),(1-rho)*c(-1/2,1,-1/2))
  tabZ <- zIni
  for(i in 1:nIt){
    tabZ <- nextGen(tabZ,tabA,tabB)}
  return(tabZ)
}

Note that there exists an upper threshold on \(\rho\) above which the output of the three similarities overlap, which is an undesirable property here. This upper threshold is denoted \(\rho_{max}\) and is found as the solution of: \[2(1-\rho_{max})\sin(\frac{\pi}{3})=2\rho_{max}\] where the left-hand side is the side length of the triangle obtained when applying the three similarities to the center of the sampling area. The solution is \(\rho_{max}=\frac{\sqrt{3}}{2+\sqrt{3}}\approx 0.46\). We show the overlapping phenomenon on three iterations below:

tabZ_noOv <- trFr(3,rho=0.3)
tabZ_ov <- trFr(3,rho=0.6)
## [1] "Fractal self-overlaps : decrease rho !"
par(mfrow=c(1,2))
plotDes(tabZ_noOv)
title(main=paste("No overlap - rho =",0.3))
plotDes(tabZ_ov)
title(main=paste("Overlap - rho =",0.6))

Avoiding overlap is the reason why we redefined a contraction parameter \(x\) in main text which is a fraction of \(\rho_{max}\).

Generating Figure 2 of main text

Using the code to generate fractal designs developped above, we can generate Figure 2 of main text :

rhoMax <- sqrt(3)/(2+sqrt(3))
nTo <- 27; vX <- c(0.33,0.66,1)
vLab <- paste("x =",vX)
lisDesFr <- lapply(vX,function(x) trFr(3,rho=x*rhoMax))
par(mfrow=c(1,3),cex.lab=1.5,cex.axis=1.1,lwd=1.6,cex=1.3)
nul <- sapply(1:3,function(j){
  lis <- lisDesFr[[j]]; plotDes(lis)
  mtext(vLab[j],3,adj=0.5,padj=-1,font=2,cex=1.8)})

## png 
##   2

Generating all the designs used in main text

In our study, we consider designs with \(27\) sampling points.

nTo <- 27

This leads to \(28\) possible values of \(p\) parameter in hybrid designs :

vecP <- (0:nTo)/nTo

For each value of \(p\) strictly above 0, we consider \(30\) replicates of the hybrid design, to account for fluctuations in the random points:

nRep <- 120
vecPRep <- rep(vecP,each=nRep)

We also consider \(28\) values of \(x\) for fractal designs, equally spaced on a log-scale between \(10^{-1.5}\) and \(1\) :

vecX <- 10^(-1.5+1.5*(0:nTo)/nTo)

We combine the parameters of all the designs tested in a single table:

tabParDes <- data.frame(
  typ=c(rep("frac",nTo+1),rep("mix",nRep*(nTo+1))),
  par=c(vecX,vecPRep))

We generate all the designs:

set.seed(446) #for replicability
lisDes <- lapply(1:nrow(tabParDes),function(iRowDes){
  if(tabParDes$typ[iRowDes]=="frac"){
    tabXY <- trFr(round(log(nTo)/log(3)),tabParDes$par[iRowDes]*sqrt(3)/(2+sqrt(3)))
  }else tabXY <- trHy(nTo,tabParDes$par[iRowDes])
  matD <- as.matrix(dist(t(tabXY)))
  des <- list(xy=tabXY,dis=matD)
  return(des)})

The total number of designs considered in our study is:

length(lisDes)
## [1] 3388

Generating environmental covariate

Each sampling point of a design is associated to a value of environmental variable. We considered three main types of environmental gradient in our analysis, a linear trend across the sampling area, rugged environmental changes within the area and an edge effect. All the variables are standardized between \(0\) and \(1\) wihtin the area.

We define a dense grid over which environmental covariates will be visualized below.

tabZ_envEx <- trGr(100)

The environmental profile is generated using a product of two sine waves along abscissas and ordinates axis respectively. Formally, this product is a two-dimensional function \(\Phi\) which associates to any position \((x,y)\) a value: \[\Phi(x,y) = \sin(2\pi x)\sin(2\pi y)\]

We visualize this function below:

PhiTest <- function(x,y) sin(2*pi*x)*sin(2*pi*y)
vecE_Phi <- apply(tabZ_envEx,2,function(v) PhiTest(v[1],v[2]))
vecE_Phi_std <- (vecE_Phi-min(vecE_Phi))/(max(vecE_Phi)-min(vecE_Phi))
drawFrame()
points(tabZ_envEx[1,],tabZ_envEx[2,],pch=24,cex=0.75,
  bg=rgb(vecE_Phi_std,0,1-vecE_Phi_std),col=rgb(0,0,0,0))
title(main=expression(Phi))

We then translate this function, considering \(\Psi\):

\[\Psi_{\phi}(x,y)=\Phi(x-\frac{\phi}{4},y+\frac{1}{4})\]

Horizontal translation is driven by a parameter \(\phi\in\left[0,1\right]\).

PsiTest <- function(x,y,phi) PhiTest(x-phi/4,y+1/4)
vecE_Psi_1 <- apply(tabZ_envEx,2,function(v) PsiTest(v[1],v[2],1))
vecE_Psi_1_std <- (vecE_Psi_1-min(vecE_Psi_1))/(max(vecE_Psi_1)-min(vecE_Psi_1))
vecE_Psi_0 <- apply(tabZ_envEx,2,function(v) PsiTest(v[1],v[2],0))
vecE_Psi_0_std <- (vecE_Psi_0-min(vecE_Psi_0))/(max(vecE_Psi_0)-min(vecE_Psi_0))
par(mfrow=c(1,2))
drawFrame()
points(tabZ_envEx[1,],tabZ_envEx[2,],pch=24,cex=0.75,
  bg=rgb(vecE_Psi_1_std,0,1-vecE_Psi_1_std),col=rgb(0,0,0,0))
title(main=paste("Psi, phi =",1))
drawFrame()
points(tabZ_envEx[1,],tabZ_envEx[2,],pch=24,cex=0.75,
  bg=rgb(vecE_Psi_0_std,0,1-vecE_Psi_0_std),col=rgb(0,0,0,0))
title(main=paste("Psi, phi =",0))

We can then apply a dilatation coefficient \(\lambda\) between \(0\) and \(1\) such that the environmental variable is ultimately defined as :

\[E_{\phi,\lambda}(x,y)=\Psi_\phi(\frac{x}{4\lambda},\frac{y}{4\lambda})\]

To obtained a standardized environmental covariate between \(0\) and \(1\), we use: \[\tilde{E}_{\phi,\lambda}(x,y) =\frac{E_{\phi,\lambda}(x,y)-E_{min}}{E_{max}-E_{min}}\] where \(E_{min}\) and \(E_{max}\) are numerically computed based on a dense grid over the sampling area. The standardized rugged environmental covariate \(\tilde{E}\) is coded as:

getEnv_rug <- function(tabZ,lam,phi,ref=NULL){
  #OBJECTIVE :
    #computing a rugged environment covariate value for each site of a sampling 
    #design
  
  #INPUTS :
    #tabZ : coordinates of sampling points; tabZ[1,i] (resp. tabZ[2,i]) is the 
      #longitude (resp. latitude) of site i
    #lam  : dilatation coefficient of the covariate (the size of a patch is 
      #lambda*sqrt(3))
  #OUTPUT : 
    #vecE : environment covariate values, vecE[i] is the environment at site i
    
    vecE <- apply(tabZ,2,function(v) PsiTest(v[1]/(4*lam),v[2]/(4*lam),phi))
    if(is.null(ref)){
      tabZ_envEx <- trGr(100)
      vecE_ref <- apply(tabZ_envEx,2,
        function(v) PsiTest(v[1]/(4*lam),v[2]/(4*lam),phi))}
    else vecE_ref <- ref
    E_min <- min(vecE_ref); E_max <- max(vecE_ref)
    vecE_std <- round((vecE-E_min)/(E_max-E_min),4)
    return(vecE_std)}

Figure 3 of main text

In Figure 3, we visualize environmental variables associated to \(\lambda,\phi\) parameters considered in main text (panels with a green frame are those analyzed in main text):

vecLam <- sqrt(3)/12*c(sqrt(3)/12,1,sqrt(12/sqrt(3)),12/sqrt(3))
vecPhi <- c(0,0.5,1)
par(mfrow=c(3,3))
nul <- sapply(vecPhi,function(phi){
  sapply(vecLam[-1],function(lam){
    vecE_rug <- getEnv_rug(tabZ_envEx,lam,phi)
    drawFrame()
    points(tabZ_envEx[1,],tabZ_envEx[2,],pch=24,cex=0.75,
      bg=rgb(vecE_rug,0,1-vecE_rug),col=rgb(0,0,0,0))
    title(main=
      paste("lambda =",round(lam,2),", phi =",round(phi,2)))
    if((phi==0)+(lam==1)>=1) box(lty=2-(phi<1)*(lam==1),lwd=8)
    })})

## png 
##   2

We compile the \(12\) combinations of parameters depicted in Figure 3 in a table :

vecLam <- sqrt(3)/12*c(sqrt(3)/12,1,sqrt(12/sqrt(3)),12/sqrt(3))
vecPhi <- c(0,0.5,1)
tabParEnv <- as.matrix(expand.grid(lam=vecLam,phi=vecPhi))

For each combination, we compute before-hand a dense grid of non-standardized environmental values that will be used to standardize the environmental covariate. This will avoid useless repetition of computation :

tabERef <- apply(tabParEnv,1,function(envPar){
  lam <- envPar[1]; phi <- envPar[2]
  tabZ_envEx <- trGr(100)
  vecE_ref <- apply(tabZ_envEx,2,
        function(v) PsiTest(v[1]/(4*lam),v[2]/(4*lam),phi))
  return(vecE_ref)})

Computing the RRMSEs for all sampling design, environmental covariate and autocorrelation range values

We build a function that computes RRMSEs associated to a couple design-environment using formulas of main text :

getTheoRRMSE <- function(des,as,sig=1){

  #OBJECTIVE : 
    #compute theoretical predictions of RRMSEs associated to
    #ML estimates for problem 1 or 2 in a single design
  
  #INPUT :
    #des : generated design with environment covariate
      #des$xy is the table of sampling points coordinates 
        #(see tabZ in functions above)
      #des$dis is the matrix of pairwise distances between sampling points
      #des$envCent is the centered vector of environment covariate 
        #at sampling points
    #as   : positive real number, value of the aurocorrelation range
    #sig  : positive real number, residual variance
    
  #OUTPUT :
    # a vector giving theoretical predictions for :
      #RRMSE_nu_p1 : RRMSE of nu in problem 1 of maintext
      #RRMSE_gam_p2 : RRMSE of gamma in problem 2 of maintext 
      #RRMSE_as : RRMSE of a_s (in either problem of main text)
  
  N <- ncol(des$xy)
  matC <- sig^2 * exp(-des$dis/as)
  dMatC <- matC*des$dis/as^2; invC <- solve(matC);  M <- invC%*%dMatC
  IF11 <- sum(invC); IF33 <- N/(2*sig^4); IF34 <- sum(diag(M))/(2*sig^2) 
  IF44 <- 1/2*sum(diag(M%*%M))
  IF21 <- sum(invC%*%des$envCent)
  IF22 <- sum((invC%*%des$envCent)*des$envCent)
  RRMSE_nu_p1 <- sqrt(1/IF11)
  RRMSE_gam_p2 <- sqrt(1/(IF22-IF21^2/IF11))
  RRMSE_as <- sqrt(IF33/(IF33*IF44-IF34^2))/as
  return(c(RRMSE_nu_p1,RRMSE_gam_p2,RRMSE_as))}

In our analysis, we consider a gradient of \(a_s\) values:

lAsMin <- -2; lAsMax <- 2; nAs <- 30
vLAs_sim <- seq(lAsMin,lAsMax,(lAsMax-lAsMin)/(nAs-1))
vAs_sim <- 10^vLAs_sim

We build a table that makes all the combinations of design, environment and autocorrelation indices:

crossTab <- as.matrix(expand.grid(
  iDes = 1:length(lisDes),iEnv =1:nrow(tabParEnv), iAs = 1:nAs))

Computing RRMSEs under all combinations of environmental covariate and autocorrelation range values

For each combination we compute the RRMSEs :

if(is.na(match("tabRRMSEs_rev1.Rdata",list.files(rD_theoScr)))){
  tabRRMSEs <- t(apply(crossTab,1,function(vIndRow){
  des <- lisDes[[vIndRow[1]]]
  tabXY <- des$xy
  envPar <- tabParEnv[vIndRow[2],]
  vecEnv <- getEnv_rug(
    tabXY,lam=envPar[1],phi=envPar[2],ref=tabERef[,vIndRow[2]])
  vecEnvCent <- vecEnv - mean(vecEnv)
  desEnv<- list(xy=tabXY,dis=des$dis,envCent = vecEnvCent)
  as <- vAs_sim[vIndRow[3]]
  return(getTheoRRMSE(desEnv,as))}))
  save(tabRRMSEs,file=file.path(rD_theoScr,"tabRRMSEs_rev1.Rdata"))
}else load(file.path(rD_theoScr,"tabRRMSEs_rev1.Rdata"))

Computing the averaged RRMSEs across replicates of sampling designs, for each environmental covariate and autocorrelation range values

We compute the average over replicates of hybrid sampling designs. To do so, we first associate a tag referring to unique autocorrelation-environment pairs (a “context”) in the combination table crossTab.

crossTab <- cbind(crossTab,paste(crossTab[,2],crossTab[,3],sep="-"))
colnames(crossTab)[4] <- "ctxtID"; vUniqID <- unique(crossTab[,4])

Then for each context, we determine the corresponding rows of crossTab and we store them in a boolean table :

tabPosCtxtId <- sapply(vUniqID,function(ID) crossTab[,4]==ID)

A duplicate of the theable of errors tabRRMSEs is cerated where all hybrid designs error values are set to NAs :

if(is.na(match("tabRRMSEs_ave_rev1.Rdata",list.files(rD_theoScr)))){
  tabRRMSEs_ave <- tabRRMSEs
  tabRRMSEs_ave[tabParDes$typ[as.numeric(crossTab[,1])]=="mix",]<- rep(NA,3)
  
  #For each combination of type of hybrid design and context, rows
  #of *crossTab* are determined and errors are averaged out :
  
  for(i in 0:nTo){ #browsing values of p parameter
    vecIDesRep <- length(vecX)+(1:nRep)+nRep*i
    vRow_crossTab <- !is.na(match(crossTab[,1],vecIDesRep))
    for(j in 1:ncol(tabPosCtxtId)){ #browsing context
      indRep <- which(vRow_crossTab*tabPosCtxtId[,j]==1)
      tabRRMSEs_ave[indRep[1],] <- apply(tabRRMSEs[indRep,],2,mean)}}
  save(tabRRMSEs_ave,file=file.path(rD_theoScr,"tabRRMSEs_ave_rev1.Rdata"))
} else{load(file.path(rD_theoScr,"tabRRMSEs_ave_rev1.Rdata"))}
crossTab_ave <- crossTab[!is.na(tabRRMSEs_ave[,1]),1:3]
tabRRMSEs_ave <- tabRRMSEs_ave[!is.na(tabRRMSEs_ave[,1]),]

Computing average RRMSEs ranks of sampling designs across autocorrelation range values, for each environmental covariate

tabRRMSEs_ave_rk <- tabRRMSEs_ave
for(iEnv in 1:nrow(tabParEnv)){
  vIEnv <- which(as.numeric(crossTab_ave[,2])==iEnv)
  crossTab_ave_env <- crossTab_ave[vIEnv,]
  tabRRMSEs_ave_env<- tabRRMSEs_ave[vIEnv,]
  tabRRMSEs_ave_env_rk <- tabRRMSEs_ave_env
  for(iAs in 1:nAs){
    vIAs <- which(as.numeric(crossTab_ave_env[,3])==iAs)
    crossTab_ave_env_as <- crossTab_ave_env[vIAs,]
    tabRRMSEs_ave_env_as<- tabRRMSEs_ave_env[vIAs,]
    tabRRMSEs_ave_env_as_rk <- apply(tabRRMSEs_ave_env_as,2,rank)
    tabRRMSEs_ave_env_rk[vIAs,] <- tabRRMSEs_ave_env_as_rk}
  tabRRMSEs_ave_rk[vIEnv,] <- tabRRMSEs_ave_env_rk}
crossTab_aveRk <- NULL; tabRRMSEs_aveRk <-NULL
for(iEnv in 1:nrow(tabParEnv)){
  vIEnv <- which(as.numeric(crossTab_ave[,2])==iEnv)
  crossTab_ave_env <- crossTab_ave[vIEnv,]
  tabRRMSEs_ave_rk_env<- tabRRMSEs_ave_rk[vIEnv,]
  vDes <- unique(crossTab_ave_env[,1])
  tabRRMSEs_aveRk_env <- matrix(NA,nrow=length(vDes),ncol=3)
  for(iDes in vDes){
    vIDes <- which(
      (as.numeric(crossTab_ave_env[,1])==iDes)*
        (abs(log10(6*vAs_sim[as.numeric(crossTab_ave_env[,3])]/sqrt(3)))<=1)==1)
        #(abs(log10(6*vAs_sim[as.numeric(crossTab_ave_env[,3])]/sqrt(3))+1)<=1)==1)
        #(abs(log10(6*vAs_sim[as.numeric(crossTab_ave_env[,3])]/sqrt(3))-1)<=1)==1)
    tabRRMSEs_ave_rk_env_des <- tabRRMSEs_ave_rk_env[vIDes,]
    vRRMSES_aveRk_env_des <- apply(tabRRMSEs_ave_rk_env_des,2,mean)
    tabRRMSEs_aveRk_env[
      match(iDes,vDes),] <- vRRMSES_aveRk_env_des}
  crossTab_aveRk <- rbind(crossTab_aveRk,
    cbind(vDes,rep(as.character(iEnv),length(vDes))))
  tabRRMSEs_aveRk <- rbind(tabRRMSEs_aveRk,tabRRMSEs_aveRk_env)}

Graphical summary of results

vecIEnvFig <- c(2:4,8,12)[c(3,4,5,2,1)]
vecIAsFig <- c(8,10,12,14)#c(8,9,10,12)

Some preliminary visualizations

Realizations of the reponse variable under various combinations of environmental covariate and autocorrelation range values

matC_1 <- exp(-as.matrix(dist(t(tabZ_envEx))))
if(is.na(match("lisTabVecResp_rev1.Rdata",list.files(rD_theoScr)))){
  cl <- makeCluster(length(vecIEnvFig))
  clusterExport(cl,
    c("tabParEnv","getEnv_rug","tabZ_envEx","PsiTest","PhiTest","trGr",
      "vecIAsFig","vAs_sim","matC_1"))
  lisTabVecResp <- parLapply(cl,vecIEnvFig,function(iEnv){
    library(mvtnorm)
    lam <- tabParEnv[iEnv,1]; phi <- tabParEnv[iEnv,2]
    vecE_rug <- getEnv_rug(tabZ_envEx,lam,phi)
    return(sapply(vecIAsFig,function(iAs){
      as <- vAs_sim[iAs]
      matC <- matC_1^(1/as)
      vecResp <- rmvnorm(1,mean=vecE_rug,sigma=matC)
      vecResp <- exp(vecResp)/(1+exp(vecResp))
      return(vecResp)}))})
  stopCluster(cl)
  save(lisTabVecResp,file=file.path(rD_theoScr,"lisTabVecResp_rev1.Rdata"))
}else{load(file.path(rD_theoScr,"lisTabVecResp_rev1.Rdata"))}
par(
  mfrow=c(length(vecIEnvFig),length(vecIAsFig)+1),cex.lab=1.5,mar=c(5,6,4,1),
  cex.axis=1.5,lwd=2)
lislisVecResp <- sapply(vecIEnvFig,function(iEnv){
  lam <- tabParEnv[iEnv,1]; phi <- tabParEnv[iEnv,2]
  vecE_rug <- getEnv_rug(tabZ_envEx,lam,phi)
  drawFrame()
  points(tabZ_envEx[1,],tabZ_envEx[2,],pch=24,cex=0.75,
    bg=rgb(vecE_rug,0,1-vecE_rug),col=rgb(0,0,0,0))
  sapply(vecIAsFig,function(iAs){
    as <- vAs_sim[iAs]
    vecResp <- lisTabVecResp[[match(iEnv,vecIEnvFig)]][,match(iAs,vecIAsFig)]
    drawFrame()
    points(tabZ_envEx[1,],tabZ_envEx[2,],pch=24,cex=0.75,
      bg=rgb(1-vecResp,vecResp,1-vecResp),col=rgb(0,0,0,0))
    if(iEnv==vecIEnvFig[1]){
      title(paste("as =",round(vAs_sim[iAs],2),sep=" "),cex.main=2)}})})

Supplementary figure 1 of main text

nul <- sapply(1:nAs,function(iAs){
    iRow <- which(
      (as.numeric(crossTab_ave[,2])==1)*(as.numeric(crossTab_ave[,3])==iAs)==1)
    subTabRRMSEs <- tabRRMSEs_ave[iRow,]
    subCrossTab <- crossTab_ave[iRow,]
    vAsFra <- subTabRRMSEs[tabParDes$typ[as.numeric(subCrossTab[,1])]=="frac",3]
    return(vAsFra)
    })
plot(vAs_sim,nul[2,],log="xy",xlab=expression(a[s]),ylab=expression(RRMSE(a[s])))

## png 
##   2

Problem 1

Figure 4

par(mfrow=c(1,length(vecIAsFig)),cex.lab=1.5,mar=c(5,6,4,1),cex.axis=1.5,lwd=2)
nul <- sapply(vecIAsFig,function(iAs){
    iRow <- which(
      (as.numeric(crossTab_ave[,2])==1)*(as.numeric(crossTab_ave[,3])==iAs)==1)
    subTabRRMSEs <- tabRRMSEs_ave[iRow,]
    subCrossTab <- crossTab_ave[iRow,]
    vNuMix <- subTabRRMSEs[tabParDes$typ[as.numeric(subCrossTab[,1])]=="mix",1]
    vAsMix <- subTabRRMSEs[tabParDes$typ[as.numeric(subCrossTab[,1])]=="mix",3]
    vNuFra <- subTabRRMSEs[tabParDes$typ[as.numeric(subCrossTab[,1])]=="frac",1]
    vAsFra <- subTabRRMSEs[tabParDes$typ[as.numeric(subCrossTab[,1])]=="frac",3]
    tabPareto <- sapply(1:(nTo+1),function(iPar){
      vNuMixMix <- 0.99*vNuMix[iPar]<vNuMix
      vAsMixMix <- 0.99*vAsMix[iPar]<vAsMix
      vNuFraFra <- 0.99*vNuFra[iPar]<vNuFra
      vAsFraFra <- 0.99*vAsFra[iPar]<vAsFra
      vNuMixFra <- 0.99*vNuMix[iPar]<vNuFra
      vAsMixFra <- 0.99*vAsMix[iPar]<vAsFra
      vNuFraMix <- 0.99*vNuFra[iPar]<vNuMix
      vAsFraMix <- 0.99*vAsFra[iPar]<vAsMix
      return(c((prod((vNuMixMix+vAsMixMix)>0)==1),
        (prod((vNuFraFra+vAsFraFra)>0)==1),
        (prod((vNuMixFra+vAsMixFra)>0)==1),
        (prod((vNuFraMix+vAsFraMix)>0)==1)))})
    plot(vNuMix,vAsMix,type="l",col=4,lty="dotted",
         xlab=expression(RRMSE(nu)),ylab=expression(RRMSE(a[s])),
         #xlim=c(min(subTabRRMSEs[,1]),max(subTabRRMSEs[,1])),
         #ylim=c(min(subTabRRMSEs[,3]),max(subTabRRMSEs[,3]))
         xlim=c(0.2,0.7),ylim=c(0.5,1)
         )
    points(vNuMix[nTo+1],vAsMix[nTo+1],pch=2,col=4,cex=2)
    points(vNuMix[1],vAsMix[1],pch=0,col=4,cex=2)
    lines(vNuMix[tabPareto[1,]],vAsMix[tabPareto[1,]],lwd=4,col=4)
    points(vNuMix[which(tabPareto[1,]*(1-tabPareto[3,])==1)],
           vAsMix[which(tabPareto[1,]*(1-tabPareto[3,])==1)],
           pch=4,col=2,cex=2)
    lines(vNuFra,vAsFra,type="l",col=2,lty="dotted")
    points(vNuFra[1],vAsFra[1],pch=2,col=2,cex=2)
    points(vNuFra[nTo+1],vAsFra[nTo+1],pch=0,col=2,cex=2)
    lines(vNuFra[tabPareto[2,]],vAsFra[tabPareto[2,]],lwd=4,col=2)
    points(vNuFra[which(tabPareto[2,]*(1-tabPareto[4,])==1)],
       vAsFra[which(tabPareto[2,]*(1-tabPareto[4,])==1)],
       pch=4,col=4,cex=2)
    title(bquote(a[s]*"="*.(round(vAs_sim[iAs],2))),cex.main=2)})

## png 
##   2

Figure 5

bigTabPareto <- sapply(1:length(vAs_sim),function(iAs){
    iRow <- which(
      (as.numeric(crossTab_ave[,2])==1)*(as.numeric(crossTab_ave[,3])==iAs)==1)
    subTabRRMSEs <- tabRRMSEs_ave[iRow,]
    subCrossTab <- crossTab_ave[iRow,]
    vNuMix <- subTabRRMSEs[tabParDes$typ[as.numeric(subCrossTab[,1])]=="mix",1]
    vAsMix <- subTabRRMSEs[tabParDes$typ[as.numeric(subCrossTab[,1])]=="mix",3]
    vNuFra <- subTabRRMSEs[tabParDes$typ[as.numeric(subCrossTab[,1])]=="frac",1]
    vAsFra <- subTabRRMSEs[tabParDes$typ[as.numeric(subCrossTab[,1])]=="frac",3]
    tabPareto <- sapply(1:(nTo+1),function(iPar){
      vNuMixMix <- 0.99*vNuMix[iPar]<vNuMix
      vAsMixMix <- 0.99*vAsMix[iPar]<vAsMix
      vNuFraFra <- 0.99*vNuFra[iPar]<vNuFra
      vAsFraFra <- 0.99*vAsFra[iPar]<vAsFra
      vNuMixFra <- 0.99*vNuMix[iPar]<vNuFra
      vAsMixFra <- 0.99*vAsMix[iPar]<vAsFra
      vNuFraMix <- 0.99*vNuFra[iPar]<vNuMix
      vAsFraMix <- 0.99*vAsFra[iPar]<vAsMix
      return(c((prod((vNuMixMix+vAsMixMix)>0)==1),
        (prod((vNuFraFra+vAsFraFra)>0)==1),
        (prod((vNuMixFra+vAsMixFra)>0)==1),
        (prod((vNuFraMix+vAsFraMix)>0)==1)))})
    return(as.vector(t(tabPareto)))})
par(mfrow=c(1,2))
image(vecX,vAs_sim,2*bigTabPareto[84+1:28,]+bigTabPareto[28+1:28,],col=c("red","orange","red","green"),breaks=c(-0.5,0.5,1.5,2.5,3.5),log="y")
abline(h=sqrt(3)/6,lwd=2,lty="dotted")
arrows(x0=vecX[2],x1=vecX[1],y0=vAs_sim[vecIAsFig],lwd=3,length=0.1,code=2)
title("Fractal")
image(vecP,vAs_sim,2*bigTabPareto[56+1:28,]+bigTabPareto[1:28,],col=c("red","orange","red","green"),breaks=c(-0.5,0.5,1.5,2.5,3.5),log="y")
title("Hybrid")
abline(h=sqrt(3)/6,lwd=2,lty="dotted")
arrows(
  x0=vecP[1]-(vecP[2]-vecP[1])/2+10^(-3),
  x1=vecP[1]-(vecP[2]-vecP[1])/2,y0=vAs_sim[vecIAsFig],lwd=3,length=0.1,code=2)

## png 
##   2

Global performance across autocorrelation values

tabRRMSEs_aveRk1 <- tabRRMSEs_aveRk[crossTab_aveRk[,2]=="1",]
vDes <- crossTab_aveRk[crossTab_aveRk[,2]=="1",1]
vIFrac <- which(tabParDes[as.numeric(vDes),1]=="frac")
vIMix <- which(tabParDes[as.numeric(vDes),1]=="mix")
plot(tabRRMSEs_aveRk1[,1],tabRRMSEs_aveRk1[,3],col=0,xlim=c(1,56),ylim=c(1,56),
     xlab=expression(paste("Average rank of ",RRMSE(nu))),
     ylab=expression(paste("Average rank of ",RRMSE(a[s]))))
points(tabRRMSEs_aveRk1[vIFrac,1],tabRRMSEs_aveRk1[vIFrac,3],bg=2,pch=21)
points(tabRRMSEs_aveRk1[vIFrac[1],1],tabRRMSEs_aveRk1[vIFrac[1],3],col=2,pch=2,cex=2)
points(tabRRMSEs_aveRk1[vIMix,1],tabRRMSEs_aveRk1[vIMix,3],bg=4,pch=21)
points(tabRRMSEs_aveRk1[vIMix[length(vIMix)],1],tabRRMSEs_aveRk1[vIMix[length(vIMix)],3],col=4,pch=2,cex=2)

Problem 2

Figure 6

par(mfrow=c(length(vecIEnvFig),length(vecIAsFig)+1),cex.lab=1.5,mar=c(5,6,4,1),
    cex.axis=1.5,lwd=2)
for(iEnv in vecIEnvFig){
  lam <- tabParEnv[iEnv,1]; phi <- tabParEnv[iEnv,2]
  vecE_rug <- getEnv_rug(tabZ_envEx,lam,phi)
  if (match(iEnv,vecIEnvFig)<=2) ltyBox <- 1 else ltyBox <- 2
  drawFrame()
  points(tabZ_envEx[1,],tabZ_envEx[2,],pch=24,cex=0.75,
    bg=rgb(vecE_rug,0,1-vecE_rug),col=rgb(0,0,0,0))
  sapply(vecIAsFig,function(iAs){
    iRow <- which(
      (as.numeric(crossTab_ave[,2])==iEnv)
      *(as.numeric(crossTab_ave[,3])==iAs)==1)
    subTabRRMSEs <- tabRRMSEs_ave[iRow,]
    subCrossTab <- crossTab_ave[iRow,]
    vGaMix <- subTabRRMSEs[tabParDes$typ[as.numeric(subCrossTab[,1])]=="mix",2]
    vAsMix <- subTabRRMSEs[tabParDes$typ[as.numeric(subCrossTab[,1])]=="mix",3]
    vGaFra <- subTabRRMSEs[tabParDes$typ[as.numeric(subCrossTab[,1])]=="frac",2]
    vAsFra <- subTabRRMSEs[tabParDes$typ[as.numeric(subCrossTab[,1])]=="frac",3]
    tabPareto <- sapply(1:(nTo+1),function(iPar){
      vGaMixMix <- 0.99*vGaMix[iPar]<vGaMix
      vAsMixMix <- 0.99*vAsMix[iPar]<vAsMix
      vGaFraFra <- 0.99*vGaFra[iPar]<vGaFra
      vAsFraFra <- 0.99*vAsFra[iPar]<vAsFra
      vGaMixFra <- 0.99*vGaMix[iPar]<vGaFra
      vAsMixFra <- 0.99*vAsMix[iPar]<vAsFra
      vGaFraMix <- 0.99*vGaFra[iPar]<vGaMix
      vAsFraMix <- 0.99*vAsFra[iPar]<vAsMix
      return(c((prod((vGaMixMix+vAsMixMix)>0)==1),
        (prod((vGaFraFra+vAsFraFra)>0)==1),
        (prod((vGaMixFra+vAsMixFra)>0)==1),
        (prod((vGaFraMix+vAsFraMix)>0)==1)))})
    box(lty=ltyBox,lwd=8)
    plot(vGaMix,vAsMix,type="l",col=4,lty="dotted",
         xlab=expression(RRMSE(gamma)),ylab=expression(RRMSE(a[s])),
         #xlim=c(min(subTabRRMSEs[,2]),max(subTabRRMSEs[,2])),
         #ylim=c(min(subTabRRMSEs[,3]),max(subTabRRMSEs[,3]))
         xlim=c(0.3,2),ylim=c(0.5,1.2)
         )
    points(vGaMix[nTo+1],vAsMix[nTo+1],pch=2,col=4,cex=2)
    points(vGaMix[1],vAsMix[1],pch=0,col=4,cex=2)
    lines(vGaMix[tabPareto[1,]],vAsMix[tabPareto[1,]],lwd=4,col=4)
    points(vGaMix[which(tabPareto[1,]*(1-tabPareto[3,])==1)],
           vAsMix[which(tabPareto[1,]*(1-tabPareto[3,])==1)],
           pch=4,col=2,cex=2)
    lines(vGaFra,vAsFra,type="l",col=2,lty="dotted")
    points(vGaFra[1],vAsFra[1],pch=2,col=2,cex=2)
    points(vGaFra[nTo+1],vAsFra[nTo+1],pch=0,col=2,cex=2)
    lines(vGaFra[tabPareto[2,]],vAsFra[tabPareto[2,]],lwd=4,col=2)
    points(vGaFra[which(tabPareto[2,]*(1-tabPareto[4,])==1)],
       vAsFra[which(tabPareto[2,]*(1-tabPareto[4,])==1)],
       pch=4,col=4,cex=2)
    if(iEnv==vecIEnvFig[1]){
      title(bquote(a[s]*"="*.(round(vAs_sim[iAs],2))),cex.main=2)}})
    box(lty=ltyBox,lwd=8)}

## png 
##   2

Figure 7

lisBigTabPareto <- lapply(vecIEnvFig[c(1,4)],function(iEnv){
  tabPareto <- sapply(1:length(vAs_sim),function(iAs){
    iRow <- which(
      (as.numeric(crossTab_ave[,2])==iEnv)
      *(as.numeric(crossTab_ave[,3])==iAs)==1)
    subTabRRMSEs <- tabRRMSEs_ave[iRow,]
    subCrossTab <- crossTab_ave[iRow,]
    vGaMix <- subTabRRMSEs[tabParDes$typ[as.numeric(subCrossTab[,1])]=="mix",2]
    vAsMix <- subTabRRMSEs[tabParDes$typ[as.numeric(subCrossTab[,1])]=="mix",3]
    vGaFra <- subTabRRMSEs[tabParDes$typ[as.numeric(subCrossTab[,1])]=="frac",2]
    vAsFra <- subTabRRMSEs[tabParDes$typ[as.numeric(subCrossTab[,1])]=="frac",3]
    bigTabPareto <- sapply(1:(nTo+1),function(iPar){
      vGaMixMix <- 0.99*vGaMix[iPar]<vGaMix
      vAsMixMix <- 0.99*vAsMix[iPar]<vAsMix
      vGaFraFra <- 0.99*vGaFra[iPar]<vGaFra
      vAsFraFra <- 0.99*vAsFra[iPar]<vAsFra
      vGaMixFra <- 0.99*vGaMix[iPar]<vGaFra
      vAsMixFra <- 0.99*vAsMix[iPar]<vAsFra
      vGaFraMix <- 0.99*vGaFra[iPar]<vGaMix
      vAsFraMix <- 0.99*vAsFra[iPar]<vAsMix
      return(c((prod((vGaMixMix+vAsMixMix)>0)==1),
        (prod((vGaFraFra+vAsFraFra)>0)==1),
        (prod((vGaMixFra+vAsMixFra)>0)==1),
        (prod((vGaFraMix+vAsFraMix)>0)==1)))})
    return(as.vector(t(bigTabPareto)))})})
par(mfrow=c(2,2))
lapply(1:2,function(iBigTabPareto){
    bigTabPareto <-lisBigTabPareto[[iBigTabPareto]] 
    image(vecX,vAs_sim,2*bigTabPareto[84+1:28,]+bigTabPareto[28+1:28,],col=c("red","orange","red","green"),breaks=c(-0.5,0.5,1.5,2.5,3.5),log="y")
  abline(h=sqrt(3)/6,lwd=2,lty="dotted")
  arrows(x0=vecX[2],x1=vecX[1],y0=vAs_sim[vecIAsFig],lwd=3,length=0.1,code=2)
  if(iBigTabPareto==2) title("Fractal - Non-Monotonic") else{
    title("Fractal - Monotonic")} 
  image(vecP,vAs_sim,2*bigTabPareto[56+1:28,]+bigTabPareto[1:28,],col=c("red","orange","red","green"),breaks=c(-0.5,0.5,1.5,2.5,3.5),log="y")
  if(iBigTabPareto==2) title("Hybrid - Non-Monotonic") else{
    title("Hybrid - Monotonic")} 
  abline(h=sqrt(3)/6,lwd=2,lty="dotted")
    arrows(
    x0=vecP[1]-(vecP[2]-vecP[1])/2+10^(-3),
    x1=vecP[1]-(vecP[2]-vecP[1])/2,
    y0=vAs_sim[vecIAsFig],lwd=3,length=0.1,code=2)})

## [[1]]
## NULL
## 
## [[2]]
## NULL
## [[1]]
## NULL
## 
## [[2]]
## NULL
## png 
##   2

Global performance across autocorrelation values

par(mfrow=c(3,3))
nul <- sapply(which(tabParEnv[,1]>vecLam[1]),function(iEnv){
  tabRRMSEs_aveRkEnv <- tabRRMSEs_aveRk[crossTab_aveRk[,2]==as.character(iEnv),]
  vDes <- crossTab_aveRk[crossTab_aveRk[,2]==as.character(iEnv),1]
  vIFrac <- which(tabParDes[as.numeric(vDes),1]=="frac")
  vIMix <- which(tabParDes[as.numeric(vDes),1]=="mix")
  plot(tabRRMSEs_aveRkEnv[,2],tabRRMSEs_aveRkEnv[,3],col=0,xlim=c(1,56),ylim=c(1,56),
       xlab=expression(paste("Average rank of ",RRMSE(gamma))),
       ylab=expression(paste("Average rank of ",RRMSE(a[s]))))
  points(tabRRMSEs_aveRkEnv[vIFrac,2],tabRRMSEs_aveRkEnv[vIFrac,3],bg=2,pch=21)
  lines(tabRRMSEs_aveRkEnv[vIFrac,2],tabRRMSEs_aveRkEnv[vIFrac,3],col=2)
  points(tabRRMSEs_aveRkEnv[vIFrac[1],2],tabRRMSEs_aveRkEnv[vIFrac[1],3],col=2,pch=2,cex=2)
  points(tabRRMSEs_aveRkEnv[vIFrac[nTo+1],2],tabRRMSEs_aveRkEnv[vIFrac[nTo+1],3],col=2,pch=0,cex=2)
  points(tabRRMSEs_aveRkEnv[vIMix,2],tabRRMSEs_aveRkEnv[vIMix,3],bg=4,pch=21)
  lines(tabRRMSEs_aveRkEnv[vIMix,2],tabRRMSEs_aveRkEnv[vIMix,3],col=4)
  points(tabRRMSEs_aveRkEnv[vIMix[length(vIMix)],2],tabRRMSEs_aveRkEnv[vIMix[length(vIMix)],3],col=4,pch=2,cex=2)
  points(tabRRMSEs_aveRkEnv[vIMix[1],2],tabRRMSEs_aveRkEnv[vIMix[1],3],col=4,pch=0,cex=2)
  title(main=tabParEnv[iEnv,])})

Figure 8 - global performance for problems 1 & 2

vPareto1 <- sapply(1:56,function(i){
  r1 <- tabRRMSEs_aveRk1[i,1]; r2 <- tabRRMSEs_aveRk1[i,3]
  return(prod(
    ((0.99*r1<tabRRMSEs_aveRk1[,1])+(0.99*r2<tabRRMSEs_aveRk1[,3]))>0))})
png(file.path(rD_theoScr,"Figure8.png"),width=1200,height=400)
par(mfrow=c(1,3),mar=c(5,8,4,1))
plot(tabRRMSEs_aveRk1[,1],tabRRMSEs_aveRk1[,3],col=0,xlim=c(1,56),ylim=c(1,56),
     xlab="",
     ylab="",
     cex.axis=2.5,cex.lab=2.5)
mtext(expression(paste("Average rank of ",RRMSE(nu))),1,padj=1.4,adj=0.5,cex=1.8)
mtext(expression(paste("Average rank of ",RRMSE(a[s]))),2,padj=-1,adj=0.5,cex=1.8)
points(tabRRMSEs_aveRk1[vIFrac,1],tabRRMSEs_aveRk1[vIFrac,3],col=2,bg=2,pch=1+20*vPareto1[vIFrac],cex=2.5)
lines(tabRRMSEs_aveRk1[vIFrac,1],tabRRMSEs_aveRk1[vIFrac,3],col=2,lwd=1)
points(tabRRMSEs_aveRk1[vIFrac[1],1],tabRRMSEs_aveRk1[vIFrac[1],3],col=2,pch=2,cex=4)
points(tabRRMSEs_aveRk1[vIFrac[nTo+1],1],tabRRMSEs_aveRk1[vIFrac[nTo+1],3],col=2,pch=0,cex=4)
points(tabRRMSEs_aveRk1[vIMix,1],tabRRMSEs_aveRk1[vIMix,3],col=4,bg=4,pch=1+20*vPareto1[vIMix],cex=2.5)
lines(tabRRMSEs_aveRk1[vIMix,1],tabRRMSEs_aveRk1[vIMix,3],col=4,lwd=1)
points(tabRRMSEs_aveRk1[vIMix[length(vIMix)],1],tabRRMSEs_aveRk1[vIMix[length(vIMix)],3],col=4,pch=2,cex=4)
points(tabRRMSEs_aveRk1[vIMix[1],1],tabRRMSEs_aveRk1[vIMix[1],3],col=4,pch=0,cex=4)
title(main="Problem 1",cex.main=2.5)
nul <- sapply(vecIEnvFig[c(1,4)],function(iEnv){
  tabRRMSEs_aveRkEnv <- tabRRMSEs_aveRk[crossTab_aveRk[,2]==as.character(iEnv),]
  vPareto2 <- sapply(1:56,function(i){
  r1 <- tabRRMSEs_aveRkEnv[i,1]; r2 <- tabRRMSEs_aveRkEnv[i,3]
  return(prod(
    ((0.99*r1<tabRRMSEs_aveRkEnv[,1])+(0.99*r2<tabRRMSEs_aveRkEnv[,3]))>0))})

  vDes <- crossTab_aveRk[crossTab_aveRk[,2]==as.character(iEnv),1]
  vIFrac <- which(tabParDes[as.numeric(vDes),1]=="frac")
  vIMix <- which(tabParDes[as.numeric(vDes),1]=="mix")
  plot(tabRRMSEs_aveRkEnv[,2],tabRRMSEs_aveRkEnv[,3],col=0,xlim=c(1,56),ylim=c(1,56),
       xlab="",
       ylab="",
      cex.axis=2.5,cex.lab=2.5)
  mtext(expression(paste("Average rank of ",RRMSE(gamma))),1,padj=1.4,adj=0.5,cex=1.8)
  mtext(expression(paste("Average rank of ",RRMSE(a[s]))),2,padj=-1,adj=0.5,cex=1.8)
  points(tabRRMSEs_aveRkEnv[vIFrac,2],tabRRMSEs_aveRkEnv[vIFrac,3],col=2,bg=2,pch=1+20*vPareto2[vIFrac],cex=2.5)
  lines(tabRRMSEs_aveRkEnv[vIFrac,2],tabRRMSEs_aveRkEnv[vIFrac,3],col=2)
  points(tabRRMSEs_aveRkEnv[vIFrac[1],2],tabRRMSEs_aveRkEnv[vIFrac[1],3],col=2,pch=2,cex=4)
  points(tabRRMSEs_aveRkEnv[vIFrac[nTo+1],2],tabRRMSEs_aveRkEnv[vIFrac[nTo+1],3],col=2,pch=0,cex=4)
  points(tabRRMSEs_aveRkEnv[vIMix,2],tabRRMSEs_aveRkEnv[vIMix,3],col=4,bg=4,pch=1+20*vPareto2[vIMix],cex=2.5)
  lines(tabRRMSEs_aveRkEnv[vIMix,2],tabRRMSEs_aveRkEnv[vIMix,3],col=4)
  points(tabRRMSEs_aveRkEnv[vIMix[length(vIMix)],2],tabRRMSEs_aveRkEnv[vIMix[length(vIMix)],3],col=4,pch=2,cex=4)
  points(tabRRMSEs_aveRkEnv[vIMix[1],2],tabRRMSEs_aveRkEnv[vIMix[1],3],col=4,pch=0,cex=4)
  if(iEnv==vecIEnvFig[1]) typ <- "Monotonic" else typ <- "Non-monotonic"
  title(main=paste("Problem 2 - ",typ,sep=""),cex.main=2.5)
})
dev.off()
## png 
##   2

Shortest spanning paths of grid and fractal designs

We use the theoretical results of the appendix to build functions computing shortest spanning paths of fractal design and grid design :

SP_frac <- function(N,L,x){
  n <- log(N)/log(3)
  rho <- sqrt(3)/(sqrt(3)+2)*x
  return(L*(2*(1-2*rho)*(1-(3*rho)^n)/(1-3*rho)+rho^n*(3^n-1)))}
SP_grid <- function(N,L){
  return(L*2*(N-1)/(sqrt(8*N+1)-3))}

We compare this length over a range of design sizes :

L <- sqrt(3); xVal <- (1/3/sqrt(3)*(2+sqrt(3)))
vX <- 10^(seq(-1.5,0,0.01))
plot(vX,sapply(vX,function(x) SP_frac(27,L,x)),log="x",ylim=c(0,SP_grid (28,L)),xlab="x",ylab="Spanning path length",pch=16-15*(vX>xVal))
abline(h=SP_grid (28,L),col=2,lwd=2)

## png 
##   2