2
$\begingroup$

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?

$\endgroup$
2
  • 1
    $\begingroup$ 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. $\endgroup$
    – Glen_b
    Commented Apr 22, 2014 at 8:25
  • $\begingroup$ Thanks Glen_b! The book gave me a clear idea of what the function does. $\endgroup$
    – Eduardo
    Commented Apr 28, 2014 at 18:12

1 Answer 1

1
$\begingroup$

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>
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.