[R] setting the random-effects covariance matrix in lme

[Rendas-Baum, Regina]

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