alternative to ANOVA or simplifying approximations for multi-factor,, batch design with binomial data (advanced beginner) this is my first post. I did already read various didactic materials and related posts here, but still had trouble figuring this one out.  Sorry this is so long but I am trying to be as specific as possible, and also mindful that I might use technical terms incorrectly, so I'm spelling things out rather literally to be sure.
Context:
I'm dealing with a rescue operation - a set of directed hypothesis test experiments were already planned, done, and  analyzed, but the statistical method (2-way ANOVA with replicates, the field's standard analysis for this type of experiment) has been called into question at the manuscript review stage. For future experiments we definitely want to know the ideal approach for such experiments, and will take the time to learn it or get a collaborator to help us. But for this case time is of the essence, so I'd also like to identify the simplest approach that is valid and adequate. The claims being made are regarding very large effects and large interactions, so it's probably not essential to maximize statistical power of a hypothesis test, just to find one that is justifiably valid or at least conservative.
These experiments have two factors, and the inference of interest is whether or not these factors interact and in which direction (whether Y enhances, inhibits, or does not influence the effect of X).  Y may or may not have an effect on its own. To give a specific example: the first factor is nominal (+ vs - intervention X), the second factor is also nominal (potential modulators Y1, Y2, or none). So there are 2x3 treatment groups. There are 3 observations per group. BUT a 2-way ANOVA with replicates is questionable choice, for a couple of reasons:
(1) The raw data are binomial. The three "observations" are percentages (% cells that are positive).  The reviewer thinks since the data are categorical we should use Fishers' exact, but I am not sold. Each of those percentages is based on scoring a large (but variable) number of randomly sampled cells (100-300). I am inclined to argue that with this N, a percentage is effectively continuous, so it's not invalid use a method that assumes a continuous quantitative measure, although we do lose information about the uncertainty of those estimates. I was initially inclined to argue that these percentages are also approximately Gaussian, such that it's not invalid to use a method that assumes normality. But a few of the treatment groups have percentages near enough to 0  to be necessarily non-Gaussian.  So this would be a reason to prefer some non-parametric alternative, or parametric one based on the underlying binomial distribution.  Unless one can argue that this particular deviation from normality could only make the ANOVA overly conservative. If that's a known theorem, then the fastest solution to our present problem may be to make that argument, and then report the ANOVAs already done.
But there's still another problem.
(2) There were 3 independent runs: observation #1 for all six conditions was made on cells derived from animal 1, obs #2 from animal #2, etc.  No specific reason to expect the animals differ, but it's a potential confound. So if I'm thinking straight, these aren't "replicates" in the sense required, and we need to do some sort of three factor or batch or nested ANOVA or multi-level hierarchical model (none of which we know how to do).  After reading several sources, I am not confident about identifying which approach(es) is/are appropriate, particularly among nonparametric ones. Also, if we don't have the chops to do the analysis that is required, nor any way to get a professional statistician's assistance on such short notice, then it would be highly desirable to simplify the problem into one we know how to do or at least is fairly easy to learn. For example, we considered pre-processing the data into "X-effect" scores that are computed by comparing + vs - X within each Y condition within-experiment. That turns the problem back into two factors (Y group x animal). But that opens up questions as to how to properly compute "effects" (do we subtract, divide, normalize, etc.). Alternatively, maybe we can do some separate analysis to show that animal identity or batch isn't a substantial source of variance, and then justify that the 2-way-ANOVA-with-replicates framework is an acceptable approximation.
(3) Important context: nobody thinks the observed effect could have arisen by chance. We just don't know how to quantify that, and industry standard in this field requires that a frequentist p value be reported. I don't have the raw count data at hand, but if I simulate binomial data using the observed percentages and a representative distribution of Ns, and plot symbols for the means and binomial confidence intervals for each of the 3 experiments, for each of the six conditions, the effects between conditions are huge compared to the binomial confidence intervals, and within each treatment group the differences between the three experiments are tiny. So I suspect the observed effects and interactions would have a very low p value even within each of the three experiments considered separately, if we knew such a test. (The raw count data can be found, but the current data file only contains the percentages).
 A: Retrieve the original count data.
Then treat this as a binomial logistic regression. In R with the glm() function, you can provide a two-column outcome specifying the number of successes and failures for each combination of predictors. You specify the values of X,Y and animal in the corresponding row of the data.
Forget trying to force this into a specific acronym like ANOVA. This is just a (generalized) linear model. You specify X and Y and their interaction(s) as predictors. I'd recommend coding the data so that the absence of intervention is the baseline category for X and "no modulator" is the baseline category of Y.
You include animal as another categorical predictor in the model. If you include an additive term for animal you allow for differences in percentages under baseline conditions among animals. In R the model then would be something like:
glm(c(success, fail) ~ X*Y + animal, data=yourData, family = binomial)

That will give a set of interaction coefficients for each of Y1 and Y2 with X, to use for testing your hypotheses.
Your sense that normal approximations could work OK in this case might be correct, but if I understand your study properly then the above approach should be unimpeachable.
