How to plot regression or LOWESS lines over data i

It is often helpful to plot a line (e.g., a fitted regression line or a non-parametric LOWESS line) over data. Likewise, when variables are confounded, it is often helpful to plot x1, y data stratified by x2. However, I can't quite figure out how to combine both using ?coplot. This ought to be quite natural in lattice, but I can't seem to figure out how to do it there, either.

The following provides some data to play with:

library(MASS)
set.seed(7422)
X <- mvrnorm(100, mu=c(5,5), Sigma=rbind(c(2.0, 1.9),
                                         c(1.9, 2.0) ) )
x1 <- X[,1];  x2 <- X[,2]
y <- 3 + .6*x1 - .4*x2 + rnorm(100)
m <- lm(y~x1+x2)
summary(m)$coefficients
#               Estimate Std. Error   t value     Pr(>|t|)
# (Intercept)  3.8281801  0.4022204  9.517619 1.476519e-15
# x1           0.5466495  0.2060269  2.653292 9.314673e-03
# x2          -0.4835804  0.2139542 -2.260205 2.604451e-02
windows()
  pairs(cbind(y, x1, x2))
xs <- seq(1.6, 9.2, by=.1)
windows()
  plot(y~x2)
  lines(xs, predict(m, data.frame(x1=5.2, x2=xs)), col="blue")
windows()
  coplot(y~x2|x1)
  ## doesn't work (lines on plot were drawn by hand)
  # lines(xs, predict(m, data.frame(x1=mean(x1[subscripts], x2=xs)), col="blue")
  # lines(lowess(y~x2, subset=x1[subscripts]), col="gray")

Scatterplot matrix of the raw data:

Here is a marginal plot of x2 vs y, with fitted the regression line and a marginal LOWESS line plotted on top of the data. This will look weird to people who aren't statistically savvy.

This is what I have in mind. It's a coplot, with the regression model and a LOWESS line fitted to each stratum of the data plotted inside each panel.

To do this with coplot, you need to define your own panel function. The default function is points which just plots the points, but you can add a regression line and a lowess curve with a simple function.

PointsWithReg = function(x, y, pch=20, col=1, ...) {
    points(x=x, y=y, col=col, pch=pch, type="p", ...)
    abline(lm(y ~ x))
    lines(lowess(x=x, y=y), col="blue")
}

coplot(y~x2|x1, panel=PointsWithReg, pch=16)