I see at least three reasonable options, although there is one I tend to do.
- Compute the difference score, $D = Pre - Post$ and then predict that. The regression equation being something like: $$\hat{D}_{i} = b_0 + b_1X_{i} + b_2Z_{i} + b_3Genderism_{i} + b_4Z_{i}*Genderism_{i}$$ One thing that is attractive about this is that it is straightforward to do.
- Predict the post scores using your model, but also controlling for pre scores. This regression equation would look something like: $$\hat{Post}_i = b_0 + b_1X_{i} + b_2Z_{i} + b_3Genderism_{i} + b_4Z_{i}*Genderism_{i} + b_5Pre_{i}$$ This tends to be my preferred approach. It does not require you to compute any new variables (which is not a big issue but can be slightly annoying). It also includes in the results an estimate of the strength of association between pre and post scores. It works because your treatment, $X$, predicts that portion of the post scores that are not explained by pre scores. The only drawback I see is that if you have a very small sample size, you lose an additional degree of freedom controlling for pre scores. That seems reasonable to me; however, because there is measurement error at both time points anyway, so it is not like computing the difference scores has zero error.
- The last approach I see is to reshape the data from wide to long, and fit a mixed effects model. This would look something like: $$\hat{Outcome}_{ij} = b_0 + u_{0i} + b_1X_{i} + b_2Z_{i} + b_3Genderism_{i} + b_4Z_{i}*Genderism_{i} + b_5time_{ij} + b_6time_{ij}*X_{i}$$ The outcome is the score at the jth time point for the ith individual. The model includes a random intercept (captured by $u_{0i}$) for each individual. The time captures the change over time, and the interaction between time and the group variable is the "treatment" effect. Although kind of cool, I think this model is far too much work in the simple case where you only have pre and post scores (if you had 3+ time points, it would make sense).
By the way, I do not particularly think that you need to center your variables prior to using them as moderators. I know many people teach that, but the models should work out the same. It can be slightly convenient as it makes the simple effects potentially more interpretable (0 = mean, therefore they are the effect of the variable at the mean of the other), but aside from that I see little gain. In more complex models, the reduction in collinearity between the variables and their product can also be helpful, but I have only ever see that matter in complicated random effects models or in some parallel latent growth models I fit once.
