Surface Plot for Logistic Regression Interactions I'm trying to visualize some different interactions from a logistic regression in R. I'd like create a surface plot of the predictive model with two predictor variables along the x and y, then the binary prediction on the z.
I've tried using plotly, geoR, persp, bplot, and a few other methods without much success. Almost had success with lattice, but I'm not sure how to plot the predictive model instead of the raw data.
Some background: below are simulated data modeled after an experiment where an observer makes a yes/no decision about a stimulus presented at 9 different settings ranging from 0.6 to 1.4 in steps of 0.1. The setting (ordinal) and an additional parameter ranging from 0.02626 to 0.14749 (continuous) were used to predict the yes/no using glmer a la:
outcome ~ (random effect) + setting * parameter
The results showed a significant negative interaction between the setting and the parameter; this is what I'm trying to visualize.
outcome <- sample(x = c(0,1), size = 4186, replace = TRUE)
setting <- sample(x = c(0.6,0.7,0.8,0.9,1.0,1.1,1.2,1.3,1.4), size = 4186, replace = TRUE)
parameter <- runif(4186, min=0.02626, max=0.14749)

The closest I've come is in using this code:
library(lattice)
library(gridExtra)

trellis.par.set("axis.line", list(col=NA,lty=1,lwd=1))
jobsat_plot <- wireframe(outcome ~ setting*parameter, data=data, 
          xlab = "setting",
          ylab = "parameter",
          main = "Logistic - setting x parameter",
          drape = TRUE,
          colorkey = TRUE,
          scales = list(arrows=FALSE,cex=.5, tick.number = 10, z = list(arrows=F), distance =c(1.5, 1.5, 1.5)),
          light.source = c(10,0,10),
          col.regions = rainbow(100, s = 1, v = 1, start = 0, end = max(1,100 - 1)/100, alpha = .8),
          screen = list(z = -60, x = -60)
)

grid.arrange(jobsat_plot, ncol=1, clip=TRUE)

Which produces something like this (with the real data; can't seem to get it to work with the simulated data):

I think this gets the point across, but I'd like the smooth predictive model to show the general trends. Any help, or pointing in the right direction is appreciated!
 A: The reason your code doesn't work is that it's plotting the raw data, not predictions made from a model.
To get what you want, you just need to make a grid of points at which to make model predictions, the predictions at these grid points, and then supply that as the data to wireframe.
library(lattice)
set.seed(42)

dat = data.frame(
    outcome=sample(x = c(0,1), size = 4186, replace = TRUE),
    setting=sample(x = c(0.6,0.7,0.8,0.9,1.0,1.1,1.2,1.3,1.4), size = 4186, replace = TRUE),
    parameter=runif(4186, min=0.02626, max=0.14749)
    )
model <- glm(outcome ~ setting*parameter, data=dat, family="binomial")

x_tilde <- expand.grid(setting=seq(0.6, 1.4, length.out=10), parameter=seq(0.02626, 0.14749, length.out=10))
x_tilde$prediction <- predict(model, x_tilde, type="response")

wireframe(prediction ~ setting + parameter, data=x_tilde,
    xlab = "setting",
          ylab = "parameter",
          main = "Logistic - setting x parameter",
          drape = TRUE,
          colorkey = TRUE,
          scales = list(arrows=FALSE,cex=.5, tick.number = 10, z = list(arrows=F), distance =c(1.5, 1.5, 1.5)),
          light.source = c(10,0,10),
          col.regions = rainbow(100, s = 1, v = 1, start = 0, end = max(1,100 - 1)/100, alpha = .8),
          screen = list(z = -60, x = -60))

Obviously you'll need to fine-tune how to display the plot, but this is the general idea.

