I am planning to simulate a mixed-effects model fitted with robustlmm::rlmer to validate if the model is correctly specified or not by using DHARMa. The model formula can be written as:
fit.robust <- rlmer(y ~ x + (1 | group)) Let us assume that there are 30 groups in the data, $g \in {1, 2, ..., 30}$ and $n$ observations.
In order to unconditionally simulate this model, I am going to set up a simulation, where the conditional mode for each group, $u_g$, is going to drawn from $N(0, \hat\sigma_u)$ for each observation and then weighted by the estimated robustness weight, $w_g$. Similarly, the Gaussian error, $\epsilon_{i}$, will be drawn from $N(0, \hat\sigma_{e} )$ for each observation and weighted by the estimated robustness weight, $w_{i, e}$. Thus, a simulated response can be expressed as follows:
$y_{i, g}$ = $\hat{b}_0$ + $\hat{b}_1x_i$ + $w_{g}$. $u_{i, g}$ + $w_{i, e}$. $\epsilon_{i}$
where $i \in {1, 2, ..., n}$, and $\hat{b}_0$ and $\hat{b}_1$ are estimated intercept and estimated fixed effect for $x$, respectively.
The vectors $w_e$ (of length n) and $w_g$ (with 30 elements) could be obtained by using getME(fit.robust, "w_e") and getME(fit.robust, "w_b"), respectively. For more details, please see the manual of robustlmm.
Would this be a correct way to carry out unconditional simulation i.e. not using estimated conditional modes, to generate response vector for calculating DHARMa residuals? I am aware that I will have to calculate a simulated-response matrix than just a simulated-response vector to get DHARMa residuals.
rlmer. You may consider using thecreateDHARMafunction to generate the residuals since it isn't already supported, as advised here and perhaps contacting Florian about adding it somehow.. $\endgroup$