[R-sig-ME] Partial regression plot for mixed models?

[Stefano Leonardi]

Hi,
I am fitting a fairly complicated mixed gaussian models with lme.
The dependent variable is a measure taken many times over
the same individuals (repeated measure).

The model involves several fixed explanatory variables, some 
interactions between them and also a random intercept and two random 
slopes. Residuals are correlated in time so I also included a
correlation matrix estimate with corARMA.

The call to lme is like this:

lme(Y ~ X1 + X2 + X3 + X4 + X5 +
     I(X1^2) + I(X5^2) + X1:X3 + X2:X4 + X3:X4,
     data = d.df,
     random = list(ID = pdDiag(~X1 + X5)),
     correlation = corARMA(p = 1, q = 1))

I obtain this model after a quite difficult simplification.

I would like to graphically show the effect of X1 and X5 (and possible 
also the other X) on Y.

I was thinking to use partial regression plots.

Is there an easy way to do this?
How should I compute regression on Y on all X except Xi and get residuals?
And what about residuals on regression on Xi on all X except Xi?

I don't know how to consider the Xi when they are included in the
random terms.

Any suggestion?
I would be very grateful. :-)

Stefano