A Course on Bioequivalence: Unit 10 - Details of Linear and Mixed Effects Models

This is Unit 10 of a complete bioequivalence analysis course using Pumas. There are 15 units in the course in total.

1 Unit overview

This unit fills in some missing details of fixed effects linear models and mixed effects (random effects) linear models. In particular we continue the discussion from Unit 5 and present further details of models used for bioequivalence.

As you can recall from Unit 5:

• With the fixed effects models, pumas_be uses Julia’s GLM package. This yields behavior compatible with SAS’s PROC GLM.
• With the mixed effects models, pumas_be uses Julia’s MixedModels package. This yields behavior compatible with SAS’s PROC MIXED.

As a general rule we use the model based on the type of design:

• If the design is a non-replicate cross over design, a fixed effects model is used.
• If the design is a replicate design (partial replicate or fully replicate), a mixed effects model is used.

Further, in replicate designs, an additional auxiliary estimator is used to estimate the within subject variability of either the reference product, the test product, or both products.

Some FDA examples present the computation of this auxiliary estimator using SAS’s PROC MIXED, but it can also be done with a fixed effects model as carried out in pumas_be. In any case, we present the details of this auxiliary estimator, which with Pumas can also be called directly via the within_subject_variability function.

In addition to these topics, this unit presents the fixed effects models used for cases of multiple formulations with multiple test products. This continues the discussion of multiple formulations in Unit 4 where we called such cases 1RnT cases (with n=2 or n=3).

Hence in summary this unit presents the following topics:

• A review and further investigation of the fixed effects model formulation used for standard \(2 \times 2\) crossover designs as well as crossover designs with multiple test products (1RnT).
• A dive into the mixed effect model formulations used for partial and fully replicate designs. This includes a brief discussion of REML (residual maximum likelihood estimation) and a brief mention of the Satterthwaite computation for such models.
• A discussion of the model used to estimate within subject variability (used with the within_subject_variability Pumas function, and also used implicitly with pumas_be).

We use the following packages.

using Pumas # available with Pumas products - also reexports Distributions
using NCA # available with Pumas products
using Bioequivalence  # available with Pumas products
using BioequivalencePower  # available with Pumas products
using PharmaDatasets # available with Pumas products
using GLM
using DataFramesMeta
using SummaryTables
using DataFrames
using DataFramesMeta
using AlgebraOfGraphics
using CairoMakie
using Random
using LinearAlgebra
using MixedModels

We have used all of these packages in previous units.

2 Fixed effects models

Mathematically (and loosely) a fixed effect model has the general form:

\[ \underbrace{y}_{\text{log endpoint}} = \underbrace{\beta^{\top} x}_{\text{linear combination of fixed effects}} + \qquad \underbrace{\varepsilon}_{\text{noise}} \]

Here \(\beta\) is a vector of coefficients that is estimated and \(x\) is composed of the predictors. Hence the response \(y\) (typically a log transformed endpoint) is a linear function of the predictors \(x\). The noise term \(\varepsilon\) is considered to be normally distributed with a \(0\) mean and a standard deviation \(\sigma\), which is also an unknown model parameter.

See Unit 5 for an introductory discussion of the linear model. You can also see more details in Nazarathy and Klok (2021), Liquet, Moka, and Nazarathy (2024), and references within.

A review of the basic bioequivalence linear model

Let’s return to a basic \(2 \times 2\) crossover design example:

data = dataset(joinpath("bioequivalence", "RT_TR", "PJ2017_3_12.csv"))
be_result = pumas_be(data)

Observation Counts
Sequence ╲ Period	1	2
RT	23	23
TR	23	23

Paradigm: Non replicated crossover design
Model: Linear model
Criteria: Standard ABE
Endpoint: AUC
Formulations: Reference(R), Test(T)
		Results(AUC)	Assessment	Criteria
R	Geometric Marginal Mean	78.95
	Geometric Naive Mean	83.08
T	Geometric Marginal Mean	86.99
	Geometric Naive Mean	86.35
	Geometric Mean T/R Ratio (%)	110.2
	Degrees of Freedom	43
	90% Confidence Interval (%)	[94.08, 129]	Fail	CI ⊆ [80, 125]
Variability	CV (%) | σ̂	46.89 | 0.4458
ANOVA	Formulation (p-value)	0.308
	Sequence (p-value)	0.1742
	Period (p-value)	0.5947

As you may recall from Unit 5, this is the formula used:

be_result.model |> formula |> println

endpoint ~ 1 + sequence + formulation + period + sequence & id_within_sequence

All of the effects 1 (the intercept), sequence (the sequence effect), formulation (the formulation effect), period (the period effect), and sequence & id_within_sequence (the subject effect), are part of the linear combination \(\beta^{\top} x\).

The bulk of the model is thus captured by the coefficients \(\beta\) which we can retrieve via the coef function:

round.(coef(be_result.model), digits = 3) |> println

