search

rvlenth / emmeans

Estimated marginal means
https://rvlenth.github.io/emmeans/
341 stars 30 forks source link

Error in emmeans for averaging object (mth$objs[[1]])(object, trms, xlev, grid, ...): Unable to match model terms #402

Closed derek-corcoran-barrios closed 1 year ago

derek-corcoran-barrios commented 1 year ago

I already posted this in stackoverflow, but I decided to add it here, I went through some other issues and stackoverflow question dealing with a similar issue, but so far non of them has worked, here is the link of the stackoverflow post

emmeans can't evaluate averaging object

The problem is that emmeans gives an error when evaluating an averaging object form the MuMIn package, even when it says it should in this link I have been trying to debug this for a couple of days with no luck. All data and code is in this repo

first we load the needed r packages and the dataset

library(emmeans)
library(lme4)
#;-) Loading required package: Matrix
library(MuMIn)
library(doParallel)
#;-) Loading required package: foreach
#;-) Loading required package: iterators
#;-) Loading required package: parallel

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)

Working with the best model works

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

Working with the average model does not works

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)
#;-) Error in (mth$objs[[1]])(object, trms, xlev, grid, ...): Unable to match model terms
Standard output and standard error ``` sh -- nothing to show -- ```
Session info ``` r sessioninfo::session_info() #;-) ─ Session info ─────────────────────────────────────────────────────────────── #;-) setting value #;-) version R version 4.2.2 Patched (2022-11-10 r83330) #;-) os Ubuntu 20.04.5 LTS #;-) system x86_64, linux-gnu #;-) ui X11 #;-) language en_US:en #;-) collate en_US.UTF-8 #;-) ctype en_US.UTF-8 #;-) tz Europe/Copenhagen #;-) date 2023-02-21 #;-) pandoc 2.19.2 @ /usr/lib/rstudio/bin/quarto/bin/tools/ (via rmarkdown) #;-) #;-) ─ Packages ─────────────────────────────────────────────────────────────────── #;-) package * version date (UTC) lib source #;-) boot 1.3-28 2021-05-03 [4] CRAN (R 4.0.5) #;-) cli 3.6.0 2023-01-09 [1] CRAN (R 4.2.2) #;-) coda 0.19-4 2020-09-30 [3] CRAN (R 4.0.2) #;-) codetools 0.2-19 2023-02-01 [4] CRAN (R 4.2.2) #;-) digest 0.6.31 2022-12-11 [1] CRAN (R 4.2.2) #;-) doParallel * 1.0.17 2022-02-07 [1] CRAN (R 4.2.1) #;-) dplyr 1.1.0 2023-01-29 [1] CRAN (R 4.2.2) #;-) emmeans * 1.8.4-1 2023-01-17 [1] CRAN (R 4.2.2) #;-) estimability 1.4.1 2022-08-05 [1] CRAN (R 4.2.1) #;-) evaluate 0.20 2023-01-17 [1] CRAN (R 4.2.2) #;-) fansi 1.0.4 2023-01-22 [3] CRAN (R 4.2.2) #;-) fastmap 1.1.0 2021-01-25 [3] CRAN (R 4.0.3) #;-) foreach * 1.5.2 2022-02-02 [1] CRAN (R 4.2.1) #;-) fs 1.6.1 2023-02-06 [3] CRAN (R 4.2.2) #;-) generics 0.1.3 2022-07-05 [1] CRAN (R 4.2.1) #;-) glue 1.6.2 2022-02-24 [3] CRAN (R 4.1.2) #;-) htmltools 0.5.4 2022-12-07 [1] CRAN (R 4.2.2) #;-) iterators * 1.0.14 2022-02-05 [1] CRAN (R 4.2.1) #;-) knitr 1.42 2023-01-25 [1] CRAN (R 4.2.2) #;-) lattice 0.20-45 2021-09-22 [4] CRAN (R 4.2.0) #;-) lifecycle 1.0.3 2022-10-07 [3] CRAN (R 4.2.1) #;-) lme4 * 1.1-31 2022-11-01 [1] CRAN (R 4.2.2) #;-) magrittr 2.0.3 2022-03-30 [3] CRAN (R 4.1.3) #;-) MASS 7.3-58.2 2023-01-23 [4] CRAN (R 4.2.2) #;-) Matrix * 1.5-3 2022-11-11 [4] CRAN (R 4.2.2) #;-) minqa 1.2.5 2022-10-19 [3] CRAN (R 4.2.1) #;-) multcomp 1.4-20 2022-08-07 [1] CRAN (R 4.2.1) #;-) MuMIn * 1.47.1 2022-09-01 [1] CRAN (R 4.2.2) #;-) mvtnorm 1.1-3 2021-10-08 [1] CRAN (R 4.2.1) #;-) nlme 3.1-162 2023-01-31 [4] CRAN (R 4.2.2) #;-) nloptr 2.0.3 2022-05-26 [3] CRAN (R 4.2.0) #;-) pillar 1.8.1 2022-08-19 [1] CRAN (R 4.2.1) #;-) pkgconfig 2.0.3 2019-09-22 [3] CRAN (R 4.0.0) #;-) purrr 1.0.1 2023-01-10 [1] CRAN (R 4.2.2) #;-) R.cache 0.16.0 2022-07-21 [1] CRAN (R 4.2.1) #;-) R.methodsS3 1.8.2 2022-06-13 [1] CRAN (R 4.2.1) #;-) R.oo 1.25.0 2022-06-12 [1] CRAN (R 4.2.1) #;-) R.utils 2.12.2 2022-11-11 [1] CRAN (R 4.2.2) #;-) R6 2.5.1 2021-08-19 [3] CRAN (R 4.1.1) #;-) Rcpp 1.0.10 2023-01-22 [3] CRAN (R 4.2.2) #;-) reprex 2.0.2 2022-08-17 [1] CRAN (R 4.2.2) #;-) rlang 1.0.6 2022-09-24 [1] CRAN (R 4.2.1) #;-) rmarkdown 2.20 2023-01-19 [1] CRAN (R 4.2.2) #;-) rstudioapi 0.14 2022-08-22 [1] CRAN (R 4.2.1) #;-) sandwich 3.0-2 2022-06-15 [1] CRAN (R 4.2.1) #;-) sessioninfo 1.2.2 2021-12-06 [3] CRAN (R 4.1.2) #;-) styler 1.8.1 2022-11-07 [1] CRAN (R 4.2.2) #;-) survival 3.5-3 2023-02-12 [4] CRAN (R 4.2.2) #;-) TH.data 1.1-1 2022-04-26 [1] CRAN (R 4.2.1) #;-) tibble 3.1.8 2022-07-22 [1] CRAN (R 4.2.1) #;-) tidyselect 1.2.0 2022-10-10 [1] CRAN (R 4.2.1) #;-) utf8 1.2.3 2023-01-31 [3] CRAN (R 4.2.2) #;-) vctrs 0.5.2 2023-01-23 [1] CRAN (R 4.2.2) #;-) withr 2.5.0 2022-03-03 [3] CRAN (R 4.1.3) #;-) xfun 0.37 2023-01-31 [1] CRAN (R 4.2.2) #;-) xtable 1.8-4 2019-04-21 [1] CRAN (R 4.2.0) #;-) yaml 2.3.7 2023-01-23 [3] CRAN (R 4.2.2) #;-) zoo 1.8-11 2022-09-17 [1] CRAN (R 4.2.1) #;-) #;-) [1] /home/au687614/R/x86_64-pc-linux-gnu-library/4.2 #;-) [2] /usr/local/lib/R/site-library #;-) [3] /usr/lib/R/site-library #;-) [4] /usr/lib/R/library #;-) #;-) ────────────────────────────────────────────────────────────────────────────── ```
rvlenth commented 1 year ago

