Multivariate/interaction prediction from glmm Previously I have asked how to calculate the predicted response for groups (split by two categorical variables) given a single continuous fixed effect, in a glmm.
Now I would like to take it one step further; to estimate a predicted response when there are two continuous variables ($X_1$ and $X_2$). Continuing from my previous question, assume the model:
mod = MCMCglmm(Life ~ X1*X2*Sex*Group,
    random = ~ Pop:Year,
    rcov = ~units,  
    family = "gaussian", 
    start = list(QUASI = FALSE), 
    data = DF1) 

where sex has two levels ($M$ and $F$) and group has two levels ($LO$ and $LL$). How do I then predict the lifespan (Y) for each group under two continuous variables? For example, how could I predict the lifespan for males from group LO, when I have known values of $X_1$ and $X_2$.

Note: The plan is to a) examine the difference between the groups at given values of $X_1$ and $X_2$ and b) plot surfaces (contour plots) for each group across a spectrum of $X_1$ and $X_2$ values - tips on implementation in R would be useful too.
 A: Where the output of summary(mod) looks something like:
                            mean   L.CI...
b0    (Intercept)           20.3 -125.5...
b1    X1                     0.1    0.9...
b2    X2                    22.3    2.3...
b3    SexM                   6.7  -32.4...
b4    GroupLO               -1.4 -127.0...
b5    X1:X2                  9.3  -18.4...
b6    X1:SexM                3.7  -16.3...
b7    X2:SexM               -8.8  -49.7...
b8    X1:GroupLO             5.6  -13.1...
b9    X2:GroupLO            21.3  -23.8...
b10   SexM:GroupLO         -10.6 -129.4...
b11   X1:X2:SexM             5.6   -2.2...
b12   X1:X2:GroupLO         -6.5  -12.8...
b13   X1:SexM:GroupLO        7.5  -28.5...
b14   X2:SexM:GroupLO        0.4  -76.5...
b15   X1:X2:SexM:GroupLO     0.1   -8.7...

such that, for example, $b_5$ is the interaction between $X_1$ and $X_2$, the predicted response for each combination of group and sex can be attained from:
$ Life_{F,LL} = b_0 + (b_1 \times X_1) + (b_2 \times X_2) + (b_5 \times X_1 \times X_2) $
$ Life_{M,LL} = (b_0 + b_3) + ((b_1 + b_6) \times X_1) + ((b_2 + b_7) \times X_2) + ((b_5 + b_{11}) \times X_1 \times X_2) $
$ Life_{F,LO} = (b_0 + b_4) + ((b_1 + b_8) \times X_1) + ((b_2 + b_9) \times X_2) + ((b_5 + b_{12}) \times X_1 \times X_2) $
$ Life_{M,LO} = (b_0 + b_3 + b_4 + b_{10}) + ((b_1 + b_6 + b_7 + b_{13}) \times X_1) + ((b_2 + b_7 + b_9 + b_{14}) \times X_2) + ((b_5 + b_{11} + b_{12} +b_{15}) \times X_1 \times X_2)$ 
The first term in the calculations is the intercept and is composed of all relevant estimates not featuring $X_1$ or $X_2$. The remaining parts of the formula are the slope effects, with the $X_1$ specific, $X_2$ specific, and then interaction between $X_1$ and $X_2$ specific terms.
These formula have been checked using the in built predict function in R.
