Visualizing multivariate multiple regression of continuous data in R I have created a multivariate multiple regression model with 3 dependent and 3 independent variables in R, and would like to generate meaningful visualizations. All variables are continuous. When working with multiple regression models with 1 dependent variable, this is fairly easy.
For example
set.seed(0)

df <- data.frame(ind1 = c(1:10),
                ind2 = runif(10,5,15),
                ind3 = runif(10,5,15))
df$dep1 <- df$ind1 * df$ind2 * df$ind3
df$dep2 <- df$ind1 * df$ind2 * df$ind3 * runif(10)
df$dep3 <- df$ind1 * df$ind2 * df$ind3 * abs(rnorm(10))

model1 <- lm(data = df, dep1 ~ ind1 + ind2 + ind3)

There are simple options for different insightful visualizations
library('ggplot2')
library('ggeffects')

ggplot(ggpredict(model1, terms = c("ind1 [1,5,10]", "ind2", "ind3")), 
       aes(x, predicted, color = group)) + geom_line() + facet_wrap(~facet)


library('car')
avPlots(model1)


library('sjPlot')
plot_model(model1, type = 'diag')





However, in a model with 3 dependent variables
model2 <- lm(data = df, cbind(dep1, dep2, dep3) ~ ind1 + ind2 + ind3)

these options seem to go out the window. I am hoping to come up with something a little more powerful than simply
hist(resid(model2))

Let me know if this topic is a better fit for R StackOverflow.

One option is to create a 3D surface plot showing the sensitivity of your dependent variables to the independent variables in your model.
First, create sequences (length = 10 here) along your independent variables.
df2 <- data.frame(ind1 = seq(from = min(df$ind1), to = max(df$ind1), by = (max(df$ind1)-min(df$ind1))/9),
                  ind2 = seq(from = min(df$ind2), to = max(df$ind2), by = (max(df$ind2)-min(df$ind2))/9),
                  ind3 = seq(from = min(df$ind3), to = max(df$ind3), by = (max(df$ind3)-min(df$ind3))/9))

Next, generate a series of predictions.
tst <- predict(model2, df2)

Now, you can plot these predicted values of dep1,dep2,anddep3 as a 3D surface.
library('plotly')

# generate a numeric matrix for use with plotly::plot_ly
mat <- matrix(ncol = ncol(tst), nrow = nrow(tst))

# populate matrix with numeric values from model prediction
for (i in 1:3) {
  mat[,i] <- tst[,i]
}

# generate 3D surface
plot_ly(type = 'surface', z = mat) %>%
  # customize axis titles
  layout(scene = list(xaxis = list(title = 'dep1'),
                      yaxis = list(title = 'dep2'),
                      zaxis = list(title = 'dep3')))

This creates an interactive plot. A still image is shown here

