Average Partial Effects from model.avg() GAMMs

Question

I used the MuMIn::model.avg() function to generate a set of top models which represented 95% of the weight of all models considered. Using this model average object, I also generated predictions from my top models.

Say my top 3 models are as follows:

catch ~ s(effort) + s(year) + s(month) + s(open_water) + s(pop_dens)
catch ~ s(effort) + s(year) + s(month) + s(forest)
catch ~ s(effort) + s(year) + s(month) + s(open_water)

where open water/forest are proportions.

I would like to see what the average partial effects of open water is in a plot, but it shows up in two of my top models and each has a different weight associated with it. Is there a straightforward way to assess this?

The current approach I am taking is generating estimates from each model, multiplying by the weight of that model, then adding the estimates from the models to get an averaged estimate across the set of models.

ex code:

# get models with variable of interest
water_mods <- top_mod_list[grep(pattern = "water", x = names(top_mod_list))]
# get normalized weights
water_weights <-  model.sel(water_mods) %>% as.data.frame() %>%
  mutate(mod = rownames(.))

# use mean values of all other predictors, only change variable of interest
pdat <- means %>% select(-water)
pdat$water <- seq(min(landing_dat$water),
                  max(landing_dat$water),
                  length.out = 200)

# predict response based on changes to variable of interest
water_predict <- 
  map(water_mods,
      ~pdat %>% 
        mutate(pred = predict(.x,
                              newdata = pdat, 
                              # type = "response",
                              newdata.guaranteed = TRUE,
               weighted_pred = pred * water_weights$weight))) 

Answer 1

Why do you want to see the average function? The estimated functions mean two different things in the models this smooth occurs. In neither model is the term s(open_water)

the effect of open water on the response.

It is the

the effect of open water on the response given (conditional upon) the other terms in the model

And as you have different terms in the models the term s(open_water) means different things which are likely incompatible; it is very likely you would be averaging an apple and an orange, not two apples.

Rather than this AIC, multimodel averaging approach, I would approach the model estimation problem here by adding extra shrinkage penalties to all or a subset of terms. The former is done by adding select = TRUE to your call:

m <- gam(catch ~ s(effort) + s(year) + s(month) + s(open_water) + s(pop_dens) + s(forest),
  data = foo,
  method = "REML",
  select = TRUE,
  family = something)

Which will never result in the effect of one or more of the smooths being exactly 0 but they could get shrunk towards 0 and may have effectively 0 (but not exactly 0) effect on the response.