73  MU Referencing

Introduction to NONMEM 7.6.0

The new methods in NONMEM are most efficiently implemented if the user supplies information on how the THETA parameters are associated arithmetically with the etas and individual parameters, wherever such a relationship holds. Calling the individual parameters phi, the relationship should be

\[\phi_i = \mu_i(\theta) + \eta(i)\]

For each parameter i that has an eta associated with it, and \(\mu_i\) is a function of THETA. The association of one or more THETAs with ETA(1) must be identified by a variable called MU_1. Similarly, the association with ETA(2) is MU_2, that of ETA(5) is MU_5, etcetera. Providing this information is as straight-forward as introducing the MU_ variables into the $PRED or $PK code by expansion of the code.

For a very simple example, the original code may have the lines

CL=THETA(4)+ETA(2)

This may be rephrased as:

MU_2=THETA(4)
CL=MU_2+ETA(2)

Another example would be:

CL=(THETA(1)*AGE**THETA(2))*EXP(ETA(5))
V=THETA(3)*EXP(ETA(3))

which would now be broken down into two additional lines, inserting the definition of a MU as follows:

MU_5= LOG(THETA(1))+THETA(2)*LOG(AGE)
MU_3=LOG(THETA(3))
CL=EXP(MU_5+ETA(5))
V=EXP(MU_3+ETA(3))

Note the arithmetic relationship identified by the last two lines, where MU_5+ETA(5) and MU_3+ETA(3) are expressed. This action does not change the model in any way.

It is better to have a linear relationship between all THETAs and MU’s (as we shall see below)

MU_5= THETA(1)+THETA(2)*LOG(AGE)
MU_3=THETA(3)
CL=EXP(MU_5+ETA(5))
V=EXP(MU_3+ETA(3))

The above parameterization would also entail log transforming initial values of THETA(1) and THETA(3).

If the model is formulated by the traditional typical value (TV, mean), followed by individual value, then it is straight-forward to add the MU_ references as follows:

TVCL= THETA(1)*AGE**THETA(2)
CL=TVCL*EXP(ETA(5))
TVV=THETA(3)
V=TVV*EXP(ETA(3))
MU_3=LOG(TVV)
MU_5=LOG(TVCL)

This also will work because only the MU_x= equations are required in order to take advantage of EM efficiency. It is not required to use the MU_ variables in the expression EXP(MU_5+ETA(5)), since the following are equivalent:

CL=TVCL*EXP(ETA(5))=EXP(LOG(TVCL)+ETA(5))=EXP(MU_5+ETA(5))

but it helps as an exercise to determine that the MU_ reference was properly transformed (in this case log transformed) so that it represents an arithmetic association with the eta.

Again, it is preferable to re-parameterize so that the MU’s are linear functions of all THETAs:

LTVCL= THETA(1)+THETA(2)*LOG(AGE)
CL=EXP(LTVCL+ETA(5))
LTVV=THETA(3)
V=EXP(LTVV+ETA(3))
MU_3=LTVV
MU_5=LTVCL

An incorrect usage of MU modeling would be:

MU_1=LOG(THETA(1))
MU_2=LOG(THETA(2))
MU_3=LOG(THETA(3))
CL=EXP(MU_1+ETA(2))
V=EXP(MU_2+MU_3+ETA(1))

In the above example, MU_1 is used as an arithmetic mean to ETA(2), and a composite MU_2 and MU_3 are the arithmetic means to ETA(1), which would not be correct. The association of MU_x+ETA(x) must be strictly adhered to.

Once one or more THETAs are modeled to a MU, the theta may not show up in any subsequent lines of code. That is, the only usage of that theta may be in its connection with MU. For example, if

CL=EXP(THETA(5)+ETA(2))

So that it can be rephrased as

MU_2=THETA(5)
CL=EXP(MU_2+ETA(2))

But later, suppose THETA(5) is used without its association with ETA(2):

...
CLZ=THETA(5)*2

