Normal or lognormal? After comparing the fit of those two models by calculating a likelihood ratio (LR), how to pool LRs from several data sets?

Question

Using the method of section 6.7.2 of Burnham and Anderson, Model selection and multimodel inference: a practical information-theoretic approach, 2nd edition, one can compute the likelihood ratio comparing the likelihood that a particular data set (numerical variable; all values are positive) would occur with sampling from a normal distribution compared to sampling from a lognormal distribution. Let's say that the likelihood ratio (LR) is 3.1.

Now repeat with a second data set (a different condition in the same experiment), and that LR is 2.5.

Assume that both data sets are measurements of the same variable in similar conditions, so both come from the same distribution family. Is it appropriate to multiply those two LR values (3.1 * 2.5 = 7.75) and conclude that all the data would be 7.75 times more likely to occur with sampling from a normal distribution than with sampling from a lognormal distribution?

It sure "feels right" to multiply those likelihood ratios, but is it mathematically correct?

...

In a comment, @Sextus Empiricus asked me to clarify the context and so clarify my question.

I don't have particular data at hand, but am asking a general question. You are about to start a series of experiments. The variable you measure is always positive, and you aren't sure if the underlying distribution (population) is normal or lognormal. So you look at existing data of the same variable.

Using normality and lognormality tests doesn't quite answer the question, and would lead to lots of questions (What value for alpha?? Which normality test? How to cope with data sets that pass both normality and lognormality tests or fail both?...). Computing likelihood ratios for each data set seems (to me) to directly answer the question. For each data set, compute the ratio of the likelihood that those data would occur with sampling from a normal distribution divided by the likelihood that the data would emerge from sampling from a lognormal distribution.

You compute that likelihood ratio for each data set of each prior experiment that measures the variable of interest (in similar conditions). To directly answer the question in the comment: this would mean fitting the mean and SD when fitting the normal distribution, and the Geometric mean and geometric SD when fitting the lognormal distribution. You assume that the underlying distribution (population) is always normal or always lognormal for every data set you analyze. But the mean and SD (and GeoMean and GeoSD) will differ among data sets.

Now you have a set of likelihood ratios. How can one combine them to come up with an overall likelihood ratio for all data sets answering the question: For all the data sets, what is the ratio of the likelihood that all the data sets were sampled from a normal distribution divided by the likelihood that all the data sets were sampled from a lognormal distribution. Assuming that that overall LR is far from 1.0, it can guide analysis of data to be collected in the future -- assume normality or assume lognormality?

Answer 1

Burnham and Anderson use a Monte Carlo simulation to evaluate selecting between the normal (N) and the log-normal (LN) models. I extend this simulation to select between the normal and the log-normal when the data comes in two batches (drawn from the same data generating process — either normal or log-normal). Empirically and under these idealized conditions, the simulation suggests that multiplying likelihood ratios across datasets is a reasonable approach.

Part I. Model selection with one dataset

I start by reproducing Table 6.14 to make sure I've set up the simulation correctly and to introduce the two evaluation measures.

(a) $\pi$ is the fraction of repeated experiments when the normal model, N, is selected over the log-normal model, LN, because it has higher likelihood: $$ \begin{aligned} \operatorname{I}\left\{\frac{L_{\text{N}}}{\operatorname{L}_{\text{LN}}} > 1\right\} = \operatorname{I}\left\{\ell_{\text{N}} > \ell_{\text{LN}}\right\} \end{aligned} $$ where $L$ is the likelihood and $\ell$ is the log-likelihood.

(b) $\operatorname{E}(w)$ is the average Akaike weight for the normal model.

This concept was new to me. The Akaike weight is defined in terms of the difference in AIC with the best candidate model, ie, the model with the smallest AIC. Here we are comparing two models with the same complexity (number of parameters), so the model with the smaller AIC has the higher likelihood.

$$ \begin{aligned} w_{\text{N}} = \frac{\exp\left\{\ell_{\text{N}} - \ell_{\text{best}}\right\}}{\exp\left\{\ell_\text{N} - \ell_{\text{best}}\right\} + \exp\left\{\ell_{\text{LN}} - \ell_{\text{best}}\right\}} = \frac{1}{1 + \exp\left\{\ell_{\text{LN}} - \ell_{\text{N}}\right\}} \end{aligned} $$

