Why should the variance equal the mean in Poisson regression? I am using Poisson regression to relate road crash data to road geometric characteristics. Why should the variance equal the mean in Poisson regression?
Thanks in advance
 A: Because it is a consequence of the functional form of the Poisson distribution that mean and variance are equal. If this condition is not met the model is inadequate and alternatives may be considered such as negative binomial regression (this is called overdispersion). See:
https://en.wikipedia.org/wiki/Poisson_distribution
A: To see this, let's consider the number of crashes for a specific road characteristic. Let's say this number follows a Poisson distribution with a mean $\mu$. This mean is for a certain number of km driven so let's introduce the rate $\lambda$, say 1 crash per km and the total number of km driven $T$. One assumption of Poisson distribution is that the rate remains constant over the total distance driven, consequently, $\mu=T \times \lambda.$
We divide the number of km driven into tiny $N$ short intervals of size $h$, so short that each sub-interval contains at most one crash. Now, the probability you see a crash in this tiny interval is like flipping a coin. We will denote this probability as $p$. This is known as a Bernoulli distribution and we will take for granted that the variance is $p \times (1 - p)$. On the other hand, we learnt earlier that the rate $\lambda$ is constant so we expect to see $\lambda \times h$ event in this sub-interval, that is, $p=\lambda \times h$. 
Now, if we assume that the probability to see a crash in this tiny sub-interval is extremely low then $1 - p$ approaches 1 (e.g. consider $h=\text{1 meter}$). We learnt earlier that the variance for Bernoulli distribution is $p \times (1-p)$ and if $p$ is extremely low, then $p \times (1-p) \simeq p=\lambda\times h.$ This is quite interesting because we just demonstrated that both mean and variance are equal to $\lambda \times h$ in this tiny sub-interval. 
If you extend this approach to $n$ consecutive intervals (like flipping coins $n$ times), you'll get something called a binomial distribution and in this case, the mean is $np$ and the variance $np(1-p)\simeq np$ when $p$ is small. To get to the point, for $N$ consecutive intervals of size $h$ with extremely low $p$, the mean and the variance are equal.
Now, in practice, this is usually not the case in observational studies. The reason is that we cannot take into account all the factors for heterogeneity in the study. For example, the mean number of accidents may differ during the day time and the night time. Yet, if we were to aggregate both without accounting for the different factors, the marginal variance may become larger than what we expect. This is called an overdispersion.
A: Poisson regression allows for inference of regression parameters - generally a parameter vector, $\mathbf{\theta}$ - under the assumption that the errors in the model are distributed according to a Poisson distribution: $\epsilon \sim \operatorname{Poisson}(\lambda)$.
This is appropriate for modelling certain data, a good example being heteroscedastic count data. Summarising an example from German Rodriguez' notes on Poisson Models for Count Data, we can consider a scenario wherein we are regressing the number of children a woman has, $\mathbf{y}$, onto predictors like her age, education level, region of residence, duration of marriage and so on, $\mathbf{X}$. In this situation, where the units of analysis are individual women, we can plot the relationship of the mean and variance of the outcome / target variable $\mathbf{y}$:


Clearly, the assumption of constant variance is not valid. Although the variance is not exactly equal to the mean, it is not far from being proportional to it.
Thus, we conclude that we can do far more justice to the data by fitting Poisson regression models than by clinging to ordinary linear models [where we regress assuming error distributed according to a singly-parametrised Gaussian].

The question then becomes why the mean, or Expectation, of the random error, $\epsilon$, is equal to its Variance, i.e. why $\mathbb{E}[\epsilon] = \mathbb{Var}[\epsilon]$ for $\epsilon \sim \operatorname{Poisson}(\lambda)$ in your regression model.
To see why these are equal, we can directly derive them.
Expectation of the Poisson errors
For the Expectation, from the definition of expectation we have from the definition of the Poisson distribution:
$$\mathrm{E}(X)=\sum_{k \geq 0} k \frac{1}{k !} \lambda^{k} e^{-\lambda}$$
Then:
$$
\begin{aligned}
\mathrm{E}(X) &=\lambda e^{-\lambda} \sum_{k \geq 1} \frac{1}{(k-1) !} \lambda^{k-1} & & \text { as the } k=0 \text { term vanishes } \\
&=\lambda e^{-\lambda} \sum_{j \geq 0} \frac{\lambda^{j}}{j !} & & \text { putting } j=k-1 \\
&=\lambda e^{-\lambda} e^{\lambda} & & \text { Taylor Series Expansion for Exponential Function } \\
&=\lambda & &
\end{aligned}
$$
Variance of the Poisson errors
For this we exploit the identity:
$$\operatorname{var}(X) =\mathrm{E}\left(X^{2}\right)-(\mathrm{E}(X))^{2}$$
So starting with the computation of $\mathrm{E}\left(X^{2}\right)$ we have:
$$
\begin{array}{rlr}
\mathrm{E}\left(X^{2}\right) & =\sum_{k \geq 0} k^{2} \frac{1}{k !} \lambda^{k} e^{-\lambda} & \text { Definition of Poisson Distribution } \\
& =\lambda e^{-\lambda} \sum_{k \geq 1} k \frac{1}{(k-1) !} \lambda^{k-1} & \text { Change of limit: term is zero when } k=0 \\
& =\lambda e^{-\lambda}\left(\sum_{k \geq 1}(k-1) \frac{1}{(k-1) !} \lambda^{k-1}+\sum_{k \geq 1} \frac{1}{(k-1) !} \lambda^{k-1}\right) \\
& =\lambda e^{-\lambda}\left(\lambda \sum_{k \geq 2} \frac{1}{(k-2) !} \lambda^{k-2}+\sum_{k \geq 1} \frac{1}{(k-1) !} \lambda^{k-1}\right) & \text { Change of limit: term is zero when } k-1=0 \\
& =\lambda e^{-\lambda}\left(\lambda \sum_{i \geq 0} \frac{1}{i !} \lambda^{i}+\sum_{j \geq 0} \frac{1}{j !} \lambda^{j}\right) & \text { putting } i=k-2, j=k-1 \\
& =\lambda e^{-\lambda}\left(\lambda e^{\lambda}+e^{\lambda}\right) & \text { Taylor Series Expansion for Exponential Function } \\
& =\lambda(\lambda+1) & \\
& =\lambda^{2}+\lambda \\
& =\lambda^{2}+\lambda-\lambda^{2} & & \\
& & &
\end{array}
$$
Then putting this together:
$$\begin{aligned}
\operatorname{var}(X) &=\mathrm{E}\left(X^{2}\right)-(\mathrm{E}(X))^{2} \\
&=\lambda^{2}+\lambda-\lambda^{2} \\
&=\lambda
\end{aligned}$$
You can see that the equality of the expectation (or mean) and variance emerges from the definition of the Poisson distribution (which by the way, emerges from taking a limiting case of the Binomial distribution when $n \rightarrow \infty$ and so necessarily $p \rightarrow 0$).
Both of the above proofs, for Expectation and Variance are available on ProofWiki.
