My supervisors and I are planning to run an ecological experiment in the very near future. They have proposed an experimental design and I would like to get some feedback on it. In particular, I am curious as to whether I will be able to use this design to calculate a cubic response surface.

We have two independent variables (named NutrientP, and Movement) and one response variable. The experiment will be replicated once. Below is how we plan on sampling the design surface (note the four values on NutrientP axis are set up in an exponential way, plus zero control).

Five points on both diagonals, plus controls on each axis for the middle point

The points were placed in a way that samples the area where we expect to have the greatest effect (near the middle). However, in the case that a quadratic fit will not suffice, I wonder if a cubic fit will be possible. This is an experiment that will last multiple months, meaning that it will not be possible to augment the design if we find that we are missing certain points to obtain a better fit.

  • $\begingroup$ Is it correct that you are referring to orders of reaction when you use the term third-order response? This is a distinct consideration from whether to fit a dynamical system with some equation order. Or are you referring to something else? $\endgroup$
    – Galen
    Jun 17, 2022 at 16:58
  • 2
    $\begingroup$ Did you intend to leave out a point near $(40, 80)$? Regardless, you show 11 points and it takes only 10 to determine a cubic function; and these are sufficiently spread out that they won't be redundant. But that leaves only one degree of freedom to estimate the error term, which hugely increases the uncertainty. Consider adding a few more points or replicating the experiment at some of them. BTW, you ought to answer this question yourself by creating a realistic random dataset and applying your intended data analysis procedure to it--before you start the experiment. $\endgroup$
    – whuber
    Jun 17, 2022 at 17:11
  • 1
    $\begingroup$ @ekatko1 Okay, alternatively called the degree of a polynomial. $\endgroup$
    – Galen
    Jun 17, 2022 at 17:23
  • 2
    $\begingroup$ @whube yes, I think that was intentional on behalf of my supervisors, but I am also wondering if should be included, along with the one at $(20, 340)$. Also, I forgot to mention that we were planning to replicate each point. Good idea, for the simulation data, I will do exactly this and perhaps post it as an answer. $\endgroup$
    – ekatko1
    Jun 17, 2022 at 17:44
  • 1
    $\begingroup$ There's a lot unsaid here. The theory of experimental design provides guidance for choosing the number of points, spacing them, selecting an appropriate form for the response surface, and making a replication strategy. At this stage it's worth questioning all those choices, because developing an optimal answer can save substantial effort and cost (or improve the results for a given budget). $\endgroup$
    – whuber
    Jun 17, 2022 at 18:13

1 Answer 1


Letting $x_1$ and $x_2$ be the planar coordinates, the most general cubic response has the form

$$\begin{aligned} z & = \beta_0 \\ & + \sum_{i=1}^2 \beta_i x_i\\ & + \sum_{1\le i \le j \le 2} \beta_{ij}x_ix_j\\ & + \sum_{1\le i\le j \le k \le 2} \beta_{ijk}x_ix_jx_k. \end{aligned}$$

The number of coefficients is $1 + 2 + 3 + 4 = 10.$ Although you propose $11$ locations, which ought to be enough to estimate all ten parameters, it's worth a check because one or more of those locations might be redundant for these purposes.

The way to check is to consider the design matrix for a single observation per location. If it is of full rank, the design matrix with multiple observations will be of full rank, too; and if it is not of full rank, replicating observations will not cure that.

Although we could do this check manually, it's so easy to get your statistical software to help that we might as well just put it to work. Begin with an array of coordinates of those 11 points:

df <- data.frame(x1 = c(0,0,0, 10,10, 20,20, 30,30, 40,40),
                x2 = c(0,80,320, 40,160, 0,80, 40,160, 0,320))

Now pass it to your ordinary least squares routine, using any response vector you like. If it returns estimates for all $10$ variables, you're good. Here I used a vector of zeros for the response:

lm(rep(0, nrow(df)) ~ x1+x2 + x1*x2+I(x1^2)+I(x2^2) + I(x1^3)+I(x1^2*x2)+I(x1*x2^2)+I(x2^3), df)
Coefficients:  (Intercept)  x1  x2  I(x1^2) I(x2^2) I(x1^3) I(x1^2 * x2) I(x1 * x2^2) I(x2^3) x1:x2  
                         0   0   0        0       0       0            0            0       0     0

Ten coefficients were estimated: your design is identifiable.

In R, you can avoid computing and writing out all the terms by using its command to form polynomials, as in

X <- cbind(1, poly(as.matrix(df), 3, raw = TRUE))

This is a matrix with ten columns, one for each variable (including the intercept). A good way to analyze it is to compute its singular values:

print(svd(X)$d, digits=2)
[1] 4.7e+07 3.0e+06 7.8e+04 5.1e+04 1.7e+04 6.0e+02 2.2e+02 4.6e+01 2.9e+00 9.9e-01

They are all nonzero--the smallest, at the end, is 0.099--but there is a huge range. That's a sign you might have difficulties estimating some coefficients (or combinations) of those coefficients accurately.

R (and many other platforms) offers a way to improve the numerical stability of the model through orthogonal polynomials. The correct analysis, as Belsley, Kuh, and Welsch demonstrate in Regression Diagnostics, is to analyze the centered variables (leaving out the intercept):

X <- cbind(1, poly(as.matrix(df), 3, raw = FALSE))
print(svd(scale(X)[,-1])$d, digits=2)
 [1] 4.63 4.10 3.65 3.33 3.12 2.89 2.72 1.33 0.25

That's much better: the range is $4.63:0.25 = 18.7.$ This is the condition number of the design matrix. It's a little large: values exceeding $5$ in an experiment suggest room for improvement. As an example of what can be accomplished, this design has a much lower condition number

df <- data.frame(x1 = c(0,0,0, 10,10, 20,20, 30,30, 40,40),
                 x2 = c(0, 160, 320, 80, 240, 0, 320, 80, 240, 0, 320))

Figure showing the point locations for this design as a scatterplot.

X <- poly(as.matrix(df), 3, raw = FALSE)
print(d <- svd(scale(X))$d, digits=2)
[1] 4.1 3.7 3.6 3.3 3.2 3.0 2.9 2.5 1.5

The singular values now run from $4.1$ down to $1.5$ for a condition number of just $2.8$. If this design meets your other criteria and is feasible to carry out, it will result in substantially more precise estimates of some coefficients for the same experimental cost (assuming this least squares model of the experimental variation is reasonable).

I am not proposing this particular redesign, but only pointing out how one can easily go about evaluating possibilities before going to the expense and trouble of doing the experiment.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.