Is it possible to use location-scale family distributions for mixed effects modeling? Is it possible to use location-scale family distributions for mixed modeling or generalized estimation equations? The only location-scale family package I know of is GAMLSS in R which is for additive models using splines. 
Is there any practical application of ls family distributons in mixed modeling? 
 A: A couple of comments:


*

*Generalized Linear Mixed Models (GLMMs) have the following general representation:
$$\left\{
\begin{array}{l}
Y_i \mid b_i \sim \mathcal F_\psi,\\\\
b_i \sim \mathcal N(0, D),
\end{array}
\right.$$
where $Y_i$ is the response for the $i$-th sample unit and $b_i$ is the vector of random effects for this unit. The response $Y_i$ conditional on the random effects has a distribution $\mathcal F$ parameterized by the vector $\psi$, and the random effects are typically assumed to follow a multivariate normal distribution with mean 0 and variance-covariance matrix $D$. Some standard GLMMs assume that the distribution $\mathcal F_\psi$ is the binomial, Poisson, negative binomial, Beta, or Gamma distribution. However, you can, in principle, define $\mathcal F_\psi$ as you like, and hence, you could use a distribution from the location-scale family. The practical question is whether your specific choice is currently available in software. In that regard, my GLMMadaptive package allows you (the user) to define your own distribution. For examples, check here.

*Generalized Estimating Equations (GEEs) are a semi-parametric approach. You need to define a model for the mean and variance of your outcome and a working correlation structure. But you do not need to define the distribution of your outcome. For example, when in geepack package in R in the geeglm() function, you specify family = poisson(), you do not actually assume a Poisson distribution for your outcome. What you do is that you specify the link function (e.g., the default is the log link), and you also specify the variance function, i.e., that the variance of your outcome is equal to the mean times an over-dispersion parameter (GEEs automatically allow for over-dispersion).

A: I don't know about the GEE but for random-effects: 

They are two functions for fitting random effects wthin a GAMLSS
  model, random() and re(). The function random() is based on the
  original random() function of Trevor Hastie in the package gam. TIn
  our version the function has been modified to allow a "local" maximum
  likelihood estima- tion of the smoothing parameter lambda. This method
  is equivalent to the PQL method of Breslow and Clayton (1993) applied
  at the local iterations of the algorithm. In fact for a GLM model and
  a simple random effect it is equivalent to glmmPQL() function in the
  package MASS see Venables and Ripley (2002). Venables and Ripley
  (2002) claimed that this iterative method was first introduced by
  Schall (1991). Note that in order for the "local" maximum likelhood
  estimation procedure to operate both argument df and lambda has to be
  NULL. The function re() is an interface for calling the lme() function
  of the package nlme. This gives the user the abilty to fit comlpicated
  random effect models while the assumtion of the normal dis- tribution
  for the response variable is relaxed. The theoretical justification
  cames again from the fact that this is a PQL method, Breslow and
  Clayton (1993).

https://cran.r-project.org/web/packages/gamlss/gamlss.pdf
For more information:
http://www.gamlss.com/wp-content/uploads/2013/01/Lancaster-booklet.pdf
