Re: RES and WRES output with Beal's M3 method
Thanks Martin, Matt and Scott for your replies, most helpful.
As a side note, my colleagues and I tried to look more into a "manual" way to compute WRES, which appears to be a non trivial task and involves the computation of partial derivatives with respect to the IIV parameters. The detailed equations in NONMEM guide 1 page 41 clearly explained why the computations of WRES simplify to RES/W when no IIV is included in the model.
Sebastien
Matt Hutmacher wrote:
> Hello all,
>
> IWRES, WRES, CWRES, or even Monte Carlo-based population residuals do not
> provide great diagnostic value unfortunately for datasets with censored
> data. BQL observations influence the fit through the censored likelihood,
> but these observations are not represented in the residual diagnostic plots
> (they are not defined for the BQL observations). In fact, the population
> residual plots should not have a mean of 0, nor a variance of 1 (if
> standardized by variance estimates), and are likely to have skewness because
> of the conditional nature of their calculation. The degree of departure
> from these typical target values depends upon the number of BQL observations
> of course. Also, these residuals will look worse using M3 method (a more
> principled approach) compared to those derived from fitting using method M1,
> which discards BQL observations during estimation. This reflects the bias
> in the estimates and predictions induced by excluding the non-ignorable
> censored observations (M1), counter to that typically expected from
> inspection of residuals. Predictive checks, as Martin suggests, are
> probably the only tool for evaluating the model when you have an influential
> number of BQL observations. Martin provided a very nice poster on VPCs for
> such situations. However, if you wanted to use residuals in a PPC as
> opposed to concentration (perhaps a less sufficient statistic), then which
> population residual you compute, WRES, CWRES, or Monte Carlo population
> residuals might not be all that influential in the model evaluation as long
> as you compute the residuals similarly for the original model fit to the
> observed data and the model fits to the simulated datasets.
>
> One other thing to consider is outliers. Residual-based determination of
> outliers cannot be applied to the observations that are BQL (perhaps the
> threshold should be determined from simulations given the expected lack of
> normality when censoring is present). However, I would argue that this
> doesn't mean that a BQL observation cannot be an outlier. Take for example
> a concentration profile that has a BQL observation reported around the time
> of TMAX and also that the concentrations before and after the BQL
> observation are a reasonable distance (perhaps 4-6 standard deviations away)
> from the lower limit of quantification. The likelihood for this BQL
> observation is PHI((QL-IPRED)/W) where PHI is the cumulative normal, and
> IPRED and W are defined as per Martin below. If IPRED is large relative to
> QL, this provides a very negative value of (QL-IPRED)/W, and PHI() will be
> near 0. When -2 X log-likelihood is computed, i.e., -2 X log(PHI()), a very
> large number will result. For example, if (QL-IPRED)/W = -6, then -2 X
> log(PHI(-6)) = 41.5. If we looked at an observation that was not censored
> that had an IWRES = -6 (assume approximate %CV = 30%), the -2 X log
> likelihood is 33.6. Thus, the QL-IPRED/W might be considered an outlier,
> and the estimation procedure might influenced by this BQL observation.
> Granted the scenario here suggests that the point might be an outlier
> without 'quantifying' its degree. Note that if (QL-IPRED)/W = 6, this
> suggest that IPRED is much less than QL. This is likely to occur in the
> elimination phase after a few other BQL observations have been observed. In
> this case, it is not an issue, because PHI(6) is near 1, and -2 X
> log(PHI(6)) is near 0 indicating a negligible contribution to the
> likelihood. In my mind, this is along the line of reasoning for proposing
> the M6 method, where QL/2 is used to impute the first BQL observation in the
> terminal phase. I am less certain about the influence of (QL-IPRED)/W = 6
> in the absorption phase, especially if Tlag is estimated.
>
> Best,
> Matt
>
Quoted reply history
> -----Original Message-----
> From: [email protected] [mailto:[email protected]] On
> Behalf Of Martin Bergstrand
> Sent: Wednesday, April 07, 2010 1:02 PM
> To: 'Sebastien Bihorel'; [email protected]
> Subject: RE: [NMusers] RES and WRES output with Beal's M3 method
>
> Dear Sebastien and NONMEM users,
>
> It is correct that NONMEM doesn't return any RES or WRES for individuals who
> have at least one observation with F_FLAG=1. It is understandable that
> NONMEM can't return residuals for the BQL observations (F_FLAG=1) since PRED
> is a likelihood in this case and not a prediction (furthermore the
> observation is an interval and can neither be used to calculate a residual).
> However, I don't understand why it doesn't do so for the observations for
> which F_FLAG=0? Can this be considered a bug?
>
> You can of course always get out individual residuals (IRES) and individual
> weighted residuals (IWRES) (below this paragraph you can see my definition
> of IWRES). What you can do if you really want to get RES and WRES is to run
> a MAXEVAL=0 run with your final estimates. The BQL data should in this case
> be omitted (or set to a fix value) and all M3 code taken out. Remember in
> this case that the IRES and IWRES and other EBE dependent variables will not
> be correct following the MAXEVAL=0 run (i.e. use only the PRED, RES and WRES
> out of this run and the rest from your original fit with the M3 code). In
> case you are using the FOCE estimation method you would preferably want to
> look at conditional weighted residuals (CWRES) rather than WRES. However I
> have yet not seen any example of a workaround to get these when using the M3
> method (Obs! CWRES will not be correctly calculated with the MAXEVAL=0
>
> trick).
>
> $ERROR
> IPRED = F
> W = THETA(11) ; SD for additive
>
> residual error ; W = THETA(11) * IPRED ; SD for proportional
>
> residual error
> ; W = SQRT(THEATA(11)**2+THETA(12)**2*IPRED**2) ; SD for combined
> residual error
> Y = IPRED + W*EPS(1)
> IRES = DV - IPRED
> IWRES = IRES/W
>
> $SIGMA 1 FIX
>
> If you use Xpose to do your diagnostic plots it can be good to know that
> Xpose by default omits all rows with WRES=0 (this is done to omit
> information on dosing row from being plotted). You can change this setting
> in Xpose by using the command inclZeroWRES=TRUE (you should then enter some
>
> other subset to omit the dose rows e.g. subset="EVID==0").
>
> My favorite type of diagnostics are VPCs (I think I share this with Prof.
> Nick Holford). VPCs can easily be adopted to handle the presence of BQL
> data. How to do this with PsN and Xpose is explained in a PAGE poster from
> 2009
> ( http://www.page-meeting.org/pdf_assets/7002-Poster_PAGE_VPC_090618_final.pd
> f).
>
> Best regards,
>
> Martin Bergstrand, MSc, PhD student
> -----------------------------------------------
> Pharmacometrics Research Group,
> Department of Pharmaceutical Biosciences,
> Uppsala University
> -----------------------------------------------
> [email protected]
> -----------------------------------------------
> Work: +46 18 471 4639
> Mobile: +46 709 994 396
> Fax: +46 18 471 4003
>
> -----Original Message-----
> From: [email protected] [mailto:[email protected]] On
> Behalf Of Sebastien Bihorel
> Sent: den 7 april 2010 17:31
> To: [email protected]
> Subject: [NMusers] RES and WRES output with Beal's M3 method
>
> Dear NONMEM users,
>
> We have noticed some problems with the table outputs of RES and WRES, when we implemented Beal's M3 method as proposed by Ahn and colleagues (J Pharmacokinet Pharmacodyn (2008) 35:401-421). While RES and WRES are correctly calculated and reported for patients without BLOQ records, these metrics are reported as being 0 for patients with at least one BLOQ sample, even for the records that were not flagged as BLOQ. This behavior seems to be common to NONMEM 6 and 7.
>
> Does anybody know about an NONMEM option or a workaround that would allow the user to access to the actual RES and WRES for the non-BLQ records?
>
> Any feedback would be greatly appreciated.
>
> Sebastien Bihorel
>
> PS: this is the $ERROR block code we used for a simple proportional RV model
>
> $ERROR (ONLY OBSERVATIONS)
>
> ;Information needed for BLQF and > BLQF samples
> LOQ=5
> SIG = THETA(3)
>
> DFLG=0 ;create a dose record flag
> IF(AMT.NE.0) DFLG=1
>
> IPRED=F+DFLG
> W=IPRED
>
> ;Computations for samples with DV > LOQ (BLQF=0)
> IF (BLQF.LT.0.25) THEN
> F_FLAG=0
> FFLG=0
> IRES=DV-IPRED
> IWRES=IRES/W
> Y=IPRED+W*SIG*EPS(1) ;NOTE: Prediction is a concentration
> ENDIF
>
> ;Computations for samples with DV <= LOQ (BLQF=1)
> IF (BLQF.GE.0.25) THEN
> F_FLAG=1
> FFLG=1
> DUM=(LOQ-IPRED)/(W*SIG)
> IPRED=PHI(DUM)
> Y=IPRED ;NOTE: prediction is a *probability*
> ENDIF