So I tried to solve this myself.
The general strategy
- Convert the profile likelihood curves into probability distribution functions by exploiting the likelihood ratio test as described here.
- Sample points from within these probability distributions and determine if one is consistently higher or lower.
The functions to do this in R are below. I'd be happy to clarify things if someone runs across this question and prods me a bit.
Generate a cdf from a likelihood profile
# Generate cdf
# Given a set of likelihood return the corresponding cumulative
# density function
cdf <- function(likelihoods) {
cdf <- (1 - pchisq(log(max(likelihoods)) - log(likelihoods), 1)) / 2
# After 0.5 (the most likely) cdf values should count up to
# 1 instead of counting back down to 0
greater_than_most_likely <- which.max(cdf):length(cdf)
cdf[greater_than_most_likely] <- 1 - cdf[greater_than_most_likely]
cdf
}
For the example I showed, this results in cdf curves as follows.

Determine the probability that two cdfs are different
# p greater
# Given two cumulative probability density curves estimate the p value
# associated with the hypothesis
# a > b
#
# Tests
# p_greater(c(0.000, 0.001), c(0.999, 1))
# # ~ 0
# p_greater(pnorm((-100:100)/100), pnorm((-100:100)/100))
# # ~ 0.50
# p_greater(pnorm((-100:100)/100, sd=2), pnorm((-100:100)/100, sd=1))
# # ~ 0.50
# p_greater(pnorm((-1000:1000)/100, mean=1), pnorm((-1000:1000)/100, mean=0))
# # ~ 1 - sum(rnorm(1000000, mean=1) >= rnorm(1000000, mean=0)) / 1000000
# # ~ 0.24
p_greater <- function(cdf_a, cdf_b, n_samples=10000) {
a_sample <- which.closest(runif(n_samples), cdf_a)
b_sample <- which.closest(runif(n_samples), cdf_b)
a_greater <- a_sample > b_sample
a_equal <- a_sample == b_sample
1 - (sum(a_greater) + sum(a_equal)/2) / n_samples
}
# p different
# Like p greater but a two sided test
p_different <- function(cdf_a, cdf_b, n_samples=10000) {
p_a_greater <- p_greater(cdf_a, cdf_b, n_samples)
p_b_greater <- p_greater(cdf_b, cdf_a, n_samples)
min(p_a_greater, p_b_greater) / 2
}
so that I can call:
p_different(cdf(likelihoods_1), cdf(likelihoods_2))
In this case I got a p-value of ~ 0.16.