[4.459, -0.13, 0.097, -0.05, -0.135, -1.743, 0.344, -1.123, 1.154, -1.032  …  -1.63, -0.345, -1.132, -0.934, 1.498, 0.672, -0.14, 2.598, 0.013, 0.0]

Note that there are 50 coefficients here, most of which are for the `sequence & id_within_sequence subject effects.

Fixed effects models for more than two formulations

As we saw at the end of Unit 4, we can also have a single model, even if there are more than two formulations. In particular we support the case of three formulations (one reference formulation and two test formulations), and four formulations (one reference formulation three test formulations). These cases are most popular in the 1RnT setting where the additional formulations are considered as test products (in the case of multiple reference products the recommendation is typically to consider different models).

When we run bioequivalence analysis for such cases, even though a separate hypothesis test is considered for each formulation, Pumas still uses a single linear model for all formulations. This allows to use all of the available data to estimate the variability and to consider the test formulation. This is most suitable for what we call the `

As an example let us focus on the case with three formulations with this dataset:

data_3_formulations = dataset("bioequivalence/RST_RTS_SRT_STR_TRS_TSR/PJ2017_4_5")
detect_design(data_3_formulations)

(design = "RST|RTS|SRT|STR|TRS|TSR",
 num_sequences = 6,
 num_periods = 3,
 num_formulations = 3,
 reference_formulation = 'R',
 test_formulations = ['S', 'T'],
 replicated = non_replicated,
 supports_rsabe = false,
 crossover = true,)

As you can see, there are 6 sequences in this design, comprised of all possible permutations of the reference R, and the test formulations S and T. To see the size of the dataset use:

size(data_3_formulations)

(186, 5)

or alternatively to see the number of rows use:

nrow(data_3_formulations)

186

Here are the first few rows of the dataset:

simple_table(first(data_3_formulations, 8))

id	sequence	period	AUC	Cmax
1	SRT	1	7260	1633
1	SRT	2	6463	1366
1	SRT	3	8759	2141
2	RTS	1	3457	776
2	RTS	2	6556	2387
2	RTS	3	4081	1355
4	TSR	1	4006	1326
4	TSR	2	4879	1028

And here are the last few rows of the dataset:

simple_table(last(data_3_formulations, 8))

id	sequence	period	AUC	Cmax
62	SRT	2	7069	1995
62	SRT	3	6530	1236
63	TRS	1	2204	495
63	TRS	2	2927	770
63	TRS	3	missing	missing
67	RST	1	4045	1025
67	RST	2	7865	2668
67	RST	3	missing	missing

As you can see, there is some missing endpoint data in the dataset. When carrying out analysis for a specific endpoint (either AUC or Cmax in this case), Pumas will throw away lines with missing values for the endpoint.

Let’s run pumas_be on both endpoints. You can see that the observation counts differ per endpoint (click the tabs to change between AUC and Cmax):

AUC

auc_result = pumas_be(data_3_formulations, endpoint = :AUC)

Observation Counts
Sequence ╲ Period	1	2	3
RST	8	8	8
RTS	11	11	11
SRT	11	11	11
STR	10	10	10
TRS	10	10	10
TSR	10	10	10

Paradigm: 3 formulations
Model: Linear model
Criteria: Standard ABE
Endpoint: AUC
Formulations: Reference(R), Test(S and T)
		Results(AUC)	Assessment	Criteria
R	Geometric Marginal Mean	3808
	Geometric Naive Mean	3824
S	Geometric Marginal Mean	5346
	Geometric Naive Mean	5412
	Geometric Mean S/R Ratio (%)	140.4
	Degrees of Freedom	114
	90% Confidence Interval (%)	[131.7, 149.7]	Fail	CI ⊆ [80, 125]
T	Geometric Marginal Mean	4426
	Geometric Naive Mean	4426
	Geometric Mean T/R Ratio (%)	116.2
	Degrees of Freedom	114
	90% Confidence Interval (%)	[109, 123.9]	Pass	CI ⊆ [80, 125]
Variability	CV (%) | σ̂	21.22 | 0.2099

Cmax

cmax_result = pumas_be(data_3_formulations, endpoint = :Cmax)

Observation Counts
Sequence ╲ Period	1	2	3
RST	9	9	8
RTS	11	11	11
SRT	11	11	11
STR	10	10	10
TRS	11	11	10
TSR	10	10	10

Paradigm: 3 formulations
Model: Linear model
Criteria: Standard ABE
Endpoint: Cmax
Formulations: Reference(R), Test(S and T)
		Results(Cmax)	Assessment	Criteria
R	Geometric Marginal Mean	832.4
	Geometric Naive Mean	837.5
S	Geometric Marginal Mean	1327
	Geometric Naive Mean	1340
	Geometric Mean S/R Ratio (%)	159.4
	Degrees of Freedom	118
	90% Confidence Interval (%)	[146.1, 173.9]	Fail	CI ⊆ [80, 125]
