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 have 16 subjects with 1-4 obs per subject.
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)
[[alternative HTML version deleted]]
On 11/16/05, Wassell, James T., Ph.D. <jtw2 at cdc.gov> wrote:> 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).So to paraphrase, your model can be written as (with the index i representing subject) y_ij = \mu + b_i + e_ij where b_i ~ N(0, \tao^2) e_ij ~ N(0, \sigma_2) and all b_i's and e_ij's are mutually independent. The model has, as you say, 3 parameters, \mu, \tao and \sigma.> I have 16 subjects with 1-4 obs per subject. > > I need a 3x3 variance-covariance matrix that includes all 3 parameters > in order to compute the variance of a specific linear combination.Can you specify the 'linear combination' that you want to estimate in terms of the model above? Deepayan
-----Original Message----- From: Wassell, James T., Ph.D. Sent: Thursday, November 17, 2005 9:40 AM To: 'Deepayan Sarkar' Subject: RE: nlme question Deepayan, Thanks for your interest. It's difficult in email but I need the variance of Kappa = mu + 1.645*tau + 1.645* sigma Just using the standard variance calculation Var(K) = Var(mu) + (1.645)^2 * Var(tau) + 1.645^2 * var(sigma) + [three covariance terms with constant multipliers]. So, I can get var(mu) [or actually the standard error] from the summary function -- and var(tau) and var(sigma) using the VarCorr function, but I still need the covariance terms. I am trying to duplicate methods in a paper by Nicas & Neuhaus, "Variability in Respiratory Protection and the Assigned Protection Factor" J. Occ & Environ Health, Feb. 2004. p 99-109. (see eqn. 12). The authors used Proc Mixed, but I can't figure out how to get covariance terms with SAS either. Thanks again. Terry Wassell -----Original Message----- From: Deepayan Sarkar [mailto:deepayan.sarkar at gmail.com] Sent: Thursday, November 17, 2005 2:52 AM To: Wassell, James T., Ph.D. Cc: r-help at stat.math.ethz.ch Subject: Re: nlme question On 11/16/05, Wassell, James T., Ph.D. <jtw2 at cdc.gov> wrote:> 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).So to paraphrase, your model can be written as (with the index i representing subject) y_ij = \mu + b_i + e_ij where b_i ~ N(0, \tao^2) e_ij ~ N(0, \sigma_2) and all b_i's and e_ij's are mutually independent. The model has, as you say, 3 parameters, \mu, \tao and \sigma.> I have 16 subjects with 1-4 obs per subject. > > I need a 3x3 variance-covariance matrix that includes all 3 parameters > in order to compute the variance of a specific linear combination.Can you specify the 'linear combination' that you want to estimate in terms of the model above? Deepayan
James, By assumption, sigma and tau are assumed uncorrelated, which I believe Deepayan noted below. Sigma is also random error so it is uncorrelated with your fixed effects. What are the covariance terms you are seeking? -----Original Message----- From: r-help-bounces at stat.math.ethz.ch [mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Wassell, James T., Ph.D. Sent: Thursday, November 17, 2005 10:29 AM To: r-help at stat.math.ethz.ch Subject: [R] nlme question -----Original Message----- From: Wassell, James T., Ph.D. Sent: Thursday, November 17, 2005 9:40 AM To: 'Deepayan Sarkar' Subject: RE: nlme question Deepayan, Thanks for your interest. It's difficult in email but I need the variance of Kappa = mu + 1.645*tau + 1.645* sigma Just using the standard variance calculation Var(K) = Var(mu) + (1.645)^2 * Var(tau) + 1.645^2 * var(sigma) + [three covariance terms with constant multipliers]. So, I can get var(mu) [or actually the standard error] from the summary function -- and var(tau) and var(sigma) using the VarCorr function, but I still need the covariance terms. I am trying to duplicate methods in a paper by Nicas & Neuhaus, "Variability in Respiratory Protection and the Assigned Protection Factor" J. Occ & Environ Health, Feb. 2004. p 99-109. (see eqn. 12). The authors used Proc Mixed, but I can't figure out how to get covariance terms with SAS either. Thanks again. Terry Wassell -----Original Message----- From: Deepayan Sarkar [mailto:deepayan.sarkar at gmail.com] Sent: Thursday, November 17, 2005 2:52 AM To: Wassell, James T., Ph.D. Cc: r-help at stat.math.ethz.ch Subject: Re: nlme question On 11/16/05, Wassell, James T., Ph.D. <jtw2 at cdc.gov> wrote:> 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).So to paraphrase, your model can be written as (with the index i representing subject) y_ij = \mu + b_i + e_ij where b_i ~ N(0, \tao^2) e_ij ~ N(0, \sigma_2) and all b_i's and e_ij's are mutually independent. The model has, as you say, 3 parameters, \mu, \tao and \sigma.> I have 16 subjects with 1-4 obs per subject. > > I need a 3x3 variance-covariance matrix that includes all 3 parameters> in order to compute the variance of a specific linear combination.Can you specify the 'linear combination' that you want to estimate in terms of the model above? Deepayan ______________________________________________ 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
Thank you for taking the time to think about my problem. The reference states: "The covariance structure must be considered, because for unbalanced data the estimates" (i.e. mu, sigma and tau hats) "are not typically independent." Page 105. It would be nice to simply assume zero covariance terms, but the authors reject this simplification.
I think the authors are mistaken. Sigma is random error, and due to its randomness it cannot be systematically related to anything. It is this ind. assumption that allows for the likelihood to be expressed as described in Pinhiero and Bates p.62. If you are finding that sigma is systematically related to a fixed effect, then your model is misspecified and there are omitted characteristics that need to be accounted for. -----Original Message----- From: Wassell, James T., Ph.D. [mailto:jtw2 at cdc.gov] Sent: Thursday, November 17, 2005 11:52 AM To: Doran, Harold; r-help at stat.math.ethz.ch Subject: RE: [R] nlme question Thank you for taking the time to think about my problem. The reference states: "The covariance structure must be considered, because for unbalanced data the estimates" (i.e. mu, sigma and tau hats) "are not typically independent." Page 105. It would be nice to simply assume zero covariance terms, but the authors reject this simplification.
Deepayan, Yes, thanks for confirming my suspicions. I know mixed models are "different" but, I did not think they were so different as to preclude estimating the var-cov matrix (via the Hessian in Maximum likelihood, as you point out). Thanks for prompting me to think about MCMC. Your suggestion to consider MCMC makes me realize that using BUGS, I could directly sample from the posterior of the linear combination of parameters - to get its variance and eliminate the extra step using the var-cov matrix. As you say, with results better than the asymptotic approximation. (Maybe I can do the same thing with mcmcsamp?, but I'm not familiar with this and will have to take a look at it.) -----Original Message----- From: Deepayan Sarkar [mailto:deepayan.sarkar at gmail.com] Sent: Thursday, November 17, 2005 2:22 PM To: Doran, Harold Cc: Wassell, James T., Ph.D.; r-help at stat.math.ethz.ch Subject: Re: nlme question On 11/17/05, Doran, Harold <HDoran at air.org> wrote:> I think the authors are mistaken. Sigma is random error, and due toits> randomness it cannot be systematically related to anything. It is this > ind. assumption that allows for the likelihood to be expressed as > described in Pinhiero and Bates p.62.I think not. The issue is dependence between the _estimates_ of sigma, tao, etc, and that may well be present. Presumably, if one can compute the likelihood surface as a function of the 3 parameters, the hessian at the MLE's would give the estimated covariance. However, I don't think nlme does this. A different approach you might want to consider is using mcmcsamp in the lme4 package (or more precisely, the Matrix package) to get samples from the joint posterior distribution. This is likely to be better than the asymptotic normal approximation in any case. Deepayan