Load packages library(tidyverse) library(plotly) Data Some sample data data(tips, package= "reshape2") Regression model lm1 <- lm(tip ~ total_bill + size, data = tips) lm1_coef <- coef(lm1) Sequence x1_seq <- seq(min(tips$total_bill), max(tips$total_bill), length.out = 25) x2_seq <- seq(min(tips$size), max(tips$size), length.out = 6) Compute grid z2 <- t(outer(x1_seq, x2_seq, function(x,y) lm1_coef[1]+lm1_coef[2]*x+lm1_coef[3]*y)) z2 #> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] #> [1,] …