The normal accumulation of bone during childhood may play an important role in avoiding or delaying osteoporosis later in life. Inadequate bone accretion during childhood can be related to lifestyle factors, such as diet and physical activity [2], chronic medical conditions with primary or secondary effects on bone [3], and concomitant medications. Thus, identifying low bone mineral density (BMD) in children may allow clinicians to make more informed decisions about treatments for children with poor bone mineral accretion.

Dual-energy X-ray absorptiometry (DXA) is the most commonly used method of measuring BMD because of its reproducibility, safety, and widespread availability, and it is recommended for clinical assessment in children [4]. However, clinical assessment of BMD in children requires special consideration of the expected patterns of change in BMD associated with growth and development [5].

The Bone Mineral Density in Childhood Study (BMDCS) is a prospective, longitudinal study with the goal of providing the necessary reference data for the clinical assessment of bone density in children. The publication of the first results of the BMDCS [6] used the LMS [7] statistical method to construct sex- and race-specific reference percentiles of BMD and bone mineral content (BMC) relative to age.

The LMS method smoothes the estimated distributions across age and adjusts the data for skewness, a common problem in biological data. The primary disadvantage of the LMS method is the inability to easily include multiple effects within the same model.

Several major classifying variables, including age, gender, and race, predict BMD. Additionally, the anthropometric variables of weight, height, and body composition (amount of lean and fat tissue in the body) correlate well with BMD [8, 9]. Several of these factors are highly collinear. As children age, their height, weight, lean body mass, and percent body fat tend to change concurrently.

Others have used general linear regression to include these anthropometric variables in their models [9–12]. Unlike the LMS method that groups subjects by age, these methods use age as a continuous variable, but like the LMS method they do not include interaction terms that would allow the model to consider the non-linear changing effects of anthropometric variables as subjects age. For instance, the effect of height on BMD may change as subjects grow older. Including interaction terms is particularly important when the observed population spans several developmental milestones. Furthermore, these methods often transform BMD by using the inverse or log functions, making interpretation of the parameter effects and the model results difficult.

Our analysis presents a strategy of modeling normal childhood BMD data by taking into account potentially relevant anthropometric variables and including the effects of the interaction of these variables with age while preserving the smoothness of the estimates as is provided by the LMS method.

Different combinations of explanatory variables resulted in different models (Table 2). The traditional model, model A, approximates the previous analysis [6], but instead of using medians and a power transformation as the LMS method does, model A uses the more standard parametric technique employing means, no power transformation and a symmetric error distribution. We find that the previous results via the LMS model and model A are almost identical. We will use model A as an analog to the output of the LMS method for comparison purposes.

The full model, model B, is the “kitchen sink” model, in which everything is included, even parameters that may not add much information to the model.

Accurately evaluating sexual maturity in practice requires careful cross-training of the assessors and may be uncomfortable for both the subject and the assessor. Measuring the percent body fat means another DXA measurement, slightly increasing radiation dose to the subject. The other parameters can be gathered quickly and agreeably, resulting in the practical model, model C.

Because of computational difficulty with orthogonalizing the nominal variable of maturity, fully smoothed models for some parameter sets were not calculated. The results for the unsmoothed model provide us with a reasonable estimate of the effectiveness of the final model. As seen in Tables 1 and ​and2,2, dropping sexual maturity as a predictor does not significantly reduce the predictive power of the model.

Statistically, model F may be optimal in terms of creating the most powerful model while adding the fewest predicting parameters, but Tukey's confidence intervals indicate that, at the black/non-black level, race still matters for some age/gender combinations. Including both weight and percent body fat will intrinsically account for lean mass in the model [9]. Model D may be optimal in terms of balancing the needs of the subjects with the most powerful model, and this model D will be investigated further.

When segmented by age, gender, and race, there are 64 simultaneous fits for each model. For the sparse model, 13 of these groups show a lack of normality in the underlying BMD values at the 0.05 level using a Shapiro-Wilk's test for normality. This lack of normality was the reason previous investigators chose to use the LMS model. For the smoothed model D, the residuals of only eight of these fits do not show normality due to outliers that make the distributions look leptokurtic. Since the outliers tend to appear in an uneven manner, the leptokurtosis leads to perceived skewness, and all eight of these groups have a moment of skewness twice the error of skewness, meaning the residuals show more than mild skewness [18]; however, selectively removing one outlier from the extreme tail of the sample allows all but one of these groups to pass a test of normality. In samples of this size, this result reinforces the notion that the underlying distribution is actually less skewed and more normal.

By examining the smoothed coefficients of each parameter, we may glean some information about the parameter's effect and how this effect changes with age. Weight is positively related to BMD (Fig. 2a), but the influence declines with increasing age. For most groups, greater percent body fat is associated with lower bone density (Fig. 2b). Since the parameters are adjusted for weight, this negative association between percent body fat and BMD indicates that for two children of the same sex, race, age, height, and weight, the child with a higher percent body fat has a lower lean body mass and lower BMD, reflecting the known strong association of lean body mass and BMD [9]. Also note that for late teen black females, the parameter reverses sign, but this is the area where data are most sparse.

Whereas weight usually has a strong positive contribution to BMD, height produces a much smaller negative correction after weight, age, and percent body fat are accounted for (Fig. 2c). This can be interpreted as an effect of the given bone mass, defined largely by body weight for a child with given sex, race, age, and percent body fat being distributed over a larger projected area in a taller child and, thus, resulting in a lower measured BMD. For children at older ages, both black males and females show a reversal of this pattern, and height contributes positively to BMD.

