Date: Sat, 05 Aug 2000 09:21:59 +1200

From: Nick Holford <n.holford@auckland.ac.nz>

Subject: [Fwd: CLIN PHAR STAT: Mixed Vs Fixed]

Stephen Senn contributed these comments to the cps@topica.com list. I thought they might be of interest to NMUSERs. Would anyone like to comment on the biased SE issue as it might apply to NONMEM? Is it possible to put "a random effect on the variances themselves" in NONMEM?

stephens@public-health.ucl.ac.uk wrote:

>

> Of course, the mixed (random effects) analysis makes an

> assumption of Normality that the fixed effects model does not

> make. There is a further technical problem. Ideally, the between

> and within patient information should be weighted by the true

> variances. In practice they are weighted by the observed variances

> and this leads to a downward bias in the estimated overall standard

> error and a consequent bias towards significance in the p-values.

> (A similar problem occurs in meta-analysis whether fixed or

> random.) This can be corrected by putting a random effect on the

> variances themselves but I would guess (I don't know) that this is

> not incorporated currently in PROC MIXED in SAS.

>

> Regards

>

> Stephen

--

Nick Holford, Dept Pharmacology & Clinical Pharmacology

University of Auckland, Private Bag 92019, Auckland, New Zealand

email:n.holford@auckland.ac.nz tel:+64(9)373-7599x6730 fax:373-7556

http://www.phm.auckland.ac.nz/Staff/NHolford/nholford.htm

*****

Date: Tue, 08 Aug 2000 07:34:04 +1200

From: Nick Holford <n.holford@auckland.ac.nz>

Subject: Re: [Fwd: CLIN PHAR STAT: Mixed Vs Fixed]

Leonid,

Thanks for this suggestion but it doesn't seem to me that THETA(10) is identifiable. I think you need some additional covariate to distinguish the components of variance for CL when using this kind of approach.

I was wondering if the undocumented Bayesian PRIOR feature of NONMEM would be appropriate here. By putting a prior on the variance of CL this is one way of putting "a random effect on the variances themselves". (Mats, Diane -- are you there? can you comment?).

Nick

"Gibiansky, Leonid" wrote:

>

> >

> > Is it possible to put "a random effect on the variances themselves" in

> NONMEM?

> >

>

> One can try

>

> VAR=THETA(10)*EXP(ETA(10)) ;assuming that the other 9 THETAs and ETAs are

> used elsewhere

>

> ET1=VAR*ETA(1)

> CL=THETA(1)*EXP(ET1)

>

> $OMEGA

> 1 FIXED ; variance of eta1

> ...

> 0.5 ; variance of random effect for variance of eta1

>

> I am not sure about convergence and whether you gain anything...

> Leonid

>

--

Nick Holford, Dept Pharmacology & Clinical Pharmacology

University of Auckland, Private Bag 92019, Auckland, New Zealand

email:n.holford@auckland.ac.nz tel:+64(9)373-7599x6730 fax:373-7556

http://www.phm.auckland.ac.nz/Staff/NHolford/nholford.htm

*****

Date: Mon, 7 Aug 2000 22:02:06 +0100 (GMT)

From: "J.G. Wright" <J.G.Wright@newcastle.ac.uk>

Subject: Re: [Fwd: CLIN PHAR STAT: Mixed Vs Fixed]

Dear Nick,

Basically, a hierarchical model has more freedom than it admits too when calculating it SEs. For fixed effects models,degrees of freedom are easy to work out. In a hierarchical random effects model, its not so clear.

Because the variance components are estimated and determine their own weighting, their hidden degrees of freedom lead to a geenral increase in freedom. The model can shuffle the variance components around to find the cheapest way to "pay" for discrpancies between data and predictions. This problem is greatly amplified by the use of joint maximum likelihood in NONMEM.

It is possible to write abbreviated code, to put a random effect on the magnitude of a random effect. The idea is then that the model will integrate over this parameters possible values rather than treating it as fixed. However, this additional random effect is another parameter to be estimated and doesn't have a natural interpretation for OMEGA in in a single population, so I guess what we basically want is a hyperprior onn a variance component. I should stress that a random effect scaling the magnitude of a further random effect is entering the likelihood integral in a highly nonlinear manner and is unlikely to be well-estimated under a linearization of this integral. However, I suppose by fixing it we could obtain such a hyperprior.

NONMEM will handle this kind of effect differently at the first and secons stage of the hierarchy - I am not quite sure which of these will "solve the problem", and indeed whether introducing an added parameter will create additional difficulties, especially in a linearised joint maximum likelihood model. I would guess (literally) that the problem Stephen Senn refers to occurs using GLS fitting procedures and, although similar considerations almost certainly apply, NONMEM will calculate its SEs from likelihood curvature. Both Mats Karlsson and Stuart Beal have done some work on random effects affecting residual error...(Cue)

Presumably, the severity of this problem decreases as the number of patients increases - anyway the properties of SEs for hypothesis testing are asymptotic by nature in nonlinear models. I think similar considerations may also apply to the degrees of freedom of the chi-squared distribution for the change in the objective function.

I hope I have understood this correctly,

James Wright

*****

Date: Mon, 07 Aug 2000 23:35:18 +0200

From: Mats Karlsson <Mats.Karlsson@biof.uu.se>

Subject: Re: [Fwd: CLIN PHAR STAT: Mixed Vs Fixed]

Dear James and Nick,

James alluded to a model for interaction between ETA's and EPS's. What he refers to is to put an interindividual random effect on the residual error. The reason you may want to do this is to acknowledge that residual error magnitude may vary between subjects (more or less model misspecification, better or worse compliers, etc). The way to do that in NONMEM is Y=F+EPS(1)*EXP(ETA(1)) where the magnitude of OM(1) would reflect the degree to which subjects differ in residual error magnitude. To use this option, you have to use FOCE INTERACTION. As opposed to interindividual variability in the ETA, which as James points out is tricky to estimate and relatively hard to grasp, interindividual variability in maginutde of residual error is often well estimated if you have a decent number of observations per subject (say >4). It can be a useful diagnostic (high POSTHOC residual variabiltiy in an individual-> poor fit) and it makes subjects whose profiles are not at all well described by the model to influence parameter estimates to a lower degree. On the topic of misspecification, I assume it would be possible also to estimate interindividual varaibility in auto-correlation between residuals, again providing a measure for varying degree of misspecification. I've never tried it or heard of anyone who has though.

Best regards,

Mats

--

Mats Karlsson, PhD

Professor of Biopharmaceutics and Pharmacokinetics

Div. of Biopharmaceutics and Pharmacokinetics

Dept of Pharmacy

Faculty of Pharmacy

Uppsala University

Box 580

SE-751 23 Uppsala

Sweden

phone +46 18 471 4105

fax +46 18 471 4003

mats.karlsson@biof.uu.se

*****

From: Stephen Senn <stephens@public-health.ucl.ac.uk>

Date: Tue, 8 Aug 2000 11:59:19 GMT0

Subject: Re: [Fwd: CLIN PHAR STAT: Mixed Vs Fixed]

Dear Nick, James and Mats

Mats is right in saying

> > Presumably, the severity of this problem decreases as the number of

> > patients increases -

in the context of meta-analysis, because here the patients are the repeated measurements on the trial effects. What we need to know is the correct variance from each trial and as the patients increase we get closer to estimating the true within-trial variance. I have a letter on this subject which will be appearing in Controlled Clinical Trials which shows that David Salsburg's recent contribution to that journal underestimates this problem and that he is wrong to recommend meta-analysis as an alternative to the linear model.

I am very ignorant on PK/PD but the analogy here would seem to be not in the number of patients but in the number of measurements per patient. In this context, there may be a bias in variance estimation and associated inferential statistics (CI, P-values (ugh!) etc) for sparse sampling. However, it depends on the way you set the model up. If you impose a common residual variance for each patient then the problem largely disappears. (THis is analogous to what is done in a fixed effects linear model analysis of a multi-centre trial.)

Of course, fully Bayesian methods put a prior on everything and so deal with this problem. (Or at least appear to deal with it.)

Regards

Stephen

Professor Stephen Senn

Department of Statistical Science &

Department of Epidemiology and Public Health

University College London

Room 126, 1-19 Torrington Place

LONDON WC1E 6BT

Tel: +44 (0) 20 7679 1698

Fax: +44 (0) 20 7383 4703

Email: stephens@public-health.ucl.ac.uk

webpage: http://www.ucl.ac.uk/~ucaksjs/

*****

Date: Wed, 09 Aug 2000 08:14:25 +1200

From: Nick Holford <n.holford@auckland.ac.nz>

Subject: Re: [Fwd: CLIN PHAR STAT: Mixed Vs Fixed]

Stephen,

Stephen Senn wrote:

>

> I am very ignorant on PK/PD but the analogy here would seem to

> be not in the number of patients but in the number of

> measurements per patient. In this context, there may be a bias in

> variance estimation and associated inferential statistics (CI, P-

> values (ugh!) etc) for sparse sampling. However, it depends on the

> way you set the model up. If you impose a common residual

> variance for each patient then the problem largely disappears.

It seemed to me that Mats had suggested a solution to this problem which was to explicitly use a DIFFERENT residual variance for each patient. The usual NONMEM model assumes a common residual error for each patient. What do you think the consequences of using different residual error for each patient might be?

> Of course, fully Bayesian methods put a prior on everything and so

> deal with this problem. (Or at least appear to deal with it.)

I had suggested using the NONMEM Bayesian prior approach earlier in this thread. NONMEM allows one to put a normal (or lognormal) prior on structural model parameters (THETA), an inverse Wishart prior on random effects parameters (OMEGA) and normal/lognormal prior on residual error parameters (SIGMA) (by indirection using THETA). How well does this approach your definition of "fully" Bayesian?

[Please reply to nmusers <nmusers@c255.ucsf.edu>]

--

Nick Holford, Division of Pharmacology & Clinical Pharmacology

University of Auckland, Private Bag 92019, 85 Park Road, Auckland, NZ

email: n.holford@auckland.ac.nz tel:+64(9)373-7599x6730 fax:373-7556

http://www.phm.auckland.ac.nz/Staff/NHolford/nholford.htm

*****

From: "Stephen Duffull" <sduffull@pharmacy.uq.edu.au>

Subject: RE: [Fwd: CLIN PHAR STAT: Mixed Vs Fixed]

Date: Wed, 9 Aug 2000 09:37:35 +1000

Nick

Here is a simple example of the prior setup for the WinBUGS code for implementing a Bayesian PK analysis:

...the model bit goes here...

z[j,i] ~ dnorm(c[j,i],tau[j]) # likelihood function

}

