2

生态学中的一个常见情况是具有二元结果（0 = 死亡，1 = 生存）的生存模型，其中个体（在此示例中，考虑鸟类的个体筑巢尝试）在他们面临死亡风险的天数方面有所不同。为了解决这个问题，我们使用修改后的逻辑回归，将曝光天数合并到链接函数中。
正如 Shaffer (2004) 所述：

“通过选择适当的预测函数，可以根据 x 对每日存活率进行建模，在我们的例子中，它应该产生介于 0 和 1 之间的值。正如在逻辑回归中所做的那样，我们使用 S 形逻辑函数：

在此处输入图像描述

我们的广义线性模型的系统分量是 [s(x)]t。接下来，我们考虑函数：

在此处输入图像描述

上述函数是关于 θ 的单调可微的，可以证明 g(θ) = β0 + β1x，满足广义线性模型中链接函数的标准。这三个组成部分：二项式响应分布、表达式 1 中给出的预测函数和表达式 2 中给出的链接函数，完全指定了我们的广义线性模型。该模型（以下简称“逻辑暴露模型”）类似于逻辑回归模型，但在链接函数的形式上有所不同。逻辑暴露链接函数在分子和分母中包含一个指数 (1/t)，而逻辑回归链接函数中不存在该指数。指数是必要的，以说明在区间中存活的概率取决于区间长度。”</p>

此链接函数的代码可在 Web 上找到，如果您键入“help(family)”，它也是 R 中描述的示例链接函数之一：

logexp <- function(days = 1)
{
linkfun <- function(mu) qlogis(mu^(1/days))
linkinv <- function(eta) plogis(eta)^days
mu.eta <- function(eta) days * plogis(eta)^(days-1) *
  .Call("logit_mu_eta", eta, PACKAGE = "stats")
valideta <- function(eta) TRUE
link <- paste("logexp(", days, ")", sep="")
structure(list(linkfun = linkfun, linkinv = linkinv,
               mu.eta = mu.eta, valideta = valideta, name = link),
          class = "link-glm")
}

它在这样的模型中工作得很好：

glm(survive ~ date, family=binomial(link=logexp(days=dat$Days)),data=dat)

我遇到的问题是尝试在 GLMER 模型中使用此自定义链接函数并添加随机效果时（我在网上找到了此方法的一个示例：http://rstudio-pubs-static.s3 .amazonaws.com/4082_51aa699bd9f041c7b3f7cf7b9252f60c.html）。

在我们的例子中，我们希望将站点包含为随机效应。模型的制定方式与之前的 GLM 相同：

glmer(survive ~ date + (1|site), family=binomial(link=logexp(days=dat$Days)),data=dat)

但是，现在我收到一条错误消息：

famType（glmFit$family）中的错误：未知链接：'logexp（3）'未知链接：'logexp（4）'未知链接：'logexp（3）'未知链接：'logexp（2）'未知链接：'logexp (3)'未知链接：'logexp(3)'未知链接：'logexp(4)'未知链接：'logexp(3)'未知链接：'logexp(2)'未知链接：'logexp(1)'未知链接：'logexp（4）'未知链接：'logexp（5）'未知链接：'logexp（4）'未知链接：'logexp（3）'未知链接：'logexp（4）'未知链接：'logexp（ 5)'未知链接：'logexp(3)'未知链接：'logexp(4)'未知链接：'logexp(3)'未知链接：'logexp(3)'未知链接：'logexp(3)'未知链接：'logexp（3）'未知链接：'logexp（3）'未知链接：'logexp（3）'未知链接：'logexp（3）'未知链接：'logexp（3）'未知链接：'logexp（3 )'未知链接：'logexp(2)'未知链接：'logexp(1)'未知链接：'logexp(3)'未知链接：'logexp(1)'未知链接：'logexp(1)'未知链接：'logexp(1)'unknown link: 'logexp(1)'unkn 另外：警告消息：在 if (!(lTyp <- match(family$link, linkNms, nomatch = 0))) stop(gettextf("unknown link: %s", : 条件长度 > 1 并且只使用第一个元素

错误消息为每一行数据列出了一个未知链接，其中数字对应于该巢访问（或数据行）的曝光天数。

例如：第一个 'logexp(3)' 对应于具有 3 个曝光天数的第一行数据。

有其他人能够在 GLMER 模型中使用此自定义链接功能吗？或者，如果没有，是否有人知道导致错误的原因？

＃＃＃＃＃＃更新＃＃＃＃＃＃

非常感谢 Ben Bolker 解决了这个问题。我更新到 3.0.2 和最新版本的 lme4 并使用了 Ben 的 R 相关帖子（https://www.rpubs.com/bbolker/logregexp）的链接功能，就是这个：

library(MASS)
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")
}
4

1 回答 1

3

您需要使用更新版本的lme4软件包，例如刚刚在 CRAN 上发布的 1.0-4 版本。早期版本不允许用户指定的链接功能。

另外，请注意，您在上面发布的代码不是最新版本中出现的代码?family，其中 (obsolete).Call("logit_mu_eta", eta, PACKAGE = "stats")被纯 R 实现替换：

logexp <- function(days = 1)
     {
         linkfun <- function(mu) qlogis(mu^(1/days))
         linkinv <- function(eta) plogis(eta)^days
         mu.eta <- function(eta) days * plogis(eta)^(days-1) * binomial()$mu_eta
         valideta <- function(eta) TRUE
         link <- paste0("logexp(", days, ")")
         structure(list(linkfun = linkfun, linkinv = linkinv,
                        mu.eta = mu.eta, valideta = valideta, name = link),
                   class = "link-glm")
     }

您在上面指定的链接实际上确实有这样一个模型的示例（但它确实需要最新版本的lme4）。

于 2013-09-25T18:16:13.387 回答