Abstract:
There are provided a system and a method therefor for analyzing pharmacokinetics in such a manner that the influence of genetic polymorphism of an individual is taken into consideration, and a system and a method therefor for allowing implementation of high-accuracy haplotype frequency estimation and diplotype configuration estimation even in a case where the number of individuals is small from which data is obtainable when making pharmacokinetic analysis. A diplotype configuration estimation step estimates diplotype configurations of individuals. Next, a pharmacokinetic model building step expresses pharmacokinetic parameters as functions of the diplotype configurations of the individuals, thereby configuring pharmacokinetic models. Moreover, a pharmacokinetic parameter estimation step estimates the pharmacokinetic parameters. Here, a diplotype configuration correction step corrects the estimation result of the diplotype configurations of the individuals, thereby implementing an enhancement in the estimation accuracy.

Description:
BACKGROUND OF THE INVENTION  
       [0001]     1. Field of the Invention  
         [0002]     The present invention relates to a pharmacokinetic analysis system and a method of the pharmacokinetic analysis for analyzing transition of in-blood drug concentration of an individual in such a manner that the influence of genetic polymorphism is taken into consideration.  
         [0003]     2. Description of the Related Art  
         [0004]     A drug dosed to humans passes through processes of liberation, absorption, distribution, metabolism, and excretion. The research field regarding transition of drug concentration in body fluid (mainly in blood) in these processes is referred to as “Pharmacokinetics: PK)”. When investigating pharmacokinetics in humans, data on the pharmacokinetics, such as in-blood drug concentration or in-urine excretion amount obtained from a clinical trial, is analyzed using statistics. After that, the pharmacokinetics is investigated based on this analysis result.  
         [0005]     When applying the pharmacokinetics analysis to the data on the in-blood drug concentration obtained from a clinical trial, the following two cases are conceivable: A case of estimating a pharmacokinetic parameter for each examinee (individual), and a case of estimating a pharmacokinetic parameter in a population out of which an individual is extracted. In the former case, a method is generally used where the measurements of the in-blood drug concentration are made for each individual at large number of points-in-time, and where the pharmacokinetic parameter is estimated for each individual on the basis of the data obtained. An early-stage clinical trial, such as clinical phase-I trial, is a typical example of this former case. The latter case is an analysis method commonly referred to as “population pharmacokinetic analysis”. In the population pharmacokinetic analysis, the pharmacokinetics in a population out of which each individual is extracted is analyzed based on the data on the individual.  
         [0006]     In general, the data on the pharmacokinetics exhibits large variations between individuals and inside an individual. Even if a constant amount of drug is dosed, the in-blood drug concentrations differ inside an individual or between individuals. Also, even if the dose is adjusted so that the in-blood drug concentrations will become the same, the effect and action of the drug will not necessarily become the same. In this way, as one of the causes by which the in-blood drug concentrations and the effect and action of the drug differ between individuals, it is conceivable that workings of genes related with the various processes described above differ depending on the individuals. As one of accomplishments of the Human Genome Project, it is conceivable that individual genetic information (i.e., the individual genetic polymorphism) is utilized in the field of Pharmacokinetics. Expectations are now getting more and more placed on this utilization of the individual genetic information as an effective method for making it possible to understand the problem of an individual difference in a drug therapy more scientifically than ever.  
         [0007]     Statistical genetics method is one of the effective methods for analyzing the relationship between the genetic polymorphism of an individual and occurrence of diseases and side effects of drugs. The statistical genetics method is a genetic name for methodologies of searching for the genetic polymorphism related with a specific trait (presence or absence of a disease, and presence or absence of side effect of a drug) by taking advantage of statistics. The statistical genetics method is the methodologies onto which attention has been focused in recent years.  
         [0008]     The statistical genetics method is classified into a linkage analysis method and an analysis method using linkage disequilibrium. The former is the analysis method using the linkage, i.e., an exception of Mendel&#39;s Independent Law. The latter is the analysis method using the linkage disequilibrium, i.e., a non-independent phenomenon of a plurality of genetic polymorphisms within a group. Family-line information is needed in the linkage analysis method; whereas the family-line information is not necessarily needed in the analysis method using the linkage disequilibrium. On account of this, in recent years, the analysis method using the linkage disequilibrium has become prevalent. In the analysis method using the linkage disequilibrium, frequencies of a genetic polymorphism on which attention is focused, or frequencies of a combination (haplotype) of a plurality of genetic polymorphisms in an area where the linkage disequilibrium exists are compared between groups of different traits (case group and control group, and responder group to a drug and non-responder group thereto). This comparison allows investigation of the relationship between a specific trait and the genetic polymorphisms or haplotype.  
         [0009]     In the analysis method using the linkage disequilibrium, in general, the analysis is made based on the haplotype or diplotype configuration (a combination of two haplotypes which an individual has). Up to the present, several proposals have been made concerning methodologies for estimating the haplotype frequency of a population and the diplotype configuration of an individual on the basis of the data on the genetic polymorphism of the individual. The representative proposals are a methodology using EM algorithm as described in Excoffier L &amp; Slatkin M: Maximum-likelihood estimation of molecular haplotype frequencies in a diploid population, Mol Biol Evol, Vol. 12, pp. 921-927, 1995, and PHASE method as described in Stephens M et al.: A new statistical method for haplotype reconstruction from population data, Am J Hum Genet, Vol. 68, pp. 978-989, 2001.  
         [0010]     As described above, it is conceivable that the genetic polymorphism of an individual is utilized in the field of Pharmacokinetics. Expectations are now getting more and more placed on this utilization of the genetic polymorphism as the effective method for making it possible to understand the problem of an individual difference in a drug therapy more scientifically than ever. However, a method for analyzing the pharmacokinetics in such a manner that the influence of the genetic polymorphism of an individual is taken into consideration has been not established yet.  
         [0011]     Also, the estimation methodologies reported up to the present for estimating the haplotype frequency of a population and the diplotype configuration of an individual require the genetic polymorphism data on a large number of individuals (a few hundreds to a few thousands of individuals in general) in order to enhance the estimation accuracy. Usually, however, the number of individuals obtainable in a clinical trial or the like is about a few tens of individuals. On account of this, acquiring the genetic polymorphism data on a large number of individuals is difficult in the pharmacokinetic analysis.  
         [0012]     Accordingly, there has existed a problem that the estimation accuracy is insufficient in estimating the haplotype frequency of a population and the diplotype configuration of an individual.  
       SUMMARY OF THE INVENTION  
       [0013]     In view of this situation, an object of the present invention is to solve the above-described problems, and to provide the method for analyzing the pharmacokinetics in such a manner that the influence of the genetic polymorphism of an individual is taken into consideration. Also, another object of the present invention is to provide a method for allowing implementation of high-accuracy haplotype frequency estimation and diplotype configuration estimation even in the case where the number of individuals is small from which the data is obtainable when making the pharmacokinetic analysis.  
         [0014]     The pharmacokinetic analysis system and method of the present invention estimate the diplotype configuration of an individual, and expresses a pharmacokinetic parameter as a function of the diplotype configuration of the individual. This processing makes it possible to configure pharmacokinetic models, and to make the pharmacokinetic analysis where the influence of the genetic polymorphism of the individual is taken into consideration. Also, the system and method of the present invention create the individual in a pseudo manner on the basis of information on the haplotype frequency and the number of individuals obtainable from documents, empirical knowledge, and the like. Then, using the data created, the system and method correct the estimation result of the diplotype configuration of the individual, thereby implementing an enhancement in the estimation accuracy. Moreover, the system and method select an optimum model from among the pharmacokinetic models, and displays a drug-concentration time course curve based on the selected model in a manner of being equipped with a confidence interval. Then, the system and method permit a user to set the confidence interval freely. This processing makes it possible to perform stratification of patients depending on a purpose.  
         [0015]     According to the pharmacokinetic analysis method of the present invention, it becomes possible to analyze the pharmacokinetics in such a manner that the influence of the genetic polymorphism of an individual is taken into consideration. Also, it becomes possible to investigate in detail the problem of an individual difference in a drug therapy.  
         [0016]     Other objects, features and advantages of the invention will become apparent from the following description of the embodiments of the invention taken in conjunction with the accompanying drawings. 
     
    
     BRIEF DESCRIPTION OF THE DRAWINGS  
       [0017]      FIG. 1  is a diagram for illustrating an example of the device configuration of a pharmacokinetic analysis device;  
         [0018]      FIG. 2  is a diagram for illustrating an example of the embodiment of a pharmacokinetic analysis program;  
         [0019]      FIG. 3  is a diagram for illustrating an example of the embodiment of the pharmacokinetic analysis program in a case of including a step of correcting the estimation result of the diplotype configuration;  
         [0020]      FIG. 4  is a diagram for illustrating an example of the screen of the pharmacokinetic analysis program;  
         [0021]      FIG. 5  is a diagram for illustrating an example of the setting screen of a population pharmacokinetic model building condition and a stratification condition;  
         [0022]      FIG. 6  is a diagram for illustrating an example of the setting screen of the haplotype frequency, sample size, and weighting coefficient;  
         [0023]      FIG. 7  is a diagram for illustrating an example of the data structure of genetic data;  
         [0024]      FIG. 8  is a diagram for illustrating an example of the data structure of the pharmacokinetic data and clinical data;  
         [0025]      FIG. 9  is a diagram for illustrating an example of the embodiment of a diplotype configuration estimation step;  
         [0026]      FIG. 10  is a diagram for illustrating an example of the embodiment of the diplotype configuration correction step;  
         [0027]      FIG. 11  is a diagram for illustrating effectiveness of the diplotype-configuration estimation result correction made by the diplotype configuration correction step;  
         [0028]      FIG. 12  is a diagram for illustrating an example of the embodiment of the pharmacokinetic model building step and the pharmacokinetic parameter estimation step;  
         [0029]      FIG. 13  is a diagram for illustrating a description example of a diagnosis-assisting knowledge database  15 ;  
         [0030]      FIG. 14  is a diagram for illustrating a 1-compartment model of oral administration;  
         [0031]      FIG. 15  is a diagram for illustrating an example of the device configuration of the pharmacokinetic analysis device in the case of including the step of correcting the estimation result of the diplotype configuration;  
         [0032]      FIG. 16  is a diagram for illustrating an example of the screen for selecting an optimum pharmacokinetic model; and  
         [0033]      FIG. 17  is a diagram for illustrating an example of the display screen of a stratification-condition satisfaction judgment result. 
     
    
     DESCRIPTION OF THE INVENTION  
       [0034]      FIG. 1  is a diagram for illustrating an example of the device configuration of a pharmacokinetic analysis device according to the present invention. The pharmacokinetic analysis device according to the present invention is mainly configured by an electronic computer such as the so-called personal computer. A processor  1 , a memory  2 , an input device  3 , a display  4 , and an external storage device  10  are connected to a system bus  5 . The external storage device  10  includes therein drug concentration data  11  on a plurality of individuals obtained from a clinical trial or the like, genetic data  12  on the plurality of individuals, clinical data  13  on the plurality of individuals, and a pharmacokinetic analysis program  14  for making the pharmacokinetic analysis according to the present invention. In addition to these, the device  10  of course includes a database and processing programs needed for implementing functions as the electronic computer.  
         [0035]     The pharmacokinetic analysis program  14  includes a data acquisition device  141  for acquiring the drug concentration data  11 , the genetic data  12 , and the clinical data  13  on a plurality of individuals, a diplotype configuration estimation device  142  for estimating diplotype configurations of the individuals on the basis of the genetic data  12  on the plurality of individuals acquired by the data acquisition device  141 , a pharmacokinetic model building device  143  for building pharmacokinetic models on the basis of the diplotype configurations of the individuals estimated by the diplotype configuration estimation device  142 , a pharmacokinetic parameter estimation device  144  for estimating a pharmacokinetic parameter on the basis of the data acquired by the data acquisition device  141  and the pharmacokinetic models built by the pharmacokinetic model building device  143 , an optimum model selection device  145  for calculating information criterion of the pharmacokinetic models on the basis of the pharmacokinetic models built by the pharmacokinetic model building device  143  and the pharmacokinetic parameter estimated by the pharmacokinetic parameter estimation device  144 , judging whether or not the pharmacokinetic models built by the pharmacokinetic model building device  143  are optimum on the basis of the information criterion calculated, and selecting an optimum pharmacokinetic model, a stratification device  146  for creating a drug-concentration time course curve on the basis of the optimum pharmacokinetic model selected by the optimum model selection device  145 , and performing stratification of the drug-concentration time course curve, and an analysis result output device  147  for outputting, as an analysis result, the data on the drug-concentration time course curve created by the stratification device  146 . Although the program like this has been described assuming that the program exists in the external storage device  10 , the program may also be installed into the present device via a medium. Hereinafter, referring to  FIG. 2  and  FIG. 3 , the explanation will be given below concerning the embodiment of the pharmacokinetic analysis program  14 .  
         [0036]      FIG. 2  is a diagram for illustrating an example of the embodiment of the pharmacokinetic analysis based on the pharmacokinetic analysis program  14 . Hereinafter, the explanation will be given below concerning each step in  FIG. 2 .  
         [0037]     A data acquisition step  210 : The data acquisition device  141  acquires the drug concentration data  11 , the genetic data  12 , and the clinical data  13  on a plurality of individuals. The data to be acquired are specified by the user via a setting screen which will be explained later using  FIG. 4 .  
         [0038]     A diplotype configuration estimation step  211 : The diplotype configuration estimation device  142  estimates diplotype configurations of the individuals on the basis of the genetic data  12  on the plurality of individuals acquired at the data acquisition step  210 . The details of the diplotype configuration estimation step  211  will be explained later using  FIG. 9 .  
         [0039]     A pharmacokinetic model building step  212 : The pharmacokinetic model building device  143  builds pharmacokinetic models on the basis of a population pharmacokinetic model building condition  21  and a stratification condition  22 , i.e., the condition for performing stratification of patients. The population pharmacokinetic model building condition  21  and the stratification condition  22  are specified by the user via a setting screen which will be explained later using  FIG. 5 .  
         [0040]     A pharmacokinetic parameter estimation step  213 : The pharmacokinetic parameter estimation device  144  estimates a pharmacokinetic parameter on the basis of the data acquired at the data acquisition step  210  and the pharmacokinetic models built at the pharmacokinetic model building step  212 .  
         [0041]     The details of the pharmacokinetic model building step  212  and the pharmacokinetic parameter estimation step  213  will be explained later using  FIG. 12 .  
         [0042]     An information criterion calculation step  214 : The optimum model selection device  145  calculates information criterion of the pharmacokinetic models on the basis of the pharmacokinetic models built at the pharmacokinetic model building step  212  and the pharmacokinetic parameter estimated at the pharmacokinetic parameter estimation step  213 .  
         [0043]     A model selection step  215 : The optimum model selection device  145  judges whether or not the pharmacokinetic models built at the pharmacokinetic model building step  212  are optimum on the basis of the information criterion calculated at the information criterion calculation step  214 , and the selection device  145  selects an optimum pharmacokinetic model.  
         [0044]     The details of the information criterion calculation step  214  and the model selection step  215  will be explained later using  FIG. 16 .  
         [0045]     A drug-concentration time course curve creation step  216 : The stratification device  146  creates a drug-concentration time course curve on the basis of the optimum pharmacokinetic model selected at the model selection step  215 .  
         [0046]     A stratification-condition satisfaction judgment step  217 : The stratification device  146  judges whether or not the stratification condition  22  is satisfied based on the drug-concentration time course curve created at the drug-concentration time course curve creation step  216 .  
         [0047]     An analysis continuation judgment step  218 : The stratification device  146  judges whether or not the pharmacokinetic analysis should be continued.  
         [0048]     The details of the drug-concentration time course curve creation step  216 , the stratification-condition satisfaction judgment step  217 , and the analysis continuation judgment step  218  will be explained later.  
         [0049]     An analysis result output step  219 : The analysis result output device  147  outputs, as an analysis result  23 , the data on the drug-concentration time course curve created at the drug-concentration time course curve creation step  216 .  
         [0050]      FIG. 3  is a diagram for illustrating an example of the embodiment of the pharmacokinetic analysis based on the pharmacokinetic analysis program  14  in a case of including a step of correcting the estimation result of the diplotype configurations using data created in a pseudo manner on the basis of given information on the haplotype frequency and the number of individuals (samples).  FIG. 15  illustrates an example of the device configuration of the pharmacokinetic analysis device in the present embodiment. The device configuration of the pharmacokinetic analysis device in the present embodiment turns out to become a one where a diplotype configuration correction device  151  for correcting the estimation result of the diplotype configurations using data created in a pseudo manner on the basis of given information on the haplotype frequency and the number of individuals (samples) is added to the device configuration illustrated in  FIG. 1 . As illustrated in  FIG. 3 , the pharmacokinetic analysis based on the pharmacokinetic analysis program  14  in the present embodiment includes each step illustrated in  FIG. 2  and a diplotype configuration correction step  310 . At the diplotype configuration correction step  310 , the diplotype configuration correction device  151  creates the data in a pseudo manner on the basis of haplotype frequency  31 , sample size  32 , and weighting coefficient  33  (indicating importance of documents or empirical knowledge) which are obtainable from documents, published database on the Web, empirical knowledge, and the like, and the correction device  151  corrects the estimation result of the diplotype configurations by using the data. The haplotype frequency  31 , the sample size  32 , and the weighting coefficient  33  are specified by the user via a setting screen which will be explained later using  FIG. 6 .  
         [0051]      FIG. 4  is a diagram for illustrating an example of the screen of the pharmacokinetic analysis program  14  of the present invention. A pharmacokinetic analysis program screen  41  in the present embodiment includes a data input unit  411 , a analysis condition setting button  412 , a diplotype configuration correction condition setting button  413 , and an analysis execution button  414 .  
         [0052]     In the data input unit  411 , the drug concentration data  11 , the genetic data  12 , and the clinical data  13  are inputted which are needed for the data acquisition step  210  in  FIG. 2  and  FIG. 3 . In the data input unit  411 , e.g., a pass or the like to a file or database is inputted where the drug concentration data  11 , the genetic data  12 , and the clinical data  13  are stored.  
         [0053]     Selecting and pressing the analysis condition setting button  412  makes it possible to set the population pharmacokinetic model building condition  21  and the stratification condition  22  needed for the pharmacokinetic model building step  212  in  FIG. 2  and  FIG. 3 . The setting screen for the population pharmacokinetic model building condition  21  and the stratification condition  22  will be explained later using  FIG. 5 .  
         [0054]     Selecting and pressing the diplotype configuration correction condition setting button  413  makes it possible to set the haplotype frequency  31 , the sample size  32 , and the weighting coefficient  33  needed for the diplotype configuration correction step  310  in  FIG. 3 . The setting screen for the haplotype frequency  31 , the sample size  32 , and the weighting coefficient  33  will be explained later using  FIG. 6 .  
         [0055]     Selecting and pressing the analysis execution button  414  causes the pharmacokinetic analysis program  14  to start the pharmacokinetic analysis.  
         [0056]      FIG. 5  is a diagram for illustrating an example of the setting screen of the population pharmacokinetic model building condition  21  and the stratification condition  22  needed for the pharmacokinetic model building step  212  in  FIG. 2  and  FIG. 3 . An analysis condition setting screen  51  in the present embodiment includes a population pharmacokinetic model building condition input unit  511 , a pharmacokinetic parameter setting unit  512 , a stratification condition input unit  513 , and a setting button  514 .  
         [0057]     In the population pharmacokinetic model building condition input unit  511 , for example, like the screen in the present embodiment, a compartment model to be used for the pharmacokinetic analysis is selected by a pull-down button or the like. For example, like the present embodiment, when a 1-compartment model with oral administration is selected, pharmacokinetic parameters related with the selected 1-compartment model with oral administration are displayed on the pharmacokinetic parameter setting unit  512  together with check buttons and the like. In the present embodiment, absorption rate K a , elimination rate K e , distribution volume V d , and in-blood drug concentration C are displayed as the pharmacokinetic parameters related with the 1-compartment model with oral administration. Of course, a pharmacokinetic parameter other than the pharmacokinetic parameters illustrated in  FIG. 5 , such as, e.g., clearance CL, may also be displayed. In the pharmacokinetic parameter setting unit  512 , which pharmacokinetic parameter of the pharmacokinetic parameters should be expressed as a function of the diplotype configuration can be set in a manner of being classified into fixed effect and random effect. In the present embodiment, the fixed effect of the elimination rate K e  is set as the function of the diplotype configuration. Also, in the pharmacokinetic parameter setting unit  512 , an error model can be set together with the pharmacokinetic parameter. For example, like the present embodiment, the error model to be used for the pharmacokinetic analysis can be selected from among an absolute error model, a relative error model, an exponential error model, and a mixed error model. In the present embodiment, the relative error model is selected as the error model. In the stratification condition input unit  513 , for example, like the screen in the present embodiment, a stratification condition for performing stratification of patients is selected by a pull-down button or the like. In the present embodiment, the stratification condition  22  is set assuming that the stratification of patients is performed when no overlapping exists in a 90-% confidence interval of the drug-concentration time course curve (i.e., no overlapping exists between lower-limit of the 90-% confidence interval of an upper-side drug-concentration time course curve and upper-limit of the 90-% confidence interval of a lower-side drug-concentration time course curve). As the stratification condition  22 , a stratification condition other than the one indicated in the present embodiment may also be set, such as, e.g., the stratification is performed when an x or more difference exists in the drug concentrations, the stratification is performed when a difference between lower-limit of the 95-% confidence interval of the upper-side drug-concentration time course curve and upper-limit of the 95-% confidence interval of the lower-side drug-concentration time course curve is y or more, the stratification is performed regarding patients whose drug concentrations attain to a critical area and do not attain thereto. The population pharmacokinetic model building condition is set by the setting button  514 .  
         [0058]      FIG. 6  is a diagram for illustrating an example of the setting screen of the haplotype frequency  31 , the sample size  32 , and the weighting coefficient  33  needed for the diplotype configuration correction step  310  in  FIG. 3 . A diplotype configuration correction condition setting screen  61  in the present embodiment includes a haplotype number input unit  611 , a haplotype frequency input unit  612 , a sample-size/weighting-coefficient input unit  613 , and a setting button  614 .  
         [0059]     In the haplotype number input unit  611 , for example, like the screen in the present embodiment, the number of haplotypes existing in a genetic area becoming the analysis target is inputted by a pull-down button or the like. For example, like the present embodiment, if the number of the haplotypes is inputted as being 3, a screen for inputting frequencies of the respective haplotypes (which are denoted by H 1  to H 3  each, for example) is displayed on the haplotype frequency input unit  612 . In the sample-size/weighting-coefficient input unit  613 , the sample size  32  and the weighting coefficient  33  are inputted, and the diplotype configuration correction condition is set by the setting button  614 . The user inputs the diplotype configuration correction condition on the basis of information obtainable from documents, published database on the Web, empirical knowledge, and the like.  
         [0060]      FIG. 7  is a diagram for illustrating an example of the data structure of the genetic data  12  of the present invention.  FIG. 7  illustrates an example of the data structure of genetic polymorphism information on each individual basis. In the example in  FIG. 7 , genetic type information on the genetic polymorphisms which are X in number (all of which are assumed to be SNP polymorphisms in the present embodiment) are stored. For example, an individual PERSON_ 1  has two A alleles for a polymorphism SNP_ 1  (i.e., is homo of the A alleles for the polymorphism SNP_ 1 ), and has two G alleles for a polymorphism SNP_ 2  (i.e., is homo of the G alleles for the polymorphism SNP_ 2 ). Similarly, an individual PERSON_ 2  has one A allele and one T allele for the polymorphism SNP_ 1  (i.e., is hetero for the polymorphism SNP_ 1 ), and has one C allele and one G allele for the polymorphism SNP_ 2  (i.e., is hetero for the polymorphism SNP_ 2 ). Although  FIG. 7  illustrates the example of the data on the SNP polymorphisms, data on another genetic polymorphism, such as ATRP polymorphism, may also be stored therein.  
         [0061]      FIG. 8  is a diagram for illustrating an example of the data structure of the drug pharmacokinetic data  11  and the clinical data  13 .  FIG. 8  illustrates an example where the dose, administration interval, elapsed time, and in-blood drug concentration on each individual basis are stored as the drug pharmacokinetic data, and where the weight on each individual basis is stored as the clinical data. Although  FIG. 7  illustrates the example where the drug pharmacokinetic data  11  and the clinical data  13  are stored into the same table, the drug pharmacokinetic data  11  and the clinical data  13  may also be stored in a separated manner.  
         [0062]     Next, referring to  FIG. 9 , the explanation will be given below concerning an example of the embodiment of the diplotype configuration estimation step  211  in  FIG. 2  and  FIG. 3 . In the present embodiment, the haplotype frequencies of a population and the diplotype configurations of individuals are estimated using the EM algorithm. Now, consider a sample group which includes n individuals. In this sample group, consider haplotypes in a plurality of marker locuses under linkage, and let frequencies of the haplotypes in the population be F=(F 1 , F 2 , . . . , F M ) Here, M is the total number of possible haplotypes. For example, if all of the marker locuses are SNP locuses, M=2L holds, letting the number of the locuses be L. Observed data on the genetic type of each individual in the plurality of marker locuses under the linkage are represented by G=(G 1 , G 2 , G n ). In many cases, G i  is incomplete data.  
         [0063]     Accordingly, the diplotype configurations corresponding to G i  are not determined into one configuration in many cases. In the case like this, posterior probability distributions (which will be referred to as “diplotype distributions”) on possible diplotype configurations are defined. The diplotype configuration corresponding to G i  for each individual i (i=1, 2, . . . , n) is represented by D ij  (j=1, 2, . . . , m i ). Here, m i  is the number of the possible diplotype configurations for G i , and maximum value of m i  is M. A storage area for storing the M possible haplotypes (which are denoted by H 1 , H 2 , . . . , H M  each) in the population and the frequencies F of the M possible haplotypes, and a storage area for storing the possible diplotype configuration D ij  for each individual and the posterior probability distributions of the possible diplotype configuration D ij  are ensured on the memory  2  on the basis of the genetic data  12  (which is equivalent to the above-described G) on the plurality of individuals acquired by the data acquisition device  141  in  FIG. 1  and  FIG. 15  at the data acquisition step  210  in  FIG. 2  and  FIG. 3 . At the diplotype configuration estimation step  211  in  FIG. 2  and  FIG. 3 , the diplotype configuration estimation device  142  in  FIG. 1  and  FIG. 15  estimates the diplotype configuration of each individual in association with the data stored on the memory  2  and in accordance with the following respective steps:  
         [0064]     An initial value setting step  91 : First, initial values F (0) of the haplotype frequencies are allocated to the M possible haplotypes (which are denoted by H 1 , H 2 , . . . , H M  each), then being stored into the storage area. The sum total of the haplotype frequencies is equal to 1. The initial value of each haplotype frequency is given at random such that the sum total of the haplotype frequencies becomes equal to 1. Otherwise, 1/M is equally given to each of the M possible haplotypes.  
         [0065]     Next, with respect to 0, 1, 2, . . . , F(t+1) is calculated from F(t) in accordance with the following diplotype distribution calculation step  92  to haplotype frequency update step  96 :  
         [0066]     The diplotype distribution calculation step  92 : Each diplotype configuration D ij  includes the two haplotypes H 1  and H m , where 1≦1≦M and 1≦m≦M. When the haplotype frequencies F(t) of the population are given, the probability with which D ij  will be obtained is given by the following Expression (1):  
             (     Expression   ⁢           ⁢   1     )                             Pr   ⁡     (     D   ij     )       =     {           F   l       (   t   )     2             l   =   m             if   ⁢           ⁢   l     =   m               2   ⁢     F   l     (   t   )       ⁢     F   m     (   t   )               l   ≠   m           if   ⁢           ⁢   l   ⁢   ⁢   m                     (   1   )             
 
         [0067]     Consequently, the posterior probability Pr(D ij |G i ) that the diplotype configuration of the individual i is D ij  under the observed data G i  on the genetic type is given by the following Expression (2) from Bayes&#39; theorem:  
             (     Expression   ⁢           ⁢   2     )                             Pr   ⁡     (       D   ij     |     G   i       )       =           Pr   ⁡     (     D   ij     )       ⁢     Pr   ⁡     (       G   i     |     D   ij       )             ∑     k   =   1       m   i       ⁢       Pr   ⁡     (     D   ik     )       ⁢     Pr   ⁡     (       G   i     |     D   ik       )             =       Pr   ⁡     (     D   ij     )           ∑     k   =   1       m   i       ⁢     Pr   ⁡     (     D   ik     )                     (   2   )             
 
         [0068]     Calculating this Expression for all j (j=1, 2, . . . , m i ) determines the diplotype distribution of the individual i. This diplotype distribution is applied to all the individuals in the sample group, then being stored into the storage area.  
         [0069]     A likelihood calculation step  93 : The determination of the diplotype distribution of the individual i makes it possible to calculate the entire likelihood in the sample group. The entire likelihood can be expressed by the following Expression (3) by connecting likelihoods of all the diplotype configurations for each individual and further, by connecting likelihoods of all the individuals. A storage area for storing the likelihood is ensured on the memory  2 , thereby storing the likelihood.  
             (     Expression   ⁢           ⁢   3     )                             L   ⁡     (     F     (   t   )       )       =       Pr   ⁡     (     G   |     F     (   t   )         )       =       ∏     i   =   1     n     ⁢           ⁢       ∑     j   =   1       m   i       ⁢     Pr   ⁡     (     D   ij     )                     (   3   )             
 
         [0070]     A haplotype-frequency expectation value calculation step  94 : The determination of the diplotype distribution of the individual i makes it possible to calculate expectation values of the haplotype frequencies of the population from the diplotype distributions of all the individuals in the sample group. The expectation values of the haplotype frequencies of the population are given by the following Expression (4):  
             (     Expression   ⁢           ⁢   4     )                             E   ⁡     [     F   i     (   t   )       ]       =       1     2   ⁢   n       ⁢       ∑     j   =   1     n     ⁢       ∑     k   =   1       m   i       ⁢       Pr   ⁡     (       D   jk     |     G   j       )       ⁢     N       D   jk     ⁢   i                       (   4   )             
 
         [0071]     Here, N Djki  is the number of H i  (i.e., any one of 0, 1, and 2) included in the diplotype configuration D jk .  
         [0072]     A convergence judgment step  95 : It is judged whether or not the value of the entire likelihood L(F) has been converged. If L(F(t+1))−L(F(t))&lt;β is satisfied, It is judged that the value has been converged, and the processing proceeds to the following diplotype distribution determination step  97 : Meanwhile, if L(F(t+1))−L(F(t))&lt;β is not satisfied, the processing returns to the diplotype distribution calculation step  92  via the following haplotype frequency update step  96 , then repeating the steps up to the convergence judgment step  95 : Here, β is a threshold value. A sufficiently small value such as, e.g., 10 −5  is set as the value of β. The haplotype frequency update step  96 : F is updated by setting F(t+1)=E[F(t)].  
         [0073]     The diplotype distribution determination step  97 : The value E[F]=F(EM) at the point-in-time when the entire likelihood has been converged is defined as maximum likelihood estimation value of the haplotype frequencies in the population. Moreover, the posterior probability Pr(D|G) at this time is defined as the diplotype distribution of the individual under the maximum likelihood estimation value of the haplotype frequencies in the population.  
         [0074]     Next, referring to  FIG. 10 , the explanation will be given below concerning an example of the embodiment of the diplotype configuration correction step  310  in  FIG. 3 . In the present embodiment, the haplotype frequencies of a population and the diplotype configurations of individuals are estimated using the EM algorithm. In the present embodiment as well, consider the sample group which includes n individuals. Then, the haplotype frequencies F in the population, the observed data G on the genetic type in a plurality of marker locuses under linkage, and the diplotype configuration D ij  (j=1, 2, . . . , m i ) corresponding to G i  for each individual i (i=1, 2, . . . , n) are set exactly as described earlier. In the present embodiment as well, the storage area for storing the M possible haplotypes (which are denoted by H 1 , H 2 , . . . , H m  each) in the population and the frequencies F of the M possible haplotypes, and the storage area for storing the possible diplotype configuration D ij  for each individual and the posterior probability distributions of the possible diplotype configuration D ij  are ensured on the memory  2  on the basis of the genetic data  12  (which is equivalent to the above-described G) on the plurality of individuals acquired by the data acquisition device  141  in  FIG. 15  at the data acquisition step  210  in  FIG. 3 .  
         [0075]     A diplotype configuration correction condition acquisition step  101 : The data on the haplotype frequency  31 , the sample size  32 , and the weighting coefficient  33  (indicating importance of documents or empirical knowledge) are acquired which are obtainable from documents, published database on the Web, empirical knowledge, and the like. The haplotype frequency  31 , the sample size  32 , and the weighting coefficient  33  are specified by the user via the setting screen explained earlier using  FIG. 6 .  
         [0076]     A replication data creation step  102 : Based on the data acquired at the diplotype configuration correction condition acquisition step  101 , samples (individuals) which are RN in number are created in a pseudo manner so that expectation value of the diplotype configuration frequency including two haplotypes H i  and H j  (whose frequencies are set at F(H i ) and F(H j ) each) becomes equal to F(H i ) 2  if i=j, and 2F(H i )F(H j ) if i         j. Here, R denotes the weighting coefficient  33 , and N denotes the sample size  32 .  
         [0077]     In the present embodiment, the data on the individuals acquired in advance (observed actually) are referred to as “original data”; while the data on the individuals created in a pseudo manner at the replication data creation step  102  are referred to as “replication data”.  
         [0078]     A data merge step  103 : The original data and the replication data are merged with each other. The genetic data  12  on the plurality of individuals acquired by the data acquisition device  141  and the replication data created in a pseudo manner at the replication data creation step  102  are added, thereby being formed into one piece of data.  
         [0079]     A step  104  to a step  1010  are steps corresponding to the initial value setting step  91  to the diplotype distribution determination step  97  in  FIG. 9 . However, at the diplotype distribution calculation step  92  in  FIG. 9 , the diplotype distributions of all the individuals are calculated. In contrast thereto, the diplotype distribution calculation step  105  in  FIG. 10  differs therefrom in the point that the diplotype distributions of the individuals of only the original data are calculated.  
         [0080]      FIG. 11  is a diagram for illustrating effectiveness of the diplotype-configuration estimation result correction made by the diplotype configuration correction step  310 . The longitudinal axis of the graph denotes accuracy of the diplotype configuration estimation result, and the transverse axis denotes the number of the samples (individuals). The accuracy is calculated by making a comparison between diplotype configuration data (correct-solution data) on individuals whose diplotype configurations are determined in advance, and diplotype configuration data created by dividing the correct-solution data into genetic type data in a plurality of marker locuses configuring the diplotype configuration, and re-estimating the diplotype configuration data from the genetic type data divided. The present invention (i.e., the proposed method) includes the diplotype configuration correction step  310 . This makes it possible to obtain the high accuracy as compared with the conventional method even in the case where the sample size is small.  
         [0081]     Next, referring to  FIG. 12 , the explanation will be given below concerning an example of the embodiment of the pharmacokinetic model building step  212  and the pharmacokinetic parameter estimation step  213  in  FIG. 2  and  FIG. 3 . The embodiment illustrated in  FIG. 12  includes a compartment model setting step  121 , a pharmacokinetic parameter setting step  122 , a likelihood calculation step  123 , and a parameter estimation step  124 . The compartment model setting step  121  and the pharmacokinetic parameter setting step  122  correspond to the pharmacokinetic model building step  212  in  FIG. 2  and  FIG. 3 . The likelihood calculation step  123  and the parameter estimation step  124  correspond to the pharmacokinetic parameter estimation step  213  in  FIG. 2  and  FIG. 3 .  
         [0082]     At the pharmacokinetic model building step  212 , based on the population pharmacokinetic model building condition  21  and the stratification condition  22  specified by the user via the setting screen explained earlier using  FIG. 5 , the pharmacokinetic model building device  143  sets a compartment model at the compartment model setting step  121 , and sets a pharmacokinetic parameter at the pharmacokinetic parameter setting step  122 . In the present embodiment, the pharmacokinetic analysis is made using the 1-compartment model with oral administration as illustrated in  FIG. 14 . In the 1-compartment model with oral administration as illustrated in  FIG. 14 , the estimation value C of the in-blood drug concentration at the point-in-time when a time t has elapsed after the administration is represented by the following Expression (5):  
             (     Expression   ⁢           ⁢   5     )                           C   =           k   a     ⁢   D         V   d     ⁡     (       k   a     -     k   el       )         ⁢     {       exp   ⁡     (       -     k   el       ⁢   t     )       -     exp   ⁡     (       -     k   a       ⁢   t     )         }               (   5   )             
 
         [0083]     Here, D is the dose, K a  is absorption rate, K el  is elimination rate, and V d  is distribution volume.  
         [0084]     The population pharmacokinetic analysis using the 1-compartment model with oral administration can be explained as follows: Now, assume that the in-blood drug concentration is measured for an individual i at m i  points-in-time t 11 , . . . , t imi , and that the resultant measurement values are x i1 , . . . , x imi . Expressing all of pharmacokinetic parameters of the individual i as a vector θ i  in batch, the measurement values x ij  of the in-blood drug concentration can be expressed as x ij =C (θ i , t ij )+ε ij . This is because the measurement values x ij  can be obtained by adding an error ε ij , i.e., a probability variable, to the estimation value C (θ 1 , t ij ) based on the 1-compartment model with oral administration. In the case of the 1-compartment model with oral administration, the pharmacokinetic parameter vector of the individual i is θ i =(k a   (i) , k el   (i) , V d   (i) ) T . Accordingly, letting the dose of the individual i be D i , the estimation value of the in-blood drug concentration at the point-in-time when a time t ij  has elapsed after the administration is represented by the following Expression (6):  
             (     Expression   ⁢           ⁢   6     )                             C   ⁡     (       θ   i     ,     t   ij       )       =           k   a     (   i   )       ⁢     D   i           V   d     (   i   )       ⁡     (       k   a     (   i   )       -     k   el     (   i   )         )         ⁢     {       exp   ⁡     (       -     k   el     (   i   )         ⁢     t   ij       )       -     exp   ⁡     (       -     k   a     (   i   )         ⁢     t   ij       )         }               (   6   )             
 
         [0085]     For simplicity, assuming that the error ε ij  follows a normal distribution of average 0 and variance σ ε   2 , probability density function representing distribution of the measurement values x i1 , . . . , x imi  at the m i  measurement points-in-time t i1 , . . . , t imi  can be represented by the following Expression (7):  
             (     Expression   ⁢           ⁢   7     )                             f   ⁡     (       x     i   ⁢           ⁢   1       ,   ⋯   ⁢           ,       x     im   i       |     k   a     (   i   )         ,     k   el     (   i   )       ,     V   d     (   i   )       ,     σ   ɛ   2       )       =         ∏     i   =   1       m   i       ⁢           ⁢       1       2   ⁢     πσ   ɛ   2           ⁢   exp   ⁢     {     -         (       x   ij     -     C   ⁡     (       θ   i     ,     t   ij       )         )     2       2   ⁢     σ   ɛ   2           }         =       ∏     j   =   1       m   i       ⁢           ⁢       1       2   ⁢     πσ   ɛ   2           ⁢   exp   ⁢     {     -         {       x   ij     -           k   a     (   i   )       ⁢     D   i           V   d     (   i   )       ⁡     (       k   a     (   i   )       -     k   el     (   i   )         )         ⁢     {       exp   ⁡     (       -     k   el     (   i   )         ⁢     t   ij       )       -     exp   ⁡     (       -     k   a     (   i   )         ⁢     t   ij       )         }         }     2       2   ⁢     σ   ɛ   2           }                   (   7   )             
 
         [0086]     For simplicity, the following normal distribution is assumed as a distribution for representing differences in the pharmacokinetic parameters between the individuals: k a   (i) ˜N(μ ka , σ ka   2 ), k el   (i) ˜N(μ kel , σ kel   2 ), V d   (i) ˜N(μ Vd , σ Vd   2 ).  
         [0087]     Although population parameters are unknown, likelihood function for the individual i is represented by the following Expression (8) under the situation that the measurement data x i1 , . . . , x imi  on the individual are given: This representation is made possible by considering the above-described probability density function as a function of the unknown population parameters. What is referred to as “the unknown population parameters” here means the variance σ ε   2  of the error, and μ ka , σ ka   2 , μ kel , σ kel   2 , μ Vd , and σ Vd   2  which specify the distribution of the pharmacokinetic parameters in the population.  
             (     Expression   ⁢           ⁢   8     )                               L   1     ⁡     (       μ     k   a       ,     σ     k   a     2     ,     μ     k   el       ,     σ     k   el     2     ,     μ     V   d       ,     σ     V   d     2     ,     σ   ɛ   2       )       =     ∫     ∫     ∫       [       ∏     j   =   1       m   i       ⁢           ⁢       1       2   ⁢     πσ   ɛ   2           ⁢   exp   ⁢     {     -         {       x   ij     -           k   a     (   i   )       ⁢     D   i           V   d     (   i   )       ⁡     (       k   a     (   i   )       -     k   el     (   i   )         )         ⁢     {       exp   ⁡     (       -     k   el     (   i   )         ⁢     t   ij       )       -     exp   ⁡     (       -     k   a     (   i   )         ⁢     t   ij       )         }         }     2       (     2   ⁢     σ   ɛ   2       )         }     ×     1       2   ⁢     πσ     k   a     2           ⁢   exp   ⁢     {     -         (       k   a     (   i   )       -     μ     k   a         )     2       2   ⁢     σ     k   a     2           }     ×     1       2   ⁢     πσ     k   el     2           ⁢   exp   ⁢     {     -         (       k   el     (   i   )       -     μ     k   el         )     2       2   ⁢     σ     k   el     2           }     ×     1       2   ⁢     πσ     V   d     2           ⁢   exp   ⁢     {     -         (       V   d     (   i   )       -     μ     V   d         )     2       2   ⁢     σ     V   d     2           }         ]     ⁢     ⅆ     k   a     (   i   )         ⁢     ⅆ     k   el     (   i   )         ⁢     ⅆ     V   d     (   i   )                         (   8   )             
 
         [0088]     If the measurement values at the m i  points-in-time t i1 , . . . , t imi  for each individual i have been found to be x i1 , . . . , x imi , logarithmic likelihood function for all the data can be represented by the following Expression (9), letting the number of all the individuals be n:  
             (     Expression   ⁢           ⁢   9     )                             l   ⁡     (       μ     k   a       ,     σ     k   a     2     ,     μ     k   el       ,     σ     k   el     2     ,     μ     V   d       ,     σ     V   d     2     ,     σ   ɛ   2       )       =     log   ⁢       ∏     i   =   1     n     ⁢           ⁢       L   i     ⁡     (       μ     k   a       ,     σ     k   a     2     ,     μ     k   el       ,     σ     k   el     2     ,     μ     V   d       ,     σ     V   d     2     ,     σ   ɛ   2       )                   (   9   )             
 
         [0089]     In the population pharmacokinetic analysis, the population parameters μ ka , σ ka   2 , μ kel , σ kel   2 , μ Vd , σ Vd   2 , and σ ε   2  which will maximize this logarithmic likelihood function are estimated. The analysis model to be used in the population pharmacokinetic analysis turns out to become a mixed effect model, which includes the fixed effect of the average values μ ka  and μ kel  of the population parameters and μ Vd , and the random effect of σ ka   2 , σ kel   2 , σ Vd   2 , and σ ε   2  whose average values are equal to zero.  
         [0090]     At the pharmacokinetic model building step  212 , the pharmacokinetic model is built by the following two steps:  
         [0091]     The compartment model setting step  121 : In accordance with the population pharmacokinetic model building condition  21  specified by the user via the setting screen explained earlier using  FIG. 5 , the 1-compartment model with oral administration is set based on the Expression (5).  
         [0092]     The pharmacokinetic parameter setting step  122 : In accordance with the population pharmacokinetic model building condition  21  specified by the user via the setting screen explained earlier using  FIG. 5 , the pharmacokinetic parameters are set. In the present embodiment, the population parameters μ ka , μ kel , and μ Vd  are expressed as functions of the diplotype configuration of each individual, and are considered to be the fixed effect of the mixed effect model, thereby configuring the population pharmacokinetic analysis model. For simplicity, in the present embodiment, it is assumed that the elimination of the drug is due to metabolism enzyme alone, and that the value of μ kel  depends on the diplotype configuration of each individual. The fixed effect and the random effect of the mixed effect model are set as the following Expression (10) and Expression (11), respectively: (Expression 10)  
             (     Expression   ⁢           ⁢   10     )                             μ     k   el       =       θ   1     D   1       ⁢     θ   2     D   2       ⁢     θ   3     D   3       ⁢     θ   4     D   4       ⁢     {               D   1     =   1     ,             D   2     =       D   3     =       D   4     =   0               (       H   1     /     H   1       )                   D   2     =   1     ,             D   1     =       D   3     =       D   4     =   0               (         H   1     /     H   2       ,       H   1     /     H   3         )                   D   3     =   1     ,             D   1     =       D   2     =       D   4     =   0               (         H   2     /     H   2       ,       H   2     /     H   3         )                   D   4     =   1     ,             D   1     =       D   2     =       D   3     =   0               (       H   3     /     H   3       )                       (   10   )               (     Expression   ⁢           ⁢   11     )                                 k   a     =       μ     k   a       ⁡     (     1   +     η     k   a         )         ,       η     k   a       ~     N   ⁡     (     0   ,     σ     η   ka     2       )           ⁢     
     ⁢         k   el     =       μ     k   el       ⁡     (     1   +     η     k   el         )         ,       η     k   el       ~     N   ⁡     (     0   ,     σ     η   el     2       )           ⁢     
     ⁢         V   d     =       μ   d     ⁡     (     1   +     η     V   d         )         ,       η     V   d       ~     N   ⁡     (     0   ,     σ     η     V   d       2       )           ⁢     
     ⁢       C   =       μ   C     ⁡     (     1   +   ɛ     )         ,     ɛ   ~     N   ⁡     (     0   ,     σ   ɛ   2       )                   (   11   )             
 
         [0093]     Here, H i  (i=1, 2, 3) denote the haplotypes, and H i /H j  (i=1, 2, 3, j=1, 2, 3) denote the diplotype configurations including the haplotypes H i  and H j . In the present embodiment, the three haplotypes are considered, then setting the fixed effect of μ kel  on each diplotype-configuration basis.  
         [0094]     At the pharmacokinetic parameter estimation step  213 , the pharmacokinetic parameter estimation device  144  calculates the likelihood at the likelihood calculation step  123 , and estimates the pharmacokinetic parameters at the parameter estimation step  124 . Hereinafter, the explanation will be given below regarding the likelihood calculation step  123  and the parameter estimation step  124 .  
         [0095]     The likelihood calculation step  123 : A storage area for storing the likelihood of each individual i to be calculated by the Expression (8) is ensured on the memory  2 . This calculation is made on the basis of the drug concentration data  11  and the clinical data  13  (which are equivalent to the above-described D, (t i1 , . . . , t imi ), and (x i1 , . . . , x imi )) on the plurality of individuals acquired by the data acquisition device  141  in  FIG. 1  and  FIG. 15  at the data acquisition step  210  in  FIG. 2  and  FIG. 3 , and the diplotype configuration data (which are equivalent to the above-described H i  and H j ) estimated by the diplotype configuration estimation device  142  in  FIG. 1  and  FIG. 15  at the diplotype configuration estimation step  211  in  FIG. 2  and  FIG. 3 . Next, the likelihood of each individual i is calculated by the Expression (8), then being stored into the storage area. Moreover, a storage area for storing the logarithmic likelihood for all the data to be calculated by the Expression (9) is ensured on the memory  2 . Furthermore, the logarithmic likelihood for all the data is calculated based on the Expression (9), then being stored into the storage area.  
         [0096]     The parameter estimation step  124 : The parameters which will maximize the logarithmic likelihood calculated at the likelihood calculation step  123  are estimated. Then, the parameters estimated are employed as the pharmacokinetic parameters in the population. The methods such as maximum likelihood method and EM algorithm are used for the parameter estimation.  
         [0097]     Next, the explanation will be given below concerning an example of the embodiment of the information criterion calculation step  214  in  FIG. 2  and  FIG. 3 . Although, in the present embodiment, Akaike Information Criterion (: AIC) is used as the information criterion, some other information criterion may also be used. Using maximum logarithmic likelihood of a model and free parameter number of the model, AIC can be calculated by the following Expression: −2* (maximum logarithmic likelihood of a model)+2*(free parameter number of the model). At the information criterion calculation step  214 , the optimum model selection device  145  calculates AIC on the basis of the pharmacokinetic models built at the pharmacokinetic model building step  212  and the pharmacokinetic parameters estimated at the pharmacokinetic parameter estimation step  213 . Using the above-described Expression, AIC is calculated from the maximum logarithmic likelihood and the pharmacokinetic parameter number calculated at the likelihood calculation step  123  and the parameter estimation step  124  in  FIG. 12 .  
         [0098]      FIG. 16  is a diagram for illustrating an example of the screen for selecting an optimum pharmacokinetic model by the optimum model selection device  145  at the model selection step  215 . An optimum model selection screen  161  illustrated at the upper portion in  FIG. 16  includes an information criterion display unit  1611  and a stratification execution button  1612 .  
         [0099]     The information criterion display unit  1611  displays calculation results of AIC of the respective models built by the pharmacokinetic model building device  143  at the pharmacokinetic model building step  212 . Using a sort button  16112  makes it possible to sort the calculation results in an ascendant or descendent order.  
         [0100]     Pressing a model detail button  16113  displays a model detail display screen  162  illustrated at the lower portion in  FIG. 16 . This screen makes it possible to know details of a model selected. The model detail display screen  162  includes a model detail display unit  1621 .  
         [0101]     The user selects the optimum model on the basis of the calculation results of AIC, then checking a selection model check box  16111  and pressing the stratification execution button  1612 . This starts the stratification by the stratification device  146 . In the present embodiment, the example has been indicated where the user selects the optimum model using the optimum model selection screen  161 . Instead of performing the interactive processing with the user, however, the optimum model selection device  145  may automatically select the model which will make AIC the minimum, and start the stratification by the stratification device  146 .  
         [0102]     Next, the explanation will be given below concerning the drug-concentration time course curve creation step  216 , the stratification-condition satisfaction judgment step  217 , and the analysis continuation judgment step  218  in  FIG. 2  and  FIG. 3 .  
         [0103]     At the drug-concentration time course curve creation step  216 , the drug-concentration time course curve is created on the basis of the optimum pharmacokinetic model selected by the optimum model selection device  145  at the model selection step  215  in  FIG. 2  and  FIG. 3 . Moreover, at the stratification-condition satisfaction judgment step  217 , it is judged whether or not the drug-concentration time course curve created satisfies the stratification condition  22  specified by the user via the setting screen explained earlier using  FIG. 5 .  
         [0104]     At the stratification-condition satisfaction judgment step  217 , the stratification device  146  judges whether or not the stratification condition  22  is satisfied based on the drug-concentration time course curve created at the drug-concentration time course curve creation step  216 . As a result of this, if the stratification condition  22  is satisfied, the processing proceeds to the analysis result output step  219 , where the analysis result output device  147  outputs the analysis result  23 . Also, at the stratification-condition satisfaction judgment step  217 , the stratification device  146  judges whether or not the stratification condition  22  is satisfied based on the drug-concentration time course curve created at the drug-concentration time course curve creation step  216 . As a result of this, if the stratification condition  22  is not satisfied, as the analysis continuation judgment step  218 , a stratification-condition satisfaction judgment result screen  171  as is illustrated in, e.g.,  FIG. 17  is displayed. The stratification-condition satisfaction judgment result screen  171  includes a stratification-condition satisfaction judgment result display unit  1711 . Pressing “Yes” button returns the processing to the pharmacokinetic model building step  212 , thereby continuing the pharmacokinetic analysis. Pressing “No” button terminates the pharmacokinetic analysis.  
         [0105]      FIG. 13  is a diagram for illustrating an example of the screen of the analysis result  23  outputted by the analysis result output step  219  in  FIG. 2  and  FIG. 3 . A pharmacokinetic estimation result output screen  131  in the present embodiment includes a drug-concentration time course curve display unit  1311  and a confidence interval setting bar  1312 .  
         [0106]     The drug-concentration time course curve display unit  1311  displays the data on the drug-concentration time course curve created at the drug-concentration time course curve creation step  216  in  FIG. 2  and  FIG. 3 . In the example of the screen in the present embodiment, graphs of drug-concentration time course curves for, e.g., four groups are displayed. Here, the four groups are an individual having the diplotype configuration of H 1 /H 1 , an individual having the diplotype configuration of H 1 /H 2  or H 1 /H 3 , an individual having the diplotype configuration of H 2 /H 2  or H 2 /H 3 , and an individual having the diplotype configuration of H 3 /H 3 . The longitudinal axis of the graph denotes in-blood drug concentration, and the transverse axis denotes elapsed time in stationary state. Each drug-concentration time course curve is displayed with the confidence interval.  
         [0107]     The confidence interval setting bar  1312  sets the confidence interval of each drug-concentration time course curve displayed on the drug-concentration time course curve display unit  1311 . The upper-limit of the confidence interval is 100%, and the lower-limit thereof is 0%. In the example of the screen in the present embodiment, the confidence interval-setting bar  1312  is set at, e.g., 90%. Accordingly, each drug-concentration time course curve (bold line) is displayed on the drug-concentration time course curve display unit  1311  together with the upper-limit and lower-limit (narrow lines) of the 90-% confidence interval.  
         [0108]     The user is permitted to freely set the confidence interval of each drug-concentration time course curve displayed as the analysis result. This feature allows the user to obtain information for performing the stratification of patients depending on a purpose.  
         [0109]     It should be further understood by those skilled in the art that although the foregoing description has been made on embodiments of the invention, the invention is not limited thereto and various changes and modifications may be made without departing from the spirit of the invention and the scope of the appended claims.