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

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

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:

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

i = 1 where C_{e}
is the intercept and C_{q} 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

constants k_{12} and k_{21}; 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:

(where V_{A} is
the volume of distribution in the 1st
compartment; CA_{o} is the amount of drug in the 1^{st}
compartment while CB_{o} denotes the
amount of drug in the 2^{nd} compartment). This leads
to the system of
differential equations:

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:

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:

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 1^{st} 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:

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.

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:

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,...H_{k+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.

- See Fletcher, Broyden,
Biggs, Greenstadt, Goldfarb, Brodlie and Huang for papers
about quasi-Newton methods and algorithms.

- See the nonlin file for implementation of
these algorithms in Min!f 3.

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 (t_{j},
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).

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:

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

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

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, Greenstadt, Goldfarb, Brodlie 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).

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 3^{rd} 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)

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:

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)

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

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)

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 = K_{10 }where V
is the volume of distribution, D
the dose administered in bolus IV and K_{10}
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)

x.2 = K

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)

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)

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)

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)

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)

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)

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)

FIGURE XV

(form of the criterion function)

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 :

g)- iterate

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

FIGURE XIX

(begin iteration processs)

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

where g is the gradient of
function f(tj,x.i) [see for
further information pharmacokinetic
sample]. 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)

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)

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)

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 **

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)

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)

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)

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)

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)

k)- global results:

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.

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)

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)

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.

____________________________________________________

©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.

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 2^{nd} 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.

As we saw in the 1^{st} 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 H_{k} which progress to
the true
Hessian asymptotically, if and only if H_{k+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).

y_{x.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., SS_{k+1}
< SS_{k}. 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 3^{rd}
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 3^{rd} 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:

with for initial estimates: x_{k} = 1 and
x_{k+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.

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.

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.