Other than literally testing each possible combination of variable(s) in a model (x1:x2 or x1*x2 ... xn-1 * xn). How do you identify if an interaction SHOULD or COULD exist between your independent (hopefully) variables?
What are best practices in attempting to identify interactions? Is there a graphical technique that you could or do use?
Cox and Wermuth (1996) or Cox (1984) discussed some methods for detecting interactions. The problem is usually how general the interaction terms should be. Basically, we (a) fit (and test) all second-order interaction terms, one at a time, and (b) plot their corresponding p-values (i.e., the No. terms as a function of $1-p$). The idea is then to look if a certain number of interaction terms should be retained: Under the assumption that all interaction terms are null the distribution of the p-values should be uniform (or equivalently, the points on the scatterplot should be roughly distributed along a line passing through the origin).
Now, as @Gavin said, fitting many (if not all) interactions might lead to overfitting, but it is also useless in a certain sense (some high-order interaction terms often have no sense at all). However, this has to do with interpretation, not detection of interactions, and a good review was already provided by Cox in Interpretation of interaction: A review (The Annals of Applied Statistics 2007, 1(2), 371–385)--it includes references cited above. Other lines of research worth to look at are study of epistatic effects in genetic studies, in particular methods based on graphical models (e.g., An efficient method for identifying statistical interactors in gene association networks).
My best practice would be to think about the problem to hand before fitting the model. What is a plausible model given the phenomenon you are studying? Fitting all possible combinations of variables and interactions sounds like data dredging to me.
Fitting a tree model (i.e. using R), will help you identify complex interactions between the explanatory variables. Read the example on page 30 here.
I'll preface this response as I entirely agree with Gavin, and if you're interested in fitting any type of model it should be reflective of the phenomenon under study. What the problem is with the logic of identifying any and all effects (and what Gavin refers to when he says data dredging) is that you could fit an infinite number of interactions, or quadratic terms for variables, or transformations to your data, and you would inevitably find "significant" effects for some variation of your data.
As chl states, these higher order interaction effects don't really have any interpretation, and frequently even the lower order interactions don't make any sense. If your interested in developing a causal model you should only include terms you believe could be pertinent to your dependent variable A priori to fitting your model.
If you believe they can increase predictive power of your model, you should look up resources on model selection techniques to prevent over-fitting your model.
How large is $n$ ? how many observations do you have ? this is crucial ...
Sobol indices will tell you the proportion of variance explained by interaction if you have a lot of observations and a few $n$, otherwise you will have to do modelling (linear to start with). You have a nice R package for that called sensitivity. Anyway the idea is quite often that of decomposing the variance (also called generalized ANOVA).
If you want to know if this proportion of variance is significant you will have to do modelling (roughly, you need to know the number of degrees of freedom of your model to compare it to the variance).
Are your variables discrete or continuous ? bounded or not really (i.e you don't know the maximum) ?
I think this is a good use case of LASSO. You can throw all the interaction terms, and LASSO can find the ones that matter, by using cross-validation to select the best regularization parameter. This doesn't need to be confined to linear interaction, as we can add much richer class of interactions (e.g. $x_m^2 * x_n$ terms)
For LASSO, you can have p >> n, so you don't have to worry about identification by having more covariates than number of observations. For a smaller dimension problem(n >> p), LASSO is pretty much the same as running t-test on each interaction terms, so I think it'd work similar to chl's answer.
