Let j = 1 , … , J index the HER2 groups. Within each group j, let y jk denote the average tumor proliferation outcome across mouse replicates under the treatment ( y jk 1 ) and control ( y jk 0 ) for patient k = 1 , … , n j . Specifically, y jk 1 is the average treatment-arm response for patient k in HER2 group j: for PDX it is the mean tumor volume at Day 21 across the 3 treatment-arm mice, and for MiniPDX it is the mean RLU at Day 10 across the 6 treatment-arm capsules (3 mice × 2 capsules each). y jk 0 is the corresponding vehicle-arm (control) average computed the same way. Each patient therefore contributes both a treatment arm and a vehicle arm, and the two arms are paired within patient: y jk 1 and y jk 0 come from the same tumor source under the same protocol, with allocation to cohorts determined by HER2 expression rather than by within-cohort randomization.

Let z jk denote the outcome corresponding to TGI or TCGI, as defined in (2.1) and (2.2), respectively. It follows that z jk = ( 1 − y jk 1 y jk 0 ) × 100 % . Here z jk is the patient-specific drug response, summarizing relative tumor inhibition for patient k in HER2 group j: z jk > 1 corresponds to net tumor regression, z jk = 1 to tumor stasis (no net growth from baseline), 0 < z jk < 1 to partial growth inhibition, z jk = 0 to no treatment effect, and z jk < 0 to tumor escape. We take z jk ≥ 0.45 as a positive-response threshold consistent with the per-tumor TCGI cutoff used in the MiniPDX literature; this corresponds to roughly 55 % of the control growth being suppressed in the treatment arm. We propose a Bayesian hierarchical model for analyzing the data { z jk }. We assume (2.3) z jk ∼ N ( z jk ; μ j , σ j 2 ) where N ( ⋅ ; μ j , σ j 2 ) is a normal distribution with mean μ j and variance σ j 2 . Moreover, to enable shrinkage and partial pooling, we specify the patient-level parameters to follow hyperpriors (2.4) ( μ j , σ j 2 ) ∼ N ( μ j ; μ 0 , τ 2 ) ⋅ TCauchy ( σ j ; a , b , ( 0 , ∞ ) ) . Here, TCauchy ( ⋅ ; a , b , ( 0 , ∞ ) ) is a half-Cauchy distribution. In addition, let (2.5) ( μ 0 , τ 2 ) ∼ N ( μ 0 ; s 0 , r 0 2 ) ) ⋅ TCauchy ( τ ; a τ , b τ , ( 0 , ∞ ) ) where s 0 , r 0 2 , a, b, a τ and b τ are specified to enable vague but realistic hyperpriors, ensuring flexibility without imposing strong assumptions.

The prior choices in equations (2.3)–(2.4) are weakly informative. We use a half-Cauchy prior on σ j and on the between-patient scale τ, following the recommendation of Gelman [4] for variance components in hierarchical models with a small number of groups, which gives a heavier-tailed alternative to a half-normal and behaves well when the number of groups is small. We use a vague normal prior on μ 0 centered at a mid-range TCGI value and with realistic spread, so that the prior is non-informative within the support of the data without sending posterior mass to implausible negative-inhibition extremes.

Let i = 1 , 2 index experimental assays (MiniPDX and PDX), j = 1 , … , J index HER2 expression groups, and k = 1 , … , n ij index patients within each assay-by-group combination. We propose a Bayesian hierarchical model designed to integrate information from both the MiniPDX and PDX assays simultaneously. This structure facilitates borrowing strength across assays, allowing the MiniPDX assay to leverage richer information from the more established PDX data. Specifically, we define assay-specific effects as α i and HER2 group-specific effects as β j , thus capturing systematic variations attributable to assay conditions and biological differences between HER2 groups.

