[R] Problems with the commands FUNCTION and DERIV to build a polynomial

[Thomas Lumley]

On Thu, 17 Sep 2009, [ISO-8859-1] Noela Sánchez wrote:

> Hi all,
>
> I need to automate a process in order to prepare a  a big loop in the future
> but I have a problem with the *command function*
>
> First I fit a model with lm
>
>>
> model1<-lm(data2[,2]~data2[,1]+I(data2[,1]^2)+I(data2[,1]^3)+I(data2[,1]^4))
>
> I extract the coefficients to build the polynomial.
>
> coef<-as.matrix(model1$coefficients)
>
> In the next step I need to define the polynomial to derive it. If I write
> the coefficients manually (writing the numbers by hand) the deriv command
> works fine!
>
> bb<-deriv(~2847.22015 -463.06063*x+ 25.43829*x^2 -0.17896*x^3, namevec="x",
> function.arg=40)
>
> I would like to automate this step by being able to extract the coefficients
> from the linear model and adding them into the polynomial (and not write
> them by hand)!
>
> But if I build the polynomial with the function(x) command calling the *
> coef* values, the numeric values are not interpreted, the command function
> does not read properly the coefficients from the linear model.
>
>>  fun<-function(x) coef[1]+coef[2]*x+coef[3]*x^2+coef[4]*x^3
>
>> fun
>
> function(x) coef[1]+coef[2]*x+coef[3]*x^2+coef[4]*x^3
>
> How can i avoid to write the values of the coefficients by hand??
>

As is often the case there are two parts to this answer: how to do what you asked, and how to do what you want.

You can use bquote() to get the value of coef[1] into the expression, so
   deriv(bquote(.(coef[1])+.(coef[2])*x+.(coef[3])*x^2+.(coef[4])*x^3),"x",
     function.arg="x")
will give the derivative.  Bill Venables would tell you to use Horner's Rule and write the polynomial as
    coef[1]+x*(coef[2]+x*(coef[3]+x*coef[4]))
to get better speed and numerical stability, but the same trick still works.

However, you really don't need deriv() to differentiate a polynomial, and so you can use ordinary lexical scope rather than substitution, for a much tidier answer

derivative <- function(coef){
     function(x) coef[2]+x*(2*coef[3]+x*3*coef[4])
}

         -thomas


Thomas Lumley			Assoc. Professor, Biostatistics
tlumley at u.washington.edu	University of Washington, Seattle