R中的apply和boxcox函数错误答案

【问题标题】：apply and boxcox function error in RR中的apply和boxcox函数错误
【发布时间】：2015-08-10 13:53:05
【问题描述】：

我在调用 boxcox 函数时遇到了使用 apply 函数的问题。

该数据集包含 35 个样本，2 天内 3 次治疗（有些人在几天内重复，但总是在治疗组内），有大约 1000 个测量值/因变量，我在这里只包括了 5 个（#1,2,4， 7,9)。

NUM day trt 1   2   4   7   9  
8145    7   L   0.986   0.423   0.86    0.648   2.031  
8169    7   L   1.013   0.512   1.157   0.633   0.55  
8201    7   L   0.236   6.144   0.604   1.759   1.181  
8212    7   L   0.455   0.707   0.693   0.972   0.615  
8168    7   L   1.261   0.618   1.138   0.943   0.868  
8193    7   L   1.754   0.273   1.224   0.719   0.895  
8145    13  L   1.257   0.36    1.505   0.729   0.735  
8169    13  L   1.549   0.577   1.53    0.771   1.116  
8201    13  L   0.489   0.378   0.842   1.538   1.676  
8212    13  L   0.991   0.34    1.03    1.157   1.076  
8168    13  L   1.355   0.391   1.416   0.953   1.479  
8193    13  L   1.029   0.308   1.027   0.902   1.934  
8214    7   C   1.224   0.298   1.113   1.445   0.218  
8139    7   C   1.277   0.554   1.443   0.895   0.74  
8151    7   C   1.312   2.025   1.197   0.675   0.791  
8160    7   C   1.555   0.405   1.432   0.826   0.501  
8196    7   C   0.938   0.717   0.917   0.801   1.462  
8213    7   C   0.835   1.863   0.942   1.967   0.739  
8139    13  C   0.958   0.275   1.273   1.351   0.842  
8151    13  C   0.864   0.517   0.98    1.368   1.865  
8160    13  C   1.516   0.895   1.318   0.551   0.779  
8239    13  C   1.071   0.194   0.955   1.87    0.68  
8196    13  C   1.299   0.594   1.14    0.877   1.873  
8213    13  C   1.601   0.733   1.375   0.738   1.273  
8231    7   H   1.401   0.483   1.001   1.052   0.548  
8232    7   H   1.292   0.574   1.634   0.641   0.464  
8219    7   H   0.785   0.396   0.886   0.903   1.734  
8177    7   H   0.525   0.252   0.563   0.914   1.174  
8143    7   H   1.398   0.266   0.947   0.94    0.781  
8219    13  H   0.903   6.225   1.109   1.218   1.073  
8143    13  H   1.086   0.435   1.4 0.922   0.925  
8167    13  H   0.83    0.574   1.09    1.338   1.563  
8231    13  H   1.498   0.482   1.375   0.855   0.719  
8232    13  H   1.055   0.811   0.887   0.606   0.634  
8177    13  H   0.857   1.324   0.954   1.635   0.675

我正在尝试在 R 的 MASS 包中运行 boxcox 函数，以评估是否需要转换数据，如果需要，则根据具有最大对数似然（或其舍入）的 lamda 值进行转换。

linear.f=function(x){lm(x~day+trt+day*trt, data=data)}
linear.multiple=apply(data[,4:ncol(data)],2,linear.f)

工作正常

单独运行 boxcox 函数也可以正常工作

boxcox(lm(x~day+trt+day*trt,data=data))$x[which.max(boxcox(lm(x~day+trt+day*trt, data=data))$y)]

将函数放入apply：

lamda.f=function(x){boxcox(lm(x~day+trt+day*trt, data=data))$x[which.max(boxcox(lm(x~day+trt+day*trt, data=data))$y)]}

但是，当尝试使用 apply 跨所有列/因变量运行 boxcox 时会出现问题：

lamda.multiple=apply(data[,4:ncol(data)],2,lamda.f)