Formally, our hierarchical model for the TCGI or TGI ( μ ij ) outcome z ijk is specified as: (2.6) z ijk ∼ N ( α i + β j , σ ij 2 ) k = 1 , … , n ij , i = 1 , 2 , j = 1 , … , J , where the mean structure explicitly partitions the expected outcome into assay-specific( α i ) and group-specific( β j ) contributions and define μ ij = α i + β j to be the total effect. Without loss of generality, we use μ ij to represent the treatment effect in both single-assay and multi-assay contexts. To ensure meaningful interpretation and identifiability, we impose the sum-to-zero constraint ∑ i = 1 I α i = 0 on the assay effects.

We further model the variance parameters hierarchically. To stabilize variance estimation and allow for flexible heterogeneity within and between assays, we introduce a hierarchical log-variance structure: (2.7) log ( σ ij 2 ) | log ( σ i 2 ) ∼ N ( log ( σ i 2 ) , τ within 2 ) , log ( σ i 2 ) | log ( σ 0 2 ) ∼ N ( log ( σ 0 2 ) , τ between 2 ) where τ within 2 and τ between 2 captures the variability within assay and between assay, and σ 0 ∼ TCauchy ( a σ 0 , b σ 0 , ( 0 , ∞ ) ).

At the top level, we place hyperpriors on the assay-level parameters α i and group-level mean parameters β j , modeling their variances directly: (2.8) α i | τ α 2 ∼ N ( α 0 , τ α 2 ) , (2.9) β j | β 0 , τ β 2 ∼ N ( β 0 , τ β 2 ) , (2.10) with constraint: ∑ i = 1 I α i = 0 .

Since we have imposed the constraint on the assay effects, we do not assign a hyperprior to α 0 , but instead fix its value at 0. To robustly model uncertainty around hyperparameters, we employ weakly informative hierarchical priors. Specifically, we set: (2.11) τ α ∼ TCauchy ( a α , b α , ( 0 , ∞ ) ) , (2.12) ( β 0 , τ β 2 ) ∼ N ( s β , r β 2 ) ⋅ TCauchy ( τ β ; a β , b β , ( 0 , ∞ ) ) where s β , r β 2 , a α , b α , a σ 0 , b σ 0 , a β and b β are specified to enable vague but realistic hyperpriors, ensuring flexibility without imposing strong assumptions.

The two-assay model imposes the sum-to-zero constraint ∑ i α i = 0 to make the assay effects identifiable under the additive coding μ ij = α i + β j . The hierarchical log-variance specification log ( σ ij 2 ) ∼ N ( log ( σ i 2 ) , τ within 2 ) and log ( σ i 2 ) ∼ N ( log ( σ 0 2 ) , τ between 2 ) stabilizes variance estimation when each assay-by-group cell has only a handful of patients. As in the single-assay model, the variance scales τ within , τ between , σ 0 , τ β carry weakly informative half-Cauchy priors per Gelman [4], and the top-level mean β 0 has the same vague-normal treatment we used for μ 0 in the single-assay model. The assay-effect SD τ α is not identifiable from the data with only I = 2 assays under the sum-to-zero constraint and is therefore fixed at τ α = 2.5 in the present analysis (see Appendix B for the sensitivity check).

We fit all Bayesian hierarchical models in R 4.4 via the rjags package (version 4.x), which interfaces to JAGS for Markov chain Monte Carlo (MCMC) sampling. For both the single-assay and two-assay models we run 3 parallel chains with an adaptation phase of 10 , 000 iterations, a burn-in of 10 , 000 iterations discarded, and 100 , 000 retained iterations per chain thinned by 10, yielding 10 , 000 effective post-thin draws per chain and 30 , 000 draws in total. Random-number seeds 1001, 1002, and 1003 are used for the three chains so that all results in this paper are reproducible. Convergence is assessed using the Gelman–Rubin potential scale reduction factor R ˆ (target < 1.05) and effective sample size, computed via coda::gelman.diag and coda::effectiveSize on the saved 3-chain mcmc.list, and trace plots are inspected for all monitored parameters. Posterior probabilities of the form Pr ( μ ij ≥ t 0 | data ) are computed as the empirical proportion of post-burn-in, thinned MCMC draws satisfying the inequality, and credible intervals are reported as 95% equal-tailed quantile intervals of the posterior draws.

