机器算法验证 - 当人口普查间隔变化时对二元结果建模 - 吾爱随笔录

当人口普查间隔变化时对二元结果建模

机器算法验证 r 物流 lme4-nlme 死亡

2022-04-12 13:43:24

对于当前的一项工作，我正在尝试模拟英国林地中山毛榉树的树木死亡概率。我记录了树木在 3 个不同的人口普查期间是活还是死，以及它们的直径和生长速度的数据。每棵树都有一个 ID 号，因此可以在每个时间间隔识别它。然而，人口普查间隔各不相同，因此一次调查和另一次调查之间的时间是 4 年、12 年或 18 年。显然，普查期越长，一棵树在下一次调查时死亡的可能性就越大。我在制作一个现实的可重现示例时遇到了问题，因此您可以在此处找到数据。

数据集中的变量是：

ID - 树的唯一 ID
Block - 树所在的 20x20m 地块的 ID
死亡 - 树的状态，死亡 (1) 或活着 (0)
GR - 上次调查的年增长率
DBH - 树的胸高直径
SL - 人口普查之间的时间间隔（以年为单位）

一旦一棵树被记录为死亡，它就会从随后的调查中消失。

理想情况下，我希望能够使用有关直径和增长率的信息来估计一棵树的年死亡率。经过一段时间的搜索，我发现逻辑暴露模型似乎能够通过使用二项式模型的 logit 链接的更改版本来解释人口普查期间的差异，正如 Ben Bolker在此处详述的那样。这最初是由 Shaffer 用来确定巢的年龄（因此暴露）不同的鸟巢存活的每日概率。我没有看到它在巢生存模型的背景之外使用，但似乎我应该能够使用它来模拟暴露不同的生存/死亡率。

require(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
  logit_mu_eta <- function(eta) {
    ifelse(abs(eta)>30,.Machine$double.eps,
           exp(eta)/(1+exp(eta))^2)
    ## OR .Call(stats:::C_logit_mu_eta, eta, PACKAGE = "stats")
  }
  mu.eta <- function(eta) {       
    exposure * plogis(eta)^(exposure-1) *
      logit_mu_eta(eta)
  }
  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")
}

目前我的模型看起来像这样，但我会在进行过程中加入更多变量：

require(lme4)
Dead<-read.csv("Stack_dead.csv",)


M1<-glmer(Dead~DBH+(1|ID),data=Dead,family=binomial(logexp(Dead$SL))) 
#I use (1|ID) here to account for the repeated measurements of the same individuals
    summary(M1)

plot(Dead$DBH,plogis(predict(M1,re.form=NA)))

我主要想知道：

我用来控制人口普查之间时间差异的统计技术是否合理？如果不是，你能想出更好的方法来解决这个问题吗？
如果第一个问题的答案是肯定的，那么使用逆 logit (plogis) 是获得表示为概率的预测的正确方法吗？

提前感谢您的帮助！

1个回答

另一种稍微简单的方法是使用条件对数对数链接 ( cloglog)，它估计对数风险而不是死亡率等结果的对数几率。从rpubs复制：

生态学（和其他地方）的一个非常常见的情况是生存/二元结果模型，其中个体（每个测量一次）的暴露量不同。解决这个问题的经典方法是使用互补的日志-日志链接。互补的 log-log 或“cloglog”功能是 $C(\mu)=\log(−\log(1−\mu))$ ; 它的倒数是 $\mu=C^{−1}(\eta)=1−\exp(−\exp(\eta))$ . 因此，如果我们预计死亡率 $\mu_0$ 在一段时间内 $\Delta t=1$ 和线性预测器 $\eta=C^{−1}(\mu_0)$ 然后

C^{- 1} (η + \log Δ t) = (1 - \exp (- \exp (η) \cdot Δ t))

$C^{−1}(\eta+\log\Delta t)=(1−\exp(−\exp(\eta) \cdot \Delta t))$ 一些代数表明这等于

1 - (1 - μ_{0})^{Δ t}

$1−(1−\mu_0)^{\Delta t}$ ，这就是我们想要的。

功能 $\exp(−\exp(x))$ 被称为Gompertz函数（它也是极值分布的 CDF ），因此用这个反向链接函数拟合模型（即拟合cloglog与生存率，而不是死亡率，概率的链接）也称为gompit (或极值）回归。

要在 R 中使用这种方法，请在公式中指定family=binomial(link="cloglog")并添加一个形式的术语offset(log(exposure))（或者，一些建模函数将偏移量作为单独的参数）。例如，

glm(surv~x1+x2+offset(log(exposure)),
    family=binomial(link="cloglog"),
    data=my_surv_data)

其中exposure是特定个体面临死亡/失败可能性的时间长度（例如，普查间隔或观察之间的时间或总观察时间）。

您可能还需要考虑检查log(exposure)作为协变量而不是偏移量包含在内的模型 - 这使得对数风险具有 $\beta_t \log(t)$ 项，或等效地使危险与 $t^{\beta_t}$ 而不是 $t$ （我相信这使得生存分布 Weibull 而不是指数分布，但我没有仔细检查这个结论）。

使用这种方法而不是 Schaffer 的幂逻辑方法的优点：

因为曝光时间包含在偏移中，而不是链接函数本身的定义中，所以 R 会更优雅地处理这个问题（例如，从原始数据集中生成具有不同曝光时间的预测会更容易）。
它稍微老一点，在统计中使用更广泛；谷歌搜索“cloglog 逻辑回归”或在CrossValidated 上搜索 cloglog将为您带来更多资源。

我能想到的唯一缺点是巢穴生存世界中的人们更习惯于谢弗的方法。对于足够大的数据集，您可能能够分辨出哪个链接实际上更适合数据（例如，适合两种方法并比较 AIC 值），但总的来说，我怀疑差异很大。

其它你可能感兴趣的问题

上一篇来自指数分布的R中的泊松过程下一篇IIA 假设：差异 logit 和 probit