Consider a simple linear model \(y(x)=\beta_0+\beta_1 x+\varepsilon(x)\) and \(x\in\mathcal{X}=\{1,2,3,4,5,6,7,8\}\). We will compare the following four exact designs of size \(N=24\):

\[ \xi_1=\left\{ \begin{array}{cccccccc} 1&2&3&4&5&6&7&8\\ 3&3&3&3&3&3&3&3 \end{array} \right\} \] \[ \xi_2=\left\{ \begin{array}{cccccccc} 1&2&3&4&5&6&7&8\\ 6&0&0&6&6&0&0&6 \end{array} \right\} \] \[ \xi_2=\left\{ \begin{array}{cccccccc} 1&2&3&4&5&6&7&8\\ 0&0&0&12&12&0&0&0 \end{array} \right\} \] \[ \xi_4=\left\{ \begin{array}{cccccccc} 1&2&3&4&5&6&7&8\\ 12&0&0&0&0&0&0&12 \end{array} \right\} \]

w1 <- rep(3,8)
w2 <- c(6,0,0,6,6,0,0,6)
w3 <- c(0,0,0,12,12,0,0,0)
w4 <- c(12,0,0,0,0,0,0,12)

Suppose that we know that the true values of the paramaters \(\beta=(\beta_0,\beta_1)'\) are \((9,8)'\).

truepar <- c(9,8)

Using the package OptimalDesign, we can easily construct the matrix of candidate regressors:

library(OptimalDesign) 
F <- Fx_cube(~x1,1,8,8)
## [1] Call of the function:
## Fx_cube(formula = ~x1, lower = 1, upper = 8, n.levels = 8)

Now, let’s simulate all four scenarios: make the required measurements as prescribed by each of the design, and then estimate the unknown parameters (this is done in the function Experiment). Do this 5000 times for each design and visualize the final estimates with respect to the true value of \(\beta\) (this is done in SimExp2).

Experiment <- function(theta,F,w){
    y <- rep(0,24)
    Ft <- matrix(0,ncol=2,nrow=24)
    u <- 1
    for (i in 1:8){
      if (w[i]!=0){
        for (j in 1:w[i]){
          y[u] <- F[i,]%*%theta+rnorm(1)
          Ft[u,] <- F[i,]
          u <- u+1
        }
      }
    }   
    lse <- lsfit(Ft,y,intercept=FALSE)
    lse$coefficients
}

SimExp2 <- function(){
    d1 <- matrix(0,ncol=2,nrow=5000)
    d2 <- matrix(0,ncol=2,nrow=5000)
    d3 <- matrix(0,ncol=2,nrow=5000)
    d4 <- matrix(0,ncol=2,nrow=5000)
    for ( i in 1:5000){
        d1[i,] <- Experiment(truepar,F,w1)
        d2[i,] <- Experiment(truepar,F,w2)
        d3[i,] <- Experiment(truepar,F,w3)
        d4[i,] <- Experiment(truepar,F,w4)
    }
    plot(d1,xlim=c(6,12), ylim=c(7.5,8.5), col="yellow", pch=".", cex=2)
    points(d2, col="blue", pch=".", cex=2)
    points(d3, col="green", pch=".", cex=2)
    points(d4, col="red", pch=".", cex=2)
    points(x=truepar[1], y=truepar[2], col="black", pch=19)
    
}

SimExp2()

We can see that the resulting scatters of points resemble ellipses. The larger the scatter, the larger the variance of the least squares estimates \(\hat{\beta}\), and, hence, the worse the initial design. If we shift these ellipses so that they are centered in the origin, we get so-called concentration ellipses.

