Calculate derivatives & curvature of a polynomial

Question

I have a data frame which records the skin temperature of a bunch of folks over time. I would like to:

1. Fit a quadratic polynomial to each id's SkinTemp over Time;
2. Calculate the curvature.

This seems much harder than it should be.

I asked about the first part in Fitting a quadratic curve for each data set that has different lengths, but I can't move forward to calculate derivatives and curvature.

df <- data.frame(Time = seq(65),
                 SkinTemp = rnorm(65, 37, 0.5),
                 id = rep(1:10, c(5,4,10,6,7,8,9,8,4,4)))

#Predict data points for each quadratic
fitted_models = df %>% group_by(id) %>% do(model =
                lm(SkinTemp ~ Time+I(Time^2), data = .))

Now I need to calculate the curvature k = y''/(1 + y' ^ 2) ^ (3 / 2), where y' and y'' are 1st and 2nd derivative of y with respect to x.

I thought I could ask the predict function to give me derivatives by passing for example deriv = 2, but it doesn't seem to work.

predQ <- lapply(unique(df$id),
                function(x) predict(deriv = 2,fitted_models$model[[x]])) 

So I amended this function, which seems to work OK but isn't there a built-in function for this task?

deriv <- function(x, y) diff(y) / diff(x)
middle_pts <- function(x) x[-1] - diff(x) / 2
second_d <- lapply(unique(df$id),
                   function(x) deriv(middle_pts(df$Time[df["id"]==x]), deriv(df$Time[df["id"]==x], df$SkinTemp[df["id"]==x])))

Answer 1

You are just hindered by derivative computations of a polynomial. How about using function g defined in my this answer?

For example, your first model has polynomial coefficients:

pc <- coef(fitted_models[[2]][[1]])
#(Intercept)        Time   I(Time^2) 
#38.36702120 -0.61025716  0.04703084 

Let's say you just want to evaluate derivatives and curvatures at observed locations:

x <- with(df, Time[id == 1])
#[1] 1 2 3 4 5

Then you can do analytical computation step by step:

## 1st derivative
d1 <- g(x, pc, 1)
#[1] -0.5161955 -0.4221338 -0.3280721 -0.2340104 -0.1399487

## 2nd derivative
d2 <- g(x, pc, 2)
#[1] 0.09406168 0.09406168 0.09406168 0.09406168 0.09406168

## curvature: d2 / (1 + d1 * d1) ^ (3 / 2)
d2 / (1 + d1 * d1) ^ (3 / 2)
#[1] 0.06599738 0.07355055 0.08069004 0.08683238 0.09136444

Isn't this much better than your finite differencing approximation?

---

Note that g can also evaluate nderiv = 0L, i.e., the polynomial itself:

g(x, pc, 0)
#[1] 37.80379 37.33463 36.95953 36.67849 36.49151

which agrees with predict.lm:

predict.lm(fitted_models[[2]][[1]], data.frame(Time = x))
#       1        2        3        4        5 
#37.80379 37.33463 36.95953 36.67849 36.49151 

The function g judges the degree of the polynomial by the length of polynomial coefficient vector pc. A length-3 vector means degree = 2. It is designed for raw polynomials, not orthogonal ones.

---

computations for all groups

To do the curvature computation for all groups, I would use Map.

polynom_curvature <- function (x, pc) {
  d1 <- g(x, pc, 1L)
  d2 <- g(x, pc, 2L)
  d2 / (1 + d1 * d1) ^ (3 / 2)
  }

pc_lst <- lapply(fitted_models[[2]], coef)
Time_lst <- split(df$Time, df$id)
result <- Map(polynom_curvature, Time_lst, pc_lst)
str(result)
#List of 10
# $ 1 : num [1:5] 0.066 0.0736 0.0807 0.0868 0.0914
# $ 2 : num [1:4] -0.106 -0.12 -0.131 -0.135
# $ 3 : num [1:10] 0.0795 0.0897 0.0988 0.1058 0.1095 ...
# $ 4 : num [1:6] -0.098 -0.107 -0.113 -0.115 -0.112 ...
# $ 5 : num [1:7] -0.0878 -0.0923 -0.0946 -0.0944 -0.0917 ...
# $ 6 : num [1:8] 0.0752 0.0811 0.0857 0.0886 0.0895 ...
# $ 7 : num [1:9] 0.0397 0.0405 0.0411 0.0414 0.0416 ...
# $ 8 : num [1:8] 0.0178 0.018 0.0182 0.0184 0.0185 ...
# $ 9 : num [1:4] -0.151 -0.161 -0.159 -0.146
# $ 10: num [1:4] 0.1186 0.1129 0.1033 0.0917