Using Fractional Polynomials for Logistic Regression Modelling in R

I am learning logistic regression modeling using the book "Applied Logistic Regression" by Hosmer.

In chpaters, he suggested using Fractional Polynomials for fitting continuous variable which does not seems to be related to logit in linear fashion. I tried the mfp package and can give exactly the same verbose as the book.

But I don't know how to write the transformed variable based on the output of fractional polynomials. The book only shows example of the transformed variable when $J=2$ with $p_1=0$ and $p_2=-0.5$ (page 101) and when $J=2$ with $p_1=2$ and $p_2=2$ (page 101), But what about the others? Currently my case is $J=2$ with $p_1=-1$ and $p_2=-1$.

I know little about fractional polynomials and the book seems not giving sufficient hits on this part. Can anyone refer me to some place which I can know how to write the polynomial? Thanks.

-

The exposition is obscure but the examples and the discussion on p. 101 make the intentions clear.

Recall that the objective (for the situation with a single continuous covariate $x$) is to generalize logistic regression from the case

$$\text{logit}(y) = \beta_0 + \beta_1 x$$

to a relatively simple nonlinear expression of the form

$$\text{logit}(y) = \beta_0 + \beta_1 F_1(x) + \cdots + \beta_J F_J(x).$$

"Fractional polynomials" [sic] are expressions of the form

$$F(x) = x^p (\log(x))^q$$

for suitably chosen powers $p$ and $q$, with $q$ a natural number and $p$ a real number close to $1$. It is intended that if a high power $q$ of the logarithm is included, then all lower powers $q-1, q-2, \ldots, 1, 0$ will also be included. To be practicable and interpretable, H&L suggest restricting the values of $p$ to the set $P$ = $\{-2, -1, -1/2, 0, 1/2, 1, 2, 3\}$ ($0$ corresponds to $\log$, as usual) and $q$ to the set $\{0,1\}$.

When we limit the fractional polynomial to just two terms ($J=2$), the only possibilities according to these rules are of the form

$$F_1(x) = x^{p_1}, F_2(x) = x^{p_2}$$

for $p_1 \ne p_2$ or

$$F_1(x) = x^p, F_2(x) = x^p\log(x).$$

(The case $p=0$ corresponds to using $F_1(x) = \log(x)$ and $F_2(x) = (\log(x))^2$.)

These possibilities can be uniquely determined by a non-decreasing sequence of $J=2$ elements of $P$. The sequence $(p_1,p_2)$ with $p_2 \gt p_1$ specifies the first kind of fractional polynomial and the sequence $(p_1,p_2) = (p,p)$ specifies the second kind. Because $P$ has eight elements, this gives $\binom{8+1}{2} = 36$ possibilities for $J=2$. For instance, your case of $(-1,-1)$ specifies the model

$$\text{logit}(y) = \beta_0 + \beta_1 \frac{1}{x} + \beta_2 \frac{\log(x)}{x}.$$

(H&L go on to recount an approximate procedure in which partial likelihood ratio tests are used to fit the best model with $J=1$ (there are just eight of these) and then the best model with $J=2$ is fit. Each contributes approximately $2J$ degrees of freedom in the resulting chi-squared test.)

Of course, to be really sure of what R is doing, you should either look at the source code, or fit the model and plot the predictions against the data, or both.

-