I am trying to build a multiple regression model using R. I have a number of predictor variables. I have some basic domain knowledge for which I am trying to build the model. To start with, I included a few predictor variables based on domain knowledge and high correlation coefficients with the response variable, while excluding some other predictors due to multicollinearity. I would like to figure out if I should include some interaction terms. But, due to large number of predictors, I am having a hard time trying to figure out which all interaction terms I should include in the model. Based on what I have read on this site about automated model selection (thanks, @gung et. al), I am trying to avoid using it.
I think you can deal with some of these issues based on your domain knowledge. 21 predictors aren't a lot with 11,000 records, if your outcome variable is some continuous measure, so the issues you face are what predictors and interactions to include and how to deal with collinearity.
For building the model, you might not want to omit any of your 21 original predictors. When you omit 1 of 2 highly correlated predictors, you are throwing out information provided by the one you omit and run the risk of your results being too closely tied to the peculiarities of those correlated variables in the particular sample that you are analyzing. Also, don't depend on correlation of independent variables with your dependent variable for choosing predictors to include. Keeping some predictors poorly correlated with the dependent variable might help improve the performance of other predictors, even in the absence of interactions.
For interactions, consider adding interactions that you think might be important based on your domain knowledge. That presumably will be a lot fewer than the 420 possible 2-way interactions among 21 predictors so that you will still have a reasonably small number of independent variables. You might even consider not including any interactions at all and seeing if the 21 predictors on their own work well enough for your purposes. Sometimes it's best to start simple, and add complexity only as needed.
One way to deal with collinearity would be based on domain knowledge: combine correlated predictors into a single predictor that captures the essential underlying phenomenon that those correlated predictors represent. That would seem to be consistent with your goal to use your model for inference. If you can combine correlated predictors in a way that's defensible based on domain knowledge, you might reduce the number of predictors in the model in a way that makes inference easier.
Alternatively, to deal with collinearity you could use an approach like ridge regression that tends to treat collinear predictors together. My impression is that ridge regression is more often used for predictive rather than for inferential models, but it does have the advantage of handling collinearity in a reasonable way. It returns coefficients for all predictors, which is either an advantage or a disadvantage depending on your perspective. Some might prefer LASSO for inference as it retains only a subset of predictors, but its particular choice among collinear predictors might be sample dependent and you would have to consider that in interpreting the results.
My guess is that a bigger problem than dealing with 21 predictor variables will be finding appropriate scaling transformations for your variables so that they work reasonably well in the approximation of a linear model.