Shows how to use coefs and their vcov matrix to get simulated means and SE- will work for any model! And this example uses a categorical variable. Using the model.matrix avoids having to understand the treatment contrasts.

sim_coefs.R

library(mvtnorm)
library(dplyr)

#make cyl a factor, like site
mtcars$cyl <- as.character(mtcars$cyl)


#fit a model
mod <- lm(mpg ~ cyl*hp, data = mtcars)


#get 10K draws of the coefficients
coef_sims <- rmvnorm(1e4, coef(mod), vcov(mod))


#10k rows by 6 cols
dim(coef_sims)

#now, let's get the site effect at the average level of HP
dat <- data.frame(cyl = unique(mtcars$cyl),
                  hp = mean(mtcars$hp),
                  mpg = 1:3) #just need some values

#ok, we need to turn this into a model frame
dat_frame <- model.matrix(mpg ~ cyl*hp, data = dat)

#the calculation!
mpg <- coef_sims %*% t(dat_frame)

#how many sims? 
dim(mpg) #10k rows, 3 cols - 1 for each cyl

#add simulated predictions to data frame
dat$mpg <- colMeans(mpg)
dat$se_mpg <- apply(mpg, 2, sd)


#plot
library(ggplot2)
ggplot(dat,
       aes(x = cyl, y = mpg, ymin = mpg-se_mpg, ymax = mpg+se_mpg)) +
  geom_pointrange()