通过 RPy 在 lme4.lmer 输出上调用 anova

Question

我正在尝试通过RPy分析用lme4.lmer（）生成的一组线性模型的偏差。这里的这个笔记本展示了一个完整的示例，我导入我的 deps、加载我的文件、运行我的 lme4.lmer（） 并且无法让 anova 在它们上运行。

为了您的方便，这里再次粘贴了失败的行，我希望看到它的工作。

compare = stats.anova(res[0], res[1], res[2])
    Error in Ops.data.frame(data, data[[1]]) : 
      list of length 3 not meaningful
    In addition: Warning message:
    In anova.merMod(<S4 object of class "lmerMod">, <S4 object of class "lmerMod">,  :
      failed to find unique model names, assigning generic names
    ---------------------------------------------------------------------------
    RRuntimeError                             Traceback (most recent call last)
    <ipython-input-47-fe0ffa3b55de> in <module>()
    ----> 1 compare = stats.anova(res[0], res[1], res[2])
    /usr/lib64/python2.7/site-packages/rpy2/robjects/functions.pyc in __call__(self, args, **kwargs)
         84                 v = kwargs.pop(k)
         85                 kwargs[r_k] = v
    ---> 86         return super(SignatureTranslatedFunction, self).__call__(*args, **kwargs)
    /usr/lib64/python2.7/site-packages/rpy2/robjects/functions.pyc in __call__(self, *args, **kwargs)
         33         for k, v in kwargs.iteritems():
         34             new_kwargs[k] = conversion.py2ri(v)
    ---> 35         res = super(Function, self).__call__(*new_args, **new_kwargs)
         36         res = conversion.ri2py(res)
         37         return res
    RRuntimeError: Error in Ops.data.frame(data, data[[1]]) : 
      list of length 3 not meaningful

此代码在 R 中完美运行，如下所示：

> mydata = read.csv("http://chymera.eu/data/test/r_data.csv")
> library(lme4)
Loading required package: lattice
Loading required package: Matrix
> lme1 = lme4.lmer(formula='RT~cat2 + (1|ID)', data=mydata, REML=FALSE)
Error: could not find function "lme4.lmer"
> lme1 = lmer(formula='RT~cat1 + (1|ID)', data=mydata, REML=FALSE)
> lme2 = lmer(formula='RT~cat2 + (1|ID)', data=mydata, REML=FALSE)
> anova(lme1,lme2)
> lme3 = lmer(formula='RT~cat2*cat1 + (1|ID)', data=mydata, REML=FALSE)
> stats::anova(lme1, lme2, lme3)
Data: mydata
Models:
lme1: RT ~ cat1 + (1 | ID)
lme2: RT ~ cat2 + (1 | ID)
lme3: RT ~ cat2 * cat1 + (1 | ID)
     Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
lme1  4 116.68 122.29 -54.342   108.68                             
lme2  4 149.59 155.19 -70.793   141.59  0.000      0          1    
lme3  6 117.19 125.59 -52.594   105.19 36.398      2  1.248e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

你能帮我让它也运行在 RPy 中吗？

Answer 1

在 R 中时，stats::anova()大概是从函数调用中未计算的表达式推断模型名称。这里是lme1、'lme2和lme3。

现在考虑在不使用变量名称的情况下重写 R 代码，因为这更接近当前实现中发生的情况，rpy2因为数据帧和拟合的模型未绑定到变量名称。这将给出以下内容（注意："更接近"而不是"等于" - 有关此的详细信息只会分散对要点的注意力）：

stats::anova(lmer(formula='RT~cat1 + (1|ID)',
                  data=read.csv("http://chymera.eu/data/test/r_data.csv"),
                  REML=FALSE),
             lmer(formula='RT~cat2 + (1|ID)',
                  data=read.csv("http://chymera.eu/data/test/r_data.csv"),
                  REML=FALSE),
             lmer(formula='RT~cat2*cat1 + (1|ID)',
                  data=read.csv("http://chymera.eu/data/test/r_data.csv"),
                  REML=FALSE))

结果是 R 中的错误。

Error in names(mods) <- sub("@env$", "", mNms) : 
  'names' attribute [6] must be the same length as the vector [3]
In addition: Warning message:
In anova.merMod(lmer(formula = "RT~cat1 + (1|ID)", data = read.csv("http://chymera.eu/data/test/r_data.csv"),  :
  failed to find unique model names, assigning generic names

这表明 R 函数lme4:::anova.meMod正在做出很容易被违反的假设，应该通知包的作者。

它还显示表达式将用于在生成的文本输出中标识模型。

以下内容可能缺乏一点优雅，但应该既是一种解决方法，也是一种保持模型标签简短的方法。

# bind the DataFrame to an R symbol
robjects.globalenv['dataf'] = dfr
# build models, letting R fetch the symbol `dataf` when it is evaluating
# the parameters in the function call
res = list()
for formula in formulae:
    lme_res = lme4.lmer(formula=formula, data=base.as_symbol("dataf"), REML='false')
    res.append(lme_res)
# This is enough to work around the problem
compare = stats.anova(res[0], res[1], res[2])
# if not happy with the model names displayed by `compare`,
# globalenv can be filled further
names = list()
for i, value in enumerate(res):
    names.append('lme%i'  % i)
    robjects.globalenv[names[i]] = value 
# call `anova`
compare = stats.anova(*[base.as_symbol(x) for x in names])

Answer 2

这是merMod对象的anova方法中的一个错误：它本质上是由于传递给 R 的对象名称太长引起的，因此当deparse() d 时，它们最终成为（意外地）具有多个元素的字符向量。 这是由 https://github.com/lme4/lme4/commit/075c78d128db9d8398f43474621e49f32fdb5bd1 修复的;现在还有一个（未记录的）参数model.names可以指定它来覆盖模型名称的解析。

您可以使用 devtools::install_github("lme4","lme4") 安装开发版本，否则您可能需要等待一段时间才能发布修补版本......除了构建调用之外，想不出解决方法，以便传递给 R 的语言对象在解析时更短......