Simultaneous multifactor Bayesian analysis (SiMBA) of PET time activity curve data

Positron emission tomography (PET) is an in vivo imaging method essential for studying the neurochemical pathophysiology of psychiatric and neurological disease. However, its high cost and exposure of participants to radiation make it unfeasible to employ large sample sizes. The major shortcoming of PET imaging is therefore its lack of power for studying clinically-relevant research questions. Here, we introduce a new method for performing PET quantification and analysis called SiMBA, which helps to alleviate these issues by improving the efficiency of PET analysis by exploiting similarities between both individuals and regions within individuals. In simulated [11C]WAY100635 data, SiMBA greatly improves both statistical power and the consistency of effect size estimation without affecting the false positive rate. This approach makes use of hierarchical, multifactor, multivariate Bayesian modelling to effectively borrow strength across the whole dataset to improve stability and robustness to measurement error. In so doing, parameter identifiability and estimation are improved, without sacrificing model interpretability. This comes at the cost of increased computational overhead, however this is practically negligible relative to the time taken to collect PET data. This method has the potential to make it possible to test clinically-relevant hypotheses which could never be studied before given the practical constraints. Furthermore, because this method does not require any additional information over and above that required for traditional analysis, it makes it possible to re-examine data which has already previously been collected at great expense. In the absence of dramatic advancements in PET image data quality, radiotracer development, or data sharing, PET imaging has been fundamentally limited in the scope of research hypotheses which could be studied. This method, especially combined with the recent steps taken by the PET imaging community to embrace data sharing, will make it possible to greatly improve the research possibilities and clinical relevance of PET neuroimaging.


Supplementary Materials S1 : Prior Specification
We made use of a multifactor model (Betancourt, 2021), with a global mean, together with partially pooled deviations from the global mean for each region across all individuals, and for each individual across all regions. We also added a final partially pooled parameter to accommodate any residual differences at the level of individual TACs after accounting for average regional and individual effects. We therefore anticipate that most of the variation in the pharmacokinetic parameters will be accounted for by individual differences and regional differences, and that TAC-level will be much smaller. For this reason, we defined wider zero-centred priors over the standard deviation across individuals and regions, but more constrainted priors over the standard deviation across TACs.
For the simulations, we did not include any covariates (other than region where we used unpooled effects for BP ND and K 1 ) for any parameters besides BP ND , which was compared between the two groups.

Global Intercepts
Below are the priors defined for the global intercepts. Note that all priors are defined over the natural logarithms of the parameters.

Individual effects
Differences between individuals were defined by specifying the primary pharmacokinetic parameters in one variance-covariance matrix, and v B and σ on their own.

Regional effects
For logBP ND and logK 1 , regional differences were defined as unpooled effects using a dummy (indicator) variable defined with reference to the dorsolateral prefrontal cortex as covariates. For simplicity, all regional differences were defined as zero-centred regularising priors with the same SD.
For the remaining parameters, regional differences were defined as pooled variables, arising from a common distribution

TAC effects
Residual TAC effects were defined for the major four pharmacokinetic parameters, but were not included for v B or σ as these were not considered to be of central importance.

Smooth sigma function
The smooth basis function makes use of a penalised regression spline. The σ spline-coefficients term refers to the standard deviation of the spline coefficients, which penalises the wiggliness of the spline. The α spline-coefficients term refers to the fixed effects term, which describes the magnitude of the influence of the smooth term around the estimated mean value.

Covariates
For the simulations, there was only one additional covariate for group effects on BP ND . For this, we used a zero-centred regularising prior. This prior conservatively assumes that no difference is most likely, and assigns 95% of its probability between differences of -48% and 48% differences between groups.

Analysis data
For the analysis, we included more covariates due to variation in the data.

Covariates
We included effects of age and sex covariates for the estimation of K 1 , V ND and BP ND . For age, we divided ages by 10, so that it was not defined by years, but by decades. This makes it easier to specify priors, and is expected to improve efficiency of estimation by estimating parameters in a similar order of magnitude.
Covariates for the measurement error included the centred natural logarithm of the average region size, as well as the centred natural logarithm of the injected radioactivity.
Because so much of the variation in the measurement error between regions would be explained by region size, we reduced the variance of the prior for the partial pooling of measurement error across regions.

Standard deviations
Parameter SD

Mean values
Here are the regional means. For simplicity of interpretation, I have added the global mean value to all of the deviations from the mean so that they are displayed as the values themselves. Region

Correlation Matrix
Parameter

Smooth sigma function
The smooth sigma function across time is depicted below.

Supplementary Materials S3: Comparisons between power estimation approaches
In Figures 1 and 2, we compare the outcomes estimated using the following datasets: 1. Large Empirical: 1000 simulated studies, with power and false positive rate calculated empirically. These datasets were only used with NLS estimation. 2. Small Empirical: 50 simulated studies, with power and false positive rate calculated empirically. 3. Small Modelled: 50 simulated studies, with power and false positive rate calculated using logspline functions applied to the confidence or credible intervals. Figure 1: Power estimated using the same and different datasets between methods, with power calculated with different methods. All methods produce similar results. Error bars are not included for NLS: t-test figures for clarity.
In Figure 3 are the fits of the logspline function to the SiMBA 95% credible intervals for the calculation of the small empirical power and false positive rates.

Supplementary Materials S4: Group difference properties in the same data
Calculated in the same set of 50 datasets for each condition, we examined the mean of the 50 estimated differences, the standard deviation of the means across these datasets, the mean standard error of these estimated means across the datasets, and the (modelled) power and false positive rates. In the Figures 4 and  5, we show plots in the style of the power and false positive plots of the standard deviation and standard error of these methods in the same datasets. Figure 4: Standard deviation of estimated group differences across simulated datasets.
The following table includes all the data together.   Figure 6: The standard deviation of estimated mean group differences for different sample sizes as a measure of consistency. With larger samples, the standard deviation of the estimates is decreased, indicating greater consistency of estimates. Linear mixed effects (LME) modelling of binding outcomes show less variation than t-tests, and the Bayesian hierarchical pharmacokinetic modelling of TACs has greatly reduced variation compared to both LME modelling and t-tests.

Supplementary Materials S6: Sensitivity to measurement error
In Figure 7, we show the power (estimated using logsplines and bootstrapped 95% confidence intervals), the mean posterior estimates of α σ , the standard deviation (SD) of the mean posterior estimates of the group differences (with 95% confidence intervals estimated using a case resampling bootstrap procedure), and the mean standard error (SE) of the differences from within each simulated study (with 95% quantiles). All values are calculated for BP ND and BP P (as these two are the same due to not estimating different V ND values between groups).
Figure 7: The SiMBA model is sensitive to measurement error, but even with the standard deviation of the measurement error equal to 20% of the mean TAC value, the power with n=10 is still over 20%.

Supplementary Materials S7: Application to data with no parameter intercorrelations
Data were simulated with the originally estimated parameter intercorrelations, as well as following sampling from univariate distributions, thereby removing the correlations between parameters. We examined the distribution of the estimated parameter intercorrelations in the original, and in the uncorrelated data.
Comparisons between correlated and uncorrelated data are shown in figures 8 and 9.
Also displayed are the estimated correlations in the original data for different sample sizes in figures 10 and 11.
In all cases, the measurement error of the simulated data is 10%, which is more than twice that of the original dataset.