While $\pi$ indicates whether we select the correct model or not, the Akaike weights scale the difference $\ell_{\text{LN}} - \ell_{\text{N}}$ to a measure between 0 and 1. While it's wrong to interpret $w_{\text{N}}$ as the "probability that the normal model is correct", the weights indicate whether it is easy or difficult to discriminate between the models. Here is example: Say the data is normal and $\ell_{\text{N}} = -11$, $\ell_{\text{LN}} = -11.2$. Since the normal likelihood is higher, we choose the normal model, correctly; its Akaike weight is 0.55. (The Akaike weight of the log-normal model is 045.) Now if we had ~10x more data and the log-likelihoods are $\ell_{\text{N}} = -110$, $\ell_{\text{LN}} = -112$, we still select the correct model by comparing the likelhoods; but now $w_{\text{N}} = 0.88$. We take the same decision but with more "certainty".

And here is my version of Table 6.14; the results are close to those in the book. It's harder to discriminate between the two models when there is less data (we can tackle this by getting more data) and when the mean is larger (the distributions are very similar, so the model choice makes little difference).

n	data	μ=5: π	μ=5: E(w)	μ=10: π	μ=10: E(w)	μ=25: π	μ=25: E(w)	μ=50: π	μ=50: E(w)
10	normal	58	54	54	51	51	50	51	50
50	normal	78	71	65	57	56	51	53	50
100	normal	87	83	72	63	59	53	54	51
500	normal	100	100	91	87	71	62	61	53
10	log-normal	38	46	44	49	48	50	49	50
50	log-normal	21	29	34	43	43	49	47	50
100	log-normal	12	18	28	37	41	47	45	49
500	log-normal	0	1	9	13	29	38	39	47

Part II. Model selection with two datasets

Now that I can run the simulation for one dataset of size $n$, it's easy to extend it to two datasets of size $n_1$ and $n_2$. For comparison with the previous results I choose the sample sizes so that $n_1 + n_2 = n$, for $n = \{25, 50, 100, 500\}$ split into ratios $9:1$, $4:1$, $2:1$ and $1:1$. [I omit $n = 10$ as it's too small to split and compute the MLEs.]

Here is how I compute the pooled measures:

(a) Exactly as suggested, $\pi$ is the fraction of repeated experiments when the normal model, N, is selected over the log-normal model, LN: $$ \begin{aligned} \operatorname{I}\left\{\frac{L_{\text{N}}^{(1)}}{\operatorname{L}_{\text{LN}}^{(1)}}\frac{L_{\text{N}}^{(2)}}{\operatorname{L}_{\text{LN}}^{(2)}} > 1\right\} = \operatorname{I}\left\{\ell_{\text{N}}^{(1)} + \ell_{\text{N}}^{(2)} > \ell_{\text{LN}}^{(1)} + \ell_{\text{LN}}^{(2)}\right\} \end{aligned} $$ where the superscripts indicate datasets $1$ and $2$.

(b) The average Akaike weight for the normal model across the two datasets is:

$$ \begin{aligned} w_{\text{N}} = \frac{\sum_{i}\exp\left\{\ell_{\text{N}}^{(i)} - \ell_{\text{best}}^{(i)}\right\}}{\sum_{i}\exp\left\{\ell_\text{N}^{(i)} - \ell_{\text{best}}^{(i)}\right\} + \exp\left\{\ell_{\text{LN}}^{(i)} - \ell_{\text{best}}^{(i)}\right\}} \end{aligned} $$

The log-likelihood is normalized separately for each dataset. So we can't compute the pooled Akaike weight from the likelihood ratio alone; we need to have the likelihood of each model separately.

