R Language – How to Approach Finding the Inverse of a Function?

Question

I'm reading a paper by Steve Horvath (2013) for the epigenetic clock he developed. He used an age calibration function before carrying out an elastic net regression analysis. The function F for transforming age is as follows (adult age = 20):

F(age)=log(age+1)-log(adult.age+1) if age<= adult.age

F(age)=(age-adult.age)/(adult.age+1) if age> adult.age

Using R, I have written the function as follows:

ageC <- function(age) {
  if (age <= 20) {
    ageC = log(age+1)-log(21)}
  if (age > 20) {
    ageC = (age-20)/(21)
  }
  return(ageC)
}

My question is how do I find the inverse of the function F? The model predicts DNA methylation age (DNAmAge) using the following:

DNAmAge = inverse.F(b0 + b1CpG1 + . . . + b353CpG353)

What is the "inverse.F" portion?

Answer 1

For brevity, let $x_0$ be the adult-age and let $x$ be the age (I will assume that the age is non-negative and adult-age is positive). Your function of interest $F: \mathbb{R}_{0+} \rightarrow \mathbb{R}$ is given by:

$$\begin{align} F(x) &= \int \limits_{x_0}^x \frac{ds}{\min(s,x_0)+1} \\[6pt] &= \begin{cases} \log(x+1) - \log(x_0+1) & & \text{for } 0 \leqslant x \leqslant x_0, \\[6pt] (x-x_0)/(x_0+1) & & \text{for } x > x_0. \\[6pt] \end{cases} \end{align}$$

This function has first derivative $F'(x) = 1/(\min(x,x_0)+1)$ so it is strictly increasing and has an inverse. (It is also simple to show that this is a continuous function, so its inverse is also continuous.) For this part I will use $r=F(x)$ to represent the age-calibration. Solving for $x$ using some simple algebra on each part gives the corresponding inverse function $F^{-1}: \mathbb{R} \rightarrow \mathbb{R}_{0+}$, which is:

$$F^{-1}(r) = \begin{cases} (x_0+1) \exp(r) - 1 & & \text{for } r \leqslant 0, \\[6pt] x_0 + r(x_0+1) & & \text{for } r > 0. \\[6pt] \end{cases}$$

---

Computation in R: It would be preferable to program these functions without hard-coding the adult-age you wish to use. If you really want to get fancy then you can also add the derivatives to the functions as gradient and hessian attributes, which is useful if you use these functions later in optimisation work. You can program these functions in R as follows:

ageC <- function(age, age.adult, gradient = FALSE, hessian = FALSE) {
  if (age < 0)          stop('Error: Input age must be non-negative')
  if (age.adult <= 0)   stop('Error: Input age.adult must be positive')
  if (age <= age.adult) { OUT <- log((age+1)/(age.adult+1)) }
  if (age > age.adult)  { OUT <- (age-age.adult)/(age.adult+1) }
  if (gradient) {
    attr(OUT, 'gradient') <- 1/(min(age, age.adult)+1) }
  if (hessian)  {
    attr(OUT, 'hessian')  <- ifelse(age < age.adult, -1/(age+1)^2, 0) }
  OUT }

ageC.inverse <- function(r, age.adult, gradient = FALSE, hessian = FALSE) {
  if (age.adult <= 0)   stop('Error: Input age.adult must be positive')
  if (r <= 0) { OUT <- (age.adult+1)*exp(r) - 1 }
  if (r > 0)  { OUT <- age.adult + r*(age.adult+1) }
  if (gradient) {
    attr(OUT, 'gradient') <- ifelse(r < 0, (age.adult+1)*exp(r), age.adult) }
  if (hessian)  {
    attr(OUT, 'hessian')  <- ifelse(r < 0, (age.adult+1)*exp(r), 0) }
  OUT }