Executive summary: this is much harder than you would expect, and there is neither a standard implementation, nor even an accepted estimator.
Let's fix some terminology. In a survey, you have primary sampling units, secondary sampling units, and so on for potentially multiple stages. In a mixed model you have top-level clusters, possibly clusters at intermediate levels, and first-level individual observations.
Stata's -mixed-, and the gllamm package before it, have a method for the (fairly common) setting where the clusters and the sampling units are the same (or at least nested).
You're talking about a setting where they aren't. When the first complex-survey multilevel models were described in 1998, Jon Rao pointed this issue out in the discussion
The authors assume that the sample is selected according to the
hierarchical structure of the model, but in multipurpose surveys the
hierarchical structure of the sample could be quite different. For
example, an educational multilevel model may include, as levels,
students, classes, schools and school boards. However, the sample of
students may be chosen by a multistage design where the stages are
geographical areas, then households and then students; for example the
Canadian National Longitudinal Survey of Children and Youth follows
such a design. Since the design clusters cut across the model clusters
for such surveys, the theory proposed may not be applicable. We are
currently studying this problem.
He didn't get very far; it's a hard problem.
One approach to this problem is to use pairwise likelihood. I'm working on software for this, but it's not ready yet; Xudong Huang worked on the theory in his PhD thesis, based on earlier work by Jon Rao and Grace Li and co-workers.
On top of the lack of user-friendly software, a straightforward implementation of pairwise likelihood gives quite inefficient estimators for the variance components -- and if you aren't interested in the variance components, you might as well be doing ordinary regression. On top of that, there is currently no approach to computing BLUPs of the random effects, which are needed for adaptive Gaussian quadrature to fit generalised linear mixed models.
Overall, I have some hopes for an implementation of this being available in R by the end of 2020, but it isn't yet.
