Suppose the data is collected from n subjects, and each of the subjects was only observed once at a random time point. Specifically, let the observed data or measurement for the ith subject be y i , and the observed time be t i , i = 1 , 2 , … , n. Since y i are from different subjects, they are independent to each other. Additional information except time are possible, such as gender or treatment group. They are represented by a p-dimensional x i for the ith subject. To introduce heteroscedasticity, we use a smooth function  g ( ν i ) to characterize the standard deviation of the regression error of the ith subject, where  v i is some index variable. The index variable can be observed such as time, or unobserved such as the mean of the corresponding subject from the linear model. The smooth function  g ( · ) areis parametrized by shape-restricted splines.

The growth pattern against time is often non-linear and present certain shape restrictions, such as increasing with time. This can be realized by using shape-restricted spline of time in the mean model. Once the degree and degree of freedom of spline bases being chosen, it can be represented using the same format as the parametric part: linear combination of coefficient and spline bases. Using spline bases in the mean pattern has been discussed extensively in the previous work [33, 20, 21, 53], and our focus is more on using the spline on the variance part. For the simplicity of expression, we choose to include the spline bases of time as part of the vector x i instead of explicitly show them separately. The process to chose the degree and degree of freedom is the same as for the splines in the error term.

Using I-splines as an example, a heteroscedastic linear model is (3.1) y i = x i ⊤ β + ε i , i = 1 , … , n , ε i ∼ N ( 0 , g 2 ( v i , θ ) ) , g ( v i ; θ ) = θ 0 + ∑ k = 1 K θ k I 2 , k ( v i ) , where  β is a p-dimensional regression coefficient vector for  x i , ε i is the normally distributed regression error with mean zero and standard deviation g ( ν i ; θ ), { I 1 , k ( · ) , k = 1 , … , K } are I-spline bases with K degrees of freedom, and θ = ( θ 0 , … , θ K ) is a ( K + 1 )-dimensional coefficient vector. The degree and degrees of freedom for each spline bases need to be selected using model selection method introduced in Section 3.3, and the internal knots are evenly spaced. The coefficients θ can be restricted to control the shape of the heteroscedasticity as a function of ν i . For example, if the variance increases with the mean (or time), the coefficients θ can be restricted to be non-negative.

If concavity or convexity is desired, the I-splines can be replaced with C-splines and a linear term of time with appropriate restrictions on the coefficients, as introduced in Section 2. Interaction terms can be introduced in the mean function to allow the covariates have time-varying coefficients [20, 21].

When more than one measures were collected from each of the n subjects, the observed data will have a clustered structure. Let the number of repeated measures on the ith subject be n i , and it might be different for each subject and sometimes might be small or even just 1. The jth observation of the ith subject is y i , j which is collected at time t i , j , where j = 1 , … , n i , i = 1 , … , n. Same as before, we still use spline bases to estimate the error variance in order to put the shape restrictions, but with a linear mixed-effects model to account for the dependence structure within the dataset. We use the matrix notation for simplicity of demonstration: y i = y i , 1 . . . y i , n i , X i = x i , 1 ⊤ . . . x i , n i ⊤ , ε i = ε i , 1 . . . ε i , n i . Again using I-Splines as an example, the model for the ith subject is (3.2) y i = X i β + Z i b i + ε i , i = 1 , … , n , ε i ∼ MVN 0 , diag ( g 2 ( v i , 1 , θ ) , … , g 2 ( v i , n i , θ ) ) , g ( v i , j , θ ) = θ 0 + ∑ k = 1 K θ k I k ( v i , j ) , j = 1 , … , n i , where x i , j is a p dimensional covariate vector for fixed effects which can include the spline bases of time, Z i is an  n i × q design matrix for random effects, β is a p dimensional fixed effects vector, b i is a q-dimensional random effects vector with covariance matrix B parameterized by vector α, and MVN is the multivariate normal distribution. Shape restrictions can be applied on the coefficient of spline bases, and K need to be selected using the model section method in Section 3.3. Other notations are the clustered analogs to those in Equation (3.1).