T	Geometric Marginal Mean	1082
	Geometric Naive Mean	1078
	Geometric Mean T/R Ratio (%)	129.9
	Degrees of Freedom	118
	90% Confidence Interval (%)	[119.1, 141.8]	Fail	CI ⊆ [80, 125]
Variability	CV (%) | σ̂	29.69 | 0.2907

Let’s focus on the AUC case. This is the model formula, as you can see it is exactly the same formula like the \(2 \times 2\) case:

auc_result.model |> formula |> println

endpoint ~ 1 + sequence + formulation + period + sequence & id_within_sequence

The difference is that there are more levels for the various factors.

Note that for brevity here we only display the first 15 rows of the coefficient table:

simple_table(first(coeftable(auc_result.model), 15), halign = :left)

Name	Coef.	Std. Error	t	Pr(>|t|)	Lower 95%	Upper 95%
(Intercept)	8.18	0.0563	145	3.02e-131	8.07	8.29
sequence: RTS	0.0481	0.0585	0.823	0.412	-0.0678	0.164
sequence: SRT	0.172	0.0585	2.93	0.00407	0.0556	0.287
sequence: STR	-0.1	0.0596	-1.68	0.0948	-0.219	0.0177
sequence: TRS	0.187	0.061	3.06	0.00273	0.0661	0.308
sequence: TSR	0.0529	0.0596	0.888	0.377	-0.0652	0.171
formulation: S	0.339	0.0386	8.79	1.8e-14	0.263	0.416
formulation: T	0.15	0.0386	3.9	0.000163	0.0741	0.227
period: 2	0.053	0.0384	1.38	0.169	-0.0229	0.129
period: 3	-0.0294	0.0388	-0.756	0.451	-0.106	0.0476
sequence: RST & id_within_sequence: 2	-0.405	0.116	-3.48	0.000702	-0.635	-0.175
sequence: RTS & id_within_sequence: 2	0.0825	0.116	0.714	0.477	-0.146	0.311
sequence: SRT & id_within_sequence: 2	-0.521	0.116	-4.51	1.59e-5	-0.75	-0.292
sequence: STR & id_within_sequence: 2	-0.125	0.115	-1.09	0.279	-0.353	0.103
sequence: TRS & id_within_sequence: 2	0.378	0.117	3.24	0.00159	0.147	0.61

You can see:

• The sequence factor (with 6 levels) has 5 coefficients.
• The period factor (with 3 levels) has 2 coefficients.
• Most importantly, the formulation factor (with 3 levels) has 2 coefficients.

Let’s extract the coefficients associated with formulation:

β_S, β_T = coef(auc_result.model)[7:8]
(β_S, β_T)

(0.3392505331133662, 0.15049814185173857)

We can then obtain the geometric mean ratio estimates (GMR) results via these:

GMR_S = round(100exp(β_S), digits = 2)
GMR_T = round(100exp(β_T), digits = 2)
(GMR_S, GMR_T)

(140.39, 116.24)

You can see (that up to rounding variations) these values agree with the pumas_be output above.

Degrees of freedom and variability estimation

The degrees of freedom (DOF) of the model’s residuals represent the number of independent pieces of information available to estimate the model’s error variance (or variability).

Intuitively, you can think of it this way: you start with a certain number of observations (n), and each time you use the data to estimate a parameter (like an intercept or a coefficient for sequence or period), you “use up” one degree of freedom. What remains are the residual degrees of freedom, used to estimate the variability. This is why, in the simplest case, they are calculated as the total number of observations minus the number of estimated parameters.

This value is critical for obtaining an unbiased estimate of the residual standard deviation and for determining the specific t-distribution used in confidence interval calculations.

In the simple case with n observations and p estimated coefficients, are obtained via n - p. For example in our case,

n = nrow(auc_result.data)
p = length(coef(auc_result.model))
n, p

(180, 70)

And hence we may expect n-p degrees of freedom:

dof = n - p

110

However, more specifically, the DOF are determined by properties of the design matrix of the data. If really needed, you can access this design matrix, X, via:

X = auc_result.model.mm.m
size(X)

(180, 70)

This matrix is internally used for least squares estimation of the model.

Now it sometimes occurs, due to missing observations that that the columns of the design matrix are not linearly independent. From a mathematical (linear algebra) point of view, this means that the column rank of the design matrix is less thank the number of columns. If really needed we can investigate it via the rank function from the LinearAlgebra package:

rank(X)

66

With this, the correct way to compute the degrees of freedom is via:

dof = n - rank(X)

114

Indeed this agrees with the output of dof_residual from the GLM package:

dof = dof_residual(auc_result.model)

114.0

Now with this, we can obtain \(\hat{\sigma}\), the estimate of the standard deviation of \(\varepsilon\) via:

σ̂ = round(sqrt(deviance(auc_result.model) / dof), digits = 4)

0.2099

Note that we can also use the dispersion function, but in this case we need to use model.model:

round(dispersion(auc_result.model.model), digits = 4)

0.2099

In any case, we can consider the value as a CV:

round(100sigma2geocv(σ̂), digits = 2)

21.22

