Ошибка при попытке предсказать без случайного эффекта от вывода bam()

Question

У меня есть набор данных, который я пытаюсь подогнать с помощью bam() в пакете mgcv. Модель имеет бинарный результат, и мне нужно указать случайные перехваты для каждого идентификатора животного. Подмножество данных приведено ниже (мои фактические данные намного больше с большим количеством ковариат):

      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()! Любая помощь будет принята с благодарностью!

РЕДАКТИРОВАТЬ

Возможно, я нашел исправление; по-видимому, если использовать в bam()модель, затем predict.bam()также использует discrete=TRUEкоторый не будет работать с отсутствующим случайным эффектом, но это работает:

      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 обойти это, вставив что- то для , не имеет значения, какое значение вы укажете (не 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