In our assignment we were asked to model the compressive strength of concrete (response variable: Strength) with predictor variables Cement (kg/m^3), Water (kg/m^3), and Coarse.Aggregate (kg/m^3). We had to use backwards stepwise model selection from a complete second order model to produce the final model.

#complete second order model
full.model <- lm(Strength ~ (Cement + Water + Coarse.Aggregate)^2 +
                I(Cement^2) + I(Water^2) + I(Coarse.Aggregate^2),
                data = concrete.df)

The statistical analysis was carried out in RStudio and below is the output from the stepwise model selection

# backwards stepwise model selection
step(full.model, direction = "backward")

# output of the final step
Step:  AIC=519.54
Strength ~ Cement + Water + Coarse.Aggregate + I(Water^2) + Water:Coarse.Aggregate

                         Df Sum of Sq   RSS    AIC
<none>                                16004 519.54
- I(Water^2)              1    691.74 16695 521.77
- Water:Coarse.Aggregate  1   2490.36 18494 532.00
- Cement                  1   3105.32 19109 535.27

I attempted to determine if this model was the best fit for the data by looking at the Anova table

# anova table
anova.lm(final.model)

# output 
Analysis of Variance Table

Response: Strength
                       Df  Sum Sq Mean Sq F value    Pr(>F)    
Cement                  1  2612.8 2612.80 15.3468 0.0001696 ***
Water                   1  2126.3 2126.27 12.4890 0.0006370 ***
Coarse.Aggregate        1   945.3  945.27  5.5522 0.0205336 *  
I(Water^2)              1    34.3   34.26  0.2012 0.6547702    
Water:Coarse.Aggregate  1  2490.4 2490.36 14.6276 0.0002355 ***
Residuals              94 16003.6  170.25                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Based on the Anova table I concluded that the I(Water^2) term was not significant given the p-value > then the threshold 0.05 and removed that term from the final model.

# final model
final.model <- lm(Strength ~ Cement + Water + Coarse.Aggregate +
                      Coarse.Aggregate:Water, data = concrete.df)

The assignment was marked and the feedback I was provided mentioned that the Anova is unnecessary and I shouldn't combine using AIC and p-value variable selection within the same model building exercise.

My question is: why is this the case? I have tried looking for resources to provide the reasoning behind this with no luck.

Note: this is a second year statistics course and it may be outside of the scope of what we are learning right now but I need to know! Thanks!

Sorry if this is not the right forum to ask, let me know if there is a better place to ask this question.