With levelplot I would like to combine fitted values and "raw" values in one plot. The surface should be based on fitted values and on top of those I would like distinguisable points based on the raw data with a colorscheme the same as the surface. # raw data x1 = rep(1:10, times = 10) x2 = rep(1:10, each = 10) y = 2 + (x1+rnorm(100, mean=0, sd=4)) + (x2+rnorm(100, mean=0, sd=2)) levelplot(y~x1+x2) # predicted data fit = lm(y ~ x1 + x2) x1.mesh = do.breaks(range(x1), 50) x2.mesh = do.breaks(range(x2), 50) grid = expand.grid(x1 = x1.mesh, x2 = x2.mesh) grid[["fit"]] = predict(fit, newdata=grid) levelplot(fit ~ x1 + x2, data=grid) So, instead of plotting both levelplots (like here) I want the first on top of the second. I tried a bunch of things (including panel.xyplot within panel) but am now doubting whether it is possible? Thanks in advance, Eelke