Perform linear regression, but force solution to go through some particular data points

I know how to perform a linear regression on a set of points. That is, I know how to fit a polynomial of my choice, to a given data set, (in the LSE sense). However, what I do not know, is how to force my solution to go through some particular points of my choosing. I have seen this be done before, but I cannot remember what the procedure was called, let alone how it was done.

As a very simple and concrete example, let us say that I have 100 points scattered on the x-y plane, and I choose to fit a whatever-order polynomial through them. I know how to perform this linear regression very well. However, let us say that I want to 'force' my solution, to go through, say, three of my data points at x-coordinates $x=3$, $x=19$, and $x=89$, (and their corresponding y-coordinates of course).

What is this general procedure called, how is it done, and are there any particular pitfalls I need to be aware of?

Edit:

I would like to add, that I am looking for a concrete way to do this. I have written a program that actually does the linear regression in one of two ways, inverting the covariance matrix directly, or through gradient descent. What I am asking is, how, exactly, step by step, do I modify what I did, such that I force the polynomial solution to go through specific points?

Thanks!

-
Why do you call it "linear" if you are using a polynomial? Every point you want it to pass through is a constraint that will reduce your degree of freedom. You can then use a constrained optimization algorithm. –  curious_cat Feb 20 '13 at 22:41
It is linear because you are finding co-efficients to a linear combination. For example, if you want to fit your data to a cubic, then you are finding the co-efficients (the $c$'s) of $y = c_0 + c_1x + c_2x^2 + c_3x^3$. –  Tarantula Feb 20 '13 at 22:44
@Mohammad: One other way to approximate what you want would be to use a weighted least squares solution, and give very large weights to the points that you want the regression line to pass through. This should force the solution to pass very closely to the points that you choose. –  Jason R Feb 21 '13 at 2:59
@JasonR Good to see you on here. Yes WLS is indeed an interesting contender. I have gone with whubers answer because of the clever polynomial factorization, and because it maintains the error structure nicely. –  Tarantula Feb 21 '13 at 18:30

The model in question can be written

$$y = p(x) + (x-x_1)\cdots(x-x_d)\left(\beta_0 + \beta_1 x + \cdots + \beta_p x^p \right) + \varepsilon$$

where $p(x_i) = y_i$ is a polynomial of degree $d-1$ passing through predetermined points $(x_1,y_1), \ldots, (x_d,y_d)$ and $\varepsilon$ is random. (Use the Lagrange interpolating polynomial.) Writing $(x-x_1)\cdots(x-x_d) = r(x)$ allows us to rewrite this model as

$$y - p(x) = \beta_0 r(x) + \beta_1 r(x)x + \beta_2 r(x)x^2 + \cdots + \beta_p r(x)x^p + \varepsilon,$$

which is a standard OLS multiple regression problem with the same error structure as the original where the independent variables are the $p+1$ quantities $r(x)x^i,$ $i=0, 1, \ldots, p$. Simply compute these variables and run your familiar regression software, making sure to prevent it from including a constant term. The usual caveats about regressions without a constant term apply; in particular, the $R^2$ can be artificially high; the usual interpretations do not apply.

(In fact, regression through the origin is a special case of this construction where $d=1$, $(x_1,y_1) = (0,0)$, and $p(x)=0$, so that the model is $y = \beta_0 x + \cdots + \beta_p x^{p+1} + \varepsilon.$)

Here is a worked example (in R)

# Generate some data that *do* pass through three points (up to random error).
x <- 1:24
f <- function(x) ( (x-2)*(x-12) + (x-2)*(x-23) + (x-12)*(x-23) )  / 100
y0 <-(x-2) * (x-12) * (x-23) * (1 + x - (x/24)^2) / 10^4  + f(x)
set.seed(17)
eps <- rnorm(length(y0), mean=0, 1/2)
y <- y0 + eps
data <- data.frame(x,y)

# Plot the data and the three special points.
plot(data)
points(cbind(c(2,12,23), f(c(2,12,23))), pch=19, col="Red", cex=1.5)

# For comparison, conduct unconstrained polynomial regression
data$x2 <- x^2 data$x3 <- x^3
data$x4 <- x^4 fit0 <- lm(y ~ x + x2 + x3 + x4, data=data) lines(predict(fit0), lty=2, lwd=2) # Conduct the constrained regressions data$y1 <- y - f(x)
data$r <- (x-2)*(x-12)*(x-23) data$z0 <- data$r data$z1 <- data$r * x data$z2 <- data$r * x^2 fit <- lm(y1 ~ z0 + z1 + z2 - 1, data=data) lines(predict(fit) + f(x), col="Red", lwd=2) The three fixed points are shown in solid red--they are not part of the data. The unconstrained fourth-order polynomial least squares fit is shown with a black dotted line (it has five parameters); the constrained fit (of order five, but with only three free parameters) is shown with the red line. Inspecting the least squares output (summary(fit0) and summary(fit)) can be instructive--I leave this to the interested reader. - whuber, this is interesting... I would be lying if I said I fully understood it just yet, but I am digesting it as we speak. If I understand correctly, here I am basically solving for the$\beta$s as usual, but they are being multiplied with$r(x)x^i$, instead of just$x^i$s as before, yes? If that is correct, then how are you computing$r(x)$exactly? Thank you. – Tarantula Feb 20 '13 at 23:34 I have added a worked example, Mohammad. – whuber Feb 20 '13 at 23:42 Oh, perfect. I will study it. Using your example, it would still be possible to force the poly to go through points that are part of the data though, right? – Tarantula Feb 20 '13 at 23:47 Absolutely that can be done: but be doubly cautious about interpreting the p-values or any other statistics, because now your constraints are based on the data themselves. – whuber Feb 20 '13 at 23:53 Your post had me up last night. I taught myself the LIP. (LIP is interesting. It is like a Fourier decomposition but with polys). – Tarantula Feb 21 '13 at 17:55 If you want to force a regression line to go through a single point, that can be done in a roundabout way. Let's say your point is$(x_i,y_i)$. You just re-center your data with that point as the origin. That is, you subtract$x_i$from every$x$-value, and$y_i$from every$y$-value. Now the point is at the origin of the coordinate plane. Then you simply fit a regression line while suppressing the intercept (forcing the intercept to be (0,0). Because this is a linear transformation, you can easily back-transform everything afterwards if you want to. If you want to force a line to go through two points in an X-Y plane, that's pretty easy to do also. Any two points can be fit with a line. You can use the point-slope formula to calculate your slope, and then use one of the points, the slope, and the equation of a line to find the intercept. Note that it may not be possible to fit a straight line through three points in a coordinate plane. However, we can guarantee that they can be fit perfectly with a parabola (that is, using both$X$and$X^2$). There is algebra for this too, but as we move up, it might be easier to just fit a model with software by including only those three (more) points in the dataset. Similarly, you could get the straight line that best approximates those three points by fitting a model that has access only to those three points. I feel compelled to mention at this point, however, that this may not be a great thing to do (unless your theory provides very solid reasons for doing so). You may also want to look into Bayesian regression, where you can allow your model to find the best combination of the information in your data, and some prior information (which you could use to strongly bias your intercept towards zero, for example, without quite forcing it). - Gung, thanks for your answer. I have modified my question a little. I did not know about Bayesian Regression, but will take a look at it. I am afraid I do not completely understand how, concretely, from an algorithmic view point, the one point and two point case you mentioned. Specifically, for the one point, I understand removing and re-adding$x_i$and$y_i$to each point before and after a block, but I do not understand how to do that block exactly. For the 2-point case, I am afraid I do not understand at all what to do there. Thanks. – Tarantula Feb 20 '13 at 22:37 Although throwing in three more points and weighting them (a la Glen_b's answer) could create such a fit, interpreting any of the statistical output would be problematic: some adjustments would be needed. – whuber Feb 21 '13 at 0:01 To add a little extra information to @gung's excellent coverage of the linear case, in the higher order polynomial case there are several ways you could do it either exactly or approximately (but pretty much as accurately as you need). First, note that the degrees of freedom for the polynomial (or indeed of any fitted function) must be at least as large as the number of "known" points. If the degrees of freedom are equal, you don't need the data at all, since the curve is completely determined. If there are more 'known' points you can't solve it. From here on, I'll just talk about when the polynomial has more d.f. than the known points (such as a cubic - with 4df and three known points). 1) "the curve must pass through this point" is a linear constraint on the parameters, resulting in constrained estimation or constrained least squares (though both terms can include other things than linear constraints, such as positivity constraints). You can incorporate linear constraints by either (a) recasting the parameterization to implicitly include each constraint resulting in a lower order model. (b) using standard tools that can incorporate linear constraints on the parameters of a least squares fit. (usually via something like the formula given at the above link) 2) Another way is via weighted regression. If you give the known points arbitrarily large weight, you can get essentially the same fit as in (1). This is often readily implemented, can be substantially quicker than reparameterizing, and can be done in packages that don't offer constrained fitting. All of @gung's caveats apply - Glen_b, I had not considered weighted regression. It might be the way to go about it. I have put it in my to do list. I believe I can teach myself that without incident. Regarding (1), may you kindly expand on this aspect of re-paramaterization? Also, what do you 'call' this that I am trying to do, where I force the polynomial to go through certain points? Part of the problem is that I do not know what to google for. If I know what this was called, I might be able to augment what you are saying with online material. Thanks. – Tarantula Feb 20 '13 at 23:11 See my edits above, which include some search terms and a link with a few more details. – Glen_b Feb 20 '13 at 23:23 +1 Weighted regression is a nice idea. Some adjustment of the output statistics, such as estimates of RMS error, might be needed. – whuber Feb 20 '13 at 23:58 @whuber +1 indeed, if the statistics (like$s^2$,$F$,$R^2\$ ... std errors etc) are intended to relate only to the not-"known" points (which is likely what one would want), aside from the parameter estimates and the fitted values, the raw stats coming out will nearly all be wrong. I originally typed a sentence relating to that but I seem to have deleted it before posting; it's important to mention it. –  Glen_b Feb 21 '13 at 0:08
Thanks for your answer Glen_b, although I have accepted @whuber 's, I still learned a lot from yours. –  Tarantula Feb 21 '13 at 18:28