Function to compute the ML/REML log-likelihood for a standard random-effects model.

ll.r

ll <- function(mu, tau2, yi, vi, method="ML") {
   k <- length(yi)
   wi <- 1/(vi + tau2)
   if (method == "ML") {
      sum(dnorm(yi, mean=mu, sqrt(vi + tau2), log=TRUE))
   } else {
      -1/2 * (k-1) * log(2*pi) + 1/2 * log(k) - 1/2 * sum(log(vi + tau2)) - 1/2 * log(sum(wi)) - 1/2 * sum(wi * (yi - mu)^2)
   }
}