Am I using the correct statistical analysis? I was wondering if anyone could help me.
I am trying to analyse the effect of body dysmorphia scores (scored on a scale of 0-4) on ratings of agency (on a scale of -100 to 100) under differing contingencies (positive negative, zero). As I have collected data from the general population, the distribution of scores in my sample is very heavily skewed, with most participants having a score of 0 (30+) on the body dysmorphic disorder questionnaire, a few with a score of 1, one with a score of 2, 3 with a score of 3 and 4 with a score of 4. A score of 4 is indicative of Body dysmorphia, so even if I was to collapse the data into positive screening vs negative screening for BDD, the positive group would still only have 4 participants in that group. Thus, I had decided on carrying out 3 separate regression analyses to test the effect of BDD on ratings under each contingency independently. I have found that BDDQ is significant predictor of ratings under each contingency, respectively, with those scoring higher on the BDDQ having more negative ratings under each contingency when tested separately. But I am unsure if this is correct, as having done some reading, and some say to use a mixed factorial ANOVA, but I have an issue with assumption violations (Box’s and Levene’s). Moreover, as I didn’t intend to recruit a high vs low BDDQ group, most of my participants fall into the low group, rather than high (high is a score of 4, which is the maximum score). The repeated measures are contingency. Each participant experiences positive, negative and zero. Each participant experiences these, but they aren’t a measure of the same variable at different time points. Does anyone one have any idea on what I should do? Are the 3 separate regressions sufficient?
Similarly, I have a separate analysis for another study, where I want to add depression as a confound to my regression models such that: ratings ~ BDDQ + BDI. Without BDI, BDDQ is not a significant predictor for ratings. When I add BDI, both BDI and BDDQ are significant predictors of ratings, with BDI being a positive predictor (one point increase in BDI leads to a the beat value of BDI increase in ratings), and BDDQ as a negative predictor (one point increase in BDDQ leads to the beta value decrease in ratings). When I add the interaction between BDI and BDDQ, only BDDQ is a significant predictor, the relationship remains negative. How do I interpret these results? The models with BDDQ + BDI and BDQ + BDI + BDDQ*BDI are significant, the interaction model has the largest R2.
Any help on these two would be very helpful.

For such response variables you need to respect the nature of the data by using a semiparametric ordinal model such as the proportional odds model which generalizes the Wilcoxon and Kruskal-Wallis tests.  See here for resources for this model.
Don't select a model by "significance" or $R^2$.  Fully pre-specify the model using subject matter knowledge.  Otherwise the standard errors and P-values will be next to meaningless.
I am unclear about the repeated measures part of the design, i.e., how many observations there are per subject.  Since you mentioned mixed models I assuming there is > 1 observation per subject.  You can deal with this using a mixed effects ordinal model using the R ordinal or mixor packages or a Bayesian proportional odds model using the rmsb or brms packages.  Bayesian methods work better when random effects are involved.
