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 <[EMAIL PROTECTED]> wrote: > > >>Sounds like the documentation should be changed. Doug? Jose? >> >>- Mary >> >>Andrzej Galecki ([EMAIL PROTECTED]) 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 ([EMAIL PROTECTED]) 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: [EMAIL PROTECTED] >> > > > > > To: Andrzej Galecki <[EMAIL PROTECTED]> >> > > > > > >> > > > > > 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 >> > > > > >> > > > > >> > > > > >> > > > > >> > > > > >> > > >> > > >> > > >> > > >> >> >> > > > > ______________________________________________ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel