# 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.