[R] nmle question

[Wassell, James T., Ph.D.]

Hello.

I have 16 subjects with 1-4 obs per subject. 

I am using the package "nlme" to fit a simple random effects (variance
components model) with 3 parameters:  overall mean (fixed effect),
between subject variance (random) and within subject variance (random).
 
I need a 3x3 variance-covariance matrix that includes all 3 parameters
in order to compute the variance of a specific linear combination.

But I can't get the 3x3 matrix.  Should I specify the formulae in lme
differently or is there some other suggestion that I might try?

Thank you very much for any advice.  my data and code follows.  

"mydata" <-
structure(list(subject = c(17, 17, 17, 17, 5, 16, 16, 8, 8, 8, 
8, 7, 7, 7, 7, 9, 9, 9, 10, 10, 11, 11, 11, 12, 12, 12, 12, 14, 
14, 14, 14, 15, 15, 15, 15, 13, 13, 13, 1, 1, 1, 2, 2, 2, 2, 
3, 3, 3, 4, 4, 4), y = c(-2.944, -5.521, -4.644, -4.736, -5.799, 
-4.635, -5.986, -5.011, -3.989, -4.682, -6.975, -6.064, -5.991, 
-8.068, -5.075, -5.298, -6.446, -5.037, -6.534, -4.828, -5.886, 
-4.025, -6.607, -5.914, -4.159, -6.757, -4.564, -5.011, -5.416, 
-5.371, -5.768, -7.962, -5.635, -4.575, -5.268, -6.975, -5.598, 
-7.669, -8.292, -7.265, -5.858, -7.003, -3.807, -5.829, -5.613, 
-3.135, -5.136, -5.394, -5.011, -5.598, -4.174)), .Names = c("subject", 
"y"), class = "data.frame", row.names = c("1", "2", "3", "4", 
"5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", 
"16", "17", "18", "19", "20", "21", "22", "23", "24", "25", "26", 
"27", "28", "29", "30", "31", "32", "33", "34", "35", "36", "37", 
"38", "39", "40", "41", "42", "43", "44", "45", "46", "47", "48", 
"49", "50", "51"))

lme.res<-lme(fixed=y~1,data=mydata,random=~1|subject,method="ML")
VarCorr(lme.res)