The intercepts of each group represent a baseline before the other effects are layered in (Fig. 2d). Again, a difference between the black and non-black groups is striking.

Building a more complete model, including anthropometric factors, should allow for narrower distributions than would be allowed by simply using age, gender, and race. We can indeed confirm that the expected distribution narrows as parameters are added; however, after age, gender, race, weight, height, and percent body fat are included in the model, sexual maturity does not contribute in narrowing the distribution further.

In addition to narrowing the estimated distributions, the inclusion of anthropometric variables may produce shifts in the expected BMD. Subject weight creates a major shift in the mean, whereas subject height has little influence (Fig. 3a). The sparse model (model A) shows a wider distribution as it includes all variations in weight and height. Gender and race produce additional shifts in the model curves (Fig. 3b).

When comparing the total root mean squared error (RMSE) of the sparse model, model A, to the total RMSE produced by a more complex model, model D, a lower total RMSE is observed for the more complex model (Table 2). Furthermore, plotting the race/gender subgroups' individual RMSEs by age (Fig. 4) shows that model D also has a consistently lower RMSE at the group level. Thus, a model including relevant anthropometric data does produce narrower distributions. Also note that, since the RMSEs are expressed in the same scale as the dependent variable, they tend to increase as the subjects grow older just as the subjects' BMD levels increase with age.

Comparing model D, the full model, and model A, the reduced model, the F statistic returns a value above 14, where the degrees of freedom in the denominator and the numerator are both above 7000. The p value below 0.001 means that adding the anthropometric parameters makes numeric sense and helps to conclude that model D is superior to model A.

Not factoring in the effect of weight as a determinant of BMD leads to the possibility that children of normal weight, based on the 50% value of the CDC growth chart, will be judged as having low BMD when that is not the case (Fig. 5a). The sparse model agrees better with individuals who are at the 75th percentile of weight and height by CDC standards (Fig. 5b) than it does with those at the 50th percentile of the CDC standards.

Manufacturers have released average BMD values starting at the age of 20 years. The method described here might be improved by including such values at the end point of the models as fixed or target values; however, without the relevant anthropometric data available, it is not clear how to account for these values.

The fact that Hologic fan-beam instruments have been used for this study may have an influence on the resulting data. Change in body weight may change the distance of the spine above the scan table, and this, in turn, would change the magnification of the bones in the image. Such magnification errors would clearly influence a geometric parameter like bone area as well as bone mineral content, which is derived by multiplying BMD with bone area. BMD is least influenced by magnification and has been shown to be relatively independent from the distance of the bone above the table [20].

The small sample size of late teen black females and the difference between the behavior of the effect coefficients that group produces, compared to the rest of the data set, calls those results into question. Combining the data and eliminating race seems attractive, but Tukey's confidence intervals for the residuals of the age/gender/weight/height/fat fits without race still indicate that race is a significant determinant for BMD.

The models presented differ from other models that take anthropometrics into account for the interpretation of DXA data. A number of authors have suggested a multi-step approach [9, 12, 21], which narrows the appropriate comparisons down by gradually including additional anthropometric variables. Our models include these variables up front, simplifying the interpretation step of the DXA data. These models, however, critically rely on the implemented smoothing approach, as the subgroups become very small even in a data set containing several thousand observations.

Pubertal stage did not have a major influence on our models. Horlick et al. made similar observations for their whole body BMC models [12], and they argued that the use of a normal population, the cross-sectional design of the study and the consideration of anthropometric parameters diminished the importance of pubertal stage. This would, however, not mean that abnormal pubertal development would not have an influence on bone density.

Whereas there are alternative approaches to correcting BMD measurements that account for skeletal size, most notably bone mineral apparent density [22], there does not appear to be clear agreement about the best way to correct for such an effect or agreement about how meaningful such a correction would be. When and if such a correction is agreed upon, the methods presented in this paper would also apply to this skeletal size-corrected parameter. However, because several anthropometric parameters were considered in the models, it is quite possible that the introduction of a new skeletal size-corrected parameter is unnecessary.

A purely statistical approach to the problem of classifying juvenile subjects with respect to BMD would be to discard age entirely as an explanatory variable. The colinearity of both weight and height with age would make such an approach seem viable. Indeed, fitting BMD by gender and race against weight, height, and percent body fat produces an impressive adjusted R² of 0.81 and an RMSE of 0.094, although not as good as model D's values of 0.83 and 0.088, respectively. By its nature, such a model would not include any interaction effects. Unfortunately, removing the percent body fat from this equation drives the R² to 0.77 and increases the RMSE to 0.103, values that are comparable to what is achieved by using just age, gender and race, the sparse model A.

In conclusion, the proposed model adjusts for height, weight, and percent body fat. A similar model using the same type of smoothing can be implemented without percent body fat if that measurement is not available. These models adjust for the differences in the distribution of body weight between this sample and the CDC growth charts.

The models presented here can be used to determine a child's Z-score adjusted not just for age, gender, and race but also for body weight, size, and composition. These Z-scores are reasonably smooth at age breaks because of the MCS² smoothing. Individuals are compared to more appropriate means, considering their anthropometric status, and placed within a narrower distribution, allowing for more sensitive clinical classification.