Nonlinear quantile mixed models
Abstract
In regression applications, the presence of nonlinearity and correlation among observations offer computational challenges not only in traditional settings such as least squares regression, but also (and especially) when the objective function is nonsmooth as in the case of quantile regression. In this paper, we develop methods for the modeling and estimation of nonlinear conditional quantile functions when data are clustered within twolevel nested designs. This work represents an extension of the linear quantile mixed models of Geraci and Bottai (2014, Statistics and Computing). We develop a novel algorithm which is a blend of a smoothing algorithm for quantile regression and a second order Laplacian approximation for nonlinear mixed models. To assess the proposed methods, we present a simulation study and two applications, one in pharmacokinetics and one related to growth curve modeling in agriculture.
Nonlinear quantile mixed models
t1Corresponding author: Marco Geraci, Department of Epidemiology and Biostatistics, Arnold School of Public Health, University of South Carolina, 915 Greene Street, Columbia SC 29209, USA. \printeade1
class=MSC] \kwd[Primary ]62F99 \kwd[; secondary ]62J99
asymmetric Laplace distribution \kwdconditional percentiles \kwdmultilevel designs \kwdrandom effects
1 Introduction
Quantile regression analysis of clustered data is a very active area of research. Since the seminal work of Koenker and Bassett (1978) on methods for crosssectional observations, there have been a number of proposals on how to accommodate for the dependency induced by clustered (e.g., longitudinal) designs. As outlined by Geraci and Bottai (2014) and then extensively reviewed by Marino and Farcomeni (2015), approaches to linear quantile regression with clustered data can be classified into distributionfree and likelihoodbased approaches. The former include fixed effects (Koenker, 2004; Lamarche, 2010; Galvao and MontesRojas, 2010; Galvao, 2011) and weighted (Lipsitz et al., 1997; Fu and Wang, 2012) approaches. The latter are mainly based on the asymmetric Laplace (AL) density (Geraci and Bottai, 2007, 2014; Yuan and Yin, 2010; Farcomeni, 2012) or other, usually flexible, parametric distributions (for example, Rigby and Stasinopoulos, 2005; Reich, Bondell and Wang, 2010; Noufaily and Jones, 2013). A different classification can be made into approaches that include clusterspecific effects (e.g., Koenker, 2004; Geraci and Bottai, 2014) and those that ignore or remove them (Lipsitz et al., 1997; Canay, 2011; Parente and Santos Silva, 2016).
In some applications, the assumption of linearity may not be appropriate. This is often the case in, for example, pharmacokinetics (Lindsey, 2001) and growth curve modeling (Panik, 2014). Contributions to statistical methods for nonlinear mean regression when data are clustered can be found in the literature of mixedeffects modeling (Lindstrom and Bates, 1990; Pinheiro and Bates, 1995, 2000) as well as generalized estimating equations (Davidian and Giltinan, 1995, 2003; Contreras and Ryan, 2000; Vonesh et al., 2002). We examined the statistical literature on parametric nonlinear quantile regression functions with clustered data (thus, in our review, we did not consider smoothing for nonparametric quantile functions). To the best of our knowledge, there seem to be only a handful of published articles. Karlsson (2008) considered nonlinear longitudinal data and proposed weighting the standard quantile regression estimator (Koenker and Bassett, 1978) with prespecified weights. Wang (2012), taking her cue from Geraci and Bottai (2007), used the AL distribution to define the likelihood of a Bayesian nonlinear quantile regression model. Huang and Chen (2016) proposed a Bayesian joint model for timetoevent and longitudinal data. An approach based on copulas is developed by Chen, Koenker and Xiao (2009). Finally, Oberhofer and Haupt (2016) established the consistency of the norm nonlinear quantile estimator under weak dependency. None of the cited papers provides an approach to modeling and estimation of nonlinear quantile functions with random effects in a frequentist framework. Such an approach is desirable when the correlation between measurements is modeled by means of random effects, but the parameters of interest are assumed to be fixed.
In this paper, we propose an extension of Geraci and Bottai’s (2014) linear quantile mixed model (LQMM) to the nonlinear case. In Section 2, we briefly outline the LQMM approach and, in Section 3, we introduce the nonlinear quantile mixedeffects model or nonlinear quantile mixed model (NLQMM) for short. Estimation is carried out using a novel algorithm which is a combination of a smoothing algorithm for quantile regression and a second order Laplacian approximation for nonlinear mixed models. We then carry out a simulation study in Section 4 to assess the performance of the proposed methods. In Section 5, we consider an application of NLQMM to pharmacokinetics and growth curves modeling. We conclude with some remarks in Section 6.
2 Linear quantile mixed models
Consider data from a twolevel nested design in the form , for and , , where is the th row of a known matrix , is the th row of a known matrix and is the th observation of the response vector for the th cluster. Throughout the paper, the covariates and are assumed to be given. The vector of ones will be denoted by , the identity matrix by , and the matrix of zeros by .
In a distributionfree approach, the linear quantile regression model for clustered (or panel) data (e.g., Koenker, 2004; Abrevaya and Dahl, 2008; Bache, Dahl and Kristensen, 2013) can be specified as
(1) 
where is the given quantile level, is a vector of coefficients common to all clusters, while the vector may vary with cluster. All the parameters are specific, although the clusterspecific effects are often specified simply as pure locationshift effects (Koenker, 2004; Lamarche, 2010). Fitting can be achieved by solving the classical norm regression problem
(2) 
where is the loss function and denotes the indicator function. The penalty on the clusterspecific effects controls the variability introduced by a large number of parameters and is usually based on the norm (Koenker, 2004; Lamarche, 2010; Bache, Dahl and Kristensen, 2013).
To mimic the minimization problem (2) in a likelihood framework, Geraci and Bottai (2014) introduced the convenient assumption that the responses , , , conditionally on a vector of random effects , independently follow the asymmetric Laplace (AL) density
with location and scale parameters given by and , respectively, which we write as . (The third parameter of the AL is the skew parameter which, in this model, is fixed and defines the quantile level of interest.) Also, they assumed that , for , is a random vector independent from the model’s error term with mean zero and variancecovariance matrix . Note that all the parameters are dependent. The random effects vector depends on through the variancecovariance matrix. If we let and , the joint density based on clusters for the th linear quantile mixed model (LQMM) is given by (Geraci and Bottai, 2014)
Geraci and Bottai (2014) proposed estimating LQMMs through a combination of Gaussian quadrature and nonsmooth optimization. They approximated the marginal (over the random effects) loglikelihood using the rule
(3) 
with , where and , , , denote, respectively, the abscissas and weights of the (onedimensional) Gaussian quadrature. In principle, one can consider different distributions for the random effects, which may be naturally linked to different quadrature rules (or penalties). For example, it is immediate to verify that the double exponential distribution confers robustness to the model and is akin to a GaussLaguerre quadrature (Geraci and Bottai, 2014). Alternatively, one can avoid parametric assumptions altogether (Alfò, Salvati and Ranallli, 2016). Throughout this paper we assume that the random effects are normally distributed and we do not explore this issue any further. However, this assumption can be relaxed and the developments that follow can be modified as appropriate.
3 Nonlinear quantile mixed models
3.1 The model
We consider the nonlinear quantile regression function
(4) 
where is a nonlinear, smooth function of the random parameter , and are two given design matrices of dimensions and , respectively, which in general contain elements of the covariates .
To stress the functional dependence of the quantiles on the fixed parameter and on the random parameter , we write . For estimation purposes, model (4) can be equivalently written as
(5) 
conditionally on , where . Moreover, we assume , independently from .
Note the similarities and dissimilarities between the proposed model (5) and the traditional nonlinear mixedeffects (NLME) model
with and . First of all, conditionally on the random effects, both models impose a restriction on the error term (Powell, 1994). However, the NLME model requires , while the ALbased specification of the error given in (5) leads to or, equivalently, . Secondly, the fixed effects can be interpreted as the average value of the clusterspecific parameters, i.e. , or as the regression parameters of the ‘zeromedian’ cluster, i.e. a cluster with a zero randomeffect vector. However, the parameter is allowed to vary with the quantile level , while in the NLME model is not (except for a location shift). Finally, in both approaches the variancecovariance matrix of the random effects gives a measure of the variability of around but, again, estimates are allowed to differ by only for the quantile mixedeffects model.
In general, neither model (5) nor the NLME model provide fixed parameters that can be interpreted as, respectively, regression quantiles or regression means for the population. This is because random effects are allowed to enter nonlinearly in the model. (Similarly, several generalized linear mixed models with nonlinear link functions lack marginal interpretability (Ritz and Spiegelman, 2004; Gory, Craigmile and MacEachern, 2016).) In contrast, the fixed effects of a linear (normal) mixed model remain the same after the random effects are integrated out, whereas, in general, this is not true for the fixed effects of the LQMMs of Geraci and Bottai (2014).
3.2 Laplacian approximation
Let be the scaled variancecovariance matrix of the random effects and define , where an an unrestricted dimensional vector, , of nonredundant parameters in . Our goal is to maximize the marginal loglikelihood
(6) 
where and .
First of all, we consider the following smooth approximation of (Madsen and Nielsen, 1993; Chen, 2007):
where and is a scalar “tuning” parameter. For , we have that . A similar approximation is given by Muggeo, Sciandra and Augugliaro (2012) who claimed that their method provides a better approximation than Chen’s (2007) algorithm. However, no analytical evidence was provided in their paper to support such claim. This point might offer scope for additional investigation but, here, it represents a secondary issue and will not be discussed any further.
Let be the vector of residuals , , for the th cluster, and define the corresponding sign vector with
(7) 
(Note that the notation above has been simplified since the ’s as well as the ’s should be written as functions of the ’s.) Then we have
(8) 
where is an diagonal matrix with diagonal elements , and are two vectors with elements
and
respectively.
We now define the function
(9) 
which is akin to a regularized, nonlinear, weighted leastsquares loss function. The gradient of with respect to is given by
(10) 
where and , while the Hessian is given by
(11) 
Moreover, let
(12) 
be the conditional mode of . For a given value of , this can be obtained by means of penalized leastsquares.
A secondorder approximation of around is given by
where , , and . Since is zero at , we have the following Laplacian approximation of the loglikelihood
(13) 
If we ignore the negligible contribution of the first term in expression (3.2) (Pinheiro and Bates, 1995), then only the firstorder partial derivatives of are required to compute (3.2).
Since does not depend on , the loglikelihood can be profiled on leading to
(14) 
where .
Estimation of the parameters can be carried out iteratively. A pseudocode of the algorithm is given in Appendix. The algorithm requires setting the starting value of and , the tuning parameter , the tolerance for the change in the loglikelihood, and the maximum number of iterations. Moreover, the modes of the random effects can be obtained by minimizing the objective function of the penalized leastsquares problem using a GaussNewton method. Let be the relative precision factor such that (Pinheiro and Bates, 2000). Then the function in (9) can be rewritten as
(15) 
where
When using the asymmetric Laplace as pseudolikelihood, inference must be restrained to point estimation (see for example Yang, Wang and He, 2016). Standard errors for nonrandom parameters can be calculated using block bootstrap, although this increases the computational cost. Bootstrap confidence intervals have been shown to have good coverage in LQMMs (Geraci and Bottai, 2014).
4 Simulation study
In this section, we perform a simulation study to evaluate the proposed methods. We start from a setting similar to the one used in Pinheiro and Bates (1995), which is ideal for normal NLME models, and then investigate scenarios more apposite for NLQMM.
In the first scenario, we simulated the data from the following logistic model
(16) 
where , , , and . The random effects are thus associated with the asymptotes ( and ) and the sigmoid’s midpoint (). Their variancecovariance matrix was defined as
In the second scenario, we used the same model (16), but we sampled the errors from a standardized chisquared distribution with degrees of freedom, i.e. .
In the third scenario, we slightly changed model (16) and used
(17) 
where , , , and . Note that the error is skewed as in the second scenario but now operates within the exponential function. In this heteroscedastic model, there is only one random effect associated with the sigmoid’s midpoint.
In the fourth and last scenario, we used the biexponential model
(18) 
where , , and , with parameters and .
In all scenarios, we used , , . Instances of replications are shown in Figure 1. For data sampled from models (16) and (17), we fitted mixedeffects logistic quantile functions with parameter , where
in the first 3 scenarios,
in the first and second scenarios, and in the third scenario. For data sampled from model (4), we fitted mixedeffects biexponential quantile functions with parameter , where .
For each scenario, we replicated datasets and fitted NLQMMs at three quantile levels using . Estimation was carried out by following the algorithm as described in Appendix. An attempt to maximize the approximated Laplacian loglikelihood (3.2) was made by using the BroydenFletcherGoldfarbShanno (BFGS) algorithm in the first instance. Upon failure of the BFGS algorithm during any iteration of the main estimation algorithm, the latter was started again and a new attempt to maximize (3.2) was made by using the NelderMead algorithm. The maximum number of iterations was set to 500, while the tolerance for the relative change in the loglikelihood was set to . Between two successive iterations, the tuning parameter was multiplied by the factor . All the parameters of the optimization algorithms in optim and nlm were set at their default values. Computations were performed using the R environment for statistical computing and graphics (R Core Team, 2016) version 3.3.2 on a desktop computer with a 3.60GHz quad core i74790 processor and 16 gigabytes of RAM.
Before we proceed with the analysis of the results, it is important to note that, in general, the nonlinearity of the models along with the presence of the random effects pose a difficulty for establishing the ‘true’ value of for quantiles other than the median (see for example the simulation study in Karlsson, 2008), even when the errors are normal. For example, in the logistic model not only the asymptotes and , but also the midpoint and the scale change with in a rather complicated way. (An exception is given by model (17) for which the lower and upper asymptotes ( and , respectively) do not change with .) We find solace in observing that such limitation brings out one of the advantages of quantilebased over momentbased modeling, since direct estimation of conditional quantiles does not require nontrivial manipulation of nonlinear relationships (Demidenko, 2013, p.435). As a reference, we can consider the corresponding results from standard nonlinear quantile regression (NLRQ) (Koenker and Park, 1996) under the assumption of independent observations. Similarity of the magnitude and direction of the estimates would support the interpretation of as regression parameters of the ‘zeromedian’ cluster, while comparing the variability of the estimates from NLQMM and NLRQ would inform us on whether accounting for clustering provides a gain in efficiency.
The average estimates and standard deviations of the estimates are reported in Tables 14. In summary, NLQMM estimates were close to NLRQ estimates in all scenarios. The variability of the estimates from NLQMM was either lower or close to that of the estimates from NLQR. Of all the results, perhaps those related to the quantile in the third scenario (Table 3) deserve more discussion. Both NLQMM and NLRQ clearly failed to provide meaningful estimates of the parameters. This is due to the fact that the range of the simulated values for was not wide enough to correctly estimate the upper asymptote at upper quantiles. This observation may have a particular relevance when modeling reference growth curves. Further, the estimated variancecovariance parameters (Table 5) and the predicted random effects obtained from (12) (Figures 25) were, in general, consistent with the parameters of the true distribution of the random effects, although, as noted before, direct comparisons are not straightforward.
We now provide basic details about the algorithm’s performance and a few recommendations. On average, it took about 7 iterations (approximately 35 seconds) to fit one model for the quantiles or , and about 6 iterations (approximately 20 seconds) for the median in the first two scenarios. In the third scenario it took between and iterations (approximately seconds on average) to fit one model for any of the three quantile levels. In the fourth scenario, the algorithm needed a similar number of iterations as in the first two scenarios but the time to convergence was, on average, twice as long. This means that, within each iteration, the number of function evaluations required by optim to fit the more complex biexponential model was greater than that needed to fit the logistic model. In the first two scenarios, the median value of the factor at the last iteration was about for all considered quantiles. In the third scenario, it was less than for the tail quantiles and for the median. In the fourth scenario, it was less than for all considered quantiles.
In a separate analysis (results not shown), the average number of iterations to convergence increased to at least 10 when was increased to . In contrast, the algorithm converged too quickly to smaller values of the loglikelihood when setting to less than 0.2. This was to be expected since controls the speed at which the smoothing parameter approaches zero. As in Chen (2007), we recommend using in most situations.
Further, in the first three scenarios the average number of iterations and the values of the estimates were not particularly sensitive to the specific algorithm used for optimizing the loglikelihood, although the BFGS algorithm did fail to converge in about of the replications, more often when estimating tail quantiles () rather than when estimating the median (). In contrast, BFGS never failed to converge in the fourth scenario. We then ran a separate analysis (results not shown) in which biexponential models were fitted exclusively using NelderMead. For , estimates were unreasonable. We recommend using BFGS as default optimization algorithm.
NLQMM  