Then THETA(5) cannot be MU modeled, because it shows up as associated with ETA(2) in one context, but as a fixed effect without association with ETA(2) elsewhere. However, if

MU_2=THETA(5)
CL=EXP(MU_2+ETA(2))
…
CLZ=CL*2

Then this is legitimate, as the individual parameter CL retains the association of THETA(5) with ETA(2), when used to define CLZ. That is, THETA(5) and ETA(2) may not be used separately in any other part of the model, except indirectly through CL, in which their association is retained.

Suppose you have:

CL=THETA(5)+THETA(5)*ETA(2)

One should see this as:

CL=THETA(5)*(1+ETA(2))

So the way to MU model this is:

MU_2=1.0
CL=THETA(5)*(MU_2+ETA(2))

Which would mean that in the end, THETA(5) is not actually MU modeled, since MU_2 does not depend on THETA(5). One would be tempted to model as follows:

MU_2=THETA(5)
CL=MU_2+MU_2*ETA(2)

But this would be incorrect, as MU_2 and ETA(2) may not show up together in the code except as MU_2+ETA(2) or its equivalent. Thus, THETA(5) cannot be MU modeled. In such cases, remodel to the following similar format:

CL=THETA(5)*EXP(ETA(2))

So that THETA(5) may be MU modeled as:

MU_2=LOG(THETA(5))
CL=EXP(MU_2+ETA(2))

Again, for EM methods, better to re-parameterize as:

MU_2=THETA(5)
CL=EXP(MU_2+ETA(2))

And log transform the initial value of THETA(5).

Sometimes, a particular parameter has a fixed effect with no random effect, such as:

Km=THETA(5)

with the intention that Km is unknown but constant across all subjects. In such cases, the THETA(5) and Km cannot be Mu referenced, and the EM efficiency will not be available in moving this Theta. However, one could assign an ETA to THETA(5), and then fix its OMEGA to a small value, such as 0.0225 = 0.15² to represent 15% CV, if OMEGA represents proportional error. This often will allow the EM algorithms to efficiently move this parameter, while retaining the original intent that all subjects have similar, although not identical, Km’s. Very often, inter-subject variances to parameters were removed because the FOCE had difficulty estimating a large parametered problem, and so it was an artificial constraint to begin with. EM methods are much more robust, and are adept at handling large, full block OMEGAs, so you may want to incorporate as many etas as possible when using the EM methods.

You should Mu reference as many of the THETAs as possible, except those pertaining to residual variance (which should be modeled through SIGMA whenever possible). If you can afford to slightly change the theta/eta relationship a little to make it MU referenced without unduly influencing the model specification or the physiological meaning, then it should be done.

When the arithmetic mean of an ETA is associated with one or more THETAs in this way, EM methods can more efficiently analyze the problem, by requiring in certain calculations only the evaluation of the MU’s to determine new estimates of THETAs for the next iteration, without having to re-evaluate the predicted value for each observation, which can be computationally expensive, particularly when differential equations are used in the model. For those THETAs that do not have a relationship with any ETA’s, computationally expensive gradient evaluations must be made to provide new estimates of them for the next iteration.

If you provide a MU reference to THETAs associated with ETAS whose OMEGA value is fixed to 0 (and if you do not turn these MU-reference off with a MUM=N designation, see below), these THETAs will also be evaluated by gradient evaluations. However, as of NONMEM 7.4, these gradients will be evaluated using analytical derivatives, which are usually faster and more accurate (see I.50 The FAST Option for use with FOCE/ITS and Differential Equation ($DES) Models (NM74) to read how this is done). This will automatically evaluate its gradient using the analytical eta derivatives, and may increase the speed of analysis, although trial and error is recommended for each case. Set MUM to N for those THETAs you want not mu referenced (for example MUM=N(4,5) means to not use mu-reference for THETAs 4 and 5), and their derivatives will then be evaluated by finite difference.

There is additional increased efficiency in the evaluation of the problem if the MU models are linear functions with respect to THETA. As mentioned in the previous examples above, we could re-parameterize such that