An initial and usually small number of tumor samples is often used to investigate therapeutics effects using MiniPDX or PDX models. After the data are observed and treatment effects estimated based on the proposed Bayesian hierarchical model, it may be the case that treatment effects for some tumors may exhibit large posterior variabilities and require additional investigation. For others, the treatment effects might be more homogeneous without further exploration. Therefore, we consider an adaptive design and sample size re-estimation (SSR) to allow additional tumor samples to be tested. We define efficacy as the posterior probability greater than a threshold η that the mean TCGI or TGI ( μ ij ) is above a pre-determined value t 0 .

The proposed design employs a Bayesian adaptive procedure to efficiently estimate the required number of additional patients. To guide this adaptive enrollment, we simulate future patient outcomes assuming the future TCGI or TGI values follow some distribution f, which reflects a realistic anticipated effect. Importantly, as new information becomes available, we can adjust this assumed distribution accordingly, underscoring the advantage of the Bayesian adaptive framework.

We iteratively update the posterior probability Pr ( μ ij ≥ t 0 | data ) with each simulated set of additional patients. Enrollment simulations continue until this posterior probability surpasses our predefined threshold η. The detailed adaptive sampling procedure is described in Algorithm 1. Algorithm 1

Sample size re-estimation based on observed data.

The dataset consists of measurements evaluating the efficacy of trastuzumab deruxtecan relative to a control in inhibiting tumor growth. Tumors are categorized into three biomarker subgroups based on HER2 expression: High, Medium, and Low. Within each subgroup, tumor tissues from different patients underwent MiniPDX and PDX assays.

The sample flow for the analyses below reflects a two-stage adaptive design and is summarized in Figure 2. In Stage 1 we enrolled 16 PDX models (Table 1, batch one), each assayed on both the PDX (TGI) and MiniPDX (TCGI) platforms. Two of the Stage 1 models have HER2-negative IHC status (HER2 −) and are not included in the Bayesian hierarchical analyses; the mechanism of action of trastuzumab deruxtecan is HER2-driven, so the no-expression cohort is not expected to share the same response distribution as the HER2-positive strata, and pooling it into the same hierarchical structure would mis-specify the borrowing across groups. We retain the HER2-negative cases in the design schematic (Figure 1) and in Table 1 for completeness of the enrolled cohort, but they do not enter equations (2.3)–(2.4) or the two-assay model. The 14 Stage 1 HER2-positive models (HER2 +, ++, +++) contribute paired PDX TGI and MiniPDX TCGI measurements (Table 2) and form the primary analysis cohort for the single-assay (Section 2.2) and two-assay (Section 2.3) hierarchical models. After the Stage 1 analysis, the MiniPDX HER2-high posterior probability of efficacy was Pr ( μ 1 ≥ 0.45 | data ) = 0.853, below our pre-specified target of 0.9, whereas the PDX HER2-high probability was already ≈ 1.00. The MiniPDX assay was therefore the bottleneck for confidence, and the adaptive sample-size re-estimation procedure described in Section 2.5 was triggered specifically for MiniPDX. Under that procedure we then prospectively enrolled 4 additional HER2-high models in Stage 2 (Table 1, batch two), assayed on the MiniPDX platform only (no PDX TGI was collected for Stage 2; Table 7 reports the MiniPDX TCGI values). These 4 Stage 2 models enter the analysis only in Section 3.4, where they illustrate the MiniPDX sample-size re-estimation under the adaptive design. Figure 2

