是否可以将 MuMIn:::predict.averaging() 与 glmmTMB 模型和栅格堆栈一起用作新数据？

Question

我在完整模型

MuMIn::dredge()

上使用了

glmmTMB::glmmTMB()

并对顶级模型进行了平均。我想获得栅格堆栈的预测，但我不断收到错误。请参阅下面使用

terra::predict()

和

MuMIn:::predict.averaging()

包装器的示例。

我认为存在多个问题......因为它没有识别参数

full

或

re.from

，因为它似乎仍在根据随机效应水平进行预测。

library(glmmTMB)
library(MuMIn)
library(tidyverse)
library(terra)
data(iris)
glimpse(iris)

#create more variables, center/scale
data<-iris%>%mutate(Species=as.factor(Species),
                    
                    Sepal.Length.2=Sepal.Length^2,
                    Sepal.Width.2=Sepal.Width^2,
                    Petal.Width.2=Petal.Width^2,
                    
                    Sepal.Length.Log=log(Sepal.Length),
                    Sepal.Width.Log=log(Sepal.Width),
                    Petal.Width.Log=log(Petal.Width))%>%
  mutate(across(c(-Species,-Petal.Length),~c(scale(.x))))

#full glmm
mod<-glmmTMB(Petal.Length~Sepal.Length+Sepal.Length.2+Sepal.Length.Log+
               Sepal.Width+Sepal.Width.2+Sepal.Width.Log+
               Petal.Width+Petal.Width.2+Petal.Width.Log+
               (1|Species),data,na.action='na.fail')
#dredge
mmodels<-dredge(mod, m.lim=c(2,3), fixed = ~(1|Species))
#get top models
confset.95p <- get.models(mmodels, subset = cumsum(weight) <= .95)
#average models
avgm <- model.avg(confset.95p)
colnames(avgm$x)
avgm$sw

#test predict on data frame
test<-avgm$x%>%as.data.frame()%>%select(-1)%>%mutate(Species=NA)
predict(avgm,test,type='link',se.fit=TRUE, re.form=~0,backtransform=FALSE, full=FALSE)#why is full ignored?

#generate raster data
sw <- sw2 <- swl <-sl <- sl2 <- sll <- pw <- rast(ncol=10, nrow=10)
set.seed(4)
values(sw) <- runif(ncell(sw), min=data$Sepal.Width, max=data$Sepal.Width)
values(sw2) <- runif(ncell(sw2), min=data$Sepal.Width.2, max=data$Sepal.Width.2)
values(swl) <- runif(ncell(swl), min=data$Sepal.Width.Log, max=data$Sepal.Width.Log)

values(sl) <- runif(ncell(sl), min=data$Sepal.Length,max=data$Sepal.Length)
values(sl2) <- runif(ncell(sl2), min=data$Sepal.Length.2,max=data$Sepal.Length.2)
values(sll) <- runif(ncell(sll), min=data$Sepal.Length.Log,max=data$Sepal.Length.Log)

values(pw) <- runif(ncell(pw), min=data$Petal.Width, max=data$Petal.Width)

spec<-(pw*0)+1

x <- c(sw, sw2, swl, sl,sl2,sll, pw, spec)
x[1] <- NA
names(x) <- c("Sepal.Width","Sepal.Width.2","Sepal.Width.Log",
              "Sepal.Length","Sepal.Length.2","Sepal.Length.Log",
              "Petal.Width","Species")

#predict on raster data
out<-terra::predict(x, avgm,type='link',se.fit=TRUE, re.form=~0,backtransform=FALSE, full=FALSE, fun = function(model, ...) MuMIn:::predict.averaging(model, ...)$se.fit)
#Error in MuMIn:::predict.averaging(model, ...) : 
#'predict' for models '26', '50', '42', '146', '82' and '274' caused errors
#In addition: There were 13 warnings (use warnings() to see them)
#Warning messages:
#1: In MuMIn:::predict.averaging(model, ...) : argument 'full' ignored
#2: In `[<-.factor`(`*tmp*`, ri, value = c(1, 1, 1, 1, 1,  ... : invalid factor level, NA generated
#3: In `[<-.factor`(`*tmp*`, ri, value = c(1, 1, 1, 1, 1,  ... : invalid factor level, NA generated
#4: In `[<-.factor`(`*tmp*`, ri, value = c(1, 1, 1, 1, 1,  ... : invalid factor level, NA generated
#5: In `[<-.factor`(`*tmp*`, ri, value = c(1, 1, 1, 1, 1,  ... : invalid factor level, NA generated
#6: In `[<-.factor`(`*tmp*`, ri, value = c(1, 1, 1, 1, 1,  ... : invalid factor level, NA generated
#7: In `[<-.factor`(`*tmp*`, ri, value = c(1, 1, 1, 1, 1,  ... : invalid factor level, NA generated
#8: In `contrasts<-`(`*tmp*`, value = aug_contrasts(c1, new_levels)) : wrong number of contrast matrix rows
#9: In `contrasts<-`(`*tmp*`, value = aug_contrasts(c1, new_levels)) : wrong number of contrast matrix rows
#10:In `contrasts<-`(`*tmp*`, value = aug_contrasts(c1, new_levels)) : wrong number of contrast matrix rows
#11:In `contrasts<-`(`*tmp*`, value = aug_contrasts(c1, new_levels)) : wrong number of contrast matrix rows
#12:In `contrasts<-`(`*tmp*`, value = aug_contrasts(c1, new_levels)) : wrong number of contrast matrix rows
#13:In `contrasts<-`(`*tmp*`, value = aug_contrasts(c1, new_levels)) : wrong number of contrast matrix rows

Answer 1

这是一个稍微简化的示例，其中

predict.averaging

与

terra::predict

使用相同的数据。结果是一样的。也许你可以用它来帮助你更加集中你的问题。

library(glmmTMB)
library(MuMIn)
library(terra)

data <- iris |> dplyr::mutate(Species=as.factor(Species),
                Sepal.Length.2=Sepal.Length^2,
                Sepal.Width.2=Sepal.Width^2,
                Petal.Width.2=Petal.Width^2,
                Sepal.Length.Log=log(Sepal.Length),
                Sepal.Width.Log=log(Sepal.Width),
                Petal.Width.Log=log(Petal.Width)) |>
        dplyr::mutate(dplyr::across(c(-Species,-Petal.Length), ~c(scale(.x))))

mod <- glmmTMB::glmmTMB(Petal.Length~Sepal.Length+Sepal.Length.2 + Sepal.Length.Log +
               Sepal.Width+Sepal.Width.2+Sepal.Width.Log +
               Petal.Width+Petal.Width.2+Petal.Width.Log +
               (1|Species), data, na.action='na.fail')

mmodels <- MuMIn::dredge(mod, m.lim=c(2,3), fixed = ~(1|Species))
confset.95p <- MuMIn::get.models(mmodels, subset = cumsum(weight) <= .95)
#Fixed terms are "cond((Int))" and "disp((Int))"
avgm <- MuMIn::model.avg(confset.95p)

test <- data.frame(avgm$x, Species=NA)[,-1]
p <- predict(avgm, test, type='link', se.fit=TRUE, re.form=~0, backtransform=FALSE)

#generate raster data
r <- terra::rast(nrow=10, ncol=15, nlyr=8, names=names(test), vals=test)
pr <- terra::predict(r, avgm, type='link', se.fit=TRUE, re.form=~0, backtransform=FALSE)
all(do.call(cbind, p) == values(pr))
# TRUE