As you can see, this value agrees with the Variability estimate in the above pumas_be output.

3 Mixed (random) effects models

Mathematically (and loosely) a mixed (random) effects model has the general form:

\[ \underbrace{y}_{\text{log endpoint}} = \underbrace{\beta^{\top} x}_{\text{linear combination of fixed effects}} + \qquad \underbrace{\zeta}_{\text{random effects}} + \qquad \underbrace{\varepsilon}_{\text{noise}} \]

Here the additional term beyond fixed effects models is denoted as \(\zeta\). It is a a random noise term which can also depend on the predictors \(x\) and on other factors.

The most basic use for it in the case of clinical trials (and bioequivalence trials) is to capture subject effects.

When considering any mixed effects model, our goal is not only to estimate \(\beta\) and the variability of \(\varepsilon\), \(\sigma\), but to also estimate the variance-covariance structure of \(\zeta\).

Two variations of mixed effects models for replicated designs

In Unit 5 we saw the use of the MixedModels package directly for a \(2 \times 2\) crossover design. This agrees with cases where one wants to use the FDA’s SAS PROC MIXED suggestion.

When we move onto replicate (partial or complete) designs, pumas_be also uses mixed models, directly under the hood. To see this consider this fully replicate dataset:

data = dataset(joinpath("bioequivalence", "RTTR_TRRT", "PJ2017_4_3"));

Let us run pumas_be on this dataset, where we also consider two options for the homogeneity argument, false for an unequal variance assumption, and true for an equal variance assumption (the default is false).

Click on this tabset to alternate between the outputs. You will see slight differences between the Degrees of Freedom, and the associated quantities such as the GMR estimate and the \(90\%\) confidence interval.

Equal variance assumption

With homogeneity = true we enforce an equal variance assumption.

result_yes_homogenous = pumas_be(data, endpoint = :Cmax, homogeneity = true)

Observation Counts
Sequence ╲ Period	1	2	3	4
RTTR	8	8	8	8
TRRT	9	9	9	8

Paradigm: Replicated crossover that supports reference scaling
Model: Mixed model (equal variance)
Criteria: Standard ABE
Endpoint: Cmax
Formulations: Reference(R), Test(T)
		Results(Cmax)	Assessment	Criteria
R	Geometric Marginal Mean	1150
	Geometric Naive Mean	1150
T	Geometric Marginal Mean	1049
	Geometric Naive Mean	1048
	Geometric Mean T/R Ratio (%)	91.21
	Degrees of Freedom	45.97
	90% Confidence Interval (%)	[83.21, 99.99]	Pass	CI ⊆ [80, 125]
Variability	CVᵣ (%) | σ̂ᵣ	21.17 | 0.2094
	CVₜ (%) | σ̂ₜ	26.87 | 0.264
	Variability Ratio (%)	126.2
ANOVA	Formulation (p-value)	0.09955
	Sequence (p-value)	0.5563
	Period (p-value)	0.4375

Unequal variance assumption

The default for linear mixed models is homogeneity = false. Hence we don’t need to use explicitly use this option.

result_non_homogenous = pumas_be(data, endpoint = :Cmax)

Observation Counts
Sequence ╲ Period	1	2	3	4
RTTR	8	8	8	8
TRRT	9	9	9	8

Paradigm: Replicated crossover that supports reference scaling
Model: Mixed model (unequal variance)
Criteria: Standard ABE
Endpoint: Cmax
Formulations: Reference(R), Test(T)
		Results(Cmax)	Assessment	Criteria
R	Geometric Marginal Mean	1150
	Geometric Naive Mean	1150
T	Geometric Marginal Mean	1049
	Geometric Naive Mean	1048
	Geometric Mean T/R Ratio (%)	91.23
	Degrees of Freedom	40.82
	90% Confidence Interval (%)	[83.18, 100.1]	Pass	CI ⊆ [80, 125]
Variability	CVᵣ (%) | σ̂ᵣ	21.17 | 0.2094
	CVₜ (%) | σ̂ₜ	26.87 | 0.264
	Variability Ratio (%)	126.2
ANOVA	Formulation (p-value)	0.1019
	Sequence (p-value)	0.5596
	Period (p-value)	0.4916

Let us see the formulas used across the models.

println("Equal variance:")
println(result_yes_homogenous.model.formula)
println("\nUnequal variance:")
println(result_non_homogenous.model.formula)

Equal variance:
endpoint ~ 1 + formulation + sequence + period + (1 | id)

Unequal variance:
endpoint ~ 1 + formulation + sequence + period + (formulation + 0 | id) + zerocorr(formulation + 0 | row)

In both cases the fixed effects are the intercept (1), and then the formulation, the sequence, and the period. These same fixed effects were in the non-random effects (non-mixed) models analyzed previously.

The difference is with the random effects. These are terms with a | character.

The equal variance model (homogeneity = true) has:

