Adding a regression plane to plot_ly 3d plot

Question

I have a multiple regression that plots well in plot_ly, however I cannot seem to get a regression line from ggpredict of a GLMM to overlay. The model I want plotted is structured as the following:

my.lm <- glmmTMB(z ~ a * b + x + y + (1|g), data = MSdf, family = beta_family()), 

where a and b are categorical binomials (0 or 1) and x and y are continuous numeric variables (density gradients).

So far I have the following which makes a nice scatterplot:

p = plot_ly(x = ~x,
            y = ~y,
            z = ~z,
            type = "scatter3d",
            size = 1)

I want to add the predicted fit:

predic = ggpredict(my.lm, terms = c('x', 'y'), type = 'fixed') 

to the above scatterplot. Add_surface seems to be the best bet, but my z value is a vector rather than a numeric matrix, which is odd becuase in the examples I see, their z value is a vector. In other words, I get this error: Error: z must be a numeric matrix. I've tried a bunch of things to structure this matrix (i.e. z = as.numeric(Emig) %>% z = as.data.frame(z); cbind(x,y,z) to name a couple) but nothing seems to work.

I've done the same with scatterplot3d and it yields another issue the same issue.

Answer 1

You are right, add_surface is the way to do this. Here is how you can achieve this. Note I just added some sample data:

library(glmmTMB)
library(ggeffects)
library(plotly)
library(dplyr)
library(tidyr)

set.seed(123)
n <- 100
MSdf <- data.frame(
  z = rbeta(n, 2, 5),
  a = sample(0:1, n, replace = TRUE),
  b = sample(0:1, n, replace = TRUE),
  x = runif(n, 0, 10),
  y = runif(n, 0, 10),
  g = factor(sample(1:10, n, replace = TRUE))
)

my.lm <- glmmTMB(z ~ a * b + x + y + (1|g), data = MSdf, family = beta_family())

x_seq <- seq(min(MSdf$x), max(MSdf$x), length.out = 100)
y_seq <- seq(min(MSdf$y), max(MSdf$y), length.out = 100)
common_g <- names(sort(table(MSdf$g), decreasing = TRUE))[1]  
grid <- expand.grid(x = x_seq, y = y_seq, a = unique(MSdf$a), b = unique(MSdf$b), g = common_g)

grid$z <- predict(my.lm, newdata = grid, type = "response")
grid_avg <- grid %>%
  group_by(x, y) %>%
  summarize(z = mean(z), .groups = 'drop')

z_matrix <- matrix(grid_avg$z, nrow = length(x_seq), ncol = length(y_seq))

p <- plot_ly(MSdf, x = ~x, y = ~y, z = ~z, type = "scatter3d", mode = "markers", marker = list(size = 2))
p <- p %>% add_surface(x = x_seq, y = y_seq, z = z_matrix)
p

Which gives