In general, we can define concentration ellipsoid of a nonsingular design \(\xi\) as \[ \mathcal{E}_{\xi}=\{z: z'M(\xi)z\leq 1\}. \]

If there are two designs \(\xi_1\) and \(\xi_2\) for which \(\mathcal{E}_{\xi_1}\subset \mathcal{E}_{\xi_2}\), we say that the design \(\xi_1\) is uniformly better than the design \(\xi_2\).

DrawEl <- function(M, cl="red"){
  eigens <- eigen(M)
  evs <- sqrt(eigens$values)
  evecs <- eigens$vectors

  alpha <- atan(evecs[ , 1][2] / evecs[ , 1][1])
  theta <- seq(0, 2 * pi, length=(1000))

  x <- evs[1] * cos(theta) * cos(alpha) - evs[2] * sin(theta) * sin(alpha)
  y <- evs[1] * cos(theta) * sin(alpha) + evs[2] * sin(theta) * cos(alpha)

  plot(x, y, type = "l", lwd=2, col = cl, xlim=c(-2,2), ylim=c(-0.5,0.5))
}

DrawEl(solve(infmat(F, w2, echo=FALSE)), cl="blue")
par(new=TRUE)
DrawEl(solve(infmat(F, w1, echo=FALSE)), cl="yellow")
par(new=TRUE)
DrawEl(solve(infmat(F, w3, echo=FALSE)), cl="green")
par(new=TRUE)
DrawEl(solve(infmat(F, w4, echo=FALSE)), cl="red")

Hence, as the picture suggest, the design \(\xi_4\) is uniformly better than all of the remaining design, while the design \(\xi_2\) is uniformly worst of all four designs.

Example: Quadratic Regression

Let’s visualize the concentration ellipsoids in a quadratic regression given by \[ y(x)=\beta_0+\beta_1 x+\beta_2 x^2+\varepsilon(x), \] where \(x\in\mathcal{X}=[-1,1]\).

Consider these three approximate designs: \[ \xi_1=\left\{ \begin{array}{ccc} -1&\cdots&1\\ \frac{1}{111} & \cdots & \frac{1}{111} \end{array} \right\} \] \[ \xi_2=\left\{ \begin{array}{ccc} -1&0&1\\ \frac{1}{3}&\frac{1}{3}&\frac{1}{3} \end{array} \right\} \] \[ \xi_3=\left\{ \begin{array}{ccc} -1&0&1\\ \frac{1}{4}&\frac{1}{2}&\frac{1}{4} \end{array} \right\} \]

F2 <- Fx_cube(~x1+I(x1^2),-1,1,111)
## [1] Call of the function:
## Fx_cube(formula = ~x1 + I(x1^2), lower = -1, upper = 1, n.levels = 111)
wq <- rep(0,111)
wq1 <- wq
wq1[1] <- 1/3; wq1[56] <- 1/3; wq1[111] <- 1/3
wq2 <- rep(1/111, 111)
wq3 <- wq
wq3[1] <- 1/4; wq3[56] <- 1/2; wq3[111] <- 1/4

For drawing the ellipsoids using the function ellipsoid3d available in the vignette for the package `rgl’, we need to compute the eigenvalues of the information matrices of the designs \(\xi_1\), \(\xi_2\) and \(\xi_3\):

X1 <- eigen(solve(infmat(F2, wq1, echo=FALSE)))$values
X2 <- eigen(solve(infmat(F2, wq2, echo=FALSE)))$values
X3 <- eigen(solve(infmat(F2, wq3, echo=FALSE)))$values
bg3d("white")
wire3d(ellipsoid3d(X1[1],X1[2],X1[3],qmesh=TRUE),col='green')
wire3d(ellipsoid3d(X2[1],X2[2],X2[3],qmesh=TRUE),col='blue')
wire3d(ellipsoid3d(X3[1],X3[2],X3[3],qmesh=TRUE),col='magenta')
grid3d(c("x", "y", "z"))

You must enable Javascript to view this page properly.

We can see that in this case, none of the designs \(\xi_2\) and \(\xi_3\) is uniformly better than the other one. In this case, the quality of the design can be measured with respect to an optimality criterion \(\Phi:\mathbb{S}^m_+\to[0,\infty)\), which is a function of the information matrix of the design (more on optimality criteria in a different vignette).

We can easily verify that the design \(\xi_2\) is better with respect to the \(D\)-optimality criterion, whereas the design \(\xi_3\) is better with respect to the \(A\)-optimality criterion.

optcrit(F2, wq1, crit="D", echo=FALSE)
## [1] 0.5291337
optcrit(F2, wq1, crit="A", echo=FALSE)
## [1] 0.3333333
optcrit(F2, wq2, crit="D", echo=FALSE)
## [1] 0.3150398
optcrit(F2, wq2, crit="A", echo=FALSE)
## [1] 0.1869022
optcrit(F2, wq3, crit="D", echo=FALSE)
## [1] 0.5
optcrit(F2, wq3, crit="A", echo=FALSE)
## [1] 0.375