# What exactly does the 'boxcox' function in R do?

I am familiar with the power transform family and I know how to estimate the MLE for $\lambda$ for given samples of a random variable. I have been using the 'boxcox' function in R for a sample of a random variable that depends linearly on a latent variable (i.e. a linear model) and it works very well. The problem is that I don't quite know how it works on such linear models.

The official reference won't give clear clues on how to answer my question. Can you provide some information on the answer to my question?

• Besides Jeromy's nice answer, you might find a few additional details in the book that MASS is the package for (Modern Applied Statistics with S, Venables and Ripley). If I was close to my copy, I'd check it now. Commented Apr 22, 2014 at 8:25
• Thanks Glen_b! The book gave me a clear idea of what the function does. Commented Apr 28, 2014 at 18:12

This is not a complete answer, but the real official reference is arguably "the code". In this case it's a little harder to find (it's a method, and it's hidden behind a namespace), but this command shows the code for the default method.

library('MASS')
getAnywhere('boxcox.default')


The result is the following (see also getAnywhere('boxcox.lm')):

A single object matching ‘boxcox.default’ was found
It was found in the following places
registered S3 method for boxcox from namespace MASS
namespace:MASS
with value

function (object, lambda = seq(-2, 2, 1/10), plotit = TRUE, interp = (plotit &&
(m < 100)), eps = 1/50, xlab = expression(lambda), ylab = "log-Likelihood",
...)
{
if (is.null(y <- object$y) || is.null(xqr <- object$qr))
stop(gettextf("%s does not have both 'qr' and 'y' components",
sQuote(deparse(substitute(object)))), domain = NA)
if (any(y <= 0))
stop("response variable must be positive")
n <- length(y)
y <- y/exp(mean(log(y)))
logy <- log(y)
xl <- loglik <- as.vector(lambda)
m <- length(xl)
for (i in 1L:m) {
if (abs(la <- xl[i]) > eps)
yt <- (y^la - 1)/la
else yt <- logy * (1 + (la * logy)/2 * (1 + (la * logy)/3 *
(1 + (la * logy)/4)))
loglik[i] <- -n/2 * log(sum(qr.resid(xqr, yt)^2))
}
if (interp) {
sp <- spline(xl, loglik, n = 100)
xl <- sp$x loglik <- sp$y
m <- length(xl)
}
if (plotit) {
mx <- (1L:m)[loglik == max(loglik)][1L]
Lmax <- loglik[mx]
lim <- Lmax - qchisq(19/20, 1)/2
dev.hold()
on.exit(dev.flush())
plot(xl, loglik, xlab = xlab, ylab = ylab, type = "l",
ylim = range(loglik, lim))
plims <- par("usr")
abline(h = lim, lty = 3)
y0 <- plims[3L]
scal <- (1/10 * (plims[4L] - y0))/par("pin")[2L]
scx <- (1/10 * (plims[2L] - plims[1L]))/par("pin")[1L]
text(xl[1L] + scx, lim + scal, " 95%", xpd = TRUE)
la <- xl[mx]
if (mx > 1 && mx < m)
segments(la, y0, la, Lmax, lty = 3)
ind <- range((1L:m)[loglik > lim])
if (loglik[1L] < lim) {
i <- ind[1L]
x <- xl[i - 1] + ((lim - loglik[i - 1]) * (xl[i] -
xl[i - 1]))/(loglik[i] - loglik[i - 1])
segments(x, y0, x, lim, lty = 3)
}
if (loglik[m] < lim) {
i <- ind[2L] + 1
x <- xl[i - 1] + ((lim - loglik[i - 1]) * (xl[i] -
xl[i - 1]))/(loglik[i] - loglik[i - 1])
segments(x, y0, x, lim, lty = 3)
}
}
list(x = xl, y = loglik)
}
<bytecode: 0x7ff6e0973120>
<environment: namespace:MASS>