MU_5=THETA(1)+THETA(2)*LOG(AGE)
CL=EXP(MU_5+ETA(5))
MU_3=THETA(3)
V=EXP(MU_3+ETA(3))

This changes the values of THETA(1) and THETA(3) such that the re-parameterized THETA(1) and THETA(3) are the logarithm of the original parameterization of THETA(1) and THETA(3). The models are identical, however, in that the same maximum likelihood value will be achieved. The only inconvenience is having to anti-log these THETAs during post-processing.

The added efficiency obtained by maintaining linear relationships between the MU’s and THETAs is greatest when using the SAEM method and the MCMC Bayesian method. In the Bayesian method, THETAs that are linearly modeled with the MU variables have linear relationships with respect to the inter-subject variability, and this allows the Gibbs sampling method to be used, which is much more efficient than the Metropolis-Hastings (M-H) method. By default, NONMEM tests MU-THETA linearity by determining if the second derivative of MU with respect to THETA is nearly or equal to 0. Those THETA parameters with 0 valued second derivatives are Gibbs sampled, while all other THETAS are M-H sampled. In the Gibbs sampling method, THETA values are sampled from a multi-variate normal conditional density given the latest PHI=MU+ETA values for each subject, and the samples are always accepted. In M-H sampling, the sampling density used is only an approximation, so the sampled THETA values must be tested by evaluating the likelihood to determine if they are statistically probable, requiring much more computation time.

As much as possible, define the MU’s in the first few lines of $PK or $PRED. Do not define MU_ values in $ERROR. Have all the MU’s particularly defined before any additional verbatim code, such as write statements. NM-TRAN produces a MUMODEL2 subroutine based on the PRED or PK subroutine in FSUBS, and this MUMODEL2 subroutine is frequently called with the ICALL=2 settings, more often than PRED or PK. The fewer code lines that MUMODEL2 has to go through to evaluate all the MU_s’ the more efficient.

MU parameters should be completely defined on every call to PK. While conditional assignments to MU are permitted, make sure to define them for all possibilities, such as:

IF(GROUP==1) THEN
MU_2=THETA(3)
ELSE
MU_2=THETA(4)
ENDIF

Time dependent covariates, or covariates changing with each record within an individual, cannot be part of the MU_ equation. For example

MU_3=THETA(1)*TIME+THETA(2)

should not be done. Or, consider

MU_3=THETA(2)*WT

Where WT is not constant within an individual, but varies with observation record (time). This would also not be suitable. However, we could phrase as

MU_3=THETA(2)
CL=WT*(MU_3+ETA(3))

where MU_3 represents a population mean clearance per unit weight, which is constant with time (observation record), and is more universal among subjects. The MU variables may vary with inter-occasion, but not with time.

Suppose we have a situation where WT has an unknown power term associated with it modeled as THETA(3) in this example:

CL=THETA(2)*WT**THETA(3)*EXP(ETA(1))

Normally, we could efficiently linear model this as follows:

MU_1=THETA(2)+THETA(3)*LOG(WT)
CL=EXP(MU_1+ETA(1))

with THETA(2) transformed into the log of clearance domain. However, if WT changes record by record within the individual, then LOG(WT) may not be in the Mu modeling. We would then remove the THETA(3)*LOG(WT) term from MU_1:

MU_1=LOG(THETA(2))
CL=WT**THETA(3)*EXP(MU_1+ETA(1))

And THETA(3) itself would not be MU modeled.

For NONMEM 7.2.0, NM-TRAN is programmed to detect some MU modeling errors. Nonetheless, the user should verify that these rules are followed.

Examples at the end of the document show examples of MU modeling for various problem types. Study these examples carefully. When transposing your own code, begin with simple problems and work your way to more complex problems.