I summarize the results in a small-multiples plot. There is a separate panel for each combination of sample size $n=\{25, 50, 100, 500\}$ and population mean $\mu=\{5, 10, 25, 50\}$. The dotted horizontal lines indicate the model selection performance (in terms of $\pi$ in blue and $\operatorname{E}(w)$ in red) when we analyze one dataset of size $n$. The solid lines indicate the "pooled" performance with $n$ samples split across two datasets according to the ratio given on the x-axis.

Figure 1: Selecting between the normal model (correct) and log-normal model (wrong) by pooling log-likelihood ratios between two datasets.

As Burnham and Anderson explain, $\pi$ reflects sampling variation (how often we select the correct model) and $\operatorname{E}(w)$ reflects inferential uncertainty (how easy it is to discriminate between the two models).

In this simulation at least, we don't lose much accuracy to select the correct model unless $\mu$ $=$ $5$ or the (total) sample size is $n$ $\leq$ $50$: the solid blue line is only slightly below the dotted blue line. But there is more uncertainty (big gap between the solid and dotted red lines) and it's worse to have the samples unevenly split between the two datasets. So we lose some information by splitting the data and the estimation into two; this is not surprising. In this simulation there is no order in how the datasets are generated. But in practice the datasets could be collected & analyzed sequentially. There may be little gain to collect a small followup sample and study it separately. It would be better to put the datasets together and analyze all the data we have, esp. if as stated in the question, the datasets are generated under similar conditions.

Figure 2: Selecting between the normal model (wrong) and log-normal model (correct) by pooling log-likelihood ratios between two datasets.

---

The R code to reproduce the simulation (mostly in base R) and the figure (with ggplot2).

R <- 100000

# Find the MLEs of the normal mean and variance; Compute the log-likelihood.
mle.norm <- function(x) {
  xbar <- mean(x)
  s2 <- mean((x - xbar)^2)
  sum(dnorm(x, xbar, sqrt(s2), log = TRUE))
}
# Find the MLEs of the log-normal mean and variance; Compute the log-likelihood.
mle.lnorm <- function(x) {
  lxbar <- mean(log(x))
  ls2 <- mean((log(x) - lxbar)^2)
  sum(dlnorm(x, lxbar, sqrt(ls2), log = TRUE))
}

