2
$\begingroup$

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?

$\endgroup$
6
  • $\begingroup$ Do you propose to fit the same models separately on the two datasets? In that case, wouldn't you expect that the MLEs of the model parameters are (at least somewhat) different as well? So the first likelihood ratio would compare models $M_A(\hat{\theta}_1)$ vs $M_B(\hat{\phi}_1)$ and the second likelihood ratio would compare $M_A(\hat{\theta}_2)$ vs $M_B(\hat{\phi}_2)$. $\endgroup$
    – dipetkov
    Jan 4 at 1:25
  • $\begingroup$ @dipetkov. You understand right. But if both data sets (or more generally, all data sets) are of the same variable in not-too-different conditions, you'd expect that one model is correct for both (all). I'm trying to figure out how to decide when to assume sampling from normal vs. lognormal. You'd need to look back at old (or preliminary) data, and see which model fits better, and then use that model for future similar studies. $\endgroup$ Jan 4 at 1:35
  • $\begingroup$ Would stock returns be a good example? Eg each data set would have a year of data, and each datum would be the S&P 500 price divided by the price on the previous day. $\endgroup$
    – Matt F.
    Jan 6 at 15:37
  • $\begingroup$ @MattF. Might be a good example. But.. A better example would be weight (or glucose concentration in the blood) of a bunch of mice. In your example, the values follow some distribution, but you collect them sequentially which adds another variable. $\endgroup$ Jan 6 at 17:26
  • $\begingroup$ @MattF. Let's make it simple. Your plan is to compare weights in mice given one of two diets. One measurement per mouse after 2 months on the diet. To plan your analysis you need to know if weight is closer to normal or lognormal. So you look back at some prior data from several prior studies. What calculations do you do to pool the data in order decide (in general) whether mouse weight is closer to normal or lognormal? $\endgroup$ Jan 6 at 20:25

2 Answers 2

1
$\begingroup$

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. Normal data

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. long-normal data


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")
$\endgroup$
4
  • $\begingroup$ Thanks. Just to clarify. For the normal distributions, you always choose a population standard deviation of 1.0. So the coefficient of variations runs from .2 to 0.02. For the lognormal distributions, the SD of the logarithms is 1.0, so the Geometric standard deviation is always exp(1) or 2.72. I understand the blue lines (solid and dotted) in your figures, but can't really understand the red lines. What are they telling us? $\endgroup$ Jan 10 at 15:32
  • $\begingroup$ I set up the simulation exactly as in Table 6.14 as the only thing I wanted to vary is the data split. From the text: "...select between the simple normal and log-normal models, variance(y) = 1 for both distributions". Of course, I also put all the code up, so feel free to vary $\sigma^2$ if that's of interest. $\endgroup$
    – dipetkov
    Jan 10 at 16:26
  • $\begingroup$ Here is my interpretation of $E(w)$ (in red) vs $\pi$ (in blue). To calculate $\pi$ you only look at the selection decision (was the normal model or the log-normal the one selected?). To calculate $E(w)$ you look at the difference between the normal and log-normal likelihoods, so you get an idea of how "difficult" it was to choose between the two models. If the likelihood of the correct model is much higher than the wrong model, it will also have higher weight (than 0.5). So the red lines are telling us that we do lose information about the model when we split the data. Hardly a surprise. $\endgroup$
    – dipetkov
    Jan 10 at 16:26
  • $\begingroup$ PS. By "lose information" I mean the gap between the solid and dotted red lines, vs. the much smaller, or in some cases, non-existent gap between the solid and dotted blue lines. $\endgroup$
    – dipetkov
    Jan 10 at 16:32
2
$\begingroup$

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.

$\endgroup$
6
  • $\begingroup$ I think this question is about multiplying the calculated likelihood ratios after maximizing the likelihoods separately on each dataset. In your answer you have $L(\theta;x_1)$ and $L(\theta;x_2)$, which are functions of $\theta$, while in the question the numbers 3.1 & 2.5 refer to $L(\hat{\theta}_1;x_1)$ and $L(\hat{\theta}_2;x_2)$ where the $\hat{\theta}$s are the dataset-specific MLEs. $\endgroup$
    – dipetkov
    Jan 4 at 10:56
  • $\begingroup$ @dipetkov clarified my question correctly. Does that change the answer? $\endgroup$ Jan 4 at 12:10
  • $\begingroup$ @HarveyMotulsky that makes that I don't understand your question anymore. What hypothesis are you testing with the likelihood ratio test if you consider two different $\hat\theta$ for the two data sets? $\endgroup$ Jan 4 at 12:46
  • $\begingroup$ Is the $-2\log$ x the likelihood ratio (approximately) Chi-squared? The two models are Normal and Log-Normal, so they are not nested? Burnham & Anderson use the likelihood ratio to select between the models (or equivalently, the AIC since the models have the same number of parameters). But they don't test a null hypothesis. $\endgroup$
    – dipetkov
    Jan 4 at 15:23
  • $\begingroup$ @dipetkov. I've lengthened my original question substantially to clarify precisely what my question is. $\endgroup$ Jan 4 at 16:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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