• The random effect (1 | id) term implies a random effect per id. This is similar to the model used in Unit 5. The | id syntax means “grouped by id”. Hence each id has its own instance of the random effect \(\zeta\), but since there are multiple occurrences of the endpoint for the subject id, they all have the same instance of subject effect.

The unequal variance model (homogeneity = false) has two variance components:

• The random effect (formulation + 0 | id) which implies that each id has a different instance of the random effect like in the case of equal model variance, but this time, the instances associated with formulations R and T are allowed to have different variability. I.e. there is one variability for R and one for T.
• The random effect zerocorr(formulation + 0 | row) is an unusual grouping in the world of mixed models, yet it is the suitable model in this case. The zerocorr() wrapper specifies that the random effects for this grouping factor have zero correlations between their elements—i.e., the random effects are uncorrelated. The formulation + 0 | row syntax means that for each row, we consider random effects based on the formulation. This is a way to allow the unexplained variability in the model to have different variability parameters for each formulation.

With these interpretations it is useful to consider the model field of the pumas_be output in both cases:

Model with equal variance assumption

result_yes_homogenous.model

	Est.	SE	z	p	σ_id
(Intercept)	6.9935	0.0747	93.56	<1e-99	0.0871
formulation: T	-0.0920	0.0547	-1.68	0.0928
sequence: TRRT	-0.0416	0.0692	-0.60	0.5471
period: 2	0.0832	0.0767	1.08	0.2783
period: 3	0.0942	0.0767	1.23	0.2194
period: 4	0.1221	0.0780	1.57	0.1172
Residual	0.2234

Model with unequal variance assumption

result_non_homogenous.model

	Est.	SE	z	p	σ_row	σ_id
(Intercept)	6.9988	0.0724	96.62	<1e-99
formulation: T	-0.0918	0.0549	-1.67	0.0942	0.1369	0.0771
sequence: TRRT	-0.0415	0.0696	-0.60	0.5507
period: 2	0.0757	0.0760	1.00	0.3194
period: 3	0.0908	0.0760	1.19	0.2321
period: 4	0.1115	0.0765	1.46	0.1449
formulation: R					0.0196	0.1004
Residual	0.2000

In the equal variance model we see the estimate σ_id (0.0871) for the random effects in addition to the residual variability estimate (0.2234).

In the unequal variance model we see that the estimate σ_id is split across formulation: T (0.0771) and formulation: R (0.1004). Further the estimate σ_row which goes on each row of the model is also split across the formulations with formulation: T (0.1369) and formulation: R (0.0196). Finally, there is also the residual variability estimate (0.2000).

This latter model can potentially capture more variability in the data and with this potentially obtain more suitable estimates for the formulation effect.

The Satterthwaite computation

You may notice the degrees of freedom appearing in the outputs of the above linear models. For the equal variance model assumption the Degrees of Freedom term is 45.97, and for the unequal variance model the Degrees of Freedom term is 40.82.

You may obtain these values (and any other value appearing in pumas_be output) via the result_table field. For example,

dof_yes_homogenous = result_yes_homogenous.result_table[:DF]

45.97

dof_non_homogenous = result_non_homogenous.result_table[:DF]

40.82

Now these degrees of freedom are computed via an approach called Satterthwaite’s method for t-statistics. A brief summary of the approach is Haubo and Christensen (2020), at this point, let us take these computations as given.

Importantly, we now calculate a confidence interval based on these degrees of freedom. We also use the stderror function from the MixedModels package. Here is the computation for each case. You can compare we obtain the same confidence intervals as those in pumas_be.

Model with equal variance assumption

β_T = coef(result_yes_homogenous.model)[2]

-0.09195391814186732

SE = stderror(result_yes_homogenous.model)[2]

0.05470239368161244

t_quantile = quantile(TDist(dof_yes_homogenous), 0.95)

1.6786829275296131

ci = [β_T - SE * t_quantile, β_T + SE * t_quantile]
round.(100exp.(ci), digits = 2) |> println

[83.21, 99.99]

Model with unequal variance assumption

β_T = coef(result_non_homogenous.model)[2]

-0.09181507690732098

SE = stderror(result_non_homogenous.model)[2]

0.05486184122988195

t_quantile = quantile(TDist(dof_non_homogenous), 0.95)

1.6830495451436747

ci = [β_T - SE * t_quantile, β_T + SE * t_quantile]
round.(100exp.(ci), digits = 2) |> println

[83.18, 100.05]

Restricted maximum likelihood (REML)

The parameter estimation in mixed models is typically initially positioned as maximum likelihood estimation (MLE). With an MLE approach, the estimation is based on a strong theoretical footing which also supports associated hypothesis tests. However, with small samples as often occurring in the context of maximum likelihood, and/or with complex variance structures of mixed models, an alternative popular approach is restricted maximum likelihood (REML). In a nutshell, such an approach presents better (less biased) estimates of variance parameters.

See Section 5.1 of Patterson and Jones (2017) for further discussion of REML in our context. Additional references are suggested in that reference.