# Individual priors

theta[j,1:2] ~

dmnorm(mu.theta[1:2],omega.theta[1:2,1:2])

tau[j] ~ dgamma(0.001,0.001)

}

# Population Priors

mu.theta[1] ~ dnorm(.05,0.1)#I(0,)

mu.theta[2] ~ dnorm(1,0.1)#I(0,)

omega.theta[1:2,1:2] ~ dwish(R[1:2,1:2],2)

In this example z[j,i] is the ith observed concentration for the jth individual, and c[j,i] is the expected concentration; there are 2 parameters (thetas say CL and V) in the model. I have used mu.theta[x] to represent the population value (sometimes in NONMEM TV is used) of the parameter. Note that tau[j] is the precision of the residual error (precision = 1/variance] and there is one for each individual. It is of course possible to put a population prior on tau[j] to allow the individual estimates of RUV to be conditional on the population.

Nick wrote:

> I had suggested using the NONMEM Bayesian prior

> approach earlier in this

> thread. NONMEM allows one to put a normal (or

> lognormal) prior on

> structural model parameters (THETA), an inverse

> Wishart prior on random

> effects parameters (OMEGA) and normal/lognormal

> prior on residual error

> parameters (SIGMA) (by indirection using THETA).

> How well does this

> approach your definition of "fully" Bayesian?

