Simulating a (simple) robust mixed-effects model to calculate DHARMa residuals

Question

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.