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

[agalecki at umich.edu]

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