Sample-flow diagram for the analyses. Stage 1 enrolled 16 PDX models, each assayed on both the PDX (TGI) and MiniPDX (TCGI) platforms (Table 1, batch one). The 14 HER2-positive Stage 1 models (HER2 +, ++, +++) form the primary analysis cohort for the single-assay and two-assay hierarchical analyses; the 2 HER2-negative Stage 1 models are excluded from the hierarchical analyses because the trastuzumab deruxtecan mechanism is HER2-driven and the no-expression cohort is not expected to share the same response distribution. After Stage 1, the MiniPDX HER2-high posterior probability of efficacy was Pr ( μ 1 ≥ 0.45 | data ) = 0.853, below the 0.9 threshold, while the PDX HER2-high probability was already ≈ 1.00; the SSR was therefore triggered specifically for the MiniPDX assay (Section 2.5), and Stage 2 prospectively enrolled 4 additional HER2-high models on MiniPDX only (no PDX TGI was collected for Stage 2; Table 1, batch two). These 4 Stage 2 models form the SSR cohort analyzed in Section 3.4. The excluded cohort is shown with a dashed border.

We applied the Bayesian hierarchical model presented in Section 2.2 to analyze the MiniPDX and PDX datasets separately. We set the parameters s 0 = 0.5, r 0 = 1, a = a τ = 0, b = b τ = 2.5 for the single assay Bayesian hierarchical model.

Based on the results from one assay: across all HER2 groups, the posterior mean TGI values derived from the PDX assay are consistently higher compared with those obtained from the MiniPDX assay. Moreover, the credible intervals for the HER2 High and Medium groups from the PDX assay are narrower, reflecting lower uncertainty. Therefore, leveraging a two assay hierarchical model, as described in Section 2.3, is particularly beneficial. After separately modeling the MiniPDX and PDX data, we apply the model described in Section 2.3 to borrow information across both groups and assays.

For the two assay model, τ within and τ between are estimated from the data (see Appendix A for details), and we set the hyperparameters as s β = a β = a σ 0 = 0, r β = 1 and b β = b σ 0 = 2.5. With only I = 2 assays, the assay-effect SD τ α is not identifiable from the data; we therefore fix τ α = 2.5 rather than sampling it under a half-Cauchy prior. The justification and a sensitivity check at τ α ∈ { 0.5 , 1 , 2.5 , 5 } are reported in Appendix B.

The results are presented in Table 5 and Figure 7. By jointly analyzing the PDX and MiniPDX data, we observe a notable reduction in the posterior variance of μ ij , reflecting that borrowing information across assays improves the precision of the inference. In particular, the posterior probability Pr ( μ 1 j ≥ 0.45 | data ) for the HER2 High group increases from 0.853 (when using only MiniPDX data) to 0.958 under the combined model, indicating substantially greater confidence in the treatment effect for this subgroup. Overall, these findings highlight the advantage of hierarchical modeling in leveraging information across assays to yield more robust and reliable estimates. Convergence diagnostics for the two-assay production fit are within standard bounds for all parameters of scientific interest ( R ˆ ≤ 1.001, ESS > 17 , 000 for α i , β j , β 0 , and μ ij ); the worst R ˆ across the model is 1.029 on a single σ ij 2 cell, well within the < 1.05 target. Table 5

Posterior inference for μ ij under different HER2 groups for MiniPDX and PDX data.

To make the efficiency contribution of the Bayesian hierarchical model concrete, we contrast its posterior estimates with simple per-group frequentist sample means and 95% Student-t confidence intervals computed from the same data shown in Table 2. Table 6 reports posterior mean, 95% credible-interval half-width (the half-length of the central 95% posterior quantile interval, computed from the saved MCMC draws), sample mean, and 95% confidence-interval half-width ( t n − 1 , 0.975 × s ∕ n ) for each HER2 group in each assay, under both the single-assay model of Section 2.2 and the joint two-assay model of this section. In the MiniPDX assay, where per-group sample sizes are small and the data are heterogeneous, the single-assay Bayesian credible intervals are already narrower than the frequentist confidence intervals across all three HER2 groups, and the two-assay model tightens them further by roughly a factor of two (for example, the HER2 High half-width drops from 0.637 to 0.296). In the PDX assay the single-assay Bayesian and frequentist widths are comparable, while the two-assay model again tightens the credible intervals across groups (the HER2 Medium half-width drops from 0.124 to 0.102). This pattern is consistent with the intuition that hierarchical shrinkage helps most when per-group information is weakest, and that joint modeling across the two assays buys further efficiency on top of within-assay borrowing. Table 6

