A simple 2^2 Full Factorial Design of Experiments (DOE): JMP vs Python

In this article we are going to review one of the simplest DOE configurations: the 2^2 full factorial design. We will see how to replicate the JMP results in Python so we can have a deeper understanding of what is going on behind the scenes.

In a 2x2 full factorial DOE we have 2 factors (e.g., X1 and X2) and each of the factors is set at two levels (e.g., -1, 1).

A full factorial design explores all the combinations of factors and levels. In our example the combinations of 2 factors with two levels are 2^2 = 4. This define a design space that can be visualized as a box in the X1-X2 space. The corners represent the four possible combinations:

The JMP way

In JMP we can use the Full factorial design under the Classical DOE interface:

For simplicity we’ll not use center points and replicates, so we have 4 experiments to run and collect the results (response Y). Here below the result with simulated Ys:

We can now fit our model. JMP automatically propose a model, associated to the DOE configuration, which will estimate the 2 main effects (X1 and X2) and their interaction (X1*X2)

This corresponds to find the betas in the following linear equation, using the Least Squares method:

?Here below the result:

The column “Estimate” reports the values for our betas. This is our model:

We have an interaction term in our model. JMP visualize this interaction in the Interaction Profiles:

The interaction profile is a simple plot that shows the behavior of the model for one variable when the other is fixed. For example, the upper-right graph shows the behavior of the model for X2 when X1 = -1 (red line) and X1 = 1 (blue line). Visually the fact that red and blue lines are not parallel indicates that there in an interaction. Mathematically this corresponds to have a X1*X2 term in the model which change the parallelism expected if only the X1 + X2 terms were present.

The Python way

Let’s now try to code all the above results in Python. We first need to load some libraries:

We can save the JMP table in a .csv file and load back in a Pandas dataframe (df):

We can now fit our model using the statsmodels OLS (Ordinary Least Square) library.The interaction term is not included by default, and we have to add it passing a formula:

The R-style syntax ‘Y ~ X1*X2’ tells to add the interaction term X1*X2 (including the individual columns too) in the model. We can now fit our model:

In the “coef” column we have the same betas estimated by JMP.

The interaction profiler

As said before we have to fix one variable and move the other. These are four 2x2 combinations matrices that we can obtain by simply filtering the dataframe columns X1 and X2:

We now need the predicted Ys vectors associated to the four X matrices above:

Let’s plot these X and Y:

And we have obtained the same JMP interaction profiles.

The Prediction Profiler

Another nice feature in JMP is the model Prediction Profiler. Here below the JMP profiler for our simple model:

The prediction profiler is a cross-section of the model at the points of the vertical red lines. For example, if we fix X2=0, the Y will vary according to the X1 graph when X1 goes from -1 to 1. Vice versa if we fix X1=0, the Y will vary accordingly with the X2 graph when X2 moves from -1 to 1.

To replicate this in Python let’s define 2 matrices with X1 and X2 set at these zero cross-section points:

Let’s predict now the two Y vectors associated:

We can now plot the X and Y:

Here we are. We obtained the JMP profiler in this specific 0,0 point. Of course, this is just a static picture and has not the nice JMP interaction capability, but this was out of scope :)