hubungan antara motivasi dan pengetahuan orang tua dengan ...

4TH International Conference on Mathematics and Statistics (ICoMS 2009) Universitas Malahayati Bandarlampung © MSMSSEA and Universitas Malahayati 2009

An Adaptive Grid Approach for modeling non-parametric data

Ibnu Gunawan1 , Mahendrawathi ER2, Nur Iriawan3 1

Post Graduate Student , Institute Teknologi Sepuluh Nopember, Indonesia 2 Department of Information System, Faculty of Information Technology, Institute Teknologi Sepuluh Nopember, Indonesia 3 Department of Statistic, Faculty of Mathematic, Institute Teknologi Sepuluh Nopember, Indonesia ____________________________________________________________________________

Abstract. Practically, data can’t always paternally be modeled as parametric form. It is due to number of data is limited and the relationship among data are frequently not known. Such patient medical data record modeling is an example of this case. Statistics non parametric, such as Non-parametric Adaptive Grid approach (NPAG), overcomes this kind of modeling and data analysis. NPAG can effectively and efficiently describe these non parametric patient medical data problem by using an adaptive grid approach. NPAG employs non parametric method couple with state-of-the-art primal-dual interior-point to maximize likelihood on each grid. This combination makes NPAG more efficient and effective and requires much less memory than the last generation of Non-parametric ExpectationMaximization (NPEM), and therefore faster computation feasible on PC. NPAG demonstrates the result as a mixture of density which lay on grids. This paper will describe NPAG algorithm and applying it to model the kidney surgery data using MATLAB. Keywords: Primal-dual interior point, Non-parametric Adaptive Grid, Likelihood

_____________________________________________________________________________ 1. Introduction Population pharmacokinetics (PPK) tries to find a model, which describes a drug course in the patients body and investigates statistical distribution of parameters of that model in the population under study. The purpose of these analysis may be purely utilitarian or cognitive (as one may wish to do with the data gathered during clinical observation). The statistical theory makes parameter estimation possible even in the case of sparse data. As an example of this case is when the number of measured drug concentrations from a given subject is less than the number of pharmacokinetic (PK) model parameters. This is a great advantage, because it makes routine clinical data usable for population analysis. However, an estimation of the model for sparse data is very difficult statistically and need a computational intensive task. A number of approaches to this problem have been investigated. Table I shows the most important of them along with the appropriate computer programs. Each of these methods constitutes significantly distinct approach to the PPK problem, both in the aspect of the assumptions as well as numeric algorithms applied. The first three approaches have their origins in the maximum likelihood method, the very fundamental statistical principle. The NONMEM programs searches for a likelihood function maximum by means of numerical algorithm. Next two methods search for likelihood maximizing parameters without necessity of explicit calculation of likelihood function. These methods differ in the way they express the parameter distribution. In the parametric EM, the parametric family of distributions (like normal or log – normal) is assumed. On the other hand, the non – parametric EM (NPEM) replaces continuous statistical distribution by their discrete approximations. This allow for a greater flexibility in a description of parameter distributions NPEM Algorithm is Nonparametric Adaptive Grid (NPAG)’s successor, used to get prior distribution from non parametric data [1]. NPAG difference from NPEM on second phase or known as maximation phase. NPAG use a primal

Wamiliana and Caccetta/4th International Conference on Mathematics and Statistics 2009

2

dual interior point method[2] while NPEM use Gaussian quadrature [3]. Because of that, we will explain NPEM as far as likelihood computation that exactly match with NPAG algorithm. Table 1. Methods and software for population pharmacokinetic analysis of sparse data. [2]

Method Nonlinear mixed-effect modelling Expectation maximization (EM) Nonparametric EM Markov Chain Monte Carlo (MCMC)

Program NONMEM (1) Thermo Kinetica (2) NPAG (3) WinBugs/PKBugs (4)

Authors/originators SheinerL., Beal S. (UCSF) MentrÈ F., Gomeni R. Jellife R., Schumitzky A., Leary R. (USC) Lunn D.J., St Maryís Hospital, London

2. NPEM The probability distribution function from a pair of data {(X1 ,Y1 ),(X2 ,Y2 ),...,(Xm ,Yn )}

,

where X i and Yi are independent, is as follow: n

f((x1 ,y1 ),(x 2 ,y 2 ),...,(x m ,y n )|θ)= f(yi |θ)f(x i |θ)

(1)

i=1

with θ consist of all distribution parameter from X and Y. Equation (1) known as complete data likelihood and ( (x1 ,y1 ),(x 2 ,y2 ),...,(x m ,yn ) ) known as complete data. When data is sparse, then parameter estimation process would be more difficult. For example, if x1 is missing, y1 can be eliminated so it will remain (n-1) pair of data. Elimination of y1 will affect y1 information missing also, eventhough y1 information can increase accuracy of estimation process. Sample x1 missing data is [4]: 

 f((x ,y ),(x ,y ),...,(x

x1 0

1

1

2

2

m

,yn )|θ)

