R help - Jun 2006 - setting the random-effects covariance matrix in lme

Dear R-users,
 
I have longitudinal data and would like to fit a model where both the
variance-covariance matrix of the random effects and the residual variance are
conditional on a (binary) grouping variable.
I guess the model would have the following form (in hierarchical notation)
Yi|bi,k ~ N(XiB+Zibi, sigmak*Ident)
bi|k ~ N(0, Dk)
K~Bernoulli(p)
 
I can obtain different sigmas (sigma0 and sigma1 based on the factor
'dx') using the weights option in the call to lme:
lme(fixed = height ~ -1 + bage + mat_ht + pat_ht + dx + time:dx + time2:dx , 

                 weights=varIdent(form=~1|dx),

                  random = ~ time + time2 |subject, data = pilot,
na.action=na.omit)

 but I cannot seem to be able to get two different matrices for the random
effects.

I'm particularly interested in obtaining,
 
Y|k ~ N(XB, Z(Dk)Z+sigmak),
 
so I need to extract D0 and D1 form the lme object.  I've looked in Pinheiro
and Bates but was unable to identify how to fit this type of covariance matrix
from the examples provided in the book - perhaps it is there expressed in a
different form?  I looked in section 4.2.2 (pdMat) of teh book but I'm still
not sure how to pass it on in the call.
 
Any help would be greatly appreciated,
 
Regina
 
Regina Rendas-Baum
Graduate Student
Brown University
Biostatistics

I haven't seen a reply to this post, so I will offer a few comments. 
  I don't know off the top of my head how to do what you want, but I 
believe it can be done.

	  Are you aware that the standard R distribution includes script files 
containing nearly all the R commands required to produce the analyses 
described in Pinheiro and Bates?  Find where R is installed on your 
computer, and then look for "~R/library/nlme/scripts".  This should
make
it vastly easier to work through the book line by line.  [It does for 
me.  For example, I keep forgetting that a "b^2" term in the book must
usually be written "I(b^2)".  Also, I've spent hours trying to 
understand why my results differed so much from the book only to find a 
stupid typographical error.  These problems disappear when using the 
script file.]

	  If this fails, please provide a simple, self-contained example, e.g., 
modifying an example in the book to show something you tried, and 
explain why it didn't work for you.

	  Hope this helps.
	  Spencer Graves
p.s.  You mentioned you looked at Section 4.2.2.  Have you considered 
Section 5.2  "Variance Functions for Modeling Heteroscedasticity". 
What
you describe sounds like "heterscedasticity" to me.

Rendas-Baum, Regina wrote:
> Dear R-users,
>  
> I have longitudinal data and would like to fit a model where both 
the variance-covariance matrix of the random effects and the residual
variance are conditional on a (binary) grouping
variable.
> I guess the model would have the following form (in hierarchical 
notation)
> Yi|bi,k ~ N(XiB+Zibi, sigmak*Ident)
> bi|k ~ N(0, Dk)
> K~Bernoulli(p)
>  
> I can obtain different sigmas (sigma0 and sigma1 based on the factor 
'dx') using the weights option in the call to
lme:
> lme(fixed = height ~ -1 + bage + mat_ht + pat_ht + dx + time:dx + time2:dx
,
> 
>    weights=varIdent(form=~1|dx),
> 
>    random = ~ time + time2 |subject, data = pilot, na.action=na.omit)
> 
>  but I cannot seem to be able to get two different matrices for the 
random effects.
> 
> I'm particularly interested in obtaining,
>  
> Y|k ~ N(XB, Z(Dk)Z+sigmak),
>  
> so I need to extract D0 and D1 form the lme object.  I've looked in 
Pinheiro and Bates but was unable to identify how to fit this type of
covariance matrix from the examples provided in the book - perhaps it
is there expressed in a different form?  I looked in section 4.2.2
(pdMat) of teh book but I'm still not sure how to pass it on in the
call.
>  
> Any help would be greatly appreciated,
>  
> Regina
>  
> Regina Rendas-Baum
> Graduate Student
> Brown University
> Biostatistics
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html