I can't look at anything because none of the RDS files are present in the github site. All that is there are two README files, an .Rproj file, a .png file, and a .gitignore file that I think lists some of the files we need.

rvlenth commented 1 year ago

One comment I would make though, while I am waiting, is that your call to emmeans() works only accidentally. The specs argument should just give the factors; it is not supposed to be a model statement. And the pairwise ~, just makes it do additional pairwise comparisons that you never use.

noise.emm <- emmeans(BestModel, ~ year * initial_habitat * treatment, at = list(year = 4), data = Data)
pairs(noise.emm, simple = "treatment")

Also, when you report issues, I really, really prefer plain and simple steps, not a series of pipes. I don't need all that stuff with |> as.data.frame |> ... |> kable(); those are outside of emmeans and have no relevance to whatever issues exist in emmeans, and including those just masks what emmeans() (or pairs() in this case) actually produced.

derek-corcoran-barrios commented 1 year ago

Thanks, that was very clumsy of me, I already updated the repo with the rds files, I will try what you said know and report back

rvlenth commented 1 year ago

Never mind -- got 'em using ReadRDS

derek-corcoran-barrios commented 1 year ago

,sorry about that, by the way thanks a lot for taking the time to debugg this, your package is really awesome

rvlenth commented 1 year ago

Here's what I found in debug mode in part of emmeans:::emm_basis.averaging:

