Bayesian SITAR model fit

Introduction

The Superimposition by Translation and Rotation (SITAR) model is a shape-invariant, nonlinear mixed-effects growth curve model that fits a population average curve to the data and aligns each individual’s growth trajectory to this underlying population average curve via a set of random effects (T. J. Cole et al., 2010). The ‘bsitar’ package (Sandhu, 2023) implements Bayesian estimation of the SITAR model and can fit a univariate model (analysis of a single outcome), univariate_by model (analysis of a single outcome separately for subgroups defined by a factor variable such as gender), and multivariate model (joint modeling of two or more outcomes). A detailed introduction to the Bayesian SITAR model, including the biological and mathematical description of the model, Markov Chain Monte Carlo (MCMC) estimation, and prior specifications, is provided in the vignette Bayesian SITAR model - An Introduction.

In this vignette, we fit a univariate Bayesian SITAR model using the ‘bsitar’ package (Sandhu, 2023) and compare the findings with the results of the non-Bayesian SITAR model (fitted via the ‘sitar’ package (T. Cole, 2022)). Hereafter, the Bayesian and non-Bayesian models estimated by the ‘sitar’ and ‘bsitar’ packages are referred to as sitar and bsitar models, respectively. The purpose is to demonstrate that both ‘sitar’ and ‘bsitar’ packages perform similarly in estimating population average and individual-specific trajectories and growth spurt parameters (timing and intensity). The aim of this exercise is to ensure that our ‘bsitar’ package performs as expected (i.e., no systematic difference in results when compared with the sitar model) before we fit a more complex model (such as a multivariate model), for which no direct comparison is possible between sitar and bsitar models. This is because the ‘sitar’ package is restricted to univariate model fitting.

We study the Berkeley data, which comprise repeated growth measurements taken on males and females from birth to adulthood. The Berkeley dataset is included in both the ‘sitar’ and ‘bsitar’ packages. To compare our findings directly with the results shown in the vignette for the sitar model Fitting models with SITAR, we analyze the exact same data, which is a subset of the Berkeley data with height measurements for females (between 8 and 18 years of age, every six months).

Model fit

If needed, install the ‘bsitar’ and ‘sitar’ packages (along with any other required packages, such as ggplot2).

# Load required packages
library(bsitar)  # To fit Bayesian SITAR model
library(sitar)   # To fit frequentist SITAR model
library(ggplot2) # For plots

Fit sitar and bsitar models using ‘sitar’ and ‘bsitar’ packages.

if(existsdata) {
  if(exists('berkeley_exdata')) data <- berkeley_exdata else stop("data does not exist")
} else {
  data <- na.omit(berkeley[berkeley$sex == "Female" & 
                         berkeley$age >= 8 & berkeley$age <= 18, 
                       c('id', 'age', 'height')])
}


# Fit frequentist SITAR model, 'model_frequ' ('sitar' package)
model_frequ <- sitar::sitar(x = age, y = height, id = id, data = data, df = 3)


# Fit Bayesian SITAR model 'model_bayes' ('bsitar' package) 
# To save time, the model is fit by running two parallel chain with 1000  
# iterations per chain. 

if(existsfit) {
  if(exists('berkeley_exfit')) model_bayes <- berkeley_exfit else stop("model does not exist")
} else {
  model_bayes <-  bsitar(x = age, y = height, id = id, data = data, df = 5,
                chains = 2, cores = 2, iter = 2000, thin = 4, refresh = 100,
                sample_prior = "no",
                seed = 123)
}

Model summary

Summarize the sitar model

SITAR nonlinear mixed-effects model fit by maximum likelihood
  Call: sitar::sitar(x = age, y = height, id = id, data = data, df = 3)
  Log-likelihood: -1112.306
  Fixed: s1 + s2 + s3 + a + b + c ~ 1 
         s1          s2          s3           a           b           c 
 27.3650000  41.5161852  24.5786864 133.7696674  -0.9158315   0.5093028 

Random effects:
 Formula: list(a ~ 1, b ~ 1, c ~ 1)
 Level: id
 Structure: General positive-definite, Log-Cholesky parametrization
         StdDev    Corr         
a        5.8650199 a      b     
b        0.8602644  0.150       
c        0.1233394  0.341 -0.672
Residual 0.5150428              

Number of Observations: 770
Number of Groups: 70 

