6. Non-linear Regression
(Back to Top)

We can also directly fit the Michaelis-Menten function to our data using non-linear regression. Remember the term "sum-of-squares" from your old regression class? When you fit a regression model, you get a "fitted value" for every data point used to fit the model. If you take the difference between the fitted value and the observed value, you get what we call a residual. Then if you square all the residuals and add them up, you get the residual sum-of-squares. The smaller that is, the better the model fits your data. You may have heard this called the least-squares method. Well, non-linear regression works the same way. With non-linear regression, we specify the form of the model we want to fit and the parameters that need to be estimated. S-PLUS then searches for parameter values that will minimize the residual sum-of-squares.

You can try non-linear regression from the pull down menus as well. From the Statistics menu, select Regression > Nonlinear... to open the dialog box. Having been through it at the command line, you should have no trouble knowing what to do in the dialog box.

7. Polynomial Regression
(Back to Top)

The last method we'll use to fit these data is polynomial regression, where the model takes on the form y = b₀ + b₁x + b₂x² + b₃x³ + × × × , etc. We're going to fit a second order polynomial like this:

> polyfit2_lm(vel~poly(conc,2), data=df)

We're using the same function lm() as linear regression, but we're also using the function poly() to generate our polynomial formula. We could just as easily fit this model this way:

> polyfit2_lm(vel~conc+conc^2, data=df)

This is fine when you have a second order polynomial, but what about when you want to fit a 10^th order, or a 20^th order polynomial? Life can be simpler sometimes using the function poly(). But there is one distinct difference between these two approaches. If you look at the summaries for both, you'll find that all the significance tests result in the same P-values, but the coefficients themselves are different. The difference is in what the poly() function does. Instead of producing a polynomial with x and x², you get an orthogonal and normalized form of the polynomial. Suffice it to say that you get appropriate results from the regression (i.e., the significance tests are correct), but the resulting coefficients are on a different scale. So before we look at the coefficients, we want to transform them back to the original scale. This is how we do that:

> coef2_poly.transform(poly(df$conc,2), coef(polyfit2)) > coef2 x^0 x^1 x^2 1.288439 25.65224 -56.50026

Now we want to draw this line on our graph. We could add it to the plot with the other two lines, but let's create a new graph and label the x and y axes:

> plot(df$conc, df$vel, xlim=c(0,0.4), ylim=c(0,5), + xlab="Substrate Concentration", ylab="Reaction Rate")

There's something new with this line. You can enter the entire command on a single line if you want, but if you hit Enter before the command is complete, you get the "+" prompt on the second line where you finish the command. NOTE: the "+" on the second line IS NOT part of the command, it is the prompt to continue. So if you enter this all together on a single line, DO NOT include the "+".

There are a couple ways to generate the y-values for the line. Perhaps the most straightforward is the following:

> y3_coef2[1] + coef2[2]*x + coef2[3]*x^2

We just plug in the coefficients and the appropriate x-values and we're done. There's another way that doesn't involve so much typing, especially when dealing with higher order polynomials. It involves matrix multiplication so hopefully you remember something about linear algebra. We're going to create a matrix of x-values and then multiply that by our coefficient vector.

> y3_cbind(1,x,x^2) %*% coef2

The function cbind() is used to bind columns together to form a matrix. The operator %*% is used for matrix multiplication. There is also a function called rbind() that binds vectors together as rows instead of columns. Add the line to the graph:

> lines(x,y3)

And it should look like this:

The last thing we're going to do is increase the polynomial order to the maximum. There are seven data points so the maximum order is six. (Why?) First we fit the new regression, then transform the coefficients, generate new y-values, and add the new line to our graph.

> polyfit6_lm(vel~poly(conc,6),data=df) > coef6_poly.transform(poly(df$conc,6),coef(polyfit6)) > y4_cbind(1,x,x^2,x^3,x^4,x^5,x^6) %*% coef6 > lines(x,y4,lty=2)

When you add the lines, you will see several warnings go by because some of the resulting y-values greatly exceed the range of the graph. The plot should now look like this:

The last thing we're going to look at is what really makes S-PLUS powerful: writing your own functions. If you do something once it's not a problem to type in each command as we've done here. But if you frequently do some particular task, it would be nice to automate it so that it could be done with a single command. In calculating the y-values for the polynomial lines for our graph, we used the cbind() function with a term for each column. If you were going to do this a lot with large order polynomials, you might want to automate the task of creating this matrix. From the file menu, select New… or click on the New toolbar button. Then select Script File from the choices and click OK. This opens a script window where you will write your function. Here is a listing of one way to write this function:

polymat_function(x, deg) { out_numeric(0) for(i in 0:deg) out_cbind(out, x^i) return(out) }

Type this into the script window exactly. Then with the script window active, press F10 or click the Run toolbar button. This will source the script you just wrote and create the function. There is a split screen at the bottom of the script window where you will see the lines of the function being executed as it's sourced. Check for any errors here. If there are errors, it's most likely that you typed something in wrong. Double check the syntax and try again. When it sources OK, use the function objects() to see that your function polymat() was really created. Here's how the function works. Give it a vector for x and the polynomial degree you want. It will return a matrix with all ones in the the first column, the second column equal to x, the third column equal to x², etc. Try it out with a simple vector:

> polymat(1:5, 3) [,1] [,2] [,3] [,4] [1,] 1 1 1 1 [2,] 1 2 4 8 [3,] 1 3 9 27 [4,] 1 4 16 64 [5,] 1 5 25 125

Note the 1:5 used here. This is just a shorthand way of getting a vector of integers between two numbers. Now you can use this function when you want to create a vector of y-values to plot. This is how you would use it as an alternative to how you created y4 previously:

> y4_polymat(x, 6) %*% coef6

Now if you really want to automate the process, let's write a function that will fit whatever order polynomial you want, transform the coefficients, and plot the graph all with one command. Open another script window and type in this function.

poly.fit_function(deg=2, data=df) { pfit_lm(data$vel~poly(data$conc,deg)) pcoef_poly.transform(poly(data$conc,deg), coef(pfit)) x_seq(0, 0.4, length=100) y_polymat(x, deg) %*% pcoef plot(data$conc, data$vel, xlim=c(0,0.4), ylim=c(0,5), xlab="Substrate Concentration", ylab="Reaction Rate") lines(x,y) invisible() }

You should recognize most of the lines in this function The invisible() function is used because nothing is being returned from poly.fit() other than a graph. Notice also that we're calling the polymat() function that we just created. You don't have to put everything in a single function. It's usually more efficient to build a library of smaller functions that do specific tasks and then pull them together as needed. Notice that the data frame is an argument to this function. This way the function is not specific to just one data set.

Remember to source the function (F10 or the Run button). Notice that both arguments have defaults so you could use the function with no arguments within the parentheses. To try a 4^th order polynomial and see the fit, use the function like this:

> poly.fit(4)

That's all there is to it. Hopefully you can see the value of being able to write your own functions this easily. Can you modify the function so that it returns the transformed coefficients?

1. search() - Use this to view your current search path. Your current working directory will likely be in position 1, meaning that's the first place S-PLUS will look for objects you're using and the location where it will save new objects.
2. attach() - Use this to attach new directories or data frames to your search path. If you want to switch working directories, use this function to attach the directory in position 1.
3. detach() - Use this to remove directories or data frames from your search path.
4. par() - Use this to set up various graphics parameters. You can change how graphs appear, put multiple graphs on a single page, etc.