The difficulty I have lies in not knowing exactly what NONMEM does when it does the NONMEM Bayesian approach (please note that I have never used this feature of NONMEM and indeed don't know how too). In contrast WinBUGS is quite transparent about both its hierarchical structure and about methodology it uses to get around nasty integration problems.

IMHO I think that if a Bayesian method is desirable then a Bayesian program should be used - but I may be biased :-).

Regards

Steve

=================

Stephen Duffull

School of Pharmacy

University of Queensland

Brisbane, QLD 4072

Australia

Ph +61 7 3365 8808

Fax +61 7 3365 1688

http://www.uq.edu.au/pharmacy/duffull.htm

*****

Date: Wed, 09 Aug 2000 12:41:14 +1200

From: Nick Holford <n.holford@auckland.ac.nz>

Subject: Re: [Fwd: CLIN PHAR STAT: Mixed Vs Fixed]

Steve,

Stephen Duffull wrote:

> Here is a simple example of the prior setup for the WinBUGS

> code for implementing a Bayesian PK analysis.

This is the nmusers list. I have deleted that WinBuggy filth :-)

> The difficulty I have lies in not knowing exactly what

> NONMEM does when it does the NONMEM Bayesian approach

> (please note that I have never used this feature of NONMEM

> and indeed don't know how too).

Here is a "simple" example of how to use NONMEM with priors on THETA and OMEGA. It requires a user defined prior subroutine which follows the control stream:

$PROB theophylline pharmacodynamics

$DATA theopd.dat IGNORE #

$INPUT ID TIME THEO AGE WT GEND RACE DIAG PEFR=DV

$ESTIM PRINT=1 POSTHOC

MSFO=theopd.msf

$COV

$THETA (0,146.,) ; E0

$THETA (0,164.,) ; EMAX

$THETA (.001,6.55,) ; EC50

$OMEGA 0.00178 ; CVE0

$OMEGA 0.239 ; CVEMAX

$OMEGA 1.46 ; CVEC50

$SIGMA 6590. ; SD 1

$SUB PRIOR=theopdB.for

;Prior normal mode for THETA

$THETA 146. FIX ; theta1

$THETA 164. FIX ; theta2

$THETA 6.55 FIX ; theta3

;Prior uncertainty on normal mode for THETA

$OMEGA BLOCK(3) FIX

122

-125 427

6.05 20.1 3.82

;Prior OMEGA

$OMEGA 0.00178 FIX ; omega1

$OMEGA 0.239 FIX ; omega2

$OMEGA 1.46 FIX ; omega3

;Degrees of Freedom for OMEGA (Nsubjects-1)

$THETA 152 FIX ; dfomega1

$THETA 152 FIX ; dfomega2

$THETA 152 FIX ; dfomega3

$PRED

e0= THETA(1)*EXP(ETA(1))

EMAX=THETA(2)*EXP(ETA(2))

EC50=THETA(3)*EXP(ETA(3))

Y = E0 + EMAX*THEO/(THEO+EC50) + ERR(1)

$TABLE ID TIME THEO AGE WT GEND RACE DIAG

E0 EMAX EC50 Y

NOPRINT ONEHEADER FILE=theopd.fit

theopdB.for

SUBROUTINE PRIOR (I,CNT,NTHP,NETP,NEPP)

DOUBLE PRECISION CNT

IF (ICALL.LE.1) THEN

NTHP=3

NETP=3

NEPP=1

ENDIF

CALL NWPRI(CNT)

RETURN

END

> In contrast WinBUGS is

> quite transparent about both its hierarchical structure and

> about methodology it uses to get around nasty integration

> problems.

Its all clear as mud.

> IMHO I think that if a Bayesian method is desirable then a

> Bayesian program should be used - but I may be biased :-).