From the perspective of pumas_be, the reml option is by default true. If wishing to carry out estimation without REML, but rather with a direct maximum likelihood approach, you may set reml = false. For example, here is the equal variance model from above, only this time with direct maximum likelihood estimation, as opposed to residual maximum likelihood.

REML estimation

be_result = pumas_be(data, endpoint = :Cmax, reml = true)

Observation Counts
Sequence ╲ Period	1	2	3	4
RTTR	8	8	8	8
TRRT	9	9	9	8

Paradigm: Replicated crossover that supports reference scaling
Model: Mixed model (unequal variance)
Criteria: Standard ABE
Endpoint: Cmax
Formulations: Reference(R), Test(T)
		Results(Cmax)	Assessment	Criteria
R	Geometric Marginal Mean	1150
	Geometric Naive Mean	1150
T	Geometric Marginal Mean	1049
	Geometric Naive Mean	1048
	Geometric Mean T/R Ratio (%)	91.23
	Degrees of Freedom	40.82
	90% Confidence Interval (%)	[83.18, 100.1]	Pass	CI ⊆ [80, 125]
Variability	CVᵣ (%) | σ̂ᵣ	21.17 | 0.2094
	CVₜ (%) | σ̂ₜ	26.87 | 0.264
	Variability Ratio (%)	126.2
ANOVA	Formulation (p-value)	0.1019
	Sequence (p-value)	0.5596
	Period (p-value)	0.4916

be_result.model |> println

Linear mixed model fit by REML
 endpoint ~ 1 + formulation + sequence + period + (formulation + 0 | id) + zerocorr(formulation + 0 | row)
 REML criterion at convergence: 13.441805854606276

Variance components:
             Column      Variance  Std.Dev.   Corr.
row      formulation: R  0.0003849 0.0196180
         formulation: T  0.0187330 0.1368686   .  
id       formulation: R  0.0100757 0.1003778
         formulation: T  0.0059463 0.0771126 +1.00
Residual                 0.0400044 0.2000111
 Number of obs: 67; levels of grouping factors: 67, 17

  Fixed-effects parameters:
───────────────────────────────────────────────────────
                     Coef.  Std. Error      z  Pr(>|z|)
───────────────────────────────────────────────────────
(Intercept)      6.99877     0.0724352  96.62    <1e-99
formulation: T  -0.0918151   0.0548618  -1.67    0.0942
sequence: TRRT  -0.0415016   0.0695515  -0.60    0.5507
period: 2        0.0756548   0.0759892   1.00    0.3194
period: 3        0.0908021   0.0759892   1.19    0.2321
period: 4        0.111534    0.0765      1.46    0.1449
───────────────────────────────────────────────────────

ML estimation

be_result = pumas_be(data, endpoint = :Cmax, reml = false)

Observation Counts
Sequence ╲ Period	1	2	3	4
RTTR	8	8	8	8
TRRT	9	9	9	8

Paradigm: Replicated crossover that supports reference scaling
Model: Mixed model (unequal variance)
Criteria: Standard ABE
Endpoint: Cmax
Formulations: Reference(R), Test(T)
		Results(Cmax)	Assessment	Criteria
R	Geometric Marginal Mean	1150
	Geometric Naive Mean	1150
T	Geometric Marginal Mean	1049
	Geometric Naive Mean	1048
	Geometric Mean T/R Ratio (%)	91.24
	Degrees of Freedom	45.06
	90% Confidence Interval (%)	[83.52, 99.67]	Pass	CI ⊆ [80, 125]
Variability	CVᵣ (%) | σ̂ᵣ	21.17 | 0.2094
	CVₜ (%) | σ̂ₜ	26.87 | 0.264
	Variability Ratio (%)	126.2
ANOVA	Formulation (p-value)	0.08819
	Sequence (p-value)	0.5352
	Period (p-value)	0.4537

be_result.model |> println

Linear mixed model fit by maximum likelihood
 endpoint ~ 1 + formulation + sequence + period + (formulation + 0 | id) + zerocorr(formulation + 0 | row)
   logLik   -2 logLik     AIC       AICc        BIC    
     4.9215    -9.8430    14.1570    19.9347    40.6133

Variance components:
             Column      Variance  Std.Dev.   Corr.
row      formulation: R  0.0001768 0.0132959
         formulation: T  0.0171759 0.1310570   .  
id       formulation: R  0.0083688 0.0914809
         formulation: T  0.0048674 0.0697665 +1.00
Residual                 0.0369367 0.1921891
 Number of obs: 67; levels of grouping factors: 67, 17

  Fixed-effects parameters:
────────────────────────────────────────────────────────
                     Coef.  Std. Error       z  Pr(>|z|)
────────────────────────────────────────────────────────
(Intercept)      6.9987      0.0686395  101.96    <1e-99
formulation: T  -0.0916953   0.0526136   -1.74    0.0814
sequence: TRRT  -0.0413819   0.0653844   -0.63    0.5268
period: 2        0.0755887   0.0728789    1.04    0.2997
period: 3        0.0907633   0.0728789    1.25    0.2130
period: 4        0.111669    0.0733674    1.52    0.1280
────────────────────────────────────────────────────────

