apparently the Bayes factor somehow uses likelihoods that represent the likelihood of each model integrated over it's entire parameter space (i.e. not just at the MLE). How is this integration actually achieved typically? Does one really just try to calculate the likelihood at each of thousands (millions?) of random samples from the parameter space, or are there analytic methods to integrating the likelihood across the parameter space?
First, any situation where you consider a term such as $P(D|M)$ for data $D$ and model $M$ is considered a likelihood model. This is often the bread and butter of any statistical analysis, frequentist or Bayesian, and this is the portion that your analysis is meant to suggest is either a good fit or a bad fit. So Bayes factors are not doing anything fundamentally different than likelihood ratios.
It's important to put Bayes factors in their right setting. When you have two models, say, and you convert from probabilities to odds, then Bayes factors act like an operator on prior beliefs:
$$ Posterior Odds = Bayes Factor * Prior Odds $$
$$ \frac{P(M_{1}|D)}{P(M_{2}|D)} = B.F. \times \frac{P(M_{1})}{P(M_{2})} $$
The real difference is that likelihood ratios are cheaper to compute and generally conceptually easier to specify. The likelihood at the MLE is just a point estimate of the Bayes factor numerator and denominator, respectively. Like most frequentist constructions, it can be viewed as a special case of Bayesian analysis with a contrived prior that's hard to get at. But mostly it arose because it's analytically tractable and easier to compute (in the era before approximate Bayesian computational approaches arose).
To the point on computation, yes: you will evaluate the different likelihood integrals in the Bayesian setting with a large-scale Monte Carlo procedure in almost any case of practical interest. There are some specialized simulators, such as GHK, that work if you assume certain distributions, and if you make these assumptions, sometimes you can find analytically tractable problems for which fully analytic Bayes factors exist.
But no one uses these; there is no reason to. With optimized Metropolis/Gibbs samplers and other MCMC methods, it's totally tractable to approach these problems in a fully data driven way and compute your integrals numerically. In fact, one will often do this hierarchically and further integrate the results over meta-priors that relate to data collection mechanisms, non-ignorable experimental designs, etc.
I recommend the book Bayesian Data Analysis for more on this. Although, the author, Andrew Gelman, seems not to care too much for Bayes factors. As an aside, I agree with Gelman. If you're going to go Bayesian, then exploit the full posterior. Doing model selection with Bayesian methods is like handicapping them, because model selection is a weak and mostly useless form of inference. I'd rather know distributions over model choices if I can... who cares about quantizing it down to "model A is better than model B" sorts of statements when you do not have to?
Additionally, when computing the Bayes factor, does one apply correction for complexity (automatically via cross-validated estimation of likelihood or analytically via AIC) as one does with the likelihood ratio?
This is one of the nice things about Bayesian methods. Bayes factors automatically account for model complexity in a technical sense. You can set up a simple scenario with two models, $M_{1}$ and $M_{2}$ with assumed model complexities $d_{1}$ and $d_{2}$, respectively, with $d_{1} < d_{2}$ and a sample size $N$.
Then if $B_{1,2}$ is the Bayes factor with $M_{1}$ in the numerator, under the assumption that $M_{1}$ is true one can prove that as $N\to\infty$, $B_{1,2}$ approaches $\infty$ at a rate that depends on the difference in model complexity, and that the Bayes factor favors the simpler model. More specifically, you can show that under all of the above assumptions, $$ B_{1,2} = \mathcal{O}(N^{\frac{1}{2}(d_{2}-d_{1})}) $$
I'm familiar with this derivation and the discussion from the book Finite Mixture and Markov Switching Models by Sylvia Frühwirth-Schnatter, but there are likely more directly statistical accounts that dive more into the epistemology underlying it.
I don't know the details well enough to give them here, but I believe there are some fairly deep theoretical connections between this and the derivation of AIC. The Information Theory book by Cover and Thomas hinted at this at least.
Also, what are the philosophical differences between the likelihood ratio and the Bayes factor (n.b. I'm not asking about the philosophical differences between the likelihood ratio and Bayesian methods in general, but the Bayes factor as a representation of the objective evidence specifically). How would one go about characterizing the meaning of the Bayes factor as compared to the likelihood ratio?
The Wikipedia article's section on "Interpretation" does a good job of discussing this (especially the chart showing Jeffreys' strength of evidence scale).
Like usual, there's not too much philosophical stuff beyond the basic differences between Bayesian methods and frequentist methods (which you seem already familiar with).
The main thing is that the likelihood ratio is not coherent in a Dutch book sense. You can concoct scenarios where the model selection inference from likelihood ratios will lead one to accept losing bets. The Bayesian method is coherent, but operates on a prior which could be extremely poor and has to be chosen subjectively. Tradeoffs.. tradeoffs...
FWIW, I think this kind of heavily parameterized model selection is not very good inference. I prefer Bayesian methods and I prefer to organize them more hierarchically, and I want the inference to center on the full posterior distribution if it is at all computationally feasible to do so. I think Bayes factors have some neat mathematical properties, but as a Bayesian myself, I am not impressed by them. They conceal the really useful part of Bayesian analysis, which is that it forces you to deal with your priors out in the open instead of sweeping them under the rug, and allows you to do inference on full posteriors.
You can compute the likelihood of a model indeed using model.logp(). As input, it requires a point. For example, the BEST model from the examples directory I can do:
np.exp(model.logp({'group1_mean': 0.1,
'group2_mean': 0.2,
'group1_std_interval': 1.,
'group2_std_interval': 1.2,
'nu_minus_one_log': 1}))
Note that this model is using transformed variables, so I have to supply these. You could then take the exp() of that and use it inside a numerical integrator, for example as provided by scipy.integrate. The problem is that even with only 5 parameters, this will be very slow.
Bayes Factors are generally very difficult to compute because you have to integrate over the complete parameter space. There are some ideas to using MCMC samples for that. See this post, and especially the comment section for more information: https://radfordneal.wordpress.com/2008/08/17/the-harmonic-mean-of-the-likelihood-worst-monte-carlo-method-ever/ The case for BIC is unfortunately similar.
If you really want to compute the Bayes Factor, you can also look at the Savage Dickey Ratio test (see e.g. http://drsmorey.org/bibtex/upload/Wagenmakers:etal:2010.pdf), but it's application is limited.
I suppose that you're trying to do model comparison which is a field with many opinions and solutions (some hard to implement, like BFs). One measure that is very easy compute is the Deviance Information Criterion. It has its downsides, although some of them can be remedied (see http://onlinelibrary.wiley.com/doi/10.1111/rssb.12062/abstract). Unfortunately we haven't ported the code pymc3 yet, but it'd be pretty easy (see here for the pymc2 implementation: https://github.com/pymc-devs/pymc/blob/895c24f62b9f5d786bce7ac4fe88edb4ad220364/pymc/MCMC.py#L410).
Kruschke favors the approach to just build the full model and let it tell you which parameters matter. You could also build variable selection into the model itself (see e.g. http://arxiv.org/pdf/math/0505633.pdf).
Finally, for a much more complete treatment, see this recent blog post: http://jakevdp.github.io/blog/2015/08/07/frequentism-and-bayesianism-5-model-selection/
Best Answer
What you are missing is that it very rarely makes sense to use a prior for which all of parameters are independent. That might make sense if the parameters were as varied and logically disconnected as, say, the set {some physical constant, some baseball player's batting average, some yeast gene's expression level}.
Most analyses look at sets of parameters that are best modeled as exchangeable, like, say, the set of all current players' batting averages or the set of all yeast genes' expression levels. For estimation, this leads to approaches like this; for testing, it leads to approaches like this.