At this point one may wonder why bother inserting MU references in your code. MU referencing only needs to be done if you are using one of the new EM or Gibbs sampling methods to improve their efficiency. The EM methods may be performed without MU references, but it will be several fold slower than the FOCE method, and the problem may not even optimize successfully. If you choose one of the new methods, and you do not incorporate MU referencing into your model, you are likely to be disappointed in its performance. For simple two compartment models, the new EM methods are slower than FOCE even with the MU references. But, for 3 compartment models, or numerical integration problems, the improvement in speed by the EM methods, properly MU modeled, can be 5-10 fold faster than with FOCE. Example 6 described at the end of the SIGL section is one example where importance sampling solves this problem in 30 minutes, with R matrix standard error, versus FOCE which takes 2-10 hours or longer, and without even requesting the $COV step. So, for complex PK/PD problems that take a very long time in FOCE, it is well worth putting in MU references and using one of the EM methods, even if you may need to rephrase some of the fixed/random (theta/eta) effects relationships. In addition, FOCE is a linearized optimization method, and is less accurate than the EM and Bayesian methods when data are sparse or when the posterior density for each individual is highly non-normal.

It cannot be stressed too much that MU referencing and using the new EM methods will take some time to learn how to use properly. It is best to begin with fairly simple problems, to understand how a particular method behaves, and determine the best option settings. When setting up a problem for the new EM methods, you should start out with some trial runs, and a limited number of iterations, and observe its behavior. Here are some starting points for the various methods:

$EST METHOD=ITS NITER=100
$EST METHOD=SAEM NBURN=500 NITER=500
$EST METHOD=IMP NITER=100 ISAMPLE=300

The convergence tests should not be used during trial runs. The convergence tests for the EM methods can be fooled into running excessively long, or ending the problem prematurely. For example, the iterations of SAEM are Markov chain dependent, and therefore, certain parameters may meander slowly. The convergence tester, if CITER and CINTERVAL are not properly set to span these meanderings, may never detect stationarity for all the parameters, and therefore may never conclude the analysis. For IMP, the parameters between iterations are less statistically correlated, and the convergence tester is a little more reliable for it.

NM-TRAN does some checking of MU statements. If you wish to turn this off (checking mu statements can take a long time for very large control stream files), then include the NOCHECKMU option on the $ABBR record:

$ABBR NOCHECKMU

1 MUM=MMNNMD

These options allow the MU reference equations for each theta to be optionally used or not used. By default, if a theta parameter is MU referenced, it will be used to facilitate theta parameter estimation. However, the user may “turn off” specific parameters so their Mu referencing is not used. M indicates that the parameter should be Mu modeled (assuming there is an association of a Mu for that theta, which the program will verify), and N indicates it should not be Mu modeled. In the above example, THETAs 1, 2, 5, 6 are MU modeled, and 3, 4 are not to be Mu modeled. D (for default) indicates you want the program to decide whether to MU model, useful for specifying back to a default option in a future $EST statement, if the present setting is N.

The MUM parameter can also be used to specify which THETAS are used in a mixture problem by marking the position with an X. For example:

MUM=DDDDX

Where THETA(5) is involved in mixture modeling (in a $MIX statement). This is only necessary for covariate dependent mixture models, such as:

$MIX
IF(KNOWGENDER==1) THEN
IF(GENDER==1) THEN
P(1)=1.0
P(2)=0.0
ELSE
P(1)=0.0
P(2)=1.0
ENDIF
ELSE
P(1)=THETA(5)
P(2)=1-THETA(5)
ENDIF

and it guarantees that the new estimation methods are aware of the proper parameters.

An alternative method for specifying MU modeled parameters is by using the following syntax:

MUM=v1(n1):v2(n2):v3(n3)...

Where v refers to a letter (N, M, D, or X), and n refers to a number list. For example, to specify THETAs 3, 5 through 8 to not be MU modeled, theta 2 is a population mixture parameter, and THETAs 6, 12 are to be MU modeled,

MUM=N(3,5-8):X(2):M(6,12)

THETAs not specified are given a default D designation.

2 GRD=GNGNNND