4 Within Subject Variability estimation

While we may potentially estimate the within subject variability directly from the model, a common approach is to setup auxiliary estimators focused directly on computing this variability. This approach is also suggested in the FDA draft guidances FDA (2011) and FDA (2012), which are further discussed in the context of reference scaling in the sequel.

The pumas_be function carries out this estimation for us, yet let us also see how to do it in steps.

Separating out R and T

Consider this example:

data = dataset(joinpath("bioequivalence", "RTTR_TRRT", "PJ2017_4_3"));

Let us separate out the R and the T data, after preprocess_be:

data_preprocessed = preprocess_be(data; endpoint = :Cmax);

data_R = @rsubset data_preprocessed :formulation == 'R'
data_T = @rsubset data_preprocessed :formulation == 'T';

Using the auxiliary estimator to estimate within subject variability

We can now use Pumas’ within_subject_variability which gives us an estimate as well as a degrees of freedom for the estimate.

within_subject_variability(data_R)

(σw = 0.2094036782727373,
 k = 15.0,)

within_subject_variability(data_T)

(σw = 0.2640679869013234,
 k = 14.0,)

Here we combine these estimates and also use sigma2geocv:

σw_estimates = combine(
    groupby(data_preprocessed, :formulation),
    subdf -> within_subject_variability(
        subdf;
        userepeatedobsonly = true,
        allownonreplicated = true,
    ),
)

@rtransform! σw_estimates :CV = round(100sigma2geocv(:σw), digits = 2)
@rtransform! σw_estimates :σw = round(:σw, digits = 4)
@rselect! σw_estimates :CV :σw :k

simple_table(σw_estimates, round_mode = nothing)

CV	σw	k
21.17	0.2094	15.0
26.87	0.2641	14.0

Observe that these same results appear with pumas_be.

pumas_be(data, endpoint = :Cmax)

Observation Counts
Sequence ╲ Period	1	2	3	4
RTTR	8	8	8	8
TRRT	9	9	9	8

Paradigm: Replicated crossover that supports reference scaling
Model: Mixed model (unequal variance)
Criteria: Standard ABE
Endpoint: Cmax
Formulations: Reference(R), Test(T)
		Results(Cmax)	Assessment	Criteria
R	Geometric Marginal Mean	1150
	Geometric Naive Mean	1150
T	Geometric Marginal Mean	1049
	Geometric Naive Mean	1048
	Geometric Mean T/R Ratio (%)	91.23
	Degrees of Freedom	40.82
	90% Confidence Interval (%)	[83.18, 100.1]	Pass	CI ⊆ [80, 125]
Variability	CVᵣ (%) | σ̂ᵣ	21.17 | 0.2094
	CVₜ (%) | σ̂ₜ	26.87 | 0.264
	Variability Ratio (%)	126.2
ANOVA	Formulation (p-value)	0.1019
	Sequence (p-value)	0.5596
	Period (p-value)	0.4916

5 Conclusion

This unit provided a deeper exploration into the statistical models that power bioequivalence analysis in Pumas. We revisited the fixed-effects linear models used for standard crossover designs, extending the discussion to handle cases with multiple test formulations within a single, unified model. We then transitioned to the more complex mixed-effects models, which are essential for analyzing replicate designs. Key concepts such as the distinction between equal and unequal variance assumptions, the role of Restricted Maximum Likelihood (REML) in estimation, and the Satterthwaite method for approximating degrees of freedom were introduced to explain the machinery behind the results.

The models and methods detailed here are not just theoretical exercises; they form the critical foundation for advanced bioequivalence assessments. The ability to correctly model subject effects and estimate within-subject variability, particularly for the reference product, is a prerequisite for the scaled bioequivalence approaches required by regulatory agencies for highly variable drugs. The next unit will build directly upon this foundation, introducing Reference-Scaled Average Bioequivalence (RSABE) and demonstrating how the outputs from these mixed models, especially the within-subject variability, are used to adjust the bioequivalence acceptance criteria.

6 Unit exercises

1. Multiple Formulations in a Fixed-Effects Model Using the data_3_formulations dataset from the unit, let’s manually verify the Geometric Mean Ratio (GMR) estimates from the underlying model.

   data_3_formulations = dataset("bioequivalence/RST_RTS_SRT_STR_TRS_TSR/PJ2017_4_5")
   auc_result = pumas_be(data_3_formulations, endpoint = :AUC)

   1. Extract the coefficient vector from the model using coef(auc_result.model).
   2. Identify the coefficients corresponding to the formulation effects for S and T (Hint: look at the coeftable). Remember that R is the reference level.
   3. Using the extracted coefficients, manually calculate the GMR as a percentage for S vs. R and T vs. R by applying 100 * exp(coefficient).
   4. Compare your calculated GMRs to the GMR (%) values reported in the pumas_be output table. Do they match?