# Generate normal data; Maximize the normal and log-normal log-likelihoods.
fit.norm <- function(n, mu) {
  sigma <- 1
  x <- rnorm(n, mu, 1)
  c(ll.a = mle.norm(x), ll.b = mle.lnorm(x))
}
# Generate log-normal data; Maximize the normal and log-normal log-likelihoods.
fit.lnorm <- function(n, mu) {
  sigma <- 1
  phi <- sqrt(mu^2 + sigma^2)
  x <- rlnorm(n, log(mu^2 / phi), sqrt(log(phi^2 / mu^2)))
  c(ll.a = mle.norm(x), ll.b = mle.lnorm(x))
}
# ---
# Do model selection using log-likelihood ratio & Akaike weights.
# NB. Since the two candidate models -- normal and log-normal -- have the same
# number of parameters, AIC(M1) - AIC(M2) = LL(M1) - LL(M2) = log{L(M1)/L(M2)}
select.model <- function(lls) {
  lls %>%
    drop_na() %>%
    mutate(
      better = ll.a > ll.b,
      Δ.a = ll.a - pmax(ll.a, ll.b),
      Δ.b = ll.b - pmax(ll.a, ll.b),
      weight = exp(Δ.a) / (exp(Δ.a) + exp(Δ.b))
    ) %>%
    summarise(
      π = mean(better),
      `E(w)` = mean(weight)
    )
}
select.model.2batches <- function(lls) {
  lls %>%
    drop_na() %>%
    mutate(
      better = (ll.a1 + ll.a2) > (ll.b1 + ll.b2),
      Δ.a1 = ll.a1 - pmax(ll.a1, ll.b1),
      Δ.b1 = ll.b1 - pmax(ll.a1, ll.b1),
      Δ.a2 = ll.a2 - pmax(ll.a2, ll.b2),
      Δ.b2 = ll.b2 - pmax(ll.a2, ll.b2),
      weight = (exp(Δ.a1) + exp(Δ.a2)) / (exp(Δ.a1) + exp(Δ.b1) + exp(Δ.a2) + exp(Δ.b2))
    ) %>%
    summarise(
      π = mean(better),
      `E(w)` = mean(weight)
    )
}
# ---
# Boilerplate code to run the simulations
simulate.norm.data <- function(n, mu, reps = R) {
  seq(reps) %>% map_dfr(~ fit.norm(n, mu))
}
simulate.lnorm.data <- function(n, mu, reps = R) {
  seq(reps) %>% map_dfr(~ fit.lnorm(n, mu))
}
simulate.norm.data.2batches <- function(n, mu, p, reps = R) {
  n1 <- round(n * p)
  n2 <- n - n1
  data1 <- simulate.norm.data(n1, mu) %>% rename_with(~ str_c(., "1"))
  data2 <- simulate.norm.data(n2, mu) %>% rename_with(~ str_c(., "2"))
  bind_cols(data1, data2)
}
simulate.lnorm.data.2batches <- function(n, mu, p, reps = R) {
  n1 <- round(n * p)
  n2 <- n - n1
  data1 <- simulate.lnorm.data(n1, mu) %>% rename_with(~ str_c(., "1"))
  data2 <- simulate.lnorm.data(n2, mu) %>% rename_with(~ str_c(., "2"))
  bind_cols(data1, data2)
}
run.simulations.save.results <- function(seed = 1234) {
  set.seed(seed)
  # Table 6.14A
  grid <- expand_grid(
    n = c(10, 25, 50, 100, 500),
    μ = c(5, 10, 25, 50)
  )
  sim.norm <- grid %>%
    mutate(
      data = pmap(list(n, μ), simulate.norm.data),
      table = map(data, select.model)
    )
  sim.lnorm <- grid %>%
    mutate(
      data = pmap(list(n, μ), simulate.lnorm.data),
      table = map(data, select.model)
    )
  bind_rows(
    norm = sim.norm,
    lnorm = sim.lnorm,
    .id = "data"
  ) %>%
    unnest(
      table
    ) %>%
    write_csv(
      "600592-Table-6.14A.csv"
    )
  # Table 6.14B
  grid <- expand_grid(
    n = c(25, 50, 100, 500),
    μ = c(5, 10, 25, 50),
    p = c(1 / 10, 1 / 5, 1 / 3, 1 / 2)
  )
  sim.norm <- grid %>%
    mutate(
      data = pmap(list(n, μ, p), simulate.norm.data.2batches),
      table = map(data, select.model.2batches)
    )
  sim.lnorm <- grid %>%
    mutate(
      data = pmap(list(n, μ, p), simulate.lnorm.data.2batches),
      table = map(data, select.model.2batches)
    )
  bind_rows(
    norm = sim.norm,
    lnorm = sim.lnorm,
    .id = "data"
  ) %>%
    unnest(
      table
    ) %>%
    write_csv(
      "600592-Table-6.14B.csv"
    )
}

# 6.7.2 A Normal Versus Log-Normal Example
library("tidyverse")

plot_selection <- function(data) {
  table.b %>%
    filter(data == {{ data }}) %>%
    ggplot(
      aes(
        p, value,
        group = name,
        color = name
      )
    ) +
    geom_line() +
    geom_hline(
      aes(
        yintercept = value,
        color = name
      ),
      linetype = 2,
      data = table.a %>% filter(data == {{ data }})
    ) +
    facet_wrap(
      ~ n + μ,
      scales = "free_y",
      labeller = "label_both"
    ) +
    scale_x_continuous(
      breaks = c(1 / 10, 1 / 5, 1 / 3, 1 / 2),
      labels = c("1:9", "1:4", "1:2", "1:1")
    ) +
    theme(
      axis.title = element_blank(),
      legend.title = element_blank(),
      panel.spacing.y = unit(0, "line")
    )
}

run.simulations.save.results()

table.a <- read_csv("600592-Table-6.14A.csv")
table.b <- read_csv("600592-Table-6.14B.csv")