Browse[2]> colnames(X)
 [1] "(Intercept)"                   "aspect"                        "elevation"                    
 [4] "I((year - 1)^2)"               "I(abs(year - 1))"              "initial_habitatMeadow"        
 [7] "initial_habitatRangeland"      "initial_habitatForest:year"    "initial_habitatMeadow:year"   
[10] "initial_habitatRangeland:year" "year:treatmentControl" 

Browse[2]> bhat
                     (Intercept)                  I((year - 1)^2)                 I(abs(year - 1)) 
                    2.155205e+00                    -1.824446e-02                     1.649255e-01 
           initial_habitatMeadow         initial_habitatRangeland treatmentPermanentExclosure:year 
                    6.357188e-01                     7.759404e-01                    -7.279999e-02 
           treatmentControl:year       initial_habitatForest:year       initial_habitatMeadow:year 
                   -1.702733e-02                    -2.934701e-02                    -3.886043e-02 
   initial_habitatRangeland:year            treatmentControl:year                        elevation 
                   -5.070798e-02                    -1.702733e-02                     7.891242e-04 
                          aspect 
                    2.044456e-05

The problem is we need to match the names of bhat with the column names of X; and in one we have treatmentControl:year and the other we have year:treatmentControl.

Going back further, we have

> names(coef(AV, full = TRUE))
 [1] "(Intercept)"                      "I((year - 1)^2)"                  "I(abs(year - 1))"                
 [4] "initial_habitatMeadow"            "initial_habitatRangeland"         "treatmentPermanentExclosure:year"
 [7] "treatmentControl:year"            "initial_habitatForest:year"       "initial_habitatMeadow:year"      
[10] "initial_habitatRangeland:year"    "treatmentControl:year"            "elevation"                       
[13] "aspect"  

> formula(AV)
richness ~ aspect + elevation + I((year - 1)^2) + I(abs(year - 
    1)) + initial_habitat + initial_habitat:year + year:treatment

Notice the last term is year:treatment which is inconsistent with the coefficients which imply it was treatment:year.

The other glaring problem is that coef(AV) includes treatmentControl:year twice -- elements 7 and 11. Clearly, something is amiss, and I wonder if the models that are being averaged are specified inconsistently. Indeed,

> ML <- attr(AV, "modelList")
> lapply(ML, formula)
$`275`
richness ~ I((year - 1)^2) + I(abs(year - 1)) + initial_habitat + 
    (1 | block_no) + treatment:year

$`403`
richness ~ I((year - 1)^2) + I(abs(year - 1)) + initial_habitat + 
    (1 | block_no) + initial_habitat:year + year:treatment

