Dependent data arise in many studies. For example, children with the same parents or living in neighboring geographic areas tend to be more alike in many characteristics than individuals chosen at random from the population at large; observations taken repeatedly on the same individual are likely to be more similar than observations from different individuals. Frequently adopted sampling designs, such as cluster, multilevel, spatial, and repeated measures (or longitudinal or panel), may induce this dependence, which , the analysis of the data needs to take into due account. In a previous publication (Geraci and Bottai, Biostatistics 2007), we proposed a conditional quantile regression model for continuous responses where a random intercept was included along with fixed-coefficient predictors to account for betweensubjects dependence in the context of longitudinal data analysis. Conditional on the random intercept, the response was assumed to follow an asymmetric Laplace distribution. The approach hinged upon the link existing between the minimization of weighted least absolute deviations, typically used in quantile regression, and the maximization of a Laplace likelihood. As a follow up to that study, here we consider an extension of those models to more complex dependence structures in the data, which are modeled by including multiple random effects in the linear conditional quantile functions. Differently from the Gibbs * Draft version: June 1, 2011. Copyright to this paper remains with the authors or their assignees. Users may produce this paper for their own personal use, but distributing or reposting to other electronic bulletin boards or archives, may not be done without the written consent of the authors. To cite this paper: Geraci, M. and Bottai, M. (1 June 2011). Linear Quantile Mixed Models. Unpublished manuscript. 1 sampling expectation-maximization approach proposed previously, the estimation of the fixed regression coefficients and of the random effects covariance matrix is based on a combination of Gaussian quadrature approximations and optimization algorithms. The former include Gauss-Hermite and Gauss-Laguerre quadratures for, respectively, normal and doubleexponential (i.e., symmetric Laplace) random effects; the latter include a modified compass search algorithm and general purpose optimizers. As a result, some of the computational burden associated with large Gibbs sample sizes is avoided. We also discuss briefly an estimation approach based on generalized Clarkes derivatives. Finally, a simulation study is presented and some preliminary results are shown.