(2)

Likelihood (2) known as incomplete data likelihood that must be maximized. If Y=(Y1 ,Y2 ,...,Ym ) incomplete data and X=(X1 ,X2 ,...,Xm ) is augmented data, then we can set it together in pair Y,Xas complete data. Therefore the relation between

g(.|θ) as a pdf of Y and the density of f(.|θ) as a pdf of a joint Y,Xthen:

g(y|θ)=  f(y,x|θ)dx

(3)

So then, L(θ|y)=g(y|θ) called uncomplete likelihood data and L(θ|y,x)=f(y,x|θ) is called as a complete likelihood data. Calculation of L(θ|y) is hard to do because x is missing data, while we must do integral for all value of x. In contrast, calculation of L(θ|y,x) is easier to do. Observation on n individu, E1 , E 2 ,..., E n, which each observation represented by ran-

dom vector  Yi ,θi  . Pair of random vector  Yi ,θi  is analogue with complete data on NPEM. For example Yi is a random vector for gentamicyn concentrate in i-th individu. The value of Yi =(i=1,2,...,n) can be measure, it is independent but non identical. Yi is in the Euclidean with n dimension, symbolized by E n . Random vector θi (=1,2,...,n) is a parametric vector for i individual. θ i is independently identically, so that θi =θ . θ value is unknown and therefore it is assumed as an un-observed variabel.

Wamiliana/4th International Conference on Mathematics and Statistics 2009

3

Yi |θ distribution is p( Yi |θ ) that can be assumed as a normal distribution, so the likelihood [2] is: LIJ =Ce-wss/2 (4) with K

WSS   ( yobs (tk )  y pred (tk ))2 /  k2

(5)

K 1

and K

C=1/(2πσ 2k )1/2

(6)

k=1

Therefore, equation (5) can be written as: K

K

LIJ =1/(2πσ ) e 2 1/2 k

-(

 (Yobs (t k )-Ypred (t k ))2 /σk2 )/2 k=1

(7)

k=1

or K

