1

我有一个数据集，我正试图与 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=TRUEwhich 将无法使用缺少随机效应，但这有效：

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
4

1 回答 1

3

tl;博士通过在 for中添加一些东西来解决这个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
于 2020-12-22T22:15:50.217 回答