Dear all,
I have a more general question. While fitting a non-linear mixed effect model in Stan, which diagnostics method could be employed apart from the plot between observed and predicted values.
I was thinking using some plots on the residuals but not sure how to best do it in Stan.
Currently, my model worked quite fine (plot between observed and predicted values shown quite a good fit). However, I am worried about other assumptions: For example: in my model
y_{ij} = mu_{ij} + \eplison_{ij}
where
\epsilon_{ij} ~ N(0, \sigma^2)
That means I assumed a constant error model. What do you suggest as the best way to do diagnostic based on residuals in Stan?
Calculate the difference between y and mu for each observation and draw from the posterior distribution. This would yield an S by N matrix where S is the number of draws from the posterior distribution. Then calculate the variance (or standard deviation) of those residuals column-wise and see how close it is to a constant.
Or just estimate another model with a Student t likelihood and estimate the degrees of freedom. If the posterior distribution of the degrees of freedom is small, then there is excess heterogeneity in the data-generating process.
Thank you for the reply, Ben.
I plotted the sd of the residuals and it looks like:
Since the definition of “close to a constant” is not straightforward, I am not quite sure in this case can we conclude that the constant variance assumption is satisfied (of course there are some observations that yielded higher variance of the residuals at around 0.55, and almost other values are quite close to each other).
I was thinking about the idea to test the assumption of constant variance of the residuals over time using the population weighted residuals (more details can be see here: http://monolix.lixoft.com/graphics/residuals/). But it is more or less like the idea of computing the studentized residuals in linear regression and plot it over time. Now in my case, when I plot absolute values of the studentized residuals against time, there is some trend (residuals decrease over time):
Moreover, according to testing the normality assumption of the residuals, I got QQ plot:
which shown clearly that the assumption is violated.
Now, since I already tried to fit the model using the log scale, log2 scale or log10 scale (model on original scale did not converge), what would be another strategy for model refining to help improve the normality and homocesdasticity of the residuals?
Sorry I know this is an old post, but have been trying to work out how to compute weighted residuals for a bayesian non-linear mixed effects model by hand, is it a simple division of your residuals by their respective standard deviations that you calculated for your plot? - Many thanks
