尝试在没有来自 bam() 输出的随机影响的情况下进行预测时出错

问题描述

我有一个数据集，我想在 mgcv 包中与 bam() 配合使用。该模型有一个二元结果，我需要为每个动物 ID 指定随机截距。下面是数据的一个子集（我的实际数据更大，协变量更多）：

dat2 <- read.csv('https://github.com/silasbergen/example_data/raw/main/dat2.csv')
dat2$Animal_id <- factor(dat2$Animal_id)
> head(dat2)
  Animal_id DEM_IA Anyrisk
1       105 279.94       0
2       105 278.68       0
3       106 329.13       0
4       106 329.93       0
5       106 332.25       0
6       106 333.52       0
> summary(dat2)
 Animal_id        DEM_IA         Anyrisk      
 105:     2   Min.   :156.3   Min.   :0.0000  
 106: 83252   1st Qu.:246.8   1st Qu.:0.0000  
 107: 22657   Median :290.1   Median :0.0000  
 108:104873   Mean   :284.8   Mean   :0.3619  
 109:142897   3rd Qu.:318.0   3rd Qu.:1.0000  
 110: 53967   Max.   :411.8   Max.   :1.0000 

我想拟合模型并预测没有随机效应的新数据：

library(mgcv)
mod <- bam(Anyrisk ~s(Animal_id,bs="re") + s(DEM_IA),data = dat2,family = "binomial",discrete=TRUE)
topred <-  data.frame(DEM_IA = c(280,320))
predict(mod,newdata = topred,exclude="s(Animal_id)",newdata.guaranteed = TRUE)

但这会引发错误：

Error in eval(predvars,data,env) : object 'Animal_id' not found

当我特别告诉它从预测中排除这个词时，它为什么需要 Animal_id？这也特别奇怪，因为我可以运行 ?random.effects mgcv 帮助文件中的类似示例，没问题，即使我修改这些示例以使用 bam() 而不是 gam()！任何帮助将不胜感激！

编辑

我可能找到了解决办法；显然，如果在 discrete=TRUE 模型中使用 bam()，那么 predict.bam() 也使用 discrete=TRUE，这将不在缺少随机效应的情况下工作，但这有效：

mod<- bam(Anyrisk ~s(Animal_id,topred,newdata.guaranteed = TRUE,discrete=FALSE)

输出：

         1          2 
-0.4451066 -0.0285989 

解决方法

tl;dr 通过将 something 放入 Animal_id 来解决此问题，您指定的值（不是 NA虽然...）

为什么？如果不深入研究代码就不能确定，但​​是……使用 model.frame(formula,newdata) 作为计算所需模型矩阵的步骤通常很方便。 （例如，可以先构建整个模型矩阵，然后将要忽略的列归零……）弄清楚可以从公式中删除哪些项可能是一个单独的、更困难的步骤。 （我不知道为什么它在 bam 和 gam 中的工作方式不同......）

这似乎工作正常：

topred <-  data.frame(DEM_IA = c(280,320),Animal_id=dat2$Animal_id[1])
predict(mod,newdata = topred,exclude="s(Animal_id)",newdata.guaranteed = TRUE)

检查您为 Animal_id 指定的内容是否真的无关紧要：

res <- lapply(levels(dat2$Animal_id),function(i) {
             dd <- transform(topred,Animal_id=i)
               predict(mod,newdata = dd,newdata.guaranteed = TRUE)
           })
do.call(rbind,res)

结果：

              1          2
[1,] -0.4451066 -0.0285989
[2,] -0.4451066 -0.0285989
[3,] -0.4451066 -0.0285989
[4,] -0.4451066 -0.0285989
[5,] -0.4451066 -0.0285989
[6,] -0.4451066 -0.0285989