$`274`
richness ~ I(abs(year - 1)) + initial_habitat + (1 | block_no) + 
    treatment:year

$`402`
richness ~ I(abs(year - 1)) + initial_habitat + (1 | block_no) + 
    initial_habitat:year + year:treatment

$`283`
richness ~ I((year - 1)^2) + I(abs(year - 1)) + elevation + initial_habitat + 
    (1 | block_no) + treatment:year

$`279`
richness ~ I((year - 1)^2) + I(abs(year - 1)) + aspect + initial_habitat + 
    (1 | block_no) + treatment:year

$`411`
richness ~ I((year - 1)^2) + I(abs(year - 1)) + elevation + initial_habitat + 
    (1 | block_no) + initial_habitat:year + year:treatment

$`282`
richness ~ I(abs(year - 1)) + elevation + initial_habitat + (1 | 
    block_no) + treatment:year

and note that these are inconsistent on including terms treatment:year and year:treatment. I think if you can re-do this in a way that the models are consistently specified, this issue will go away. Moreover, I suspect that the inconsistent model specifications are screwing up the model-averaging algorithm.

derek-corcoran-barrios commented 1 year ago

Hey Russell Thanks I have been working on this for a couple of days, and I made a far simpler example with my top 2 best models made by hand:

library(MuMIn)
library(lme4)
library(emmeans)

Model1 <- glmer(formula = richness ~ I((year - 1)^2) + I(abs(year - 1)) +  initial_habitat + (1 | block_no) + treatment:year, data = Data, family = poisson, control = glmerControl(optimizer = "bobyqa",  optCtrl = list(maxfun = 1e+05)))

Model2 <- glmer(formula = richness ~ I((year - 1)^2) + I(abs(year - 1)) + initial_habitat + (1 | block_no) + initial_habitat:year + treatment:year, data = Data, family = poisson, control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+05)))

TestList <- list(Model1 = Model1, 
                 Model2 = Model2)

AV <- model.avg(TestList, fit = TRUE)

noise.emm <- emmeans(AV, ~ year * initial_habitat * treatment, at = list(year = 4), data = Data)

And I am getting exactly the same error.

Error in (mth$objs[[1]])(object, trms, xlev, grid, ...) : 
  Unable to match model terms

Now looking at the coefficients as you did I found that again the same parameter (interaction between treatment:year) is giving problems, as you can see here:

coef(AV)

results in 

                     (Intercept)                  I((year - 1)^2) 
                     2.154985923                     -0.029031801 
                I(abs(year - 1))            initial_habitatMeadow 
                     0.192109142                      0.619439556 
        initial_habitatRangeland treatmentPermanentExclosure:year 
                     0.805688309                     -0.112639199 
           treatmentControl:year       initial_habitatForest:year 
                    -0.005549537                     -0.071254762 
      initial_habitatMeadow:year    initial_habitatRangeland:year 
                    -0.096462928                     -0.127065079 
           treatmentControl:year 
                    -0.005549537

As you can see the 7th coeficient is treatmentControl:year and the 11th parameter is the same, this is so weird, I guess this is more a MuMIn issue than a emmeans issue, do you have any clue on what is going on here?

rvlenth commented 1 year ago

That is strange indeed. My initial thought is it's something that happened in pdredge -- sort of a "right hand doesn't know what the left hand is doing" situation. But you've demonstrated it is not that, and in fact that the same thing happens when the models both specify treatment:year in that order.

I do think it'd be a good idea to report it to the MuMIn site. Maybe you can construct a very simple, small example where the same thing happens.

I do see some unusual characteristics of this model specifications: there is no main effect of either year or treatment, yet we have the interaction of these factors. It seems like it'd solve the problem to manually create those interaction variables, e.g. yt1 <- year*(treatment=="Control") and include those variables instead of treatment:year. Then it'd be impossible to reverse the names.

bbolker commented 1 year ago

