If I understand your request correctly, you want to use something like
"weights=varIdent(...)" as an argument to lme(). varIdent and the other
varFunc constructors have an argument "fixed" that allow you to specify
values for some or all of the coefficients of the variance function.
See ?varIdent. The actual error variance will be varFunc() * sigma^2,
where sigma^2 is estimated.
R. Woodrow Setzer, Jr. Phone:
(919) 541-0128
Experimental Toxicology Division Fax: (919)
541-5394
Pharmacokinetics Branch
NHEERL MD-74; US EPA; RTP, NC 27711
|---------+------------------------------>
| | "J.R. Lockwood" |
| | <lockwood at rand.org>|
| | Sent by: |
| | owner-r-help at stat.m|
| | ath.ethz.ch |
| | |
| | |
| | 08/29/02 02:00 PM |
| | |
|---------+------------------------------>
>------------------------------------------------------------------------------------------------------|
| |
| To: r-help at stat.math.ethz.ch |
| cc: "J.R. Lockwood" <lockwood at rand.org> |
| Subject: [R] lme() with known level-one variances |
>------------------------------------------------------------------------------------------------------|
Greetings,
I have a meta-analysis problem in which I have fixed effects
regression coefficients (and estimated standard errors) from identical
models fit to different data sets. I would like to use these results
to create pooled estimated regression coefficients and estimated
standard errors for these pooled coefficients. In particular, I would
like to estimate the model
\beta_{i} = \mu + \eta_{i} + \epsilon_{i}
\eta_{i} ~ iid N(0,\tau^2) and independent of the \epsilon_{i}, the
latter themselves being independent with variances assumed known and
equal to the squared standard errors reported in the regression
output.
I would like to use lme() to estimate \tau^2 by REML, and also get a
sensibly weighted estimate for \mu from the fixed effects output. I
am not sure how to do this. I have tried
lme(fixed=beta~1,random=~1|group,weights=~beta.v)
where "beta" are my coefficients, "group" is a trivial factor
indicating that each observation is its own group, and "beta.v" are
the squared standard errors. Whatever I get out of this doesn't make
sense to me, and I suspect that I have specified the model
incorrectly.
Incidentally, if I just run the simple unidentifiable model
lme(fixed=beta~1,random=~1|group)
lme() somehow manages to produce estimates of the two variance
components, although the estimated confidence intervals are huge and
contain zero. If I square and sum the estimated variance components,
I do get the sample variance of my regression coefficients, but why
that particular parceling of variance was chosen as opposed to any
other with the same property eludes me.
Here are my specs:
platform i686-pc-linux-gnu
arch i686
os linux-gnu
system i686, linux-gnu
status
major 1
minor 5.1
year 2002
month 06
day 17
language R
Thanks in advance for your help -- I've learned a ton of statistics
and computing on this list.
J.R. Lockwood
412-683-2300 x4941
lockwood at rand.org
http://www.rand.org/methodology/stat/members/lockwood/
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.
-.-.-.-
r-help mailing list -- Read
http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._.
_._._._
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
lme() with known level-one variances
3 messages · Setzer.Woodrow@epamail.epa.gov, Thomas Lumley, J.R. Lockwood
On Fri, 30 Aug 2002 Setzer.Woodrow at epamail.epa.gov wrote:
If I understand your request correctly, you want to use something like "weights=varIdent(...)" as an argument to lme(). varIdent and the other varFunc constructors have an argument "fixed" that allow you to specify values for some or all of the coefficients of the variance function. See ?varIdent. The actual error variance will be varFunc() * sigma^2, where sigma^2 is estimated.
That's the problem. As happens in meta-analysis as well, the problem is to estimate a model with a variance component fixed. Not fixed up to a scale parameter. Fixed. In meta-analysis the model is that within each trial a treatment effect parameter is constant, and as the trial is large the variance of the estimated treatment effect is very accurately known conditional on the true treatment effect for that trial. The unconditional variance is then the known conditional variance plus an unknown variance. It doesn't seem that lme() is designed for this, and last time I tried to do it I gave up and changed the model more or less as you suggest. -thomas -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html Send "info", "help", or "[un]subscribe" (in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
If I understand your request correctly, you want to use something like "weights=varIdent(...)" as an argument to lme(). varIdent and the other varFunc constructors have an argument "fixed" that allow you to specify values for some or all of the coefficients of the variance function. See ?varIdent. The actual error variance will be varFunc() * sigma^2, where sigma^2 is estimated.
That's the problem. As happens in meta-analysis as well, the problem is to estimate a model with a variance component fixed. Not fixed up to a scale parameter. Fixed. In meta-analysis the model is that within each trial a treatment effect parameter is constant, and as the trial is large the variance of the estimated treatment effect is very accurately known conditional on the true treatment effect for that trial. The unconditional variance is then the known conditional variance plus an unknown variance. It doesn't seem that lme() is designed for this, and last time I tried to do it I gave up and changed the model more or less as you suggest.
Thanks to you both for your helpful input. Thomas is correct that I want to treat them as fixed (no free parameters) rather than fixed up to a scale factor. Indeed I have not yet been able to keep lme() from estimating two variance components. If Thomas gave up, that's pretty much enough evidence for me that I need to look elsewhere. I do not think that for this particular problem, writing a routine to do the REML/ML estimation of the inter-study variance component would be too difficult. I will add this to the growing list of items I would like to contribute to this great project. Thanks again and best regards, J.R. Lockwood 412-683-2300 x4941 lockwood at rand.org http://www.rand.org/methodology/stat/members/lockwood/ -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html Send "info", "help", or "[un]subscribe" (in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._