Larger-than-expected variance in a Poisson I have data set that contains values between 10-25.
c(21, 19, 12, 20, 25, 22,14, 14, 22, 11, 15, 14,
12, 13, 16, 23, 20, 16, 17, 16)

We have been asked to give reasons that might cause a larger-than-expected variance.
I've said that if we were to get an outlier of 40 in the data set, that this will cause the probability of a number occurring to be higher, meaning that the variance will also be higher.
Is this a correct statement to make? 
 A: The mean of the $n = 20$ observations is 17.10 and their variance is 17.04.
Thus if the data are Poisson, it seems reasonable to estimate the
population mean $\lambda$ as $\hat \lambda = 17.1.$
I did a chi-squared goodness-of-fit test with the five 'categories'
0-13, 14-16, 17-19, 20-22, and 23+, which have counts $X = (4, 7, 2, 5, 2)$
and under the hypothetical distribution $\mathsf{Pois}(17.1),$ they
have probabilities $p = (.194, .264, .270, .172, .100).$ 
Thus the
expected counts are $E = np = (3.88, 5.28, 5.40, 3.44, 2.00).$
Finally, the chi-squared statistic is
$$Q = \sum_{i=1}^5 \frac{(X_i-E_i)^2}{E_i} = 3.41.$$
Under the null hypothesis, $Q$ is approximately distributed as $\mathsf{Chisq}(5-2 = 3),$
which has 95th percentile (5% critical value) 7.815. 
So the
hypothesis that the data are from $\mathsf{Pois}(17.1)$ is not
rejected. [That is not saying a lot, because the power of this goodness-of-fit
test with so little data is not very good.]
A histogram of the data along with Poisson probabilities (dots)
is shown below. Of course, one cannot expect a really good fit of the histogram to the
exact distribution with only $n = 20$ observations, but the match is
not so bad that the null hypothesis is rejected.

Sampling from a Poisson distribution with mean 17.1, a value as large as 40 would be extremely unlikely and would make one wonder
whether the resulting data could come from that null distribution: For $X \sim \mathsf{Pois}(17.1),$
$P(X \ge 40) \approx 1.6 \times 10^{-6}.$
1 - ppois(39, 17.1)
## 1.626658e-06

Certainly, changing the observation 25 to 40 would increase the sample
mean by a little (from 17.1 to 17.85) and the sample variance by a lot
(from 17.04 to 40.77), as you say.
mean(x1); var(x1)
[1] 17.1
[1] 17.04211
x2 = c(x1[1:19], 40)
mean(x2); var(x2)
[1] 17.85
[1] 40.76579

A: Statistical tests on your dataset will not tell you what caused the observed values to be what they are, and so if you get an outlier in your data, and you want to known the cause, that is something you will have to examine by looking at the source of your data.  Notwithstanding this limitation, one way of examining over-dispersion within a Poisson model is to generalise this to a negative binomial model and then undertake explicit hypothesis testing with regard to the over-dispersion parameter using a likelihood-ratio test.  To do this, we posit the generalised distribution:
$$X_1, ..., X_n \sim \text{IID NB} \Big( \frac{\lambda}{\phi}, \frac{\phi}{1+\phi} \Big),$$
where $\lambda > 0$ is the mean parameter and $\phi > 0$ is a dispersion parameter.  With this formulation we have $\mathbb{E}(X_i) = \lambda$ and $\mathbb{V}(X_i) = (1 + \phi) \lambda$, and in the special case where $\phi \rightarrow 0$, the distribution reduces to the Poisson distribution.  Thus, for completeness, we can take $\phi=0$ to be a valid input for the parameter, and the distribution is the Poisson distribution in this case.

LR test for over-dispersion: Suppose we observe a data vector $\mathbf{x} = (x_1,...,x_n)$ and we use this to calculate the full and restricted MLEs:
$$(\hat{\lambda}, \hat{\phi}) \equiv \underset{\lambda > 0, \phi \geqslant 0}{\arg \max} \text{ } \text{ } \ell_\mathbf{x}(\lambda, \phi) \quad \quad \quad \hat{\lambda}_0 \equiv \underset{\lambda > 0}{\arg \max} \text{ } \text{ } \ell_\mathbf{x}(\lambda, 0).$$
We can use a standard likelihood-ratio test to test the hypotheses:
$$H_0: \phi = 0 \quad \quad \quad H_\text{A}: \phi \neq 0.$$
Under the null hypothesis we have $LR(\mathbf{x}) \equiv 2 \cdot (\ell_\mathbf{x}(\hat{\lambda}, \hat{\phi}) - \ell_\mathbf{x}(\hat{\lambda}_0, 0)) \sim \text{Chi-Sq}(1)$, with higher values of the test-statistic being more conducive to the alternative hypothesis.  Hence, the p-value for the test is:
$$p(\mathbf{x}) = \int \limits_{LR(\mathbf{x})}^\infty \text{Chi-Sq}(r|1) dr.$$

Application to your data: You can implement the above test on your data by using the MASS package in R.  The required code is as follows:
#Input the data
DATA <- c(21, 19, 12, 20, 25, 22,14, 14, 22, 11, 15, 14, 12, 13, 16, 23, 20, 16, 17, 16);

#Calculate the LR statistic
library(MASS)
LOGLIKE_NEGBIN  <- fitdistr(DATA, "Negative Binomial")$loglik;
LOGLIKE_POISSON <- fitdistr(DATA, "Poisson")$loglik;
LR_STATISTIC    <- 2*max(LOGLIKE_NEGBIN - LOGLIKE_POISSON, 0);

#Output the LR-statistic and p-value of LR test
LR_STATISTIC
[1]  0

1 - pchisq(LR_STATISTIC, 1);
[1]  1

As you can see, in this particular case the log-likelihood is maximised under the restricted model for the Poisson distribution, so there is no evidence of over-dispersion in the data (p-value of one).  We can repeat this analysis with the addition of one more data point of 40, to see the effect this has on the test.
#Input the updated data
DATA2 <- c(DATA, 40);

#Calculate the LR statistic
library(MASS)
LOGLIKE_NEGBIN  <- fitdistr(DATA2, "Negative Binomial")$loglik;
LOGLIKE_POISSON <- fitdistr(DATA2, "Poisson")$loglik;
LR_STATISTIC    <- 2*max(LOGLIKE_NEGBIN - LOGLIKE_POISSON, 0);

#Output the LR-statistic and p-value of LR test
LR_STATISTIC
[1]  6.466906

1 - pchisq(LR_STATISTIC, 1);
[1]  0.01099017

As you can see, with this additional data point, there is now reasonably strong evidence of over-dispersion in the model (p-value a little over one-percent).  The test does not allow you to identify the cause of the observed data value, but it does show that the addition of this 'outlier' adds evidence of over-dispersion.
A: 
I've said that if we were to get an outlier of 40 in the data set,
  that this will cause the probability of a number occurring to be
  higher, meaning that the variance will also be higher.

Possibilities:


*

*Your Poisson statistical assumption is invalid. Maybe something like negative binomial with overdispersion is better. The true variance might be higher than what you would expect.

*Your Poisson statistical assumption is valid, you are just lucky to get the number

*Your Poisson statistical assumption is valid, but you should mix it with another model (mixed modelling)

