MIN!f: A software package for optimization practice on
Microsoft® Excel spreadsheet


© AstroResearch


revised on January 12, 2005

SUMMARY

Min!f is a high nonlinear regression program. It can solve all the optimization problems that appear in the field of biosciences, pharmacokinetics, pharmacodynamics, enzyme-kinetic, etc. The user can use more than 10 algorithms, all derived from the quasi-Newton method. We show that, with Min!f, the need to have "good" first estimates of the parameters is no longer necessary, that is the user can enter even very poor estimates without non convergence occurs. Min!f can solve more than 50 row data with more than 20 adjustable parameters to fit, or a set of more than 10 simultaneous equations, including special cases for scalar or vectorial functions (like array functions). The user can define and store his own functions or even, his own algorithms by use of a special formula due to Huang. So, Min!f appears as a very user-friendly software and many options (graphs, statistics and report data) are easily available by only choosing the appropriate bar menu. Finally, it can run either on Apple® computer or IBM PC.

Key-words: Macintosh, IBM PC, spreadsheet, optimization, biosciences, quasi Newton, statistics, algorithms.

Download Min!f 3.5

Plan : 1. introduction - 2. Computational method and theory - 3. Description of the program - 4. Results : data from Hazelrig - data from Kaplan - 5. Example [ a. data - b. model - c. variables - d. initial graph - e. criterion function - f. compute quasi Newton algorithm - g. iterations - h. statistics - i. report - j. comparisons - k. final results] - 6. Discussion - conclusion - bibliography -



Min!f: A software package for optimization practice on Microsoft® Excel spreadsheet

1. Introduction

Many software packages [1 - 7] are available in the fields of bio sciences, pharmacokinetics, pharmacodynamics but none of them can be considered really as user-friendly, because they don't allow the user to choose in many ways the options of the procedure, especially those which are written in Basic [6, 7], Fortran,... Our aim was to write in a simple way a program devoted to the search of values of multi variable nonlinear functions without forcing the user to go always in the same way or the same sub routine. However, the goal was to make a full program which allows all the options that are desirable to make an entire search of the optimum including reports, graphics and statistics. Considering this requirement, we have used the well known Microsoft® Excel spreadsheet:

- It can run either on PC (OS/2®, Windows™) or Apple® computer without modifying anything on the program.
- The program (MACRO language) is very simple to learn and to write [with version 4.0 and 5.0 ; higher versions use a special Basic code].
- There are great capacities for the graphic options and exchanges among file (that is, the user can import data from other folders or import data from Lotus 1.2.3® or Resolve® file very easily.).
- The "Wysiwyg" model and the combinations with text, graphics and equations in a single document makes the software very user-friendly.

2. Computational methods and theory

The aim of all the optimization methods is to minimize a criteria function [8, 9], that is, e.g., in least squares, a sum of square like:

(1)

where j = 1,...n defines the number of data plots (in pharmacokinetics, a pair of tj,Cj with tj for the time and Cj for an observed concentration in the blood), w denotes the weighting factor for the jth observation (one can take often w = 0 -no weight-, w = 1 or w = 2.) x.i are the parameters to find. Some authors have defined other methods like Extended least squares [10] (with Min!f, the user can define all the weights he wants, see hereunder). f(tj,x.i) is the model equation defined by the user (in pharmacokinetics, often a sum of exponentials [11] like

(2)

i = 1 where Ce is the intercept and Cq is the slope). Such models can be also described in terms of linear system of differential equations [12] if we consider the compartment analysis model: for instance, if we take a two compartment open model, we have (FIGURE I) a central compartment (A) connected with 1 peripheral compartment (B) by a system of micro-rate


FIGURE I

constants k12 and k21; we can consider also that at the time t = 0, if a drug is administered intravenously in bolus, all the drug is carried in the central compartment (A), so that the amount of the drug is:

(3a)  (3b)

(where VA is the volume of distribution in the 1st compartment; CAo is the amount of drug in the 1st compartment while CBo denotes the amount of drug in the 2nd compartment). This leads to the system of differential equations:

(4)

With Min!f, the user can define the criteria function in term of equations like (2) or system like (4) with initial conditions as above (3a and 3b). The explicit solution of (4) is given by:

(5)

with  (5a) and   (5b)

Once the model has been defined, the user must define explicitly the objective function to be solved for the variables x.i, in the equation (1). So, let x.i = {1,...,q} be a column vector of several variables that define the right hand of (1). If we take our model in (2), we see that the function is nonlinear in its 2 parameters and the vector of variables can be defined as:

 (6)

The solution can't be reached in a full step but in a certain number of iterations. At each iteration, namely k+1, the actual values of the estimates x.i are replaced by new values which must minimize the objective function f(x.i). This is done by a nonlinear procedure of fitting, namely a nonlinear optimization method. There are many numerical methods which are devoted to this procedure:

