---
title: "A Course on Bioequivalence: Unit 2 - Log transformed data, geometric means, coefficient of variation, and the log-normal distribution"
execute:
  error: false
  echo: true
  enabled: true
  daemon: 300
engine: julia
author:
  - Yoni Nazarathy
  - Andreas Noack
  - Vijay Ivaturi
format:
  html:
    include-in-header:
      text: |
        <style>
          .header-container {
            width: 100%;
          }
          .logo {
            display: block;
            width: 150px;
            margin-left: -2.5%;
          }
        </style>
        <style>
          .no-arrow::before {
            content: '';
          }
          .no-arrow > summary {
            list-style: none;
          }
          .no-arrow > summary::-webkit-details-marker {
            display: none;
          }
        </style>
        <link rel = "shortcut icon" href = "https://pumas.ai/favicon.ico" />
    include-before-body:
      text: |
        <div class="header-container">
          <img src="https://pumas-assets.s3.amazonaws.com/CompanyLogos/PumasAI/RGB/PNG/Pumas+AI+Primary%404x.png" alt="Pumas Logo" class="logo">
        </div>
    include-after-body:
      text: |
    css: style.css
    self-contained-math: true
    anchor-sections: true
    theme: default
    toc: true
    toc-depth: 4
    toc-expand: 2
    toc-location: left
    toc-title: Contents
    number-sections: true
    number-depth: 2
    code-summary: Show/Hide Code
bibliography: ../references.bib
fig-format: svg
fig-width: 8
fig-height: 6
license: CC BY-SA 4.0
profile-type: 
  no-intermediate-text: false #set to true to hide the intermediate text.
---

