I'm modeling the effect of a treatment on a population of flies. For each fly, I have the following covariates:

• Treatment (control or treated)
• Sex
• Cage in which the fly was reared. I have 3 cages of treated flies and 3 cages of control flies.
• 96-well plate on which the fly was sequenced. Each plate is entirely treatment or entirely control flies, from a mix of the 3 cages of that treatment. (i.e., plate is nested within treatment, and cage is nested within treatment, but plate and cage are not nested in each other).

and I also have my response variable, which is continuous and determined by the sequencing.

I'd think to model this using a mixed-effects model, treating treatment and sex as fixed effects and cage and plate as random:

Response ~ treatment + sex + 1|cage + 1|plate

I have no issues with this so far, but I'm interested in applying a permutation test to this data and am not sure how best to do so. I want to build a null distribution of test statistics for the treatment fixed effect. As I understand it, I could just shuffle treatment labels to get a null distribution:

Response ~ treatment_permuted + sex + 1|cage + 1|plate

However, I'm concerned that the random effects here will be based on the variance associated with the true treatment, and will not be truly permuted. Alternatively, I could use the permuted treatment labels in making new cage:treatment and plate:treatment groups for the random effects, but then I'd have more of these groups than in my true case, since plate and cage will no longer be nested with permuted treatment.