A special choice for the index variable v i j is the mean of the response variable. For a mixed-effects model, this mean can be conditional on the random effects or not. If random effects are conditioned on, the mean is (3.3) μ i , c = E [ y i ∣ b i ] = X i β + Z i b i ; otherwise it is (3.4) μ i , m = E [ y i ] = X i β . For ease of referencing, we call them conditional mean and marginal mean, respectively. When the error variance changes with the marginal mean, the response variable still has a multivariate normal distribution. If the conditional mean is in the error variance structure, there is dependence between the random effects and the error term, and the response variable no longer has a multivariate normal distribution. To calculate the likelihood function for this situation, as needed in AIC and BIC calculations, numerical integration is needed. See details in Section 3.3.

Same as for the indepdent data, C-splines can be used for concavity or convexity shape restriction, and interaction terms in the mean function can allow time-varying coefficients for the covariates.

The maximum likelihood method can be used to get parameter estimates in theory as long as appropriate restrictions on the coefficients are imposed to enforce the shape restrictions. To obtain the maximum likelihood estimator, we propose an iteratively reweighted fitting procedure that takes advantage of existing software packages for linear mixed-effects models allowing weights. This method is flexibile to deal with different scenarios including the error variance changing with conditional or marginal mean. It can also be easily computed since no closed-form solutions need to be derived. The steps are summarized in Algorithm 1 for clustered data when the error variance is changing with the mean. We use  μ i in the algorithm to represent either the conditional mean or the marginal mean, and  μ ˆ i for the estimated value of the mean. Algorithm for independent data is similar and simpler, and will not be repeated here. The shape restrictions on the heteroscedasticity (and the mean model) can be imposed with a constrained optimizer, such as the constrOptim() function in R.

To construct reference quantiles in a GCA, we suggest using parametric bootstrap. This is very similar to the resampling-subject bootstrap (RSB) method [20, 53] since the bootstrap sample is generated on subject level to maintain the cluster structure. The main difference is that our bootstrap sample is generated with estimated parameter, instead of the residuals. This is because the error variance in our model is estimated with spline bases. The bootstrap method avoids deriving the likelihood function and the Hessian matrix [39, p. 227], which is challenging when the error variance changes with the conditional mean. This is of particular importance when some of the estimated θ’s are on the boundaries of the constrained parameter space [2]. Since we focus on the final fitted curve instead of the basis coefficients, bootstrap provides a natural solution for the quantitles of the fitted curve regardless of whether some of the estimated θ’s are on the boundaries.

Algorithm 2 summarizes the steps to get one parametric bootstrap sample for clustered data. Same as before, the algorithm for independent data will be similar and will not be displayed here. Repeating this process gives a large sample of bootstrap copies of the point estimates of the model parameters. Their empirical standard deviations are then used as the standard errors of the model parameter estimates. For mixed-effects models, parametric bootstrap method has better performance compared with bootstrap methods that only re-sample observations or residuals, as it produces more accurate standard deviation of estimated parameters, and closer-to-nominal coverage rates for confidence intervals [11, 54, 47].

After the model fitting process, model checking can be done using the residuals. The standardized residuals, i.e., the residuals divided by the estimated error standard deviation, should follow a standard normal distribution, and their normality can be checked visually using a normal Q-Q plot, or other normality testing methods.

With the internal knots evenly spaced, the values of degree and degree of freedom are needed to generate the spline bases. They should depend on the sample size n and the number of observations each subject has n i for clustered data. Additionally, for clustered data, candidate models can either have the conditional mean or the marginal mean in the error variance function in Model (3.2).