a)- The gradient methods [13] which need function values and gradient values at each iteration (that includes the Raphson-Newton algorithm [14], the Levenberg-Marquard algorithm [15,16] and all the quasi-Newton methods [17]). All of these
methods need the first partial derivatives of the function with respect to each parameter x.i and the second partial derivatives (for the Raphson-Newton algorithm [18]).

b)- The direct search methods [19] which need only function values, like the simplex algorithm. The 1st method suffers from two serious drawbacks:

- At first, it is difficult to compute and inverse the Jacobian matrix: we recall that the Jacobian is the matrix M of dimension (jxq) of all the first partial derivatives of x.i (so that, in our example, we have 4 columns formed by each of the x.i's). We can approach the analytic values of these derivatives by first differences [20], as in many applications, the user is unwilling or unable to provide expressions of these. Min!f uses the most simple expression for the evaluations, so:

(7)

Our experience shows that we can gain, so, a lot of time of computation without damages in the results.

Second, the fact that the method requires modifications [21] until it frequently fails to converge. Namely, and it is a always a great problem with the optimization methods, a good initial estimate of the solution may be provided [22], unless meaningless "solutions" and even divergence may occur: that is, a local minimum is found instead of the true optimum. The FIGURE II shows a function of 2 variables with a local minimum

(the Rosenbrock's function: see equation 8 (8)) plotted in one dimension.


FIGURE II

In order to overcome these difficulties, some suggestions has been made:

Scaling the Jacobian matrix [15, 16] by a factor which increases each aii, so that the Hessian matrix, which is the inverse matrix of second derivatives becomes:

(9)

where T denotes a transposition. This method is closed to the classical steepest-descent [23] when J is great (J>1) (where the Hessian matrix is the identity matrix I) and to the Raphson-Newton when J is small (J<1).

Finding the Hessian matrix scale by scale, without inversing the cross product of M: it is the basis for all the methods called quasi-Newton [24, 25] and, indeed, they have many advantages:

- The inversion matrix is not needed
- If the first approximation of the inverse of the Hessian, namely H-1 is a positive definite matrix, so will be the further estimates H2,...Hk+1 (where k is the actual iteration).
- Only gradients (found by equation 7) and function values must be needed. One can compute both the gradients and function values in a single pass. Moreover, a direction search is needed to ensure that a minimum is found (i.e.,). This search can be made with a relative accuracy, so that usually, it is said that this is non exact: however, we will show in the further section that the speed of the iterations is improved and that it seams to give as good final values as a classical and entire procedure of exact search. The TABLE I shows a flowchart of Min!f for computing a quasi-Newton algorithm. 


TABLE I

- See Fletcher, Broyden, Biggs, GreenstadtGoldfarbBrodlie and Huang for papers about quasi-Newton methods and algorithms.
- See the nonlin file for implementation of these algorithms in Min!f 3.

3. Description of the program

Min!f solves scalar or vectorial functions and allows the user to enter even far away values from the initial estimates, as we will see in the further section. The program is composed by a MACRO file called 'MQNI' and by a tutorial file called 'HelpMin'. By opening 'MQNI', one can show a new bar menu, with the different bar control of the MACRO drives. If, for example, we want to enter the equation (1) with equation (5) as f (tj, x.i), we can choose ENTER FUNCTION from the
FUNCTION menu [see the demo file for exemples]. The cell $A$2 is selected and we can write:

{=x. 1 *exp(-x.2*t. 1 )+x.3*exp(-x.4*t. 1)} (10)  
[Note carefully stake enter hook of the equation, written under matrix shape, see demo file
]


The response is NAME? because we have not yet named {x.i} and the independent variable t.1 (the time). Now, we invoke the bar called DATA and we select ENTER t.i. The screen shows another region of the document where we can write the data: first, the time in column A, then the concentration in column B, but one might have entered other t.i columns and the corresponding observed concentrations (obs.1, obs.2,...,obs.i or under matrix shape obs.a). We can now by selecting RETURN TO f enter the x.i's and corresponding values. Once the initial values has been computed, a dialog box appears to ask the user for computing a quasi-Newton algorithm or the Levenberg-Marquard (LM) algorithm (only available for vectorial functions). The user can also cancel the procedure if he wishes to change of quasi-Newton method or of 6 value for the derivative. For that, he can invoke the menu OPTION and then choose PREFERENCE: a new convergence or failure (3 warnings are available scaling incidents of procedure, so that the program is really interacting with the user). Once the convergence has been achieved, the final estimates of the parameters and the estimate of the Hessian matrix are computed. Then, some options are available (including: report of all the procedure, statistic tests [35] with residuals, run test, Akaike's information criterion [36], ±SD of variables and test F for testing several models). The user can enter:

- More than 50 row data and more than 20 parameters to estimate
- A set of more than 10 simultaneous equations to fit
- Equalities or inequalities constraints (x.i <= a or a < x.i < b)
- Special functions, including quartic functions on the kind:

 (11) where A can be a diagonal matrix or other matrix... or a set of more than 50 simultaneous linear equations which may be of the kind:  (12) where the t.i's can be scalar or vectorial arguments. It is possible, too, to minimize a function with only 1 variable.
- A function to be maximized with equalities and inequalities constraints like the well known following problem:

Maximize f(x.i) = x.1*x.2*x.3 subject to the constraints: 

with: x.4 = X.1 + 2*X. 2 + 2*X.3

as described by Rosenbrock [38] and Box [37]. This problem is solved by Min!f in 3 iterations (the solution is: f(x.i) = 3300 with: x.1 = 20, x.2 = 11 and x.3 = 15).

4. Results

We have considered 2 examples, the first taken from Hazelrig and al. [39], the other taken from Giltinan and Ruppert [40].

a)- In the first example, Hazelrig illustrates the MFIT library by a set of data taken from Sheppard and Martin: they have proposed that the following equation could be appropriate to fit the exchange of potassium between plasma and erythrocytes in vitro:

