Bug in getVarCov.gls method (PR#9752)

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