错误：

 Error in model.frame.default(formula = x ~ day + trt + day * trt, data = data,  : 
  invalid type (list) for variable 'x'

我在做一些明显不正确的事情吗？因为 boxcox 将在每列数据上一个一个地工作..

我也尝试过使用 for 循环而不是 apply，但没有成功

任何建议或帮助将不胜感激。

【问题讨论】：

标签： r apply

【解决方案1】：

该问题与apply() 调用无关。这与您在公式中使用 x 的事实有关，但您传递给 lm() 的数据集没有 x 列。

lm() 函数有点“宽松”，也许是不可取的。如果在给定的数据集中找不到给定公式中的符号，则lm() 允许它绑定到可以在公式的闭包环境链中找到的任何变量。来自https://stat.ethz.ch/R-manual/R-devel/library/stats/html/model.frame.html：

公式、子集和 ... 中的所有变量首先在数据中查找，然后在公式环境中查找（有关详细信息，请参阅 formula() 的帮助）并收集到数据框中。

...

只有类型为原始、逻辑、整数、实数、复数或字符的变量可以包含在模型框架中：这包括分类变量，例如因子（其基础类型为整数），但不包括列表。

现在，实际上，您的linear.f() 和lamda.f() 函数都有一个函数参数x，这允许lm() 调用在两个函数中成功。 IOW，当lm() 在这些函数中运行时，x 无法绑定到给定数据集 (data) 中的任何列，但随后绑定到公式的闭包环境中的函数参数（即当前功能评估的评估环境）。因此，lm() 调用成功。

更详细地说：评估函数时，其参数始终存储在为对该特定函数进行特定评估而创建的评估环境中。由于您已在函数内部逐字定义公式，因此它围绕当前评估环境关闭，因此当 lm() 运行时，在未能绑定到数据集中的列后，搜索符号 @987654338 的目标@ 首先命中该环境并绑定到环境中的函数参数条目。

我注意到您表示lm()+boxcox() 的一次性顶级调用成功。这很奇怪，因为没有要绑定的函数参数。我只能猜测，在你运行那行代码时，你碰巧在全局环境中有一个变量x，它对模型有效。如果没有，那么它应该会失败，因为公式中的 x 符号在数据集或闭包环境链中都没有绑定。注意 lm() 应该失败；我什至不在这里谈论boxcox()：

if (exists('x')) rm(x); ## remove global x
lm(x~day+trt+day*trt,data=data);
## Error in eval(expr, envir, enclos) : object 'x' not found

现在，在您运行第二次 apply 调用时，我们可能猜到您已经用一个列表覆盖了全局 x 变量，这将导致您得到的确切错误消息：

x <- list();
lm(x~day+trt+day*trt,data=data);
## Error in model.frame.default(formula = x ~ day + trt + day * trt, data = data,  :
##   invalid type (list) for variable 'x'

但是这里肯定有其他事情发生。回想一下我是如何指出linear.f() 和lamda.f() 都有一个函数参数x，所以lm() 应该绑定到它，因为这是一个有效的预测向量，它应该会成功。

在我看来，boxcox() 正在自动调用lm()，并且在这种情况下以某种方式阻止了闭包链搜索。我相信我可以用下面的代码证明这个猜想，它使用了IIFEs（最初是一个 JavaScript 术语，但适用于任何具有一流可表达函数的语言，如 R）：

if (exists('x')) rm(x); ## remove global x
data2 <- data.frame(y=1:3); ## no x
(function(x) lm(y~x,data2))(1:3); ## x binds to function parameter
##
## Call:
## lm(formula = y ~ x, data = data2)
##
## Coefficients:
## (Intercept)            x
##           0            1
##
(function(x) boxcox(lm(y~x,data2)))(1:3);
## Error in eval(expr, envir, enclos) : object 'x' not found
traceback();
## 14: eval(expr, envir, enclos)
## 13: eval(predvars, data, env)
## 12: model.frame.default(formula = y ~ x, data = data2, drop.unused.levels = TRUE)
## 11: stats::model.frame(formula = y ~ x, data = data2, drop.unused.levels = TRUE)
## 10: eval(expr, envir, enclos)
## 9: eval(mf, parent.frame())
## 8: lm(formula = y ~ x, data = data2, y = TRUE, qr = TRUE)
## 7: eval(expr, envir, enclos)
## 6: eval(call, parent.frame())
## 5: update.default(object, y = TRUE, qr = TRUE, ...)
## 4: update(object, y = TRUE, qr = TRUE, ...)
## 3: boxcox.lm(lm(y ~ x, data2))
## 2: boxcox(lm(y ~ x, data2)) at #1
## 1: (function(x) boxcox(lm(y ~ x, data2)))(1:3)

如您所见，boxcox() 似乎正在运行lm()，但与直接的lm() 调用不同，它无法绑定到函数参数x，即使它应该能够绑定，因为传递给boxcox() 的公式参数包含闭包环境指针。我想，我们可以将此归结为 boxcox() 函数的一个弱点。

无论如何，就解决方案而言，一个好的方法是不依赖闭包环境系统来解析未绑定的公式符号，而是始终确保给定的数据集包含所有列公式中的符号。简单例子：

data2 <- data.frame(y=1:3,x=10:12);
boxcox(lm(y~x,data2)); ## succeeds

但是，在您的情况下，这并不是那么简单，因为您有一个动态的结果变量。一个好的方法是在将结果变量传递给 lm() 调用之前用结果变量补充您的 data.frame，例如cbind(data,x)。

不幸的是，令人惊讶的是，我收到以下错误：

lamda.f <- function(x) { data.x <- cbind(data,x); boxcox(lm(x~day+trt+day*trt,data=data.x))$x[which.max(boxcox(lm(x~day+trt+day*trt,data=data.x))$y)]; };
lamda.multiple <- apply(data[,4:ncol(data)],2,lamda.f);
## Error in is.data.frame(data) : object 'data.x' not found

我什至尝试将lm() 调用与boxcox() 调用分开，以防boxcox() 对其参数进行某种疯狂的非标准评估，从而尝试在上下文中评估data.x这将阻止与lamda.f() 评估环境的绑定：

lamda.f <- function(x) { data.x <- cbind(data,x); m <- lm(x~day+trt+day*trt,data=data.x); b <- boxcox(m); b$x[which.max(b$y)]; };
lamda.multiple <- apply(data[,4:ncol(data)],2,lamda.f);
## Error in is.data.frame(data) : object 'data.x' not found

在我看来，boxcox() 在很大程度上依赖于全球环境。在函数范围内运行 boxcox() 并依赖于任何局部变量只会破坏它。我怀疑它正在检查存储在模型对象上的调用（例如m$call）并尝试直接读取符号。无论如何，这真的很奇怪。

我认为您可以通过将补充的 data.frame 存储在全局环境中并确保它在 boxcox() 运行时存在于那里来解决这个问题。您可以为此使用超赋值运算符：

lamda.f <- function(x) { data.x <<- cbind(data,x); m <- lm(x~day+trt+day*trt,data=data.x); b <- boxcox(m); b$x[which.max(b$y)]; };
lamda.multiple <- apply(data[,4:ncol(data)],2,lamda.f);
lamda.multiple;
##         X1         X2         X4         X7         X9
##  1.2323232 -0.6666667  0.7474747 -0.6666667  0.2222222

【讨论】：