# logistic regression multivariable fractional ploynomials stata vs. R

I a going through Hosmer, Lemenshow and Sturdivant's (HLS) Applied Logistic Regression (2013) and trying to interpret the difference between what STATA is doing and what R is doing. Concerning the fit of the model using multivariable fractional ploynomials (MFP), HLS looks at the glow500 study where the dependent variable is fracture and there are both dichotomous and continuous independent variable. You do not need the book to consider these questions:

Here is the data and the R call to the MFP function

library(mfp)
library(data.table)
colnames(mydat) = c("SUB_ID","SITE_ID","PHY_ID","PRIORFRAC","AGE","WEIGHT","HEIGHT","BMI","PREMENO","MOMFRAC","ARMASSIST","SMOKE",
"RATERISK","FRACSCORE", "FRACTURE","V16","V17")
m = mfp(FRACTURE ~ PRIORFRAC + as.factor(RATERISK) + MOMFRAC  + ARMASSIST + SMOKE + PREMENO +
fp(AGE, df=4, select = .15, alpha = .05) +
fp(BMI, df=4, select = .15, alpha = .05) +
fp(WEIGHT, df=4, select = .15, alpha = .05) +
fp(HEIGHT, df=4, select = .15, alpha = .05) ,
family = binomial, data= mydat, verbose = TRUE)
summary(m)


In HLS's book on page 142 there is table that is Cycle 2 of the MFP procedure. Cycle 2 they says contains the "partial liklihood ratio tests that are based, not on the full 11 covariate model, but on a 5 covariate model (from Cycle 1)." NOTE: HERE HE SAYS THERE ARE 5 VARIABLES that are significant from cycle 1.

Here is the MFP function call in R using the same data he used in STATA.

I have a couple questions.

QUESTION 1

On page 141 in the book BMI is not significant so it is not included in cycle 2. You can see the same non-significance in the R output. If you look at Cycle 1 you see

           BMI
506.098
503.837     1
501.743     0.5
501.574     -2 1


So the liklihood ratio test is

1-pchisq(506.098 - 501.574,4) = p-value of .34 > .15 so DO NOT INCLUDE in cycle 2


This is the Deviance of the NULL Model excluding BMI compared to the deviance of the best 2-term polynomial fit powers (-2,1) the p-value = .34 which is > .15 so BMI SHOULD NOT BE INCLUDED in CYCLE 2. But it is. You can look under Cycle 2 of the output and you see:

                 BMI
506.674
505.865     1
504.843     -2
504.224     -2 -2


So now in R - in Cycle 2 you see 6 (NOT 5 like he says in the book) variables in the model in cycle 2. BMI again is not significant

1-pchisq(506.674 - 504.224,4) =  p-value = .65 so again not significant in cycle 2 comparing the NULL model in cycle 2 vs. the best 2 term polynomial (-2,2)


So my first question is why is BMI included in Cycle 2 in R? Looking at BMI in cycle 1 it is NOT significant.

QUESTION 2 Under the Fractional Polynomials part of the output I see:

Fractional polynomials
df.initial select alpha df.final power1 power2
PRIORFRAC                     1   1.00  0.05        1      1      .
AGE                           4   0.15  0.05        1      1      .
as.factor(RATERISK)2          1   1.00  0.05        1      1      .
as.factor(RATERISK)3          1   1.00  0.05        1      1      .
MOMFRAC                       1   1.00  0.05        1      1      .
BMI                           4   0.15  0.05        0      .      .
WEIGHT                        4   0.15  0.05        0      .      .
ARMASSIST                     1   1.00  0.05        1      1      .
HEIGHT                        4   0.15  0.05        1      1      .
SMOKE                         1   1.00  0.05        1      1      .
PREMENO                       1   1.00  0.05        1      1      .


In the documentation (https://cran.r-project.org/web/packages/mfp/vignettes/mfp_vignette.pdf) df.final is defined as

"a vector containing the degrees of freedom of each covariate at convergence of the backfitting  algorithm (m=4 for second degree FP, m=2 for first degree FP, m=1 for untransformed variable, m=0 if covariate was excluded)."


can you explain how the df changes for weight and BMI so that they are excluded from the final model?

• I'm voting to close this question as off-topic because I think it better fits to stats.stackexchange (for statistics) or stack.overflow (as it is about software).
– Ric
Jul 6 '16 at 8:14