Bayesian Priors

Setting up Prior Information

Prior information is important for Bayesian analysis. By default, priors to THETAS are multivariate normal, and to OMEGAS and SIGMAS are inverse Wishart distributed. Alternatively, a residual variance in the form of its square root, may be modeled via THETA (a sigma-like Theta parameters is set up in example 2). See $PRIOR for prior specifications. Here we discuss the setup for most Bayesian analysis purposes.

To set up the $PRIOR NWPRI statement such as

1

$PRIOR NWPRI NTHETA=4, NETA=4, NEPS=1 NTHP=4, NETP=4, NEPP=1

one should note

NTHETA=number of THETAs to be estimated;
NETA=number of Etas (OMEGAs) to be estimated (and is to be described by an NETAxNETA OMEGA matrix);
NEPS=number of EPSILONs (SIGMAs) to be estimated (and is to be described by an NEPSxNEPS SIGMA matrix);
NTHP=number of THETAs which have a prior;
NETP=number of OMEGAs with prior;
NEPP=number of SIGMAs with prior (NM73). Before NM73, the NEPP option was ignored, as supplying priors for SIGMAs was not activated, so that the $THETA records list the parameters, in order, the following:
NTHETA of initial THETAs;
NTHP of Priors to THETAS;
Degrees of freedom to each OMEGA block Prior;
Degrees of freedom to each SIGMA block Prior;

the $OMEGA records list the variances, in order, the following:

NETAxNETA of initial OMEGAS;
NTHPxNTHP of variances of Priors to THETAS;
NETPxNETP of priors to OMEGAS, matching the block pattern of the initial OMEGAS;

the $SIGMA records list the variances, in order, the following:

NEPSxNEPS of initial SIGMAS
NEPPxNEPP of priors to SIGMAS, matching the block pattern of the initial SIGMAS (NM73).

The following snippet shows details of setting up a THETA, OMEGA, and SIGMA.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44


$THETA 2.0 2.0 4.0 4.0 ; Initial THETAs
$OMEGA BLOCK(4) ; Initial Parameters for OMEGA
0.4
0.01 0.4
0.01 0.01 0.4
0.01 0.01 0.01 0.4
$SIGMA 0.1

$PRIOR NWPRI NTHETA=4, NETA=4, NEPS=1, NTHP=4, NETP=4, NEPP=1

; Prior information of THETAS (NTHP=4 of them)
$THETA (2.0 FIX) (2.0 FIX) (2.0 FIX) (2.0 FIX)

; Variance to prior information of THETAS (NTHPxNTHP=4x4 of them).
; Because variances are very large, this means that the prior
; information to the THETAS is highly uninformative.  Note that the
; order of $THETA values among the THETA records, and the order
; of $OMEGA values among the OMEGA records, is very important,
; But $THETAs and $OMEGAs can be interspersed.
$OMEGA BLOCK(4)
10000 FIX
0.00 10000
0.00  0.00 10000
0.00  0.00 0.0 10000

; Prior to OMEGA (NETPxNETP=4x4 if them)
$OMEGA BLOCK(4)
0.2 FIX
0.0  0.2
0.0  0.0 0.2
0.0  0.0 0.0 0.2

; Set degrees of freedom of OMEGA Prior (one value per OMEGA block)
; Uninformative Omega prior is designated by having a DF that is equal to
; the dimension size of the Omega block.
$THETA (4 FIX)

; Prior to SIGMA (NEPPxNEPP=1x1 if them)
$SIGMA 0.05 FIX

; Set degrees of freedom of SIGMA Prior (one value per SIGMA block)
; Uninformative SIGMA prior is designated by having a DF that is equal to
; the dimension size of the Sigma block.
$THETA (1 FIX)

By default, the number of prior experiments is 1. However, perhaps you have more than one previous study, and you wish to average their contribution, forming a composite average set of prior parameters to influence the present analysis. In this case, add NEXP=n to the $NWPRI record above, where n is the number of experiments. Then add the prior information of each additional study with additional $THETA, $OMEGA, and $SIGMA statements. The order is as following.

$THETA records list the parameters, in order, the following:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10


      NTHETA of initial THETAs
      Exp 1:
      NTHP of Priors to THETAS
      Degrees of freedom to each OMEGA block Prior
      Degrees of freedom to each SIGMA block Prior
      Exp 2:
      NTHP of Priors to THETAS
      Degrees of freedom to each OMEGA block Prior
      Degrees of freedom to each SIGMA block Prior
      ...

The $OMEGA records list the variances, in order, the following:

