Biostatistics Advance Access originally published online on April 24, 2006
Biostatistics 2007 8(1):140-154; doi:10.1093/biostatistics/kxj039
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Quantile regression for longitudinal data using the asymmetric Laplace distribution
Department of Epidemiology and Biostatistics, University of South Carolina, 800 Sumter Street, Columbia, SC 29208, USA geraci{at}gwm.sc.edu
* To whom correspondence should be addressed.
| SUMMARY |
|---|
|
|
|---|
In longitudinal studies, measurements of the same individuals are taken repeatedly through time. Often, the primary goal is to characterize the change in response over time and the factors that influence change. Factors can affect not only the location but also more generally the shape of the distribution of the response over time. To make inference about the shape of a population distribution, the widely popular mixed-effects regression, for example, would be inadequate, if the distribution is not approximately Gaussian. We propose a novel linear model for quantile regression (QR) that includes random effects in order to account for the dependence between serial observations on the same subject. The notion of QR is synonymous with robust analysis of the conditional distribution of the response variable. We present a likelihood-based approach to the estimation of the regression quantiles that uses the asymmetric Laplace density. In a simulation study, the proposed method had an advantage in terms of mean squared error of the QR estimator, when compared with the approach that considers penalized fixed effects. Following our strategy, a nearly optimal degree of shrinkage of the individual effects is automatically selected by the data and their likelihood. Also, our model appears to be a robust alternative to the mean regression with random effects when the location parameter of the conditional distribution of the response is of interest. We apply our model to a real data set which consists of self-reported amount of labor pain measurements taken on women repeatedly over time, whose distribution is characterized by skewness, and the significance of the parameters is evaluated by the likelihood ratio statistic.
Keywords: Asymmetric Laplace distribution; Clinical trials; Markov Chain Monte Carlo; Quantile regression; Random effects
| 1. INTRODUCTION |
|---|
|
|
|---|
In the last few years, quantile regression (QR) has become a more widely used technique to describe the distribution of a response variable given a set of explanatory variables. Median regression is a special case. In their seminal work, Koenker and Bassett (1978)
The foundations of the methods for independent data are now consolidated, and some computational commands for the analysis of such data are provided by most of the available statistical software (e.g. R/S-PLUS, SAS, and Stata).
We explore the use of QR for the analysis of longitudinal data. These kinds of data occur in a wide variety of studies, such as clinical trials, panel studies, and epidemiological surveys. The within-subject variability due to the measurements on the same subject should be accounted for to avoid bias in the parameter's estimate. As motivating example, in Section 4 we consider a clinical trial study on the effectiveness of a medication for relieving labor pain. The marginal distribution of the response (Figure 2) is characterized by a marked skewness. In addition, the skewness of the conditional distribution changes sign over time. In such a case, the appropriateness of mean regression should be seriously questioned: it would be highly inefficient to estimate the location, and, more important, entirely inadequate to estimate the quantiles of the conditional distribution of the response.
|
We propose a method to estimate the parameter of conditional quantile functions with subject-specific, location-shift random effects. This approach is analogous to estimating mean regressions with random intercepts. We introduce a likelihood based on the asymmetric Laplace density (ALD), which provides an automatic choice of the optimal level of penalization. In a simulation study, our strategy outperformed some of the best competing models.
Several works have appeared in the literature that propose methods for repeated measurements data. Koenker (2004)
considered the penalized interpretation of the classical random-effects estimator in order to estimate quantile functions with subject-specific fixed effects. The variability induced in the estimation process by these effects was controlled by some forms of shrinkage, or regularization, whose degree and type required a suitable choice of a penalization term. The application's areas include reference charts in medicine (Heagerty and Pepe, 1999
), (Cole and Green, 1992
) and clinical trials (Robins and others, 1995), (Lipsitz and others, 1997
). Jung (1996)
proposed a quasi-likelihood approach for median regression estimation, where the dependency structure is modeled by the covariance matrix of the indicator functions. Lipsitz and others (1997) made use of Jung's work and described the use of QR methods to analyze longitudinal data on CD4 cell counts, and following Robins and others (1995) accounted for missing at random dropouts by weighting the generalized estimating equations. They estimated the regression parameters' standard errors by using the bootstrap.
The remainder of this paper proceeds as follows: in Section 2, we briefly present the conditional quantile functions for independent data and the ALD. We then propose a model for repeated measurements along with related inferential procedures. In Section 3, we present the results of a simulation. In Section 4, we apply the methodology to clinical trial data, and, in Section 5, we conclude the paper with some final remarks.
| 2. MODEL AND METHOD |
|---|
|
|
|---|
Consider data in the form
, for
, where
are independent scalar observations of a continuous random variable with common cumulative distribution function
, whose shape is not exactly known, and
are row p-vectors of a known design matrix X.
The linear conditional quantile functions are defined as
|
| (2.1) |
where
,
, and
is a column vector of length p with unknown fixed parameters.
We refer to the
th regression quantile as any solution
,
, to the minimization problem
![]() | (2.2) |
In order to highlight the
-distributional dependency, the parameter
and the solution
should be indexed by
, i.e.
and
. For sake of simplicity, however, we will omit this notation in the remainder of the paper.
A possible parametric link between the minimization of the sum of the absolute deviates in (2.2) and the maximum likelihood theory is given by the ALD. This skewed distribution appeared in the papers by Koenker and Machado (1999)
and Yu and Moyeed (2001)
, among others. We say that a random variable Y is distributed as an ALD with parameters
,
, and
, and we write it as
, if the corresponding probability density is given by
|
| (2.3) |
where
is the loss function,
is the indicator function,
is the skewness parameter,
is the scale parameter, and
is the location parameter. The support of the random variable Y is the real line. It should be noted that the loss function
assigns weight
or
to the observations greater or, respectively, less than
and that
. Therefore, the distribution splits along the scale parameter into two parts, one with probability
to the left, and one with probability
to the right. See Yu and Zhang (2005)
for further properties and generalizations of this distribution.
Set
and
. Assuming that
, then the likelihood for N independent observations is, bar a proportionality constant,
![]() | (2.4) |
If we consider
a nuisance parameter, then the maximization of the likelihood in (2.4) with respect to the parameter
is equivalent to the minimization of the objective function
in (2.2). Thus, the ALD proves useful as a unifying bridge between the likelihood and the inference for QR estimation. Koenker and Machado (1999)
introduced a goodness-of-fit process for QR and related inference processes, and they also considered the likelihood ratio statistic under the parametric assumption of a Laplacean distribution for the error term. For a full Bayesian approach, see Yu and Moyeed (2001)
, among others.
We now consider repeated measurements data in the form
, for
and
, where
are row p-vectors of a known design matrix
and
is the jth scalar measurement of a continuous random variable on the ith subject.
A possible penalized approach is offered by Koenker (2004)
, who proposed an
penalty term and considered
![]() | (2.5) |
where the weight w controls the relative influence of the
th quantile on the estimation of the fixed subject-specific effects
, and
is the penalization parameter.
The parameter
represents a fundamental issue. Its value needs to be arbitrarily set, yet its choice heavily influences the inference for the parameter
, as confirmed by the simulation results in Section 3. Besides, following a penalized approach, estimating a
-distributional individual effect would be impracticable. Also, as the sample size increases, so does the dimension of the parameter. We circumvented these drawbacks by letting the location-shift effects be random and the parameter of their distribution be dependent on
.
We define the linear mixed quantile functions of the response
|
| (2.6) |
where
is the inverse of the cumulative distribution function of the response conditional on a location-shift random effect
.
We assume that
, conditionally on
, for
and
, are independently distributed according to the ALD
|
| (2.7) |
where
is the linear predictor of the
th quantile. In our model, the conditional quantile to be estimated,
, is fixed and known. The random effects induce a correlation structure among observations on the same subject. We assume that the
are identically distributed according to some density
characterized by a
-dependent parameter
(i.e.
), and they are mutually independent. Also, we assume that
are independent, and
and
are independent of one another. In Section 2.4, we elaborate on the distribution of the random effects.
Let
and
be the density for the ith subject conditional on the random intercept
. The complete-data density of
, for
, is then given by
|
| (2.8) |
where
is the density of
and
is the parameter of interest.
If we let
and
, the joint density of
based on N subjects is given by
![]() |
We obtain the marginal density of
by integrating out the random effects, leading to
. Thus, the inference about the parameter
should be based on the marginal likelihood,
. The latter involves an integral that, in general, does not have a closed form solution.
To estimate
, we propose a Monte Carlo EM algorithm, which has been widely applied to generalized linear mixed models (LMMs) (see for example, Meng and Van Dyk, 1998
; Booth and Hobert, 1999). The presence of intractable integrals in the E-step can be overcome, at some computational expenses, by using a sampler. Handling missingness due to nonresponse also represents an attractive feature of such an algorithm (see for example, Ibrahim and others, 1999, 2001).
We write the E-step for the ith subject at the
th EM iteration as
![]() | (2.9) |
where
,
is the log likelihood for the ith observation, and the expectation is taken with respect to the "missing" data
, conditioned on the observed data
.
The approximation of the expected value of the log likelihood can be carried out with a suitable Monte Carlo technique, which depends on the target density to be sampled.
If the random effects were known and the model in (2.6) is correctly specified, we could define
![]() | (2.10) |
where
, and apply a linear programming algorithm to solve
![]() | (2.11) |
since the
values would be no longer dependent.
Since the random effects are unobserved, we propose to take a sample
of size
from the density of the random effects conditional to the response
|
| (2.12) |
and to approximate the expected value in (2.9) with
|
|
To obtain the maximum likelihood estimate of the parameter
for the
th quantile model, we applied the following iterative procedure:
Step 1: Initialize the parameter
, set
, and plug
in (2.12);
Step 2: Draw a sample
from
, for
, independently, by using a sampling technique, e.g. via the Gibbs sampler with the adaptive rejection sampling algorithm (Gilks and others, 1995);
Step 3: Solve
![]() |
where
, set the solution of the minimization problem equal to
, and calculate
![]() |
and
|
|
where
is the maximum likelihood estimate (estimator) of
, conditional on
;
Step 4: Set
and repeat Steps 13 until the parameter
reaches the convergence.
The minimization problem in Step 3 is equivalent to a weighted QR with an offset and thus can be easily solved. Then, at each iteration,
is the maximum likelihood estimate of
, given
.
The sampling step, obviously, requires evaluating the Markov chain's convergence (see for example, Cowles and Carlin, 1996
) and the sensitivity of the parameter's estimates to different burn-in and sample sizes. Alternative strategies for choosing the size of the Gibbs samples are available. One can consider a constant size
, for all i and each EM iteration, or increase
at each EM iteration.
The availability of a parametric likelihood function allowed us to evaluate a likelihood ratio statistic to test the null hypotheses
,
. A bootstrap method can also be used to obtain an estimate of the standard error of the estimates of
and it can be implemented for longitudinal data by regarding the matrix
as the basic resampling unit (Lipsitz and others, 1997). Buchinsky (1995)
studied the performance of various estimators for the asymptotic variance matrix for QR models and found that the bootstrap procedure yielded the best results, performing well for relatively small samples and also when the data exhibit heteroscedasticity.
In Section 1, we argued that the usual choice of the normal distribution for the error model is sometimes questionable and, most important, leaves the estimation out to the influence of the outliers. The use of a joint normal distribution within the mixed models simplifies considerably the analytical form of the expected value of the log likelihood when using the EM algorithm and thus reduces the time to convergence. Robustness issues, however, apply to the random effects as well (see for example, Lange and others, 1989; Richardson and Welsh, 1995
; Richardson, 1997
). For our model, we considered the normal distribution, although it is possible to consider a vast range of symmetric as well as asymmetric distributions such as the asymmetric Laplace. Simulation results are reported in Section 3.
Assume that
, identically and independently for
. The joint density of
for the ith subject is given by
![]() | (2.13) |
Note that the conditional distribution of the random effect
is log-concave in
, where the definition of log-concavity here is derivative-free (Gilks, 1995
), i.e. a function g, which is a density with respect to Lebesgue measure and its domain D is an interval of the real line, is defined to be log-concave if
|
|
Thus, the Gibbs sampling in conjunction with the adaptive rejection sampling can be applied and the Monte Carlo E-step carried out by sampling from
, for each subject i,
.
The maximum likelihood estimator of
,
, at Step 3 of the iterative procedure presented in Section 2.3, is then given by
![]() |
Assuming that
, the joint density of
for the ith subject would be
![]() | (2.14) |
The Gibbs sampling can be applied, since the conditional distribution
derived from the joint density in (2.14) is still log-concave in
. Thus, the maximum likelihood estimator of
is given by
![]() |
Note also that it is possible to consider
and to estimate the degree of skewness
of the random effects in the case the assumption of symmetry of u is untenable.
Consider the median regression for example, i.e.
. If we set
, the density in (2.14) turns into
![]() |
where the argument of the exponential has a similar form as in the penalized quantile regression (PQR) in (2.5).
| 3. SIMULATION STUDY |
|---|
|
|
|---|
We conducted a simulation study to evaluate the performance of the proposed method for two different quantile functions,
and
. We used
subjects and
, for
, subject measurements. To generate the data, we used the location-shift model
|
| (3.1) |
The
values were set equal to the values
, constant throughout the simulation study.
Different probability distributions were considered for the error term: the standard normal,
, the Student's t distribution with 3 degrees of freedom,
, and the chi-square with 3 degrees of freedom,
. The random effects were generated by using a standard normal. For the data
,
, generated using the model in (3.1), we assumed the joint density in (2.13) and we estimated the parameter
using the method proposed in Section 2.3. We simulated 1000 replications from each distribution of the error term.
For each replication, the Monte Carlo E-step was carried out via the Gibbs sampler (see Section 2.4) with a constant sample size of
,
, and a burn-in size of 2500.
We compared the quantile regression with random effects (QRRE) with the PQR reported in (2.5), using
,
, and
. The latter is equivalent to the ordinary QR, in which the fixed effects are fully shrunk to zero. The choice of these values was made empirically, since there are no theoretically based or practically accepted guidelines as to how to choose the level of shrinkage. The ratio of the estimated standard deviations of
and
, when using the normal and the Student's t distributions, has been sometimes proposed as reference value for
.
We also compared the proposed model with the ordinary LMM when the error term distribution was symmetric (i.e. Normal and t-distribution). When the location model was of interest, comparing the performance of the least squares with the least absolute estimators proved to be useful to evaluate the robustness of our model.
In Table 1, for each estimator's component
,
, and for each model, we reported the estimated relative bias averaged over the simulations,
![]() |
where
is the estimated quantile for the rth replication and
is the corresponding true quantile of the conditional distribution
. We also reported the estimated relative efficiency, averaged over the simulations, using the Monte Carlo variance of the QRRE model in the denominator
![]() |
where
and
.
|
The results show clearly that our method generally performs better than the other across the simulated scenarios. The
penalization model performs as well as ours for some values of the penalization parameter, but poorly for others, which is an undesirable situation in practical applications.
The relevance of the bias in the mean squared error is particularly evident for the estimates of the first quartile. A poor choice for the penalization parameter may lead to large bias. In fact, in our simulations the relative bias ranged from
to
. It can be noted that, as
increases, the bias of the penalized regression decreases from large positive to large negative values. This suggests that there might be a reasonable value of the penalization parameter that ensures no bias. It can be argued that our model automatically provides a nearly optimal degree of shrinking. When estimating the first quartile, the loss of efficiency of the penalized regressions, with respect to our model, ranged from 12 to
and from 17 to
for the intercept and the slope parameters, respectively.
The gap in terms of mean squared error between our model and the penalized regression decreased when estimating the median. We observed a smaller loss of efficiency and a leveling off of the relative bias for the symmetric distributions. The bias associated with the random effects QR, except for the case in which the error was
distributed, turned out to be slightly higher than the others. The reason might be due to the Monte Carlo error that is less conspicuous when the estimated bias is high as in the case of the first quartile.
As expected, the least-squares estimation outperforms the median regression when the error term is normally distributed, but it yields highly inefficient estimates when the distribution has heavy tails.
A small simulation was also performed on the quantile
, and the results (not shown here) were equivalent to the ones obtained for the first quartile.
| 4. LABOR PAIN DATA: COMPARISONS |
|---|
|
|
|---|
The labor pain data were reported by Davis (1991)
women in labor, of which 43 were randomly assigned to a pain medication group and 40 to a placebo group. The response was measured every 30 min on a 100-mm line, where 0 means no pain and 100 means extreme pain. A nearly monotone pattern of missing data was found for the response variable and the maximum number of measurements for each woman was six. In Figure 1, a boxplot shows the median pain over time for the placebo and the medication groups. Some degree of distributional dependence on the temporal course of the quartiles of the response is evident, especially for the placebo group.
|
Jung (1996)
|
| (4.1) |
where
is the amount of pain for the ith patient at time j,
is the treatment indicator taking 1 for placebo and 0 for treatment, and
is the measurement time divided by 30 min, and
is a zero-median error term.
Let
![]() |
where
is the last measurement time for the ith woman.
Figure 2 reports the triangular kernel estimates of the marginal density of the response in the entire sample, in the pain medication, and in the placebo group. As optimal level of smoothing, the bandwidth of the kernel function was set equal to twice its standard deviation. Other kernel functions (e.g. Gaussian, rectangular, Epanechnikov, biweight, cosine) and different bandwidths were explored; they all showed similar pictures. The resulting shapes are clearly far from the symmetric, bell-shaped normal distribution. The estimated conditional densities for the pain medication and the placebo group appear to be very different from one another. The former presents a high density peak toward the minimum of the pain score scale, indicating that the treatment is effective. The estimated mean and median pain score are equal to 18.87 and 9.00, respectively. The latter is approximately uniform, with estimated mean and median pain score equal to 49.51 and 45.50, respectively.
Given the longitudinal nature of the observations, we should expect some degree of correlation between the measurements pertaining to the same subject. The application of random-effects regression models to such data has a very widespread success. For the case at hand, however, the estimation of a LMM would drive to misleading, inappropriate inferential results. In fact, the assumption of normality is clearly unreasonable. Further, and most importantly, the median of the pain score might be more informative than the mean for such a skewed distribution, for the mean is highly influenced by unusually large observations, which may occur with substantial probability in skewed populations.
We then estimated the median regression with random effects
|
| (4.2) |
where
is a random intercept and, for
, we assumed the density given in (2.13). According to the results of the simulation in Section 3, the model in (4.17) and the use of the normal Laplace joint distribution should be superior to the LMM. The advantage is two-fold:
- the quantile model provides an estimate for the median, which is a more informative measure of central tendency than the mean for such a skewed distribution;
- the quantile model makes no assumptions about the distribution of the residual error which allows correct inference about other quantiles.
The Monte Carlo E-step was carried out via the Gibbs sampler with a constant sample size of
, and a burn-in size of 5000.
We also estimated the corresponding penalized median regression with fixed
values. Lacking any practical guidelines as to how to choose the penalizing constant,
, we set it at two arbitrarily chosen values, 0.01 and 5.
The estimates of the parameter
are reported in Table 2, along with the results of the quasi-likelihood approach showed by Jung (1996)
and the least-squares estimates of the random intercept model. The significance of the parameter estimates was obtained by using a likelihood ratio statistic. According to our fitted model, the baseline pain does not differ significantly between the two groups, while time and its interaction with treatment affect the amount of pain. We also evaluated approximated
confidence intervals by using a bootstrap technique with a sample of size 499, and no differences were observed.
|
The median regression with random effects and the random-effects regression tally in terms of the magnitude and the significance of the estimate of
. To some extent, this result supports the use of random effects. Moreover, the rest of the QRs detected the significance of
only, which might be due to small power. By examining the plot in Figure 1, it is evident that a first-degree polynomial model does not adequately describe the median pain score over time. We added square and cubic terms for both the time, T, and the interaction between time and treatment, RT, to model (4.16). We also estimated the first and the third quartiles and we used the Akaike information criterion (AIC) to evaluate the fit of each model. The possibility to describe more generally the conditional distribution of the response through the estimation of its quantiles, while accounting for the dependence among the observations, represents a great advantage of our models with respect to normal mixed regressions. The latter, in fact, cannot realistically approximate the quantiles of a skewed or heavy-tailed distributions.
Figure 3 presents a summary of the QR results. The magnitude and the sign of the skewness of the pain score distribution for the placebo group change over time. In fact, the skewness remains positive in the interval between 30 and 90 min, circa, and it turns negative after 90 min. This result suggests that there is a modest fraction of women that is more sensitive to pain or that is not influenced by the placebo in the first halftime of the treatment, and a modest fraction of women that are less sensitive to pain or influenced by the placebo in the second half. The skewness for the pain medication group is always positive and its magnitude increases over time, suggesting that the medication is effective constantly over time for at least half of the group. A more detailed analysis of the distributional composition of the two groups could be extended to other quantiles. In general, however, the study of the tails of the distribution requires a suitable number of observations to ensure stability of the estimates.
|
For the labor pain data, the lowest AIC value was obtained when adding the terms
,
, and
for both the first quartile and the median. For the third quartile, the AIC was minimized when adding the terms
,
, and
. On one hand, the level of pain for the most sensitive women of the placebo group seems to increase rapidly over time and then to stabilize at the median value, suggesting that they did not experience a placebo effect or they did it to a lesser extent. On the other hand, the third quartile of the medication group seems still on a rise at the end of the observational time. The estimate and the significance of the coefficients (not shown here) of a polynomial of degree higher than one must be interpreted with caution in this case. The value of the self-measured pain is, in fact, constrained to be less than or equal to 100 and the curvature of the upper quartile for the placebo group might be unnatural. These results, therefore, are purely an indication of the inadequacy of the model in (4.1).
We monitored convergence of the Gibbs sampler through the diagnostic techniques suggested by Cowles and Carlin (1996)
and we did not observe any lack of convergence. Lag-1 autocorrelation was always within
for the first and the third quartiles and within
for the median model.
| 5. FINAL REMARKS |
|---|
|
|
|---|
Extending the QR methods to include random effects is currently an active and promising area of research. We proposed a novel linear QR model with a subject-specific random intercept that accounts for within-group correlation. In order to estimate the model parameter, we applied an iterative procedure. At each step, a sample from the conditional density of the random effects was drawn via a Gibbs sampler with an adaptive rejection metropolis sampling algorithm. Then, the predicted random effects were subtracted from the response variable and a simple linear programming algorithm was applied to estimate the regression quantiles. Finally, we used the likelihood ratio statistic to test the null hypotheses that the fixed effects are equal to zero. Alternatively, a bootstrap technique can be used to estimate the variance of the
's estimator. In our experience, the ALD seems to be very promising within the regression quantile models, since it clearly represents a suitable error law for the least absolutes estimator. Also, the regression median with random effects is seen to be more efficient than the least-squares estimator in the LMM for any distribution for which the median is more efficient than the mean in the location model.
Further studies should be aimed at alleviating the computational burden related to the presence of the multiple integral with respect to the random effects. The combined use of generalized directional derivatives and numerical integration techniques (e.g. the Clark derivative and the adaptive quadrature) represents a possible solution to the maximization of the marginal likelihood, which could be easily implemented with standard statistical languages such as R. The extensions to conditional quantile functions with more than a single random effect and to the use of different link functions also represent interesting areas of development.
| ACKNOWLEDGMENTS |
|---|
Conflict of Interest: None declared.
| REFERENCES |
|---|
|
|
|---|
-
Booth JG and Hobert JP. (1999) Maximizing generalized linear mixed model likelihoods with an automated Monte Carlo EM algorithm. Journal of the Royal Statistical Society, Series B 61:26585.[CrossRef]
Buchinsky M. (1995) Estimating the asymptotic covariance matrix for quantile regression models: a Monte Carlo study. Journal of Econometrics 68:30338.[CrossRef]
Cole TJ and Green PJ. (1992) Smoothing reference centile curves: the LMS method and penalized likelihood. Statistics in Medicine 11:130519.[Web of Science][Medline]
Cowles MK and Carlin BP. (1996) Markov chain Monte Carlo convergence diagnostics: a comparative review. Journal of the American Statistical Association 91:883904.[CrossRef]
Davis CS. (1991) Semi-parametric and non-parametric methods for the analysis of repeated measurements with applications to clinical trials. Statistics in medicine 10:195980.[Web of Science][Medline]
Gilks WR. (1995) Derivative-free adaptive rejection sampling for Gibbs sampling. In Bernardo JM, Berger JO, David AP, Smith AFM (Eds.). Bayesian Statistics 4(Clarendon, Oxford) pp. 6419.
Gilks WR, Best NG, Tan KKC. (1995) Adaptive rejection metropolis sampling within Gibbs sampling. Applied Statistics 44:45572.[CrossRef]
Heagerty PJ and Pepe MS. (1999) Semiparametric estimation of regression quantiles with application to standardizing weight for height and age in US children. Journal of the Royal Statistical Society, Series C 48:53351.[CrossRef]
Ibrahim JG, Chen MH, Lipsitz SR. (1999) Monte Carlo EM for missing covariates in parametric regression models. Biometrics 99:10411.
Ibrahim JG, Chen MH, Lipsitz SR. (2001) Missing responses in generalized linear mixed models when the missing data mechanism is nonignorable. Biometrika 88:55164.
Jung S. (1996) Quasi-likelihood for median regression models. Journal of the American Statistical Association 91:2517.[CrossRef]
Koenker R. (2004) Quantile regression for longitudinal data. Journal of Multivariate Analysis 91:7489.[CrossRef]
Koenker R. (2005) Quantile Regression(Cambridge University Press, New York).
Koenker R and Bassett G. (1978) Regression quantiles. Econometrica 46:3350.[CrossRef]
Koenker R and Machado J. (1999) Goodness of fit and related inference processes for quantile regression. Journal of the American Statistical Association 94:1296309.[CrossRef]
Lange KL, Little RJA, Taylor JMG. (1989) Robust statistical modeling using the t-distribution. Journal of the American Statistical Association 84:88196.[CrossRef]
Lipsitz SR, Fitzmaurice GM, Molenberghs G, Zhao LP. (1997) Quantile regression methods for longitudinal data with drop-outs: application to CD4 cell counts of patients infected with the human immunodeficiency virus. Journal of the Royal Statistical Society, Series C 46:46376.[CrossRef]
Meng XL and Van Dyk D. (1998) Fast EM-type implementations for mixed effects models. Journal of the Royal Statistical Society, Series B 60:55968.[CrossRef]
Richardson AM. (1997) Bounded influence estimation in the mixed linear model. Journal of the American Statistical Association 92:15461.[CrossRef]
Richardson AM and Welsh AH. (1995) Robust restricted maximum likelihood in mixed linear models. Biometrics 51:142939.[CrossRef]
Robins JM, Rotnitzky A, Zhao LP. (1995) Analysis of semiparametric regression models for repeated outcomes in the presence of missing data. Journal of the American Statistical Association 90:10621.[CrossRef]
Yu K and Moyeed RA. (2001) Bayesian quantile regression. Statistics and Probability Letters 54:43747.[CrossRef]
Yu K and Zhang J. (2005) A three-parameter asymmetric Laplace distribution and its extension. Communications in StatisticsTheory and Methods 34:186779.[CrossRef]
Received November 11, 2005; revised April 18, 2006; accepted for publication April 21, 2006.
![]()
CiteULike
Connotea
Del.icio.us What's this?
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||

















) and efficiency (
) for different error term distributions and quantile models: PQR with penalization parameter taking values
= 0.01, 5, and 

