Two-dimensional Gibbs sampling method based on hierarchical Bayesian model

Abstract

This invention discloses a two-dimensional Gibbs sampling method based on a hierarchical Bayesian model. This method constructs a P-function that only requires calculation of two parameters based on a known hierarchical Bayesian model, then uses two-dimensional Gibbs sampling to obtain the target parameters, and finally provides an individualized drug dosage. This significantly improves computational efficiency: since only two parameters, CL and V, need to be calculated, the calculation steps are simplified, resulting in a substantial improvement in computational efficiency. No function calculation is required; the sampling step only needs to select the average values ​​of CL and V in the population as initial values ​​and substitute them into the P-function, making the sample distribution after sampling closer to the target distribution, further improving sampling efficiency.

Description

Technical Field

[0001] This invention relates to the field of drug sampling analysis, specifically a two-dimensional Gibbs sampling analysis method based on a hierarchical Bayesian model. Background Technology

[0002] Markov Chain Monte Carlo (MCMC) is an increasingly popular method for obtaining information about distributions, particularly for estimating posterior distributions in Bayesian inference. MCMC is very useful in Bayesian inference because it focuses on the posterior distribution, which is often difficult to handle through analytical tests. Bayesian inference uses information provided by observational data about one or more parameters, in the form of likelihood, to update prior states of beliefs about one or more parameters to become posterior states of beliefs about the set of one or more parameters.

[0003] Among various pharmacokinetic methods, population pharmacokinetics (PPK) is widely used, but its application in personalized dosing studies still has limitations:

[0004] 1. The constructed PPK model has at least three parameters, making its application in personalized drug delivery computationally complex. Because many factors need to be considered during the construction of the PPK model, the final model typically contains numerous parameters, leading to extensive computation in practical applications and severely impacting computational efficiency.

[0005] 2. Excessive parameters require a large amount of patient information, which can easily lead to information leakage and poor data confidentiality. Current personalized drug administration studies do not fully utilize the information from PPK models and require substantial patient data, thus compromising patient privacy and security.

[0006] 3. The PPK model suffers from uncertainty in constructing the posterior distribution and poor data reliability, resulting in less than optimistic model accuracy. While compatible with various datasets, the PPK model's poor reliability and the potential introduction of erroneous data further complicate the accuracy of the resulting model. Summary of the Invention

[0007] To address the shortcomings of the prior art, this invention provides a two-dimensional Gibbs sampling method based on a hierarchical Bayesian model. This method constructs a P-function that only requires the calculation of two parameters based on a known hierarchical Bayesian model, then uses two-dimensional Gibbs sampling to obtain the target parameters, and finally provides an individualized drug dosage.

[0008] The technical solution to achieve the objective of this invention is:

[0009] A two-dimensional Gibbs sampling method based on a hierarchical Bayesian model includes the following steps:

[0010] (1) The P function is established based on the Bayesian hierarchical model. The P function is as follows:

[0011] p(η|y)∝p(y|η)·p(η)=P(y|η)p(τ)p(β|μ,∑ -1 )p(μ)p(∑ -1 )

[0012] In the formula: P(y|η) is the likelihood function of blood drug concentration under parameter η; p(τ) is the prior distribution of τ; p(μ) is the prior distribution of μ; p(∑ -1 ) is ∑ -1 The prior distribution;

[0013] p(β|μ,∑ -1 ) is a given parameter μ, ∑ -1 Density functions of CL and V under given conditions;

[0014] Where η=(β) i ,τ,μ,∑ -1 ); β i Let CL, V, and τ be the values ​​of y for the i-th individual. ij The accuracy of , μ is the mean of CL and V, and ∑ is the covariance matrix of CL and V;

[0015] (2) Sample CL and V using the two-dimensional Gibbs sampling method. The steps of Gibbs sampling are as follows:

[0016] 1) First, a stationary distribution π(CL, V) is given;

[0017] 2) At t=0, a random initial state (CL) is generated. (0) V (0) );

[0018] 3) In the conditional probability distribution P(V|CL) (0) Sampling in CL (0) V (1) The conditional distribution function in the acceptance / rejection sampling at this point is P(V|CL). (0) ), can take q(x) ~ U(a, b) and auxiliary uniform distribution U(0, 1);

[0019] 4) In the conditional probability distribution P(CL|V) (1) Sampling in CL (1) V (1) The conditional distribution function in the acceptance / rejection sampling at this point is P(CL|V). (1) ), can take q(x)~U(a′,b′) and auxiliary uniform distribution U(0,1);

[0020] 5) Rotate the coordinate axes until the target number of samples is reached.

[0021]

[0022] (3) Combining the P function established in step (1), the Gibbs sampling in step (2) is applied to the two-dimensional target parameters for sampling:

[0023] 1) From P(V|CL) (k) Sample V in y) (k+1) ;

[0024] 2) From P(CL|V (k+1) ,y) sampling CL (k+1) ;

