RE: omega matrix - friendly suggestion

________________________________________ From: Stephen Duffull [[email protected]] Sent: Thursday, June 14, 2012 1:51 PM To: Bob Leary; Bauer, Robert; Matt Hutmacher; 'Nick Holford'; [email protected] Subject: RE: [NMusers] 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 > -----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 > > > > >