Using the Geyer-Thompson algorithm, you can compute the ratio of two likelihoods (needed to do a likelihood ratio test or to compute an information criterion like AIC or BIC). This is useful when comparing two hierarchical models using data cloning, which uses the Bayesian machinery to compute ML estimates. However, you need to compute the average of the ratio on the data rather than the log scale.
## Generate some data and create a linear model to play with
set.seed(123)
n <- 1000
beta0 <- 20
beta1 <- 5.92
beta2 <- -1.2
x <- runif(n=n, min=0, max=5)
mu <- beta0 + beta1 * x + beta2 * x**2
sd <- 2
y <- rnorm(n=n, mean=mu, sd=sd)
data <- data.frame(x=x, y=y)
# plot(data$x, data$y)
m1 <- glm(formula = y ~ poly(x,1), family=gaussian(link="identity"), data = data)
m2 <- glm(formula = y ~ poly(x,2), family=gaussian(link="identity"), data = data)
m3 <- glm(formula = y ~ poly(x,3), family=gaussian(link="identity"), data = data)
bbmle::ICtab(m1,m2,m3,mnames=c("linear", "quadratic", "cubic"), type = "AIC",
weights = T, delta = T, base = T, logLik = T)
# logLik AIC dLogLik dAIC df weight
# quadratic -2112.7 4233.4 414.8 0.0 4 0.58
# cubic -2112.0 4234.0 415.5 0.6 5 0.42
# linear -2527.5 5060.9 0.0 827.6 3 <0.001
# sjPlot::plot_model(m2, type = "pred", show.data = T)
Using data cloning, you compute the ratio of m1 and m2 on the log scale (pp 358 in Ponciano et al 2009) and then take the average of the exponent of this ratio (mean(exp(lnm1 - lnm2))), i.e.,
#### LRT
likRATIO1 <- exp(logLik(m1)-logLik(m2)) |> as.vector()
likRATIO2 <- exp(logLik(m1)-logLik(m3)) |> as.vector()
-2* log( likRATIO1 ) + 2*( 3-4 ) ## m1 vs m2
-2* log( likRATIO2 ) + 2*( 3-5 ) ## m1 vs m3
-2* log( 1/(likRATIO1/likRATIO2) ) + 2*( 4-5 ) ## m2 vs m3
m1$aic - m2$aic # [1] 827.5576
m1$aic - m3$aic # [1] 826.9519
m2$aic - m3$aic # [1] -0.6057001
When the difference in the log-likelihood is a large negative number, the exponent is zero and the likelihood is -Inf, e.g.,
> -2* log( likRATIO1 ) + 2*( 3-4 ) ## m1 vs m2
[1] Inf
> -2* log( likRATIO2 ) + 2*( 3-5 ) ## m1 vs m3
[1] Inf
> -2* log( 1/(likRATIO1/likRATIO2) ) + 2*( 4-5 ) ## m2 vs m3
[1] NaN
Is there a trick here to computing the log( . ) that doesn't result in numerical issues?
Thompson, E. A. 1994. Monte Carlo likelihoods in genetic mapping. Statistical Science 9:355–366.
Lele, S. R., Dennis, B., & Lutscher, F. (2007). Data cloning: Easy maximum likelihood estimation for complex ecological models using Bayesian Markov chain Monte Carlo methods. Ecology Letters, 10(7), 551–563. https://doi.org/10.1111/j.1461-0248.2007.01047.x
Ponciano, J. M., Taper, M. L., Dennis, B., & Lele, S. R. (2009). Hierarchical models in ecology: confidence intervals, hypothesis testing, and model selection using data cloning. Ecology, 90(2), 356–362.
You haven't quite given us an example of your problematic case, but suppose we want to compute mean(exp(lnm1 - lnm2)) when lnm1 and lnm2 are vectors with some pairs of elements very different from each other:
lnm1 <- c(0, 1000, 0)
lnm2 <- c(1000, 0, 1000)
exp(lnm1-lnm2) ## 0 Inf 0
matrixStats::logSumExp(lnm1-lnm2) ## 1000 = log(sum(exp(lnm1-lnm2)))
In this case the answer is sort of trivial because the second element will dominate the first and third. We can check that this works correctly when the values are computable either way:
x <- 1:10
y <- 10:1
all.equal(matrixStats::logSumExp(x-y), log(sum(exp(x-y)))) ## TRUE
(Computing the mean is left as an exercise ...)
There are many packages that include equivalents of logsumexp (e.g. limma::logsumexp, MatrixGenerics::<various>, mclust::logsumexp)