尝试在没有 bam() 输出随机效应的情况下进行预测时出错

Question

我有一个数据集，我正在尝试将其与 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

，这将 not 与丢失随机效应一起工作，但这有效：

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,topred, exclude="s(Animal_id)",newdata.guaranteed = TRUE,discrete=FALSE)

输出：

         1          2 
-0.4451066 -0.0285989 

Answer 1

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, 
                       exclude="s(Animal_id)",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

Answer 2

当我在 newdata 中包含一些 NA 时，即使选择 na.action 过程也会出现错误：

顶红<- data.frame(DEM_IA = c(280,320, NA)) predict(mod, topred, exclude="s(Animal_id)", newdata.guaranteed = TRUE, discrete=FALSE, na.action = na.pass)

Predict.matrix.tprs.smooth(object, dk$data) 中的错误： 外部函数调用中的 NA/NaN/Inf（参数 1）

gam() 和 bam() 都会发生这种情况。