I'm trying to plot the predictions of a linear model with multiple numeric predictors as a collection of plots such that $\hat Y$ is on the y-axis, each $X_j$ in turn is on the x-axis, and the remaining $X_{i\ne j}$ are all set to plausible values given $X_j$. Here is some example R code (syntax might look strange, it's that way for conciseness, but it does parse correctly in R):
mdata<-data.frame(x1=x1<-rnorm(100,1.5),
x2=x2<-x1+2*rnorm(100,0.5),
x3=x3<-x1*2+rnorm(100,sd=.2),
y=x1-x2+x3-0.5*x1*x2+0.5*x2*x3-0.5*x1*x3+rnorm(100,sd=.1));
mfit<-lm(y~x1*x2*x3,mdata);
Note that in the above dataset, x2 and x3 are both correlated with X1, and all three are predictors for y.
The widespread approach of setting each $X_{i\ne j}=0$ or to the population mean or median is misleading if there is the possibility of multicollinearity among the $X$s because no distinction is made between commonly-occurring combinations of variables and ones that are completely out of the support range of the data. Example using above objects:
predictAtX1<-seq(min(mdata$x1),max(mdata$x1),len=100);
## Do you trust this? I don't.
plot(predict(mfit,newdata=data.frame(x1=predictAtX1,
x2=mean(mdata$x2),
x3=mean(mdata$x3)))~predictAtX1,
xlab='x1',ylab='y-hat');
## Ditto.
plot(predict(mfit,newdata=data.frame(x1=predictAtX1,
x2=median(mdata$x2),
x3=median(mdata$x3)))~predictAtX1,
xlab='x1',ylab='y-hat');
## Ditto, even worse, and yet the most common value at which to fix x_js.
plot(predict(mfit,newdata=data.frame(x1=predictAtX1,
x2=0,
x3=0))~predictAtX1,
xlab='x1',ylab='y-hat');
I thought of fitting each of the $X_j$ to a given $X_i$ and using the predicted values for a given sequence of $X_i$....
plot(predict(mfit,newdata=data.frame(x1=predictAtX1,
x2=predict(lm(x2~x1,mdata),
newdata=list(x1=predictAtX1)),
x3=predict(lm(x3~x1,mdata),
newdata=list(x1=predictAtX1))))~predictAtX1,
xlab='x1',ylab='y-hat')
...and that doesn't look as obviously nonsensical as the plots for globally fixed values of $X_2$ and $X_3$ but given that there are additional interactions involving $X_2$ and $X_3$ that are being ignored, these predictor values are not representative of the population either.
But. For a given sequence of values for $X_1$ it should be possible to construct a smoothed path through the variable space corresponding to the most likely values of $X_2$ and $X_3$. I can think of brute-force optimization approaches, but given how obvious and universal this problem is for any real-life regression model with multiple predictors, there must be an already-implemented and/or simpler approach to dealing with it.
So, my revised question is, does anybody know the name of such a method or the name of the subfield of statistics that deals with it, or the name of an R package that implements it?
Thanks.