[Rd] Bug in getVarCov.gls method (PR#9752)

From: <agalecki_at_umich.edu>
Date: Mon, 25 Jun 2007 22:17:12 +0200 (CEST)


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
>> > > > >
>> > > > >
>> > > > >
>> > > > >
>> > > > >
>> > >
>> > >
>> > >
>> > >
>>
>>
>>
>
>
>
>



R-devel_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel Received on Tue 26 Jun 2007 - 07:28:20 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Wed 27 Jun 2007 - 07:35:57 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-devel. Please read the posting guide before posting to the list.