Different values of spline degree and degree of freedom, and the choice of using conditional mean or marginal mean, will significantly impact the model fitting results. The popular model selection criteria for such problems are AIC and BIC [33] as they consider both the fitting accuracy and the model complexity. It can give similar results to the ‘deleting subject cross validation’ method, and is faster to compute [33, 20]. Details are as follow: AIC = − 2 log L + 2 P , BIC = − 2 log L + P log n , where L is the likelihood function of the fitted model, P is the number of parameters, and n is the sample size. Models with smaller AIC or BIC are preferred.

The numerical integration is needed for the scenario when the model has error variance changing with the conditional mean (3.3). From the definition in Model (3.2), the error variance now depends on the random effects, and the likelihood function should be: L = ∏ i = 1 n ∫ − ∞ + ∞ f ( y i ∣ b i ) f ( b i ) d b i , b i ∼ MVN 0 , B ˆ , y i ∣ b i ∼ MVN X i β ˆ + Z i b i , g 2 ( v i , θ ˆ ) , where  f ( y i ∣ b i ) is the probability density function (pdf) of the response vector  y i conditioning on the random effects  b i , and  f ( b i ) is the pdf of the random effects  b i . The distributions are showed in the equation. There is no closed-form result for this integral, but it can be numerically calculated.

This study mimics a scenario where each subject has one measure; see the fetal pancreas length application in Section 5.1. The data generating model was y i = β 0 + β 1 x 1 i + β 2 x 2 i + ε i , i = 1 , … , n , where  y i is ith response, x 1 i is a Bernoulli ( 0.5 ) variable, x 2 i is a Uniform ( 0 , 2 ) variable, ( β 0 , β 1 , β 2 ) = ( 1 , 1 , 1 ), and ε i is a zero-mean normally distributed error term with heteroscedasticity. Specifically, the standard deviation of ε i is  g ( μ i ), where μ i = β 0 + β 1 x 1 i + β 2 x 2 i . Three functional forms were considered for g: g 1 ( μ ) = 0.25 ( μ − 0.9 ), g 2 ( μ ) = 0.02 ( μ 3 + 1.2 ), and g 3 ( μ ) = 0.1 ( 5 Φ ( μ − 2 ) / 0.3 ) + 1 ), where  Φ ( · ) is the cumulative distribution function of the standard normal distribution. The specific parameter values in these functions were chosen such that the function values were positive, the ratio of the maximum value over the minimum value was around 30, and the resulting signal-to-noise ratio in the linear regression model was around 3. Two sample sizes were considered n ∈ { 100 , 200 }.

For each configuration, 1000 datasets were generated. For comparison, both the naive linear regression model with heteroscedasticity ignored and the proposed GCA with shape-restricted heteroscedasticity were fitted to each dataset. For the naive model, the variances of the regression coefficients were obtained using the robust estimator to account for the heteroskedasticity [52, 25] as implemented in the R package sandwich [59, 57, 58]. In the proposed method, I-spline bases were used to enforce monotonicity. As more degrees of freedom were needed for fitting more complex patterns with satisfying accuracy, we picked quadratic I-spline basis with 2, 3, and 7 degrees of freedom when fitting the three g ( · ) patterns, respectively, with evenly spaced internal knots. Parametric bootstrap was used to calculate the standard deviation of the estimates and to construct 95% confidence interval of the regression coefficients and the error variance curve. The number of bootstrapping replicates was 1000.

Table 1 summarizes the empirical bias, empirical standard error (se), estimated standard error ( se ˆ), and the coverage percentage (CP) of the 95% confidence intervals of parameter estimates. The bias from the proposed method is close to zero under all the scenarios. The estimated standard deviations from parametric bootstrap se ˆ is close to the empirical value from the proposed method, and the CP is close to the nominal level 95%. The point estimate of β 2 has lower variation than that for β 1 , which is expected because the continuous covariate  x 2 provides more information than the binary  x 1 . As sample size increases, all standard errors decreases. In comparison with proposed method, the naive method leads to much higher standard deviations in the regression coefficient estimation. Although the 95% confidence intervals from the naive method seem to have appropriate coverage percentage, they are much longer than those from the proposed method as evident from the standard errors.

