Structural parameters and QRPEM PML example

QRPEM (Quasi-Random Parametric Expectation) is an importance sampling-based EM engine in Phoenix NLME that is based on the general EM paradigm of sampling the empirical Bayesian posterior of each subject given the current estimate of the fixed and random effect parameters, and then updating these parameters based on simple statistics computed from the samples. In particular, the fixed effect updates are based on posterior sample means of the random effect samples. In the most ideal case, each random effect is associated with a unique set of one or more fixed effects, and the updates to the fixed effects can be made from simple means or a simple linear regression model on the means of the corresponding sampled random effects.

Due the specific fixed effect update algorithm used by QRPEM, the correspondences between fixed and random effects in the PML stparm statements that define structural parameters in the model must be clear and restricted to a few acceptable types that are compatible with the fixed effect update methodology. Thus some stparm statements in models that are acceptable for other engines cannot be used directly with QRPEM. For example, non-linear covariate models within a PML stparm statement are perfectly acceptable for other engines, but will not work with QRPEM. Attempting to execute such a model with QRPEM will result in the error message (in the Core Status file in the Results tab)

  FATAL ERROR IN QRPEM/IMPEM!
   Model not suitable for QRPEM analysis
   Possibly non linear covariate model or
   some other unimplemented feature

Similarly, time varying covariates inside a stparm statement can be used freely with other engines but not QRPEM. Attempting to run such a model will result in an error message that the selected engine cannot work with covariate effects where the covariate(s) have multiple values within a given subject, and will identify the variable covariate.

However, in such cases there is almost always an easy workaround that is based on splitting the stparm statement or block into two parts, a conforming part that remains inside the stparm and is acceptable to QRPEM, and a second nonconforming part, which is inserted into the body of the model outside the stparm and contains the offending features.

This section describes the splitting technique and gives some examples. It also points to two ‘ready-to-go’ working example projects (Remifen.phxproj and QRPEM_with_time_varying_covariate.phxproj) that illustrate the method. All such example NLME projects are in the directory <Phoenix_install_dir>\application\Examples\NLME. You can save a copy of the Examples directory (installed with Phoenix) to your Phoenix project directory via the Project Settings in the Phoenix Preferences dialog.

The QRPEM method works most naturally and efficiently when all structural parameters in a model can be expressed simply in terms of fixed effects and random effects in such a way that either the structural parameter or its log is a normal random variable, and the mean of the normal random variable is a linear function of fixed effects and random effects. This is essentially the same concept as “mu-modeling” in NONMEM. When all the structural parameters are expressed this way, a very simple and efficient EM update formula allows the fixed effects to be updated in by estimating the means of the empirical Bayesian posterior distribution of each subject’s random effects by random or quasi-random importance sampling. A simple linear regression of the estimated means with the fixed effects (including any linear covariate models) as predictors provides the update.

For example, consider a simple volume of distribution structural parameter V. The most common way to model this in an stparm statement (because V is inherently nonnegative) is one of the two alternative, but mathematically equivalent, lognormal forms.

a) stparm(V=tvV*exp(etaV))

b) stparm(V=exp(tvlogV+etaV))

But also the simple additive parameterization:

c) stparm(V=tvV+etaV)

fits the ‘mu-modeled’ framework of acceptable QRPEM structural parameters, although technically this additive parameterization might be more suitable for a structural parameter that can be positive or negative.

In both a) and b), log(V) is the normally distributed structural parameter, where in c) V is normally distributed. Clearly b) is mathematically equivalent to a) with tvV=exp(tvlogV). In the log normal cases, either formulation a) or b) is can be used with QRPEM, which does the log transforms automatically when it encounters case a). This contrasts with NONMEM, where only b) can be used in the mu-modeled case.

A simple extension allows the mean of the structural parameter to include a linear function of covariates. For example, suppose V includes a covariate on weight (here we assume a mean of 75 on weight, and use wt-75 as the centered covariate):

d) stparm(V=exp(tvlogv+coefwt*(weight-75)+etaV))

