在R 3.0.+中使用.Call()出现错误

Question

我正在尝试使用一个函数来模拟鸟类巢穴成功率，该函数使用逻辑暴露链接函数。当我在R 3.0.0或3.0.1中使用上面的示例代码运行此函数时，会出现错误：

Error in .Call("logit_mu_eta", eta, PACKAGE = "stats") : 
  "logit_mu_eta" not available for .Call() for package "stats"

然而，在R 2.15.3中它运行良好。

我希望它在更新的R版本中也能够运行良好，因为我需要使用这些版本来进一步分析输出。如果有人有任何建议、解决方案或更正意见，我将非常乐意尝试。

- MAHE

Answer 1

这并不是技术上的漏洞，因为该函数使用了一个现在已经改变位置的内部函数。我在https://rpubs.com/bbolker/logregexp上发布了一个可行的示例...关键在于将logit_mu_eta更改为stats:::C_logit_mu_eta，如下所示。当然，这仍然会对内部的未来更改产生不稳定性...

logexp <- function(exposure = 1)
{
    linkfun <- function(mu) qlogis(mu^(1/exposure))
    ## FIXME: is there some trick we can play here to allow
    ##   evaluation in the context of the 'data' argument?
    linkinv <- function(eta)  plogis(eta)^exposure
    mu.eta <- function(eta) exposure * plogis(eta)^(exposure-1) *
      .Call(stats:::C_logit_mu_eta, eta, PACKAGE = "stats")
    valideta <- function(eta) TRUE
    link <- paste("logexp(", deparse(substitute(exposure)), ")",
                   sep="")
    structure(list(linkfun = linkfun, linkinv = linkinv,
                   mu.eta = mu.eta, valideta = valideta, 
                   name = link),
              class = "link-glm")
}

- Ben Bolker

---

使用 SAS 中的聊天数据，截距为2.6973。使用旧的 R 函数时，截距也是2.6973。然而，使用这个新的 R 函数和 parastat + patsize 时，截距为2.747。我明天会再看一下。 - Mark Miller

如果我只是将.Call(stats:::C_logit_mu_eta, eta, PACKAGE = "stats")插入旧函数中，仍然会得到2.747。我想知道这是否意味着C_logit_mu_eta已经改变了？但我需要明天再回来看看。 - Mark Miller

你能否发布聊天数据的链接（或其他可重现的示例）？上面的链接已经失效了... - Ben Bolker

谢谢，Ben。这是我找到数据的链接：http://zhenglei-gao.github.io/StatRMemo/posts/logexp_exposure.html 问题似乎是有人重新定义了分类变量的级别。这导致R使用不同的默认级别，从而导致参数估计值不同（截距不同，斜率符号相反）。很抱歉在发布评论之前没有发现这样一个微不足道的问题。现在我已经使用您的函数使R返回与SAS相同的估计值。 - Mark Miller

我刚刚注意到你帖子中的链接已经失效了。我只是好奇想看一下你提到的工作示例。不管怎样，再次感谢你。 - Mark Miller

Answer 2

这很奇怪。可能是一个bug。你可以将此发送到R邮件列表中，看看他们的想法。作为一个粗略的解决方法，你可以在R中重新编写该函数。以下是它的C代码，来自src/library/stats/src/family.c文件：

SEXP logit_mu_eta(SEXP eta)
{
    SEXP ans = PROTECT(duplicate(eta));
    int i, n = LENGTH(eta);
    double *rans = REAL(ans), *reta = REAL(eta);

    if (!n || !isReal(eta))
    error(_("Argument %s must be a nonempty numeric vector"), "eta");
    for (i = 0; i < n; i++) {
    double etai = reta[i];
    double opexp = 1 + exp(etai);

    rans[i] = (etai > THRESH || etai < MTHRESH) ? DOUBLE_EPS :
        exp(etai)/(opexp * opexp);
    }
    UNPROTECT(1);
    return ans;
}

THRESH 被定义为 30。因此，看起来您可以用这个函数替换外部调用：

logit_mu_eta<-function(x){
  ex<-exp(x)
  ans<-ex/(1+ex)^2
  ans[abs(x)>30]<-.Machine$double.eps
  ans
}

接着，你需要修改调用该函数的任何函数。

- mrip