How to run backward elimination in $R$ with both categorical data and numeric data? [closed]

Usually in backward elimination, we start with the full model of all covariates, check the $$p$$-value of $$t$$-statistic for each covariate (which is compared between the full model and the model minus the given covariate), and remove the covairate with the highest $$p$$-value until the maximum of $$p$$-value is all below some threshold $$\alpha$$.

What about when the data has both categorical (of multiple levels) and numeric data? I know the most natural thing is to do a similar procedure, but now with the $$t$$-statistic replaced with the $$F$$-statistic, and compare the entire group of covariates against the full model to eliminate (we remove the categorical variable of all levels all at once). But how can I easily implement this in $$R$$ and if there is already some package that integrates the scenario. In the usual case of all numeric covariates, I can perform something like

run_backward_elimination = function(alpha){
S = 1:p
while(TRUE){
pvalues = summary(lm(highway_mpg~[,S]))$coefficients[-1,4] if(max(pvalues) <= alpha){ break } remove_ind = S[which.max(pvalues)] S = setdiff(S,remove_ind) } XS = X[,S,drop=FALSE]; colnames(XS) = S print(summary(lm(Y~XS))$coefficients)
}


Then it should be fine, as I can directly check the $$t$$-static at the summary table of linear regression. But with the Anova table, the $$F$$-statistic is compared among the inner models, not against the full model like $$t$$-static does. So I have to manually run run $$p$$ many regressions to remove one covariate, which has a rather high complexity and is also tedious to implement.