R help - Jul 2014 - Box-cox transformation

Hi,

I am trying to do Box-Cox transformation, but I am not sure how to do it
correctly.  Here is an example showing what I am trying:



# example from MASS

require(MASS)
boxcox(Days+1 ~ Eth*Sex*Age*Lrn, data = quine,
       lambda = seq(-0.05, 0.45, len = 20))

# Here is My attempt at getting the profile likelihood for the Box-Cox parameter
lam <- seq(-0.05, 0.45, len = 20)
dev <- rep(NA, length=20)

for (i in 1:20) {
a <- lam[i]
ans <- glm(((Days+1)^a-1)/a ~ Eth*Sex*Age*Lrn, family=gaussian, data = quine)
dev[i] <- ans$deviance
}

plot(lam, dev, type="b", xlab="lambda",
ylab="deviance")

I am trying to create the profile likelihood for the Box-Cox parameter, but
obviously I am not getting it right.  I am not sure that ans$deviance is the
right thing to do.

I would appreciate any guidance.

Thanks & Best,
Ravi




	[[alternative HTML version deleted]]

Hi Ravi,

Deviance is the SS in this case, but you need a normalizing constant
adjusted by the lambda to put them on the same scale.  I modified your
example below to simplify slightly and use the normalization (see the
LL line).

Cheers,

Josh

######################################

require(MASS)

myp <- function(y, lambda) (y^lambda-1)/lambda


lambda <- seq(-0.05, 0.45, len = 20)
N <- nrow(quine)
res <- matrix(numeric(0), nrow = length(lambda), 2, dimnames list(NULL,
c("Lambda", "LL")))

# scaling contant
C <- exp(mean(log(quine$Days+1)))

for(i in seq_along(lambda)) {
  r <- resid(lm(myp(Days + 1, lambda[i]) ~ Eth*Sex*Age*Lrn, data = quine))
  LL <- (- (N/2) * log(sum((r/(C^lambda[i]))^2)))
  res[i, ] <- c(lambda[i], LL)
}

# box cox
boxcox(Days+1 ~ Eth*Sex*Age*Lrn, data = quine, lambda = lambda)
# add our points on top to verify match
points(res[, 1], res[,2], pch = 16)

######################################



On Mon, Jul 7, 2014 at 11:33 AM, Ravi Varadhan <ravi.varadhan at jhu.edu>
wrote:
> Hi,
>
> I am trying to do Box-Cox transformation, but I am not sure how to do it
correctly.  Here is an example showing what I am trying:
>
>
>
> # example from MASS
>
> require(MASS)
> boxcox(Days+1 ~ Eth*Sex*Age*Lrn, data = quine,
>        lambda = seq(-0.05, 0.45, len = 20))
>
> # Here is My attempt at getting the profile likelihood for the Box-Cox
parameter
> lam <- seq(-0.05, 0.45, len = 20)
> dev <- rep(NA, length=20)
>
> for (i in 1:20) {
> a <- lam[i]
> ans <- glm(((Days+1)^a-1)/a ~ Eth*Sex*Age*Lrn, family=gaussian, data =
quine)
> dev[i] <- ans$deviance
> }
>
> plot(lam, dev, type="b", xlab="lambda",
ylab="deviance")
>
> I am trying to create the profile likelihood for the Box-Cox parameter, but
obviously I am not getting it right.  I am not sure that ans$deviance is the
right thing to do.
>
> I would appreciate any guidance.
>
> Thanks & Best,
> Ravi
>
>
>
>
>         [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.


-- 
Joshua F. Wiley
Ph.D. Student, UCLA Department of Psychology
http://joshuawiley.com/
Senior Analyst, Elkhart Group Ltd.
http://elkhartgroup.com
Office: 260.673.5518