Summarize the bsitar model

 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: height ~ SITARFun(age, a, b, c, s1, s2, s3) 
         a ~ 1 + (1 | 1 | gr(id, dist = "gaussian"))
         b ~ 1 + (1 | 1 | gr(id, dist = "gaussian"))
         c ~ 1 + (1 | 1 | gr(id, dist = "gaussian"))
         s1 ~ 1
         s2 ~ 1
         s3 ~ 1
   Data: structure(list(id = structure(c(1L, 1L, 1L, 1L, 1L (Number of observations: 770) 
  Draws: 2 chains, each with iter = 2000; warmup = 1000; thin = 5;
         total post-warmup draws = 400

Multilevel Hyperparameters:
~id (Number of levels: 70) 
                             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(a_Intercept)                  6.11      0.53     5.10     7.19 1.01      107      182
sd(b_Intercept)                  0.88      0.08     0.75     1.04 1.00      142      229
sd(c_Intercept)                  0.13      0.01     0.11     0.15 1.00      272      294
cor(a_Intercept,b_Intercept)     0.17      0.12    -0.06     0.40 1.00      155      194
cor(a_Intercept,c_Intercept)     0.32      0.11     0.09     0.52 1.00      171      239
cor(b_Intercept,c_Intercept)    -0.65      0.07    -0.77    -0.50 1.01      273      365

Regression Coefficients:
             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
a_Intercept    133.52      1.03   131.44   135.58 1.01      119      239
b_Intercept     -0.93      0.10    -1.12    -0.74 1.00      155      191
c_Intercept      0.50      0.03     0.43     0.56 1.00      227      280
s1_Intercept     3.87      0.09     3.70     4.04 1.01      333      348
s2_Intercept    -0.32      0.01    -0.35    -0.29 1.00      454      407
s3_Intercept    31.85      0.71    30.54    33.24 1.00      353      382

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.52      0.02     0.49     0.55 1.00      291      390

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

The figure below shows growth curves adjusted via the random effects. As described earlier (see Introduction), the random effects align individual-specific growth trajectories toward a common population average (mean) trajectory. The three growth patterns that SITAR adjusts via random effects are size, timing, and intensity:

• Size is an individual’s mean height compared to the average (measured in cm).
• Timing is the individual’s age at peak height velocity (i.e., the age when they are growing fastest) compared to the average (measured in years).
• Intensity is their growth rate compared to the average (measured as a fraction, or percentage divided by 100).

Adjusted curves aligned with the mean curve indicate that most of the inter-individual variability has been captured. This is achieved by: 1) shifting the high curves down and the low curves up to adjust for the size, 2) shifting the early curves later and the late curves earlier (translation) to adjust for the timing, and 3) steep curves are made shallower by stretching them along the age scale whereas shallow curves are made steeper by shrinking them along the age scale (rotation) to adjust for the intensity. This explains the name SITAR - Super Imposition by Translation And Rotation. The curve adjustment involves translation (shifting the curves up/down and left/right) and rotation (making them steeper or shallower), and the effect is to superimpose the curves. The complexity of the spline curve’s shape is determined by the degrees of freedom (df) argument in the sitar and bsitar call. In our example, the models are fit using five degrees of freedom.

Observed curves along with the superimposed adjusted curves (via random effects) for the sitar and bsitar models. The observed curves are shown as solid light grey lines whereas the adjusted curves for the sitar model are displayed as solid black lines. The adjusted curves for the bsitar model are shown as dotted black lines.

Random effects (size, timing and intensity)

The effect of SITAR adjustment in individual curves via random effects is shown more clearly in the figure below, which displays observed and adjusted curves for the tallest and shortest individuals along with the population average (mean) curve. The corresponding estimates of random effects (size, timing, and intensity) are provided in the tables below. Note that intensity c, multiplied by 100, is a percentage. Results for both models are very similar to each other. Results for the sitar model are presented below.

For the shortest individual, size is 14.6 cm smaller than the average size, timing of the spurt is -1.1 year earlier than average timing, and the growth rate (intensity) is 27.3 % faster than average growth rate. In contrast, the size for tallest individual is -9.4 cm larger than the average size and the spurt occurs 1.2 year later than average timing with growth rate -14.9 % slower than average intensity.

For a broader comparison between sitar and bsitar models, random effects for the first ten individuals estimated are shown below (see Tables). The association (correlation) between random effects (across all individuals) is shown as a scatterplot in the figures, and corresponding estimates are provided in the tables. Results show that timing is negatively correlated with intensity, indicating that early puberty is associated with faster growth.

Observed and adjusted curve (via random effects) curves along with the population averag curves for sitar and bsitar models. The observed curves for all individuals are shown as solid light grey lines whereas curves for the tallest and shortest individuals are displayed using the solid black lines. The adjusted curves for the tallest and shortest individuals are shown as dotted black lines, and the population average curve is shown as solid black line.

Correlation scatterplot of random effects for the sitar model

Correlation scatterplot of random effects for the bsitar model

Growth curves (distance and velocity)

Population average (mean) distance and velocity curves

Population average distance curves for the sitar and bsitar models

The estimated distance curve (see figure below) shows that the increase in height follows a sigmoid pattern, with maximum gain in height occurring around 12 years of age. This is confirmed by the velocity curve shown in the figure (below), which shows a clear spurt in growth rate that peaks at around 11.7 years of age (see Population average (mean) estimates).

Population average velocity curves for the sitar and bsitar models

Individual specific distance and velocity curves

Individual specific distance curves for the sitar and bsitar models

Individual specific velocity curves for the sitar and bsitar models

Individual-specific distance and velocity curves (shown below) follow a pattern similar to the population average distance and velocity curves described earlier. Individual velocity curves (shown below) show that all individuals experience a pubertal growth spurt, but the timing of the spurt varies between 10 and 14 years, with varying growth patterns among individuals. Some individuals are consistently taller than average, while others are consistently shorter. While some start relatively short and become taller, others have an opposite growth pattern, i.e., they are taller at an early age and achieve less height in later years compared to those who were short at an early age.

Growth parameters (timing and intensity)

References