The goal of emmeansProblem is to generate a reproducible example to solve the issues within our model averaging within the emmeans package
first we load the needed r packages and the dataset
library(emmeans)
library(lme4)
library(MuMIn)
library(doParallel)
Data <- readRDS("Data.rds")
Now we fit a general model:
Model <- glmer(richness ~ aspect + elevation +
initial_habitat +
I(abs(year - 1)) +
I((year-1)^2) +
slope +
treatment:initial_habitat +
year:initial_habitat +
year:treatment +
year:treatment:initial_habitat +
(1 | block_no), family = poisson, data = Data, control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 100000)))
And we make a model selection (skip this step since it takes a while, the “SelectRichness.rds” file is in this github)
options(na.action = "na.fail")
library(doParallel)
cl <- makeCluster(4)
registerDoParallel(cl)
clusterEvalQ(cl, library(lme4))
#> [[1]]
#> [1] "lme4" "Matrix" "stats" "graphics" "grDevices" "utils"
#> [7] "datasets" "methods" "base"
#>
#> [[2]]
#> [1] "lme4" "Matrix" "stats" "graphics" "grDevices" "utils"
#> [7] "datasets" "methods" "base"
#>
#> [[3]]
#> [1] "lme4" "Matrix" "stats" "graphics" "grDevices" "utils"
#> [7] "datasets" "methods" "base"
#>
#> [[4]]
#> [1] "lme4" "Matrix" "stats" "graphics" "grDevices" "utils"
#> [7] "datasets" "methods" "base"
clusterExport(cl, "Data")
Select <- MuMIn::pdredge(Model, extra = list(R2m = function(x) r.squaredGLMM(x)[1, 1], R2c = function(x) r.squaredGLMM(x)[1, 2]),fixed = ~YEAR:Treatment, cluster = cl)
stopCluster(cl)
saveRDS(Select, "SelectRichness.rds")
And now we Select the best models, I will do this twice, since the
outcome of the subset
function will be used to show how the best model
does not have issues and the averaged model from get.models
which is
the result I need is not working
Select <- readRDS("SelectRichness.rds")
Selected <- subset(Select, delta < 2)
SelectedList <- get.models(Select, delta < 2)
As specified above the goal is to find if the treatments do yieald differences by year 4, based on the model. So first we will show this with the best model
BestModel <- get.models(Selected, 1)[[1]]
noise.emm <- emmeans(BestModel, pairwise ~ year + initial_habitat + initial_habitat:year + year:treatment, at = list(year = 4), data = Data)
pairs(noise.emm, simple = "treatment") |>
as.data.frame() |> dplyr::filter(p.value < 0.05) |>
dplyr::arrange(initial_habitat, estimate) |>
dplyr::select(-SE, -df, -z.ratio) |>
knitr::kable()
contrast | year | initial_habitat | estimate | p.value |
---|---|---|---|---|
PermanentExclosure - Control | 4 | Forest | -0.2193592 | 3.06e-05 |
PermanentExclosure - Control | 4 | Meadow | -0.2193592 | 3.06e-05 |
PermanentExclosure - Control | 4 | Rangeland | -0.2193592 | 3.06e-05 |
This does not work
AV <- model.avg(SelectedList, fit = TRUE)
noise.emm_av <- emmeans(AV, pairwise ~ year + initial_habitat + initial_habitat:year + year:treatment, at = list(year = 4), data = Data)
pairs(noise.emm_av, simple = "treatment") |>
as.data.frame() |> dplyr::filter(p.value < 0.05) |>
dplyr::arrange(initial_habitat, estimate) |>
dplyr::select(-SE, -df, -z.ratio) |>
knitr::kable()
Giving the following error
Error in (mth$objs[[1]])(object, trms, xlev, grid, ...) :
Unable to match model terms