85 Nested Random Levels
1 Adding Nested Random Levels Above Subject ID (NM75)
The following section provides an improved method of modeling random levels above subject ID that differs from NONMEM 7.4 and earlier. This method will give similar results, but will provide standard errors to THETAs that reflect the variability in the highest random level, rather than at the lowest, subject level. It does this by pointing the Mu referencing to the eta of the highest random level, rather than at the individual level, and preventing level etas summing to 0 (centering), by setting $EST LEVCENTER=0.
Suppose you wish to model inter-site variability, or inter-trial variability, so that several subjects belong to a trial. An easy, method would be to use the $LEVEL feature. Consider the following control stream fragment, which in addition to inter-subject variability eta(1) for clearance (CL), there is inter-site variability eta(5):
$PK
MU_5=THETA(1)
MU_2=THETA(2)
MU_3=THETA(3)
MU_4=THETA(4)
CL=DEXP(MU_5+ETA(5)+ETA(1))
V1=DEXP(MU_2+ETA(2))
Q=DEXP(MU_3+ETA(3))
V2=DEXP(MU_4+ETA(4))
S1=V1
...
$LEVEL
SID=(5[1])
…
$EST METHOD=ITS INTERACTION PRINT=1 NSIG=2 NITER=500 SIGL=8 FNLETA=0 NOABORT
CTYPE=3 MCETA=0
LEVCENTER=0 ; Do not force level etas to sum to 0
$EST METHOD=1 INTERACTION PRINT=5 NSIG=2 NBURN=1000 NITER=500 SIGL=10
FNLETA=0 NOHABORT SLOW NONINFETA=1 MCETA=20
$COV MATRIX=R UNCONDITIONAL SIGL=10Note that LEVCENTER=0 so as not to force level etas to sum to 0. Notice also that FNLETA=0 is set, so that the etas reflect what were used in the estimation. If FNELTA=0 is not set, super ID eta values outputted using $TABLE will incorrectly differ with each subject, rather than averaged for each LEVEL item value.
Let us suppose that the data item named SID is the site ID. NONMEM needs to know that SID is to be associated with eta(5), and in turn eta(1) is nested within eta(5). The data file need not be sorted for super ID values. The $LEVEL record gives this information:
$LEVEL
SID=(5[1])such that SID is a super ID data item associated with eta(5) (inter-site eta), and eta(1) nests within eta(5) (5[1]). NONMEM will then perform appropriate summary statistics for eta(5), and make the appropriate constraints on eta(5), so eta(5) changes by site, that is, by every SID value change, and not by every ID value change. You may have additional parameters having site variability etas and their suitable nesting etas, such as for V1, Q, and V2:
$PK
MU_5=THETA(1)
MU_6=THETA(2)
MU_7=THETA(3)
MU_8=THETA(4)
CL=DEXP(MU_5+ETA(5)+ETA(1))
V1=DEXP(MU_6+ETA(6)+ETA(2))
Q=DEXP(MU_7+ETA(7)+ETA(3))
V2=DEXP(MU_8+ETA(8)+ETA(4))
S1=V1
...
$LEVEL
SID=(5[1],6[2],7[3],8[4])Perhaps in addition to SID, you have country ID, let’s call that data item CID. Perhaps there are several sites belonging to one country, some other sites belonging to another country, etc. This would provide a nesting level of 2 above that of ID, and is expressed as follows, for example (..\examples\superid2_*.ctl):
$PK
MU_9=THETA(1)
MU_10=THETA(2)
MU_11=THETA(3)
MU_12=THETA(4)
CL=DEXP(MU_9+ETA(9)+ETA(5)+ETA(1))
V1=DEXP(MU_10+ETA(10)+ETA(6)+ETA(2))
Q=DEXP(MU_11+ETA(11)+ETA(7)+ETA(3))
V2=DEXP(MU_12+ETA(12)+ETA(8)+ETA(4))
S1=V1
...
$LEVEL
SID=(5[1],6[2],7[3],8[4])
CID=(9[5],10[6],11[7],12[8])Thus, for clearance, eta(9) is the country variability that has nested in it the site variability eta(5), which in turn has nested in it the subject variability (the standard ID data) eta(1). When performing FOCE with $LEVEL, you must use the SLOW option in $EST, and MATRIX=R for the covariance step $COV should be selected. Notice in these examples the THETAs are mu referenced to the highest random level.
Sometimes you may wish to have inter-site variability, but not inter-subject variability. In the following example, residual SD is to change only with each site, and the inter-subject variance is set to 0:
$PK
MU_1=THETA(1)
MU_2=THETA(2)
MU_4=THETA(3)
CL=DEXP(MU_1+ETA(1))
V=DEXP(MU_2+ETA(2))
W=DEXP(MU_4+ETA(4)+ETA(3))
S1=V
$ERROR
IPRED=F
Y = F + F*W*EPS(1)
; Initial values of THETA
$THETA 0.5 1.0 -.5
;INITIAL values of OMEGA
$OMEGA BLOCK(2)
0.02
0.001 0.02
$OMEGA
(0.0 FIXED)
$OMEGA
0.2
$SIGMA
1.0 FIXED
$LEVEL
SID=(4[3])As of nm74, instead of making a nesting connection to an eta of 0 variance, you can reference eta 0, which tells NONMEM that there is no variability for that level underneath. Thus the above model can be re-written to:
$PK
MU_1=THETA(1)
MU_2=THETA(2)
MU_3=THETA(3)
CL=DEXP(MU_1+ETA(1))
V=DEXP(MU_2+ETA(2))
W=DEXP(MU_3+ETA(3))
S1=V
$ERROR
IPRED=F
Y = F + F*W*EPS(1)
; Initial values of THETA
$THETA 0.5 1.0 -0.8
;INITIAL values of OMEGA
$OMEGA BLOCK(2)
0.02
0.001 0.02
$OMEGA
0.2
$SIGMA
1.0 FIXED
$LEVEL
SID=(3[0])with no need to have an extra ETA as a dummy variance. Thus, NONMEM is informed that eta(3) is a random effect to W, associated with changes in SID, rather than the usual ID, and no other eta shares the variability for that parameter.
As of NM74, the syntax to giving the level mapping may use the TO(:)/BY pattern method. For example,
SID=(5[1],6[2],7[3],8[4])
CID=(9[5],10[6],11[7],12[8])may be expressed as:
SID=(5:8[1])
CID=(9:12[5])So, 5:8[1] means etas 5 to 8 nest into etas beginning at 1. The “to” value in the brackets is not necessary. However, you may wish to express a by step:
SID=(5:11by3[12])means etas 5 to 11, by steps of 3, nest into etas beginning at 12, also by steps of 3. So, etas 5, 8, 11 map into etas 12, 15, 18, respectively. If the by step value must differ between the two lists, then the by value must be given in both lists:
SID=(5:11by3[12by4])So, etas 5, 8, 11 nest into 12, 16, 20, respectively. Negative steps may also be used:
SID=(11:5by-3[20])means etas 11, 8, 5 map into etas starting at 20, also by steps of -3 (20, 17, 14).
Commas should only be used to separate x[y] structures, but commas may not appear within x or within y. Also, x may have a to(:)/BY pattern and y need only have a value indicating a starting number of the nested etas, as given in the examples above, but the converse may not be true. That is, the y may imply its to/by pattern from x, but x may not imply its to/by pattern from y.
Nesting below the subject ID as for previous versions of NONMEM, as shown for inter-occasion variability, example 7.
This method using $LEVEL is approximate to modeling data such that the SID data item is renamed as ID, and the the original ID data item is renamed as OCC (for occasion). Then, the 16 SID levels become 16 “super” subjects, each super subject having 50 occasions to represent the original 50 subjects within each level. This of course places considerable computational demand, and executes much more slowly than the $LEVEL method. The following example (superid3_20.ctl) shows the layout of how this would be done. Notice the $OMEGA … SAME(49) line, to allow for 49 additional occasions (formerly subjects) within each subject (formerly SID) sharing the repeated 49 omega blocks to represent “inter-occasion” (formerly inter-subject) variability.
$SIZES LVR=160
$PROB RUN#
$ABBR DERIV2=NO
$INPUT C OCC TIME DV AMT RATE EVID MDV CMT ROWNUM ID
$DATA superid3_occ.csv IGNORE=C
$SUBROUTINES ADVAN2 TRANS2
$ABBR REPLACE ETA(OCC_KA)=ETA(,4 TO 151 BY 3)
$ABBR REPLACE ETA(OCC_CL)=ETA(,5 TO 152 BY 3)
$ABBR REPLACE ETA(OCC_V)=ETA(,6 TO 153 BY 3)
$PK
MU_1=THETA(1)
MU_2=THETA(2)
MU_3=THETA(3)
KA=DEXP(MU_1+ETA(1)+ETA(OCC_KA))
CL=DEXP(MU_2+ETA(2)+ETA(OCC_CL))
V=DEXP(MU_3+ETA(3)+ETA(OCC_V))
S2=V
$ERROR
IPRE=F
Y = IPRE + IPRE*EPS(1)
$THETA 1.79067E-01 -5.31302E+00 -3.08220E+00
$OMEGA BLOCK(3) 5.23413E-02
3.14841E-02
-5.55179E-03 2.64664E-02
2.73335E-02 -6.08426E-03
$OMEGA BLOCK(3) 9.72772E-03
1.00860E-02
1.48793E-04 9.80457E-03
5.38548E-04 6.57755E-04
$OMEGA BLOCK(3) SAME(49)
$SIGMA ;[P]
2.99892E-03
$EST METHOD=ITS INTERACTION PRINT=1 NSIG=2 NITER=0 SIGL=8 FNLETA=0
NOABORT CTYPE=3 MCETA=0
$COV MATRIX=R UNCONDITIONAL SIGL=10The equivalent model using $LEVEL would be as follows (superid3_21.ctl):
$PROB RUN#
$INPUT C ID TIME DV AMT RATE EVID MDV CMT ROWNUM SID
$DATA superid3.csv IGNORE=C
$SUBROUTINES ADVAN2 TRANS2
$PK
MU_4=THETA(1)
MU_5=THETA(2)
MU_6=THETA(3)
KA=DEXP(MU_4+ETA(4)+ETA(1))
CL=DEXP(MU_5+ETA(5)+ETA(2))
V=DEXP(MU_6+ETA(6)+ETA(3))
S2=V
$ERROR
IPRE=F
Y = IPRE + IPRE*EPS(1)
$THETA 0.2 -4 -2
$OMEGA BLOCK(3)
0.1
0.001 0.1
0.001 0.001 0.1
$OMEGA BLOCK(3)
0.3
0.001 0.3
0.001 0.001 0.3
;Initial value of SIGMA
$SIGMA
0.1 ;[P]
$LEVEL
SID=(4[1],5[2],6[3])
$EST METHOD=ITS INTERACTION PRINT=1 NSIG=2 NITER=500 SIGL=8 FNLETA=0 NOABORT
CTYPE=3 MCETA=0
$EST METHOD=1 INTERACTION PRINT=5 NSIG=2 NBURN=1000 NITER=500 SIGL=10
FNLETA=0 NOHABORT SLOW NONINFETA=1 MCETA=20
$COV MATRIX=R UNCONDITIONAL SIGL=101.1 LEVWT=0 (default) (NM74)
By default, LEVWT=0, and weights each level value equally, regardless of number of subjects per level value. If you wish to weight according to number of subjects for that value, set LEVWT=1 on the $EST record.
1.2 LEVCENTER=[0|1] (no default) (NM75)
If LEVCENTER=1, this ensures the etas of super ID random levels sum to 0. In earlier versions of NONMEM, the default (and only) action was to ensure the etas of super ID random levels sum to 0. To obtain similar results as earlier versions of NONMEM, set LEVCENTER=1 on the $EST record.
1.3 LEVOBJTYPE=0 (default) (NM75)
If LEVOBJTYPE=0, then for IMP/ITS/FOCE/Laplace, the reported objective function (and what is used for estimation) is the average of two cycles of optimization. The first optimization floats all etas, and the second optimization fixes the super-ID etas to the appropriate average of etas per super-id group from the first optimization, and then floats the ID (and below) level etas. If LEVOBJTYPE=1, then the OBJ reported (and used for estimation) is only that of the second optimization cycle. The LEVOBJTYPE=0 provides the most consistent results between FOCE/Laplace and the EM and Bayes methods, particularly with super-ID levels higher than 1 (so, site ID and country ID, for example). The LEVOBJTYPE=1 provides an OBJ more consistent with theory.
By default, NONMEM sets the estimation process to SLOW when a $LEVEL option is used. As of nm751, you may set NOSLOW or FAST. While there is an increase in speed, but there may not be as good optimization progress as with the SLOW option. It may be best to use SLOW for the $EST step, and use NOSLOW or FAST for the $COV step.
2 Adding Nested Random Levels Above Subject ID in Versions NM73 and NM74
In previous versions of NONMEM providing the $LEVEL feature, the paradigm used to model nested random levels above subject ID differed from the method described in the previous section in two ways. First, the MU referencing pointed to the individual (ID) level rather than the top random level. Secondly, the etas to random levels above ID were constrained to sum to 0, with the OMEGA estimation corrected for the loss of 1 degree of freedom. This older method applies some additional constraints on the model, and offers greater stability, especially for EM/BAYES methods, if the number of distinct items in a level are few (<5).
This original method of $LEVEL modeling was designed to replicate the following paradigm. Consider defining separate THETAs for each value pertaining to a super ID data item, so that theta is shared only by the subjects with the particular SID value. This is suitable if there are not too many distinct values of the super ID data item, otherwise, the number of THETAs can become very large, and the analysis may take a considerable amount of time. This analysis method could be performed in earlier versions of NONMEM, but the many THETAs that needed to be mapped with the different levels could make the NM-TRAN code quite large and tedious to write. Fortunately since NM73, a series of substitution variable techniques and short-hand entries for initial values can be used, and this method is now easier to program in NM-TRAN.
Here is an example to code using separate THETAs pertaining to each value of the SID data item (example superid3_6):
$SIZES LTH=60
$PROB RUN#
$INPUT C ID TIME DV AMT RATE EVID MDV CMT ROWNUM SID TYPE L2
$DATA superid3_6.csv IGNORE=C
$SUBROUTINES ADVAN2 TRANS2
$ABBR REPLACE THETA(SID_KA)=THETA(,4 to 19)
$ABBR REPLACE THETA(SID_CL)=THETA(,20 to 35)
$ABBR REPLACE THETA(SID_V)=THETA(,36 to 51)
$ABBR DECLARE DOWHILE I
$ABBR DECLARE INTEGER NSID
$PK
MU_1=THETA(1)
MU_2=THETA(2)
MU_3=THETA(3)
NSID=16
THSUM_KA=0.0
THSUM_CL=0.0
THSUM_V=0.0
I=1
DO WHILE (I<=NSID)
THSUM_KA=THSUM_KA-THETA(I+3)
THSUM_CL=THSUM_CL-THETA(I+19)
THSUM_V=THSUM_V-THETA(I+35)
I=I+1
ENDDO
IF(SID<NSID) THEN
KA=DEXP(MU_1+ETA(1)+THETA(SID_KA))
CL=DEXP(MU_2+ETA(2)+THETA(SID_CL))
V=DEXP(MU_3+ETA(3)+THETA(SID_V))
ELSE
; for the last SID level, NSID, use the negative sum of the THETAs of the other SID levels,
; so that the sum of all THETAs is 0, that is, the super-nested average theta is 0.
KA=DEXP(MU_1+ETA(1)+THSUM_KA)
CL=DEXP(MU_2+ETA(2)+THSUM_CL)
V=DEXP(MU_3+ETA(3)+THSUM_V)
ENDIF
S2=V
$ERROR
IPRE=F
IF(TYPE==0) Y = IPRE + IPRE*EPS(1)
IF(TYPE==1.AND.SID<NSID) Y=THETA(SID_KA)+EPS(2) ; The fitting of the pseudo-data (TYPE>0)
IF(TYPE==1.AND.SID==NSID) Y=THSUM_KA+EPS(2) ; constrains the SID level THETAs to be
IF(TYPE==2.AND.SID<NSID) Y=THETA(SID_CL)+EPS(3) ; constrained, and modeled using extra
IF(TYPE==2.AND.SID==NSID) Y=THSUM_CL+EPS(3) ; Sigma variances 2-4.
IF(TYPE==3.AND.SID<NSID) Y=THETA(SID_V)+EPS(4)
IF(TYPE==3.AND.SID==NSID) Y=THSUM_V+EPS(4)
$THETA 0.2 -4 -2
(0.1)x15 (0.0 FIXED)
(0.1)x15 (0.0 FIXED)
(0.1)x15 (0.0 FIXED)
$OMEGA BLOCK(3) VALUES(0.1,0.001)
$SIGMA
0.1 ;[P]
$SIGMA BLOCK(3) VALUES(0.3,0.001) ; This is the inter-SID variance.
$EST METHOD=1 INTERACTION PRINT=1 NSIG=2 SIGL=10 FNLETA=0 NOHABORT NONINFETA=1 MCETA=20
$COV MATRIX=R UNCONDITIONAL SIGL=10Notice the use of variable replacement mapping ($ABBR REPLACE), short-hand entries for initial THETAs, OMEGAs, and SIGMAs, and that the sum of the THETAs to the SID data item are fixed to 0 by constraining the theta pertaining to the highest SID value (NSID) to be the negative sum of the THETAs to the other SID values (1 through NSID-1) using a DOWHILE loop.
For this method, some pseudo-data must be added to the data file.
The idea in doing this is to cause the following term to be added to the objective function:
\[\sum_{i=1}^{\text{NSID}} \left[\theta_i^T \Sigma^{-1} \theta_i + \log(|\Sigma|)\right]\]
Where \(\theta_i\) is the vector of SID THETAs, and \(\Sigma\) is the variance among the SID THETAs. For the above example, \(\theta_i\) is a 3x1 vector, one element each for KA (TYPE=1), CL (TYPE=2), and V (TYPE=3), for i=1 to NSID, where NSID is the number of possible values of SID, which in this example NSID=16. The \(\Sigma\) matrix is the 3x3 block matrix to Epsilons 2, 3, and 4. NONMEM is fooled into constructing the above term by use of the additional data records for which \(DV_{ij}=0\) (or nearly so), for which are modeled \(IPRED_{ij} = \theta(3+(TYPE-1)\cdot j+i)\), for i=1 to 16 SID values, and j=1 to 3 TYPE values. NONMEM thus adds, for each TYPE>0 data record, objective function value terms \((DV_i - IPRED_i)^T \Sigma^{-1} (DV_i - IPRED_i)\) that evaluates to \(\theta_i^T \Sigma^{-1} \theta_i\), and the control stream file places a dependency of the last \(\theta_i\) of each element (that is, each of the three TYPE’s) such that \(\sum_{i=1}^{\text{NSID}} \theta_i = 0\). The L2 data item allows NONMEM to assess correlation (hence off-diagonal elements to the SIGMA block) between the three TYPEs, within a given SID. Thus for the added data portion, NONMEM sees 16 “subjects”, one for each of the SID values, each of which have 3 “data points”, one for each PK parameter (TYPE).
The above problem can alternatively be coded more easily using the $LEVELS mapping of etas as follows (example superid3_1), without needing to add pseudo data to the data file:
$PROB RUN#
$INPUT C ID TIME DV AMT RATE EVID MDV CMT ROWNUM SID
$DATA superid3.csv IGNORE=C
$SUBROUTINES ADVAN2 TRANS2
$PK
MU_1=THETA(1)
MU_2=THETA(2)
MU_3=THETA(3)
KA=DEXP(MU_1+ETA(1)+ETA(4))
CL=DEXP(MU_2+ETA(2)+ETA(5))
V=DEXP(MU_3+ETA(3)+ETA(6))
S2=V
$ERROR
IPRE=F
Y = IPRE + IPRE*EPS(1)
; Initial values of THETA
$THETA 0.2 -4 -2
;INITIAL values of OMEGA
$OMEGA BLOCK(3)
0.1
0.001 0.1
0.001 0.001 0.1
$OMEGA BLOCK(3) ; Inter-SID variance
0.3
0.001 0.3
0.001 0.001 0.3
;Initial value of SIGMA
$SIGMA
0.1 ;[P]
$LEVEL
SID=(4[1],5[2],6[3])
$EST METHOD=ITS INTERACTION PRINT=1 NSIG=2 NITER=500 SIGL=8 FNLETA=0 NOABORT CTYPE=3 MCETA=0
LEVCENTER=1
$EST METHOD=IMP INTERACTION PRINT=1 NSIG=2 NITER=500 SIGL=8 FNLETA=0 NOABORT CTYPE=3 MCETA=0
ISAMPLE=300 MAPITER=0
$EST METHOD=SAEM INTERACTION PRINT=10 NSIG=2 NITER=100 SIGL=8 FNLETA=0 NOABORT CTYPE=3 MCETA=0
ISAMPLE=2 CONSTRAIN=0
$EST METHOD=IMP EONLY=1 INTERACTION PRINT=1 NSIG=2 NITER=5 SIGL=8 FNLETA=0 NOABORT CTYPE=3
MCETA=0 ISAMPLE=300 MAPITER=0
$EST METHOD=BAYES INTERACTION PRINT=10 NSIG=2 NBURN=1000 NITER=500 SIGL=8 FNLETA=0
NOABORT CTYPE=3
$EST METHOD=1 INTERACTION PRINT=5 NSIG=2 NBURN=1000 NITER=500 SIGL=10 FNLETA=0 NOHABORT
SLOW NONINFETA=1 MCETA=20
$COV MATRIX=R UNCONDITIONAL SIGL=10Notice that the MU referencing points to the individual level (MU 1,2,3), rather at the SID level as was done in the previous section (superid_21.ctl, MU 4,5,6). Also, LEVCENTER=1 is set in the present example (superid3_1.ctl), to center (constrain) the super ID etas to 0 (this is done by default in earlier versions of NONMEM, but in nm75, you must set LEVCENTER=1 to obtain similar results as in earlier versions).
Both methods (LEVCENTER=0 or LEVCENTER=1) provide similar estimates of population parameters. However, the standard errors of the THETAs in superid3_1 are smaller than from superid3_21. The standard errors in superid3_1 represent what one would get from information at the 800 subjects level, which is equiavalent to the error of THETA(1), THETA(2), and THETA(3) obtained in superid3_6, with the “difference of SIDS from overall mean” THETAs 4 to 51 (or the OMEGA(4,4), OMEGA(5,5), OMEGA(6,6)) absorbing the larger standard error based on the fewer (only 16 of them) SID levels. In contrast, the standard errors of the THETAs in superid3_21 incorporate the uncertainty due to variability between sites, and that there are only 16 of these sites. To get this standard error associated with site uncertainty from superid3_1, you can obtain this by dividing the site OMEGA by number of sites, and square-rooting it. For example, SQRT(OMEGA(4,4)/16)=site standard error of theta(1).