Moderation in R: Plot moderation for regression model with one interaction and controll variable I'm testing a moderation using R and would like to plot it.
This is my regression model (2 predictors + interaction):
PBC.lm <- lm(IntAll ~ cPB01*cSK, data=data)

For plotting, normally I use this function of QuantPsyc library:
par(mfrow = c(1,1), mar = c(4, 4, 2, 1))
simpleSlope <- sim.slopes(mod=PBC.lm, z=data$cSK, zsd=1, mcz=T)
graph.mod(simpleSlope, x=cPB01, y=IntAll, data=data, 
          ylab="Intention" , xlab="centered PBC",
          title = "Relation of PBC and Intention")

However, I want to do the same, cotrolling for one further predictor. Thus, my regression model is:
PBC.lm2 <- lm(IntAll ~ cPB01*cSK + gcKnowledgeScore, data=data)

My question now is, can I use the same command (below) when including a third predictor? If not, what would be the correct alternative to display the moderation?
par(mfrow = c(1,1), mar = c(4, 4, 2, 1))
simpleSlope <- sim.slopes(mod=PBC.lm2, z=data$cSK, zsd=1, mcz=T)
graph.mod(simpleSlope, x=cPB01, y=IntAll, data=data, 
          ylab="Intention" , xlab="centered PBC",
          title = "Relation of PBC and Intention")

Thank you!
 A: alternative manual approach
It doesn't seem like sim.slopes is able to generate simple slopes for 3-way or higher order interactions. You could do this manually.

# create data
x  <- rep(1:10,4)
z1 <- c(rep(0,20),rep(1,20))
z2 <- c(rep(0,10),rep(1,10),rep(0,10),rep(1,10))
y <- 1 + x + z1 + z1*x + 3*z2 + rnorm(40,0,1)
data <- as.data.frame(cbind(y, x, z1, z2))

# model
m <- lm(y ~ 1 + x*z1 + z2, data = data)

# plot points
colours = hsv((1+z1+2*z2)/6,rep(1,40),rep(0.5,40)+(1+z1+2*z2)/8)
plot(x, y, pch=21, bg=colours, col = colours, ylim = c(0,25))

# custom function to plot simple lines 
# thus uses predict to get to the estimates of the line coefficients
#     the se is here not just the se of the trend lines 
#     but also the se of the lines *plus* error of the estimate in the mean y
sslineplot <- function(ix,iz1,iz2) {
  n <- length(ix)
  colours1 = rep(hsv((1+iz1+2*iz2)/6,1,0.5+(1+iz1+2*iz2)/8),n)
  colours2 = rep(hsv((1+iz1+2*iz2)/6,1,0.5+(1+iz1+2*iz2)/8,0.25),n)
  yp <- predict(m, newdata = as.data.frame(cbind(x = ix, z1=rep(iz1,n), z2=rep(iz2,n), y=rep(0,n))),se.fit = TRUE)
  lines(ix, yp$fit, col = colours1)
  lines(ix, yp$fit + 2*yp$se.fit, col = colours1, lty=2)
  lines(ix, yp$fit - 2*yp$se.fit, col = colours1, lty=2)
  polygon(c(rev(ix), ix), c(rev(yp$fit + 2*yp$se.fit), yp$fit - 2*yp$se.fit), col = colours2, border = NA)

}

sslineplot(seq(0,12,0.1),0,0)
sslineplot(seq(0,12,0.1),0,1)
sslineplot(seq(0,12,0.1),1,0)
sslineplot(seq(0,12,0.1),1,1)

Note that this function uses the predict function. This predicts the $y$ values instead of the slope parameters, but underlying those predictions for $y$ is that the prediction of slope parameters.
what does sim.slopes do?
If you type the function sim.slopes directly into the console then you get to see the source code and you will be able to evaluate the way in which the slopes are calculated. I am not sure whether the way to calculate those slopes is actually standard. At least I think that it can be done in different ways and the approach here is not necessarily the best.


*

*Sim.slopes calculates slopes/lines for particular values of $z$. This works by re-ordering the equation:
$$ y_i = a + bx_i + c z_i + d(x_i z_i) = (a + c z_i) + (b+d z_i)  x_i$$
the error of the coefficients $(a + c z_i)$ and $(b + d z_i)$ are calculated using the standard error for the coefficients $a$, $b$, $c$, $d$ and their correlations.

*Sim.slopes computes the standard error for those lines only for the slope parameter but not for the intercept parameter. 

*Sim.slopes plots three lines (one centered and two extremes) based on the variation of the $z$ parameter (which may be useful only when $z$ is not categorical) but not on the variation of the error of the lines. You won't be able to do this for two continuous variables. Sim.slopes chooses a high and low value for $z$ giving an idea about the variation of the slope for that parameter. But with varying two parameters together you do not have two outer ends (what you could do is find out the mutual direction of change that changes the slope the most and choose a low and high point on that direction)
It seems to me that you can vary a lot within this theme.
> sim.slopes
function (mod, z, zsd = 1, mcz = FALSE) 
{
    if (!mcz) {
        z <- z - mean(z, na.rm = TRUE)
    }
    else {
        z <- z
    }
    zhi <- mean(z, na.rm = TRUE) + zsd * sd(z, na.rm = TRUE)
    zlo <- mean(z, na.rm = TRUE) - zsd * sd(z, na.rm = TRUE)
    zme <- mean(z, na.rm = TRUE)
    b0 <- summary(mod)$coef[1, 1]
    b1 <- summary(mod)$coef[2, 1]
    b2 <- summary(mod)$coef[3, 1]
    b3 <- summary(mod)$coef[4, 1]
    x.zhi <- (b1 + b3 * zhi)
    x.zlo <- (b1 + b3 * zlo)
    x.zme <- (b1 + b3 * zme)
    int.zhi <- (b0 + b2 * zhi)
    int.zlo <- (b0 + b2 * zlo)
    int.zme <- (b0 + b2 * zme)
    seb.zhi <- sqrt(vcov(mod)[2, 2] + 2 * zhi * vcov(mod)[2, 
        4] + zhi^2 * vcov(mod)[4, 4])
    seb.zlo <- sqrt(vcov(mod)[2, 2] + 2 * zlo * vcov(mod)[2, 
        4] + zlo^2 * vcov(mod)[4, 4])
    seb.zme <- sqrt(vcov(mod)[2, 2] + 2 * zme * vcov(mod)[2, 
        4] + zme^2 * vcov(mod)[4, 4])
    td <- qt(0.975, df = summary(mod)$df[2])
    zhi.u <- x.zhi + td * seb.zhi
    zhi.l <- x.zhi - td * seb.zhi
    zlo.u <- x.zlo + td * seb.zlo
    zlo.l <- x.zlo - td * seb.zlo
    zme.u <- x.zme + td * seb.zme
    zme.l <- x.zme - td * seb.zme
    mat <- matrix(NA, 3, 5)
    colnames(mat) <- c("INT", "Slope", "SE", "LCL", "UCL")
    rownames(mat) <- c("at zHigh", "at zMean", "at zLow")
    mat[1, ] <- c(int.zhi, x.zhi, seb.zhi, zhi.l, zhi.u)
    mat[2, ] <- c(int.zme, x.zme, seb.zme, zme.l, zme.u)
    mat[3, ] <- c(int.zlo, x.zlo, seb.zlo, zlo.l, zlo.u)
    return(data.frame(mat))
}
<environment: namespace:QuantPsyc>

