How to investigate a 3-way interaction? There is a significant 3-way interaction in a data-set I'm working with.
The interaction involves both categorical and quantitative variables.
I have been pointed towards simple slopes and this website but I find the explanations lacking. I only have a basic background in statistics, and Googling for other examples has not been very helpful. 
Any insight as to how to begin to understand this interaction will be very welcome, especially using R rather than SPSS. Thank you very much in advance, from a lower-year college student who has suddenly found himself over his head!
 A: +1 to comments above about graphing interaction.
An initial approach to thinking about a three-way interaction is that it is saying that the pattern of results contained in the interaction between A and B depends upon the level/value of C. The following is framed in a linear regression kind of framework, but conceptually similar to e.g. logistic regression.
This assumes that you're happy with thinking about two-way interactions, of course :-)
Let's say there are three predictors, A (continuous), B (categorical), and C (categorical).
A two way A*B interaction would indicate that the slope of A (relating to the outcome) depends on the level of B to which an individual belongs.
In this context, a three way A*B*C interaction would indicate that the A*B interaction (previously discussed differential slopes of A on the outcome according to group of B) depends on the group of C to which one belongs... 
Things are much more complex if there are two continuous/ordinal predictors in the interaction (usually too complex to conceptualise: I would be interested to hear others' opinions on this.)
A: I was afraid this was going to be a continuous $\times$ continuous $\times$ categorical interaction... OK, here goes - first, we create some toy data (foo is a binary predictor, bar and baz are continuous, dv is the dependent variable):
set.seed(1)
obs <- data.frame(foo=sample(c("A","B"),size=100,replace=TRUE),
  bar=sample(1:10,size=100,replace=TRUE),
  baz=sample(1:10,size=100,replace=TRUE),
  dv=rnorm(100))

We then fit the model and look at the three-way interaction:
model <- lm(dv~foo*bar*baz,data=obs)
anova(update(model,.~.-foo:bar:baz),model)

Now, for understanding the interaction, we plot the fits. The problem is that we have three independent variables, so we would really need a 4d plot, which is rather hard to do ;-). In our case, we can simply plot the fits against bar and baz in two separate plots, one for each level of foo. First calculate the fits:
fit.A <- data.frame(foo="A",bar=rep(1:10,10),baz=rep(1:10,each=10))
fit.A$pred <- predict(model,newdata=fit.A)
fit.B <- data.frame(foo="B",bar=rep(1:10,10),baz=rep(1:10,each=10))
fit.B$pred <- predict(model,newdata=fit.B)

Next, plot the two 3d plots side by side, taking care to use the same scaling for the $z$ axis to be able to compare the plots:
par(mfrow=c(1,2),mai=c(0,0.1,0.2,0)+.02)
persp(x=1:10,y=1:10,z=matrix(fit.A$pred,nrow=10,ncol=10,byrow=TRUE),
  xlab="bar",ylab="baz",zlab="fit",main="foo = A",zlim=c(-.8,1.1))
persp(x=1:10,y=1:10,z=matrix(fit.B$pred,nrow=10,ncol=10,byrow=TRUE),
  xlab="bar",ylab="baz",zlab="fit",main="foo = B",zlim=c(-.8,1.1))

Result:

We see how the way the fit depends on (both!) bar and baz depends on the value of foo, and we can start to describe and interpret the fitted relationship. Yes, this is hard to digest. Three-way interactions always are... Statistics are easy, interpretation is hard...
Look at ?persp to see how you can prettify the graph. Browsing the R Graph Gallery may also be inspirational.
A: Often it helps to plot the relationships and see how things change when you change one of the variables.  A couple of tools in R that help with these plots are Predict.Plot and TkPredict functions in the TeachingDemos package.
