How exactly does one marginalize over parameters in an N-dimensional likelihood? I see no equations for the following, so I'm not sure exactly what they are talking about:
"For each model, we determine the best fit parameters from
the peak of the N-dimensional likelihood surface. For each parameter in the model we also compute its one dimensional likelihood function by marginalizing over all other parameters."
How do you obtain a one-dimensional likelihood function by marginalizing over all other parameters? 
 A: You should always reference your quotes.
The only reference on the Internet I found with this quote 

"For each model, we determine the best fit parameters from the peak of
  the N-dimensional likelihood surface. For each parameter in the model
  we also compute its one dimensional likelihood function by
  marginalizing over all other parameters." (p.3)

is an cosmology arxiv paper on WMAP by Spergel et al. (2003). If you turn the page after this quote, you will find an equation defining the expectation under the "marginal likelihood" as
$$<α_i> = \int d^N α\mathcal{L}(α)α_i\,.$$ This means that the "marginal likelihood" is
$$\int d^{N-1} α_{-i}\mathcal{L}(α)\,.$$
If I quote from the previous sentences of the paper as well

"For each model studied in the paper, we use a Monte Carlo Markov Chain to explore the likelihood surface. We assume flat priors in our basic parameters, impose positivity constraints on the matter and baryon
  density (these limits lie at such low likelihood that they are
  unimportant for the models. We assume a flat prior in τ , the optical
  depth, but bound τ < 0.3. This prior has little effect on the fits but
  keeps the Markov Chain out of unphysical regions of parameter space. (p.3)"

this means that the authors take a Bayesian stance with a flat prior in their parametrisation. 
A: Say you have the state of Information $I$, some observations $\{y_i\}$ and some parameters $\theta_p$ where $p \in \{1,2, \ldots n\}$ then in the continuous case you get the marginal likelihoods from the joint likelihood $p(\theta_1,\theta_2, \ldots \theta_n |\{y_i\},I)$ by integration:
\begin{align*}
p(\theta_k|\{y_i\},I) &=\int_{-\infty}^{+\infty} \ldots \int_{-\infty}^{+\infty} p(\theta_1,\theta_2, \ldots \theta_n |\{y_i\},I)\,d\theta_{p_1} \cdots \,d\theta_{p_{n-1}} 
\end{align*}
where $p_j \in \{1,2, \ldots n\}\setminus k$
In the discrete case you will have sums instead of integrals as you can see here.