68.00 (0.73)  12.66 (0.40)  3.06 (0.10)  8.79 (0.34)  
70.23 (0.38)  9.99 (0.28)  3.04 (0.05)  9.70 (0.22)  
73.47 (0.56)  7.28 (0.38)  3.15 (0.10)  9.76 (0.52)  
NLRQ  
68.47 (2.76)  12.77 (0.58)  3.04 (0.27)  9.54 (0.53)  
69.79 (0.92)  9.99 (0.33)  3.00 (0.18)  10.11 (0.71)  
72.26 (0.85)  7.19 (0.53)  3.11 (0.32)  9.53 (2.75) 
NLQMM  

69.28 (0.68)  12.70 (0.40)  3.05 (0.08)  10.25 (0.22)  
71.39 (0.35)  9.97 (0.27)  3.05 (0.05)  10.64 (0.19)  
74.67 (0.55)  7.26 (0.38)  3.15 (0.10)  10.98 (0.58)  
NLRQ  
69.67 (2.53)  12.77 (0.55)  3.03 (0.24)  10.80 (0.46)  
71.03 (0.90)  9.98 (0.31)  3.01 (0.18)  11.25 (0.72)  
73.49 (0.82)  7.17 (0.54)  3.11 (0.32)  10.73 (2.86) 
NLQMM  