(13)

The data appear on TABLE II in the two first columns.


TABLE II

As pointed out before, the user has only to enter equation (13) like:

{=x.2*exp(x.1*t.1)+x.3} (17)

and the data (the screen shows exactly the same look that the TABLE II); then the procedure is exactly the same that has been depicted above (for the entry of the independent variable t.i) and the definition of the objective function. In order to test fully the performances of all the algorithms, we have tried the following initial estimates:

x.1 = x.2 = x.3 = 0

and it can be pointed out that Hazelrig and coll. don't specify their own initial estimates.


TABLE III
-
(see Fletcher, Broyden, Biggs,
GreenstadtGoldfarbBrodlie and Huang for papers about quasi-Newton methods and algorithms. )

The TABLE III shows the results for the most important algorithms implemented in Min!f. It will be noted that the most powerful method, namely the Broyden's formula (known also as BFGS), which has been quite unanimously recommended for the past 15 years, has failed with our initial estimates of {x.i} but we note also that the rank one formula and the Oren's method has been successful; this will be discussed hereunder. The TABLE IV shows a typical output from Min!f (with the data of Sheppard and Martin).


TABLE IV

b)- In the second example, we have worked with the data of Giltinan and Ruppert, taken from Kaplan. These data appear on the TABLE V. We don't take care at first to the initial values of the parameters and they are arbitrary taken at {1 ;1 ;...;1} for both rate constants and intercepts; then, we simply try a 1, 2 and 3 open model compartment and give the statistics which permit to compare ad equally the different models. We see at the TABLE VI that the data are in good agreement with a two compartment open model (see above for the definition): both the AIC and the F test show us that a 3rd compartment is not necessary, in despite of a decrease of the SS. We see, moreover, that it is no longer necessary to take care of precise values for the first estimates. The Broyden's algorithm was used to perform these tests and the TABLE VII shows the results with all the algorithms of Min!f.


TABLE V (bolus IV)


5. Example

We are going to take for example the data of the TABLE V to give an idea of Min!f's possibilities. Procedure consists in:

a)- seizing data:

It is necessary to introduce data from the cell $A$132 [ R132C1] of the spreadsheet; dependent variable at first [time], then in the following column, the data.


FIGURE III
(entry data)


b)- entry of the model:

One can begin with a model in 1 compartment. The equation of this model is:

{ =x.1*EXP (-x.2*t.1)}  (18)

It is arranged in the cell $A$2 [ R2C1] of the spreadsheet, by activating the dialog box giving access to the bookshop of pharmacocinetik. The {...} denotes a matrix form input.


FIGURE IV
(activating libraries)



FIGURE V
(bottom of the dialog box : pharmacokin. library of build-in functions)

What gives the following aspect on the sheet :


FIGURE VI
(the model equation)

Where one can see equation in the bar of formula, its name [f] and names of the variable. Here, x.1 = D/V and
x.2 = K10
where V is the volume of distribution, D the dose administered in bolus IV and K10 the constant of elimination [ Kel].

c)- entering variables:

It is necessary to begin by entering values for the dependent variables x.1 and x.2. With the majority of programs, the values which we are going to take would lead to « data error. » Due to the algorithm developed by Oren and Spedicato, the minimum of the function criterion f will be reached quickly [although data make that function is not quadratic]. In our case, we are going to propose values x.1 = x.2 = 1. This is done by selecting the ENTRY DATA of the bar colon


FIGURE VII
(entering  x.i)


What leads to the seizure box:


FIGURE VIII
(entering number os parameters x.i)


