Why is my simulated Poisson regression mean not equal to the variance? I simulated a simple poisson regression model as follows
$x \sim Normal(2, 1)$
$\mu = \exp(0.6 + 0.9 x)$
$y \sim \text{Poisson}(\mu)$
set.seed(123)
n = 1001
x = rnorm(n, 2)
mu = exp(0.6 + 0.9 * x)
sim <- rpois(n, lambda = mu)

However, the mean and variance are not equal
mean(sim)
# 16.68731
var(sim)
# 349.5331

but the mean is quite close to the standard deviation:
sd(sim)
# 18.69581

I expected the mean to equal the standard deviation. Why is that not occurring here? How can I fix it?
I eventually want to be able to simulate a Poisson regression with multiple continuous covariates and at least 2 categorical covariates, Possibly with hierarchical effects. So any code to extend to that case would be great as well.

Your simulated data sim are not Poisson. That is, they are not unconditionally Poisson, but only conditionally on the predictor. Thus, equidispersion does not hold - or, more precisely, equidispersion will hold conditionally on the predictor: if you simulate many times for a fixed value of x (or equivalently, mu), the result will be equidispersed. But then, that will not come as a surprise.
What you have is a case of a compound distribution, specifically, a Poisson-lognormal distribution: the parameter to a Poisson distribution is itself lognormally distributed.
