Surface Plot for Logistic Regression Interactions [closed]

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!

closed as off-topic by Sycorax, Michael Chernick, Jeremy Miles, Peter Flom♦Jun 9 '18 at 14:22

This question appears to be off-topic. The users who voted to close gave this specific reason:

• "This question appears to be off-topic because EITHER it is not about statistics, machine learning, data analysis, data mining, or data visualization, OR it focuses on programming, debugging, or performing routine operations within a statistical computing platform. If the latter, you could try the support links we maintain." – Michael Chernick, Jeremy Miles, Peter Flom
If this question can be reworded to fit the rules in the help center, please edit the question.

• I dislike 3-D graphs and there is evidence (see the work of William S. Cleveland) that they are very hard to interpret correctly and easy to get completely wrong. – Peter Flom Jun 9 '18 at 14:22

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. 