While researching this problem (before I realized you'd gone much farther than I had already) I made the following change in emm_basis.averaging to issue a marginally more informative error message: would you like to incorporate it, or to have a pull request?

    missing.terms = colnames(X)[is.na(perm)]
    if (length(missing.terms) > 0) {
        stop ("In emm_basis.averaging: Unable to match model terms ",
              paste(missing.terms, collapse = ", "),
              call. = FALSE)
    }
rvlenth commented 1 year ago

Ben,

Thanks so much. This is a good idea. Why don't I just paste this in, rather than make you create a PR?

rvlenth commented 1 year ago

This does look much better, and doesn't show us messages about my custom method-dispatching quagmire:

> noise.emm <- emmeans(AV, ~ year * initial_habitat * treatment, at = list(year = 4), data = Data)
Error: In emm_basis.averaging: Unable to match model term(s):
year:treatmentControl
rvlenth commented 1 year ago

OK, I've got something that works: For the recent example with just two models averaged:

> emmeans(AV, ~ year * initial_habitat * treatment, data = Data)
 year initial_habitat treatment          emmean     SE  df lower.CL upper.CL
 1.99 Forest          PermanentExclosure   2.25 0.1777 234     1.90     2.60
 1.99 Meadow          PermanentExclosure   2.85 0.2080 234     2.44     3.26
 1.99 Rangeland       PermanentExclosure   3.01 0.1412 234     2.73     3.29
 1.99 Forest          Control              2.24 0.1382 234     1.97     2.51
 1.99 Meadow          Control              2.84 0.1752 234     2.49     3.18
 1.99 Rangeland       Control              3.00 0.0835 234     2.83     3.16

Results are given on the log (not the response) scale. 
Confidence level used: 0.95

What I did was use strsplit(, ":") to separate the interacting factors with each set of names, and then matched those sets of names. I coded it to still end up with the same perm vector and the same error message if we don't get a match. You'll see a message here when I've pushed up the revised code.

rvlenth commented 1 year ago

I guess I referred to this as a "bug fix", though I don't really think it's my bug that I fixed!

bbolker commented 1 year ago

Something like your strsplit hack occurred to me, although I think the problem would be better diagnosed and handled upstream (might be even more effort, though, and definitely not your problem ...)

rvlenth commented 1 year ago

Just got a chance to look at this more closely. Recal in our example that coef's at indices 7 and 11 have the same name and identical values:

> coef(AV, full=TRUE)
                     (Intercept)                  I((year - 1)^2)                 I(abs(year - 1)) 
                     2.154968701                     -0.029030951                      0.192107028 
           initial_habitatMeadow         initial_habitatRangeland treatmentPermanentExclosure:year 
                     0.619482054                      0.805708818                     -0.061322102 
           treatmentControl:year       initial_habitatForest:year       initial_habitatMeadow:year 
                    -0.005549882                     -0.032463194                     -0.043947716 
   initial_habitatRangeland:year            treatmentControl:year 
                    -0.057890015                     -0.005549882

However, the covariances are not identical. Here are the 7th and 11th rows of the covariance matrix.

> vcov(AV, full=TRUE)[7,]
                     (Intercept)                  I((year - 1)^2)                 I(abs(year - 1)) 
                   -2.511656e-03                    -8.423893e-05                    -3.325564e-05 
           initial_habitatMeadow         initial_habitatRangeland treatmentPermanentExclosure:year 
                    1.416865e-03                     3.070248e-03                     3.419387e-03 
           treatmentControl:year       initial_habitatForest:year       initial_habitatMeadow:year 
                    3.604075e-03                    -2.056235e-03                    -2.781465e-03 
   initial_habitatRangeland:year            treatmentControl:year 
                   -3.657338e-03                     3.262311e-03 

> vcov(AV, full=TRUE)[11,]
                     (Intercept)                  I((year - 1)^2)                 I(abs(year - 1)) 
                   -2.345281e-03                     7.198211e-07                    -9.524718e-06 
           initial_habitatMeadow         initial_habitatRangeland treatmentPermanentExclosure:year 
                    1.436469e-03                     3.086873e-03                     3.204092e-03 
           treatmentControl:year       initial_habitatForest:year       initial_habitatMeadow:year 
                    3.262311e-03                    -2.026883e-03                    -2.743934e-03 
   initial_habitatRangeland:year            treatmentControl:year 
                   -3.614440e-03                     3.262311e-03

So for starters, the variances for treatmentControl:year do not match ( 3.604075e-03 vs. 3.262311e-03). In the 7th row, we have the covariance between the two versions of the coefficient unequal to the variance. And generally, the covariances between the two versions of treatmentControl:year and other coefficients ((Intercept), etc.) are not the same. Sometimes they are in the same ballpark, sometimes they differ by more than an order of magnitude (e.g., with I((year - 1)^2), though these are very small values).

So what's clear to me is that my "bug fix" yields a set of estimates, SEs, etc, but not the estimates.

If I run temmeans call in debug mode and force it to select the 11th coefficient instead of the 7th, I get

year initial_habitat treatment          emmean     SE  df lower.CL upper.CL
 1.99 Forest          PermanentExclosure   2.25 0.1777 234     1.90     2.60
 1.99 Meadow          PermanentExclosure   2.85 0.2080 234     2.44     3.26
 1.99 Rangeland       PermanentExclosure   3.01 0.1412 234     2.73     3.29
 1.99 Forest          Control              2.24 0.1380 234     1.97     2.51
 1.99 Meadow          Control              2.84 0.1755 234     2.49     3.18
 1.99 Rangeland       Control              3.00 0.0843 234     2.83     3.16

Results are given on the log (not the response) scale. 
Confidence level used: 0.95

The estimates are the same, because the coefs are; but the last 3 SEs are a little bit different.

rvlenth commented 1 year ago

@derek-corcoran-barrios -- have you reported this error to MuMIn?

derek-corcoran-barrios commented 1 year ago

Yes I have, but no response yet, I can't find a github for the report, so I sent him an email to the one in the package, I did this 5 days ago I think

bbolker commented 1 year ago

There is an r-forge site for MuMIn with a bug tracker (you may need to have an r-forge account/be logged in order to post to it ...) (I see very recent commits to the SVN version control system there, so development is actively happening there ...)

derek-corcoran-barrios commented 1 year ago

yt1 <- year*(treatment=="Control")

Dear Russell I have kept working on this, I wanted to try this approach, however I am a bit confused on how to deal with this.

I tried to fit a main model with:

Model <- glmer(richness ~ aspect + elevation + initial_habitat + I(abs(year - 1)) + I((year-1)^2) + slope + treatment:initial_habitat + year:initial_habitat + year*(treatment=="Control") + year*(treatment=="PermanentExclosure") + year:treatment:initial_habitat + (1 | block_no), family = poisson, data = Mols2, control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 100000)))

and that works, however I cannot dredge that one, I feel like a novice here.

derek-corcoran-barrios commented 1 year ago

Hey @bbolker and @rvlenth I made it work!!!

So there was a fix that was made to the MuMIn package, but also there was a better call that I made to the dredge model averaging funcion, so when I made it:

AV <- model.avg(SelectedList, fit = TRUE, random = TRUE)

It fixed the call to emmeans, so of course as we all suspected it was a MuMIn issue, but this might be useful to all of us, so after doing that everything works

rvlenth commented 1 year ago

Glad you got it to work.

I think we may have this issue pretty well completed, so will close it. That does not prevent comments; and it can even be re-opened if necessary.

rvlenth commented 1 year ago

Just a note to @bbolker. This issue focused on correctly matching coefficients to a model matrix. But #409 turns this all on its head, and there I overhauled the code so we construct a model-like matrix that matches the coefficient names, rather than vice versa. This is needed because some models may come out parameterized differently due to quirks in the ways that R excludes the first level of a factor or not.