# 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.

• This is confused on several levels. The main problem is already explained in the first answer. Here is one more: Even for an unconditional Poisson, the mean does not equal the SD unless exceptionally the SD and the variance coincide at 1. Commented Jul 27, 2022 at 16:13

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.