LIJ = k=1

 1 K  (Y (t )-Y (t ) 2  obs k pred k      2 k=1  k   



1 2πσ 2k

e

(8)

that’s likelihood calculation process in NPEM as the same as in NPAG method. 3. NPAG In NPAG, equation (8) must be change to equation (9) bellow. Then, this equation will

perform and compatible with each grid.

M     LIJ PJ  I=1  J=1  N

Likelihood =

(9)

to maximize equation (9) , it would be simpler when we convert to log function, so equation (9) will become:

M  ln    LIJ PJ  I=1  J=1  N

Log Likelihood =

(10)

subject to constraint: N

 P =1 J

J=1

PJ  0 to maximized equation (9) or to get result from equation (10), we can use primal dual interior point method that substitute gaussian quadrature in NPEM. [2] Primal-dual InteriorPoint used in NPAG is Newton Raphson method. Newton – Raphson is proposed here to


4

minimize the steps to reach a faster convergence of the likelihood computation. Equation (11) is Newton – Rapshon basic equation for f(x) with the initial value x 0 :

f ' (x n ) x n+1 =x n - '' f (x n )

(11)

Where f ' (x n ) is the first derivative function n in data x n and f '' (x n ) is the second order derivative function n on data x n . Because of that, we must find the first and second derivative from equation (8). To make it simple we can use a natural logarithmic function, by multiply equation (8) with natural logarithmic (ln) like this: 2

n 1 J  yobs (t k )-ypred (t k )  lnLIJ =- ln(2π)-    -nlnσ X 2 2 I=1  σk 

(12)

then we substitute yobs (t k )-ypred (t k ) with xk  k so equation (12) will become 2

n 1 J  x -μ  lnLIJ =- ln(2π)-   k k  -nlnσ X 2 2 I=1  σ k 

(13)

We must find first and second derivative from equation (13), namely L  ln Lij , with respect to  is as follow

L 1  2  

 x   

(14)

and

 1    2  x      L    2    2

2 L n  2 2   

(15)

while derivative with respect to  2 would be:

L n 1 2    3  (x   ) 2   

(16)


5

and second order derivative from  2 is





 

1 3  2 L  n 2      (x   )  2  

2 L n  x    2 3 2   4



2

(17)

So then we can substitute equation (14) and (15) to (11):

  2 1 n1  n   2  n 

n1  n  



 x    

x   

(18)

n

The same as  , by substituting equation (16) and (17) to (11) we will get the Newton – Raphson generation for  as follows:

 

  3n    x    2 2  n 1   n     2n  3  x    

     

2

2

(19)

2 F F 1 Finally for the PJ , the first and second derivative will be  and 2  2 P 2 , and  P P P therefore the Newton-Raphson generation for P would be Pn1  Pn  P with subject to conN

straint

 P =1 and P J

J

0

J=1

4.

Implementation Issue using Matlab These section describe how to implement the NPAG in mathlab. For example, equation (18) will be implemented using kidney surgery data from Rumah Sakit Syaiful Anwar (RSSA), Malang – East Java, 1997. The data is about Gentamicyn dose for urology patients, shown table 1.


6

Table 1: Gentamicyn dose for urology patients in RSSA Malang (1997)

Dose

No

Patien’s Name

1 SUD 2 HAR 3 SOF 4 SUM 5 ISM 6 AHM 7 ROC 8 DJU 9 MUL 10 URI 11 SOE 12 SUW 13 ADI

Dose Per Body Given weight Time mg/kg duration (mg) 80 1,19 5,00 80 1,57 4,33 80 1,60 5,75 80 1,01 5,97 80 1,48 4,75 80 1,51 5,83 80 1,82 4,00 80 1,82 4,00 80 1,43 4,17 80 1,86 4,58 80 1,45 5,00 80 1,60 4,83 80 1,33 6,67

Time for taking sample*

First sample ConcentraAver- Time tion** age dura(ug/ml) (ug/ml) tion 1,70 1,72 1,710 7,75 2,34 2,48 2,410 7,00 2,71 2,73 2,720 7,75 0,95 0,99 0,970 7,80 2,25 2,28 2,265 7,83 1,25 1,26 1,255 8,00 1,82 1,81 1,815 6,92 2,02 1,98 2,000 6,92 2,04 1,95 1,995 7,00 1,30 1,24 1,270 7,17 1,61 1,48 1,545 7,00 3,38 3,33 3,355 7,00 1,46 1,45 1,455 8,00

Second sample Concentration Aver** age (ug/ml) (ug/ml) 0,92 0,92 0,920 1,61 1,59 1,600 2,34 2,35 2,345 0,75 0,78 0,765 1,40 1,38 1,390 0,87 0,86 0,865 0,92 0,91 0,915 1,12 1,15 1,135 1,29 1,20 1,245 0,88 0,91 0,895 1,06 1,00 1,030 3,14 3,09 3,115 0,95 0,94 0,945

* Time for taking sample is calculated from IV gentamicyn ** Sample is taken 2 times

From table 1, we extract data from data konsentrasi cuplikan 1 and cuplikan 2 in 13 patien represented in variable named “data1” and “data2” in matlab. “Data1” and “data2” in matlab representing x in equation (18), while  variabley is random data representing by “ B “. So the source code is: data1 = [1.710, 2.41, 2.72, 0.97, 2.265, 1.255, 1.815, 2, 1.995, 1.270, 1.545, 3.355, 1.455]; data2 = [0.92, 1.61, 2.34, 0.75, 1.40, 0.87, 0.92, 1.12, 1.29, 0.88, 1.06, 3.14, 0.95]; B1(1,1)= rand(1,1); for i=1:500 B1((i+1),1)=B1(i,1)+(((data1(1)-B1(i,1))+(data2(1)-B1(i,1)))/2); if (B1((i+1),1)-B1(i,1)) < 0.1 break end B1(i,1)=B1(i+1,1); end

If we run these source code, it can convergence in 3 iteration, with resulting  :1.3150. We can compare it with NPEM that convergence in 175 iteration. [5].

5. Conclusion From experiment we can make conclusion that, newton raphson method can make significant performance improvement than NPEM. NPAG proves effectively and efficiently describe non – parametric patient kidney surgery data using an adaptive grid approach. Coupling between non-parametric method and state-of-the-art primal – dual interior – point using Newton-Raphson method prove maximization of likelihood on each grid is faster than NPEM. This combination makes NPAG more efficient and effective and requires much less memory than the last generation of Non-parametric Expectation Maximization


7

(NPEM), and therefore faster computation feasible on PC. NPAG is already demonstrate the result as a mixture of density of kidney surgery data which lay on grids. References [1] Dempster, A.P., Laird, N.M., dan Rubin, D.B. (1977), “Maximum Likelihood from Incomplete Data via the EM Algorithm”,Journal Royal Statistical, Society, Series B: Vol. 39, hal.1-21. [2] Leary, R. (2001) An Adaptive Grid Non Parametric Approach to Pharmacokinetic and Dynamic (PK/PD) Population Models, Proceedings of the Fourteenth IEEE, Symposium on Computer-Based Medical Systems. [3] Schumitzky, A. (1991), “Nonparametric EM Algorithm for Estimating Prior Distribution”, Applied Mathematic and Computation, Vol. 45, hal.143-157. nd [4] Casella, G. dan Berger, R.L. (2002), Statistical Inference, 2 Edition, Duxbury, New York. [5] Prastyo, D. (2008). Pemodelan Farmakokinetika Populasi dan Individu menggunakan Algoritma EMNonParametrik dan Analisis Bayesian, ITS, hal 58.