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.

EXAMPLE

## 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 approximate the ratio as the 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 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.