SODF2: IV Dissolution similarity of tablets made at two manufacturing sites
Introduction
In this second example, we compare the dissolution profiles of a tablet product (SODF2) made at two different manufacturing sites, referred to as sites 1 and 2. We treat site 1 as the approved (or reference) site and site 2 as an unapproved (or test) site. We illustrate Bayesian model dependent tests of similarity based on both \(F_2\) and multi-parameter space similarity regions.
Data description
We index the SODF2 manufacturing site with \(m\) = 1 or 2. \(NL1\)= 8 lots were obtained from manufacturer 1 and \(NL2\) = 5 lots from manufacturer 2. Twelve units (tablets) from each lot were tested ( \(NU1\)= 96 and \(NU2\) = 60). Percent dissolution measurements were obtained from each unit at 1, 2, 4, 8, 16, 32, and 45 minutes. Thus, there were \(N1\) = 8x12x7 = 672 and \(N2\) = 5x12x7 = 420 total observations from these sites. The vectors \(lotm_n\) , \(unitm_n\) , \(ym_n\) , and \(xm_n\) identify the unique lot, unique unit, dissolution time (0 to 45), and observed % dissolution, respectively, for observation \(n\) from site \(m\).
Model description
Some key assumptions of this model
Conditional independence of sampling and prior events (i.e., all “~”s are iid).
Homogeneity of error variance, \(\sigma^2\) as a function of time and mean % dissolution. This is an important assumption that is partially justified in this example. However, usually variance is lower near zero and 100% dissolution and higher near 50% dissolution.
Lot and unit (tablet) effects can be approximated as linear effects on the 3 underlying Weibull coefficients (\(M\), \(\log(T)\), and \(\log(b)\)).
Covariance of lot effects are negligible. Actually this assumption has little justification. however, with only 6 lots the covariances are not well determined. So we assume they are zero.
The dissolution profile of all units (tablets) can be adequately described by a 3-parameter Weibull model. The Weibull model does have a good success history. However, there are other models that may be appropriate in various applications.
Hierarchical likelihood sampling model
The model is identical to that used in the accompanying SODF1 example, except that data from each manufacturing site are included. The data from each manufacturing site are fitted independently except that the same analytical standard deviation, \(\sigma\), is assumed for both since the dissolution testing laboratory was the same for both products.
The variable names have been necessarily modified from the SODF1 example and are slightly different than those described in BMIPR Chapter 23. The variable names include an “\(m\)” suffix which identifies the manufacturing site. The model for each site is given below.
\({{ym}_{n}}\overset{iid}{\mathop{\sim{\ }}}\,N\left( {{wbm}_{n}},{{\sigma }^{2}} \right)\)
\({{wbm}_{n}}={{\theta m }_{1,uni{{tm}_{n}}}}\cdot \left[ 1-\exp \left\{ {{\left( -\frac{{{xm}_{n}}}{\exp \left( {{\theta m}_{2,uni{{tm}_{n}}}} \right)} \right)}^{\exp \left( {{\theta m}_{3,uni{{tm}_{n}}}} \right)}} \right\} \right]\)
\({{\left( \begin{matrix} Mm \\ \log \left( Tm \right) \\ \log \left( bm \right) \\ \end{matrix} \right)}_{uni{{t}_{n}}}}={{\left( \begin{matrix} {{\theta m}_{1}} \\ {{\theta m}_{2}} \\ {{\theta m}_{3}} \\ \end{matrix} \right)}_{uni{{t}_{n}}}}=\left( \begin{matrix} {{\mu m}_{1}} \\ {{\mu m}_{2}} \\ {{\mu m}_{3}} \\ \end{matrix} \right)+{{\left( \begin{matrix} {{\beta lm }_{1}} \\ {{\beta lm }_{2}} \\ {{\beta lm }_{3}} \\ \end{matrix} \right)}_{lo{{tm}_{n}}}}+{{\left( \begin{matrix} {{\beta um }_{1}} \\ {{\beta um }_{2}} \\ {{\beta um }_{3}} \\ \end{matrix} \right)}_{uni{{tm}_{n}}}}\)
\({{\mathbf{\theta m}}_{uni{{tm}_{n}}}}=\mathbf{\mu m}+{{\mathbf{\beta lm }}_{lo{{tm}_{n}}}}+{{\mathbf{\beta um }}_{uni{{tm}_{n}}}}\)
\({{\mathbf{\beta lm }}_{i}}\sim MV{{N}_{3}}\left( \mathbf{0},\left[ \begin{matrix} {{\sigma lm}_{1}} & 0 & 0 \\ 0 & {{\sigma lm}_{2}} & 0 \\ 0 & 0 & {{\sigma lm}_{3}} \\ \end{matrix} \right] \right),i=1:{{NLm}}\)
\({{\mathbf{\beta um }}_{i}}\sim MV{{N}_{3}}\left( \mathbf{0},\Sigma_{unitm} \right),i=1:{{NUm}}\)
\(\Sigma unitm=\left[ \begin{matrix} \sigma u{{m}_{1}} \\ \sigma u{{m}_{2}} \\ \sigma u{{m}_{3}} \\ \end{matrix} \right]\cdot \left[ \begin{matrix} 1 & C{{m}_{12}} & C{{m}_{13}} \\ C{{m}_{12}} & 1 & C{{m}_{23}} \\ C{{m}_{13}} & C{{m}_{23}} & 1 \\ \end{matrix} \right]\cdot \left[ \begin{matrix} \sigma u{{m}_{1}} & \sigma u{{m}_{2}} & \sigma u{{m}_{3}} \\ \end{matrix} \right]=\mathbf{\sigma }um\cdot \mathbf{Cm}\cdot \mathbf{\sigma }u{{m}^{T}}\)
Prior model
Initially, more informative priors were employed, but these were gradually relaxed. The vague priors seemed to have little impact on the posterior estimates. Ultimately it was found that essentially “flat” priors gave essentially the same estimates as vague priors. It is possible that informative priors would be more appropriate scientifically and also might lead to better convergence properties.
\(\sigma \sim Unif\left( 0,+\infty \right)\)
\({{\mu m}_{i}}\sim Unif\left( -\infty ,+\infty \right),i=1:3\)
\({{\sigma lm}_{i}}\sim Unif\left( 0,+\infty \right),i=1:3\)
\({{\sigma um}_{i}}\sim Unif\left( 0,+\infty \right),i=1:3\)
\(\mathbf{Cm }\sim LKJcor\left( 1 \right)\)
Stan model code
# Create model for Stan
modelString = "
data {
int<lower=1> K;
int<lower=0> N1;
int<lower=1> NL1;
int<lower=1> NU1;
int<lower=1,upper=NL1> lot1[N1];
int<lower=1,upper=NU1> unit1[N1];
vector[N1] y1;
vector[N1] x1;
int<lower=0> N2;
int<lower=1> NL2;
int<lower=1> NU2;
int<lower=1,upper=NL2> lot2[N2];
int<lower=1,upper=NU2> unit2[N2];
vector[N2] y2;
vector[N2] x2;
}
transformed data {
vector[K] zero;
zero = rep_vector(0., K);
}
parameters {
real<lower=0> sigma;
cholesky_factor_corr[K] L1;
vector<lower=0>[K] sig_l1;
vector<lower=0>[K] sig_u1;
vector[K] beta_l1[NL1];
vector[K] beta_u1[NU1];
vector[K] mu1;
cholesky_factor_corr[K] L2;
vector<lower=0>[K] sig_l2;
vector<lower=0>[K] sig_u2;
vector[K] beta_l2[NL2];
vector[K] beta_u2[NU2];
vector[K] mu2;
}
transformed parameters {
vector[N1] wb1;
vector[K] theta1;
vector[N2] wb2;
vector[K] theta2;
matrix[K,K] C1 = L1 * L1';
matrix[K,K] C2 = L2 * L2';
matrix[K,K] Sigma_unit1 = quad_form_diag(C1, sig_u1);
matrix[K,K] Sigma_unit2 = quad_form_diag(C2, sig_u2);
for (n in 1:N1){
theta1 = mu1 + beta_l1[lot1[n]] + beta_u1[unit1[n]];
wb1[n] = theta1[1]*(1 - exp(-(x1[n]/exp(theta1[2]))^exp(theta1[3])));
}
for (n in 1:N2){
theta2 = mu2 + beta_l2[lot2[n]] + beta_u2[unit2[n]];
wb2[n] = theta2[1]*(1 - exp(-(x2[n]/exp(theta2[2]))^exp(theta2[3])));
}
}
model {
// Priors
//sigma ~ normal(0, 1);
//mu1[1] ~ normal(100, 320);
//mu1[2] ~ normal(0, 64);
//mu1[3] ~ normal(0, 64);
//sig_l1 ~ normal(0, 64);
//sig_u1 ~ normal(0, 64);
L1 ~ lkj_corr_cholesky(1);
//mu2[1] ~ normal(100, 320);
//mu2[2] ~ normal(0, 64);
//mu2[3] ~ normal(0, 64);
//sig_l2 ~ normal(0, 64);
//sig_u2 ~ normal(0, 64);
L2 ~ lkj_corr_cholesky(1);
// Likelihood
for (i in 1:NL1) beta_l1[i] ~ normal(0, sig_l1);
beta_u1 ~ multi_normal_cholesky(zero, diag_pre_multiply(sig_u1,L1));
y1 ~ normal(wb1, sigma);
for (i in 1:NL2) beta_l2[i] ~ normal(0, sig_l2);
beta_u2 ~ multi_normal_cholesky(zero, diag_pre_multiply(sig_u2,L2));
y2 ~ normal(wb2, sigma);
}
" # close quote for modelStringMCMC execution time
The code is set up to request 16,000 MCMC draws ( 4 chains of 4,000 draws each after warm up). This can take as long as 2 hours on a 4 core machine. Considerably fewer draws may be adequate. Some parameters appear poorly estimated and result in disappointing n_eff and Rhat (e.g., \(\sigma u1_3\)). However, the model is likely somewhat over-specified. However I feel that convergence was adequate with reliable estimation and prediction uncertainty.
I recommend modifying the code to execute the stan() function and then save the StanFit object to disk. Then in subsequent executions (if needed) to adjust the code to read the StanFit object from disk.
As always, patience, persistence, perseverance, and due-diligence are in order to assure MCMC convergence. I have found it is helpful to purchase a rubber ducky and keep it nearby as you work. It can be a source of comfort and inspiration in times of stress.
[1] "2022-09-23 16:39:25 BST"[1] "2022-09-23 16:39:30 BST"Posterior statistical summary
[1] "Table 23.5: Posterior summary for two manufacturers of SODF2." mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
sigma 2.28017066 7.025449e-04 6.206215e-02 2.163255698 2.23742922 2.27871696 2.32099255 2.40476285 7803.78790 1.001072
sig_l1[1] 2.41475494 1.065102e-02 9.251326e-01 1.280218921 1.80139038 2.21653073 2.77658711 4.73453700 7544.41175 1.000619
sig_l1[2] 0.08747982 7.597702e-04 4.146739e-02 0.028806398 0.06055267 0.08070320 0.10552064 0.18841095 2978.85166 1.000249
sig_l1[3] 0.04950760 2.180364e-04 2.108883e-02 0.022171220 0.03542551 0.04539870 0.05846073 0.10212481 9355.07447 1.000367
sig_l2[1] 5.44518040 1.504359e-01 4.471154e+00 2.089350362 3.27963694 4.40631096 6.19202696 14.97376988 883.35754 1.004873
sig_l2[2] 0.39509036 8.253042e-03 3.268176e-01 0.152495805 0.23531278 0.31329146 0.44060345 1.13635348 1568.13228 1.003097
sig_l2[3] 0.06345259 1.007974e-03 5.697027e-02 0.004187357 0.03025603 0.05088474 0.07904362 0.20203897 3194.46283 1.001520
sig_u1[1] 1.35736328 4.726255e-03 2.195613e-01 0.939765040 1.20944367 1.35082682 1.49913433 1.79824680 2158.12957 1.002936
sig_u1[2] 0.15717111 1.324626e-04 1.256546e-02 0.134534303 0.14826410 0.15647380 0.16523001 0.18367945 8998.50843 1.000126
sig_u1[3] 0.02102480 1.281161e-03 1.136480e-02 0.001411619 0.01254504 0.02087509 0.02897995 0.04398608 78.68937 1.046081
sig_u2[1] 2.57945803 2.820568e-03 3.209238e-01 2.017568358 2.35490036 2.55983303 2.78010014 3.27259280 12945.85528 1.000201
sig_u2[2] 0.14196343 1.380837e-04 1.534317e-02 0.114786489 0.13130494 0.14076951 0.15153301 0.17492592 12346.54103 1.000104
sig_u2[3] 0.08866669 2.035069e-04 1.537214e-02 0.059796450 0.07825088 0.08823686 0.09844477 0.12028818 5705.72502 1.000291
C1[1,1] 1.00000000 NaN 0.000000e+00 1.000000000 1.00000000 1.00000000 1.00000000 1.00000000 NaN NaN
C1[1,2] 0.70624704 1.913693e-03 1.016140e-01 0.487182324 0.64336911 0.71382389 0.77706867 0.88431367 2819.43746 1.003651
C1[1,3] 0.14853170 9.950780e-03 3.834284e-01 -0.599322329 -0.12749011 0.15211476 0.43214202 0.85330828 1484.75306 1.002639
C1[2,1] 0.70624704 1.913693e-03 1.016140e-01 0.487182324 0.64336911 0.71382389 0.77706867 0.88431367 2819.43746 1.003651
C1[2,2] 1.00000000 4.578939e-19 5.131565e-17 1.000000000 1.00000000 1.00000000 1.00000000 1.00000000 12559.43236 0.999750
C1[2,3] 0.45586639 7.418190e-03 3.243305e-01 -0.358520931 0.28393129 0.50176243 0.69114552 0.93554032 1911.52397 1.001038
C1[3,1] 0.14853170 9.950780e-03 3.834284e-01 -0.599322329 -0.12749011 0.15211476 0.43214202 0.85330828 1484.75306 1.002639
C1[3,2] 0.45586639 7.418190e-03 3.243305e-01 -0.358520931 0.28393129 0.50176243 0.69114552 0.93554032 1911.52397 1.001038
C1[3,3] 1.00000000 4.784209e-19 5.283211e-17 1.000000000 1.00000000 1.00000000 1.00000000 1.00000000 12194.82910 0.999750
C2[1,1] 1.00000000 NaN 0.000000e+00 1.000000000 1.00000000 1.00000000 1.00000000 1.00000000 NaN NaN
C2[1,2] 0.36807151 1.197018e-03 1.336909e-01 0.085009391 0.28188831 0.37566302 0.46306393 0.60617986 12473.90228 1.000010
C2[1,3] -0.14239190 2.094871e-03 1.911487e-01 -0.497730810 -0.27689751 -0.14864541 -0.01537035 0.24669587 8325.83606 1.000362
C2[2,1] 0.36807151 1.197018e-03 1.336909e-01 0.085009391 0.28188831 0.37566302 0.46306393 0.60617986 12473.90228 1.000010
C2[2,2] 1.00000000 6.909260e-19 8.677163e-17 1.000000000 1.00000000 1.00000000 1.00000000 1.00000000 15772.20630 0.999750
C2[2,3] -0.51004753 1.404582e-03 1.401675e-01 -0.748983319 -0.61065500 -0.52228036 -0.42029909 -0.20380040 9958.65003 1.000241
C2[3,1] -0.14239190 2.094871e-03 1.911487e-01 -0.497730810 -0.27689751 -0.14864541 -0.01537035 0.24669587 8325.83606 1.000362
C2[3,2] -0.51004753 1.404582e-03 1.401675e-01 -0.748983319 -0.61065500 -0.52228036 -0.42029909 -0.20380040 9958.65003 1.000241
C2[3,3] 1.00000000 9.048801e-19 7.444731e-17 1.000000000 1.00000000 1.00000000 1.00000000 1.00000000 6768.86636 0.999750
mu1[1] 95.57517761 1.668159e-02 9.558253e-01 93.715785423 95.01577926 95.56262236 96.12219737 97.48783404 3283.08443 1.002009
mu1[2] 1.88329995 5.539943e-04 3.701418e-02 1.808017780 1.86141766 1.88330704 1.90494729 1.95778940 4464.01413 1.001386
mu1[3] 0.54360908 2.900954e-04 2.056420e-02 0.502194135 0.53178729 0.54368020 0.55568525 0.58442917 5025.06724 1.000151
mu2[1] 95.75226750 1.241873e-01 3.212510e+00 89.807536471 94.36059805 95.72847901 97.08067699 101.22728362 669.16771 1.006042
mu2[2] 1.56043740 7.020145e-03 2.266532e-01 1.132841981 1.46427397 1.55927301 1.65740267 1.99766039 1042.39302 1.001799
mu2[3] 0.50412983 7.496008e-04 3.943456e-02 0.425603647 0.48464568 0.50360980 0.52302921 0.58509991 2767.54094 1.001870Comparison of observed and predicted % dissolution
The dissolution profiles for each unit are shown in Figure 23.5 below. Each of the 13 panels includes dissolution profiles for 12 tablets with symbols indicating individual observed dissolution values and lines joining corresponding Weibull model based predicted posterior mean values.
Posterior probability of similarity relative to the spheroidal \(F_2\) similarity region
As described in the text, only the first \(N_{times}=5\) time points (1, 2, 4, 8, and 16 minutes) are employed for the \(F_2\) similarity test. The condition of similarity is
\({{F}_{2}}=50\cdot {{\log }_{10}}\left( \frac{100}{\sqrt{1+\tfrac{1}{{{N}_{times}}}\sum\limits_{i=1}^{{{N}_{times}}}{{{\left( wb{{2}_{i}}-wb{{1}_{i}} \right)}^{2}}}}} \right)>50.\)
Letting \(\Delta wb_i=wb{{2}_{i}}-wb{{1}_{i}}\), this condition is equivalent to
\(R=\sqrt{\sum\limits_{i=1}^{{{N}_{times}}}{\Delta wb _{i}^{2}}}<\sqrt{99\cdot {{N}_{times}}}\)
where \(R\) is the radius of an \(N_{times}\) dimensional hypersphere centered at the origin (no difference in true mean dissolution profiles).Notice that the definition of similarity now depends critically on the number of dissolution time points included in the calculation with larger values of \(N_{times}\) allowing greater differences at some time points to be compensated by smaller differences at others. In this sense a definition of similarity based on \(F_2\) seems incoherent. Never-the-less, it continues to serve as a regulatory standard in leu of bio-equivalence.
Figure 23.7 provides a kernal density estimate of the posterior distribution of \(F_2\). The estimated posterior probability of similarity, the event \(F_2 > 50\), is thus 0.396, obtained as a simple counting exercise.
There is no universal regulatory standard for the minimum posterior probability of similarity, not should their be. From a scientific and risk management perspective, this minimum should depend on the consequences of non-similarity and decision errors in any particular application. However, this result seems low. Figure 23.6 below shows that the grand mean dissolution difference between sites at 5 minutes is greater than 10% which is consistent with the low posterior probability of similarity.
Posterior probability of similarity relative to a hyper-rectangular similarity region
An alternative definition of similarity is based on the proposal of Sathe et al (1996) who defined a hyper-rectangular region \(SR1\) of similarity as.
\(SR1=\left( -3{{{\hat{S}}}_{1}}\le \Delta {{\mu }_{1}}\le 3{{{\hat{S}}}_{1}} \right)\cap \left( -3{{{\hat{S}}}_{2}}\le \Delta {{\mu }_{2}}\le 3{{{\hat{S}}}_{2}} \right)\cap \left( -3{{{\hat{S}}}_{3}}\le \Delta {{\mu }_{3}}\le 3{{{\hat{S}}}_{3}} \right)\)
where \({{\hat{S}}}_{p}\) are univariate standard deviation estimates that include both between and within lot variation estimated for the pth parameter from reference site 1. We can estimate these limits from the posterior draws as follows:
\({{\hat{S}}_{p}}=\tfrac{1}{{{N}_{draws}}}\sum\limits_{i=1}^{{{N}_{draws}}}{\sqrt{{{\left( \sigma l1_{p}^{[i]} \right)}^{2}}+{{\left( \sigma u1_{p}^{[i]} \right)}^{2}}}},p=1,2,3\)
Here, the superscript \([i]\) indicates the \(i^{th}\) posterior draw of the between-lot or within-lot variances respectively for the \(p^{th}\) parameter from the reference site 1.
Again, using the simple counting exercise to determine the proportion of posterior draws lying within the hyper-rectangular similarity region \(SR1\), we find the posterior probability of similarity based on the \(SR1\) definition to be 0.8870625. This differs considerably from the posterior probability based on the \(F_2\) definition, but there is no reason to expect agreement because these definitions of similarity are different. In particular, the factor of 3 used in the \(SR1\) definition is quite liberal.
A criticism of the \(SR1\) definition of similarity is that it is based on estimated (plug-in) values of \({{\hat{S}}}_{p}\). Thus the definition of similarity will be different for each new set of data. This does not meet the normative requirements for a coherent similarity definition outlined by Eaton et al (2003). A definition of similarity based on the required properties of underlying model parameters would be preferable from a perspective of safety, efficacy, and life-cycle quality control. A simple modification of the above procedure would be to set fixed marginal limits on each of the three Weibull parameters based on the fit-for-purpose requirements of the pharmaceutical product. In principal, the above procedure can accommodate any similarity region whose shape can be expressed algorithmically. This is because the posterior probability of similarity is estimated using a simple counting exercise. It is only necessary to determine whether each MCMC draw is, or is not, contained within the defined region.