Not biased -- just different set of priors.

--

Nick Holford, Division of Pharmacology & Clinical Pharmacology

University of Auckland, Private Bag 92019, 85 Park Road, Auckland, NZ

email: n.holford@auckland.ac.nz tel:+64(9)373-7599x6730 fax:373-7556

http://www.phm.auckland.ac.nz/Staff/NHolford/nholford.htm

*****

Date: Tue, 08 Aug 2000 17:48:43 -0700

From: LSheiner <lewis@c255.ucsf.edu>

Subject: Re: [Fwd: CLIN PHAR STAT: Mixed Vs Fixed]

Nick Holford wrote:

>

> Stephen,

>

> Stephen Senn wrote:

> >

> > I am very ignorant on PK/PD but the analogy here would seem to

> > be not in the number of patients but in the number of

> > measurements per patient. In this context, there may be a bias in

> > variance estimation and associated inferential statistics (CI, P-

> > values (ugh!) etc) for sparse sampling. However, it depends on the

> > way you set the model up. If you impose a common residual

> > variance for each patient then the problem largely disappears.

>

> It seemed to me that Mats had suggested a solution to this problem which

> was to explicitly use a DIFFERENT residual variance for each patient.

> The usual NONMEM model assumes a common residual error for each patient.

> What do you think the consequences of using different residual error for

