4

i have the following data and created a model with the package glmmTMB in R for plant diameters ~ plant density (number of plants) with a random plot effect:

d <- data.frame (diameter = c(17,16,15,13,11, 19,17,15,11,11, 19,15,14,11,8),
                      plant_density = c(1000,2000,3000,4000,5000, 1000,2000,3000,4000,5000, 1000,2000,3000,4000,5000),
                      plot = c(1,1,1,1,1, 2,2,2,2,2, 3,3,3,3,3))

glmm.model <- glmmTMB(diameter ~  plant_density + (1|plot),
                        data = d,
                        na.action = na.omit,
                        family="gaussian",
                        ziformula = ~ 0)

My intention was to create a plot with predicted diameter data for different plant densities with an included random plot effect. So i tried to predict the data:

new.dat <- data.frame(diameter= d$diameter,
                      plant_density = d$plant_density,
                      plot= d$plot) 

new.dat$prediction <- predict(glmm.model, new.data = new.dat, 
                              type = "response", re.form = NA)

Unfortunately I get an output for every plot but wanted a generalized prediction for the diameter ~ plant density.

My goal is to create a plot like here, but with a regression model from glmmTMB which consider the random effect.

Thanks for ur help!

Schneiderhansl
  • 89
  • 1
  • 8
  • when you say "consider the random effect"- do you mean you want different curves for different levels of d$plot ? – SushiChef Jan 06 '22 at 13:23
  • Hi! No, I wanted to create new data according to the measured ones (f.e. for 1 to 8000 plants) so I can plot one line to show the general context. The plots just are three areas (sampling method) where I measured my data. – Schneiderhansl Jan 06 '22 at 13:29
  • 1
    is this what your looking for: plot(ggeffects::ggpredict(glmm.model)) – SushiChef Jan 06 '22 at 13:36
  • Thanks for that tip. Thats what I'm looking for. I think the random effects are considered :-) – Schneiderhansl Jan 06 '22 at 14:41
  • One last question: Can you tell me how to change the limits of the x axis from 0 to 20? I don't know how to include it in the code? – Schneiderhansl Jan 06 '22 at 14:44
  • do you mean the y-axis? there's nothing in that range for your x-axis – SushiChef Jan 06 '22 at 15:36
  • Let us [continue this discussion in chat](https://chat.stackoverflow.com/rooms/240797/discussion-between-schneiderhansl-and-sushichef). – Schneiderhansl Jan 06 '22 at 16:46

1 Answers1

4

The ggeffects package makes this type of thing very easy to implement and customize.

For example

library('ggplot2')
library('glmmTMB')
library('ggeffects')
d <- data.frame (diameter = c(17,16,15,13,11, 19,17,15,11,11, 19,15,14,11,8),
                 plant_density = c(1000,2000,3000,4000,5000, 1000,2000,3000,4000,5000, 1000,2000,3000,4000,5000),
                 plotx = as.factor( c(1,1,1,1,1, 2,2,2,2,2, 3,3,3,3,3)))

glmm.model <- glmmTMB(diameter ~  plant_density + (1|plotx),
                      data = d,
                      family="gaussian")

# basically what your looking for
plot(ggpredict(glmm.model, terms = "plant_density"))

# with additional a change of limits on the y-axis
plot(ggpredict(glmm.model, terms = "plant_density")) + 
scale_y_continuous(limits = c(0, 20))

You can really do whatever you'd like with it from there, changing colors, themes, scales, the package has some nice vignettes as well.

SushiChef
  • 588
  • 3
  • 12
  • Thanks for ur help!! At least I wanted to predict values from the model. I measured densities from 1000 to 5000 which helped me to create a model. Now I wanted to predict new plant diameters for density values between 0 and 10000. I tried to change the x axis but R shows only the regression for the given densities. Is there any way to solve this problem with that package? – Schneiderhansl Jan 06 '22 at 17:41
  • I'm not sure I follow you, if your asking for a prediction interval for new levels: I don't think its quite as straight forward, but is easy enough even with predict(new_dat, type = interval), it might be answered somewhere else on SO or you can ask another question – SushiChef Jan 06 '22 at 18:34