# RE: omega matrix - friendly suggestion

From: Stephen Duffull <stephen.duffull>
Date: Thu, 14 Jun 2012 18:51:00 +0000

It turns out that there are many symmetric matrices where off-diagonals are=
between -1 and 1 that have negative determinants. Just to add my 2c there =
are also non-positive definite symmetric matrices that have positive determ=
inants. Some are obvious and some are not. So if you were to parameterise=
using SD and corr you would need some care. I am pleased I don't have to c=
ode search algorithms that require positive definite matrices as a search s=
tep :-)

Sampling algorithms don't have this issue as there are good priors that gua=
rantee "legal" matrices.

Steve

> -----Original Message-----
> From: owner-nmusers
On
> Behalf Of Bob Leary
> Sent: Friday, 15 June 2012 5:11 a.m.
> To: Bauer, Robert; Matt Hutmacher; 'Nick Holford'; nmusers
> Subject: RE: [NMusers] omega matrix - friendly suggestion
>
> Easy countereaxample to Nick's hypothesis that a covariance matrix with a=
n
> underlying correlation matrix with all elements of magnitude <=1 is pos=
itive
> definite
>
> Let all sd's be 1, and let the correlation matrix be c =
> 1.0000 0.9900 -0.9900
> 0.9900 1.0000 0.9900
> -0.9900 0.9900 1.0000
>
> This is in fact an impossible correlation matrix, since it is not poitive
> definite:
>
> > det(c)
> ans =
> -3.8809
>
> (THe determinant of a positive definite matrix must be positive)
>
>
> The two most common ways to parameterize a symmetric pos def matrix that
> guarantee positive definiteness are the cholesly factor method, and the
> matrix exponential of arbitrary symmetric matrix . Neither one can be use=
d to
> force an arbitrary element of a full block matrix to a given value while
> preserving postive definiteness.
> _______________________________________
> From: owner-nmusers
lf Of
> Bauer, Robert [Robert.Bauer
> Sent: Wednesday, June 13, 2012 4:05 PM
> To: Matt Hutmacher; 'Nick Holford'; nmusers
> Subject: RE: [NMusers] omega matrix - friendly suggestion
>
> Thank you all for your suggestions. The more flexible OMEGA element entr=
ies
> will not be easy to implement in NONMEM's structure, but I will look into=
it.
> A couple of things need to be noted.
>
> 1) Regarding Matt's suggestions, to enter initial Omega elements in Chole=
sky
> format in NM 7.2, \$OMEGA CHOLESKY may be used.
> 2) The traditional methods (FO/FOCE/Laplace) actually vary the omega elem=
ents
> in the Cholesky format for estimation, than transform back in the varianc=
e
> format after each iteration. This assures that as NONMEM arbitrarily sea=
rches
> throughout the Cholesky space (that is, as each cholesky element is
> arbitrarily varied in search of the minimum OFV), the variance omega is
> guaranteed to be positive definite. Thus, initially fixing an off-diagon=
al
> variance element to a non-zero value within a block of variances cannot b=
e
> honored by the search algorithm past the first iteration. Additionally, =
to my
> best understanding, arbitrarily fixing an off-diagonal variance element t=
o 0
> does not correspond to fixing a particular Cholesky off-diagonal element =
to 0,
> except for banded matrix patterns, which is already allowed in NONMEM.
> 3) EM algorithms such as Imp and SAEM update the Omega variance in the
> variance domain, and also, allow a simple constraint filter to be super-
> imposed on the resulting matrix after each iteration. In NONMEM, this ca=
n be
> done using constraint.f90, wherein a person may code in fixing any omega
> element to 0, or any other value, over-riding the EM algorithm's update. =
This
> is okay in EM, as the omega elements are not primary search parameters, b=
ut
> are constructed as a consequence of the expectation step, which provide
> conditional means and variances from which the omega elements are evaluat=
ed.
> User beware, should this constraint cause lack of positive definiteness. =
The
> constraint.f90 routine is not used for the traditional update methods, as=
this
> would upset their gradient search pattern by providing a dis-continuous
> intrusion.
>
>
>
> Robert J. Bauer, Ph.D.
> Vice President, Pharmacometrics, R&D
> ICON Development Solutions
> 7740 Milestone Parkway
> Suite 150
> Hanover, MD 21076
> Tel: (215) 616-6428
> Mob: (925) 286-0769
> Email: Robert.Bauer
> Web: www.iconplc.com
>
> -----Original Message-----
> From: owner-nmusers
On
> Behalf Of Matt Hutmacher
> Sent: Wednesday, June 13, 2012 2:35 PM
> To: 'Nick Holford'; nmusers
> Subject: RE: [NMusers] omega matrix - friendly suggestion
>
> Dear Paolo, Nick, all,
>
> A way to ensure a nonnegative definite (NND) variance-covariance matrix i=
s to
> use a Cholesky decomposition of the matrix (Ref 1). I use this often and=
it
> greatly stabilizes estimation such that I can often get better convergenc=
e and
> improved \$COV step success for larger OMEGA matrix. It also is great if =
you
> sample from the \$COV for simulation as an NND \$OMEGA will be achieved wit=
h
> each sampling. I currently implement this by hand such that all \$OMEGA
> becomes a Diagonal matrix with 1 FIXED as all of its entries.
> Note that I have not evaluated whether NONMEM computes the appropriate WR=
ES or
> CWRES for this parameterization (as I have calculated these myself - Bob =
may
> be able to comment on whether this is a concern - ie putting the variance=
s
> into linear combinations of thetas affecting the residuals). It is not en=
ough
> to just model variances and correlations; such a parameterization will no=
t
> guarantee NND OMEGA matrices.
>
> With respect to the original question which was to constrain certain
> covarainces to 0, one might (I have not tried this) be able to do this by
> writing the covariance matrix as wished, then solving for the Cholesky
> components, and modeling with the derived Cholesky parameterization. Upo=
n
> completion, one can then compute the \$OMEGA matrix from the Cholesky
> parameterization. I do not know how easy this would be to generalize in =
a
> software application, but if possible, this could perhaps be done behind =
the
> scenes as an option (may make standard errors a bit tricky) and would pro=
vide
> a lot of flexibility for analysis (as opposed to hand coding discussed be=
low).
>
> There is also a modified Cholesky that can be used, which has the nice fe=
ature
> that one can simply, conveniently remove all the variance components rela=
ted
> to a certain random effect (Ref 2). This was recently used in the attach=
ed
> reference to aid in selection of the \$OMEGA structure with the fixed effe=
cts.
>
> Best, Matt
>
>
> Ref.1: Pinhiero, J. and Bates, D. (1996). Unconstrained parameterizations=
for
> variance-covariance matrices. Statistics and Computing 6, 289-286.
>
> Ref.2: Bondell HD, Krishna A, Ghosh (2010). Joint Variable Selection for =
Fixed
> and Random Effects in Linear Mixed-Effects Model. Biometrics 66, 1069-10=
77.
>
>
> -----Original Message-----
> From: owner-nmusers
On
> Behalf Of Nick Holford
> Sent: Wednesday, June 13, 2012 1:20 PM
> To: nmusers
> Subject: Re: [NMusers] omega matrix - friendly suggestion
>
> Paolo,
>
> You may remember this discussion about the history of the correlation fea=
ture
> for OMEGA (and SIGMA) blocks.
> http://www.cognigencorp.com/nonmem/current/2009-July/1804.html
>
> IMHO the only sensible way for a human to write a covariance matrix is us=
ing
> the SD and correlation style introduced in NONMEM 7. If I understand thin=
gs
> correctly (this kind of math is not my strong point) it should be impossi=
ble
> to write non-positive definitive matrices provided the correlations are
> -1<=0<=+1 (although NONMEM might complain at -1 or +1).
>
> Can anyone provide any good reason to prefer the variance-covariance
> parameterization over the SD-correlation parametization?
>
> The other issue brought up by Pavel is a different one -- how to specify =
more
> simply a zero correlation|covariance between random effects.
>
> Monolix has a 0|1 matrix that allows a user to fix any covariance matrix
> element to zero so presumably the math exists somewhere to allow NONMEM h=
ow to
> understand a matrix to be recognized in this format.
>
> Bob (Bauer) -- do you know how to do this? Perhaps you could talk to Ma=
rc
> Lavielle and consider adding something similar to NM-TRAN. I'd prefer the
> style shown by Pavel to having to write a separate matrix of
> 0|1 elements as used in the Monolix GUI.
>
> Best wishes,
>
> Nick
>
>
>
> On 13/06/2012 5:16 p.m., Paolo Denti wrote:
> > Dear Pavel and Ken,
> > I also share the same pain, especially when coding correlations for
> > between-occasion variability ETAs. This implies reorganizing them in
> > blocks and renumbering everything. A real chore, and a very
> > error-prone process.
> >
> > Maybe one can think of using the OMEGA in the correlation format,
> > which should make it easier to write "legitimate" OMEGA matrices. Or
> > NONMEM can check the positive definiteness of the initial estimate and
> > complain if necessary (I believe it already does so).
> >
> > So, dear NONMEM developers, please count a +1 in the survey for this
> > feature. :)
> >
> > Thank you and ciao,
> > Paolo
> >
> > On 2012/06/12 17:08, Ken Kowalski wrote:
> >>
> >> Dear Pavel,
> >>
> >> I certainly feel your pain but you have to be careful how you fix
> >> certain elements in Omega to ensure that you have a valid positive
> >> definite covariance matrix. The starting values in your \$OMEGA block
> >> do not give rise to a valid covariance matrix. Note in particular
> >> that the covariance between ETA3 and ETA4 is too large relative to
> >> the variances for ETA3 and ETA4 such that the correlation is greater
> >> than 1.0, i.e.,
> >>
> >> Corr(ETA3,ETA4) = 0.03/[SQRT(0.0166)*SQRT(0.0166)]=1.807 > 1
> >>
> >> You also have the same problem for Corr(ETA1,ETA3) > 1.
> >>
> >> Ken
> >>
> >> Kenneth G. Kowalski
> >>
> >> President & CEO
> >>
> >> A2PG - Ann Arbor Pharmacometrics Group, Inc.
> >>
> >> 110 Miller Ave., Garden Suite
> >>
> >> Ann Arbor, MI 48104
> >>
> >> Work: 734-274-8255
> >>
> >> Cell: 248-207-5082
> >>
> >> Fax: 734-913-0230
> >>
> >> ken.kowalski
> >>
> >> www.a2pg.com <http://www.a2pg.com>
> >>
> >> *From:*owner-nmusers
> >> [mailto:owner-nmusers
> >> *nonmem
> >> *Sent:* Tuesday, June 12, 2012 10:24 AM
> >> *To:* nmusers
> >> *Subject:* [NMusers] omega matrix - friendly suggestion
> >>
> >> Hello NONMEM Community,
> >>
> >> Sometimes it takes more time to choose the best omega matrix than to
> >> develop a PD model. Selecting the omega matric is a tedious, time
> >> consuming and less creative part of the model development. I hope
> >> you feel my pain. Will it be helpful to rewrite the NONMEM software
> >> so that any element of the omega matrix can be fixed to any value?
> >> It may look like this:
> >>
> >> \$OMEGA BLOCK (4) 3.60E-02 FIX
> >>
> >> 0.01 3.23E-02
> >> 0.03 0 FIX 1.66E-02
> >> 0.01 0 FIX 0.03 FIX 1.66E-02
> >>
> >> This change can make many NONMEM users happy.
> >>
> >> Thanks!
> >>
> >> Pavel
> >>
> >
> > --
> > ------------------------------------------------
> > Paolo Denti, PhD
> > Junior Lecturer
> > Division of Clinical Pharmacology
> > Department of Medicine
> > University of Cape Town
> >
> > K45 Old Main Building
> > Groote Schuur Hospital
> > Observatory, Cape Town
> > 7925 South Africa
> > phone: +27 21 404 7719
> > fax: +27 21 448 1989
> > email:paolo.denti
> > ------------------------------------------------
>
> --
> Nick Holford, Professor Clinical Pharmacology
>
> First World Conference on Pharmacometrics, 5-7 September 2012 Seoul, Kore=
a
> http://www.go-wcop.org
>
> Dept Pharmacology& Clinical Pharmacology, Bldg 505 Room 202D University =
of
> Auckland,85 Park Rd,Private Bag 92019,Auckland,New Zealand
> tel:+64(9)923-6730 fax:+64(9)373-7090 mobile:+64(21)46 23 53
> email: n.holford
> http://www.fmhs.auckland.ac.nz/sms/pharmacology/holford
>
>
>
>
> ICON plc made the following annotations.
> -------------------------------------------------------------------------=
-----
> This e-mail transmission may contain confidential or legally privileged
> information that is intended only for the individual or entity named in t=
he e-
> mail address. If you are not the intended recipient, you are hereby notif=
ied
> that any disclosure, copying, distribution, or reliance upon the contents=
of
> this e-mail is strictly prohibited. If you have received this e-mail
> transmission in error, please reply to the sender, so that ICON plc can
> arrange for proper delivery, and then please delete the message.
> Thank You,
> ICON plc
> Leopardstown
> Dublin 18
> Ireland
> Registered number:
> 145835_________________________________________________________________
> NOTICE: The information contained in this electronic mail message is inte=
nded
> only for the personal and confidential use of the designated recipient(s)
> named above. This message may be an attorney-client communication, may be
> protected by the work product doctrine, and may be subject to a protectiv=
e
> order. As such, this message is privileged and confidential. If the reade=
r of
> this message is not the intended recipient or an agent responsible for
> delivering it to the intended recipient, you are hereby notified that you=
have
> received this message in error and that any review, dissemination,
> distribution, or copying of this message is strictly prohibited. If you h=
ave