Testing one SNP at a time does not fully realise the potential of genome-wide association studies to identify multiple causal variants of small effect, which is a plausible scenario for many complex diseases. Moreover, many simulation studies assume a single causal variant and so more complex realities are ignored. Analysing large numbers of variants simultaneously is now becoming feasible, thanks to developments in Bayesian stochastic search methods. We combine Bayesian shrinkage methods together with a local stochastic model search to identify complex interactions, both local and distal. Our approach can analyse up to 10,000 SNPs simultaneously, and finds multiple potential disease models each with an associated probability. We illustrate its power in comparison with a range of alternative methods, in simulations that incorporate multiple causal loci, acting singly or in interacting pairs, among 4,000 SNPs in a 20Mb region. We argue that, implemented in a two-stage procedure, our hybrid Bayesian analysis can provide a powerful solution to the problem of extracting maximal information from genome-wide association studies.