如何在R中使用Box-Cox幂变换

Question

我需要将一些数据转换为“正常形状”，而我读到 Box-Cox 可以确定用于转换数据的指数。

据我所了解

car::boxCoxVariable(y)

在线性模型中，被用作响应变量。

MASS::boxcox(object)

对于公式或拟合模型对象。由于我的数据是数据框的变量，我发现唯一可以使用的函数是：

car::powerTransform(dataframe$variable, family="bcPower")

这正确吗？或者我漏掉了什么吗？

第二个问题是关于在我获得之后该做什么。

Estimated transformation parameters
dataframe$variable
0.6394806

我应该简单地将变量乘以这个值吗？ 我已经这样做了：

aaa = 0.6394806
dataframe$variable2 = (dataframe$variable)*aaa

然后我运行了Shapiro-Wilks正态性检验，但是我的数据似乎仍不符合正态分布：

shapiro.test(dataframe$variable2)
data:  dataframe$variable2
W = 0.97508, p-value < 2.2e-16

- dede

---

我发现手册中的数据转换章节提供了清晰的R代码和示例，对于其他转换也同样适用。该手册名为用R进行扩展计划评估的总结与分析。 - Valentin_Ștefan

@Valentin在提到的书中提供了非常好的解释。非常感谢！ - Rohit parihar

Answer 1

Box和Cox（1964）提出了一系列变换，旨在减少线性模型中误差的非正态性。这样做通常也会减少非线性。

这里有一份原始工作及其后续工作的简要概述：http://www.ime.usp.br/~abe/lista/pdfm9cJKUmFZp.pdf

但是您会注意到，控制lambda幂变换选择的对数似然函数依赖于基础模型的残差平方和（SO上没有LaTeX--请参见参考文献），因此不能在没有模型的情况下应用任何变换。

典型的应用如下：

library(MASS)

# generate some data
set.seed(1)
n <- 100
x <- runif(n, 1, 5)
y <- x^3 + rnorm(n)

# run a linear model
m <- lm(y ~ x)

# run the box-cox transformation
bc <- boxcox(y ~ x)

(lambda <- bc$x[which.max(bc$y)])
[1] 0.4242424

powerTransform <- function(y, lambda1, lambda2 = NULL, method = "boxcox") {

  boxcoxTrans <- function(x, lam1, lam2 = NULL) {

    # if we set lambda2 to zero, it becomes the one parameter transformation
    lam2 <- ifelse(is.null(lam2), 0, lam2)

    if (lam1 == 0L) {
      log(y + lam2)
    } else {
      (((y + lam2)^lam1) - 1) / lam1
    }
  }

  switch(method
         , boxcox = boxcoxTrans(y, lambda1, lambda2)
         , tukey = y^lambda1
  )
}


# re-run with transformation
mnew <- lm(powerTransform(y, lambda) ~ x)

# QQ-plot
op <- par(pty = "s", mfrow = c(1, 2))
qqnorm(m$residuals); qqline(m$residuals)
qqnorm(mnew$residuals); qqline(mnew$residuals)
par(op)

正如您所看到的，这不是万能解决方案——只有一些数据可以有效地转换（通常λ小于-2或大于2是不应使用该方法的信号）。与任何统计方法一样，在实施之前请谨慎使用。

要使用两个参数的Box-Cox变换，请使用geoR包找到λ值：

library("geoR")
bc2 <- boxcoxfit(x, y, lambda2 = TRUE)

lambda1 <- bc2$lambda[1]
lambda2 <- bc2$lambda[2]

修改：根据@Yui-Shiuan指出的，Tukey和Box-Cox实现的混淆已经修复。

- mlegge

---

你可能会指出，使用MASS中的模型lm(y ~ 1)（在这种情况下，bc <- boxcox(variable ~ 1, data=dataframe)）也可以得到相同的答案。 powerTransform()给出了“正确”的λ值，但数据中存在一些问题，使得仅使用Box-Cox无法强制正态性。 - Sam Dickson

很棒的回答！我可以问一下为什么你在回答开头强调了“错误”吗？这是因为一个人应该转换响应变量（例如 y 在 y ~ x_1 + x_2 中），而不是协变量（x_1 或 x_2），还是协变量也可以进行转换？ - Helen

Answer 2

根据文献《Box,George E. P.; Cox,D.R.(1964). "An analysis of transformations"》中的Box-cox变换公式，我认为mlegge的帖子需要稍作修改。转换后的y应为(y^(lambda)-1)/lambda而不是y^(lambda)。（实际上，y^(lambda)被称为Tukey变换，它是另一种不同的变换公式。） 因此，代码应该为：

(trans <- bc$x[which.max(bc$y)])
[1] 0.4242424
# re-run with transformation
mnew <- lm(((y^trans-1)/trans) ~ x) # Instead of mnew <- lm(y^trans ~ x) 

更多信息

• R中boxcox()函数正确实现Box-Cox变换公式：
  https://www.r-bloggers.com/on-box-cox-transform-in-regression-models/

• Box-Cox变换和Tukey变换的绝佳比较：http://onlinestatbook.com/2/transformations/box-cox.html

• 在维基百科上也可以找到Box-Cox变换公式： en.wikipedia.org/wiki/Power_transform#Box.E2.80.93Cox_transformation

如果我有误，请纠正我。

- Sean Yun-Shiuan Chuang

---

谢谢您指出这一点（并提供了优秀的文档！）。我已经更新了我的答案，试图解决这个问题。 - mlegge

Answer 3

如果我只想转换响应变量y，而不是指定x的线性模型，例如我想要转换/规范化数据列表，我可以将x设置为1，然后该对象就成为一个线性模型：

library(MASS)
y = rf(500,30,30)
hist(y,breaks = 12)
result = boxcox(y~1, lambda = seq(-5,5,0.5))
mylambda = result$x[which.max(result$y)]
mylambda
y2 = (y^mylambda-1)/mylambda
hist(y2)

- Hongzhen Tian

---

有4个不同的变量，所有4个直方图都显示为非正态分布，这个解决方案帮助我将它们单独转换为正态分布。 - Arvind Reddy

Answer 4

目前可以使用geoR软件包，对数据应用BoxCox变换而无需任何底层模型。具体来说，您可以使用boxcoxfit()函数找到最佳参数，然后使用BCtransform()函数预测转换后的变量。

- Mike.forest