在 lme4 中使用 cloglog 链接进行 PIRLS 步进减半

这个问题从根本上说是由 cloglog 反向链接函数的上尾细引起的，该函数与参数的默认起始值交互;默认情况下，固定效应系数设置为零，随机效应（缩放）方差默认设置为 1。这通常是有道理的，但这就是为什么当我们有一个大的偏移量和一个 cloglog 链接时它不起作用的原因：

cc <- make.link("cloglog")$linkinv
par(las=1, bty = "l")
curve(1-cc(x), from = -5, to = 10, log = "y", n= 501,
          ylab = "1-Prob(x)", xlab = "Linear predictor")

因此，线性预测变量的任何值大于约 3.6 时，将给出一个常数值。1-.Machine$double.eps

如果所有偏移量都明显高于 3.6（即总丰度大于约 = 36.6），则所有预测值都将位于上层“架子”上;此外，逆链接函数的梯度在此区域中也为零，尽管 R 将其设置为可能下溢的区域。exp(3.6).Machine$double.eps

目前尚不清楚为什么这会破坏PIRLS算法;它确实取决于导数（请参阅此处的等式 13，但此处描述的场景并没有明显打破这个等式。（权重矩阵中使用的逆方差会很大，但同样，这不一定是问题。1/.Machine$double.eps

library(lme4)
simfun <- function(seed = NULL) {
    if (!is.null(seed)) set.seed(seed)
    dd <- data.frame(sp = factor(rep(1:30, each = 5)),
                     tot_ab = rnbinom(150, mu = 250, size = 50))
    dd$pres <- suppressMessages(simulate(~ (1|sp) + offset(log(tot_ab))
                      , family = binomial(link = "cloglog")
                      , newdata = dd
                      , newparam = list(beta=-8, theta = 1))[[1]])
    dd
}

解决方案：使截距（本例中唯一的固定效应参数）的起始值小得多，使预测的起始值在更合理的范围内：

glmer(pres ~ (1|sp) + offset(log(tot_ab)),
      family = binomial(link = "cloglog"),
      data = dd,
      start = list(fixef = -log(min(dd$tot_ab))))

拟合有效（尽管我们确实收到一些收敛警告）。

定义自己的链接函数以不同方式设置下限和上限可能会起作用，但这可能会遇到其他数值问题......cloglog

另一种解决方案：与其使用总丰度的对数作为偏移量，不如使用 ;也就是说，不是对每个采样的个体生物体的发生概率进行建模，而是对最小样本大小的样本中发生的概率进行建模。（您还应该能够选择与最小样本数量相同数量级的整数。下面的代码演示了在对模拟数据进行的 1000 次运行中，拟合总是失败，但总是成功（尽管可能存在奇异拟合、收敛警告等），并且具有最小缩放偏移量。log(tot_ab)/min(tot_ab)offset(log(tot_ab))

fitfun <- function(dd, fix_offset = FALSE) {
    if (!fix_offset) {
        form <- pres ~ (1|sp) + offset(log(tot_ab))
    } else {
        form <- pres ~ (1|sp) + offset(log(tot_ab)/min(tot_ab))
    }
    res <- try(
        suppressMessages(
            glmer(
                form
              , family = binomial(link = "cloglog")
              , data = dd)), silent = TRUE
    )
    if (!inherits(res, "try-error")) return("OK")
    as.character(res)
}

适合原始偏移量：

set.seed(101)
res <- replicate(1000, fitfun(simfun()))
table(res)

Error : (maxstephalfit) PIRLS step-halvings failed to reduce deviance in pwrssUpdate
                                                                                     8 
                                            Error : cannot generate feasible simplex
                                                                                   983 
                          Error : pwrssUpdate did not converge in (maxit) iterations

                                                                                     9 

现在拟合最小缩放偏移量：

res <- replicate(1000, fitfun(simfun(), fix_offset = TRUE))
table(res)

res
  OK 
1000