[0025] 3) The initial value during simulation is θ. (0) = (0, 0), n+k samples are collected. According to the properties of Markov chains, to ensure that all samples are taken from a stationary distribution, the samples from the first k iterations are discarded, that is, the first k samples are pre-burned, and finally the target sample θ(n) = (CL) is obtained. (n) V (n) ).

[0026] The P function described in step (1) comprises three levels, the specific forms of which are as follows:

[0027] 1) P(y|η), p(τ): y obj For the observed blood drug concentration, y pre The blood drug concentration predicted by sample CL and V, i.e., y pre =function(CL, V), calculated from the room model;

[0028] P(y|η) is y obj The density function, where y obj ~N(y) pre , σ 2 )

[0029]

[0030] τ ~ Gamma(0.0001, 0.0001)

[0031] σ 2 The conjugate prior distribution of τ theoretically follows an inverted Gamma distribution, therefore τ follows a Gamma distribution, and the choice of its parameter provides a very flat, information-free distribution for τ.

[0032] 2)p(β|μ,∑ -1 The parameter vector follows a multivariate normal distribution.

[0033] β~MVN p (μ, ∑)

[0034] Among them, μ=(μ1, μ2). MVN p Let represent a p-dimensional multivariate normal distribution (where p = 2); μ1 and μ2 are the means of the parameters to be estimated, CL and V, respectively; ∑ is the variance-covariance matrix of CL and V, with the variances of CL and V on the main diagonal and the correlation coefficients of CL and V on the secondary diagonal.

[0035] 3) p(μ): μ usually follows a multivariate normal distribution, and the parameter distribution is set according to prior information.

[0036]

[0037]

[0038] in, The mean of the population's CL and V can be used for estimation; Ω is the covariance matrix of μ, which can be estimated using the covariance matrices of the population's CL and B.

[0039] p(∑ -1 ): ∑ -1 It typically follows a Wessat distribution; the parameter distribution is set based on prior information.

[0040] ∑ -1 ~Wishart(∑0,ρ)

[0041] ∑0 is the variance-covariance matrix obtained from prior information; ρ is the degree of freedom, which is usually taken as the number of parameters;

[0042] Finally, the complete joint posterior distribution of the parameters is obtained:

[0043]

[0044]

[0045] Where η=(β) i ,τ,μ,∑); It is a multivariate Gamma function (where p = 2), trace(∑ -1 ∑) is the trace of the matrix.

[0046] The advantages of this invention are:

[0047] 1. Two-dimensional Gibbs sampling significantly improves computational efficiency: Since only two parameters, CL and V, need to be calculated, the calculation steps are simpler and more convenient, resulting in a significant improvement in computational efficiency. It eliminates the need for further calculations. objIn the calculation of the function, the sampling step only needs to select the average value of CL and V in the population as the initial value, and then substitute it into the P function to make the sample distribution after sampling closer to the target distribution, thereby further improving the sampling efficiency.

[0048] 2. This method achieves highly effective confidentiality. It does not require comprehensive information about a specific patient group; it only needs to utilize a P-function established by a Bayesian hierarchical model to perform two-dimensional Gibbs sampling on this group, and can accurately predict blood drug concentrations and dosages.

[0049] 3. The model uses acquired data and prior information to derive posterior information for statistical inference, resulting in high accuracy. The P function explicitly provides the construction of the posterior distribution, demonstrating a strong theoretical foundation and high reliability. Detailed Implementation

[0050] The present invention will be further illustrated below through examples.

[0051] Example:

[0052] The article "Development of a Vancomycin Personalized Dosing Decision Support System" included studies in neonatal, pediatric, adult (including elderly), and neurosurgical patients (PPK studies). All included studies used a first-order elimination one-compartment model as the pharmacokinetic structural model. This example uses an adult patient population for the vancomycin study.

[0053] 1. Establish the P function

[0054] A function P is constructed from a hierarchical Bayesian model, containing two variables, CL and V. The function P is as follows:

[0055] p(η|y)∝p(y|η)·p(η)=P(y|η)p(τ)p(β|μ,∑ -1 )p(μ)p(∑ -1 )

[0056] P(y|η) is the likelihood function of blood drug concentration under parameter η; p(τ) is the prior distribution of τ; p(μ) is the prior distribution of μ; p(∑ -1 ) is ∑ -1 The prior distribution of p(β|μ,∑ -1 ) is a given parameter μ, ∑ -1 The density functions of CL and V under the given conditions; where η = (β i ,τ,μ,∑ -1 ); β i Let CL, V, and τ be the values ​​of y for the i-th individual. ij The accuracy is given by μ, where μ is the mean of CL and V, and ∑ is the covariance matrix of CL and V.

[0057] It is important to emphasize that the prior distributions p(τ), p(μ), and p(∑) in the P function at this point are... -1 Use the typical values ​​for the population given in the above references.

