[R] lme() with known level-one variances (SOLUTION)

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