0.95 (0.06)  4.18 (0.20)  0.92 (0.08)  0.00 (0.00)  
1.00 (0.03)  3.19 (0.08)  0.85 (0.05)  0.00 (0.01)  
0.25 (0.93)  2.09 (1.34)  1.16 (1.40)  1.11 (0.36)  
NLRQ  
1.01 (0.03)  3.96 (0.06)  1.01 (0.02)  0.00 (0.00)  
1.00 (0.04)  3.18 (0.10)  0.87 (0.06)  0.00 (0.01)  
0.48 (0.57)  1.24 (1.65)  2.19 (0.98)  0.40 (0.22) 
NLQMM  

1.93 (0.19)  1.06 (0.12)  0.46 (0.13)  ()  
2.04 (0.15)  0.72 (0.10)  1.03 (0.11)  3.19 (0.17)  
2.05 (0.17)  0.61 (0.15)  1.82 (0.12)  2.48 (0.11)  
NLRQ  
1.86 (0.26)  0.98 (0.19)  0.61 (0.14)  ()  
2.04 (0.18)  0.70 (0.14)  1.01 (0.13)  3.27 (0.24)  
2.17 (0.21)  0.61 (0.19)  1.55 (0.14)  2.32 (1.89) 
First scenario  

1.73  0.97  1.68  
3.30  1.83  3.19  
6.61  3.67  6.36  
Second scenario  
1.76  0.99  1.73  
3.34  1.85  3.25  
6.62  3.65  6.38  
Third scenario  
0.01  0.04  0.00  
Fourth scenario  
0.16  0.22  0.04  
0.24  0.25  0.08  
0.16  0.18  0.17  
0.09  0.17  0.01 
5 Applications
5.1 Pharmacokinetics
We begin with the analysis of a dataset taken from an old pharmacokinetics study (Kwan et al., 1976), often used as a toy example in nonlinear regression modeling (Davidian and Giltinan, 1995; Pinheiro and Bates, 2000). The data consists of 11 measurements of plasma concentrations of indomethacin which was injected intravenously in 6 subjects. We used the biexponential model
where denotes the th measurement of drug concentration (mcg/ml), , on the th subject, , and is the time (hr) of the measurement since injection (given that the dataset is balanced, does not depend on ). We modeled the variancecovariance of the random effects using the diagonal matrix (variance components). Note that the regression model above includes 3 random effects, one for each of the first 3 fixed effects. In a separate analysis (results not shown), we found that the random effect associated with , , had nearzero variance (see also Pinheiro and Bates, 2000, p.283).
In this twocompartment model, the first exponential term determines the initial, rapidly declining distribution phase of the drug. The elimination of the drug is the predominant process during the second phase and is primarily determined by the second exponential term. Besides the average rates at which the drug is distributed and then eliminated, it might be of interest to assess the change over time of concentrations that are higher or lower than the mean. The left plot of Figure 6 shows the boxplots of indomethacin concentration at each measurement occasion. It appears that the scale and possibly even the shape of the distribution are changing over time. This would mean that the rates are not similar across the quantiles of the conditional distribution.
2.31 (0.48)  0.99 (0.16)  0.30 (0.13)  1.19 (0.57)  
2.55 (0.28)  0.58 (0.19)  0.44 (0.17)  1.33 (0.23)  
3.73 (0.52)  0.75 (0.35)  0.69 (0.34)  1.49 (0.37)  
NLME  2.83 (0.26)  0.77 (0.11)  0.46 (0.11)  1.35 (0.23) 
0.78  0.06  0.02  
0.59  0.08  0.02  
1.34  0.05  0.06  
NLME  0.33  0.03  0.01 
The fitted biexponential curves for are given in the right plot of Figure 6, while estimates of the fixed effects and their standard errors are reported in Table 6. The average rate (NLME) at which the drug decreases during the distribution phase was comparable to that of the 90th centile. However, the decrease was about faster at the lower 10th centile but about slower at the median as compared to the mean. During the second phase, the rate of decrease was, again, greatest at the 10th centile. However, the average rate was similar to that of the median and greater than that of the 90th centile. One implication is that the distribution of the response becomes increasingly rightskewed as time passes. This is an advantage of NLQMM over NLME as the latter cannot obviously model changes in the shape of the distribution.
Finally, there was substantial heterogeneity among subjects regarding starting concentration levels during the distribution phase, especially at the 90th centile (Table 7).
5.2 Growth curves
In this section, we analyze data on growth patterns of two genotypes of soybeans: Plant Introduction (P), an experimental strain, and Forrest (F), a commercial variety (Davidian and Giltinan, 1995). The response variable is the average leaf weight of 6 plants chosen at random from each experimental plot and measured at approximately weekly intervals, between two and eleven weeks after planting. The experiment was carried out over three different planting years: 1988, 1989, and 1990. Eight plots were planted with each genotype in each planting year, giving a total of 48 plots in the study (Pinheiro and Bates, 2000).
Figure 7 shows the temporal trajectories of the average leaf weight for individual plots. It is apparent that the experimental strain yielded heavier leaves that the F variety, at least on average. There also seem to be differences between planting years, with a wider spread of the curves in 1989. Given that previous analyses of these data focused on the average growth curves (Davidian and Giltinan, 1995; Pinheiro and Bates, 2000), we set out to investigate growth in the tails of the distribution. For our analysis, we used the same logistic model as that in Pinheiro and Bates (2000, p.293) which was selected over a number of alternative models, that is
where denotes the average leaf weight (g) observed on occasion , , in the th plot, , and is the time (days) of the measurement after planting. The number of measurements per plot ranged between 8 and 10. The design matrices of the parameter were given by
and . Thus, is a vector. The covariates in the matrix are dummy variables for year of planting, and , and genotype, . The baseline is represented by year 1988 and variety F. The only random effect included in the model was associated with the asymptote.
The three plots in Figure 8 show the th centile, th centile, and mean predicted growth curves by variety and planting year, while the estimates and standard errors of the fixed effects are reported in Table 8. For the sake of brevity, we confine our discussion to the genotypic effect on the asymptote. In 1988, the experimental strain had an advantage over the commercial variety but only at the 95th centile of the leaf weight distribution, with an estimated asymptote difference of g. In the following year, there was a statistically significant interaction between variety and year at the 5th centile, corresponding to an estimated overall effect equal to g. The estimated overall effect of variety P on the asymptote at the 95th centile was 10.67 g, thus greater than the effect at the 5th centile and at the mean (7.19 g). Finally, the estimated differences between asymptotes of the growth curves for the experimental and commercial strains in 1990 were negligible: 0.58 g (5th centile), 2.81 g (95th centile) and 0.77 g (mean). In summary, the experimental strain did yield heavier leaves than the F variety, not only in 1989 as estimated by NLME, but also in 1988, and the magnitude of the genotypic effect was dependent on the quantile level.
NLME  
17.49 (1.47)  21.43 (2.34)  19.43 (0.95)  
7.99 (1.53)  7.02 (2.30)  8.84 (1.07)  
0.66 (2.06)  1.67 (2.49)  3.71 (1.18)  
1.64 (2.01)  6.31 (1.99)  1.62 (1.04)  
8.59 (1.93)  4.36 (2.41)  5.57 (1.17)  
2.22 (2.05)  3.50 (2.01)  0.15 (1.18)  
56.16 (1.13)  53.71 (2.57)  54.81 (0.75)  
3.30 (2.11)  0.86 (2.85)  2.24 (0.97)  
1.94 (2.48)  3.14 (2.79)  4.97 (0.97)  
2.50 (1.70)  0.51 (0.97)  1.30 (0.41)  
8.11 (0.32)  8.63 (0.79)  8.06 (0.15)  
0.29 (0.51)  0.76 (0.85)  0.90 (0.20)  
0.40 (0.49)  0.44 (0.91)  0.67 (0.21) 
6 Discussion
Mixedeffects modeling has a long tradition in statistical applications. There is a vast number of applications of mixed models to the analysis of clustered data in the social, life and physical sciences (Pinheiro and Bates, 2000; Demidenko, 2013). Linear quantile mixed models (Geraci and Bottai, 2007, 2014) have too been used in a wide range of research areas, including marine biology (Muir et al., 2015; Duffy et al., 2015; Barneche et al., 2016), environmental science (Fornaroli et al., 2015), cardiovascular disease (Degerud et al., 2014; Blankenberg et al., 2016), physical activity (Ng et al., 2014; Beets et al., 2016), and ophthalmology (Patel et al., 2015; Patel, Geraci and CortinaBorja, 2016). The present paper provides a novel and valuable contribution to the modeling of nonlinear quantile functions which broadens the applicability of quantile regression for clustered data.
NLQMMs represent a flexible alternative to nonlinear mixed models for the mean as they allow direct estimation of conditional quantile functions without imposing normal assumptions on the errors. As shown in two real data examples, NLQMMs reveal nonlinear relationships that may be quantitatively and qualitatively different at different quantiles. Also, changes in location, scale, and shape of the response distribution determined by the covariates are naturally brought into light by examining central and tail quantiles (Geraci, 2016). As compared to nonlinear quantile regression for independent data, our nonlinear estimators are more efficient and they provide additional information about the heterogeneity among clusters.
Acknowledgements
This research has been supported by the National Institutes of Health – National Institute of Child Health and Human Development (Grant Number: 1R03HD08480701A1).
Appendix
The estimation algorithm for NLQMM is based on a set of decreasing values of . This optimization approach has the appealing advantage of reducing the original problem to an approximated problem. The pseudocode is given below.
Smoothing Algorithm with Laplacian Approximation for Nonlinear Quantile Mixed Models