This is **Unit 2** of a [complete bioequivalence analysis course using Pumas](https://tutorials.pumas.ai/#module-16-a-course-on-bioequivalence).
There are **15 units** in the course in total.

## Unit overview

In the [previous unit](unit-1-what-is-be.html) we learned what is a bioequivalence study.
In this unit we lay some of the basic groundwork with regards to quantities and tools of interest.
In particular we consider:

  - Log transformed data
  - The log-normal distribution
  - Geometric means
  - The geometric mean ratio (GMR)
  - Bounds on the GMR
  - Coefficient of variation
  - Maximum likelihood estimation for normal and log-normal data

::: {.content-hidden when-meta="profile-type.no-intermediate-text"}
In this unit we use the following packages, of which the first 4 are Pumas packages (available only with a Pumas installation).
:::

```{julia}
using Pumas # available with Pumas products - also reexports Distributions
using NCA # available with Pumas products
using Bioequivalence  # available with Pumas products
using PharmaDatasets # available with Pumas products
using AlgebraOfGraphics
using CairoMakie
using SummaryTables
using DataFramesMeta
using StatsBase
```

::: {.content-hidden when-meta="profile-type.no-intermediate-text"}
We have used some of these packages in the previous unit. New to this unit is the `StatsBase` package which provides additional basic statistical functionality (not available in the in-built `Statistics` package). Similarly, the `Distributions` package provides support of all common probability distributions. While we don't explicitly have a `using Distributions` statement, that package is re-exported by `Pumas`.

:::{.callout-warning}
This unit has a bit more equations and algebra than other units. If you find the notation overwhelming, do not let that deter you from advancing. Keep in mind that some of the mathematical details are only there for further insight and may be skipped on a first reading.
:::
:::

## From Normal to Log-Normal

Central to statistics is the [normal distribution](https://en.wikipedia.org/wiki/Normal_distribution) (or Gaussian distribution). Mathematically it is a probability distribution with density,

$$
f_N(u) = \frac{1}{\sqrt{2 \pi\sigma^2}} e^{-\frac{(u-\mu)^2}{2\sigma^2}}.
$$

The parameters $\mu$ and $\sigma^2$ are the mean and the variance respectively.

::: {.content-hidden when-meta="profile-type.no-intermediate-text"}
A random value distributed according to this distribution can take any value, positive, or negative. Practically it is almost certain to fall in the interval $[\mu - 3 \sigma, \mu + 3 \sigma]$.

Many statistical inference models and procedures assume that the data (or some component of the data) is distributed according to a normal distribution, often with $\mu$ and $\sigma^2$ unknown, and estimated. This includes many of the bioequivalence procedures.

In a statistical model, one almost always makes assumptions about distributions. However, certain quantities, such as endpoints, Cmax, and AUC, **are not** distributed according to a normal distribution and thus making such assumptions won't work. That is for such endpoints, assuming a normal distribution would typically be a highly invalid assumptions.

For example let us look at the `R` and `T` data from some arbitrary dataset (this happens to be a dataset of a four period fully replicate design - yet we simply separate out `T` and `R`). Here we fit [kernel density](https://en.wikipedia.org/wiki/Kernel_density_estimation) plots to the distribution to get a rough feel for the shape of the distribution. The subtle wiggles in the curve are simply because we are estimating the distribution. But importantly the distributions only take on values on the positive numbers and are very **asymmetric**.
:::

```{julia}
df = dataset(joinpath("bioequivalence", "RTTR_TRRT", "SLTGSF2020_DS16"))
@rtransform!(df, :treatment = :sequence[:period])
specs =
    data(df) * mapping(:PK, color = :treatment) * AlgebraOfGraphics.density(bandwidth = 0.5)
draw(specs, axis = (; limits = ((-1, 5), nothing)))
```

::: {.content-hidden when-meta="profile-type.no-intermediate-text"}
Note that `PK` just stands for an arbitrary endpoint like Cmax or AUC. Also note that we plot the `R` and `T` data separably.

It is clear from the above that the distributions are only of positive values, and that they are asymmetric. In particular these distributions are far from being normal distributions.
:::

As a remedy for the situation, we carry out a **log transformation** of the endpoint data.

::: {.content-hidden when-meta="profile-type.no-intermediate-text"}
With this remember that the logarithm function takes values that are strictly positive and returns values that are either positive or negative. With such a log transformation, the shape of the distribution changes and is much more similar to a normal distribution. Note that the base of the logarithm we use is not of great importance because logarithms of different bases are the same up to a constant factor. Still it is common to use the natural logarithm. Here we just use Julia's `log` function which is the natural logarithm, sometimes denoted in mathematics via $\ln$, but in our case denoted via $\log$.
:::

```{julia}
@rtransform!(df, :logPK = log(:PK))
specs =
    data(df) *
    mapping(:logPK, color = :treatment) *
    AlgebraOfGraphics.density(bandwidth = 0.5)
draw(specs, axis = (; limits = ((-3, 3), nothing)))
```

::: {.content-hidden when-meta="profile-type.no-intermediate-text"}
Mathematically this transformation can be viewed as a relationship between two random quantities (random variables), $X_{L}$ and $X_N$, where $X_{L}$ represents the original endpoint and $X_N$ represents the log transformed endpoint. Specifically,
:::

$$
X_N = \log(X_L),
\qquad \qquad
\text{or similarly}
\qquad \qquad
e^{X_N} = X_L.
$$

::: {.content-hidden when-meta="profile-type.no-intermediate-text"}
Here it is $X_N$, the transformed data, which is assumed to have a normal distribution with density as stated above and is determined its mean $\mu$ and standard deviation $\sigma$. However, $X_L$, the original observed endpoint values, are distributed differently, they have a probability density function,
:::

$$
f_L(u) = \frac{1}{u\sqrt{2 \pi\sigma^2}} e^{-\frac{(\log(u)-\mu)^2}{2\sigma^2}},
\qquad
\text{for positive}~u.
$$

This is the [lognormal distribution](https://en.wikipedia.org/wiki/Log-normal_distribution).

::: {.content-hidden when-meta="profile-type.no-intermediate-text"}
Hence $X_L$ denotes "lognormal" while $X_N$ denotes "normal". Note that the derivation of the density $f_L(u)$ follows from that the density $f_N(u)$ based on some probability theory and elementary calculus. We omit the details.

As an illustration of the relationship between $X_N$ (normal) and $X_L$ (lognormal), consider this numerical experiment on synthetic data.
:::

```{julia}
μ, σ = 0.4, 1.2
X_N_points = rand(Normal(μ, σ), 10^6)
X_L_points = exp.(X_N_points)

df = DataFrame(
    value = vcat(X_N_points, X_L_points),
    variable = vcat(
        fill("X_N_points", length(X_N_points)),
        fill("X_L_points", length(X_L_points)),
    ),
)

plt =
    data(df) *
    mapping(:value, color = :variable) *
    AlgebraOfGraphics.density(npoints = 10^4, bandwidth = 0.1)
draw(plt, axis = (; limits = ((-3, 8), nothing)))
```

::: {.content-hidden when-meta="profile-type.no-intermediate-text"}
You can see the distribution of the normal random values ($X_N$) and compare it with the transformed log normal values ($X_L$). Note that in an endpoint data bioequivalence setting we go the other way around, we convert from $X_L$ to $X_N$.

Note that by transforming from $X_N$ to $X_L$ we get a different mean and variance. That is for the random quantity $X_L$, the parameters $\mu$ and $\sigma^2$ are no longer the mean and variance, but are rather just parameters of the distribution. This is evident by considering the sample mean and sample standard deviation as such:
:::

```{julia}
@show mean(X_N_points), std(X_N_points)
@show mean(X_L_points), std(X_L_points);
```

::: {.content-hidden when-meta="profile-type.no-intermediate-text"}
Mathematically, one can compute the mean and variance of the log normal distribution, in particular, these formulas hold:
:::

$$
{\mathbb E}[X_L] = e^{\mu + \frac{\sigma^2}{2}},
\qquad
\text{and}
\qquad
\text{Var}(X_L) = (e^{\sigma^2}-1)e^{2\mu + \sigma^2}.
$$

::: {.content-hidden when-meta="profile-type.no-intermediate-text"}
You can see that these theoretical values are very much near the estimated values (from the synthetic data):
:::

```{julia}
@show exp(μ + σ^2 / 2)
@show √((exp(σ^2) - 1) * exp(2μ + σ^2));
```

::: {.content-hidden when-meta="profile-type.no-intermediate-text"}
Hence if we have **endpoint data** distributed as $X_L$ (log-normal), we can transform it to $X_N$ (normal), carry out statistical inference to estimate $\mu$ and $\sigma$, with the estimates denoted $\hat{\mu}$ and $\hat{\sigma}$. Then if we wanted to, we would get back to the log-normal distribution with these formulas.

In practice, considering a representation using Geometric means, and the coefficient of variation is common. We outline this approach now.
:::

## Geometric Means, Arithmetic Means, and Logs

::: {.content-hidden when-meta="profile-type.no-intermediate-text"}
In general, when we have a collection of values, $u_1,\ldots,u_n$, there many ways summarize these values via a measure of centrality (or average). In our case, these two ways are important:
:::

$$
\text{mean}(u_1,\ldots,u_n) = \frac{1}{n} \sum_{i=1}^n u_i,
\qquad
\qquad
\text{geomean}(u_1,\ldots,u_n) = \sqrt[n]{\prod_{i=1}^n u_i}.
$$

Note that "mean" here is just short for the **arithmetic mean** and "geomean" is short for the **geometric mean** (the big symbol $\sum$ is for summation and the big symbol $\prod$ is for product).

::: {.content-hidden when-meta="profile-type.no-intermediate-text"}
Each of these means is suitable for their own context. Arithmetic means, the most popular, are good for additive quantities and geometric means are good for multiplicative quantities. (Note also that with geometric means we typically assume that $u_i > 0$ while for arithmetic means there is no such restriction.)

Here is an example illustrating the use of each of these means:
:::

```{julia}
println("Arithmetic mean is good for totals:")
weights = [10, 20]
@show mean(weights)
total_weight = mean(weights) * 2
@show total_weight
@show sum(weights)

println("\nGeometric mean is good for rates:")
rates_of_growth_per_year = [1.1, 1.2]
@show geomean(rates_of_growth_per_year)
percentage_growth_after_two_years = geomean(rates_of_growth_per_year)^2
@show percentage_growth_after_two_years
@show prod(rates_of_growth_per_year);
```

Now in our case, we have **endpoint data** (e.g. AUC). Let's denote it as,

$$
\text{endpoint data} = \{x_1,\ldots,x_n\}.
$$

When we transform the data with a log transformation, we have the data,
$$
\text{log transformed endpoint data} = \{\log(x_1),\ldots, \log(x_n)\}.
$$

::: {.content-hidden when-meta="profile-type.no-intermediate-text"}
Now look at the arithmetic mean of the log transformed endpoint data:
:::

$$
\begin{align*}
\text{mean}(\log(x_1),\ldots, \log(x_n)) &= \frac{1}{n} \sum_{i=1}^n \log(x_i) \\
&= \frac{1}{n} \log \Big(\prod_{i=1}^n x_i  \Big)\\
&= \log \Big(\sqrt[n]{\prod_{i=1}^n x_i} \Big)\\
&= \log \Big(\text{geomean}(x_1,\ldots,x_n)\Big).
\end{align*}
$$

::: {.content-hidden when-meta="profile-type.no-intermediate-text"}
So we see:
:::

**The arithmetic mean of the log transformed data is the log of the geometric mean of the original data.**

::: {.content-hidden when-meta="profile-type.no-intermediate-text"}
Similarly we have,
:::

$$
e^{\text{mean}(\log(x_1),\ldots, \log(x_n))} = \text{geomean}(x_1,\ldots,x_n).
$$

::: {.content-hidden when-meta="profile-type.no-intermediate-text"}
Hence we also see:
:::

**The geometric mean of the original data is the exponent of the arithmetic mean of the log transformed data.**

::: {.content-hidden when-meta="profile-type.no-intermediate-text"}
Now from a statistical perspective, because the log transformed data of endpoints such as AUC and Cmax can often be (more) safely assumed to follow a normal distribution, we almost always carry out the bioequivalence statistical analysis on the log transformed data. Yet when present the results, and summarize the data, we want to speak about the original data.

For the original data, we then use the geometric mean as the description. This means, that we take an exponent of the quantity of the original data. In certain cases this quantity is a direct arithmetic mean, but in other cases, it can be a statistical point estimate, a confidence interval boundary or similar (more on this in the sequel).

To see this let us create the same NCA based dataset that we analyzed in unit 1 with a parallel study, this time focusing only on Cmax.
:::

```{julia}
#| output: false
df = dataset("nca/crossover_study.csv")

@subset!(df, :period .== 1)
@rtransform!(df, :sequence = :sequence == "RT" ? "R" : "T")

pop = read_nca(df; observations = :concentration, group = [:sequence, :period])
nca_result = run_nca(pop; parameters = [:cmax])
be_analysis_data = nca_result.reportdf;
```

::: {.content-hidden when-meta="profile-type.no-intermediate-text"}
Here is the dataset.
:::

```{julia}
simple_table(be_analysis_data)
```

::: {.content-hidden when-meta="profile-type.no-intermediate-text"}
Now let us repeat the application of `pumas_be` (this same output appeared in unit 1):
:::

```{julia}
pumas_be(be_analysis_data, endpoint = :cmax)
```

::: {.content-hidden when-meta="profile-type.no-intermediate-text"}
As you can see, we report a `Geometric Naive Mean` of `7735` and and `7610` for the reference and test respectively (the phrase `Naive` simply says it is based on the data and not on some model estimate, as may occur with more complex models as discussed in the sequel).

In this case of this simple model, let us see that these values are in agreement with the arithmetic means of the log transformed data. What `pumas_be` does is first log transform the data and computes the arithmetic means. These are the values (not visible as output):
:::

```{julia}
arithmetic_mean_log_R = mean(log.(be_analysis_data.cmax[1:10]))
arithmetic_mean_log_T = mean(log.(be_analysis_data.cmax[11:20]))

@show arithmetic_mean_log_R
@show arithmetic_mean_log_T;
```

::: {.content-hidden when-meta="profile-type.no-intermediate-text"}
Then when presenting the summary, the model exponentiates the arithmetic means, and as described above presents them as geometric means.
:::

```{julia}
@show exp(arithmetic_mean_log_R)
@show exp(arithmetic_mean_log_T);
```

::: {.content-hidden when-meta="profile-type.no-intermediate-text"}
This is a common bioequivalence pattern, where we work with the log transformed data internally for statistical computations, and then present results in the original scale.
:::

## The Geometric Mean Ratio (GMR)

::: {.content-hidden when-meta="profile-type.no-intermediate-text"}
At the end of the day, our goal with bioequivalence is to compare the test and the reference. In the common approach we wish to compare means, and this is called **average bioequivalence**. Now comparing of quantities may be done in two common ways: One way is taking the difference and seeing if it is greater or less than $0$. Another way is taking the ratio and seeing if it is greater or less than $1$.

For geometric means, the ratio is natural, and the ratio of the two geometric means is called the `Geometric Mean T/R Ratio` (or GMR - **not to be confused with the "Geometric mean of R"**). This value in percent is ideally as close to 100% as possible since 100% means that both geometric means are the same.

But observe also,
:::

$$
\text{GMR} =
\frac{\text{geomean}(x_1^T,\ldots,x_{n_T}^T)}{\text{geomean}(x_1^R,\ldots,x^R_{n_R})} =
e^{\text{mean}(\log(x_1^T),\ldots,\log(x_{n_T}^T)) ~-~ \text{mean}(\log(x_1^R),\ldots,\log(x^R_{n_R}))}.
$$

So,

$$
\log(\text{GMR}) = \text{mean}(\log(x_1^T),\ldots,\log(x_{n_T}^T)) ~-~ \text{mean}(\log(x_1^R),\ldots,\log(x^R_{n_R})).
$$

Here $x_i^T$ for $i=1,\ldots,n_T$ are the test observations and similarly $x_i^R$ for $i=1,\ldots,n_R$ are the reference observations.

::: {.content-hidden when-meta="profile-type.no-intermediate-text"}
Now indeed if we consider the inner workings of the model, we can inspect that statistical model via,
:::

```{julia}
be_out = pumas_be(be_analysis_data, endpoint = :cmax)
be_out.model[1][1]
# This use of [1][1] is not done typically,
# but is rather a way of seeing the model output in this case.
```

::: {.content-hidden when-meta="profile-type.no-intermediate-text"}
All of the details of this model are not typically important and this is why `pumas_be` only presents a neater table output as a summary.
:::

Importantly, you can see a `point estimate` in the model with the line `point estimate: -0.0161906`.

Let us also compute this same value manually,

```{julia}
round(arithmetic_mean_log_T - arithmetic_mean_log_R, digits = 7)
```

And here it is back in the original (geometric mean ratio) scale, presented as a rounded percentage, just like in the output of `pumas_be`:

```{julia}
round(100exp(arithmetic_mean_log_T - arithmetic_mean_log_R), digits = 2)
```

::: {.content-hidden when-meta="profile-type.no-intermediate-text"}
These same transformations also take place for quantities such as confidence interval bounds. We focus more deeply on confidence intervals in the sequel, yet in the context of the log transformations, still consider them here.
:::

Let us inspect an output of the `pumas_be` computation using the `result` field:

```{julia}
simple_table(be_out.result)
```

In particular, you may see that `lnLB = -0.203` and `lnUB = 0.17`.

::: {.content-hidden when-meta="profile-type.no-intermediate-text"}
These are bounds of a 90% confidence interval (as a side note observe that the 90% interval is narrower than the 95% confidence interval appearing in the model output above). As you can see, we can retrieve the 90% confidence interval on the original scale by exponentiation:
:::

```{julia}
println(round.(100exp.([only(be_out.result.lnLB), only(be_out.result.lnUB)]), digits = 2))
```

::: {.content-hidden when-meta="profile-type.no-intermediate-text"}
Compare this result with the `pumas_be` output above.
:::

## Bounds on the GMR

::: {.content-hidden when-meta="profile-type.no-intermediate-text"}
Much of the bioequivalence analysis that we do is **average bioequivalence**, where we want the geometric mean of `R` and the geometric mean of `T` to be close. For this we need to define what we mean by "close".
:::

Closeness is sometimes characterized by a percentage such as $20\%$ or $10\%$, where we want that the difference will **"not be more than 20%"**, or **"not be more than 10%"**, etc.

::: {.content-hidden when-meta="profile-type.no-intermediate-text"}
Let's consider the most common value, $20\%$, and see how it can be used with the geometric mean ratio (GMR). For this let's call $\text{GM}_T$ the geometric mean of the test formulation and $\text{GM}_R$ the geometric mean of the reference formulation. Hence,
:::

$$
\text{GMR} =
\frac{\text{GM}_T}{\text{GM}_R}.
$$

Now one natural way (**not used today**) to quantify a maximal difference of $20\%$ would be to require that,

$$
0.8\text{GM}_R \le \text{GM}_T \le 1.2\text{GM}_R,
$$

and by dividing the inequality by $\text{GM}_R$ we see that this is equivalent to,

$$
0.8 \le \text{GMR} \le 1.2.
$$

While this may seem natural, observe what happens if we take logarithms as required for statistical analysis. In particular $\log(0.8) \approx -0.223144$ and $\log(1.2) \approx 0.182322$. Hence requiring this from the GMR is equivalent to,

$$
-0.223144 \le \text{mean}(\log(x_1^T),\ldots,\log(x_{n_T}^T)) ~-~ \text{mean}(\log(x_1^R),\ldots,\log(x^R_{n_R})) \le 0.182322.
$$

As you can see we get a criterion that is **not** symmetric on the log scale.

::: {.content-hidden when-meta="profile-type.no-intermediate-text"}
The center of the interval $[-0.223144,  0.182322]$ is at about $-0.020411$ and one would expect this center to be the most probable case (in fact, once you visit the concept of power you may see it is the point with maximal power). Taking the exponent of this center we get $0.979796$, and this means that with such bounds, the "most powerful center" is at about $97.98\%$ and not at $1$!
:::

Due to this asymmetry on the log-scale, while early bioequivalence specifications used this approach, today the approach that is recommended is a bound of the form,

$$
0.8 \le \text{GMR} \le 1.25.
$$

Observe that $1/0.8 = 1.25$.

While these bounds are not symmetric on the non-log-transformed scale, they still can be interpreted as a 20% bound since $0.8 \le \text{GMR} \le 1.25$ is equivalent to,
$$
0.8\text{GM}_R \le \text{GM}_T
\qquad
\text{and}
\qquad
0.8\text{GM}_T \le \text{GM}_R.
$$

Importantly, they imply that
$$
-0.223144 \le \text{mean}(\log(x_1^T),\ldots,\log(x_{n_T}^T)) ~-~ \text{mean}(\log(x_1^R),\ldots,\log(x^R_{n_R})) \le 0.223144,
$$
which is centered at $0$. This puts $\text{GMR} = 1$, as the central point from a statistical perspective.

Finally note that if we want the more stringent $10\%$ difference and not the $20\%$ difference the bounds are
$$
0.9 \le \text{GMR} \le 1.1111,
$$

where we observe that $1/0.9 = 1.1\bar{1}$. On the log scale these bounds translate to

$$
-0.105361 \le \log(\text{GMR}) \le 0.105361.
$$

::: {.content-hidden when-meta="profile-type.no-intermediate-text"}
All of these above values are common in bioequivalence analysis. We finally mention that in the sequel, where we expand beyond average bioequivalence to reference scaling methods, we'll see that the concepts of bounds can also take variability of the product into account.
:::

## Coefficient of variation

::: {.content-hidden when-meta="profile-type.no-intermediate-text"}
The variability of a distribution is very naturally described by the **variance** or alternatively, by its square root, the **standard deviation**. The standard deviation is useful because it is in the same units of the quantity measured, and in case of a normal distribution, we know to expect almost all observations to be within $\pm 3$ standard deviations
:::

::: {.content-hidden when-meta="profile-type.no-intermediate-text"}
However, it is sometimes convenient to use a unit free, scaled version of the standard deviation which is called the **coefficient of variation**. In particular if the mean (this is the mathematical mean or expectation) of a random quantity $X$, is ${\mathbb E}[X]$ and the variance is $\text{Var}(X)$, then the coefficient of variation is,
:::

$$
\text{CV}(X) = \frac{\sqrt{\text{Var}(X)}}{{\mathbb E}[X]} = \frac{\text{standard deviation}}{\text{mean}}.
$$

::: {.content-hidden when-meta="profile-type.no-intermediate-text"}
The usefulness in using the coefficient of variation instead of the standard deviation is that it is comparable across scales. Pharmacology experts then get a sense of what are high CV values and what are low ones, independently of the units used.
:::

For example, as a general rule, in Bioequivalence analysis we can expect CV values to range from 0.05 to 1.5, with any values greater than 0.3 typically interpreted as exhibiting a high level of variability.

::: {.content-hidden when-meta="profile-type.no-intermediate-text"}
Now if we return to the log normal distribution and its formulas for the mean and variance we may compute the CV in terms of the parameters $\mu$ and $\sigma$. In fact, in this case the CV is sometimes called the **geometric CV**, yet we just call it "CV" in most cases.
:::

::: {.content-hidden when-meta="profile-type.no-intermediate-text"}
In particular for the log normal distribution using the formulas for the mean and variance from above,
:::

$$
\text{CV}(X) = \frac{\sqrt{\text{Var}(X)}}{{\mathbb E}[X]}= \frac{\sqrt{(e^{\sigma^2}-1)e^{2\mu + \sigma^2}}}{e^{\mu + \frac{\sigma^2}{2}}} = \sqrt{e^{\sigma^2}-1}.
$$

::: {.content-hidden when-meta="profile-type.no-intermediate-text"}
Interestingly we see that the CV does not depend on $\mu$ but only on $\sigma$. This is a desirable attribute. We can then easily transform from $\sigma$ to $\text{CV}$:
:::

$$
\text{CV} =\sqrt{e^{\sigma^2}-1}.
$$

::: {.content-hidden when-meta="profile-type.no-intermediate-text"}
Going the other way around from $\text{CV}$ to $\sigma$, we can solve the above equation to obtain.
:::

$$
\sqrt{\log(\text{CV}^2 + 1)} = \sigma.
$$

::: {.content-hidden when-meta="profile-type.no-intermediate-text"}
As a convenience, these transformations are also available with the `Bioequivalence` package via the `geocv2sigma` and `sigma2geocv` functions. For example, in the output for the parallel design, the output in the last row is ``Variability	CV (%) | σ̂	23.73 | 0.234`. This means that the CV is estimated at `23.73`percent and has a corresponding estimate for $\sigma$ which is`0.234`. You can see that these values match as follows:
:::

```{julia}
round(geocv2sigma(23.73 / 100), digits = 3)
```

```{julia}
100round(sigma2geocv(0.234), digits = 4)
```

Note very minor differences due to variations of rounding.

## Maximum likelihood estimation for the log-normal distribution

::: {.content-hidden when-meta="profile-type.no-intermediate-text"}
Statistical inference, both in bioequivalence and more generally, often requires estimating the parameters of the underlying distribution. Here let's discuss estimation for the log-normal distribution. A fundamental estimation method is **maximum likelihood estimation** (MLE), which chooses the parameter values that make the observed data most probable under the assumed model.

Recall, if our endpoint data for a treatment (`T` or `R`), denoted $x_1, \ldots, x_n$, are assumed to be independent and log-normally distributed, then their logarithms $\log(x_i)$ are independent and normally distributed with mean $\mu$ and variance $\sigma^2$.

In MLE, we begin by writing down the **likelihood function**, which expresses the probability (density) of observing the actual data as a function of the distribution’s parameters. For the log-normal case, this means finding the likelihood of observing our endpoint data given unknown parameters $\mu$ and $\sigma^2$, which fully determine the shape and location of the distribution. The goal is to choose the values of $\mu$ and $\sigma^2$ that maximize this likelihood, making the observed data as "likely" as possible under the model. This approach leads directly to concrete formulas for estimating these parameters.
:::

The likelihood function for the observed data, as a function of $\mu$ and $\sigma^2$, is given by:
$$
L(\mu, \sigma^2\,|\,x_1, \ldots, x_n) = \prod_{i=1}^n \frac{1}{x_i \sqrt{2\pi\sigma^2}} \exp\left(-\frac{(\log x_i - \mu)^2}{2\sigma^2}\right).
$$

::: {.content-hidden when-meta="profile-type.no-intermediate-text"}
With maximum likelihood we now have an **optimization problem** of finding $\mu$ and $\sigma^2$ that maximize this likelihood function. In certain cases, this optimization is solved numerically, but in other cases such as here we can find a formula (closed form) solution.

In any case, it is often common to consider the **log-likelihood** which is the logarithm of the likelihood function. We can optimize the log-likelihood instead of the likelihood and the optimizers will be the same. It is denoted via, $\ell(\mu, \sigma^2)$, and is,
:::

$$
\ell(\mu, \sigma^2) = -n \log(\sigma) - n \log(\sqrt{2\pi}) - \sum_{i=1}^n \log x_i - \frac{1}{2\sigma^2} \sum_{i=1}^n (\log x_i - \mu)^2.
$$

Now with some calculus, it is possible to take derivatives of the log-likelihood, equate to zero, and solve for $\mu$ and $\sigma^2$. It is also possible to consider the second order considerations to see it is a maximum (we omit these details). Importantly, after such a calculation, the MLE for the parameters are:

::: {.content-hidden when-meta="profile-type.no-intermediate-text"}
The parameter $\mu$ is estimated as the arithmetic mean of the log-transformed data:
:::

$$
\hat\mu = \frac{1}{n} \sum_{i=1}^n \log x_i
$$

::: {.content-hidden when-meta="profile-type.no-intermediate-text"}
The parameter $\sigma^2$ is estimated as the sample variance of the log-transformed data:
:::

$$
\hat\sigma^2 = \frac{1}{n} \sum_{i=1}^n (\log x_i - \hat\mu)^2
$$

::: {.content-hidden when-meta="profile-type.no-intermediate-text"}
(Note: For unbiasedness, sometimes $1/(n-1)$ is used in place of $1/n$.)

These estimators are direct consequences of the normal likelihood for the transformed data.
:::

### Geometric mean and the MLE

::: {.content-hidden when-meta="profile-type.no-intermediate-text"}
Recall from the previous section that the geometric mean of $x_1,\ldots,x_n$ is related to $\hat\mu$:
:::

$$
\text{Geometric Mean} = \left( \prod_{i=1}^n x_i \right)^{1/n} = \exp\left( \frac{1}{n} \sum_{i=1}^n \log x_i \right) = \exp(\hat\mu).
$$

Thus, the **MLE for the location parameter in the log-normal model corresponds to the logarithm of the geometric mean** of the data.

::: {.content-hidden when-meta="profile-type.no-intermediate-text"}
This tight relationship further underscores why geometric means—and not arithmetic means—are the sensible notion of "average" for log-normally distributed quantities such as pharmacokinetic endpoints.
:::

## Conclusion

This unit of the bioequivalence analysis course introduces foundational statistical concepts essential for bioequivalence studies, specifically focusing on the use of log-transformed data, the log-normal distribution, geometric means, geometric mean ratios (GMR), and the coefficient of variation (CV).

::: {.content-hidden when-meta="profile-type.no-intermediate-text"}
It explains why pharmacokinetic endpoints like Cmax and AUC are better modeled by log-normal rather than normal distributions and demonstrates how log transformation allows for more appropriate statistical analysis. The unit details how geometric means (derived from log-transformed arithmetic means) provide meaningful central tendency for multiplicative data, and clarifies the interpretation and calculation of the GMR and its confidence intervals—central to average bioequivalence assessment. In addition, it covers the significance of bounds for GMR (such as the equivalence range 0.8–1.25), and introduces CV as a unit-free measure of variability, highlighting its role and calculation in the context of log-normal data. Throughout, practical examples using Pumas and relevant packages reinforce the theoretical discussion by analyzing both real and synthetic datasets.
:::

## Unit exercises

 1. **Log-Transformation and Normality**
    
    Using the provided pharmacokinetic endpoint dataset (e.g., `be_analysis_data` for Cmax), log-transform the data for both Test and Reference treatments. Create histogram or kernel density plots for:
    
    (a) the original data, and
    (b) the log-transformed data.
    
    Visually assess and comment on the distributional shape before and after transformation. Is the log-transformed data approximately normal? Note that in standard bioequivalence practice, such "checking" of normality is not carried out.

 2. **Geometric Mean Calculation**
    
    Given the following data for a Test group:
    $x = [6200, 7800, 9900, 8500, 7400, 6700, 8000, 9200, 8800, 7600]$
    (a) Compute the arithmetic mean of the original data.
    (b) Compute the geometric mean (you may use Julia’s `geomean()` or compute it manually).
    (c) Explain in which context each mean is most meaningful for pharmacokinetic endpoints like Cmax.
 3. **MLE for Log-Normal Parameters**
    
    Implement the MLE estimators for $\mu$ and $\sigma^2$ for the log-normal distribution using the data from Exercise 2. Provide the code and final estimated values of $\hat{\mu}$ and $\hat{\sigma}^2$.
 4. **GMR and Log Difference**
    
    Let the sample means of log-transformed Cmax values for Test and Reference be $\bar{x}_T = 8.948$ and $\bar{x}_R = 8.970$, respectively.
    (a) Compute the geometric mean ratio (GMR) of Test to Reference.
    (b) Express the GMR as a percentage.
    (c) What does this value mean in a bioequivalence analysis?
 5. **CV and $\sigma$, Conversion Practice**
    
    For a log-normal endpoint, suppose the estimated $\sigma$ from the log-transformed data is $0.21$.
    (a) Calculate the coefficient of variation (CV) as a percentage.
    (b) If the estimated CV was 36%, what would be the corresponding $\sigma$?
    (You may use the formulas or Julia functions provided in the unit.)
 6. **Equivalence Bounds and Interpretation**
    
    Explain why the commonly used bioequivalence bounds for GMR are $[0.8, 1.25]$, rather than $[0.8, 1.2]$ as might seem natural. Show by computation (in Julia or by calculator) that $\log(0.8)$ and $\log(1.25)$ are equidistant from zero, and interpret the practical meaning of this in the context of statistical analysis on the log scale.

## References {.unnumbered}