2. Degrees of Freedom and Model Rank For a fixed-effects model, the residual degrees of freedom (DOF) are crucial for inference. They are calculated as the number of observations minus the rank of the design matrix. Let’s investigate this with the auc_result model from Exercise 1, which has missing data.

   1. Determine the number of observations used in the model by finding the number of rows in auc_result.data.
   2. Determine the number of coefficients in the model using length(coef(auc_result.model)).
   3. Calculate the “naive” degrees of freedom by subtracting the number of coefficients from the number of observations.
   4. Now, find the true degrees of freedom by first getting the design matrix X = auc_result.model.mm.m and then calculating n - rank(X), where n is the number of observations. Compare this to the value from part (c).
   5. Verify your result from part (d) by calling dof_residual(auc_result.model). Why is using the matrix rank more robust than simply counting coefficients, especially when data is missing?

3. Comparing Mixed Models: Equal vs. Unequal Variances This unit showed two different mixed models for replicate designs based on the homogeneity argument. Let’s explore the difference in their variance components. Use the RTTR_TRRT dataset.

   data = dataset(joinpath("bioequivalence", "RTTR_TRRT", "PJ2017_4_3"));
   result_yes_homogenous = pumas_be(data, endpoint = :Cmax, homogeneity = true)
   result_non_homogenous = pumas_be(data, endpoint = :Cmax, homogeneity = false)

   1. Print the model summary for the equal variance case: println(result_yes_homogenous.model). Identify the estimated standard deviations (σ) for the random effects (id) and the residual.
   2. Now print the model summary for the unequal variance case: println(result_non_homogenous.model). Identify all the estimated standard deviations (σ).
   3. Based on the model outputs, explain what the homogeneity assumption means in practice. How do the random effect structures differ between the two models? Which model allows for different subject variability on the Test vs. Reference product?

4. Manual CI Calculation with Satterthwaite’s Method The confidence interval for a mixed model in pumas_be depends on the coefficient estimate, its standard error, and the Satterthwaite degrees of freedom. Let’s manually reconstruct the 90% CI for the unequal variance model (result_non_homogenous from Exercise 3).

   1. Extract the coefficient for formulation: T from result_non_homogenous.model.
   2. Extract the corresponding standard error for that coefficient using stderror(result_non_homogenous.model).
   3. Extract the degrees of freedom from the pumas_be output table: result_non_homogenous.result_table[:DF].
   4. Calculate the critical t-value using quantile(TDist(dof), 0.95), where dof is the value from part (c).
   5. Construct the 90% CI on the log scale: [coefficient ± t_value * std_error].
   6. Convert the interval to the original scale as a percentage and verify that it matches the 90% CI in the pumas_be output.

5. Estimating Within-Subject Variability pumas_be uses an auxiliary estimator to estimate the within-subject variability for each formulation in a replicate design. We can replicate this using the within_subject_variability function. Use the same RTTR_TRRT dataset.

   1. First, preprocess the data for the Cmax endpoint: data_preprocessed = preprocess_be(data; endpoint = :Cmax).
   2. Create two new DataFrames: one containing only observations for formulation R (data_R) and one for formulation T (data_T).
   3. Run within_subject_variability(data_R) and within_subject_variability(data_T). This will return a NamedTuple containing σw and the degrees of freedom k.
   4. For both R and T, take the σw value and convert it to a CV% using the sigma2geocv function.
   5. Compare your calculated CV% values to the Within-Subject Variability (CV%) for R and T reported in the output of pumas_be(data, endpoint = :Cmax). Do they match?

References

FDA. 2011. “FDA, Draft Guidance on Progesterone. Recommended Apr 2010; Revised Feb 2011.” 2011. https://www.accessdata.fda.gov/drugsatfda_docs/psg/Progesterone_caps_19781_RC02-11.pdf.

———. 2012. “FDA, Draft Guidance on Warfarin Sodium. Recommended Dec 2012.” 2012. https://www.accessdata.fda.gov/drugsatfda_docs/psg/Warfarin_Sodium_tab_09218_RC12-12.pdf.

Haubo, Rune, and B Christensen. 2020. “Satterthwaite’s Method for Degrees of Freedom in Linear Mixed Models.” 2020. https://github.com/runehaubo/lmerTestR/blob/master/pkg_notes/Satterthwaite_for_LMMs.pdf.

Liquet, Benoit, Sarat Moka, and Yoni Nazarathy. 2024. Mathematical Engineering of Deep Learning. Chapman; Hall/CRC. https://deeplearningmath.org/.

Nazarathy, Yoni, and Hayden Klok. 2021. Statistics with Julia: Fundamentals for Data Science, Machine Learning and Artificial Intelligence. Springer Nature. https://statisticswithjulia.org/.

Patterson, Scott D, and Byron Jones. 2017. Bioequivalence and Statistics in Clinical Pharmacology. Chapman; Hall/CRC. https://www.routledge.com/Bioequivalence-and-Statistics-in-Clinical-Pharmacology/Patterson-Jones/p/book/9780367782443.