2

I am running a glmm using glmmTMB and using predict() to calculate predicted means. I want to know if it is possible to calculate predicted means based on specific fixed effects from the model.

run the model

model_1 <- glmmTMB(step.rate ~ Treatment*Week + Logger.ID + (1|Animal.ID), 
           data = data.df, family = nbinom1)

I am currently using the following code to calculate predicted means for all fixed effects:

new data

new.dat <- data.frame(Treatment = data.df$Treatment, Week = data.df$Week, Logger.ID = 
           data.df$Logger.ID, Animal.ID = data.df$Animal.ID)

predicted means

new.dat$predicted <- predict(model_1, new.data = new.dat, type = "response", re.form = NA)

I don't want Logger.ID to be included when calculating predicted means and want to know if it possible to ignore Logger.ID and how to do this.

alex
  • 153
  • 1
  • 10
  • Hi alex. Do you mean you want to know what the mean would be _without_ logger ID? Why not just re-run the model without logger ID? Your model already "ignores" other variables, such as animal age, by simply not including them. If you want to ignore a variable, just don't put it in your model. Or am I misunderstanding what you're asking? – Allan Cameron Aug 17 '20 at 13:54
  • Hi, yes I would like to calculate the mean without logger ID. Logger ID has to stay in the model as a fixed effect. I want to know how to calculate the predicted means from specific fixed effects in the model while ignoring other fixed effects in the model. I know how to tell predict to ignore all random effects, but not how to know ignore some of the fixed effects. – alex Aug 17 '20 at 14:09
  • But alex, ignoring a fixed effect is the same as leaving it out of the model. What you are saying is "I want to leave logger ID out of my model, but still have it included in my model". The predictions for the other fixed effects with logger ID taken out of your model is logically identical to "ignoring" logger ID. – Allan Cameron Aug 17 '20 at 14:14
  • I did wonder if that would be the case, but wanted to know for sure if it was possible and how to do it if so... Thank you for clearing that up! – alex Aug 17 '20 at 14:32
  • Is it possible to make a prediction for treatment and week only from this model? I had to include logger ID as a fixed effect in the model originally because one logger was a lot more sensitive to step rate than the other loggers used during the trial. I wanted to explain the variance caused by this logger rather than control for it if I was to fit logger ID as a random effect. – alex Aug 17 '20 at 14:39
  • If it's very skewed compared to the other loggers than I guess you could just choose the median logger as the basis for predictions. This does seem like you are effectively treating the logger as a random effect (which you would imagine is how it should be treated). You could try predicting with the median logger, removing the logger ID altogether or adding it as a random effect. You may find it doesn't make too much difference. – Allan Cameron Aug 17 '20 at 15:14

1 Answers1

2

You cannot "remove" a predictor from the predicted means (when you use predict()), because this would return just NA values. Thus, I would recommend the same as Allen and use a sensible / meaningful value at which you hold Logger.ID constant. Here's an example from the ggeffects package:

library(ggeffects)
library(glmmTMB)
data("Salamanders")
m <- glmmTMB(count ~ spp * mined + sample + (1 | site), family = nbinom1, data = Salamanders)

# hold "sample" constant at its mean value
ggpredict(m, c("spp", "mined"))
#> 
#> # Predicted counts of count
#> # x = spp
#> 
#> # mined = yes
#> 
#> x    | Predicted |   SE |       95% CI
#> --------------------------------------
#> GP   |      0.04 | 1.01 | [0.01, 0.27]
#> PR   |      0.11 | 0.60 | [0.03, 0.36]
#> DM   |      0.38 | 0.36 | [0.19, 0.78]
#> EC-A |      0.11 | 0.60 | [0.03, 0.36]
#> EC-L |      0.32 | 0.38 | [0.15, 0.68]
#> DF   |      0.56 | 0.32 | [0.30, 1.04]
#> 
#> # mined = no
#> 
#> x    | Predicted |   SE |       95% CI
#> --------------------------------------
#> GP   |      2.27 | 0.20 | [1.53, 3.37]
#> PR   |      0.46 | 0.33 | [0.24, 0.88]
#> DM   |      2.42 | 0.20 | [1.64, 3.58]
#> EC-A |      0.89 | 0.27 | [0.53, 1.50]
#> EC-L |      3.20 | 0.19 | [2.21, 4.65]
#> DF   |      1.85 | 0.21 | [1.22, 2.81]
#> 
#> Adjusted for:
#> * sample = 2.50
#> *   site = NA (population-level)
#> Standard errors are on the link-scale (untransformed).

# predicted means when sample is set to "0"
ggpredict(m, c("spp", "mined"), condition = list(sample = 0))
#> 
#> # Predicted counts of count
#> # x = spp
#> 
#> # mined = yes
#> 
#> x    | Predicted |   SE |       95% CI
#> --------------------------------------
#> GP   |      0.04 | 1.02 | [0.00, 0.27]
#> PR   |      0.11 | 0.62 | [0.03, 0.36]
#> DM   |      0.38 | 0.38 | [0.18, 0.80]
#> EC-A |      0.11 | 0.61 | [0.03, 0.36]
#> EC-L |      0.32 | 0.40 | [0.15, 0.69]
#> DF   |      0.54 | 0.34 | [0.28, 1.06]
#> 
#> # mined = no
#> 
#> x    | Predicted |   SE |       95% CI
#> --------------------------------------
#> GP   |      2.22 | 0.24 | [1.40, 3.52]
#> PR   |      0.45 | 0.36 | [0.22, 0.90]
#> DM   |      2.37 | 0.24 | [1.49, 3.78]
#> EC-A |      0.87 | 0.30 | [0.49, 1.58]
#> EC-L |      3.14 | 0.22 | [2.04, 4.81]
#> DF   |      1.81 | 0.25 | [1.11, 2.95]
#> 
#> Adjusted for:
#> * site = NA (population-level)
#> Standard errors are on the link-scale (untransformed).

Created on 2020-09-14 by the reprex package (v0.3.0)

Daniel
  • 7,252
  • 6
  • 26
  • 38
  • I have done this and it works fine, thank you for your help! I also ran the model with Logger.ID as a random effect and there wasn't a huge difference between the output compared to the output of the model when predicting with a median logger. – alex Sep 16 '20 at 18:07