1
2
3
4
5
6
7
8


  NETAxNETA of initial OMEGAS
  Exp 1:
  NTHPxNTHP of variances of Priors to THETAS
  NETPxNETP of priors to OMEGAS, matching the block pattern of the initial OMEGAS
  Exp 2:
  NTHPxNTHP of variances of Priors to THETAS
  NETPxNETP of priors to OMEGAS, matching the block pattern of the initial OMEGAS
  ...

The $SIGMA records list the variances, in order, the following:

1
2
3
4
5
6


  NEPSxNEPS of initial SIGMAS
  Exp 1:
  NEPPxNEPP of priors to SIGMAS, matching the block pattern of the initial SIGMAS
  Exp 2:
  NEPPxNEPP of priors to SIGMAS, matching the block pattern of the initial SIGMAS
  ...

As of NM73, you can use more informative names as follows:

$THETAP for theta priors;
$THETAPV for variance to theta priors;
$OMEGAP for omega priors;
$OMEGAPD for degrees of freedom (or dispersion factor) for omega priors;
$SIGMAP for SIGMA priors;
$SIGMAPD for degrees of freedom (or dispersion factor) for SIGMA priors.

This allows you to intersperse these records at will in the control stream files, but it also gives NMTRAN an alternative source for values to NTHETA, NETA, NTHT, NETP, NEPS, and NEPP that is typically given in the $PRIOR NWPRIOR record. However, if these values are also listed in ~$PRIOR NWPRI~, then these values are chosen over what is surmised from the informatively labeled THETA/OMEGA/SIGMA records. Thus, the above control stream file could be structured as follows, with the various records in any order, and a shortened $PRIOR record (in the following example uninformative priors are used):

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40


$PRIOR NWPRI

; Prior information of THETAS (NTHP=4 of them)
$THETAP (2.0 FIX) (2.0 FIX) (2.0 FIX) (2.0 FIX)

$THETA 2.0 2.0 4.0 4.0 ; Initial THETAs
$OMEGA BLOCK(4) ; Inital Parameters for OMEGA
0.4
0.01 0.4
0.01 0.01 0.4
0.01 0.01 0.01 0.4

; Set degrees of freedom of SIGMA Prior (one value per SIGMA block)
$SIGMAPD (1 FIX)

;intial parameters to sigma
$SIGMA 0.1



; Set degrees of freedom of OMEGA Prior (one value per OMEGA block)
$OMEGAPD (4 FIX)

; Prior to OMEGA (NETPxNETP=4x4 if them)
$OMEGAP BLOCK(4)
0.2 FIX
0.0  0.2
0.0  0.0 0.2
0.0  0.0 0.0 0.2


; Variance to prior information of THETAS (NTHPxNTHP=4x4 of them).
$THETAPV BLOCK(4)
10000 FIX
0.00 10000
0.00  0.00 10000
0.00  0.00 0.0 10000

; Prior to SIGMA (NEPPxNEPP=1x1 if them)
$SIGMAP 0.05 FIX

Informative prior information may come from a previous study. Typically, they are used as follows. The theta priors for the present analysis are obtained from the estimates of THETAs from the previous study. For example, in the report file of the previous study:

1
2
3
4


FINAL PARAMETER ESTIMATE
 THETA - VECTOR OF FIXED EFFECTS PARAMETERS
         TH 1      TH 2      TH 3      TH 4
         1.64E+00  1.57E+00  7.58E-01  2.35E+00

would be placed in the present study control stream file as:

1
2


$THETAP   (1.64 FIXED) (1.57 FIXED) (0.758 FIXED)
          (2.35 FIXED)

The variance-covariance to theta priors of the present analysis are obtained from the variance-covariance submatrix pertaining to the theta estimates from the previous study. For example, the information in the report file of the previous study:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10


COVARIANCE MATRIX OF ESTIMATE
            TH 1      TH 2      TH 3      TH 4
 TH 1
+        2.33E-03
 TH 2
+        4.76E-04  2.86E-03
 TH 3
+        7.87E-04  1.27E-04  5.35E-03
 TH 4
+        7.80E-05  2.36E-04  1.76E-03  2.98E-03

would be placed in the control stream file of the present study as:

1
2
3
4
5


$THETAPV BLOCK(4)
        2.33E-03 FIXED
        4.76E-04  2.86E-03
        7.87E-04  1.27E-04  5.35E-03
        7.80E-05  2.36E-04  1.76E-03  2.98E-03

The OMEGA priors of the present analysis are obtained from the estimates of OMEGAs from the previous study. For example, from the report file of the previous study:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10


OMEGA - COV MATRIX FOR RANDOM EFFECTS - ETAS
            ETA1      ETA2      ETA3      ETA4
 ETA1
+        1.75E-01
 ETA2
