What exactly needs to be independent in GLMs?

Question

I am reading up on Poisson Regression.

Say I have an input vector $X^{T} = (X_1, X_2, \ldots, X_p)$ and I want to predict an output $Y$ which is correlated to these inputs.

Then the mathematical form of a Poisson Regression Model is:

$$\log(Y) = \alpha + \beta_{1}x_1 + \beta_{2}x_2 + \ldots + \beta_{p}x_p$$

where $\alpha, \beta_{i}$ are numeric coefficients.

Now, one of the assumptions is that the observations must be independent of one another. I'm not quite sure what this means?

Does it mean that each individual response variable $x_1, x_2, \ldots , x_p$ have to be independent, or does it mean that each $X_{i}$ vector in the data set has to be independent of one another?

Answer 1

What is actually required is conditional independence of the response variable. Conditional on the regressors, that is. A Poisson regression model - for independent data - is no different from an ordinary least squares model in this assumption except that the OLS conveniently expresses the random error as a separate parameter.

In a Poisson regression model, the specific form of the conditional response is debatable due to disagreements in the fields of probability and statistics. One popular option would be $Y/\hat{Y}$, which are still random variables albeit not exactly Poisson distributed - the value is considered an ancillary statistic like a residual, which does not depend on the estimated model parameters. An interesting thing to observe here is that, if the linear model is misspecified (such as omitted variable bias), this may induce a kind of "dependence" on the residue of the fitted component of the model that is not correctly captured.

Independence has a very specific probabilistic meaning, and most attempts to diagnose dependence with diagnostic tests are futile. This is mostly complicated by the important mathfact that independence implies covariance is zero, but zero covariance does not imply independence.

Poisson regression and OLS are special cases of "generalized linear models" or GLMs so we can conveniently deal with the independence of observations in GLMs. The classic OLS residual plot, which we use to detect heteroscedasticity, is effective to visualize a covariance structure, but additional assumptions are needed to declare independence. In a GLM such a Poisson, a non-trivial mean-variance relationship is an expected feature of the model; a standard residual plot would be useless. So we rather consider the Pearson residuals versus fitted as a diagnostic test. In R, a simple poisson GLM(x <- seq(-3, 3, by=0.1); y<-rpois(length(x), exp(-3 + x)); f <- glm(y ~ x, family=poisson). The resulting graph shows a tapering curve of residuals and, arguably a funnel shape, with a LOESS smoother showing a mostly constant and 0-expectation mean residual trend.

As an example, we may describe the distribution of the $x$ as the design of the study. While many examples here and in textbooks deal with $X$ as random, that's merely a convenience and not that reflective of reality. The design of experiment is to assess viral load of infected mouse models treated with antivirals at a sequence of dose concentrations, say, control, $X = (0 \text{ control}, 10, 50, 200, 1,000,$ and $5,000)$ mg/kg. For an effective ART, the sequence of viral loads is expected to be descending because of a dose-response relationship. The outcome might look $Y = (10^5, 10^5, 10^4, 10^3, LLOD, LLOD)$. This response vector is not unconditionally independent, there is a strong "autoregressive" trend induced by the design. Trivially, when the fitted effect is estimated through the regressions, the conditional response is completely mutually independent.

A more involved but real life example is detailed in Agresti's categorical data analysis second edition in Chapter 3 covering poisson regression. This deals with the issue of estimating the number "satellite crabs" in a horseshoe crab nest, a sort of interesting polyamory. Data analyses can be found here.

Answer 2

A general form of expressing a model is

$$Y_i = f(\textbf{X}_i,\boldsymbol\beta,\epsilon_i)$$

The function $f$ expresses the random variable $Y_i$ in terms of a random latent variable $\epsilon_i$, a fixed/known regressor variable $\textbf{X}_i$, and some distribution parameters $\boldsymbol\beta$.

Note:

• Here the subscript $i$ refers to the observation within the sample. In your question you have a subscript relating to the elements in the vector $\textbf{X}^{T} = (X_1, X_2, \ldots, X_p)$. You could also write it with two subscripts $\textbf{X}_i^{T} = (X_{1i}, X_{2i}, \ldots, X_{pi})$ where $\textbf{X}_i$ is a matrix and the row index relates to the observation and the column index relates to the regressor/feature.

• It is the random part $\epsilon_i$ that is assumed to be independent. (often, more complex models can assume some dependency between the different $\epsilon_i$)

• If you perform an experiment then the regressor variables $\textbf{X}_i$ can be dependent. For example you might repeat a measurement with the same values for several $\textbf{X}_i$. But what needs to be independent is the random part $\epsilon_i$.

Example:

Say we have

$$Y_i = Q(\mu = a + b x_i, p= \epsilon_i)$$

where we use $\epsilon_i$ a uniform distributed variable, $Q(\mu,p)$ is the quantile function of the Poisson distribution with $\mu$ the mean and $p$ the quantile.

Then some simulated data, with parameters $a = 10$ and $b=1$, could look like

       x    epsilon  y
 [1,] 10 0.26550866 17
 [2,] 10 0.37212390 18
 [3,] 10 0.57285336 21
 [4,] 10 0.90820779 26
 [5,] 20 0.20168193 25
 [6,] 20 0.89838968 37
 [7,] 20 0.94467527 39
 [8,] 20 0.66079779 32
 [9,] 30 0.62911404 42
[10,] 30 0.06178627 31
[11,] 30 0.20597457 35
[12,] 30 0.17655675 34

An r-code to compute the above numbers is:

set.seed(1)
a = 10
b = 1
x = c(10,10,10,10,20,20,20,20,30,30,30,30)
epsilon = runif(12) # generate the nois part based on uniform distribution 
y = qpois(epsilon,a + b * x) # transform to the y variable using the quantile function
plot(x,y)
cbind(x,epsilon,y)