生态学中的常见情况是具有二元结果(0=死亡,1=生存)的生存模型,其中个体(在本例中考虑鸟类的个体筑巢尝试)在暴露于死亡风险的天数方面存在差异。为了解决这个问题,我们使用了修改后的逻辑回归,它将暴露天数合并到链接函数中。
正如 Shaffer (2004) 所描述的:
“通过选择适当的预测函数,以 x 为单位对每日生存率进行建模,在我们的例子中,该函数应该产生 0 到 1 之间的值。正如逻辑回归中所做的那样,我们使用 S 形逻辑函数:
我们的广义线性模型的系统分量是 [s(x)]t。接下来,我们考虑这个函数:
上述函数对于 θ 是单调且可微的,可以证明 g(θ) = β0 + β1x,满足广义线性模型中链接函数的标准。这三个组成部分:二项式响应分布、表达式 1 中给出的预测函数以及表达式 2 中给出的链接函数完全指定了我们的广义线性模型。该模型(以下简称“逻辑暴露模型”)与逻辑回归模型类似,但链接函数的形式有所不同。 Logistic 暴露链接函数的分子和分母中包含指数 (1/t),而 Logistic 回归链接函数中不存在该指数。指数对于解释区间幸存概率取决于区间长度这一事实是必要的。”
此链接函数的代码可在网络上找到,如果您输入“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)'未知链接: 'logexp(1)' 未知
另外:警告消息:
在 if (!(lTyp 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")
}