How to construct an interaction plot

Question

I am debating how to construct an interaction plot with my supervisor.

We have a dataset comprising 8 independent variables. We are trying to analyse the effect of 2 of the 8 independent variables on the dependent variable.

My supervisor is suggesting that to draw an interaction plot, we first fit a full model using all of the 8 independent variables we have, and for every possible combination of the levels of the 2 independent variables that we are particularly interested in (let's call them var1 and var2), calculate their predicted value based on the full model that we constructed earlier. However when applying this method I was running into a problem because in order for my statistical software to make a prediction, I had to assign values for the 8-2=6 variables that are left in the dataset, which are undetermined. So I suggested to my supervisor that instead of relying on the full model for calculating predicted values, I fit a model like the one below:

y = var1 + var2 + var1*var2

(i.e. instead of y = var1 + var2+ var3 + var4 + var5 + var6 + var7 + var8 + var1*var2)

My supervisor, however, disagrees with my view and is telling me to go on by using the mean values of var3, var4, var5, var6, var7, var8, which I can calculate from our original dataset, to come up with predictions.

Is there something wrong with my method of analysing the effect of interaction? I prefer my method because the interaction plot looks much better with my method. However, if my method is theoretically wrong then I guess I have to stick with what is said by my supervisor.

Answer 1

I'm not sure that I completely understand your supervisor's suggestion, but the principle that I use when choosing how to create a graph is the make sure that the graph represents the analysis that I'm reporting in my paper. Based on this principle, I would use whatever model to create your graph that you're reporting in your paper. Thus, if you are reporting the following model:

$y = var1 + var2 + var1 * var2$

then I would use this model to obtain the predicted values that you plot on your graph. On the other hand, if you are reporting the following model:

$y = var1 + var2 + var3 + var4 + var5 + var6 + var7 + var8 + var9 + var1 * var2$

then I would plot the $var1 * var2$ interaction from this model, mean-centering var3 through var9 when you obtain the predicted values for your graph.

Assuming that the model with your control variables is the one that you are reporting in your paper, I have included some R code simulating data and creating a graph using those data below. You may want to consider plotting your $y$ points marginalized for your various control variables; if you do not know how to do this, I describe how to accomplish this here.

# Set the seed
set.seed(2314)

# Create the data
dat <- matrix(NA, nrow = 200, ncol = 9)
colnames(dat) <- paste0("var", 1:9)
dat <- data.frame(dat)
for(i in 1:9)
{
  dat[, paste0("var", i)] <- rnorm(200, sd = 1)
}
dat$y <- .5 * dat$var1 + .5 * dat$var2 + .5 * dat$var1 * dat$var2 + rnorm(200, sd = 1)

# Fit the model
mod <- lm(y ~ var1 * var2 + var3 + var4 + var5 + var6 + var7 + var8 + var9, data = dat)

# Create a matrix of desired predicted values for the model.  I am holding the control variables
# constant at their means
pX <- expand.grid(var1 = seq(min(dat$var1), max(dat$var1), by = .1), 
                  var2 = c(mean(dat$var2) - sd(dat$var2), mean(dat$var2) + sd(dat$var2)),
                  var3 = mean(dat$var3),
              var4 = mean(dat$var4),
                  var5 = mean(dat$var5),
              var6 = mean(dat$var6),
                  var7 = mean(dat$var7),
              var8 = mean(dat$var8),
                  var9 = mean(dat$var9)
                  )

# Get the predicted values
pY <- predict(mod, pX)

# Create a plotting space
plot(dat$var1, dat$y, frame = F, type = "n", xlab = "var1", ylab = "y")

# Plot the points.  Points for var1 below the median on var2 are plotted in red, 
# points for var1 above the median on var2 are plotted in blue
points(dat[dat$var2 < median(dat$var2), "var1"], dat[dat$var2 < median(dat$var2), "y"], pch = 16, cex = .5, col = "red")
points(dat[dat$var2 >= median(dat$var2), "var1"], dat[dat$var2 >= median(dat$var2), "y"], pch = 16, cex = .5, col = "blue")

# Plot the lines. Lines are colored to be consistent with the points
lines(pX[pX$var2 == mean(dat$var2) - sd(dat$var2), "var1"], pY[pX$var2 == mean(dat$var2) - sd(dat$var2)], col = "red", lwd = 2)
lines(pX[pX$var2 == mean(dat$var2) + sd(dat$var2), "var1"], pY[pX$var2 == mean(dat$var2) + sd(dat$var2)], col = "blue", lwd = 2)

Answer 2

Might be worth saying explicitly what's wrong with your proposal & why you should follow the advice given in @Patrick's answer:

First, if the model you're using involves other predictors besides the two involved in the interaction, you clearly need to specify values for all of them to make a prediction using the model.

Second, even if you're only interested in showing the form of the expected response $\operatorname{E} Y$ against two predictors, $x_1$ & $x_2$, consider what happens when the full model is

$$\operatorname{E} Y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \beta_3 x_3 + \ldots + \beta_8 x_8 + \beta_{12} x_1x_2$$

& you fit a reduced model

$$\operatorname{E} Y = \beta_0^* + \beta_1^* x_1 + \beta_2^* x_2 + \beta_{12}^*x_1x_2$$

Does $\beta_1=\beta_1^*$, $\beta_2=\beta_2^*$, & $\beta_{12}=\beta_{12}^*$?

Answer:–

Not in general—only if you have taken pains in the experiment design to ensure orthogonality. So the interaction plot could look completely different for the two models.