Outside of linear regression, what does it mean to say that the residuals should be independent? I understand how to think about residuals in the context of linear regression: we have a dataset of N rows and columns x, y, of linear data $y = mx + \varepsilon$, we pick a linear model $f(x) = \hat m x$ and estimate the $\hat m$ which minimizes the error function (squared residuals). Then for a given $x_i$ the prediction of your model is $f(x_i)$ and the residuals is the set of $e_i = y_i - f(x_i)$. Then the requirement "the residuals should be indenpendent" refers to this set. Great.
Now consider non-linear regression. Specifically I'm trying to get a Poisson regression to work. Here the error function is the maximum likelihood function and the parameter that you're fitting is the usual Poisson $\lambda$. An important point is that in my dataset the $y$ is a count variable (i.e. {0, 1, 2, ....} e.g. number of squirrels you saw that day) and $x$ is a variable that I believe could be a predictor. I incorporate this by imposing that Poisson's $\lambda$ should depend on the predictor in a specific way: $\log \lambda = mx$, and now I'm fitting $m$. After the procedure above we end up with $\lambda(x)$ which then goes into the Poisson if we need to sample. The specifics of the model don't really matter, whether or not these are good ideas is another question, I'm just telling you this to give you context.
Now here's the question: In this context, that does it mean to say independence of residuals? Related question: what are predictions of this model? You know that in the Poisson distribution the mean coincides with the $\lambda$ parameter. For a given $x$, is a prediction of the model the $\lambda$ that we fitted? Or is a prediction one specific instance of sampling from a Poisson with the $\lambda$ that fitted? This matters because in the first case the residuals would be the difference between $y_i$ and $\lambda(x_i)$, whereas in the second case I'm not sure, maybe you would have to compute a histogram and compute the differnce btween that and the Poisson PMF, or something like that?
 A: Error terms in multiple linear regression: The assumptions of linear regression can be framed in the same way as a GLM, without specifying the error terms at all.  To do this, assume that the response variables are conditionally independent given the explanatory variables, with the following conditional distributional form:
$$Y_i | \boldsymbol{x} \sim \text{N}(\mu_i, \sigma^2) \quad \quad \quad \mu_i \equiv \mathbb{E}(Y_i | \boldsymbol{x}) =  \beta_0 + \sum_{i=1}^m \beta_k x_{ik}.$$
This is one way to fully specify the form of the multiple linear regression model.  This distributional form can then be reframed by specifying it in terms of the distribution of $\varepsilon_i \equiv Y_i - \mu_i$, which has the conditional distributional form $\varepsilon_i | \boldsymbol{x} \sim \text{IID N}(0, \sigma^2)$ (i.e., the error terms are independent of the explanatory variables, and are uncorrelated normal values with constant variance and zero mean).  The reason that the error terms are interesting in linear regression is that they are the "deviance residuals" that arise from this model (i.e., they correspond to a measure of the signed change in the log-likelihood).  They have a nice simple distribution because removing the mean from a normal random variable still gives a normal random variable.
Error terms in generalised linear models (GLMs): In the broader class of GLMs you also specify a conditional distribution for the response variables, and a mean function that relates them to the explanatory variable (in a GLM this is a linear function transformed by a 'link' function).  It would be something like this:
$$Y_i | \boldsymbol{x} \sim \text{Dist}(\mu_i, \cdots ) \quad \quad \quad \mu_i \equiv \mathbb{E}(Y_i | \boldsymbol{x}) = g \Big(  \beta_0 + \sum_{i=1}^m \beta_k x_{ik} \Big).$$
This also yields an implicit distribution for the "error terms" $\varepsilon_i \equiv Y_i - \mu_i$, but this is generally not very interesting, since these terms do not always correspond to the deviance residuals in the more general case.  (Their distribution is easy to obtain, but it is not very illuminating.)  Within the context of GLMs we can define analogous quantities:
$$\varepsilon_i \equiv \text{sgn}(y_i - \mu_i) \cdot 2(\log \hat{p}( y_i ) - \log p( y_i | \mu_i, \cdots) ),$$
where $\hat{p}$ is the density for the saturated model (with one parameter per observation).  These latter quantities are the true deviance terms of the model, which generalise the notion of the "error term" in a linear regression.  They have a different form depending on the distribution and link function in the particular model, and their distributional form may be complicated.

Regardless of whether you are working with linear regression or some broader GLM, the true deviance values for the data points are conditionally independent, given the underlying model parameters.  This is a direct result of the assumption that the response variables are conditionally independent.  Of course, the deviance residuals (as opposed to the true deviances) will have some slight correlation owing to the estimation of model parameters, but this will tend to dissipate as $n \rightarrow \infty$.  In linear regression this is manifested in the fact that the covariance matrix for the residual vector depends on the design matrix.  In GLMs there is a similar phenomenon, but the distributions are more complex.
