My goal is to write a function in R that accepts coefficients for a fractional polynomial (FP) and returns a vectorized function which evaluates the specified FP for given input numbers. The FP definition has two important rules:
- x^0 is defined as log(x)
- powers can have multiple coefficients, where the 2nd coefficient for power p adds a factor of log(x) to the additive term (x^p*log(x)), the 3rd adds a factor of log(x)^2 (x^p*log(x)^2), and so on
My current solution below builds the FP-function as a string, parses the string and returns a function which evaluates the expression. My question is if there is a better/faster way that avoids eval(parse())
- possibly using some substitute()
magic.
The function must deal with having the number of coefficients per power not known in advance, but specified when being called. The final FP evaluation needs to be fast as it is called very often.
It would be nice not to be limited to the standard powers -2, -1, -0.5, 0, 0.5, 1, 2, 3. Ideally, the desired function would do two steps at once: accept FP-coefficients as well as a vector of numbers and return the FP-values for the input while still being fast.
getFP <- function(p_2, p_1, p_0.5, p0, p0.5, p1, p2, p3, ...) {
p <- as.list(match.call(expand.dots=TRUE)[-1]) # all args
names(p) <- sub("^p", "", names(p)) # strip "p" from arg names
names(p) <- sub("_", "-", names(p)) # replace _ by - in arg names
## for one power and the i-th coefficient: build string
getCoefStr <- function(i, pow, coefs) {
powBT <- ifelse(as.numeric(pow), paste0("x^(", pow, ")"), "log(x)")
logFac <- ifelse(i-1, paste0("*log(x)^", i-1), "")
paste0("(", coefs[i], ")*", powBT, logFac)
}
onePwrStr <- function(pow, p) { # for one power: build string for all coefs
coefs <- eval(p[[pow]])
pwrStr <- sapply(seq(along=coefs), getCoefStr, pow, coefs)
paste(pwrStr, collapse=" + ")
}
allPwrs <- sapply(names(p), onePwrStr, p) # for each power: build string
fpExpr <- parse(text=paste(allPwrs, collapse=" + "))
function(x) { eval(fpExpr) }
}
An example would be -1.5*x^(-1) - 14*log(x) - 13*x^(0.5) + 6*x^0.5*log(x) + 1*x^3
which has specified powers (-1, 0, 0.5, 0.5, 3) with coefficients (-1.5, -14, -13, 6, 1).
> fp <- getFP(p_1=-1.5, p0=-14, p0.5=c(-13, 6), p3=1)
> fp(1:3)
[1] -13.50000000 -14.95728798 0.01988127