Here, you enter the number « 2 » [x.1 and x.2]. Then, you can enter the initial estimates of the parameters x.1 and x.2. A new seizure box appears where you enter the values [it will appear twice because we have two variables ; in our case, the choice is « 1 » twice].


FIGURE IX
(entering the initial estimates of x.1 and x.2)


At this moment, values entered the sheet and a dialog box appears, who suggests entering algorithm. Click on " cancel ".


FIGURE X
(entry of model, data and initial guesses values of parameters)


Why? Because, in this stage, equation 18 is not a criterion function; criterion function should be a sum  of squares to be minimized; we are still not it there. We are going at first to build a graph, associating data and model.

d)- building a simple graph:

We have to see that program recognizes the columns of dependent and independent variables; in other words, we have to build registers. For it, we shall call " T.I REGISTER " command from the DATA menu.


FIGURE XI
(t.i and o.i registering)


A new dialog box displays which one suggests to you entering the number of columns of dependent and independent variables: in this case, one seizes « 1 » and « 1 » every time. One takes charge not of the bottom where is shown " matrix pooled dated ".


FIGURE XII
(editing t.i's and o.i's vectors)


By clicking on OK, a new column is created with an estimation according to the initial values of dependent variable's x.i. A graph is constructed, in semi-log plots. For reasons connected to incompatibilities of graph between Excel 4.0 and Excel 98 and 2000, the reader is invited to choose the mode of representation himself.


FIGURE XIII
(semi-log plots of data - see figure III)


One can see that equation is indicated in the top of the graph and that the labels of x and y are there arranged according to the data of the figure XII. [labels of t.i and o.i].

e)- To build criterion function:

You have merely to order " TYPE of  F " of the menu " FUNCTION ".


FIGURE XIV
(type of criterion function)


You reach another dialog box where a choice is proposed; in this case [noisy data], you choose option 3.


FIGURE XV
(form of the criterion function)


Now, it is necessary to seize the factor of weighting schema, according to the equation   (19), with w = 0 , 1 or 2 or other expression [e.g. extended least squares of Sheiner (45)] one can choose ... For the purpose, we suggest to enter the value 2, which is classical. So that the criterion function is :

{=SUM(IF(o.1=0;0;1/o.1^2*(o.1-(x.1*EXP(-x.2*t.1)))^2))} (20)

where o.1 is the data vector (column B) ; t.1 is the independent vector of time ; and the model equation is still :

x.1*EXP(-x.2*t.1) (21, see equation 18)

so that 21 has the form of equation 1. The conditional test IF(0.1=0;0...)  must be computed in order to take care of missing values.

f)- compute the quasi Newton algorithm:

Now, we can compute the algorithm, according to our other in-line papers (1, 2, 3, 4). This is done by QUASI NEWTON of the bar menu SET.


FIGURE XVI
(compute algorithm)


so that it leads with the dialog box :


FIGURE XVII
(choose your method)


By selecting quasi-Newton, you will find the dialog box for the quasi Newton algorithms :


FIGURE XVIII
(quasi newton methods)


Although these algorithms are very successful, we advise you to choose by default Oren (33). It has the form :

Hi+1 = [Hi-{Hi.g.gT.Hi}/{gT.Hi.g}+t.{v.vT}].gama+{d.dT}/{dT.g}

with:                   v = {gT.Hi.g}^0,5.[d/{dT.g}-{Hi.d}/{gT.Hi.g}

and:             gama = {dT.g}/{gT.Hi.g}.(1-p)+{gT.d}/{gT.Hi.g}.p

where g is the gradient of function f(tj,x.i) [see for further information pharmacokinetic sample].

g)- iterate

This is done by the command « ITERATE » of the bar menu FIND.



FIGURE XIX
(begin iteration processs)


A warning informs one when the end appends. In this case [with an one model compartment IV bolus] :


FIGURE XX
(end of iteration procedure)


The same graph as before [figure XIII] can be observed :


FIGURE XXI
(semi-log plots - final ; see figure XIII et III)


The final estimated values for x.1 and x.2 are


FIGURE XXII
(SS and estimated values for x.i)


x.1 = 128.24 g/l ; x.2 = 0.089 h-1.

h)- statistics:

For more developments, see the section pharmacokin sample. Having calculated all the tests, a table of that kind appears:


FIGURE XXIII
(statistical reports for the one compartment open model)


We can see that the model equation does not fit « at best » the data : C2 = 59.95 [p < 0.001].

i)- report:

_____________________________________________
©Min!f, AstroResearch Program. Version 3.5, 2/2003           
        AstroResearch   
        one compartment open model bolus IV   
        24/2/2003   
        function f{x.i}   
        SUM(IF(o.1=0,0,1/o.1^2*(o.1-(x.1*EXP(-x.2*t.1)))^2))   
        initial value of f{x.i}   
        11,96450362   
        final value of f{x.i}   
        0,572916828   
        initial {x.i} values   
        1    1
        final {x.i} values   
        128,2437092    0,08981458
        ± sd of {x.i}   
        1,858113206    0,004351
        +/-IC (0.95)       
        1,23691783    0,002896   
        Hessian matrix       
        62,35138615    0,051456   
        0,051455539    0,000342   

        algorithm:    Oren   
        iterations       
        23       
        time (s)       
        1       
        step analysis       
n° iter    function        variables   
               
1    11,95874        1,084021    0,924179
2    11,92958        1,510219    0,724623
3    11,77091        2,89718    0,339362
4    11,63496        3,602345    0,221483
5    11,59733        3,771857    0,20173
6    11,57987        3,85047    0,193966
7    10,61341        5,48162    0,055841
8    9,747068        5,801647    0,029348
9    9,670625        6,101085    0,005019
10    9,266756        5,998246    0,013453
11    9,258899        6,080822    0,0124
12    8,366047        10,86115    0,022672
13    5,678725        58,48175    0,141963
14    4,553134        116,8627    0,280341
15    2,341504        73,51424    0,082949
16    2,167239        69,04179    0,067982
17    2,152391        69,99789    0,070844
18    0,813099        103,9629    0,080558
19    0,574582        126,1408    0,08891
20    0,572918        128,187    0,089793
21    0,572918        128,1882    0,089788
22    0,572917        128,2154    0,089801
23   
0,572917        128,2436    0,089813
_____________________________________________________________________

j)- Comparison with the other algorithms - model with two compartments.

We give here only comparison with a model with two compartments. The reader will have only to redo all the stages, since b)-. Model to be entered corresponds to the equation 10. Procedure remains exactly the same. Initial values for the vector of variable x.i is {1; 1; 1; 1}.  We give below the graph and the statistics for this model. The final values of parameters can be read to the table VI.



FIGURE XXIV
(semi-log plots with a two compartment open model)


As we see it, equation corresponding to the model 17 was replaced by that of the model 10.We see now statistical tests.


FIGURE XXV
((statistical reports for the two compartments open model)

As we see it, SS is of 0.57 ± 0.058 on one side and of 0.06 ± 0.007 on the other side. A test t can be made on these data [t16 = 8.72, with p < 0.001 ; df = (10+8)-2].


FIGURE XXVI
(t test)


The figure XXVI shows a test t comparing the value of SS of two models.

k)- global results:

The table VI summarizes all the results from 1 to 3 compartments.


TABLE VI


l)- working with Excel's solver

Excel has a remarkable solver. Nevertheless, it does not arrange any option allowing to have access to results fathomed as those that we saw in the paragraph i)-.It is very easy to work from Min!F with theExcel's solver. If we resume the data of Kaplan with a system two compartments's system, we owe:

-To define criterion function as seen in e)-. Here is an example of screen with the window of the solveur opened.


