BACKGROUND: I'm conducting a hierarchical linear regression using R (specifically R studio, Version 4.1.3). I want to use robust linear models (using the rlm function, MM-estimator) for each of my step, instead of a traditional OLS model (lm function). This is because I have some influential outliers.

1. For example, here is an example of my "STEP 1"/"Model 1" code:

Note: used the f.robftest function from the sfsmisc package to get p-values from beta-coefficients:

model1_controlsRLM = rlm(ER40_CR ~ Age_in_Yrs + Gender2, data = datasetfinal, method = c("MM"))
f.robftest(model2_EFcontrolsRML, var = "Age_in_Yrs")
f.robftest(model2_EFcontrolsRML, var = "Gender21")

2. Example of "STEP 2"/MODEL 2 code:

model2_EFcontrolsRML = rlm(ER40_CR ~ Age_in_Yrs + Gender2 + CardSort_AgeAdj, data = datasetfinal, method = c("MM"))
f.robftest(model2_EFcontrolsRML, var = "Age_in_Yrs")
f.robftest(model2_EFcontrolsRML, var = "Gender21")
f.robftest(model2_EFcontrolsRML, var = "CardSort_AgeAdj)

3. And then comparing models using an ANOVA

anova(model1_controlsRLM, model2_EFcontrolsRML)'

[Missing p-values here, and also not sure how to compare changes in R^2 or an equivalence metric that would communicate similar information]

MY QUESTIONS

• 1. Is it allowed to use an RLM (instead of standard OLS) for each step/model in a hierarchical linear regression?
• 2. If so, how can I calculate changes in adjusted R^2 or an equivalent metric appropriate for RLMs? [I'm open to using a different estimator instead of MM (e.g., M-estimator) if this makes a difference].
• 3. How do I get the p-values of my ANOVA comparing the two RLM models (comparing step 1 vs. step 2)

Thank you for your help in advance.