What is the correct way to fit sigma.formula in gamlss() to add a per-dependent variable variance?

I have data in the format:

 | organism | treatment | replicate | response |
| ---      | ---       | ---       | ---      |
| A        | X         | 1         | 20       |
| A        | X         | 2         | 5        |
| A        | Y         | 1         | 10       |
| A        | Y         | 2         | 12       |
| A        | Z         | 1         | 8        |
| A        | Z         | 2         | 6        |
| B        | X         | 1         | 20       |
| B        | X         | 2         | 5        |
| B        | Y         | 1         | 10       |
| B        | Y         | 2         | 12       |
| B        | Z         | 1         | 8        |
| B        | Z         | 2         | 6        |


Where 'organism', 'treatment', and replicate are all categorical variables. There is also a 'length' variable which represents the exposure rate for each organism. The 'response' variable appears to follow a negative binomial distribution, where the variance is not equal to the mean. I have previous modelled this data using a multilevel model where 'organism' is the random effect:

glmer.nb(response ~ treatment + (1|organism) + offset(log(length)) - 1, data = ...)


I removed the intercept in the model, because it is uninterpretable to have an expected response with no organism or treatment.

This model however has large amounts of heteroscedasticity where the fitted versus standard residuals shows an obvious trend. Checking with dharma also indicates overdispersion in the residuals (if this is the correct way to describe the two-sided uniform test outlined in the dharma manual?). Plotting residuals versus treatment shows differing amounts of residual variance for each treatment, suggesting that variance should also be modelled based on the dependent variable.

I have followed up with the gamlss package and fit a model with dependent-variable specific variance as:

gamlss(formula       = response ~ treatment + random(organism) + offset(log(length)) - 1,
sigma.formula = ~ organism + treatment,
data = ...)


Comparing sigma.formula = ~ 1 and sigma.formula = ~ organism + treatment models using BIC suggests I should accept the latter model as being more parsimonious. The heteroscedasticity has also disappeared, where there residual versus fitted plot now appears normally distributed along the y-axis 0.

Is this the correct way to fit variance as a dependent variable in a model? I've only had previous experience with glms and mixed effect models, so I don't have intuition about what seems sensible when using gamlss.