0

我有一个函数可以在多元线性模型中对协方差矩阵的相等性进行 Box 的 M 检验。我想把它变成一个带有公式方法的S3泛型函数，这是最自然的接口。

完整的当前代码位于https://gist.github.com/friendly/749b5a69a067e02b87dd。我可以把它全部粘贴在这里，但也许那个链接就足够了。

我不了解访问模型对象组件的函数中使用的很多魔法。leveneTest我将在包中找到的代码用作模板car，它解决了单变量模型的类似问题。

这是使用默认方法的快速测试boxM.default：

data(iris)
res <- boxM(iris[, 1:4], iris[, "Species"])
res

这给出了预期的结果：

>     data(iris)
>     res <- boxM(iris[, 1:4], iris[, "Species"])
>     res

        Box's M-test for Homogeneity of Covariance Matrices

data:  iris[, 1:4]
Chi-Sq (approx.) = 140.94, df = 20, p-value < 2.2e-16
>

当我尝试boxM.formula直接调用公式方法时，它也可以工作，给出与上面相同的输出。

boxM( cbind(Sepal.Length, Sepal.Width, Petal.Length, Petal.Width) ~ Species, data=iris)

但是，该boxM.lm方法的此测试失败：

> iris.mod <- lm(cbind(Sepal.Length, Sepal.Width, Petal.Length, Petal.Width) ~ Species, data=iris)
> boxM(iris.mod)
Error in cbind(Sepal.Length, Sepal.Width, Petal.Length, Petal.Width) : 
  object 'Sepal.Length' not found
> traceback()
8: cbind(Sepal.Length, Sepal.Width, Petal.Length, Petal.Width)
7: eval(expr, envir, enclos)
6: eval(predvars, data, env)
5: model.frame.default(form, data)
4: model.frame(form, data) at boxM.R#59
3: boxM.formula(formula(y), data = model.frame(y), ...) at boxM.R#76
2: boxM.lm(iris.mod) at boxM.R#2
1: boxM(iris.mod)
>

我想我理解它为什么会失败——这与在 中查找变量的环境有关model.frame()，但与如何纠正它无关。

有人可以帮忙吗？

4

2 回答 2

1

您设计的boxM函数可以将lm对象作为输入。该实现尝试从 中提取公式和 model.framelm并将它们与boxM.formula.

似乎这没有成功的原因是因为model.frame(iris.mod)不返回原始data.frame数据，而是返回一个 2 列 data.frame，其中第一列包含左侧变量的矩阵，第二列包含右侧的向量. 您可以通过以下方式检查

class(model.frame(iris.mod))
dim(model.frame(iris.mod))
names(model.frame(iris.mod))
model.frame(iris.mod)[,1]
model.frame(iris.mod)[,2]

由于model.frame(iris.mod)已经将数据解析为可计算格式，因此您可以应用boxM.default而不是在输入对象boxM.formula时应用lm。例如，这似乎有效：

boxM.default(Y = model.frame(iris.mod)[,1], 
             group = model.frame(iris.mod)[,2])

#    Box's M-test for Homogeneity of Covariance Matrices

#data:  model.frame(iris.mod)[, 1]
#Chi-Sq (approx.) = 140.94, df = 20, p-value < 2.2e-16
于 2016-02-10T16:59:30.053 回答
0

解决这个问题的同事说，“你被不标准的评估咬了。”

这是一个可行的解决方案，并且更普遍地符合模型对象的 S3 方法。它data在模型公式的环境中找到。

boxM.lm <- function(y, ...) {
  data <- getCall(y)$data
  y <- if (!is.null(data)) {
    data <- eval(data, envir = environment(formula(y)))
    update(y, formula(y), data = data)
  }
  else update(y, formula(y))

  boxM.formula(formula(y), data=eval(data, envir = environment(formula(y))), ...)
}
于 2016-02-10T17:47:26.823 回答