table.a %>%
  pivot_wider(
    names_from = μ,
    names_glue = "{μ},{.value}",
    values_from = c(π, `E(w)`),
    values_fn = ~ round(100 * .)
  ) %>%
  select(
    n, data,
    starts_with("5,"),
    starts_with("10,"),
    starts_with("25,"),
    starts_with("50,")
  ) %>%
  knitr::kable()

table.a <- table.a %>%
  filter(n >= 25) %>%
  pivot_longer(c(π, `E(w)`))
table.b <- table.b %>%
  filter(n >= 25) %>%
  pivot_longer(c(π, `E(w)`))

plot_selection("norm")
plot_selection("lnorm")

Answer 2

Case of simple hypotheses

Yes the likelihood for independent events is multiplicative.

$$\mathcal{L}(\theta ; x_1, x_2)= \mathcal{L}(\theta ; x_1) \cdot \mathcal{L}(\theta ; x_2)$$

And as a consequence the likelihood ratio is multiplicative as well.

$$\frac{\mathcal{L}(\theta_1 ; x_1, x_2)}{\mathcal{L}(\theta_2 ; x_1, x_2)}= \frac{\mathcal{L}(\theta_1 ; x_1)}{\mathcal{L}(\theta_2 ; x_1)} \cdot \frac{\mathcal{L}(\theta_1 ; x_2)}{\mathcal{L}(\theta_2 ; x_2)}$$

---

In the case of log likelihood and log of the likelihood ratio, the multiplicative property changes into an additive property. $\log(ab) = \log(a) + \log(b)$

---

Be aware that the multiplication is only true for independent events. Multiple measurements are not always independent. For instance, the second measurement can be made as a consequence of the first measurements, and this leads to the stopping problem.

Case of a composite alternative hypothesis 1

It might be that you are using the likelihood ratio test for testing a null hypothesis $H_0: \theta = \theta_0$ against an alternative composite hypothesis $H_a: \theta \neq \theta_0$ and you are using the likelihood ratio

$$\lambda_x = \frac{\mathcal{L}(\theta_0;x)}{\mathcal{L}(\hat\theta_{ML};x)}$$

Here you can again multiply the likelihood ratio values (or sum the log values), but you need to consider that you computed two times parameters and the distribution of the ratio will be different.

Typically the distribution of the logarithm of the likelihood ratio, multiplied with -2, is considered to be (approximately) chi-squared distributed with $\nu$ degrees of freedom (where $\nu$ is the number of free parameters that are used to fit) and now you will have the sum of two such values which will be chi-squared distributed with $2\nu$ degrees of freedom.

Using the product of the likelihood ratios is similar to Fisher's method of combining p-values. (Which involves multiplication of the p-values and comparing with a chi-squared distribution)

Case of a composite alternative hypothesis 2

In your explanation you have a composite hypothesis for both hypotheses. You can still multiply likelihood functions and ratios. You can consider the two different datasets as if it was a single dataset that you fit with twice as many parameters.

A question is however what use your likelihood ratio has. Often you use it for a hypothesis test and you can know some distribution for the likelihood distribution given that null hypothesis is true. In your example, hypotheses of lognormal versus normal, it may not be easy to figure out this distribution. But, this is a problem that is unrelated to combining ratios and already occurs with a single dataset.

and conclude that all the data would be 7.75 times more likely to occur with sampling from a normal distribution than with sampling from a lognormal distribution?

This is not the interpretation of the likelihood ratio when it is computed with composite hypotheses. In the case of composite hypotheses there are nuisance parameters where you fit the cases that maximise the likelihood.

Compare for instance with a Bayesian analysis where those parameters could have a prior distribution. In that case you would have a likelihood ratio that can be interpreted as the relative probabilities for the data to occur with the two different hypotheses.

Case of information interpretation.

In the comments a situation is described that seems like model selection and one might use the likelihood as part of the Akaike information criterion (AIC).

In this case, AIC, it might make sense to add up the average log likelihood, the log likelihood divided by the number of observations. It is the average log likelihood that approaches the cross entropy and it might make more sense to add up the entropy rather than the likelihood (which differ by a constant, the number of observations).

If the number of observations are different for the different datasets then using this way of combining the likelihoods will be of influence.