Lmer violating residuals' normality assumption: What should I do? When "enough data" is enough? I'm trying to plot the following lmer:
mod1 <- lmer(SCORE ~ X1_c * X2 + (1|PARTICIPANT), data = data)

With THIS dataset. (this is a Git link)
However, I can't go past the residual-normality assumption, since no matter what I do, I don't get a normal-distributed residual distribution.
qqnorm(resid(mod1))
qqline(resid(mod1))


### Model's residuals test (1)
shapiro.test(resid(mod1))

Shapiro-Wilk normality test

data:  resid(mod1)
W = 0.95616, p-value = 0.0001235

### Model's residuals test (2)

shapiro_test(resid(mod1))

variable    statistic  p.value
  <chr>           <dbl>    <dbl>
1 resid(mod1)     0.956 0.000123

** It's important to note that: **

*

*I've already tried to exclude outliers

*I've already tried to exclude Participants with Cook's Distance highier than 1

*I've seen MANY MANY posts here on this issue, but I hope this isn't a duplicate since I couldn't apply the other posts' answers to my data

*I've seen some people say that Shapiro-Wilk test don't matter for "large samples", even Andy Field himself in his book and in a video . But what is "big enough" ? I have a total of 148 non-NA paired values. If I can ignore Shapiro-Wilk's results, on WHAT basis can I do that? How could I justify this decision?

*I've seen some people talk about robust modeling with this package or OLS , but I'm a beginner, I'm not familiar with these tests and though I'm really willing to learn about them, I need to return this analysis asap, so I can't do that for the present moment (but I will, any suggestions will be much appreciated)

*Log-transformations of Y or X1_c, in this particular case, wouldn't make much of a sense theorically since both variables are exam/tests scores

*The random intercept normality is ok:

### Random intercept normality test:
r_int<- ranef(mod1)$PARTICIPANT$`(Intercept)`
qqnorm(r_int)
qqline(r_int)
shapiro.test(r_int)

p > 0.05


*

*This is the plot(fitted(mod1), resid(mod1))


*

*I've been trying to decide if I should/can continue with this model within the last 72hours and, honestly, I just wanna to sit down and cry. I really mean that any thoughts would be much appreciated. So if you think that you can help me, I already thank you in advance!

Obs: the dataset I've put on git is the raw data-set with all outliers and Nas. The only thing is that X1_c is a centered predictor around its own mean.

You could consider log-transforming your outcome variable. This would imply a multiplicative relationship between your predictors and $score$. I.e., changing your predictor variables changes score by a certain % amount, rather than by a constant (e.g. 5 points). This can often make sense but it does depend on the scientific context and your beliefs about the manner in which the variables relate to each other.
Alternatively, you could ignore the non-normality of the residuals. According to Gelman and Hill, normality of errors is usually the least important of the regression assumptions. Violations of this assumption, especially of the quite modest type which you show above, have little or no influence on the estimation of the regression coefficients. Only if you wish to use the model to make predictions does this assumption become important.