By default, if a theta parameter has a Mu associated with it, and its relationship to its Mu is sufficiently linear (the program tests this by evaluating the partial second derivative of MU with respect to theta), then the program will use Gibbs sampling for that parameter. However for Mu modeled parameters, the user can override these decisions made by the program, and force a given parameter to be Gibbs sampled (G), or Metropolis-Hastings sampled (N). In the above example, THETAs 1 and 3 are to be Gibbs sampled, and the other THETAs are M-H sampled. If the parameter is not Mu modeled, or its Mu modeling is turned off by an MUM option setting, the program performs an M-H sampling. D (for default) specifies you want the program to decide whether to use Gibbs sampling.

For SIGMA parameters, if a particular SIGMA is associated with only one data point type, and conversely, the data point type has only that one SIGMA parameter defining its residual error, and that data point type is not linked by an L2 item with any other data point types, then that SIGMA will by default be Gibbs sampled with a chi-square distribution. Otherwise, that SIGMA parameter will be sampled by Metropolis-Hastings. You can force Metropolis-Hastings by specifying an N. The first m letters of GRD refer to the m THETAs. Then, the m+1th letter refers to SIGMA(1,1), m+2 refers to SIGMA(2,2), etc (going along the diagonal of SIGMA). Not all THETAs and SIGMAs need to be designated. If just the THETAs are designated, for example then the designations for SIGMA are assumed to be D.

For example, for

Y=IPRED + (CMT-1)*IPRED**GAMMA*EPS(1) +(2-CMT)*IPRED*EPS(2)

And with no correlation set between SIGMA(1,1) and SIGMA(2,2), then both SIGMA(1,1) and SIGMA(2,2) will be Gibbs sampled.

Mixed homoscedastic/heteroscedastic residual errors are not Gibbs sampled:

Y=IPRED + IPRED*EPS(1) +EPS(2)

2.1 GRD=DDDDDDSSN

The S and D specification are used only for Monte Carlo EM methods. The S specification is optional, and can improve the speed of IMP, IMPMAP, and SAEM methods. Sometimes, users model parameters that could have been a Sigma parameter, but model them as Theta parameters instead, such as:

Y=IPRED+THETA(7)*IPRED*EPS(1)+THETA(8)*EPS(2)

These theta parameters are therefore “Sigma-like”, and are typically not MU referenced. To have the S designation, these THETAs are not allowed to be involved in evaluating the predicted function F, or compartment values A(x). Specifying theta parameters 7 and 8 as “sigma-like” in this example (note 7th and 8th position of S in the GRD option setting) indicates to the program that when it evaluates forward difference partial derivatives to these THETAs (which it must when etas are not associated with theta parameters), it does not have to re-evaluate the predicted function, which can be computationally expensive, especially if one of the differential equation solver ADVAN’s are used.

Another example in which the theta can be designated S, is something like the following:

$PK
EMAX=THETA(4)
EC50=THETA(5)
…
$ERROR
IPRED=EMAX*F/(EC50+F)
Y=IPRED+IPRED*EPS(1)
…
$EST … GRD=TS(4-5)

Note that EMAX and EC50 are not MU modeled, and they are not involved in the evaluation of F or A(x), so they would benefit from an S designation in terms of efficient evaluation in a Monte Carlo EM estimation, especially for differential equation problems.

An alternative method for specifying GRD modeled parameters is by using the following syntax:

GRD=t1v1(n1): t2v2(n2): t3v3(n3)...

Where t refers to a parameter type (T for theta, S for SIGMA), v refers to a letter (S, D, or N), and n refers to a number list. For example, to specify THETAs 3, 5 through 8 to be Gibbs samples, theta 4 is sigma-like, and SIGMAs 1-3 are to be Metropolis-Hastings processed,

GRD=TG(3,5-8):TS(4):SN(1-3)

THETAs and SIGMAs not specified are given a default D designation. The SN() designation is also used by EM methods to not determine the derivatives of the objective function with respect to the SIGMAs analytically (which is faster), but numerically.