Hello, I am using R2.5 under Windows. Looks like the following statement vars <- (obj$sigma^2)*vw in getVarCov.gls method (nlme package) needs to be replaced with: vars <- (obj$sigma*vw)^2 With best regards Andrzej Galecki Douglas Bates wrote:>I'm not sure when the getVarCov.gls method was written or by whom. To >tell the truth I'm not really sure what it should do. > >Does anyone have recommendations about what to do? The simplest >thing, if the method is giving incorrect results, is to eliminate the >method entirely. Any objections to my doing that? > >The method is currently defined as > >getVarCov.gls <- > function(obj, individual = 1, ...) >{ > S <- corMatrix(obj$modelStruct$corStruct)[[individual]] > if (!is.null( obj$modelStruct$varStruct)) > { > ind <- obj$groups==individual > vw <- 1/varWeights(obj$modelStruct$varStruct)[ind] > } > else vw <- rep(1,nrow(S)) > vars <- (obj$sigma^2)*vw > result <- t(S * sqrt(vars))*sqrt(vars) > class(result) <- c("marginal","VarCov") > attr(result,"group.levels") <- names(obj$groups) > result >} > > >On 2/23/06, Mary Lindstrom <lindstro at biostat.wisc.edu> wrote: > > >>Sounds like the documentation should be changed. Doug? Jose? >> >>- Mary >> >>Andrzej Galecki (agalecki at umich.edu) writes: >> > Hi Mary, >> > >> > Thank you for your response. Note that we followed R documentation, >> > which states that getVarCov() can be used with gls() objects. >> > >> > Thanks again. >> > >> > Andrzej >> > >> > Mary Lindstrom wrote: >> > >> > >Sorry once again for the slow response. When I wrote getVarCov I was >> > >not thinking of using it for gls objects. Sounds like it doesn't work >> > >for them. Either the code should be changed to reject the request when >> > >the object has class gls or it should be rewritten to work correctly. >> > > >> > >- Mary >> > > >> > >Andrzej Galecki (agalecki at umich.edu) writes: >> > > > Hi Mary, >> > > > >> > > > Our question is about variance-covariance matrix, hat(sigma_i), >> > > > returned by getVarCov() function when applied to model fit generated by >> > > > gls( ) function. It appears that at least in our example the returned >> > > > matrix was incorrect. Are you aware of any problem with getVarCov () >> > > > when applied to an object generated by gls()? >> > > > >> > > > To be more specific our impression is that getVarCov() expects to get >> > > > variances as an input >> > > > and instead it gets standard deviations and corresponding ratios. >> > > > >> > > > Sorry we did not state our question more precisely. >> > > > >> > > > Thank you >> > > > >> > > > Andrzej >> > > > >> > > > PS. Here is our original question to Jose. His response was that the >> > > > gls() syntax was correct and he directed us to you or to D.Bates >> > > > >> > > > We used gls() to fit a marginal model with unstructured >> > > > covariance matrix. Could you verify gls() syntax, especially part >> > > > specifying correlation matrix and variance function? Assuming that our >> > > > gls() code is correct, we think that getVarCoV() returns incorrect >> > > > marginal covariance matrix in this example. Are you aware of any >> > > > problems with getVarCov(), when used with a model fit returned by >> > > > gls()? If you think that getVarCov() should work fine in the context of >> > > > this model we are ready to send you a relevant code, in which we extract >> > > > information from gls() fit and 'manualy' calculate marginal >> > > > variance-covariance matrix and cross-check it with the output from SAS. >> > > > Our impression is that getVarCov() expects to get variances as an input >> > > > and instead it gets standard deviations and corresponding ratios. We >> > > > would really like to get to the bottom of this issue. >> > > > >> > > > >> > > > Mary Lindstrom wrote: >> > > > >> > > > >Hi Andrzej >> > > > > >> > > > > > -------- Original Message -------- >> > > > > > Subject: Re: [Fwd: Mixed Models book.] >> > > > > > Date: Thu, 19 Jan 2006 13:06:09 -0500 >> > > > > > From: jose.pinheiro at novartis.com >> > > > > > To: Andrzej Galecki <agalecki at umich.edu> >> > > > > > >> > > > > > The getVarCov function was actually written by Mary Lindstrom and, as a >> > > > > > far as I know, adapted to R >> > > > > > by Douglas Bates. I did think about including it in the S-PLUS version >> > > > > > as well, but never got around >> > > > > > to it. My understanding is that it should also work with gls objects, >> > > > > > but I believe the main motivation >> > > > > > was to have it for lme objects. I am not sure if the example is >> > > > > > triggering some possible bug in the >> > > > > > code, but perhaps you can ask Mary or Douglas if an updated version may >> > > > > > be available. I will not have >> > > > > > time to look at the code and get back to you before the deadlines you >> > > > > > have for submitting the book. >> > > > > > Have you tried using getVarCov with simpler gls models (e.g., compound >> > > > > > symmetry with no variance functions) >> > > > > > to see if you get the same type of problem? >> > > > > >> > > > >The point of getVarCov was to let students compute and see hat(D) the >> > > > >estimate of the Variance matrix of the random effects and >> > > > >hat(Sigma_i), the estimate of the marginal variance of y_i. If you are >> > > > >fitting a General Linear Model (not a generalized mixed effects linear >> > > > >model just to be clear) then there are no random effects so computing >> > > > >hat(D) makes no sense. You could however still compute out >> > > > >hat(Sigma_i). >> > > > > >> > > > >Does this answer your question? - Mary >> > > > > >> > > > > >> > > > > >> > > > > >> > > > > >> > > >> > > >> > > >> > > >> >> >> > > > >
On Mon, 25 Jun 2007, agalecki at umich.edu wrote:> I am using R2.5 under Windows.I presume you mean 2.5.0 (there is no R2.5: see the posting guide). But which version of nlme, which is the relevant fact here? The R posting guide suggests showing the output of sessionInfo() to establish the environment used.> Looks like the following statement > > vars <- (obj$sigma^2)*vw > > in getVarCov.gls method (nlme package) needs to be replaced with: > > vars <- (obj$sigma*vw)^2We need some evidence! Please supply a reproducible example and give your reasoning. For example, the FAQ says The most important principle in reporting a bug is to report _facts_, not hypotheses or categorizations. It is always easier to report the facts, but people seem to prefer to strain to posit explanations and report them instead. If the explanations are based on guesses about how R is implemented, they will be useless; others will have to try to figure out what the facts must have been to lead to such speculations. Sometimes this is impossible. But in any case, it is unnecessary work for the ones trying to fix the problem. It is very useful to try and find simple examples that produce apparently the same bug, and somewhat useful to find simple examples that might be expected to produce the bug but actually do not. If you want to debug the problem and find exactly what caused it, that is wonderful. You should still report the facts as well as any explanations or solutions. Please include an example that reproduces the problem, preferably the simplest one you have found. It should be easily possible to cross-check an example with one of the many other ways available to do GLS fits in R. [...] -- Brian D. Ripley, ripley at stats.ox.ac.uk Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UK Fax: +44 1865 272595
Two attachments: 1. getVarCovBugReport.R - Rcode with an example illustrating the problem and how to fix it 2. getVarCovBugReportSession.txt contains code and R session results. Thank you Andrzej Galecki Prof Brian Ripley wrote:> On Mon, 25 Jun 2007, agalecki at umich.edu wrote: > >> I am using R2.5 under Windows. > > > I presume you mean 2.5.0 (there is no R2.5: see the posting guide). > But which version of nlme, which is the relevant fact here? The R > posting guide suggests showing the output of sessionInfo() to > establish the environment used. > >> Looks like the following statement >> >> vars <- (obj$sigma^2)*vw >> >> in getVarCov.gls method (nlme package) needs to be replaced with: >> >> vars <- (obj$sigma*vw)^2 > > > We need some evidence! Please supply a reproducible example and give > your reasoning. For example, the FAQ says > > The most important principle in reporting a bug is to report _facts_, > not hypotheses or categorizations. It is always easier to report the > facts, but people seem to prefer to strain to posit explanations and > report them instead. If the explanations are based on guesses about > how > R is implemented, they will be useless; others will have to try to > figure > out what the facts must have been to lead to such speculations. > Sometimes this is impossible. But in any case, it is unnecessary work > for the ones trying to fix the problem. > > It is very useful to try and find simple examples that produce > apparently the same bug, and somewhat useful to find simple examples > that might be expected to produce the bug but actually do not. If you > want to debug the problem and find exactly what caused it, that is > wonderful. You should still report the facts as well as any > explanations or solutions. Please include an example that reproduces > the > problem, preferably the simplest one you have found. > > It should be easily possible to cross-check an example with one of the > many other ways available to do GLS fits in R. > > [...] > >
Thank you, will be fixed in the next release of nlme (not yet scheduled, as everythign is code frozen for the release of R 2.5.1 tomorrow). On Tue, 26 Jun 2007, Andrzej Galecki wrote:> > Two attachments: > > 1. getVarCovBugReport.R - Rcode with an example illustrating the problem and > how to fix it > 2. getVarCovBugReportSession.txt contains code and R session results. > > Thank you > > Andrzej Galecki > > > > Prof Brian Ripley wrote: > >> On Mon, 25 Jun 2007, agalecki at umich.edu wrote: >> >>> I am using R2.5 under Windows. >> >> >> I presume you mean 2.5.0 (there is no R2.5: see the posting guide). But >> which version of nlme, which is the relevant fact here? The R posting >> guide suggests showing the output of sessionInfo() to establish the >> environment used. >> >>> Looks like the following statement >>> >>> vars <- (obj$sigma^2)*vw >>> >>> in getVarCov.gls method (nlme package) needs to be replaced with: >>> >>> vars <- (obj$sigma*vw)^2 >> >> >> We need some evidence! Please supply a reproducible example and give your >> reasoning. For example, the FAQ says >> >> The most important principle in reporting a bug is to report _facts_, >> not hypotheses or categorizations. It is always easier to report the >> facts, but people seem to prefer to strain to posit explanations and >> report them instead. If the explanations are based on guesses about how >> R is implemented, they will be useless; others will have to try to figure >> out what the facts must have been to lead to such speculations. >> Sometimes this is impossible. But in any case, it is unnecessary work >> for the ones trying to fix the problem. >> >> It is very useful to try and find simple examples that produce >> apparently the same bug, and somewhat useful to find simple examples >> that might be expected to produce the bug but actually do not. If you >> want to debug the problem and find exactly what caused it, that is >> wonderful. You should still report the facts as well as any >> explanations or solutions. Please include an example that reproduces the >> problem, preferably the simplest one you have found. >> >> It should be easily possible to cross-check an example with one of the many >> other ways available to do GLS fits in R. >> >> [...] >> >> > >-- Brian D. Ripley, ripley at stats.ox.ac.uk Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UK Fax: +44 1865 272595