boxcox transformation method compare.(R, python)

user_13267574 0
Viewed 43

I have a simple question: I want to make a boxcox transformation on a series of data. I tried different method, then got different results. Why I don't have consistent result?

Data:(in R form)

x <- c(112,118,132,129,121,135,148,148,136,119,104,118,115,126,141,135,125,149,170,170,158,133,114,140,145,150,178,163,172,178,199,199,184,162,146,166,171,180,193,181,183,218,230,242,209,191,172,194,196,196,236,235,229,243,264,272,237,211,180,201,204,188,235,227,234,264,302,293,259,229,203,229,242,233,267,269,270,315,364,347,312,274,237,278,284,277,317,313,318,374,413,405,355,306,271,306,315,301,356,348,355,422,465,467,404,347,305,336,340,318,362,348,363,435,491,505,404,359,310,337,360,342,406,396,420,472,548,559,463,407,362,405,417,391,419,461,472,535,622,606,508,461,390,432)

Plan A: R script from boxcox definition

par2 <- 200
par3 <- 100
numlam <- par2 + par3 + 1
n <- length(x)
c <- array(NA,dim=c(numlam))
l <- array(NA,dim=c(numlam))
mx <- -1
mxli <- -999
for (i in 1:numlam)
{
  l[i] <- (i-200-1)/100
  if (l[i] != 0)
  {
    x1 <- (x^l[i] - 1) / l[i]
  } else {
    x1 <- log(x)
  }
  c[i] <- cor(qnorm(ppoints(x), mean=0, sd=1),sort(x1))
  if (mx < c[i])
  {
    mx <- c[i]
    mxli <- l[i]
    x1.best <- x1
  }
}
print(paste("best lambda:", mxli))
hist(x)
hist(x1.best)

Plan B: R package

library(MASS)
bc= boxcox(lm(x~x))
lmbd2 = bc$x[which(bc$y == max(bc$y))]
print(paste("best lambda:", lmbd2))

Plan C: python package

from scipy.stats import boxcox
import matplotlib.pyplot as plt
x = [112,118,132,129,121,135,148,148,136,119,104,118,115,126,141,135,125,149,170,170,158,133,114,140,145,150,178,163,172,178,199,199,184,162,146,166,171,180,193,181,183,218,230,242,209,191,172,194,196,196,236,235,229,243,264,272,237,211,180,201,204,188,235,227,234,264,302,293,259,229,203,229,242,233,267,269,270,315,364,347,312,274,237,278,284,277,317,313,318,374,413,405,355,306,271,306,315,301,356,348,355,422,465,467,404,347,305,336,340,318,362,348,363,435,491,505,404,359,310,337,360,342,406,396,420,472,548,559,463,407,362,405,417,391,419,461,472,535,622,606,508,461,390,432]
plt.hist(x, bins=10)
y,lmda = boxcox(x)
print(lmda)
plt.hist(y, bins=13)

I got 3 different for lambda: 0.22, 0.1414, 0.1480

r python
1 Answers

Your script optimizes correlation. Both MASS::boxcox() in R and scipy.stats.boxcox() in Python use maximum likelihood estimation, instead.

The difference between R and Python comes from the fact that MASS::boxcox() uses a fairly sparse grid search by default. By using a finer grid, we can get the results to agree:

> with(boxcox(x ~ 1, lambda = seq(0.1, 0.2, by = 0.001)), x[which.max(y)])
# [1] 0.148

With scipy you can also optimize correlation, to match your script’s result:

>>> scipy.stats.boxcox_normmax(x, method='pearsonr')
# 0.21082925698127897

Further, in R, you can also manually optimize for the log-likelihood function (following scipy):

llf <- function(x, lambda) {
  y <- (x^lambda - 1) / lambda
  n <- length(x)
  (lambda - 1) * sum(log(x)) - n / 2 * log(sum((y - mean(y))^2) / n)
}

optimize(function(lambda) -llf(x, lambda), interval = c(-2, 2))
#> $minimum
#> [1] 0.1480227
#> 
#> $objective
#> [1] 679.5431
user_45506950
Related