RE: omega matrix - friendly suggestion
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 determinants.
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 code search
algorithms that require positive definite matrices as a search step :-)
Sampling algorithms don't have this issue as there are good priors that
guarantee "legal" matrices.
Steve
Quoted reply history
> -----Original Message-----
> From: [email protected] [mailto:[email protected]] On
> Behalf Of Bob Leary
> Sent: Friday, 15 June 2012 5:11 a.m.
> To: Bauer, Robert; Matt Hutmacher; 'Nick Holford'; [email protected]
> Subject: RE: [NMusers] omega matrix - friendly suggestion
>
> Easy countereaxample to Nick's hypothesis that a covariance matrix with an
> underlying correlation matrix with all elements of magnitude <=1 is positive
> 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 used to
> force an arbitrary element of a full block matrix to a given value while
> preserving postive definiteness.
> _______________________________________
> From: [email protected] [[email protected]] On Behalf Of
> Bauer, Robert [[email protected]]
> Sent: Wednesday, June 13, 2012 4:05 PM
> To: Matt Hutmacher; 'Nick Holford'; [email protected]
> Subject: RE: [NMusers] omega matrix - friendly suggestion
>
> Thank you all for your suggestions. The more flexible OMEGA element entries
> 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 Cholesky
> format in NM 7.2, $OMEGA CHOLESKY may be used.
> 2) The traditional methods (FO/FOCE/Laplace) actually vary the omega elements
> in the Cholesky format for estimation, than transform back in the variance
> format after each iteration. This assures that as NONMEM arbitrarily searches
> 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-diagonal
> variance element to a non-zero value within a block of variances cannot be
> honored by the search algorithm past the first iteration. Additionally, to my
> best understanding, arbitrarily fixing an off-diagonal variance element to 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 can 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, but
> are constructed as a consequence of the expectation step, which provide
> conditional means and variances from which the omega elements are evaluated.
> 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: [email protected]
> Web: www.iconplc.com
>
> -----Original Message-----
> From: [email protected] [mailto:[email protected]] On
> Behalf Of Matt Hutmacher
> Sent: Wednesday, June 13, 2012 2:35 PM
> To: 'Nick Holford'; [email protected]
> Subject: RE: [NMusers] omega matrix - friendly suggestion
>
> Dear Paolo, Nick, all,
>
> A way to ensure a nonnegative definite (NND) variance-covariance matrix is 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 convergence 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 with
> 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 WRES 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 variances
> into linear combinations of thetas affecting the residuals). It is not enough
> to just model variances and correlations; such a parameterization will not
> 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. Upon
> 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 provide
> a lot of flexibility for analysis (as opposed to hand coding discussed below).
>
> There is also a modified Cholesky that can be used, which has the nice feature
> that one can simply, conveniently remove all the variance components related
> to a certain random effect (Ref 2). This was recently used in the attached
> reference to aid in selection of the $OMEGA structure with the fixed effects.
>
> 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-1077.
>
>
> -----Original Message-----
> From: [email protected] [mailto:[email protected]] On
> Behalf Of Nick Holford
> Sent: Wednesday, June 13, 2012 1:20 PM
> To: [email protected]
> Subject: Re: [NMusers] omega matrix - friendly suggestion
>
> Paolo,
>
> You may remember this discussion about the history of the correlation feature
> 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 using
> the SD and correlation style introduced in NONMEM 7. If I understand things
> correctly (this kind of math is not my strong point) it should be impossible
> 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 how to
> understand a matrix to be recognized in this format.
>
> Bob (Bauer) -- do you know how to do this? Perhaps you could talk to Marc
> 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
> >>
> >> [email protected] <mailto:[email protected]>
> >>
> >> www.a2pg.com http://www.a2pg.com
> >>
> >> *From:*[email protected]
> >> [mailto:[email protected]] *On Behalf Of
> >> *[email protected]
> >> *Sent:* Tuesday, June 12, 2012 10:24 AM
> >> *To:* [email protected]
> >> *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:[email protected]
> > ------------------------------------------------
>
> --
> Nick Holford, Professor Clinical Pharmacology
>
> First World Conference on Pharmacometrics, 5-7 September 2012 Seoul, Korea
> 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: [email protected]
> http://www.fmhs.auckland.ac.nz/sms/pharmacology/holford
>
>
>
>
>