Skip to content

lme() with known level-one variances

2 messages · Setzer.Woodrow@epamail.epa.gov, J.R. Lockwood

Original
Setzer.Woodrow@epamail.epa.gov
# 
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:

            

      
      
      
    

        
other

            

      
      
      
    

        
specify

            

      
      
      
    

        
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._

  
  


  


    

    

      
      

        


  

        

      J.R. Lockwood
    

    

      
      #
    

  


  

      
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._