I can't decide what is the best way to assess the stability of a higher order fractional polynomial model. To use an example I have been working on, I am analyzing a dataset with panel data selected from 5 different universities. I am modelling the non-linear changes of bone mass during lifetime. Using F-test and log-likelihood tests the results suggest that the best-fitting model consists of three powers,
Age^2 and a linear component of
Age. I am then fitting the selected model in Stata:
fp <Age>, center fp(3 2 1) replace: reg Bone_Mass <Age> i.Sex i.Uni1 i.Uni2 i.Uni3 /* */ i.Uni4 i.Uni5 estimates store eq_1 predict fit1
Uni are entered as dummy variables in the model. Because this identified model is more complex than regression models with one or two powers I wanted to some further verification that it is actually a stable model. Research has extensively described the
mfpboot command to resample residuals and make sure about variable inclusion in multivariate fractional polynomial models (e.g. here and here). However, I am actually fitting a univariate model.
My first idea was to leave a cohort out and try to replicate the model selected. However, if all cohorts have a significant effect on the best-fitting model why would I be expecting to replicate the solution? Another idea would be to actually bootstrap the residuals for the identified regression using the separate powers as predictors and see whether they all remain significant in this model. Something that looks like this for example:
gen Age3=Age^3 gen Age2=Age^2 bootstrap, reps(1000): reg Bone_Mass Age*_b[Age_1] Age2_b*[Age_2] Age3*_b[Age_3]/* */i.Sex i.Uni1 i.Uni2 i.Uni3 i.Uni4 i.Uni5
_b[x] represent the regression coefficients from the final fitted fractional polynomial model
I don't, however, see how multicollinearity would not cause issue in the interpretation of the standard errors of these predictors.