How to determine which independent variables to add in a multivariable regression model when sample size is small As a post-doc I am working on some data that I did not collect myself. 
The central question I'm trying to answer: is there a difference in the fat mass of neonates born to mothers with gestational diabetes mellitus (GDM), versus mothers with normal glucose tolerance? 
I thus want to create a model where fat mass is the dependent variable, and GDM (coded 0 for mothers without GDM, and 1 for those with) is the independent variable. I also need to control for other variables such as gestational weight gain, pre-pregnancy BMI, gestational age, etc., as these could potentially confound the relationship between neonatal fat mass and GDM.
Sample size is 72 for GDM mothers/babies and 211 for non-GDM mothers/babies. I'd also like to run the model taking out all missing data for comparison, which would leave me with 68 GDM mothers/babies and 112 non-GDM mothers/babies. 
I know I need to consider the sample size when determining how many IVs to add to the model, and I've searched online and found some general rules of thumb for how many IVs are appropriate in multivariable regression with a given sample size. It appears that the data set I'm working on is too small to accommodate the number of potential confounders, especially for the smaller diabetes sub-sample.
Therefore, I'm stuck, as I cannot do what I would like to do - collect more data - and I have about 10 IVs that I think should be added into the model based on their potential to confound the association of interest. 
Another thing to note: I recently read a paper discussing how it is inappropriate to test bivariable associations of potential IVs with your outcome of interest in order to determine whether they should be added into the model (based on their bivariable significance). So I've already ruled that out as an option for determining which IVs should make it into the model. 
Does anyone have any advice on how to determine which/how many IVs to add that goes beyond what I've already considered? 
Thanks so much.  
 A: I would say your sample isn't that small and can definitely accommodate 10 predictors if entered without interactions with the treatment. If this really concerns you, you can estimate a propensity score as the predicted probability of GMD given the predictors and simply fit a model regressing the outcome on GMD and the propensity score (though it would be important to allow this model to be flexibly estimated). 
If you're worried some of your predictors are not true confounders, you can use the group lasso with doubly robust estimation described in Koch, Vock, & Wolfson (2018), which estimates a doubly robust treatment effect (robust to misspecification of either the propensity score model or outcome model) and uses lasso to select only the relevant covariates (so the models aren't overly saturated). This was specifically designed to select only the covariates from a large pool of potential confounders that are helpful in estimating the treatment effect.
For missing data, it would inappropriate to drop cases. You should use multiple imputation or FIML estimation to retain your original sample size. It does not significantly complicate the problem to use either of these methods (though FIML cannot be used with propensity scores).