Figure 1 displays the estimated heteroscedasticity form and the averaged point-wise 95% confidence intervals from the proposed method using parametric bootstrap. The averaged estimates represented by the solid lines are close to the true curves in dashed lines, and the true curves lie within the 95% point-wise confidence intervals. The wider intervals near the right boundary, especially in the case of g 3 , reflect that this is a challenging situation; the functional form is convex first and concave later, which needs more degrees of freedom to capture.

For clustered data, we considered a setting where each subject has repeated measures; see the chicken weight application in Section 5.2. The data generating model was a mixed-effects model with a random effect at the subject level, y i j = β 0 + β 1 x 1 i + β 2 x 2 i j + b i + ε i j , i = 1 , … , n , j = 1 , … , 5 , where y i j is the response variable for the jth observation from the ith subject, x 1 i is a subject-level covariate generated from Bernoulli ( 0.5 ), x 2 i j is an observation-level continuous covariate generated from Uniform ( 0 , 5 ), ( β 0 , β 1 , β 2 ) = ( 1 , 1 , 1 ), b i is the subject-level random effect generated from N ( 0 , σ b 2 ) with σ b = 0.1, and  ε i j is a zero-mean normally distributed error term with heteroscedasticity. The standard deviation of ε i j was set to be g ( ν i j ), where g ∈ { g 1 , g 2 , g 3 } is same as in the independent data setting, ν i j is either the conditional mean or the marginal mean of the jth observation of the ith subject defined in Equation (3.3) and (3.4), respectively. The total number of subjects was set to be again  n ∈ { 100 , 200 }.

Three different estimation methods were compared. The first is the naive method that fits a linear mixed-effects model with constant error variance. To better capture the standard error, the robust sandwich estimator is used [48]. The second method uses the function lme() from the R package nlme, which provides some preset forms of heteroscedasticity. For settings where g 1 and g 2 are used and μ is the conditional mean, we used the true setup to specify the error variance when fitting the model. For other settings, the correct setup is not available in the lme() function, but we still used a power function of the conditional mean to specify the heteroscedasticity. The third method is the proposed method, where quadratic I-spline bases were used with degrees of freedom 5 for  g 1 and g 2 , and 7 for  g 3 . Internal knots of spline bases were chosen to be evenly spaced. The number of replications for parametric bootstrap was 1000. For each configuration, results were obtained for 1000 datasets.

Tables 2 and 3 summarize the simulation results for the variance as a function of marginal mean and conditional mean, respectively. All three methods seem to give unbiased point estimates for the regression coefficients, as they all have correct specification of the regression model. Their differences are in their uncertainty levels and coverage percentages of the 95% confidence intervals. The naive method has the worst performance since it did not consider the heteroscedasticity at all. The method with function lme() performs better than the naive method, but is still not satisfactory due to the misspecified heteroscedasticity form through the limited choices offered by lme(). Even in settings where the heteroscedasticity is correctly specified by the power function  g 2 ( ) with conditional mean as the index variable, see Table 3, its coverage percentages for β 2 are much lower than 95%. In contrast, in all three settings, the proposed method gives point estimates with lower standard errors as well as confidence intervals with coverage percentage close to 95%. As the sample size n increases, the performance becomes better as expected.

Also reported in Tables 2 and 3 are estimates of the standard deviation σ b of the random effects. It is known that the standard errors of the random-effects variance parameters are hard to get and confidence intervals constructed from profile likelihood or parametric bootstrap are preferred [4]. So here we focus on the point estimate of  σ b . When the heteroscedasticity takes more complicated forms such as g 2 and g 3 , the proposed method has much smaller bias than the lme() method. In some settings such as g = g 2 ; the lme() estimates of σ b have bias but lower variation compared to those from the proposed method. This echos that caution is needed when using standard errors of the random effect variance. The empirical standard errors of the point estimates from the proposed method decrease as the sample size increases, but apparently not at the rate of 1 / n , suggesting that a larger sample size is needed for the asymptotic properties of the random-effect variance estimator to hold.

