0

Thanks in advance for all help.

I would like to produce a plot similar to the below from my model averaged results. This plot was produced from a single model using the effect() function in the effects package.

Desired plot outcome

To my knowledge this cannot be achieve from model averaged results from the model.avg() function so I have tried to achieve a similar result by first predicting my model averaged results and then creating a plot.

I have built two models from this data:

igm_20 <- glmmTMB(igm_pres ~  fRHDV2_arrive_cat + fseason + sage + save_ajust_abun + fseason*fRHDV2_arrive_cat + (1 | fsite), data = test, family = binomial)
igm_21 <- glmmTMB(igm_pres ~  fRHDV2_arrive_cat + fseason + sage + save_ajust_abun + fseason*fRHDV2_arrive_cat + sage*fRHDV2_arrive_cat + (1 | fsite), data = test, family = binomial)

I can average these models like so:

mod_ave_list_1 <- list(igm_20, igm_21)
mod_ave_1 <- model.avg(mod_ave_list_1, rank = AICc)
s1 <- summary(mod_ave_1)

I have then tried to make predictions from the model averaged results for each combination of fseason*fRHDV2_arrive_cat from the dataset provided, test:

a <- as.data.frame(c("Summer", "Autumn", "Winter", "Spring", "Summer", "Autumn", "Winter", "Spring"))
a$fRHDV2_arrive_cat <- c("Pre-arrival", "Post-arrival", "Pre-arrival", "Post-arrival", "Pre-arrival", "Post-arrival", "Pre-arrival", "Post-arrival")
colnames(a) <- c("fseason", "fRHDV2_arrive_cat")

predict(mod_ave_1, full = TRUE, backtransform = TRUE, newdata = a, se.fit = TRUE)

The predict function runs fine if I do not include newdata = a, consequently I believe the data frame I supply to newdata is not of the appropriate structure.

If someone is able to help me to structure newdata properly so that I can produce predictions for each combination of fseason*fRHDV2_arrive_cat it would be much appreciated? I beleive I can produce the desired plot once I have the predictions.

My situation is very similar to that described in another post (here), which remains unanswered. Above I have described my attempt to achieve a similar plot via another means. I also note that there are other similar posts, such as here; I have not found these useful to my situation.

Pat Taggart
  • 321
  • 1
  • 9

1 Answers1

0

Your newdata must include all the explanatory variables used by the model, which it currently does not.

Variables required for prediction:

all.vars(formula(mod_ave_1)[-2]) # extracts all names from formula minus response
# [1] "fRHDV2_arrive_cat" "fseason"           "sage"              "save_ajust_abun"

Variables in the new data:

colnames(a)
# [1] "fseason"           "fRHDV2_arrive_cat"
Kamil Bartoń
  • 1,482
  • 9
  • 10
  • Thanks Kamil. I can now get the ```predict()``` function to work. However, I am interested in how one decides what value to include in ```newdata``` for random effects in the model, in my instance ```fsite```? Ultimately I am interested in making predictions for each combination of ```fseason*fRHDV2_arrive_cat```. I can essentially remove other numeric variables in my model from model averaged predictions by assigning a value of ```0``` in ```newdata```, but I cannot do the same for ```fsite``` as it is a factor. What would one do in this situation, assign the reference level to ```fsite```? – Pat Taggart Jan 02 '20 at 06:20
  • @PatTaggart you can predict at the population level, i.e. setting the random effects to zero. This is done via the `re.form` argument to `predict` (although `?predict.glmmTMB` says it is "not yet implemented", but it does not clarify if the predictions include RE or not). It is implemented in `lmer` though, and for older models such as `lme` the syntax was `predict(..., level = 0)`). – Kamil Bartoń Jan 02 '20 at 19:57
  • Thank you very much – Pat Taggart Jan 02 '20 at 21:31