[0058] 2. Sampling Scheme Design

[0059] The two-dimensional Gibbs sampling method is used to sample CL and V. The steps of Gibbs sampling are as follows:

[0060] 1) First, a stationary distribution π(CL, V) is given;

[0061] 2) At t=0, a random initial state (CL) is generated. (0) V (0) );

[0062] 3) In the conditional probability distribution P(V|CL) (0) Sampling in CL (0) V (1) The conditional distribution function in the acceptance / rejection sampling at this point is P(V|CL). (0) ), can take q(x) ~ U(a, b) and auxiliary uniform distribution U(0, 1);

[0063] 4) In the conditional probability distribution P(CL|V) (1) Sampling in CL (1) V (1) The conditional distribution function in the acceptance / rejection sampling at this point is P(CL|V). (1) ), can take q(x)~U(a′,b′) and auxiliary uniform distribution U(0,1);

[0064] 5) Rotate the coordinate axes until the target number of samples is reached.

[0065]

[0066] 3. Application

[0067] Based on the P function established above, the Gibbs sampling in step (2) is applied to the two-dimensional target parameters for sampling:

[0068] 1) From P(V|CL) (k) Sample V in y) (k+1) ;

[0069] 2) From P(CL|V (k+1) ,y) sampling CL (k+1) ;

[0070] 3) The initial value during simulation is θ. (0)= (0, 0), n+k samples are collected. According to the properties of Markov chains, to ensure that all samples are taken from a stationary distribution, the samples from the first k iterations are discarded, that is, the first k samples are pre-burned, and finally the target sample θ is obtained. (n) =(CL) (n) V (n) ).

[0071] 4. Sampling effect

[0072] The posterior distribution and two-dimensional Gibbs sampling scheme designed in this paper were simulated using MATLAB software. First, the parameter values ​​of each prior distribution were given during the simulation. Second, iterative sampling was performed on CL and V, and the collected samples were pre-burned to ensure that all samples were obtained from stationary distributions. Descriptive indicators such as mean and variance were calculated for the collected samples. Finally, combined with the known compartmental model from the references, the predicted blood drug concentration was obtained, and the effectiveness of the established function model was evaluated, providing corresponding indicators and statistical test results.

Claims

1. A two-dimensional Gibbs sampling method based on hierarchical Bayesian model, comprising the following steps: (1) establishing a P function by a hierarchical Bayesian model, wherein the P function contains two variables CL and V, and the P function is as follows: parameters likelihood function of the blood concentration under the condition; parameters prior distribution of the parameters parameters prior distribution of the parameters parameters prior distribution of the parameters for a given parameter under the condition density function; in For the first individual , for accuracy, for The mean, for The covariance matrix; The P function contains three levels, and the three levels are as follows: 1) : Cpred is the observed blood concentration, Cpred is the blood concentration predicted from the sample CL and V, i.e. Cpred is the blood concentration predicted from the sample CL and V, i.e. is a density function of wherein wherein The conjugate prior distribution of p(x | θ) is theoretically inverse Gamma distributed, so The choice of the parameters of the Gamma distribution, which provides a very flat, uninformative distribution for τ; 2) : The parameter vector is subject to a multivariate normal distribution, as follows: wherein, . denote a multivariate normal distribution (in this case = 2); are the mean values of the parameters to be estimated; is the covariance matrix with the main diagonal being the variances of and and the secondary diagonal being the correlation coefficients of and ; 3) : Generally subject to a multivariate normal distribution, the parameter distribution is set according to prior information, as follows: wherein, The mean of CL and V of the available population is estimated; is The covariance matrix of and of the available population is estimated; and The covariance matrix of and of the available population is estimated; : Generally subject to a Wishart distribution; parameters distributed according to prior information, as follows: is a variance-covariance matrix obtained from prior information, is the degrees of freedom, usually the number of parameters; Finally, the complete joint posterior distribution of the parameters is obtained: wherein , is a polygamma function (in this case ), is the trace of a matrix; (2) sampling CL and V by using a two-dimensional Gibbs sampling method, and the steps of the Gibbs sampling are as follows: 1) Given a stationary distribution ​ 2) in case of a random initial state ; 3) sampling from conditional probability distribution The conditional distribution function in the rejection sampling is and the auxiliary uniform distribution ;​​ 4) sampling from conditional probability distribution ; the conditional distribution function in the rejection sampling at this time is , which can be taken and auxiliary uniform distribution ;​ 5) rotating the coordinate axis, and sampling samples until the target number is reached ； (3) applying the Gibbs sampling of step (2) to the two-dimensional target parameters for sampling by combining the P function established in step (1). 1) sampling from ;​ 2) from ;​ 3) Initial value taken in simulation , production samples, according to the properties of Markov chain, to ensure that the sample is taken from the stationary distribution, discard the first iteration times of samples, that is, the first sample pre-burn, the final target sample .

Images