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.

Hypothesis tests in multivariate multiple regression models can be visualized using hypothesis-error (HE) plots, implemented in the heplots package, https:friendly.github.io/heplots/
These show data ellipses for the residuals (E) in relation to the predictors (H).
For 1 df regression terms, the H "ellipses" collapse to a line.  They have the property that the H ellipse projects outside the E ellipse if and only if the term is significant by Roy's test.
For this example, we get:
library(heplots)
heplot(model2, fill=TRUE, fill.alpha=0.1)


Here, only ind1 is significant.
> car::Anova(model2, test="Roy")

Type II MANOVA Tests: Roy test statistic
     Df test stat approx F num Df den Df Pr(>F)  
ind1  1      6.14     8.19      3      4  0.035 *
ind2  1      1.43     1.91      3      4  0.270  
ind3  1      1.45     1.94      3      4  0.265  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

A scatterplot matrix version is obtained using pairs.mlm
pairs(model2, , fill=TRUE, fill.alpha=0.1)

The angles of the H lines for the predictors ind1, ind2, and ind3 wrt the dep variables show the correlations of the former with the latter.

For more complex models, involving factors (df > 1 terms), these results can be visualized in "canonical space" using the candisc package, candisc
