Getting adjusted average treatment effect when factors specified as random effects

I am not familiar with the term ATE, but one thing that I find weird is that you seem to ignore the interaction (1| treatment:gender) term in your estimate of the effect. This might be OK for some use cases, but I am not 100% sure it is always appropriate. Using the model you’ve given, I can look at the mean of the interaction coefficients:

inter_coef <- coef(random, summary = FALSE)$`treatment:gender`
hist(rowMeans(inter_coef[,,1]))

And I see that for many samples the average is quite far from zero, so the interaction effects may take a bit of the weight from the main effect… (i.e. they all may be positive for a given treatment so taking only the main effect of treatment would underestimate the total effect you model). This might cancel out on average, but this is AFAIK not guaranteed to happen.

I usually find it easier to make inferences in such cases via the predictions of the model, as this makes it explicit which predictors do I use and in what way.

I’ve written about those issues recently for other people at Inferences marginal of random effects in multi-level models and Multivariate model interpreting residual correlation and group correlation so check if that helps (but feel free to ask for clarifications if that’s not clear).

Best of luck with your model!