5
$\begingroup$

I want to compare two profile likelihood curves and determine if they are significantly different from one another.

For example are the following curves significantly different from one another:

two profile likelihood curves

I realize I can find a 95% confidence interval for a value from a profile likelihood curve using likelihood ratio testing like so. With this method I can determine if a single value is outside of the 95% confidence range of one profile likelihood curve. But I want to compare two curves.

|cite|improve this question
$\endgroup$
  • $\begingroup$ What is your purpose in comparing the curves? What is the substantial question? $\endgroup$ – kjetil b halvorsen Feb 28 '15 at 10:59
  • $\begingroup$ @kjetilbhalvorsen I want to to know if the estimate of theta from my first model (black) is significantly different from the estimate of theta from my second model (red). $\endgroup$ – Ben Haley Mar 1 '15 at 15:52
  • $\begingroup$ Why not K-S test or AIC? $\endgroup$ – corey979 Jun 6 '15 at 22:46
  • $\begingroup$ @corey979 the K-S test looks very promising, thanks! Not sure how I would apply AIC to this problem... $\endgroup$ – Ben Haley Jun 7 '15 at 23:02
  • $\begingroup$ @BenHaley Ah, sorry, I thought you have an underlying dataset. So not AIC. $\endgroup$ – corey979 Jun 9 '15 at 6:13
2
$\begingroup$

So I tried to solve this myself.

The general strategy

  1. Convert the profile likelihood curves into probability distribution functions by exploiting the likelihood ratio test as described here.
  2. 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.

enter image description here

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.

|cite|improve this answer
$\endgroup$
  • 4
    $\begingroup$ Your solution does not appear to address the substantive question, nor does it appear to have any theoretical support. As far as I can tell, you are evaluating the chance that a random variable drawn from the distribution corresponding to the peak of one profile curve would exceed a random variable drawn from the distribution corresponding to the peak of another profile curve. What justifies using that as a p-value? It's not one. $\endgroup$ – whuber♦ Oct 22 '14 at 18:44
  • $\begingroup$ I'm not trying to draw from the peak of the likelihood profile, but anywhere along the likelihood profile with probability proportional to the how probable that value is given the likelihood. This is why I convert the the profile likelihood into a cdf for that number. I guess that's not clear. If you tell me what made you think I was only looking at peaks I can try to clarify. $\endgroup$ – Ben Haley Oct 22 '14 at 19:36
  • 2
    $\begingroup$ A likelihood does not describe a distribution of its parameters--not by any means: stats.stackexchange.com/questions/2641/…. $\endgroup$ – whuber♦ Oct 22 '14 at 19:38
  • $\begingroup$ Yeah that is true generally. But is is possible to derive a confidence boundary for profile likelihoods. For example [1]. And hence I propose we can derive the full distribution function. Are you talking about profile likelihoods in particular or likelihoods in general? [1]: stats.stackexchange.com/a/78361/7494 $\endgroup$ – Ben Haley Oct 22 '14 at 19:53
  • 2
    $\begingroup$ I am still baffled by why treating a profile likelihood as if it were a CDF (of what random variable?) ought to produce a quantity that can be interpreted as a p-value. I'm not saying it isn't, but only that the connection hasn't been made. A p-value would reflect how the profile likelihoods would vary under the null hypothesis that both were generated by the same underlying distribution (and therefore differed only by chance). It would be based on a statistic that quantifies the discrepancies between two such curves and involves a Neyman-Pearson critical region. $\endgroup$ – whuber♦ Oct 22 '14 at 20:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.