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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
See the rmeta package on CRAN (the metaDSL() function would do what you are looking for). BTW, in your code note: i)In the model estimation, the estimated coefficients beta have to be weighed by the inverse of their variance (the smaller the variance, the more accurate the estimated beta, the more important its influence on the pooled estimate) ii) As far as I know, it makes no sense to set no.group=no.observations; otherwise you could not estimate the intra-group variance. best, vito ----- Original Message ----- From: "J.R. Lockwood" <lockwood at rand.org> To: <r-help at stat.math.ethz.ch> Cc: "J.R. Lockwood" <lockwood at rand.org> Sent: Thursday, August 29, 2002 8:00 PM 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 -- Readhttp://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 _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
Dear Vito, Thank you for your response. I still have some concerns:> > See the rmeta package on CRAN (the metaDSL() function would do what you are > looking for).It is not clear to me that this function does what I need. I apologize for casting my original question in terms of a meta-analysis; it is a bit of a red herring. My more general question is whether it is possible to tell lme() to treat some variance components as known, and to estimate others conditional on them.> > BTW, in your code note: > > i)In the model estimation, the estimated coefficients beta have to be > weighed by the inverse of their variance (the smaller the variance, the more > accurate the estimated beta, the more important its influence on the pooled > estimate)Unless I am misinterpreting, the "weights" argument that I passed to lme() is inverted. But after further consideration I do not think that weighting is what I am seeking.> ii) As far as I know, it makes no sense to set no.group=no.observations; > otherwise you could not estimate the intra-group variance. >I am not sure to what function the arguments you reference belong (they do not appear to be part of meta.DSL()). But the larger point is that I do not want to estimate the intra-group variances -- I want to treat them as known. best, 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 _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
Setzer.Woodrow@epamail.epa.gov
2002-Aug-30 15:35 UTC
[R] lme() with known level-one variances
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
Setzer.Woodrow@epamail.epa.gov
2002-Aug-30 16:28 UTC
[R] lme() with known level-one variances
Right. After thinking about the problem a few minutes longer, I
realized my mistake. Chapter 5 ("Inference Based on Individual
Estimates") of Davidian and Giltinan (1995) : "Nonlinear Models for
Repeated Measurement Data" can be used to do exactly this. It turns out
to be quite simple for a one-dimensional parameter, and I coded it up.
This is my interpretation of their "Global two-stage method":
"estmixed" <-
function (mi, sei)
{
if (length(mi) == 1) {
c(muhat = mi, muhat.se = sei, sdhat = NA)
}
else {
vari <- sei^2
start <- c(mean(mi), log(var(mi)))
out <- try(optim(par = start, fn = mll, method = "BFGS",
mi = mi, vari = vari))
if (inherits(out, "try-error")) {
print(cbind(mi, sei))
print(out)
stop("Problem in optim")
}
if (out$convergence == 0)
c(muhat = out$par[1], muhat.se = sqrt(1/sum(1/(vari +
exp(out$par[2])))), sdhat = sqrt(exp(out$par[2])))
else c(muhat = NA, muhat.se = NA, sdhat = NA)
}
}
On input, mi is the vector of parameter estimates, and sei is the vector of
their standard errors.
On output, muhat is the global estimate, muhat.se is an estimate of its standard
error, and sdhat is an estimate of the among-group standard
deviation.
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
Thomas Lumley
<tlumley at u.washin To: Woodrow
Setzer/RTP/USEPA/US at EPA
gton.edu> cc: "J.R.
Lockwood" <lockwood at rand.org>, r-help at stat.math.ethz.ch
Subject: Re: [R] lme() with
known level-one variances
08/30/02 12:11 PM
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
Can you do this using gls() with control = glsControl( sigma = yourValue ) ? --Matt -----Original Message----- From: Thomas Lumley [mailto:tlumley at u.washington.edu] Sent: Friday, August 30, 2002 9:12 AM To: Setzer.Woodrow at epamail.epa.gov Cc: J.R. Lockwood; r-help at stat.math.ethz.ch Subject: Re: [R] lme() with known level-one variances 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 _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._. _._ -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- 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 _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
Setzer.Woodrow@epamail.epa.gov
2002-Aug-30 18:34 UTC
[R] lme() with known level-one variances
That argument is not documented in the current Windows version of nlme, and gives an error when I try it.:> glsControl(sigma=1)Error in glsControl(sigma = 1) : unused argument(s) (sigma ...)> package.description("nlme")[c("Version","Date")]Version Date "3.1-29" "2002/08/10"> version_ platform i386-pc-mingw32 arch i386 os mingw32 system i386, mingw32 status major 1 minor 5.1 year 2002 month 06 day 17 language R 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 "Austin, Matt" <maustin at amgen.co To: 'Thomas Lumley' <tlumley at u.washington.edu>, Woodrow Setzer/RTP/USEPA/US at EPA m> cc: "J.R. Lockwood" <lockwood at rand.org>, r-help at stat.math.ethz.ch Subject: RE: [R] lme() with known level-one variances 08/30/02 01:56 PM Can you do this using gls() with control = glsControl( sigma = yourValue ) ? --Matt -----Original Message----- From: Thomas Lumley [mailto:tlumley at u.washington.edu] Sent: Friday, August 30, 2002 9:12 AM To: Setzer.Woodrow at epamail.epa.gov Cc: J.R. Lockwood; r-help at stat.math.ethz.ch Subject: Re: [R] lme() with known level-one variances 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 theother> varFunc constructors have an argument "fixed" that allow you tospecify> 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 _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._. _._. _._ -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- 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 _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
I apologize for the confusion, you are correct this is not implemented in the current R version of nlme, it's in the current S-Plus version. --Matt -----Original Message----- From: Setzer.Woodrow at epamail.epa.gov [mailto:Setzer.Woodrow at epamail.epa.gov] Sent: Friday, August 30, 2002 11:35 AM To: Austin, Matt Cc: r-help at stat.math.ethz.ch Subject: RE: [R] lme() with known level-one variances That argument is not documented in the current Windows version of nlme, and gives an error when I try it.:> glsControl(sigma=1)Error in glsControl(sigma = 1) : unused argument(s) (sigma ...)> package.description("nlme")[c("Version","Date")]Version Date "3.1-29" "2002/08/10"> version_ platform i386-pc-mingw32 arch i386 os mingw32 system i386, mingw32 status major 1 minor 5.1 year 2002 month 06 day 17 language R 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 "Austin, Matt" <maustin at amgen.co To: 'Thomas Lumley' <tlumley at u.washington.edu>, Woodrow Setzer/RTP/USEPA/US at EPA m> cc: "J.R. Lockwood" <lockwood at rand.org>, r-help at stat.math.ethz.ch Subject: RE: [R] lme() with known level-one variances 08/30/02 01:56 PM Can you do this using gls() with control = glsControl( sigma = yourValue ) ? --Matt -----Original Message----- From: Thomas Lumley [mailto:tlumley at u.washington.edu] Sent: Friday, August 30, 2002 9:12 AM To: Setzer.Woodrow at epamail.epa.gov Cc: J.R. Lockwood; r-help at stat.math.ethz.ch Subject: Re: [R] lme() with known level-one variances 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 theother> varFunc constructors have an argument "fixed" that allow you tospecify> 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 _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._. _._. _._ -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- 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 _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
Thanks to the help of R. Woodrow Setzer, Jr. and Thomas Lumley, the
problem has been solved. A brief recap:
The model of interest is
\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. We
would like to estimate \tau^2 and \mu, using either ML or REML. The
motivating application was a meta-analysis of regression coefficients
(the \beta_{i}) with the fixed variances equal to the squared standard
errors reported in the regression output.
The conclusion was that lme() is not suitable for this estimation, or
at least it is not obvious how to trick lme() into doing it. The
solution is that the problem is easy enough to attack with brute
force. I modified the code very kindly posted by Mr. Setzer to handle
either ML or REML estimation; it is below.
-------------------------------------------
estmixed <- function (mi, sei, method="ML"){
## mi are individual estimates, sei are their standard errors
mll<-function (par, mi, vari){
## calculate -2 * log likelihood
## par[1] is the grand mean
## par[2] is the log of the between-group variance component
sum(log(vari + exp(par[2]))) + sum((mi - par[1])^2/(vari +
exp(par[2])))
}
mll.reml<-function (par, mi, vari){
## calculate -2 * log REML likelihood (see Lindstrom and Bates)
sum(log(vari + exp(par[2]))) + sum((mi - par[1])^2/(vari +
exp(par[2]))) + log(sum(1/(vari + exp(par[2]))))
}
if(method!="ML" & method!="REML"){
stop("Invalid value of method")
}
if (length(mi) == 1) {
c(muhat = mi, muhat.se = sei, sdhat = NA)
}
else {
vari <- sei^2
start <- c(mean(mi), log(var(mi)))
objfun <- if(method == "ML") mll else mll.reml
out <- try(optim(par = start, fn = objfun, method
= "BFGS",control=list(maxit=500),mi = mi, vari = vari))
if (inherits(out, "try-error")) {
print(cbind(mi, sei))
print(out)
stop("Problem in optim")
}
if (out$convergence == 0)
c(muhat = out$par[1], muhat.se = sqrt(1/sum(1/(vari +
exp(out$par[2])))), sdhat = sqrt(exp(out$par[2])))
else {
c(muhat = NA, muhat.se = NA, sdhat = NA)
}
}
}
-------------------------------------------------------
An example application follows:
beta.est<-c(0.40,-8.02,-1.60,-8.50,-3.20)
beta.se<-c(3.64,3.73,2.35,3.14,2.05)
estmixed(beta.est,beta.se,"ML")
## -3.6835705 1.2241980 0.0700697
estmixed(beta.est,beta.se,"REML")
## -3.806156 1.428062 1.513580
Thanks again--
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._