florianhartig / DHARMa

Diagnostics for HierArchical Regession Models
http://florianhartig.github.io/DHARMa/
200 stars 21 forks source link

Understanding fittedPredictedResponse and simulations #387

Closed lumen-roses closed 10 months ago

lumen-roses commented 10 months ago

Hello Florian,

I am trying to understand how does DHARMa generates the fittedPredictedResponse values and simulations matrices for a Poisson Regression model. I am still learning a lot about GLMMs, and I still haven't scratched the surface of some basic, important concepts, that is why I try to replicate everything DHARMa generates as I work through my data. I will work through a specific example. Consider the following model:

$$ Y{ijk} \sim \text{Pois}(\Lambda{ijk}) $$

$$ \log(\Lambda_{ijk}) = X_i \beta + U_j + V_k + \log(\Delta t_i)$$

$$ U_j \sim \mathcal{N} (0,\sigma_u) $$

$$ V_k \sim \mathcal{N} (0,\sigma_v) $$

$$ U_j + V_k \sim \mathcal{N}(0, \sqrt{\sigma_u^2 + \sigmav^2}) = \mathcal{N}(0,\sigma{uv})$$

Here $X_i$ represents one of the rows of a matrix $X$ (model matrix, no intercept column) and $\beta$ is a (column) vector of coefficients. $U_i$ and $V_k$ are normal random effects. The term $\log(\Delta ti)$ is a time offset. I know that the expected value of the response $Y{ijk}$, conditional on zero random effect values is:

$$ E[Y_{ijk}|U_i = 0, V_j = 0] = \exp(X_i \beta + \log(\Delta t_i))$$

This is the part where I think my reasoning might be wrong. In order to test the structure of the model as a whole, I think that the fittedPredictedResponse should be the total expected value of the response $Y{ijk}$, that is $E[Y{ijk}]$. That means obtaining the unconditional expectation of the response, I think:

$$ E[Y_{ijk}] = \exp(X_i \beta + \log(\Delta ti) + \sigma{uv}^2/2) $$

Good simulations from the model would be obtained by generating Poisson rate parameters from a lognormal distribution, and then integers from a Poisson distribution. The marginal distribution for $Y_{ijk}$ (Poisson lognormal) is the average over all possible random effect values.

What is the definition of "unconditional" residuals and predictions as seen in #43 and in the main vignette? I would think that the only way to obtain unconditional predictions or simulations is by simulating ALL levels in the model, and such simulations could be obtained by marginalizing (if possible) out some levels of the hierarchy. If I don't do this, and condition on the fitted values or zero, wouldn't I get a pattern in the residuals?

I am currently fitting models like this one using glmmTMB and Stan.

Is this reasoning correct? Where could I be wrong? Thanks for any support you can give me.

florianhartig commented 10 months ago

Hello @lumen-roses,

first of all, in all DHARMa documentation I use those words in the following meaning:

It is true that for a nonlinear model, the expectation value of the unconditional simulation is not equal to the unconditional prediction, and thus you could additionally define what I would call a "marginal prediction". Actually, if you create a DHARMa object via createDHARMa with STAN simulations, this it the default, i.e. median posterior simulations are used as predictions.

That being said, I don't see a particular reason to prefer this definition of the current one - the main point of #43 is that you don't want to plot residuals against fitted REs, because the latter spuriously correlate with the residuals. Whether you use marginal or unconditional predictions should roughly result in the same patterns.

lumen-roses commented 10 months ago

I see. This makes sense. Thank you!

With that clarified, I will close this issue.