e) stparm(V=tvV*exp(coefwt*(weight-75)+etaV))

Then clearly log(V)~N(tvlogv+coefwt*(weight-75), OmegaV), and the mean is a linear function of the fixed effects tvlogV and coefwt, so the conditions for simple EM updates apply.

But there are a number of cases where this simple update formula breaks down.

Case 1 — Bare fixed effects

If a structural parameter enters the model as a ‘bare’ fixed effect that is not paired or associated with a random effect, for example stparm(V=tvV), with no associated random effect, then that bare fixed effect tvV cannot be estimated from the mean of a sampled empirical Bayesian posterior distribution, since there is no posterior associated with V. Such bare fixed effects must be estimated from a direct numerical optimization of a log likelihood function involving those parameters. Details are not presented here, but it is simply noted that such bare fixed effects are allowed within QRPEM ‘stparm’ statements, which automatically sets up the necessary log likelihood optimizations.

However, there are other cases that also cannot be handled, even if all the fixed effects are paired with random effects. These require some restructuring to make them work in QRPEM. The two most important are a) nonlinear covariate models, and b) covariate models with time-varying covariates, whether linear or not in the fixed effects in the covariate model.

Case 2 — Nonlinear covariate models

Suppose a covariate model of the form:

f) stparm(V=tvV*(1+coefwt*(weight-75))*exp(etaV))

is desired rather than one of the forms given in d) or e) above. Now the mean of log(V) is log(tvV+coefwt*(weight-75)), which cannot be expressed as a linear function of fixed effects. This is a common choice of covariate model, but cannot be accommodated within the simple “mu-modeled” EM update framework.

All engines, other than QRPEM, allow such non-linear covariate models to be included within a ‘stparm’ statement. However, if QRPEM encounters such a ‘stparm’ statement, it will exit with an error message as described above. But the above nonlinear covariate model can be run in QRPEM after a simple restructuring. The non-linear covariate model must be broken into two parts, a standard log normal ‘mu-modeled’ part:

stparm(Vbase=tvV*exp(etaV))

and a non-linear part in the body of the model outside the stparm statement:

V=Vbase*(1+coefwt*(weight-75))

Note this reformulated ‘split’ covariate model is obviously mathematically equivalent to the original model. Also, the common similar alternative version:

g) stparm(V=(tvV+coefwt*(weight-75))*exp(etaV))

can be rewritten as:

stparm(V=tvV*(1+coefwt/tvV)*(weight-75)*exp(etaV))

and the ‘split version’ looks like:

stparm(Vbase=tvV*exp(etaV))
V=Vbase*(1+coefwt2*(weight-75))

where now, coefwt2=coefwt/tvV.

An example of this splitting technique is included in the ‘ready-to-go’ project Remifen.phxproj in the examples directory. The stparm block of the original model in ‘ELS_ FOCE_ Formulation’ in this project has the form:

stparm(
  V1TV=theta1-theta7*(AGE-40)+theta12*(LBM-55)
  V2TV=theta2-theta8*(AGE-40)+theta13*(LBM-55)
   V3TV=theta3
  CL1TV=theta4-theta9*(AGE-40)+theta14*(LBM-55)
  CL2TV=theta5-theta10*(AGE-40)
  CL3TV=theta6-theta11*(AGE-40)
  V1=V1TV*exp(eta1)
  V2=V2TV*exp(eta2)
  V3=V3TV*exp(eta3)
  CL1=CL1TV*exp(eta4)
  CL2=CL2TV*exp(eta5)
  CL3=CL3TV*exp(eta6)
  K10=CL1/V1
  K12=CL2/V1
  K13=CL3/V1
  K21=CL2/V2
  K31=CL3/V3
)

This will run correctly with the FOCE ELS engine, but not QRPEM since clearly this has multiple instances of type of nonlinear covariate model of the type discussed in g) above. Another difficulty is the inclusion of the assignment statements defining the rate constants K10, K12, K13, K21, and K31, which are prohibited in stparm statements in correctly formulated QRPEM models.

