# How can I fit a polynomial given some prior knowledge but fewer observations than coefficients?

I want to fit a second degree polynomial to data from a fractional factorial design in five dimensions. Because of the cost of trials, I am initially only going to do eight runs. So I have fewer observations than the full number of first and second degree terms, including interaction terms.

My goal is to find a local optimal point or region in the vicinity of my trials (if there is one) or direction of greatest increase (if there isn’t). I am saving budget for a few (I hope) confirmatory tests.

For some, but less than half, of the total number of polynomial coefficients I have a very rough idea of their magnitude from a literature review. For others I know the sign, but not the value.

It is my belief that I can estimate the coefficients using some kind of penalized regression. But I don’t know how to incorporate my prior knowledge, and I don’t know whether I should alter my experimental design in some way based on that knowledge (e.g. by structuring is so that the aliased interaction terms are ones where I have some prior information).

I have done a lot of kinds of regression over the years, but I have not looked at anything Bayesian since I was an undergraduate, lo these many years ago. So I’d like a response that’s a little more of a how-to than “Oh, use a Bayesian LASSO” or whatever. Extra karma points for an answer with some pointers toward an R implimentation.

• Optimize some objective function, possibly including some penalty terms, subject to constraints. If you know the sign of a coefficient, make that a constraint. If you have a range of possible values for a coefficient, impose lower and upper bounds on it. If you don't have hard bounds, include penalty factor in the objective which penalizes the objective (maybe not quadratic) if the coefficient's too far from a reasonable range. You might be able to formulate and solve as Quadratic Programming (QP) problem, i.e., quadratic objective and linear constraints. There are QP solver functions in R. – Mark L. Stone Aug 18 '16 at 1:43
• You use the term "runs" for your experiments. Am I right in assuming that you are performing "computer experiments", i.e., running a costly simulation code for each point in your design? If so, then your results are deterministic, i.e., running twice the same design point (replication) leads to exactly the same response, which is of course not the case for a physical experiment. In this case it would be better to use a space filling design (optimal LHS, orthogonal, uniform designs, etc.) than classic DOE designs such as fractional/full factorial, central composite, etc. – DeltaIV Aug 18 '16 at 7:43
• Hi DeltaIV! No, "runs" was just a poor word choice. These are experiments -- in fact, they are experiments on myself, of some drug/supplement combinations. I've got my doctor's go-ahead to adjust within a limited range. He's sort of interested in the idea, I think. The "cost" is mainly in time to completion. – andrewH Aug 18 '16 at 8:36
• you might want to look at eg Bayesian optimisation - which aims to update the function as you do experiments (AFAIK) en.wikipedia.org/wiki/Bayesian_optimization – seanv507 Sep 5 '16 at 13:45

In $d=5$ dimensions, a second degree polynomial has $\frac{(d+1)(d+2)}{2}=21$ coefficients, which you want to estimate with 8 runs. With such a tight budget, I definitely wouldn't go for a training/test set approach. With cross-validation, you will be able to use all your runs for fitting the final model. This way, you may be able to increase the number of runs a little bit (10, maybe?). In your situation, anything helps.

Having said that, fitting a polynomial with a number of coefficient larger than the size of your design, or in general fitting an underdetermined linear model, is really the perfect job for LASSO, ridge regression or elastic net (in case you don't know whether your "ideal" model is likely to be sparse or not). The glmnet function of the R package glmnetallows you to set lower.limits and upper.limits on each coefficient: this allows you to easily impose sign constraints. For example, suppose you have four coefficients in your model: let the first one be the intercept $\beta_0$, and assume that $\beta_1$ is known to be positive, $\beta_3$ is known to be negative while we have no prior knowledge on $\beta_2$. The intercept is never shrinked in shrinkage methods, nor it can be constrained (at least with the glmnet function). Thus you can only specify constraints for $\beta_i,\ i=1,\dots 3$. You will set

lower.limits <- c(0, -Inf, -Inf)
upper.limits <- c(Inf, Inf, 0)


For the order of magnitude constraints, you could introduce lower/upper bounds as before. Another option would be to specify different penalty factors for each coefficient with the parameter penalty.factor. A larger penalty factor will shrink the coefficient towards a smaller absolute value (or even have it dropped out of the model, for LASSO). Suppose you have again 4 coefficients: as noted before, the intercept is not included in the penalty term and thus cannot be shrinked. This leaves us with three coefficients for which you can specify penalty.factor. Since the default value is 1, then if you want the last coefficient to be larger in module than the other ones, you can speficy

penalty.factor <- c(1, 1, 0.1)


There are two drawbacks with this approach: the first one is that, even assuming that all predictors have similar mean values and standard deviations, a penalty.factor of 0.1 doesn't mean that $\beta_3$ will be 10 times as larger as the other two coefficients. This depends on the data set and also on whether you're using ridge regression (which penalizes the squared coefficients), LASSO (which penalizes the absolute values of the coefficients) or elastic net (which is a combination of both). If we drop the assumption that all predictors have similar means and standard deviations, then the relationship between penalty.factor and the coefficient size is even more complicated, because by default glmnet standardizes all predictors before applying shrinkage. The reason is explained in paragraph 6.2.1 of An Introduction to Statistical Learning by James, Witten, Hastie and Tibshirani. This means that the penalty is applied to the coefficients of the standardized predictors, but of course glmnet returns the coefficients in the original (unstandardized) scales. Thus, unless your predictors have similar ranges initially, a penalty.factor as in my example will not necessarily return a $\beta_3$ bigger than $\beta_1$ and $\beta_2$.

Hope this helps: if you want a more "hands-on" answer, please provide a reproducible example.

• Wow! Excellent answer! Love all the practical advice. – andrewH Aug 19 '16 at 2:28