Figure 2 displays the fitted heteroscedasticity from the proposed method. The two panels show the results with the index variables being the marginal mean and the conditional mean, respectively. Similar to the independent data scenarios, the estimated curve is close to the true curve, and is within the averaged  95 % point-wise confidence intervals. The intervals are narrower than those in the independent data scenarios as the repeated measures provide more information.

Fetus pancreatic dysplasia and hypertrophy, i.e., abnormality of fetal pancreas, are associated with congenital malformations [31, 19]. The fetal pancreas growth curve is a critical tool for prenatal screening for disorders. The growth pattern of fetal pancreas’ lengths during the prenatal period has been investigated [12, 22, 15], but not the changing variation of the measurements. The proposed method allows capturing the fetal pancreas’ growing patterns in both the mean level and the variation level.

One of the authors of this paper, Dr. Wenling Song from the Second Hospital of Jilin University, collected the healthy pancrea length data and provided it for analysis in this paper. The data were collected from 44 pregnant women at different stages of pregnancy who visited the Second Hospital of Jilin University in China during April to July of 2012. The data is provided by Dr.Wenling Song, who is also the author of this paper. No patient showed any external pathology or anomaly. The dataset contains a single measure of the fetal pancreas’ length from each patient. Figure 3 (I) shows the pancreas’ lengths in millimeters versus the pregnant duration in days. It is reasonable to assume that the growth speed slows down as the fetus matures, in which case the growth curve would be increasing and concave. It is also reasonable to assume that the variation increases with time. Therefore, we use C-splines to model the mean growth level and use I-splines for the heteroscedasticity. Specifically, the model is y i = β 0 + β 1 t i + ∑ k = 1 K 1 β k + 1 C k ( t i ) + ε i , ε i ∼ N 0 , θ 0 + ∑ k = 1 K 2 θ k I k ( t i ) 2 , where y i is the i t h length measurement, t i is the corresponding time with linear coefficient β 1 , C k ( t i ) is the kth C-splines basis evaluated at t i with coefficient β k + 1 , k = 1 , … , K 1 , ε i is the independent error term with heteroscedasticity, and I k ( t i ) is the kth I-spline basis evaluated at t i with coefficient θ k , k = 1 , … , K 2 . The spline bases were selected using the method stated in Section 3.4. The internal knots of both spline bases were evenly spaced. The degrees of freedom K 1 and K 2 for the C-splines and I-splines, respectively, were both chosen to be 4 by BIC.

Figure 3 shows fitted results for the growth curve of the pancreas length. By shape restrictions, the fitted curve is increasing and concave while the variance is increasing over time. The estimated curve and the point-wise confidence intervals in panel (I) accurately capture the mean and variation pattern. The residual plot with fitted  95 % confidence intervals in panel (II) shows good performance on estimating the heterogeneous pattern of error standard deviation. The Q-Q plot of the standardized residuals in panel (III) shows no alarming deviation from the normality.

The fitted results show accurate estimate of the pancrea growth curve along with pregnancy days, including the mean and the quantiles. This can provide a better guidance of screening abnormality specifically for babies in that area.

The chicken weight data, which is available in R package datasets, is a classic example of clustered data for GCA [14, p. 4]. It contains the body weights in grams of 50 chicks measured at birth date and every other day thereafter until day 20, plus an additional measurement on day 21. These 50 chicks were divided into four groups to have different protein diets, and the scatter plots of weight versus time for each group are shown in Figure 4. We removed Chick No. 24 as an outlier from the original dataset since its weight stopped increasing after day 6. All four groups show increasing trends on both the mean the variation level. Previous studies focused on building regression models on weight gain [14], i.e., the weight change, but no work has been done to directly capture the heteroscedasticity in GCA.