+        8.33E-03  1.51E-01
 ETA3
+        2.98E-02  1.74E-02  2.41E-01
 ETA4
+       -8.05E-03  1.84E-02  5.14E-02  1.62E-01

and we can use it in the control stream of the present study:

1
2
3
4
5


$OMEGAP BLOCK(4)
    1.75E-01 FIXED
    8.33E-03  1.51E-01
    2.98E-02  1.74E-02  2.41E-01
   -8.05E-03  1.84E-02  5.14E-02  1.62E-01

Similarly for Sigma priors, the results of the previous study:

1
2
3
4


SIGMA - COV MATRIX FOR RANDOM EFFECTS - EPSILONS  ***
            EPS1
 EPS1
+        5.28E-02

is used in the present study as:

1

$SIGMAP (5.28E-02 FIXED)

Degrees of freedom

The degrees of freedom to the omega priors of the present analysis are at most the total number of subjects in the previous study. We use the following formula for selecting degrees of freedom (justification is described in appendix G of the NONMEM7 Technical Guide), which is appropriate when the previous analysis was an NLME (maximum likelihood) form:

$$ \text{DF}=2\left[ \frac{\hat{\Omega}_{\text{previous analysis}}}{\text{SE}(\hat{\Omega}_{\text{previous analysis}})} \right]^2 $$

$$ \text{DF}=2\left[ \frac{\hat{\Omega}_{\text{previous analysis}}}{\text{SE}(\hat{\Omega}_{\text{previous analysis}})} \right]^2 + 1 $$

to adjust for degrees of freedom loss in the estimate of OMEGA of the previous study.

For an OMEGA block, use the smallest DF calculated among the OMEGA diagonal estimates in that block.

A similar formula would apply for SIGMA priors, with the proviso that the DF be no larger than the total number of data points that apply for that sigma in the previous study (for example, if there are two SIGMAs, one for PK data, and another for PD data, then the sigma for PK data gets no more than total number of PK data points in the previous study).

If the previous analysis is Bayesian, the estimated OMEGAs (or SIGMAs) are expanded relative to the modal values from a MLE, and the ratio of estimated OMEGAs to the SE of OMEGAs from a Bayesian analysis is smaller that from an MLE (using Fisher information). As a result, the recommended transformation when incorporating results from a previous Bayesian analysis are as follows (where p is the dimension of the OMEGA, Q the covariance matrix from the previous analysis):

OMEGAPD: $$ \text{DF}=2\left[ \frac{\hat{\Omega}_{\text{previous analysis}}}{\text{SE}(\hat{\Omega}_{\text{previous analysis}})} \right]^2 + p + 3 $$

OMEGAP: $$ R=((\text{DF}-p-1)/\text{DF})*Q $$

where are is the matrix that you would place as OMEGAP. The justification of the above is based on the mean and variance formulas for an inverse Wishart (see appendix G of the NONMEM7 Technical Guide).

The SCALE option for $OMEGA(P) and $SIGMA(P) allows the user to conveniently input a re-scaling of the input matrix without having to hand calculate every element.

For convenient transfer of the entire information from a previous analysis to prior information for a subsequent analysis, see section I.101 priorget: Transfer Results of an Analysis to NMTRAN Prior Information (NM75).

As yet another alternative, the mean of the inverse OMEGAS or SIGMAS from MCMC Bayesian samples and their variances are such that they are similar to those obtained from an nlme analysis. Therefore, the inverse of the mean inverse OMEGA can be used as the input $OMEGAP without re-scaling, and the effective DF for $OMEGAPD would be similar to that of an NLME analysis:

$$ \text{DF}=2\left[ \frac{\bar{\Omega}^{-1}_{\text{previous analysis}}}{\text{SE}(\bar{\Omega}^{-1}_{\text{previous analysis}})} \right]^2 $$ Here $\bar{\Omega}$ denotes average.

The priorget utility will calculate the inverse of the average omega inverse from the non-negative iterations in the .ext file of an MCMC Bayesian analysis, and provide the appropriate prior information from these.

In NM74+, the degrees of freedom to the inverse Wishart algorithms for OMEGAs and SIGMAs may be any positive real number. Thus the inverse Wishart distribution can substitute for the inverse Gamma matrix distribution: the parameter beta is the rate parameter (with inverse units of the variable), and alpha is the shape parameter, for the Gamma distribution. This distribution to the inverse residual variance can be expressed with an equivalent Wishart distribution to the inverse residual variance: set 2alpha for the $SIGMAPD, and beta/alpha for $SIGMA. The Gamma distribution to the inverse residual variance is equivalent to the inverse Gamma distribution of the residual variance.