Biostatistics Advance Access originally published online on June 16, 2007
Biostatistics 2008 9(1):159-171; doi:10.1093/biostatistics/kxm020
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Nonlinear growth generates age changes in the moments of the frequency distribution: the example of height in puberty
Centre for Paediatric Epidemiology and Biostatistics, UCL Institute of Child Health, 30 Guilford Street, London WC1N 1EH, UK tim.cole{at}ich.ucl.ac.uk
Department of Social Medicine, University of Bristol, Bristol, UK
Statistical Laboratory, University of Cambridge, Cambridge, UK
Centre for Paediatric Epidemiology and Biostatistics, UCL Institute of Child Health, 30 Guilford Street, London WC1N 1EH, UK
* To whom correspondence should be addressed.
| SUMMARY |
|---|
|
|
|---|
Higher moments of the frequency distribution of child height and weight change with age, particularly during puberty, though why is not known. Our aims were to confirm that height skewness and kurtosis change with age during puberty, to devise a model to explain why, and to test the model by analyzing the data longitudinally. Heights of 3245 Christ's Hospital School boys born during 1927–1956 were measured twice termly from 9 to 20 years (
). Treating the data as independent, the mean, standard deviation (SD), skewness, and kurtosis were calculated in 40 age groups and plotted as functions of age t. The data were also analyzed longitudinally using the nonlinear random-effects growth model
, with
the cross-sectional data,
the individual mean curve, and
and
subject-specific random effects reflecting variability in age and height at peak height velocity (PHV). Mean height increased monotonically with age, while the SD, skewness, and kurtosis changed cyclically with, respectively, 1, 2, and 3 turning points. Surprisingly, their age curves corresponded closely in shape to the first, second, and third derivatives of the mean height curve. The growth model expanded as a Taylor series in
predicted such a pattern, and the longitudinal analysis showed that adjusting for age at PHV on a multiplicative scale largely removed the trends in the higher moments. A nonlinear growth process where subjects grow at different rates, such as in puberty, generates cyclical changes in the higher moments of the frequency distribution.
Keywords: Growth; Height; Kurtosis; Random effects; Skewness
| 1. INTRODUCTION |
|---|
|
|
|---|
The construction of growth references for height or weight involves modeling the frequency distribution as a function of age, using techniques such as the
-
-
(LMS) method (Cole and Green, 1992
Puberty is a time when these distributional changes are most marked. The mean cross-sectional height curve is steeper in puberty than before or after, reflecting the greater mean annual height increment, and the SD is increased due to the variation in growth tempo. The distribution also becomes transiently skew, as Henry Bowditch (1891)
showed over 100 years ago. Franz Boas (1892a
, b
) provided a mathematical explanation for this, but "couched in obscure prose and almost impenetrable algebra" (Tanner, 1981
), it has been little cited.
The pubertal height growth spurt in individual children is dramatic, the height increment amounting to 30% of adult height in boys (Tanner, 1978
). Peak height velocity (PHV) can reach 10 cm/year at a mean age of around 14 years, though the timing of the spurt varies between individuals reflecting differences in their rate of maturation or "tempo of growth" (Tanner, 1962
). The population distribution of the growth spurt is summarized by the mean and SD of PHV and age at PHV. Figure 1 illustrates height and height velocity for 3 typical boys with early, average, and late age at PHV; PHV and age at PHV are seen to be inversely correlated (Tanner, 1962
).
|
The variability in age at PHV is usually assumed to operate on an additive scale, yet the underlying variability is in the rate of maturation, implying a multiplicative age scale (Beath, 2007). Which is the more appropriate scale is one of the questions the paper addresses.
Also of interest is the question whether, and how well, cross-sectional data can estimate the parameters of PHV usually estimated from longitudinal data. Margaret Merrell (1931)
showed formally that during puberty, the population mean height curve is less steep than the height curve as seen in the average individual, due to differences in growth tempo. This "Merrell bias" in the population curve underestimates PHV and, coupled with an absence of information about variation in the age at PHV, it has led to the view that longitudinal information is essential to quantify the growth spurt.
The motivation for the current study was access to a uniquely large data set of heights in boys from Christ's Hospital (CH) School followed longitudinally through puberty. The main aims of the study were 3-fold: (i) to quantify changes with age in the first 4 moments of the cross-sectional height distribution, (ii) to develop a mathematical model to explain why age-related changes in the moments are the same shape as derivatives of the mean height curve, and (iii) to reanalyze the data longitudinally to test the assumptions underlying the model.
The aims taken together tell a story, so the statistical methods have been split into 3 parts and are presented with their corresponding results (Sections 3–5).
| 2. SUBJECTS AND DATA |
|---|
|
|
|---|
The school CH was established in the 16
century as a poor school for orphans and beggars, and by the 20
century it accepted boys from a wide variety of social backgrounds. The school has a long tradition of measuring the heights and weights of its pupils. Dr Gerald Edward Friend (1935)
The data consist of serial heights in 3245 boys born during 1927–1956 who attended CH during 1936–1969. The boys were measured twice a term between the ages of 9 and 20 years, a mean of 40 heights per boy (median 42, range 1–63) and 129 700 measurements altogether. Other recent papers based on data from the same study have focussed on the 3175 boys who were subsequently traced via the National Health Service Central Register (Sandhu and others, 2006a
,b). For more information about the data set, see the dedicated CH Study Web site http://www.epi.bris.ac.uk/chs/.
| 3. CROSS-SECTIONAL HEIGHT DISTRIBUTION |
|---|
|
|
|---|
For the initial analysis, the height data were treated as cross-sectional and were sorted and split into 40 age groups of equal size. This used the entire data set to be compatible with the later longitudinal analysis and ensured that on average each group contained just one measurement per boy. In practice, averaged over the central 30 groups, 66% of boys per group had 1 measurement and 17% had 2.
The first 4 moments of height: mean, SD, skewness, and kurtosis were derived for each group and plotted against the group mean age. Defined as in Appendix (see supplementary material available at Biostatistics online, http://www.biostatistics.oxfordjournals.org), skewness and kurtosis are zero for a normal distribution. The age trends in each moment plot were summarized by a cubic smoothing spline curve with 10 effective degrees of freedom (edf), chosen to give a good fit and a smooth first derivative.
Prior to deriving the moments, height was adjusted for age within each group using linear regression. Throughout, the data have been analyzed both as height and log
height. The main analyses use log height, which was closer than height to a normal distribution. However, 2 figures are back-transformed to the height scale.
The data consisted of 129 700 heights from 3245 boys aged between 9.0 and 20.5 years. A total of 192 heights (0.15%) were excluded as being grossly inconsistent with others for the same boy, leaving 129 508 heights for analysis, an average of 39.9 per boy. The data were sorted by age and split into 40 groups of size 3237 or 3238, with a median age range within each group of 0.19 years.
Figure 2 (top row) shows the first 4 moments of log height for each age group plotted against age, with cubic smoothing spline curves superimposed. The steep slope of the mean curve mid-puberty shows the most rapid growth at that age. The SD curve for log height corresponds to the coefficient of variation (CV) of height and is maximal at about age 14 years. The skewness curve starts and ends at zero, indicating that log height is not skew outside puberty. During puberty, the skewness initially becomes more positive, then switches rapidly to negative skewness and returns to zero postpuberty. The kurtosis curve shows slight (positive) leptokurtosis for most of puberty except for a short period around 14 years when it is clearly platykurtic.
|
The peaked shape of the SD curve is reminiscent of height velocity during the same period (see Figure 1). Cross-sectional log height "velocity" is obtained by differentiating the cubic spline fitted to the mean log height curve, shown in Figure 2 below and to the right of the mean curve and labeled "D mean," where D indicates the derivative. Also shown in the second row of Figure 2 are the first derivatives of the SD curve and the skewness curve, similarly labeled.
The first derivative of the mean curve is indeed similar in shape to the SD curve above it, with both peaking at 14 years, but surprisingly the same is true of the other 2 derivative curves as well. The D SD curve shows the same rise, steep fall then rise as the skewness curve, and the D skewness curve has the same dip at 14 years as the kurtosis curve.
With suitable differentiation, the relationships between the moments can be summarized as follows:
![]() |
where
means "has the same shape as" and
and
indicate second and third derivatives. Thus, the SD, skewness, and kurtosis curves correspond in shape to the first, second, and third derivatives, respectively, of the mean curve.
Section 4 shows that the explanation for this lies in the tempo-shifted pattern of pubertal growth.
| 4. GENERALIZED GROWTH PROCESS |
|---|
|
|
|---|
Consider a boy whose age and height at PHV are average for the population and whose longitudinal height measurements are summarized by the smooth curve
plotted against age t. This is the growth tempo-adjusted individual mean curve (Tanner and others, 1966
, where mean PHV
is given by
. The form of
is assumed to be the same for all boys but with its timing varying according to the boy's tempo of growth as expressed by his age at PHV. This forms the basis for a simple generalized growth process
|
| (4.1) |
where
and
are subject-specific random variables with mean 0, variances
and
, and correlation
, indicating the boy's offset in age and height relative to the mean at PHV. Population samples from the growth process
are then height by age as observed cross-sectionally (with the distribution seen in Figure 2), and the mean of
is the population mean curve. It is in effect a nonlinear random-effects model with
and
as random effects and
a nonparametric smooth curve.
The form of (4.1) implies that growth tempo operates on an additive age scale, with
measured in units of years. However, a multiplicative scale is equally plausible, with rate of maturation from birth inversely proportional to age at PHV, in which case growth tempo
is additive on the log age scale. The 2 alternatives correspond to t in (4.1) being age and log age, respectively. Equally (4.1) is assumed to apply to height or log height. The following algebra applies generally to linear or log age and height scales.
The aims here are (i) to show that the moments of
depend on the derivatives of
, (ii) to compare the means and slopes of
and
at PHV so as to confirm the Merrell (1931)
bias, and (iii) to use the moments of
to estimate
.
The moments of
are expressed in terms of
and its derivatives by expanding (4.1) as a Taylor series (see Appendix, supplementary material available at Biostatistics online http://www.biostatistics.oxfordjournals.org). The leading terms of the first 4 moments are given by
![]() | (4.2) |
where
and
are the mean and SD of
;
,
, and
are the first (velocity), second (acceleration), and third derivatives of
, and
. The SD is arranged to show the excess variance due to
. The second terms for the mean and kurtosis are small (and zero when
, and the square root term for the SD is close to unity, so to first order, for
, the
moment of
is proportional to
. Thus, variation in growth tempo (i.e.
alters the higher moments of the underlying distribution of
in a way that depends on the form of
.
Differentiating
in (4.2) twice gives
![]() | (4.3) |
The mean age at PHV
is defined as the age when height velocity
is at a peak, so by definition
and
. Furthermore, if
is assumed to be symmetric about
then
. Then, from (4.2) and (4.3), at
,
![]() | (4.4) |
This proves that
peaks at
. In words, at the mean age at PHV, population mean height and individual mean height are equal (4.4) and the slopes of both curves are maximal, but for
the population mean curve is shallower in slope than the individual mean curve (4.4). This confirms and generalizes Merrell's (1931)
bias.
From (4.2),
peaks when
is maximal at
, while the kurtosis and D skewness curves are approximately proportional to
which is minimal at
. Thus, the 4 moment curves of
each provide an independent estimate of
—the ages when (i) the mean curve is steepest (positive), (ii) the SD curve is maximal, (iii) the skewness curve is steepest (negative), and (iv) the kurtosis curve is minimal. | 4.2 Results |
|---|
|
|
|---|
Estimates of
were obtained by identifying the relevant turning points in the cubic spline curves for the moments in Figure 2 (D mean, SD, D skewness, and kurtosis). They are summarized in Table 1, for height and log height separately, and are consistently between 14.0 and 14.4 years.
|
| 5. LONGITUDINAL HEIGHT GROWTH |
|---|
|
|
|---|
The data were reanalyzed longitudinally i) to see if growth tempo operates on an additive or multiplicative age scale, ii) to compare the individual and population mean curves, and iii) to test the assumptions underlying the generalized growth process.
The nonlinear mixed-effects model (4.1), along with an implied within-subject residual term, was fitted using the nlme function in R (Beath, 2007; Pinheiro and Bates, 2000
) with
estimated as a fixed-effect regression spline with 11 degrees of freedom (df). To compare the additive and multiplicative age effects, (4.1) was fitted separately for
and
.
To compare the population and individual mean curves,
was also estimated as a cubic regression spline with 11 df fitted to the data treated cross-sectionally. The spline curves for
and
were then back-transformed to give curves of height versus age; these in turn were differentiated, and estimates of mean PHV and mean age at PHV recorded. The key assumption of the growth process (4.1) is that adjusting for growth tempo makes the distribution around the individual mean curve
normal with variance
at all ages. This follows from (4.2), where if
then all the age-varying terms vanish. To adjust for growth tempo, maturational age T was derived such that for child i on occasion j,
(additive) or
(multiplicative). The distribution of log height was then investigated by maturational age group based on
and summarized as before with cubic smoothing splines of 10 edf.
| 5.2 Results |
|---|
|
|
|---|
Model (4.1) for log height with
age provided a dramatically better fit to the data than with
. The within-subject residual SD was 0.0073 for the multiplicative model as against 0.0079 for the additive model, that is, an 8% reduction. Except where stated the following results are for the multiplicative model.
Figure 3 illustrates the estimated subject-specific random effects, with intercept
(as exp(
)) plotted against offset
(as exp(
)) for 3245 boys. The correlation
between
and
was 0.39, showing that boys who were older at PHV were also taller at PHV. The SD of
was
, lower than the mean SD of log height by age in Figure 2. The SD of
was
log years, while for the additive model (4.1) the SD was 1.1 years.
|
Figure 4 compares the cross-sectional (dotted) and longitudinal (solid) spline curves
and
(left) and their first derivatives (right). Their ages at PHV
were both about 14.3 years, but PHV
for the cross-sectional curve (7.4 cm/year) was appreciably less than for the longitudinal curve (9.1 cm/year). Based on log height and log age, PHV was 0.65 units cross-sectionally as against 0.80 units longitudinally.
|
The generalized growth process predicts the moments of
to be approximately proportional to
and its derivatives (4.2). The upper row of Figure 5 shows
and its derivatives, while the lower row shows the predicted
moment curves based on
(4.2). The upper and lower curves are very similar in shape, both to each other and to the upper curves in Figure 2, justifying the statement that the derivatives define the moments. However, Figure 5 predicts a wider range of moments than observed in Figure 2 for the SD and particularly the kurtosis.
|
The distributions of
and
are compared in Figure 6. The first 4 moments of log height are shown grouped by chronological age t (top row, copied from Figure 2) and by maturational age T (bottom row), using the same axis scales. Grouped by T, the 3 higher moments are strikingly more constant across age than the corresponding curves by t. This confirms that the age trends in the higher moments of
seen in Figure 2 are generated by the variability in age at PHV and that adjusting for growth tempo on a multiplicative scale removes the trends and shows the underlying homoscedastic normal distribution about
.
|
| 6. DISCUSSION |
|---|
|
|
|---|
The study has shown that the distribution of height changes cyclically during puberty. Figure 2 shows 1, 2, and 3 turning points, respectively, in the age-varying SD, skewness, and kurtosis. The age when mean height increases fastest is also the age when skewness decreases fastest, which is in turn the age when height variability is greatest and kurtosis least. This age of 14.3 years represents in 4 distinct ways the mean age of PHV for the population and is typical of a boys population measured half a century ago (Tanner, 1978
The age-related skewness pattern in Figure 2 confirms Boas's (1892a
, b
) observation that skewness increases when puberty starts. It subsequently falls sharply and then rises again, returning to where it started. The kurtosis pattern in Figure 2 has not previously been reported. Yet, the Tanner–Whitehouse (Tanner and others, 1966
), British 1990 (Freeman and others, 1995
) and Dutch 1997 (Fredriks and others, 2000
) references all failed to detect skewness or kurtosis and treated height as normally distributed. By contrast, the US Centers of Disease Control 2000 reference (Kuczmarski and others, 2000
) did adjust for skewness, but the age-related changes were complex and showed nothing of the cyclical pattern in puberty seen here. This demonstrates the very large sample sizes needed to detect skewness in height, which probably relates to its small variability (Pan and Cole, 2004
) (log
in Figure 2, i.e.
). The impact of skewness on height centiles is relatively modest, less than 0.5% on the 2
or 98
centile, and the dip in kurtosis at the age at PHV (Figure 2) is even harder to detect and has a tiny effect on the extreme centiles. This justifies the decision to treat height as normally distributed in growth references.
However, the cyclical change in the height distribution has wider implications. A remarkable feature of the age-related changes in the moments is that each corresponds in shape to the derivative of the previous moment. So the rise and fall in the SD are synchronized with the rise and fall in height velocity and similarly for skewness and SD "velocity" and for kurtosis and skewness "velocity." These observations together show that age-related changes in the variability, skewness, and kurtosis of height are similar to those for height velocity, height acceleration, and the third derivative of height, respectively. Not only does this confirm and generalize Boas's observation that skewness changes, but also predicts that skewness changes in line with height acceleration.
The pattern can be explained mathematically with a simple growth process that treats height growth as a nonlinear function where individuals vary in growth tempo. Fitting this to the longitudinal data as a mixed-effects model shows clearly that growth tempo varies on a multiplicative not an additive age scale, inversely proportional to each child's age at PHV. The finding generalizes in the following way: in any nonlinear process, where individuals "grow" following the same pattern but at different rates, the shape of the underlying individual mean "growth" curve defines time changes in the moments of the cross-sectional distribution.
If, for example, the growth curve had nonzero first and second derivatives but zero higher derivatives, that is, the growth curve was linear or quadratic, then acceleration would be zero and skewness constant. It would not necessarily be zero as the growth process superimposes changes on the underlying distribution, which like height may already be skewed before puberty. So the age changes in height SD, skewness, and kurtosis arise from the complex shape of the mean individual height curve.
The reason why the multiplicative (log) age model fits so much better than the additive model is at least partly because it accounts for the known correlation between PHV and age at PHV seen in Figure 1 (Tanner, 1962
). Figure 7 illustrates predicted height and height velocity for individuals whose age at PHV is early, average, and late, that is,
,
and
. Back-transforming to the age scale has led to a steeper height curve for early maturers compared to late and a corresponding inverse trend in PHV.
|
The effect of a multiplicative scale for age at PHV. Back-transforming the log age scale generates steeper predicted height curves for early maturers, resulting in the known inverse correlation between PHV and age at PHV (Tanner, 1962
The growth process also clarifies the relationship between the population mean curve and the individual mean curve. Merrell (1931)
showed in a mathematical tour de force that assuming a logistic growth pattern, the population mean curve is flattened relative to the individual mean curve, due to the time shifts in individual curves. Equation (4.4) shows the same relationship for a smooth growth curve of unspecified shape, where the population mean curve and individual mean curve coincide at the mean age and mean height at PHV, but the population curve is shallower than the individual curve unless
is zero. This generalizes Merrell's finding, as seen in Figure 4 with the present data.
Apart from their theoretical interest, these moment-derivative relationships have a practical application in the construction of growth references. The LMS method, for example, summarizes growth reference centiles in terms of curves representing the median M, coefficient of variation S, and skewness L by age (Cole and Green, 1992
). The S curve is closely similar to the log SD curve, and the L curve, being the Box–Cox power to transform to normality, is essentially the inverse of the skewness curve. As shown here, derivatives of the M curve contain information about the likely shapes of the S and L curves. The M curve is the most precise of the 3 to estimate, being based on the first moment of the distribution, whereas the S curve refers to the second moment and the L curve to the third. This means that a relatively modest sample can be used to estimate the M curve, then the (underpowered) S curve can "borrow strength" from the M velocity curve and similarly for the (even more underpowered) L curve from the M acceleration curve. It could be done by using each derivative of the M curve as a covariate with the GAMLSS package (Rigby and Stasinopoulos, 2005
), which extends the LMS method by including covariates and optionally modeling kurtosis as well. It would be a potentially useful way of improving the precision of the S and L curves and/or reducing the required sample size.
In conclusion, a very large sample of heights has shown cyclical changes in the higher moments of the distribution through puberty, the patterns of change corresponding in shape to derivatives of the mean height curve. This relationship is explained by a tempo-adjusted growth process. As Major Greenwood stated prophetically in his appendix to Friend's (1935)
book describing the CH heights, "I think that Dr. Friend has provided material which will be of service to workers for some time to come."
| FUNDING |
|---|
|
|
|---|
Funding to pay the Open Access publication charges for this article was provided by a private research account.
| ACKNOWLEDGMENTS |
|---|
We thank George Davey Smith for his initiative in following up the CH cohort and making the data available and Yoav Ben-Shlomo for pointing us to the quotation by Major Greenwood. Comments by 2 anonymous referees considerably improved the paper. Research at the University College London Institute of Child Health and Great Ormond Street Hospital for Children National Health Service (NHS) Trust benefits from Research and Development funding received from the NHS Executive. Tim J. Cole is funded by Medical Research Council program grant G9827821. Conflict of Interest: None declared.
| REFERENCES |
|---|
|
|
|---|
-
Beath KJ. Infant growth modelling using a shape invariant model with random effects. Statistics in Medicine (2007) 26:2547–2564.[CrossRef][Medline]
Boas F. The growth of children. Science (1892a) 19:256–257.[Medline]
Boas F. The growth of children. Science (1892b) 19:281–282.[CrossRef][Medline]
Bowditch HP. The growth of children studied by Galton's percentile grades. In: 22nd Annual Report of the State Board of Health of Massachusetts (1891) Boston, MA: Wright and Potter.
Cole TJ, Freeman JV, Preece MA. British 1990 growth reference centiles for weight, height, body mass index and head circumference fitted by maximum penalized likelihood. Statistics in Medicine (1998) 17:407–429.[CrossRef][ISI][Medline]
Cole TJ, Green PJ. Smoothing reference centile curves: the LMS method and penalized likelihood. Statistics in Medicine (1992) 11:1305–1319.[ISI][Medline]
Fredriks AM, Van Buuren S, Burgmeijer RJ, Meulmeester JF, Beuker RJ, Brugman E, Roede MJ, Verloove-Vanhorick SP, Wit JM. Continuing positive secular growth change in The Netherlands 1955–1997. Pediatric Research (2000) 47:316–323.[ISI][Medline]
Freeman JV, Cole TJ, Chinn S, Jones PRM, White EM, Preece MA. Cross-sectional stature and weight reference curves for the UK, 1990. Archives of Disease in Childhood (1995) 73:17–24.[Abstract]
Friend GE. The Schoolboy: A Study of His Nutrition, Physical Development and Health (1935) Cambridge, UK: Heffer.
Kuczmarski RJ, Ogden CL, Grummer-Strawn LM, Flegal KM, Guo SS, Wei R, Mei Z, Curtin LR, Roche AF, Johnson CL. National Center for Health Statistics. In: CDC Growth Charts: United States (2000) Hyattsville, MD.
Merrell M. The relationship of individual growth to average growth. Human Biology (1931) 3:37–70.
Pan H, Cole TJ. A comparison of goodness of fit tests for age-related reference ranges. Statistics in Medicine (2004) 23:1749–1765.[CrossRef][ISI][Medline]
Pinheiro JC, Bates DM. Mixed-Effects Models in S and S-PLUS (2000) New York: Springer.
Rigby RA, Stasinopoulos DM. Generalized additive models for location, scale and shape (with discussion). Applied Statistics (2005) 54:507–544.
Sandhu J, Ben-Shlomo Y, Cole TJ, Holly J, Davey Smith G. The impact of childhood body mass index on timing of puberty, adult stature and obesity: a follow-up study based on adolescent anthropometry recorded at Christ's Hospital (1936–1964). International Journal of Obesity (2006a) 30:14–22.[CrossRef][ISI][Medline]
Sandhu J, Davey Smith G, Holly J, Cole TJ, Ben-Shlomo Y. Timing of puberty determines serum insulin-like growth factor-1 in late adulthood. Journal of Clinical Endocrinology and Metabolism (2006b) 91:3150–3157.
Tanner JM. Growth at Adolescence (1962) Oxford: Blackwell.
Tanner JM. Foetus into Man: Physical Growth from Conception to Maturity (1978) London: Open Books.
Tanner JM. A History of the Study of Human Growth (1981) Cambridge, UK: Cambridge University Press.
Tanner JM, Whitehouse RH, Takaishi M. Standards from birth to maturity for height, weight, height velocity, and weight velocity: British children, 1965 parts I and II. Archives of Disease in Childhood (1966) 41:454–471. 613–635.[ISI][Medline]
Received July 28, 2006; revised January 15, 2007; revised March 16, 2007; accepted for publication April 25, 2007.
![]()
CiteULike
Connotea
Del.icio.us What's this?
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||