> each patient might be?

>

What Mats suggests is to use a random effect on the residual eror - this means that you think of the error variance as being different for each person, but just like individual random effects ("etas" to us), this does not add a degree of freedom for each indivdiual, but only one df in total, for the estimated variance (of the variances).

LBS.

--

_/ _/ _/_/ _/_/_/ _/_/_/ Lewis B Sheiner, MD (lewis@c255.ucsf.edu)

_/ _/ _/ _/_ _/_/ Professor: Lab. Med., Bioph. Sci., Med.

_/ _/ _/ _/ _/ Box 0626, UCSF, SF, CA, 94143-0626

_/_/ _/_/ _/_/_/ _/ 415-476-1965 (v), 415-476-2796 (fax)

*****

From: Stephen Senn <stephens@public-health.ucl.ac.uk>

Date: Wed, 9 Aug 2000 11:46:26 GMT0

Subject: Re: [Fwd: CLIN PHAR STAT: Mixed Vs Fixed]

Nick,

It is interesting to compare meta-analysis to the general linear model

A. In a general linear model analysis of a data set consisting of two treatments and many strata you (usually) impose the assumption that the within-cell variances are constant. This then weights the contribution from the respective strata according to the design information: it depends only on the number of observations and their distribution between treatment and control. The danger is that if heteroscedasticty applies this will be inefficient.

B. In a meta-analysis you use the calculated (and hence estimated) variances for each within-stratum estimate as the means of combining them. The danger is that differences in the weights reflect random differences rather than true differences.

Because we often have little information in clinical work the dangers inherent in B are often more serious than those in A. Hence a constrained common variance is often preferable. Typically, frequentist solutions end up at one of two extremes. A Bayesian approach permits compromises between these two.

Now, turning to PK/PD work. Here you have non-linear models. This means that evn if you constrain a residual variance to be equal the estimated variances of parameters will not be fixed by the design alone. It depends on the response as well. The fact that you have a hierarchy on the parameters may help with this problem but I suspect that in any case it means that approach A veres towards approach B (but I am guessing). However, I would also guess that freeing the variances will make this problem worse rather than better.

Regards

Stephen