Interpretation of interaction effects in linear mixed models with numeric and factorial IVs I hope my question is not too basic for this community. I can*t figure out how to interpret the output of my linear mixed model,  especially the interaction effects.
I do my analysis in R using lme from the nlme package. My model consists of 3 factorial IVs (exp.type [2 levels - field vs. lab], org.type [3 levles - autotroph, heterotroph, mixed], system [2 levels - marine vs. freshwater]) and 2 numeric IVs (even, duration). 

lme1 = lme(reist ~ even + exp.type + even:exp.type + org.type + even:org.type + system + 
        duration, random =~1|authors.year, data = data)


My question is: How do I interpret the interaction effects? How do I get the slopes of the linear regressions for each level of the IVs that have significant interactions
 A: When there is an interaction, there is no unique slope for any of the main effects; those change based on other variables in the interaction. 
I have found that the simplest method of interpreting interactions is visually: make graphs of your DV at different levels of your IVs. Since you have two numeric IVs and 3 categorical ones, you could make 6 graphs. For each, plot one numeric IV on the x axis and the predicted value of the DV on the Y axis, with a line for each level of one of the categorical DVs (you will have to make some assumptions about the value of the other IVs - the mean for the continuous and the mode for the categorical may be sensible). 
Alternatively, you could make a lattice plot (using the lattice package) or use faceting (with the ggplot2 package). 
Another way to go is to make a table of the predicted values of the DV for various common combinations of the IVs. (E.g. the quartiles of the two continuous IVs and all the values of the DVs - which would give a table with 3*3*3*2*2 = 108 rows.
As to what the interactions mean - an interaction means that the relationship between one IV and the DV is different at different levels of the other IV. 
A: Short update: I was wrong in my interpretation, allthough it fits nicely. I did the same analysis using a simple artificial data set w/o mixed effects (allowing me to fit individual lm for each level of both factors:
a = c(rep(1:10, 4))
b = c(10,20,30,40,50,60,70,80,90,100, 5,8,10,14,17,22,27,35,42,50, 90,82,73,64,56,48,40,33,25,18, 5,6,8,10,12,14,17,20,23,26)
c = c(rep("male", 10), rep("female", 10), rep("male", 10), rep("female", 10))
d = c(rep("low", 20), rep("high", 20))
e = data.frame(yval = b, xval = a, sex = c, education = d)

scatterplot(yval~xval | education, smooth = F, grid = T, spread = F, reg.line = T, data = e, xlab = "x", ylab = "y") 

lm2 = lm(yval~xval+sex+xval:sex+education+xval:education, data = e)
summary(lm2)

It turns out that the interaction for the levels of factor A are calculated while ignoring the interaction between the levels of factor B (seeing that I remember that I learned something like that long ago for simple ANOVAs... :).
Strangely, the way I calculated the slopes above in my first post still works, giving the exact numbers from the table output, but completely wrong slopes :)
Calculating:
new_data_lm = expand.grid(xval = c(1:10), sex = c("female", "male"), education = c("high", "low"))
new_data_lm$pred = predict(lm2, new_data_lm, level = 0)
View(new_data_lm)

gives a table with predicted y values for given x values.
When calculating the slopes of these predicted values for each level ob the interactions:
lm_f_h = lm(pred~xval, data = subset(new_data_lm, sex == "female"& education == "high"))
summary(lm_f_h)
lm_f_l = lm(pred~xval, data = subset(new_data_lm, sex == "female"& education == "low"))
summary(lm_f_l)    

gives us the slope -1.503 (for female_high) and 8.794 (for female_low). Both slopes are completely wrong but when added together, they give the exact number for the interaction "xval:educationlow" (10.297) from the output of:
summary(lm2)

The same holds true when calculating the difference in slopes between male_high and female_high giving us -2.667, which is exactly what summary(lm2) gives us for the interaction xval:sex_male.
So all in all I managed to find out how to read the table but I don't understand what I did wrong and why it still fits. Have to think about that very thoroughly.
When predicting the y values from the model by predict(lme2) (so without the additional "new_data" argument), the values of the slopes are correct, but when adding a "new_data" argument with x-values to predict from, the slopes dont match, the interaction terms howver, do?!
So whats the difference between these two approaches? What did I miss?