We fitted Model (3.2) to this dataset. The model includes fixed effects consisting of linear I-spline bases of time and their interactions with the diet, as well as a chick-level random effects on the slope of time. Heteroscedasticity is characterized by an I-spline of an index variable, which is either the marginal mean or the conditional mean of the linear effects model. Specifically, the model is y i j = β 0 + ∑ k = 1 K 1 ( β k I 1 k ( t i j ) + β k ( 2 ) I 1 k ( t i j ) D i ( 2 ) + β k ( 3 ) I 1 k ( t i j ) D i ( 3 ) + β k ( 4 ) I 1 k ( t i j ) D i ( 4 ) ) + b i t i j + ε i j , b i ∼ N ( 0 , σ b 2 ) , ε i j ∼ N 0 , θ 0 + ∑ k = 1 K 2 θ k I 2 k ( ν i j ) 2 , where y i j is the weight of chick i at time t i j , { I 1 k ( t i j ) : k = 1 , … , K 1 } is a set of I-splines bases with K 1  degrees of freedom used in the mean model, ( D i j 2 , D i j ( 3 ) , D i j ( 4 ) ) are the dummy variables of diet using diet 1 as the reference level, b i is a normally distributed chick-level random effect with variance σ b 2 , the error term ε i j is normal with mean zero and variance changing with index variable  ν i j (either marginal or conditional mean), { I 2 k ( μ i j ) : k = 1 , … , K 2 } is a set of I-spline bases with K 2 degrees of freedom used in the heteroscedasticity model, and the regression coefficients to be estimated are { β k , β k ( 2 ) , β k ( 3 ) , β k ( 4 ) : k = 1 , … , K 1 } and  { θ k : k = 1 , … , K 2 }. The internal knots of spline bases were always chosen to be evenly spaced.

We first need to decide whether to use the marginal mean or the conditional mean as the index variable in the heteroscedasticity model. For both situations, the BIC chose the same number of degrees of freedom of the I-splines. The I-spline bases for the mean pancreas length had degree 0 with 3 degrees of freedom. The I-spline bases for the variance model had degree 1 and 9 degrees of freedom. Since both the marginal mean model and the conditional mean model had the same number of parameters, we can choose the best model by only comparing the log-likelihood. The conditional mean model had a log-likelihood of − 2164.899, which is significantly higher than that of the marginal mean model, − 2199.183. Both of them were much higher than − 2315.220, the log-likelihood of the naive model that did not consider heteroscedasticity.

Figure 5 shows the diagnostic plots for the three models. For the naive model, the residual plot suggests increasing variance as fitted value increases and the Q-Q plot suggest heavier tail than the normal distribution. After considering heteroscedasticity with the proposed method, the model with variance changing with the marginal mean model still has heavy tail problem as seen from the Q-Q plot. The model with variance changing with the conditional mean, however, shows no obvious deviation from the normal distribution in the Q-Q plot of the standardized residuals. These diagnostics are consistent with the model comparison results in terms of log-likelihood.

Point-wise quantiles are of important practical value in GCA since they can be used as reference to check if an individual is in the normal range. For this linear mixed-effects model, we approximated the quantiles by generating 10,000 individuals using the random effects and for each of them, simulating their grow curve using the fitted model. The upper and lower 5% quantiles will give us 90% confidence interval, and the upper and lower 2.5% quantiles will form the 95% confidence interval. Finally, the fitted growth curves for the four diet groups along with their point-wise 90% and 95% confidence intervals from the model with residual variance changing with the conditional mean are shown in Figure 6. Also overlaid are the fitted individual curves for all the chicks in each diet group. It indicates that diet 3 works best for the chicks with fastest growing speed and highest weight at the end of the experiment, while diet 1 is the least favorable. The estimated results along with other information such as diet cost, can help the farmer decide which diet to use to generate the highest profit. The quantile estimation can also provide guidance on unhealthy chicken screening, and stop potential disease spead at early stage.