FIGURE XXVIII
(Excel's solver)


One can see that the cell to be defined is named f, the variable cells under shape of the vector x.i; value to be achieved is the minimum of the criterion function, perceived in the bar of formula and identical to the equation 10. Here is the result which shows the Excel's solver.


FIGURE XXIX
(solver's results)


As one sees it, only the final results of the estimations of x.i are noted; no detail of the procedure of analysis with regard to the paragraph i)-. In comparison, here are results obtained with Min!F:

____________________________________________________
©Min!f, AstroResearch Program. Version 3.5, 2/2003                   
        _          
        _           
        25/2/2003           
        function f{x.i}           
        SUM(IF(o.1=0,0,1/o.1^2*(o.1-(x.1*EXP(-x.2*t.1)+x.3*EXP(-x.4*t.1)))^2))           
        initial value of f{x.i}           
        11,92911707           
        final value of f{x.i}           
        0,060579317           
        initial {x.i} values           
        1    1    1    1
        final {x.i} values           
        64,27664719    0,066106128    108,746132    0,304062173
        ± sd of {x.i}           
        14,83426237    0,005741591    14,59591849    0,094212375
        +/-IC (0.95)           
        9,078439293    0,003513804    8,93257492    0,057657153
        Hessian matrix           
        10255,18862    3,909051322    -4324,714349    57,98257933
        3,909051322    0,001664048    -1,736165625    0,021513682
        -4324,714349    -1,736165625    5985,297225    -15,84062334
        57,98257933    0,021513682    -15,84062334    0,370064249
                   
        algorithm:    Oren       
        iterations           
        28           
        time (s)           
        8           
        step analysis           
n° iter    function        variables       
                   
1    11,78443555        1,805215404    0,518880442    1,805215404    0,518880442
2    11,56757029        2,455380823    0,281292438    2,455380823    0,281292438
3                       
4    9,50278208        2,504039082    0,007646889    2,504039082    0,007646889
5    9,478978823        2,504043783    0,010265967    2,504043783    0,010265967
6    9,478919068        2,504181634    0,010131698    2,504181634    0,010131698
7    9,352201244        3,037104154    0,020383399    3,037104154    0,020383399
8    8,284830383        6,237269441    0,038158198    6,237271132    0,038158198
9    7,186056667        9,625647672    0,029284789    9,625645005    0,032489206
10    4,79801018        18,87235797    0,039358985    18,86757979    0,077986627
11    4,344869602        23,3857923    0,035096837    23,31778141    0,100882099
12    1,356457102        45,08360933    0,060352387    45,46193428    0,157170677
13    0,172565219        77,01745671    0,067292755    77,43703651    0,274774875
14    0,115459176        79,5140322    0,072369442    79,28327762    0,277121635
15    0,108424421        80,44267142    0,072051634    85,12471633    0,293180506
16    0,063413385        68,00570613    0,066979951    104,7057722    0,328389591
17    0,062004916        68,15800799    0,067736884    106,021571    0,330879156
18    0,061836463        68,16059613    0,06751719    106,4982399    0,331870751
19    0,061758929        68,18091451    0,067588217    107,2296844    0,333071267
20    0,060640324        65,19619441    0,066430242    108,3510494    0,308870868
21    0,060617127        65,04468521    0,06643171    108,1571048    0,307658489
22    0,060598556        64,6985848    0,066260991    108,512552    0,30725669
23    0,060596341        64,65172706    0,066258542    108,6524347    0,307260371
24    0,060579777        64,33934269    0,066129774    108,7395016    0,304604851
25    0,060579686        64,32551662    0,066127734    108,7255752    0,304476046
26    0,060579342        64,29273187    0,06611349    108,7372088    0,304167377
27    0,060579329        64,28178523    0,066107889    108,7441543    0,304124415
28    0,060579325        64,27728892    0,06610668    108,7503972    0,304103052
_________________________________________________________________________________

We give below an example of graph of evolution of the minimization of the criterion function.


FIGURE XXX
(minimisation of f : 28 iterations)


Note that missing value corresponds to a Hessian's reset due to a singularity.

6. Discussion

We have found in the literature a number of algorithms for the least squares estimation or for the minimization or the search of stationary values for several functions. All the authors have pointed out the necessity of beginning the procedure of optimization with good start estimates. With Min!f, the user can almost take merely faraway values without taking care of some drawbacks:

- Non convergence: this occurs when the initial values are far from correct; for instance, with our 2nd example, it is obvious that negative values of ai (i.e. : x.2 and x.4) are irreconcilable with normal linear pharmacokinetic modeling. So, one must take at least positive initial values for both slopes and intercepts.

-Singularity of matrices: It is possible, too, that the inverse of the actual estimation of the Hessian becomes singular in one iteration: a warning shows then 3 possibilities:

1 —> resetting the Hessian (i.e. replacing the actual G-1 array by the identity matrix I)
2 —> trying another algorithm (with the option of resetting (y/n) together G-1)
3 —> stop.

We must say that, in our 2 examples, no reset was necessary with the Oren algorithm, while 2 resets were done for the Broyden's method. Last, the selected equation does not reasonably fit the data: the user can change the equation during the optimization procedure without loosing, for instance, the actual estimates of x.i.

-Slow convergence: this occurs when the convergence criterion is too stringent; with Microsoft® Excel, Min!f works with absolute precision (that is 15 digits). However, as the nonlinear regression procedures nearly always move downhill, so it is possible that it converges at the bottom of one valley just over a ridge, as seen in the FIGURE III.


FIGURE XXVII

As we saw in the 1st example, the BFGS algorithm has failed contrary to the rank one formula (SR1): Conn et al have pointed out that the SR1 formula generates more accurate Hessian approximations in a number of circumstances [34]. Moreover, we have seen, that, with poor initial estimates, the Levenberg-Marquard method [13, 15, 16] fails very often after some iterations: the reason of why is, precisely, that it is the Hessian itself which is computed by this kind of method: in the quasi-Newton methods, the inverse of the Hessian is first approximated by any definite positive matrix (often the identity matrix, I) and, then, we obtain a sequence of matrices Hk which progress to the true Hessian asymptotically, if and only if Hk+1 remains definite positive: this is done internally by the program by an optimal choice of the parameter a (see above) by the "line search direction": this search direction is not exact, that is we don't locate exactly the minimum between xi and yi, in equation (14) and (15).

(14)

(15)

yx.ik denotes a first essay which will show if the line direction search is "good", in the way of the objective function (e.g., a sum of squares: SS) is minimized, i.e., SSk+1 < SSk. It has been pointed by BOX [13] that failure in convergence to locate the true minimum stems from the fact that, indeed, many methods tend to locate the minimum along a line too precisely. This point is very important and it can explain that, in despite of very poor estimates (e.g., our 1st example from Hazelrig), the convergence has been reached without any problem. So, we must point out once more that the choice of initial estimates of parameters sufficiently closed to the optimal values is not longer necessary with Min!f.

Other difficulties can occur frequently in nonlinear regression with the ill conditioning problem. It arises often when fitting sums of exponentials, when the first partial derivatives (see equation 7) are collinear: we have to care care of this drawback when the estimates of the parameters are substantially different with large standard errors for the components of x.i. We see that, with the 3rd compartment model taken from Kaplan's data, we have been able to detect differences of less than 30% in intercepts (see TABLE VII), even if, at last, the AIC was not in agreement with a 3rd compartment model (we recall that, usually in pharmacokinetics, it is not convenient to take care of parameters which exhibit less than 30% of difference compared with the others). Nevertheless, Min!f seems to detect even little differences among a number of parameters (we report the reader to Jennrich and Ralston [8] for further comments about the ill conditioning problem).

We want now to point out some of the definite advantages which occur with the use of a spreadsheet:

-the user can always see both the data and the optimization's procedure "in progress";
-the program is inter-active and the user can stop the procedures without spending data or the actual parameter values; the algorithm can be change always and the actual estimate of the inverse of the Hessian matrix can be reinitialized if singularity occurs, without holding all the procedure;
- a number of constraints can be computed, including boundary values on the parameters themselves (see above the function f(x.i)=x.1*x.2*x.3 from Rosenbrock [38]).
-the program is speedy: 1 iteration takes less than 1 s for 3 at 4 parameters; we have tried a 16 parameters model like:

(16)

with for initial estimates: xk = 1 and xk+1 = -1.2. The reader can find this function in [33]; the solution is f(x.i)=0, corresponding to x.i = 1;...1. The program takes 4 at 5 s for 1 iteration and, using Oren's algorithm, the convergence occurs in 25 iterations.

Conclusion

We must say that all the tests were performed on an Apple® computer, with CPU 68030 at 40 Mhz, FPU 6882 and 10 Mb of RAM [it was in 1994...].
Finally, it is astonishing to note how infrequent are the works devoted to the spreadsheets to fit data in despite of all their great advantages: we have found only 3 papers [41, 43, 44] but none about nonlinear regression fitting [42].
Our program can be downloaded, free of charge, by clicking here: it contains the program itself "MQNI' and a full didactic file (HelpMin). The user must possess Microsoft® Excel version 4.0 or later. All the program takes about 185 Kb on the hard-disk.

Bibliography

1. D.W.A. Bourne, Multi-forte, a microcomputer program for modeling and simulation of pharmacokinetic data, Comput. Methods Prog. Biomed. 23 (1986) 277-281

2. D.W.A. Bourne, Boomer, a simulation and modeling program for pharmacokinetic and pharmacodynamic data analysis, Comput. Methods Prog. Biomed. 29 (1989) 191-195

3. G.D. Alien, Modfit: a pharmacokinetics computer program, Biopharm. Drug Dispos. 11 (1990) 477-498

4. R.G. Duggleby, a nonlinear regression program for small computers, Anal. Biochem. 110(1981) 9-18

5. A.H.J. Scaf, pharmacokinetic analysis with Rugfit: an interactive pharmacokinetic computer program, Biopharm. Drug Dispos. 9 (1988) 415-446

6. K. Yamaoka, T. Nakagawa, H. Tanaka et al, a nonlinear multiple regression program multi2 (Bayes), based on bayesien algorithm for microcomputer, J. Pharm. Dyn. 8 (1985) 246-256

7. K. Yamaoka, T. Nakagawa, a nonlinear least squares program based on differential equations, multi (Runge), for microcomputers, J. Pharm. Dyn. 6 (1983) 595-606

8. R.I. Jennrich, M.L. Ralston, fitting nonlinear models to data, Ann. Rev. Biophys. Bioeng. 8 (1979) 195-238

9. J.M. Chambers, fitting nonlinear models: numerical techniques, Biometrika 60 (1973) 1-13

10. A.H. Thomson, A.W. Kelman, B. Whiting, evaluation of nonlinear regression with Extended Least Squares: simulation study, J. Pharm. Sci. 12 (1985) 1327-1330

11. J. Powers, statistical analysis of pharmacokinetic data, J. Vet. Pharmacol. Therap. 13(1990) 113-120

12. R.I. Jennrich, P.B. Bright, fitting systems of linear differential equations using computer generated exact derivatives, Technometrics 4 (1976) 385-392

13. M.J. Box, a comparison of several current optimization methods, and the use of transformations in constrained problems, Computer J. 7 (1966) 67-77

14. H.O. Hartley, the modified Gauss-Newton method for the fitting of nonlinear regression functions by least squares, Technometrics 3 (1961) 269-280

15. K. Levenberg, a method for the solution of certain nonlinear problems in least squares, Quart. Appl. Math. 2 (1944) 164-168

16. D.W. Marquardt, an algorithm for least squares estimation of nonlinear parameters, J. Soc. Ind. Apll. Math. 11 (1963) 431-441

17. M.J.D. Powell, recent advances in unconstrained optimization, Math. Programming 1 (1971) 26-57

18. N.R. Draper, H. Smith, an introduction to nonlinear estimation, in Applied Regression Analysis, pp. 458-474 (Wiley, New-York, 1968, 1981)

19. J.A. Nelder, R. Mead, a simplex method for minimization, Computer J. 7 (1965)308-313

20. G.W. Stewart, a modification of Davidon's minimization method to accept differences approximations of derivatives, J. Ass. Comp. Mach. 14 (1967) 72-83

21. C.G. Broyden, a class of methods for solving nonlinear simultaneous equations, Math. Comp. 19 (1965) 577-593

22. H.J. Motulsky, LA. Ransnas, fitting curves to data using nonlinear regression: a practical and non mathematical review, F.A.S.E.B. J. 1 (1987) 365-374

23. H.B. Curry, the method of steepest descent for nonlinear minimization problems, Quart. Apll. Math. 2 (1944) 258-261

24. C.G. Broyden, quasi-Newton methods and their application to function minimization, Math. Comp. 21 (1967) 368-381

25. W.C. Davidon, variable metric methods for minimization, A.E.C. Research and Development Report, Argonne National Laboratories (Rev. TID-4500, 14th ed)(1959)

26. R. Fletcher, M.J.D. Powell, a rapidly convergent method for minimization, Computer J. 6 (1963) 163-168

27. C.G. Broyden, the convergence of a class of double-rank minimization algorithms. 2. The new algorithm, J. Inst. Maths Applies 6 (1970) 222-231

28. M.C. Biggs, minimization algorithms making use of non quadratic properties of the objective function, J. Inst. Maths Applies 8 (1971) 315-327

29. J. Greenstadt, variations on variable metric methods, Math. Comp. 24 (1970) 1-22

30. D. Goldfarb, a family of variable metric methods derived by variational means, Math. Comp. 24 (1970) 23-26

31. K.W. Brodlie, an assessment of two approaches to variable metric methods, Math. Programming 12 (1977) 344-355

32. H.Y. Huang, unified approach to quadratically convergent algorithms for function minimization, J.O.T.A. 5 (1970) 405-423

33. S.S. Oren, self-scaling variable metric (SSVM) algorithms, Manag. Sci. 5 (1974) 863-874

34. A.R. Conn, N.I.M. Gould, Ph.L Toint, convergence of quasi-Newton matrices generated by the symmetric rank one update, Math. Programming 50 (1991) 177-195

35. H.G. Boxenbaum, S. Riegeiman, R. Elashoff, statistical estimations in pharmacokinetics, J. Pharmacokinet. Biopharm. 2 (1974) 123-148

36. H. Akaike, a new look at the statistical model identification, I.E.E.E. Trans. Aut.Contr. 19(1973)716-723

37. M.J. Box, a new method of constrained optimization and a comparison with other methods, Computer J. 8 (1965) 42-52

38. H.H. Rosenbrock, an automatic method for finding the greatest or least value of a function, Computer J. 3 (1960) 175-184

39. J.B. Hazelrig, M.E. Turner Jr, E.Ackerman, a function minimization package (MFIT) for nonlinear parameter estimation providing readily accessible maximum likelihood estimates, Comput. Biomed. Res. 11 (1978) 51-64

40. D.M. Giltinan, D. Ruppert, fitting heteroscedastic regression models to individual pharmacokinetic data using standard statistical software, J. Pharmacokinet. Biopharm. 5 (1989) 601-614

41. R.C. Shumaker, H. Boxenbaum, G.A. Thompson, Absplots: a Lotus 1.2.3 spreadsheet for calculating and plotting drug absorption rates, Pharmac. Research 4 (1988) 247-248

42. H. Delboy, nonlinear fitting in pharmacokinetics with Microsoft® Excel spreadsheet, Int. J. Biomed. Comp, 37, 1-14, 1994

43. L.P. Brion, use of a Lotus 1-2-3 spreadsheet to plot the confidence and the prediction intervals in linear regression analysis

44. C. Sundqvist, K. Enkvist, the use of Lotus 1-2-3 in statistics, Comput. Biol. Med. 17(1987)395-399

45. Sheiner, L.B. and Beal, S.L. (1985), "Pharmacokinetic Parameter Estimates from Several Least Squares Procedures: Superiority of Extended Least Squares," Journal of Pharmacokinetics and Biopharmaceutics, 13, 185 -201.