Generalized linear model: variance vs mean

Question

I have just read this blog: Linear Models, ANOVA, GLM etc

The author tries to explain in which situations it is better to use a generalized linear model instead of linear regression. At a certain point, he takes some data that represents the count of warms in a field, under different treatments $T_i$. The author applies a generlized linear model with a Poisson family (exponential link) and gets some results. Until here everything is fine.

Than he says:

One of the assumptions of the Poisson distribution is that its mean and variance have the same value.

and he checks the mean of all the observations vs the variance of all the observations. These turn out to be different and the author states that a quasipoisson family is better.

And here is where I don't agree. Isn't the mean the same as the variance given a single factor? I mean, if I take treatment $T_i$ and I compare the mean number of warms, then I would expect it to be the same as the variance. But the mean across all treatments must not necessarily be the same as the variance across all treatments. Am I missing something?

Answer 1

I think you are right in that the Poisson regression (like any other regression) is defined conditionally on the covariates (and thus so are the implied assumptions).
If $Y$ is the outcome and $X$ a set of covariates, the assumption is that $Y \mid X$ follows a Poisson distribution.

Then, given $X$, the Poisson assumption implies that $\mathbb E(Y \mid X) = {\rm Var}(Y \mid X)$.

But it's not because $$\mathbb E(Y \mid X) = {\rm Var}(Y \mid X)$$ that we necessarily have $$ \mathbb E(Y) = {\rm Var}(Y). $$ This is because $$ \mathbb E(Y) = \mathbb E [ \mathbb E(Y\mid X) ] $$ while $$ {\rm Var}(Y) = { \rm Var }[ \mathbb E(Y \mid X)] + \underbrace{\mathbb E[ {\rm Var }(Y \mid X )]}_{=\mathbb E(Y)} $$ And thus $\mathbb E(Y) \neq {\rm Var}(Y) $ whenever ${ \rm Var }[ \mathbb E(Y \mid X)] > 0$.

---

For example if $X$ is a simple binary covariate such that $\mathbb P(X=1) = 0.5$ and $$ Y \mid X=0 \sim \mathcal{P}(\lambda=1) $$ $$ Y \mid X=1 \sim \mathcal{P}(\lambda=2) $$ then, by assumption, $$ \mathbb E(Y \mid X) = {\rm Var}(Y \mid X). $$ However we have $$\mathbb E(Y) = 0.5 \times 1 + 0.5 \times 2= 1.5$$ while \begin{align*} {\rm Var}(Y) &= \underbrace{0.5 \times 0.25 + 0.5 \times 0.25 }_{{ \rm Var }[ \mathbb E(Y \mid X)]} + 1.5 \\ &=1.75 \end{align*}

data<-data.frame(do.call(rbind,lapply(1:1e5,function(i){
X<-rbinom(1,1,0.5)
Y<-rpois(1,lambda=ifelse(X==0,1,2))
return(c(X,Y))
})))
mean(data[,2])  # around 1.5
var(data[,2])   # 1.75