**From Michael_Looby@sandwich.pfizer.com Fri Oct 4 03:14:19 1996
**

I have encountered a strange error while checking the nmtran pred code listed below. When I run this code (without the estimation step) on a VMS Alpha, I get the following NONMEM error:

PROGRAM TERMINATED BY OBJ, ERROR IN ELS

WITH INDIVIDUAL 1 (IN INDIVIDUAL RECORD ORDERING)

VAR-COV OF DATA FROM INDIVIDUAL RECORD ESTIMATED TO BE SINGULAR

MESSAGE ISSUED UPON EVALUATION OF OBJ. FUNCTION

*** Line 1 *** which calculates A3 from the negative sum of A1 and A2 seems to be causing the problem. When I replace LINE 1 with LINE 2 (which gives the result of the calculation that should be performed in LINE 1) the code runs without error

Can anyone explain this?

Best wishes

Mick Looby

-------------------------------------------------

$PROBLEM TEST

$INPUT ID TIME CON=DV ND

$DATA TEST.DAT

$PRED

A1 = THETA(1)

A2 = THETA(2)

LA1 = THETA(3)

LA2 = THETA(4)

LA3 = THETA(5)

; *** LINE 1 ***

A3 = -(A1 + A2)

; *** LINE 2 ***

;A3 = -9.95

TAU = 8

TI = TIME - (ND - 1) * TAU

EX1 = EXP(-LA1 * TI)

EX2 = EXP(-LA2 * TI)

EX3 = EXP(-LA3 * TI)

EX11 = 1 - EXP(-LA1 * TAU * ND)

EX22 = 1 - EXP(-LA2 * TAU * ND)

EX33 = 1 - EXP(-LA3 * TAU * ND)

EX111 = 1 - EXP(-LA1 * TAU)

EX222 = 1 - EXP(-LA2 * TAU)

EX333 = 1 - EXP(-LA3 * TAU)

AA1 = A1 * EX1 * EX11/EX111

AA2 = A2 * EX2 * EX22/EX222

AA3 = A3 * EX3 * EX33/EX333

F = AA1 + AA2 + AA3

Y = F * EXP(ERR(1))

$THETA

6.75

3.2

10.7

0.136

9.17

$OMEGA .2

;$ESTIM

------------------------------------

Data to be simulated:

1 0 0 1

1 1 0 1

1 4 0 1

1 8 0 2

1 9 0 2

-----------------------------------

****

**From alison Fri Oct 4 10:46:26 1996
**

This interesting puzzle arises from the impossibility of storing decimal fractions in a binary computer; i.e., it is a rounding problem.

Although 9.95 and 6.75 and 3.2 are exact decimal fractions, they have an infinite number of binary digits as binary fractions, and cannot be expressed exactly in a digitial computer.

The important lines of the control stream are:

A1 = THETA(1)

A2 = THETA(2)

; *** LINE 1 ***

A3 = -(A1 + A2)

; *** LINE 2 ***

;A3 = -9.95

...

Y = F * EXP(ERR(1))

$THETA

6.75

3.2

It so happens that, in the data file, all records are observation events; there are no records with MDV=1. The first record, at TIME=ND=0, should result in F=0. But with proportional error model, predictions of 0 are not allowed and give rise to the VAR-COV error message. So NONMEM has done the appropriate thing.

I changed from proportional to additive error model Y=F+ERR(1) and was able to run the problem with both versions of A3. I also included a table, displaying PRED values.

With A3=-(A1 + A2), the first line of the table shows

ID TIME ND TI PRED RES WRES

.00E+00 0.00E+00 1.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00

With A3=-9.95, the first line of the table shows

ID TIME ND TI PRED RES WRES

.00E+00 0.00E+00 1.00E+00 0.00E+00 4.77E-08 -4.77E-08 -1.07E-07

Because PRED is not exactly 0, the proportional error model does not give rise to a NONMEM error message - but it is still the wrong model! It should not be used with tiny predictions, because the weights of these predicitions become excessive.

I did a run with the following code, and displayed TEMP3 in the table:

A3 = -(A1 + A2)

C9 = 9.95

TEMP = A3 + C9

Table output:

C9 A3 TEMP CON PRED RES WRES

9.95E+00 -9.95E+00 -4.77E-08 0.00E+00 0.00E+00 0.00E+00 0.00E+00

They still look identical.

I displayed values with the generated PRED using the debugger:

(dbx) p c9

c9 = 9.949999999999999

(dbx) p a3

a3 = -9.950000047683716

(dbx) p temp3

temp3 = -4.768371653085524e-08

Somehow, the value of 9.95 has been computed differently by code within NMTRAN (in which it is a double precision constant C9=9.95D0), and within NONMEM. In NONMEM, the values from the $THETA record are read by FORTRAN I/O routines and stored in the double precision vector THETA. When converting a constant to double precision, the compiler has "rounded down", but the I/O routines appear to have "rounded up". (I put "rounded" in quotes because I doubt that this is rounding as-such, but rather a low order bit that is either included or not; rounding in the world of binary fractions, perhaps.)

To summarize:

The code A3=-(A1 + A2) is correct ; the proportional error model should not be used.