Single-assay Bayesian, two-assay Bayesian, and per-group frequentist estimates of the mean response μ ij in each HER2 group (HER2 High/Medium/Low), for the MiniPDX and PDX assays. The two-assay Bayesian rows use the joint hierarchical model presented above. Sample sizes are n = 4 for MiniPDX High, MiniPDX Medium, and PDX High, and n = 5 for MiniPDX Low, PDX Medium, and PDX Low. The frequentist columns are identical in the single-assay and two-assay blocks because they depend only on the raw per-group data.

We confirmed the robustness of the two-assay inference along two prior axes: the variance multiplier on ( τ within , τ between ) swept over { 1 , 2 , 5 , 10 }, and the fixed assay-effect SD τ α swept over { 0.5 , 1 , 2.5 , 5 }. Across both sweeps the cell-mean posterior estimates μ ij shift by at most ≈ 0.03 on the TGI/TCGI scale, the rank ordering High > Medium > Low is preserved within each assay, and the HER2-High threshold-exceedance probability Pr ( μ 1 , High ≥ 0.45 | data ) remains ≥ 0.95 in MiniPDX and ≈ 1.00 in PDX. Full numerical tables and the per-run prior scales are reported in Appendix B.

After completing the one & two assays analyses, we conducted a sample size re-estimation for the MiniPDX data as a demonstration. Based on initial findings, the drug’s efficacy is most promising in the HER2 High subgroup for the MiniPDX assay, although this subgroup also shows the largest variability given the unusual value of Tumor 1. The posterior probability Pr ( μ 1 ≥ 0.45 | data ) = 0.853 at the end of stage 1, and we set a target posterior probability of 0.9 for concluding that the drug is efficacious in this subgroup. Therefore, following Algorithm 1, we conducted adaptive sample size re-estimation and the parameters were established as f = InvGamma ( 2 , 5 ) − 1, with η = 0.9, t 0 = 0.45, and n sim = 1000. We then calculate the values of { n ( 1 ) , … , n ( n sim ) } for the HER2 High group. Figure 8 shows that the simulated extra number n ∗ of subjects needed is highly skewed, ranging from 1 to 138, with the kernel-density mode at 1.49. Because n ∗ is a discrete sample size that takes only integer values ≥ 1, the kernel-density mode at 1.49 is a smoothing artifact and is not directly interpretable as a sample-size recommendation; we therefore summarize the simulated distribution by its empirical quantiles. The median of n ∗ is 3, the 25th and 75th percentiles are 1 and 9, the 90th percentile is 18, and the empirical mass at small integers is Pr ( n ∗ = 1 ) = 0.292, Pr ( n ∗ ≤ 2 ) = 0.419, and Pr ( n ∗ ≤ 3 ) = 0.506. The distribution is dominated by the smallest discrete values with a long right tail (95th percentile at 25, only 0.4 % of simulations exceed 50). The median value n ∗ = 3 is the most defensible single-number recommendation under this design; we also note that the modal integer is n ∗ = 1, while a slightly more conservative integer choice n ∗ = 2 adds a small safety buffer over the mode without committing to the long-tail upper quantiles. In the actual study we prospectively collected 4 additional HER2-high patients (Table 7), which sits between the median and the upper quartile of the simulated distribution; these four observations were collected after the stage-1 analysis indicated that additional data were needed under the adaptive design, and they are not retrospectively chosen for illustration. We acknowledge that summarizing n ∗ by quantiles, rather than by a kernel-density mode of a discrete distribution, is the standard reporting we would recommend in any prospective application of this procedure. Figure 8

Density for n(sim) using Inverse Gamma prior.