Set the maximum number of iterations ; the factor for reducing the tuning parameter ; the tolerance for the change in the loglikelihood; and . Estimate the starting values as follows:

obtain an estimate for using nonlinear quantile regression (Koenker and Park, 1996). See, for example, the function nlrq in the R package quantreg (Koenker, 2016) which supports selfstarting models such as SSlogis. If the nonlinear quantile regression algorithm fails, consider the estimate of the fixed effects from the NLME model in step (1.b) below or, if the latter fails too, a standard nonlinear least squares estimate (Bates and Watts, 1988);

obtain an estimate for from an NLME model. See, for example, the function nlme in the R package (Pinheiro et al., 2014). If the NLME algorithm fails, provide an arbitrary value ;

obtain an estimate for . For example, this can be estimated as the mean of the absolute residuals from step (1.a) above;

provide a starting value (see, for example, Chen, 2007, p.143);

using , , and , solve the penalized leastsquares problem (12) to obtain , . See, for example, the R function nlm.


Update and , .
References
 Abrevaya and Dahl (2008) {barticle}[author] \bauthor\bsnmAbrevaya, \bfnmJ.\binitsJ. \AND\bauthor\bsnmDahl, \bfnmC. M.\binitsC. M. (\byear2008). \btitleThe effects of birth inputs on birthweight: Evidence from quantile estimation on panel data. \bjournalJournal of Business and Economic Statistics \bvolume26 \bpages379397. \endbibitem
 Alfò, Salvati and Ranallli (2016) {barticle}[author] \bauthor\bsnmAlfò, \bfnmMarco\binitsM., \bauthor\bsnmSalvati, \bfnmNicola\binitsN. \AND\bauthor\bsnmRanallli, \bfnmM. Giovanna\binitsM. G. (\byear2016). \btitleFinite mixtures of quantile and Mquantile regression models. \bjournalStatistics and Computing \bvolume27 \bpages547570. \endbibitem
 Bache, Dahl and Kristensen (2013) {barticle}[author] \bauthor\bsnmBache, \bfnmStefanHolstMilton\binitsS., \bauthor\bsnmDahl, \bfnmChristianMÃ¸ller\binitsC. \AND\bauthor\bsnmKristensen, \bfnmJohannesTang\binitsJ. (\byear2013). \btitleHeadlights on tobacco road to low birthweight outcomes. \bjournalEmpirical Economics \bvolume44 \bpages15931633. \endbibitem
 Barneche et al. (2016) {barticle}[author] \bauthor\bsnmBarneche, \bfnmD. R.\binitsD. R., \bauthor\bsnmKulbicki, \bfnmM.\binitsM., \bauthor\bsnmFloeter, \bfnmS. R.\binitsS. R., \bauthor\bsnmFriedlander, \bfnmA. M.\binitsA. M. \AND\bauthor\bsnmAllen, \bfnmA. P.\binitsA. P. (\byear2016). \btitleEnergetic and ecological constraints on population density of reef fishes. \bjournalProceedings of the Royal Society B: Biological Sciences \bvolume283 \bpages18. \endbibitem
 Bates and Watts (1988) {bbook}[author] \bauthor\bsnmBates, \bfnmDouglas M\binitsD. M. \AND\bauthor\bsnmWatts, \bfnmDonald G\binitsD. G. (\byear1988). \btitleNonlinear regression analysis and its applications. \bpublisherWiley, \baddressHoboken, NJ. \endbibitem
 Beets et al. (2016) {barticle}[author] \bauthor\bsnmBeets, \bfnmMichael W.\binitsM. W., \bauthor\bsnmWeaver, \bfnmRobert Glenn\binitsR. G., \bauthor\bsnmTurnerMcGrievy, \bfnmGabrielle\binitsG., \bauthor\bsnmMoore, \bfnmJustin B.\binitsJ. B., \bauthor\bsnmWebster, \bfnmCollin\binitsC., \bauthor\bsnmBrazendale, \bfnmKeith\binitsK. \betalet al. (\byear2016). \btitleAre we there yet? Compliance with physical activity standards in YMCA afterschool programs. \bjournalChildhood Obesity \bvolume12 \bpages237246. \endbibitem
 Blankenberg et al. (2016) {barticle}[author] \bauthor\bsnmBlankenberg, \bfnmStefan\binitsS., \bauthor\bsnmSalomaa, \bfnmVeikko\binitsV., \bauthor\bsnmMakarova, \bfnmNataliya\binitsN., \bauthor\bsnmOjeda, \bfnmFrancisco\binitsF., \bauthor\bsnmWild, \bfnmPhilipp\binitsP., \bauthor\bsnmLackner, \bfnmKarl J.\binitsK. J. \betalet al. (\byear2016). \btitleTroponin I and cardiovascular risk prediction in the general population: the BiomarCaRE consortium. \bjournalEuropean Heart Journal \bvolume37 \bpages24282437. \endbibitem
 Canay (2011) {barticle}[author] \bauthor\bsnmCanay, \bfnmI. A.\binitsI. A. (\byear2011). \btitleA simple approach to quantile regression for panel data. \bjournalEconometrics Journal \bvolume14 \bpages368386. \endbibitem
 Chen (2007) {barticle}[author] \bauthor\bsnmChen, \bfnmColin\binitsC. (\byear2007). \btitleA finite smoothing algorithm for quantile regression. \bjournalJournal of Computational and Graphical Statistics \bvolume16 \bpages136164. \endbibitem
 Chen, Koenker and Xiao (2009) {barticle}[author] \bauthor\bsnmChen, \bfnmXiaohong\binitsX., \bauthor\bsnmKoenker, \bfnmRoger\binitsR. \AND\bauthor\bsnmXiao, \bfnmZhijie\binitsZ. (\byear2009). \btitleCopulabased nonlinear quantile autoregression. \bjournalEconometrics Journal \bvolume12 \bpagesS50S67. \endbibitem
 Contreras and Ryan (2000) {barticle}[author] \bauthor\bsnmContreras, \bfnmMartha\binitsM. \AND\bauthor\bsnmRyan, \bfnmLouise M.\binitsL. M. (\byear2000). \btitleFitting nonlinear and constrained generalized estimating equations with optimization software. \bjournalBiometrics \bvolume56 \bpages12681271. \endbibitem
 Davidian and Giltinan (1995) {bbook}[author] \bauthor\bsnmDavidian, \bfnmM.\binitsM. \AND\bauthor\bsnmGiltinan, \bfnmD. M.\binitsD. M. (\byear1995). \btitleNonlinear models for repeated measurement data. \bpublisherChapman and Hall/CRC, \baddressBoca Raton, FL. \endbibitem
 Davidian and Giltinan (2003) {barticle}[author] \bauthor\bsnmDavidian, \bfnmMarie\binitsM. \AND\bauthor\bsnmGiltinan, \bfnmDavid M.\binitsD. M. (\byear2003). \btitleNonlinear models for repeated measurement data: An overview and update. \bjournalJournal of Agricultural, Biological, and Environmental Statistics \bvolume8 \bpages387. \endbibitem
 Degerud et al. (2014) {barticle}[author] \bauthor\bsnmDegerud, \bfnmEirik\binitsE., \bauthor\bsnmLøland, \bfnmKjetil H\binitsK. H., \bauthor\bsnmNygård, \bfnmOttar\binitsO., \bauthor\bsnmMidttun, \bfnmØivind\binitsØ., \bauthor\bsnmUeland, \bfnmPer M\binitsP. M., \bauthor\bsnmSeifert, \bfnmReinhard\binitsR. \betalet al. (\byear2014). \btitleVitamin D status was not associated with ‘oneyear’ progression of coronary artery disease, assessed by coronary angiography in statintreated patients. \bjournalEuropean Journal of Preventive Cardiology \bvolume22 \bpages594602. \endbibitem
 Demidenko (2013) {bbook}[author] \bauthor\bsnmDemidenko, \bfnmE.\binitsE. (\byear2013). \btitleMixed models: Theory and applications with R, \beditionSecond ed. \bpublisherJohn Wiley & Sons, \baddressHoboken, NJ. \endbibitem
 Duffy et al. (2015) {barticle}[author] \bauthor\bsnmDuffy, \bfnmL. M.\binitsL. M., \bauthor\bsnmOlson, \bfnmR. J.\binitsR. J., \bauthor\bsnmLennertCody, \bfnmC. E.\binitsC. E., \bauthor\bsnmGalvánMagaña, \bfnmF.\binitsF., \bauthor\bsnmBocanegraCastillo, \bfnmN.\binitsN. \AND\bauthor\bsnmKuhnert, \bfnmP. M.\binitsP. M. (\byear2015). \btitleForaging ecology of silky sharks, Carcharhinus falciformis, captured by the tuna purseseine fishery in the eastern Pacific Ocean. \bjournalMarine Biology \bvolume162 \bpages571593. \endbibitem
 Farcomeni (2012) {barticle}[author] \bauthor\bsnmFarcomeni, \bfnmAlessio\binitsA. (\byear2012). \btitleQuantile regression for longitudinal data based on latent Markov subjectspecific parameters. \bjournalStatistics and Computing \bvolume22 \bpages141152. \endbibitem
 Fornaroli et al. (2015) {barticle}[author] \bauthor\bsnmFornaroli, \bfnmR.\binitsR., \bauthor\bsnmCabrini, \bfnmR.\binitsR., \bauthor\bsnmSartori, \bfnmL.\binitsL., \bauthor\bsnmMarazzi, \bfnmF.\binitsF., \bauthor\bsnmVracevic, \bfnmD.\binitsD., \bauthor\bsnmMezzanotte, \bfnmV.\binitsV. \betalet al. (\byear2015). \btitlePredicting the constraint effect of environmental characteristics on macroinvertebrate density and diversity using quantile regression mixed model. \bjournalHydrobiologia \bvolume742 \bpages153167. \endbibitem
 Fu and Wang (2012) {barticle}[author] \bauthor\bsnmFu, \bfnmLiya\binitsL. \AND\bauthor\bsnmWang, \bfnmYouGan\binitsY.G. (\byear2012). \btitleQuantile regression for longitudinal data with a working correlation model. \bjournalComputational Statistics & Data Analysis \bvolume56 \bpages25262538. \endbibitem
 Galvao (2011) {barticle}[author] \bauthor\bsnmGalvao, \bfnmA. F.\binitsA. F. (\byear2011). \btitleQuantile regression for dynamic panel data with fixed effects. \bjournalJournal of Econometrics \bvolume164 \bpages142157. \endbibitem
 Galvao and MontesRojas (2010) {barticle}[author] \bauthor\bsnmGalvao, \bfnmA. F.\binitsA. F. \AND\bauthor\bsnmMontesRojas, \bfnmG. V.\binitsG. V. (\byear2010). \btitlePenalized quantile regression for dynamic panel data. \bjournalJournal of Statistical Planning and Inference \bvolume140 \bpages34763497. \endbibitem
 Geraci (2016) {barticle}[author] \bauthor\bsnmGeraci, \bfnmM.\binitsM. (\byear2016). \btitleQtools: A collection of models and other tools for quantile inference. \bjournalR Journal \bvolume8 \bpages117138. \endbibitem
 Geraci and Bottai (2007) {barticle}[author] \bauthor\bsnmGeraci, \bfnmM.\binitsM. \AND\bauthor\bsnmBottai, \bfnmM.\binitsM. (\byear2007). \btitleQuantile regression for longitudinal data using the asymmetric Laplace distribution. \bjournalBiostatistics \bvolume8 \bpages140154. \endbibitem
 Geraci and Bottai (2014) {barticle}[author] \bauthor\bsnmGeraci, \bfnmM.\binitsM. \AND\bauthor\bsnmBottai, \bfnmM.\binitsM. (\byear2014). \btitleLinear quantile mixed models. \bjournalStatistics and Computing \bvolume24 \bpages461479. \endbibitem
 Gory, Craigmile and MacEachern (2016) {barticle}[author] \bauthor\bsnmGory, \bfnmJ. J.\binitsJ. J., \bauthor\bsnmCraigmile, \bfnmP. F.\binitsP. F. \AND\bauthor\bsnmMacEachern, \bfnmS. N.\binitsS. N. (\byear2016). \btitleMarginally interpretable generalized linear mixed models. \bjournalArXiv eprints \bvolume1610.01526. \endbibitem
 Huang and Chen (2016) {barticle}[author] \bauthor\bsnmHuang, \bfnmYangxin\binitsY. \AND\bauthor\bsnmChen, \bfnmJiaqing\binitsJ. (\byear2016). \btitleBayesian quantile regressionbased nonlinear mixedeffects joint models for timetoevent and longitudinal data with multiple features. \bjournalStatistics in Medicine \bvolume35 \bpages56665685. \endbibitem
 Karlsson (2008) {barticle}[author] \bauthor\bsnmKarlsson, \bfnmA.\binitsA. (\byear2008). \btitleNonlinear quantile regression estimation of longitudinal data. \bjournalCommunications in StatisticsSimulation and Computation \bvolume37 \bpages114131. \endbibitem
 Koenker (2004) {barticle}[author] \bauthor\bsnmKoenker, \bfnmR.\binitsR. (\byear2004). \btitleQuantile regression for longitudinal data. \bjournalJournal of Multivariate Analysis \bvolume91 \bpages7489. \endbibitem
 Koenker (2016) {bmanual}[author] \bauthor\bsnmKoenker, \bfnmR.\binitsR. (\byear2016). \btitlequantreg: Quantile Regression \bnoteR package version 5.29. \endbibitem
 Koenker and Bassett (1978) {barticle}[author] \bauthor\bsnmKoenker, \bfnmR.\binitsR. \AND\bauthor\bsnmBassett, \bfnmG.\binitsG. (\byear1978). \btitleRegression quantiles. \bjournalEconometrica \bvolume46 \bpages3350. \endbibitem
 Koenker and Park (1996) {barticle}[author] \bauthor\bsnmKoenker, \bfnmR.\binitsR. \AND\bauthor\bsnmPark, \bfnmB. J.\binitsB. J. (\byear1996). \btitleAn interior point algorithm for nonlinear quantile regression. \bjournalJournal of Econometrics \bvolume71 \bpages265283. \endbibitem
 Kwan et al. (1976) {barticle}[author] \bauthor\bsnmKwan, \bfnmK. C.\binitsK. C., \bauthor\bsnmBreault, \bfnmG. O.\binitsG. O., \bauthor\bsnmUmbenhauer, \bfnmE. R.\binitsE. R., \bauthor\bsnmMcMahon, \bfnmF. G.\binitsF. G. \AND\bauthor\bsnmDuggan, \bfnmD. E.\binitsD. E. (\byear1976). \btitleKinetics of indomethacin absorption, elimination, and enterohepatic circulation in man. \bjournalJournal of Pharmacokinetics and Biopharmaceutics \bvolume4 \bpages25580. \endbibitem
 Lamarche (2010) {barticle}[author] \bauthor\bsnmLamarche, \bfnmC\binitsC. (\byear2010). \btitleRobust penalized quantile regression estimation for panel data. \bjournalJournal of Econometrics \bvolume157 \bpages396498. \endbibitem
 Lindsey (2001) {bbook}[author] \bauthor\bsnmLindsey, \bfnmJ. K.\binitsJ. K. (\byear2001). \btitleNonlinear models for medical statistics. \bpublisherOxford University Press, \baddressNew York, NY. \endbibitem
 Lindstrom and Bates (1990) {barticle}[author] \bauthor\bsnmLindstrom, \bfnmM. J.\binitsM. J. \AND\bauthor\bsnmBates, \bfnmD. M.\binitsD. M. (\byear1990). \btitleNonlinear mixed effects models for repeated measures data. \bjournalBiometrics \bvolume46 \bpages673687. \endbibitem
 Lipsitz et al. (1997) {barticle}[author] \bauthor\bsnmLipsitz, \bfnmS. R.\binitsS. R., \bauthor\bsnmFitzmaurice, \bfnmG. M.\binitsG. M., \bauthor\bsnmMolenberghs, \bfnmG.\binitsG. \AND\bauthor\bsnmZhao, \bfnmL. P.\binitsL. P. (\byear1997). \btitleQuantile regression methods for longitudinal data with dropouts: Application to CD4 cell counts of patients infected with the human immunodeficiency virus. \bjournalJournal of the Royal Statistical Society C \bvolume46 \bpages463476. \endbibitem
 Madsen and Nielsen (1993) {barticle}[author] \bauthor\bsnmMadsen, \bfnmK.\binitsK. \AND\bauthor\bsnmNielsen, \bfnmH. B.\binitsH. B. (\byear1993). \btitleA finite smoothing algorithm for linear estimation. \bjournalSiam Journal on Optimization \bvolume3 \bpages223235. \endbibitem
 Marino and Farcomeni (2015) {barticle}[author] \bauthor\bsnmMarino, \bfnmMaria Francesca\binitsM. F. \AND\bauthor\bsnmFarcomeni, \bfnmAlessio\binitsA. (\byear2015). \btitleLinear quantile regression models for longitudinal experiments: an overview. \bjournalMETRON \bvolume73 \bpages229247. \endbibitem
 Muggeo, Sciandra and Augugliaro (2012) {barticle}[author] \bauthor\bsnmMuggeo, \bfnmV. M. R.\binitsV. M. R., \bauthor\bsnmSciandra, \bfnmM.\binitsM. \AND\bauthor\bsnmAugugliaro, \bfnmL.\binitsL. (\byear2012). \btitleQuantile regression via iterative least squares computations. \bjournalJournal of Statistical Computation and Simulation \bvolume82 \bpages15571569. \endbibitem
 Muir et al. (2015) {barticle}[author] \bauthor\bsnmMuir, \bfnmP. R.\binitsP. R., \bauthor\bsnmWallace, \bfnmC. C.\binitsC. C., \bauthor\bsnmDone, \bfnmT.\binitsT. \AND\bauthor\bsnmAguirre, \bfnmJ. D.\binitsJ. D. (\byear2015). \btitleLimited scope for latitudinal extension of reef corals. \bjournalScience \bvolume348 \bpages11351138. \endbibitem
 Ng et al. (2014) {barticle}[author] \bauthor\bsnmNg, \bfnmS. W.\binitsS. W., \bauthor\bsnmHoward, \bfnmA. G.\binitsA. G., \bauthor\bsnmWang, \bfnmH. J.\binitsH. J., \bauthor\bsnmSu, \bfnmC.\binitsC. \AND\bauthor\bsnmZhang, \bfnmB.\binitsB. (\byear2014). \btitleThe physical activity transition among adults in China: 19912011. \bjournalObesity Reviews \bvolume15 \bpages2736. \endbibitem
 Noufaily and Jones (2013) {barticle}[author] \bauthor\bsnmNoufaily, \bfnmAngela\binitsA. \AND\bauthor\bsnmJones, \bfnmM. C.\binitsM. C. (\byear2013). \btitleParametric quantile regression based on the generalized gamma distribution. \bjournalJournal of the Royal Statistical Society C \bvolume62 \bpages723740. \endbibitem
 Oberhofer and Haupt (2016) {barticle}[author] \bauthor\bsnmOberhofer, \bfnmWalter\binitsW. \AND\bauthor\bsnmHaupt, \bfnmHarry\binitsH. (\byear2016). \btitleAsymptotic theory for nonlinear quantile regression under weak dependence. \bjournalEconometric Theory \bvolume32 \bpages686713. \endbibitem
 Panik (2014) {bbook}[author] \bauthor\bsnmPanik, \bfnmM. J.\binitsM. J. (\byear2014). \btitleGrowth curve modeling: Theory and applications. \bpublisherJohn Wiley and Sons, \baddressHoboken, NJ. \endbibitem
 Parente and Santos Silva (2016) {barticle}[author] \bauthor\bsnmParente, \bfnmP. M. D. C.\binitsP. M. D. C. \AND\bauthor\bsnmSantos Silva, \bfnmJ. M. C.\binitsJ. M. C. (\byear2016). \btitleQuantile regression with clustered data. \bjournalJournal of Econometric Methods \bvolume5 \bpages116. \endbibitem
 Patel, Geraci and CortinaBorja (2016) {barticle}[author] \bauthor\bsnmPatel, \bfnmD. E.\binitsD. E., \bauthor\bsnmGeraci, \bfnmM.\binitsM. \AND\bauthor\bsnmCortinaBorja, \bfnmM.\binitsM. (\byear2016). \btitleModelling normative kinetic perimetry isopters using mixedeffects quantile regression. \bjournalJournal of Vision \bvolume16 \bpages16. \endbibitem
 Patel et al. (2015) {barticle}[author] \bauthor\bsnmPatel, \bfnmD. E.\binitsD. E., \bauthor\bsnmCumberland, \bfnmP. M.\binitsP. M., \bauthor\bsnmWalters, \bfnmB. C.\binitsB. C., \bauthor\bsnmRussellEggitt, \bfnmI.\binitsI., \bauthor\bsnmCortinaBorja, \bfnmM.\binitsM. \AND\bauthor\bsnmRahi, \bfnmJ. S.\binitsJ. S. (\byear2015). \btitleStudy of optimal perimetric testing in children (OPTIC): Normative visual field values in children. \bjournalOphthalmology \bvolume122 \bpages17111717. \endbibitem
 Pinheiro and Bates (1995) {barticle}[author] \bauthor\bsnmPinheiro, \bfnmJC\binitsJ. \AND\bauthor\bsnmBates, \bfnmDM\binitsD. (\byear1995). \btitleApproximations to the loglikelihood function in the nonlinear mixedeffects model. \bjournalJournal of Computational and Graphical Statistics \bvolume4 \bpages1235. \endbibitem
 Pinheiro and Bates (2000) {bbook}[author] \bauthor\bsnmPinheiro, \bfnmJ. C.\binitsJ. C. \AND\bauthor\bsnmBates, \bfnmD. M.\binitsD. M. (\byear2000). \btitleMixedeffects models in S and SPLUS. \bpublisherSpringer Verlag, \baddressNew York. \endbibitem
 Pinheiro et al. (2014) {bmanual}[author] \bauthor\bsnmPinheiro, \bfnmJ.\binitsJ., \bauthor\bsnmBates, \bfnmD.\binitsD., \bauthor\bsnmDebRoy, \bfnmS.\binitsS., \bauthor\bsnmSarkar, \bfnmD.\binitsD. \AND\bauthor\bsnmR Core Team (\byear2014). \btitlenlme: Linear and nonlinear mixed effects models \bnoteR package version 3.1117. http://CRAN.Rproject.org/package=nlme. \endbibitem
 Powell (1994) {binbook}[author] \bauthor\bsnmPowell, \bfnmJames L.\binitsJ. L. (\byear1994). \btitleEstimation of semiparametric models In \bbooktitleHandbook of Econometrics \bvolumeVolume 4 \bchapter41, \bpages24432521. \bpublisherElsevier. \endbibitem
 Reich, Bondell and Wang (2010) {barticle}[author] \bauthor\bsnmReich, \bfnmBrian J.\binitsB. J., \bauthor\bsnmBondell, \bfnmHoward D.\binitsH. D. \AND\bauthor\bsnmWang, \bfnmHuixia J.\binitsH. J. (\byear2010). \btitleFlexible Bayesian quantile regression for independent and clustered data. \bjournalBiostatistics \bvolume11 \bpages337352. \endbibitem
 Rigby and Stasinopoulos (2005) {barticle}[author] \bauthor\bsnmRigby, \bfnmRA\binitsR. \AND\bauthor\bsnmStasinopoulos, \bfnmDM\binitsD. (\byear2005). \btitleGeneralized additive models for location, scale and shape. \bjournalJournal of the Royal Statistical Society C \bvolume54 \bpages507554. \endbibitem
 Ritz and Spiegelman (2004) {barticle}[author] \bauthor\bsnmRitz, \bfnmJohn\binitsJ. \AND\bauthor\bsnmSpiegelman, \bfnmDonna\binitsD. (\byear2004). \btitleEquivalence of conditional and marginal regression models for clustered and longitudinal data. \bjournalStatistical Methods in Medical Research \bvolume13 \bpages309323. \endbibitem
 R Core Team (2016) {bmanual}[author] \bauthor\bsnmR Core Team (\byear2016). \btitleR: A Language and Environment for Statistical Computing \bpublisherR Foundation for Statistical Computing, \baddressVienna, Austria. \endbibitem
 Vonesh et al. (2002) {barticle}[author] \bauthor\bsnmVonesh, \bfnmE. F.\binitsE. F., \bauthor\bsnmWang, \bfnmH.\binitsH., \bauthor\bsnmNie, \bfnmL.\binitsL. \AND\bauthor\bsnmMajumdar, \bfnmD.\binitsD. (\byear2002). \btitleConditional secondorder generalized estimating equations for generalized linear and nonlinear mixedeffects models. \bjournalJournal of the American Statistical Association \bvolume97 \bpages271283. \endbibitem
 Wang (2012) {barticle}[author] \bauthor\bsnmWang, \bfnmJing\binitsJ. (\byear2012). \btitleBayesian quantile regression for parametric nonlinear mixed effects models. \bjournalStatistical Methods & Applications \bvolume21 \bpages279295. \endbibitem
 Yang, Wang and He (2016) {barticle}[author] \bauthor\bsnmYang, \bfnmYunwen\binitsY., \bauthor\bsnmWang, \bfnmHuixia Judy\binitsH. J. \AND\bauthor\bsnmHe, \bfnmXuming\binitsX. (\byear2016). \btitlePosterior inference in Bayesian quantile regression with asymmetric Laplace likelihood. \bjournalInternational Statistical Review \bvolume84 \bpages327â344. \bnotedoi:10.1111/insr.12114. \endbibitem
 Yuan and Yin (2010) {barticle}[author] \bauthor\bsnmYuan, \bfnmY.\binitsY. \AND\bauthor\bsnmYin, \bfnmG.\binitsG. (\byear2010). \btitleBayesian quantile regression for longitudinal studies with nonignorable missing data. \bjournalBiometrics \bvolume66 \bpages105114. \endbibitem