I have points on a 2D graph. I want to find the best 3rd polynomial that fits this model and get its first derivative. But I can't get D function working. Here is simple example:
a <- 0:10
b <- c(2, 4, 5, 8, 9, 12, 15, 16, 18, 19, 20)
plot(a, b)
m1 <- lm(b ~ a + I(a ^ 2) + I(a ^ 3))
s <- coef(m1)
## try to get 1st derivative of the regression polynomial
D(expression(s[1] + s[2] * a + (a ^ 2) * s[3] + (a ^ 3) * s[4]), "a")
Error in D(expression(s[1] + s[2] * a + (a^2) * s[3] + (a^3) * s[4]), :
Function '
[' is not in the derivatives table
I want to avoid computing numerical derivative by differencing. Thanks for help!
The error message you see "Function '[' is not in the derivatives table" is because D can only recognize a certain set of functions for symbolic operations. You can find them in ?D:
The internal code knows about the arithmetic operators ‘+’, ‘-’,
‘*’, ‘/’ and ‘^’, and the single-variable functions ‘exp’, ‘log’,
‘sin’, ‘cos’, ‘tan’, ‘sinh’, ‘cosh’, ‘sqrt’, ‘pnorm’, ‘dnorm’,
‘asin’, ‘acos’, ‘atan’, ‘gamma’, ‘lgamma’, ‘digamma’ and
‘trigamma’, as well as ‘psigamma’ for one or two arguments (but
derivative only with respect to the first). (Note that only the
standard normal distribution is considered.)
While the "[" is actually a function in R (read ?Extract or ?"[").
To demonstrate the similar behaviour, consider:
s <- function (x) x
D(expression(s(x) + x ^ 2), name = "x")
# Error in D(expression(s(x) + x^2), name = "x") :
# Function 's' is not in the derivatives table
Here, even though s has been defined as a simple function, D can do nothing with it.
Your problem has been solved by my recent answers for Function for derivatives of polynomials of arbitrary order (symbolic method preferred). Three methods are provided in three of my answers, none of which are based on numerical derivatives. I personally prefer to the one using outer (the only answer with LaTeX math formula), as for polynomials everything is exact.
To use that solution, use the function g there, and specify argument x by values where you want to evaluate the derivative (say 0:10), and pc by your polynomial regression coefficients s. By default, nderiv = 0L so the polynomial itself is returned as if predict.lm(m1, newdata = list(a = 0:10)) were called. But once nderiv is specified, you get exact derivatives of your regression curve.
a <- 0:10
b <- c(2, 4, 5, 8, 9, 12, 15, 16, 18, 19, 20)
plot(a, b)
m1 <- lm(b ~ a + I(a ^ 2) + I(a ^ 3))
s <- coef(m1)
#(Intercept) a I(a^2) I(a^3)
# 2.16083916 1.17055167 0.26223776 -0.02020202
## first derivative at your data points
g(0:10, s, nderiv = 1)
# [1] 1.1705517 1.6344211 1.9770785 2.1985237 2.2987568 2.2777778 2.1355866
# [8] 1.8721834 1.4875680 0.9817405 0.3547009
Other remark: It is suggested that you use poly(a, degree = 3, raw = TRUE) rather than I(). They do the same here, but poly is more concise, and make it easier if you want interaction, like in How to write interactions in regressions in R?