I've used MuMIn::dredge() on a full model glmmTMB::glmmTMB() and averaged the top models. I would like to get predictions for a raster stack, but I keep getting an error. See example below using terra::predict() and a wrapper for MuMIn:::predict.averaging().
I think there are multiple issues...as it is not recognizing the arguments full or re.from as it seems to still be making predictions depending on random effect level.
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
Here is a somewhat simplified example that uses the same data for predict.averaging as for terra::predict. The results are the same. Perhaps you can use this to help focus your question a bit more.
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