The equivalent QRPEM-compatible ‘split version’ is:

stparm(
  V1qr= theta1*exp(eta1)
  V2qr= theta2*exp(eta2)
  V3=   theta3*exp(eta3)
  CL1qr=theta4*exp(eta4)
  CL2qr=theta5*exp(eta5)
  CL3qr=theta6*exp(eta6)
)
  V1=V1qr*(1-theta7*(AGE-40)+theta12*(LBM-55))
  V2=V2qr*(1-theta8*(AGE-40)+theta13*(LBM-55))
  CL1=CL1qr*(1-theta9*(AGE-40)+theta14*(LBM-55))
  CL2=CL2qr*(1-theta10*(AGE-40))
  CL3=CL3qr*(1-theta11*(AGE-40))
  K10=CL1/V1
  K12=CL2/V1
  K13=CL3/V1
  K21=CL2/V2
  K31=CL3/V3

The second main type of example where splitting is useful is time varying covariates where not all instances of a covariate within a given individual have the same value. This violates some of the basic assumptions that allow simple linear regression based updates even if the covariate model is linear. Consider for example the linear covariate model in d) above:

stparm(V=exp(tvlogv+coefwt*(weight-75)+etaV)

If the covariate weight has several different values within a single individual, then attempting to run this otherwise acceptable model with QRPEM will fail with an error message describing the specific offending time varying covariate.

However, the same simple basic splitting technique where the time-varying covariate is placed outside the stparm statement allows it to run in QRPEM:

stparm(Vbase=exp(tvlogv+etaV)
V=Vbase*exp(coefwt*(weight-75))

An example of this is given in the example project QRPEM_with_time_varying_covariate.phxproj. Here the basic model is a Cl/V parameterization of a simple one-compartment IV-bolus case, where log(Cl) has a linear dependency on a covariate ‘scr’:

test(){
    cfMicro(A1, Cl/V)
    dosepoint(A1)
    C=A1/V
    error(CEps=1)
    observe(CObs=C*(1+CEps))
    stparm(V=tvV*exp(nV))
    stparm(Cl=tvCl*exp(scr*dCldscr+nCl))
    fcovariate(scr)
    fixef(tvV=c(, 2,))
    fixef(tvCl=c(, 0.5,))
    fixef(dCldscr(enable=c(0))=c(, 1,))
    ranef(diag(nV, nCl)=c(1, 1))
}

The model as written is acceptable for QRPEM if scr is not time-varying but rather has a fixed value for each individual. The QRPEM_with_fixed_cov and ELSFOCE_with_fixed_cov models in this project are set up to run QRPEM and FOCE ELS, respectively, on exactly this model for such a dataset with fixed scr covariate. A second data set with a time-varying scr is also provided. Here exactly the same model as above will run FOCE ELS correctly (model ELSFOCE_with_var_cov), but the stparm for Cl must be split for QRPEM to work in this case (workflow (QRPEM_with_var_cov):

test(){
    cfMicro(A1, Cl/V)
    dosepoint(A1)
    C=A1/V
    error(CEps=1)
    observe(CObs=C*(1+CEps))
    stparm(V=tvV*exp(nV))
    stparm(Clbase=tvCl*exp(nCl))
    #following moves time-varying covariate out of stparm statement
    #and now model is acceptable for QRPEM
    Cl=Clbase*exp(scr*dCldscr)
    fcovariate(scr)
    fixef(tvV=c(, 1,))
    fixef(tvCl=c(, 0.5,))
    fixef(dCldscr(c(, 1,))
    ranef(diag(nV, nCl)=c(1, 1))
}

Finally, note that the split version can be run with the fixed covariate dataset, but it is more efficient to use the un-split version when feasible. The splitting in effect forces the use of the less efficient log likelihood optimization method to update the fixed effect dCldscr in the Cl covariate model, whereas this can be updated with the more efficient linear regression technique in the un-split case.