I used the glmulti package in R to generate a linear model with interactions. The response variable has 21 observations, set against 5 input variables, one of which is a binary factor while the rest are continuous. The "best model" produced by glmulti, as judged according to aic, had 6 predictor terms (4 main effects and two interactions). The model fit beautifully in terms of fit (R2, all terms are significant p<0.01) and diagnostics, but the degrees of freedom are a problem.
Each of the 21 observations is the average of about 4 observations. This is to say that there were 21 sites with four sub-samples per site. I then re-fit the same model with the sub-samples nested within sites with a GLMM with a gaussian distribution. The same results were produced as before. The coefficients still only have 13 degrees of freedom, however the number of observations for the model in its entirety has jumped up to around 75.
Do I bypass the problem of having too few degrees of freedom by opting for a mixed model? I reviewed some of the other answers concerning df in GLMM but have yet to discern a clear answer to this problem.