#  lqmix: an R package for longitudinal data analysis via linear quantile
mixtures

Marco Alfò Dipartimento di Scienze Statistiche, Sapienza Università di Roma
Maria Francesca Marino Dipartimento di Statistica, Informatica, Applicazioni,
Università degli Studi di Firenze<EMAIL_ADDRESS>Maria
Giovanna Ranalli Dipartimento di Scienze Politiche, Università degli Studi di
Perugia Nicola Salvati Dipartimento di Economia e Management, Università di
Pisa

###### Abstract

The analysis of longitudinal data gives the chance to observe how unit
behaviors change over time, but it also poses series of issues. These have
been the focus of a huge literature in the context of linear and generalized
linear regression moving also, in the last ten years or so, to the context of
linear quantile regression for continuous responses. In this paper, we present
lqmix, a novel R package that helps estimate a class of linear quantile
regression models for longitudinal data, in the presence of time-constant
and/or time-varying, unit-specific, random coefficients, with unspecified
distribution. Model parameters are estimated in a maximum likelihood
framework, via an extended EM algorithm, and parameters’ standard errors are
estimated via a block-bootstrap procedure. The analysis of a benchmark dataset
is used to give details on the package functions.

Keywords: Quantile regression, random effects, EM algorithm, unobserved
heterogeneity, R, panel data.

##  1 Introduction

Quantile regression (Koenker and Bassett,, 1978) has become quite a popular
technique to model the effect of observed covariates on the conditional
quantiles of a continuous response of interest. This represents a well
established practice for the analysis of data when the focus goes beyond the
(conditional) mean. Comprehensive reviews on this topic can be found, among
others, in Koenker and Hallock, (2001), Yu et al., (2003), Koenker, (2005),
and Hao et al., (2007).

With longitudinal data, dependence between observations recorded on the same
statistical unit needs to be taken into account to avoid bias and inefficiency
in parameter estimates. Different modeling alternatives are available in the
literature for handling such a dependence (for a general presentation, see
e.g., Diggle et al.,, 2002; Fitzmaurice et al.,, 2004). A possibility consists
in including in the quantile regression model unit-specific random
coefficients that are meant to describe dependence between repeated
observations coming from the same unit, giving rise to a conditional model
specification (see Marino and Farcomeni,, 2015, for a review of available
alternatives). In this framework, a standard way of proceeding is based on
specifying a parametric distribution for the random coefficients, as
suggested, e.g., by Geraci and Bottai, (2007, 2014). An alternative consists
in leaving the distribution unspecified and estimating it from the observed
data, by using a finite mixture specification. This can result from
considering discrete, unit-specific, random coefficients having unspecific
distribution that remains constant or evolves over time according to a hidden
Markov chain, as proposed by Alfó et al., (2017) and Farcomeni, (2012),
respectively. A more flexible specification based on combining both time-
constant (TC) and time-varying (TV), unit-specific, discrete, random
coefficients may also be adopted, as proposed by Marino et al., (2018). When
compared to fully parametric alternatives, this semi-parametric approach
offers a number of specific advantages as it helps (i) avoid unverifiable
assumptions on the random coefficient distribution; (ii) account for extreme
and/or asymmetric departures from the homogeneous model; (iii) avoid integral
approximations, considerably reducing the computational effort in parameter
estimation.

In this paper, we describe the R package (R Core Team,, 2019) lqmix (Marino et
al.,, 2023), available from the Comprehensive R Archive Network (CRAN) at
https://CRAN.R-project.org/package=lqmix, which is intended to provide maximum
likelihood (ML) estimates for TC and/or TV mixtures of linear quantile
regression models for longitudinal data. An indirect estimation approach,
based on an extended Expectation-Maximization algorithm (EM - Dempster et
al.,, 1977) is employed.

The package at hand shares features with some relevant R packages available
for the analysis of longitudinal data on the CRAN repository. It is related to
packages lme4 (Bates et al.,, 2015) and lqmm (Geraci,, 2014), tailored,
respectively, to the analysis of general clustered observations via mixed
models for the (conditional) mean and the (conditional) quantiles of a
response. Here, TC random coefficients following a parametric distribution are
considered to model dependence between observations and a ML approach is
employed to derive parameter estimates. Further, it is strongly related to the
LMest R package (Bartolucci et al.,, 2017) developed to model the mean of
longitudinal observations via hidden Markov models (Zucchini and MacDonald,,
2009; Bartolucci et al.,, 2013). In this case, TV random coefficients with
unspecific distribution that evolves over time according a Markov chain are
considered to account for dependence between measures from the same unit, and
a ML approach is adopted to derive parameter estimates.

Other R packages that can be fruitfully related to lqmix are rqpd (Koenker and
Bache,, 2014) and npmlreg (Einbeck et al.,, 2018). The former allows for the
estimation of linear quantile regression models for panel data by considering
either a penalized fixed effect estimation method (Koenker,, 2004) or a
correlated random-effect method (Abrevaya and Dahl,, 2008). In both cases,
parameters are estimated by minimizing an extended quantile loss function. On
the other side, the npmlreg R package entails mixtures of (generalized) linear
models for clustered observations by employing a ML approach for parameter
estimation. Last, it is worth mentioning the quantreg R package (Koenker,,
2022) for estimation and inference in linear models for conditional quantiles.

The modeling of longitudinal data via quantile regression is also possible by
considering packages and environments out of the R world, although none allows
to consider random coefficients and the possibility of modeling them
dynamically. The Stata modules xtqreg (Machado and Santos Silva,, 2021) and
xtmdqr (Pons and Melly,, 2022) allow for the estimation of linear quantile
regression models for longitudinal data based on the use of fixed effects.
These modules differ in the way model parameters are estimated: the former
considers the method of moments introduced by Machado and Santos Silva,
(2019), while the latter builds up on the minimum distance approach described
by Galvao and Wang, (2015). The qregpd module (Baker,, 2016) is also worth to
be mentioned; it implements the quantile regression estimator developed in
Graham et al., (2015) for panel data. Last, the qreg Stata module (Machado et
al.,, 2021) allows for the estimation of linear quantile regressions under the
assumption of independent observations. Similarly does also the quantreg SAS
procedure.

The paper is structured as follows. Different proposals entailing finite
mixtures of linear quantile regression models are provided in Section 2, where
the case of TC, TV, and a combination of TC and TV random coefficients are
discussed. ML estimation is introduced in Section 3, where details of the EM
algorithm and the bootstrap procedure for standard error estimation are
described. The proposed R package is introduced in Section 4, where the
analysis of a benchmark dataset is discussed. The last section contains some
concluding remarks.

##  2 Modeling alternatives

A longitudinal study is planned to record the value of a continuous response
variable $Y$ and a set of $p$ explanatory variables $x_{1},\dots,x_{p}$ on a
sample of $n$ statistical units, at $T$ different measurement occasions. This
is a balanced study, with common (at least in unit-time) and equally spaced
occasions, in discrete time. Let $Y_{it}$ denote the value of the response for
unit $i=1,\dots,n,$ at occasion $t=1,\dots,T,$ $y_{it}$ the corresponding
realization, and $\mbox{\boldmath$x$}_{it}=(x_{it1},\dots,x_{itp})^{\prime}$ a
$p$-dimensional vector of covariates, for $i=1,\dots,n,$ and $t=1,\dots,T$ .
Further, let $\mbox{\boldmath$Y$}_{i}=(Y_{i1},\dots,Y_{iT})^{\prime}$,
$\mbox{\boldmath$y$}_{i}=(y_{i1},\dots,y_{iT})^{\prime}$, and
$\mbox{\boldmath$X$}_{i}=(\mbox{\boldmath$x$}_{i1},\dots,\mbox{\boldmath$x$}_{iT})^{\prime}$
be the $T$-dimensional vector of responses, its realization, and the $(T\times
p)$-dimensional design matrix for unit $i$. A frequent issue to address when
dealing with longitudinal studies is that of missingness. This may entail both
the response and the explanatory variables. Here, we assume that the latter
are fully observed, while an ignorable missingness (Rubin,, 1976) may affect
the outcome. That is, some units in the sample may present incomplete response
sequences due to either monotone or non-monotone missing data patterns (Little
and Rubin,, 2002). In this sense, a varying number of measures $T_{i}$ may be
available for each unit, even though missingness is assumed to be independent
from unobserved responses – M(C)AR assumption.

Random coefficient models represent a standard approach to analyze the effect
of observed explanatory variables on a response $Y$ that is repeatedly
observed over time. This also holds in the quantile regression framework,
where the interest is on modeling the (conditional) quantiles of the outcome
distribution as a function of fixed and random coefficients. To ensure
flexibility and avoid unverifiable parametric assumptions on the random
coefficient distribution, a specification based on finite mixtures represent a
viable strategy to adopt. For this purpose, we developed the lqmix R package.

In this section, we introduce the methodology underlying the proposed package,
and present some alternative formulations of quantile regression models to
deal with longitudinal data. Discrete random coefficients with unspecific
distribution are added to the model with the aim of capturing unit-specific
sources of unobserved heterogeneity due to omitted covariates. According to
the chosen specification, these may remain constant and/or evolve over time,
leading to a model based on TC and/or TV random coefficients, respectively.

In Section 2.1, we discuss the linear quantile regression model with TC random
coefficients; in Section 2.2, such a model is extended by considering TV
random coefficients while, in Section 2.3, the linear quantile regression
model with both TC and TV random coefficients is introduced.

###  2.1 Linear quantile mixtures with TC random coefficients

For a given quantile level $q\in(0,1)$, let $\mbox{\boldmath$\beta$}_{q}$
denote a quantile-dependent, $p$-dimensional, vector of parameters associated
with the (design) vector
$\mbox{\boldmath$x$}_{it}=(x_{it1},\dots,x_{itp})^{\prime}$. Also, let
$\mbox{\boldmath$z$}_{it}=(z_{it1},\dots,z_{itr})^{\prime}$ denote a further
set of $r\geq 1$ explanatory variables (not included in
$\mbox{\boldmath$x$}_{it}$) associated with a vector of unit- and quantile-
specific random coefficients
$\mbox{\boldmath$b$}_{i,q}=(b_{i1,q},\dots,b_{ir,q})^{\prime}$. These latter
account for unobserved heterogeneity that is not captured by the elements in
$\mbox{\boldmath$x$}_{it}$ and may be used to describe dependence between
repeated measurements from the same unit. In this sense, conditional on
$\mbox{\boldmath$b$}_{i,q}$, the longitudinal responses coming from the same
unit, $Y_{i1},\dots,Y_{it}$, are assumed to be independent (local independence
assumption) for a given $q$. For a given quantile level $q\in(0,1)$ and
conditional on the vector $\mbox{\boldmath$b$}_{i,q}$, the response $Y_{it}$
is assumed to follow an Asymmetric Laplace density (ALD - e.g., Yu and
Moyeed,, 2001)

$f_{y\mid
b}(y_{it}\mid\mbox{\boldmath$b$}_{i,q};q)=\left[\frac{q(1-q)}{\sigma_{q}}\right]\exp\left\\{-\rho_{q}\left[\frac{y_{it}-\mu_{it,q}}{\sigma_{q}}\right]\right\\}.$

Here, $\rho_{q}(\cdot)$ denotes the quantile asymmetric loss function (Koenker
and Bassett,, 1978), while $q,\sigma_{q}$, and $\mu_{it,q}$ denote the
skewness, the scale, and the location parameter of the distribution. The ALD
is a working model that is used to recast estimation of parameters for the
linear quantile regression model in a maximum likelihood framework. In detail,
the location parameter $\mu_{it,q}$ is modeled as

$\mu_{it,q}=\mbox{\boldmath$x$}_{it}^{\prime}\mbox{\boldmath$\beta$}_{q}+\mbox{\boldmath$z$}_{it}^{\prime}\mbox{\boldmath$b$}_{i,q}.$
(1)

The modeling structure is completed by the mixing distribution
$f_{b,q}({\mbox{\boldmath$b$}}_{i,q};{\mbox{\boldmath$\Sigma$}}_{q})$, i.e.,
the distribution of the random coefficients ${\mbox{\boldmath$b$}}_{i,q}$,
where $\mbox{\boldmath$\Sigma$}_{q}$ identifies a (possibly quantile-
dependent) covariance matrix.

Rather than specifying such a distribution parametrically as in, e.g., Geraci
and Bottai, (2014), Alfó et al., (2017) proposed to leave it unspecified and
estimate it directly from the observed data via a NonParametric Maximum
Likelihood approach (NPML – Laird,, 1978; Lindsay, 1983a, ; Lindsay, 1983b, ).
This is known to lead to the estimation of a (quantile-specific) discrete
mixing distribution defined over the set of locations
$\\{\mbox{\boldmath$\zeta$}_{1,q},\dots,\mbox{\boldmath$\zeta$}_{G_{q},q}\\}$,
with mixture probabilities
$\pi_{g,q}=\Pr(\mbox{\boldmath$b$}_{i,q}=\mbox{\boldmath$\zeta$}_{g,q})$,
$i=1,\dots,n,\,g=1,\dots,G_{q}$, and $G_{q}\leq n$. Under this approach, for
$\mbox{\boldmath$b$}_{i,q}=\mbox{\boldmath$\zeta$}_{g,q}$, the location
parameter of the ALD in equation (1) becomes

$\mu_{itg,q}=\mbox{\boldmath$x$}_{it}^{\prime}\mbox{\boldmath$\beta$}_{q}+\mbox{\boldmath$z$}_{it}^{\prime}\mbox{\boldmath$\zeta$}_{g,q},$
(2)

while the model likelihood is defined by the following expression:

$L(\cdot\mid
q)=\prod_{i=1}^{n}\sum_{g=1}^{G_{q}}\left[\prod_{t=1}^{T_{i}}f_{y\mid
b}(y_{it}\mid\mbox{\boldmath$b$}_{i,q}=\mbox{\boldmath$\zeta$}_{g,q};q)\right]\pi_{g,q}.$
(3)

As it can be easily noticed, equation (3) resembles the likelihood of a finite
mixture of linear quantile regression models with TC, discrete, random
coefficients.

###  2.2 Linear quantile mixtures with TV random coefficients

While the model specification introduced in the previous paragraph accounts
for unit-specific omitted covariates, it may not be able to handle time
variations in such unobserved heterogeneity. To solve this issue Farcomeni,
(2012) introduced a linear quantile regression model with TV, discrete, random
intercepts. To ensure parsimony and interpretability of the model, such TV
unobserved heterogeneity is not left free to vary, but it is rather modeled
via a discrete-time Markov chain. While such a proposal entails unit-specific
intercepts only, a broader specification with general TV, discrete, random
coefficients can also be considered.

As before, let $\mbox{\boldmath$x$}_{it}=(x_{it1},\dots,x_{itp})^{\prime}$
denote a set of $p$ explanatory variables associated with the vector of
parameters $\mbox{\boldmath$\beta$}_{q}$. Further, let
$\mbox{\boldmath$w$}_{it}=(w_{it1},\dots,w_{itl})^{\prime}$ denote a set of
$l\geq 1$ explanatory variables not included in $\mbox{\boldmath$x$}_{it}$,
associated with a vector of unit-, time-, and quantile-specific random
coefficients
$\mbox{\boldmath$\alpha$}_{it,q}=(\alpha_{it1,q},\dots,\alpha_{itl,q})^{\prime}$.
These are assumed to evolve over time according to a homogeneous, first order,
hidden Markov chain $\\{S_{it,q}\\}$ that depends on the specified quantile
level $q$. In detail, $\\{S_{it,q}\\}$ is defined over the finite state space
$\mathscr{S}_{q}=\\{1,\dots,m_{q}\\}$, with initial and transition
probabilities $\delta_{s_{i1},q}$ and $\gamma_{s_{it}|s_{it-1},q}$, defined
as:

$\delta_{s_{i1},q}=\Pr(S_{i1,q}=s_{i1,q})$

and

$\gamma_{s_{it,q}|s_{it-1,q},q}=\Pr(S_{it,q}=s_{it,q}\mid
S_{it-1,q}=s_{it-1,q}).$

As for the model described in Section 2.1, local independence holds, together
with the convenient assumption of conditional ALD for the responses. For a
given quantile level $q\in(0,1)$, the joint conditional distribution of the
observed $T_{i}$-dimensional response vector
$\mbox{\boldmath$y$}_{i}=(y_{i1},\dots,y_{iT_{i}})^{\prime}$ is given by

$f_{y\mid s}(\mbox{\boldmath$y$}_{i}\mid\mbox{\boldmath$s$}_{i,q};q)=f_{y\mid
s}(y_{i1}\mid{s_{i1,q}};q)\prod_{t=2}^{T_{i}}f_{y\mid s}(y_{it}\mid
s_{it,q};q).$

Here, $\mbox{\boldmath$s$}_{i,q}=(s_{i1,q},\dots,s_{iT_{i},q})^{\prime}$ is
the vector of states visited by the $i$-th statistical unit and $f_{y\mid
s}(y_{it}\mid s_{it,q};q)$ denotes the ALD with skewness $q$, scale
$\sigma_{q}$, and location parameter $\mu_{its_{it,q},q}$ that is defined by

$\mu_{its_{it,q},q}=\mbox{\boldmath$x$}_{it}^{\prime}\mbox{\boldmath$\beta$}_{q}+\mbox{\boldmath$w$}_{it}^{\prime}\mbox{\boldmath$\alpha$}_{s_{it,q},q}.$

According to the assumptions above, the likelihood function is

$L(\cdot\mid
q)=\prod_{i=1}^{n}\sum_{\mbox{\footnotesize\boldmath{$s$}}_{i,q}}\left[\delta_{s_{i1},q}\prod_{t=2}^{T_{i}}\gamma_{s_{it,q}|s_{it-1,q}}\right]\left[\prod_{t=1}^{T_{i}}f_{y\mid
s}(y_{it}\mid{s_{it,q}};q)\right],$ (4)

where $\sum_{\mbox{\footnotesize\boldmath{$s$}}_{i}}$ is a shorthand for
$\sum_{s_{i1}}^{m_{q}}\cdots\sum_{s_{iT_{i}}=1}^{m_{q}}$. When looking at the
above expression, we may recognize it is a dynamic extension of the likelihood
reported in equation (3), as it represents the likelihood of a dynamic finite
mixture of linear quantile regression models with TV, discrete, random
coefficients. The coefficients are assumed to vary as a function of the hidden
states visited by the observed units over the observed time window.

###  2.3 Linear quantile mixtures with both TC and TV random coefficients

In some real data applications, both TC and TV sources of unit-specific
unobserved heterogeneity may be present and influence the response
distribution. We report in Figure 1 the longitudinal trajectories representing
the evolution of a given response measured over (at most) 6 time occasions for
a sample of 25 statistical units. From this figure, it is clear that both
sources of unobserved heterogeneity affect the response. TC unobserved
features may be responsible of differences between units in terms of baseline
levels and systematic temporal trends; TV unobserved features may instead
explain the sudden temporal shocks characterizing individual trajectories.
These latter may be rather difficult to capture via TC random coefficients or
unit-specific random slopes associated to a time variable. In these
situations, the linear quantile mixtures described above are no longer
appropriate, as they may account for one source at a time only. To model
empirical cases where both sources of between and within-unit variation are
present, Marino et al., (2018) introduced a linear quantile regression model
where TC and TV random coefficients may be jointly present in the linear
predictor.

Figure 1: An example of longitudinal trajectories affected by both time-
constant and time-varying unobserved heterogeneity

Let $\mbox{\boldmath$x$}_{it}=(x_{it1},\dots,x_{itp})^{\prime}$ indicate the
$p$-dimensional (design) vector of explanatory variables associated with the
parameters $\mbox{\boldmath$\beta$}_{q}$. Furthermore, let
$\mbox{\boldmath$z$}_{it}=(z_{it1},\dots,z_{itr})^{\prime}$ and
$\mbox{\boldmath$w$}_{it}=(w_{it1},\dots,w_{itl})^{\prime}$ denote two
disjoint vectors of $r\geq 1$ and $l\geq 1$ explanatory variables (not
included in $\mbox{\boldmath$x$}_{it}$), respectively. The former is
associated with the vector of unit- and quantile-specific random coefficients
$\mbox{\boldmath$b$}_{i,q}=(b_{i1,q},\dots,b_{ir,q})^{\prime}$ taking value in
the in the set
$\\{\mbox{\boldmath$\zeta$}_{1,q},\dots,\mbox{\boldmath$\zeta$}_{G_{q},q}\\}$
with probability
$\pi_{g,q}=\Pr(\mbox{\boldmath$b$}_{i,q}=\mbox{\boldmath$\zeta$}_{g,q}),g=1,\dots,G_{q}$.
The latter is associated to the vector of unit-, time-, and quantile-specific
random coefficients
$\mbox{\boldmath$\alpha$}_{it,q}=(\alpha_{it1,q},\dots,\alpha_{itl,q})^{\prime}$,
which evolves over time according to a homogeneous, first order, and quantile-
dependent, hidden Markov chain $\\{S_{it,q}\\}$. As before, this is defined
over the finite state space $\mathscr{S}_{q}=\\{1,\dots,m_{q}\\}$ and is fully
described by means of the initial probability vector
$\mbox{\boldmath$\delta$}_{q}=(\delta_{1,1},\dots,\delta_{m_{q},q})^{\prime}$
and the transition probability matrix $\mbox{\boldmath$\Gamma$}_{q}$, with
generic element $\gamma_{s_{it,q}\mid s_{it-1,q},q}.$ The reader must note
that the intercept term is included either in $\mbox{\boldmath$z$}_{it}$ or in
$\mbox{\boldmath$w$}_{it}$ (or in none of these vectors), but not in both.
Further, random coefficients $\mbox{\boldmath$\zeta$}_{g,q}$ and
$\mbox{\boldmath$\alpha$}_{s_{it,q},q}$ are assumed to be independent.

Conditional on $\mbox{\boldmath$b$}_{i,q}$ and
$\mbox{\boldmath$\alpha$}_{it,q}$, longitudinal responses recorded on the same
unit are assumed to be independent (local independence) and to follow an AL
distribution with skewness, scale, and location parameter denoted by $q$,
$\sigma_{q}$, and $\mu_{itgs_{it},q}$, respectively. In detail, for a given
quantile level $q\in(0,1)$ and
$\mbox{\boldmath$b$}_{i,q}=\mbox{\boldmath$\zeta$}_{g,q}$, the location
parameter is defined by the following regression model:

$\mu_{itgs_{it,q},q}=\mbox{\boldmath$x$}_{it}^{\prime}\mbox{\boldmath$\beta$}_{q}+\mbox{\boldmath$z$}_{it}^{\prime}\mbox{\boldmath$\zeta$}_{g,q}+\mbox{\boldmath$w$}_{it}^{\prime}\mbox{\boldmath$\alpha$}_{s_{it,q},q}.$

Based on such assumptions, the likelihood function is

$\displaystyle L(\cdot\mid q)=\prod_{i=1}^{n}$
$\displaystyle\sum_{g=1}^{G_{q}}\sum_{\mbox{\footnotesize\boldmath{$s$}}_{i}}\left[\delta_{s_{i1,q}}\prod_{t=2}^{T_{i}}\gamma_{s_{it,q}\mid
s_{it-1,q}}\right]\left[\prod_{t=1}^{T_{i}}f_{y\mid
b,s}(y_{it}\mid\mbox{\boldmath$b$}_{i,q}=\mbox{\boldmath$\zeta$}_{g,q},{s_{it,q}};q)\right]\pi_{g,q},$

where, as before, $\sum_{\mbox{\footnotesize\boldmath{$s$}}_{i}}$ is a
shorthand for $\sum_{s_{i1}}^{m_{q}}\cdots\sum_{s_{iT_{i}}=1}^{m_{q}}$ and
$f_{y\mid
b,s}(y_{it}\mid\mbox{\boldmath$b$}_{i,q}=\mbox{\boldmath$\zeta$}_{g,q},{s_{it,q}};q)$
denotes the AL density. In this framework, unobserved unit-specific features
that remain constant over time are captured by the random coefficients
$\mbox{\boldmath$\zeta$}_{g,q}$. Similarly, sudden temporal shocks in the
unit-specific profiles, due to TV sources of unobserved heterogeneity, are
captured by the random coefficients $\mbox{\boldmath$\alpha$}_{s_{it,q},q}$.

To conclude this section, it is worth noticing that, when a single hidden
state ($m_{q}=1$) is considered and
$\mbox{\boldmath$b$}_{i,q}=\mbox{\boldmath$\zeta$}_{g,q}$, the location
parameter $\mu_{itgs_{it,q},q}$ simplifies to $\mu_{itg,q}$ and the model
reduces to the linear quantile mixture with TC random coefficients only. Also,
when a single mixture component is considered ($G_{q}=1$), the location
parameter $\mu_{itgs_{it,q},q}$ simplifies to $\mu_{its_{it,q},q}$ and the
model above reduces to the dynamic finite mixture of linear quantile
regressions with TV random coefficients only. Last, when both $G_{q}$ and
$m_{q}$ are equal to 1, the model reduces to a standard linear quantile
regression model, without random coefficients. These properties make this
latter specification more flexible and general than the alternatives described
so far, at the cost of a higher computational complexity.

##  3 Model estimation and inference

In this section, we describe the algorithm for ML estimation of model
parameters in the linear quantile mixture models detailed in Section 2. We
focus on the specification based on both TC and TV random coefficients
(Section 2.3), as the simpler alternatives detailed in Sections 2.1 and 2.2
can be derived as special cases. We provide details of an extended EM
algorithm for parameter estimation in Section 3.1; the procedure for deriving
standard errors and choosing the optimal number of mixture components and/or
states of the hidden Markov chain is described in Section 3.2.

###  3.1 Maximum likelihood estimation

Let $\mbox{\boldmath$\theta$}_{q}$ denote the global set of free model
parameters for a given quantile level $q\in(0,1)$. To derive an estimate for
such a vector, we may rely on indirect maximization of the likelihood function
via an extended version of the EM algorithm (Dempster et al.,, 1977). Let
$u_{i,q}(g)$ and $v_{it,q}(h)$ be the indicator variables for
$\mbox{\boldmath$b$}_{i,q}=\mbox{\boldmath$\zeta$}_{g,q}$ and $S_{it,q}=h$,
respectively, and let $v_{it,q}(h,k)=v_{it-1,q}(h)\times v_{it,q}(k)$ denote
the indicator variable for unit $i$ moving from state $h$ at occasion $t-1$ to
state $k$ at occasion $t$, with $g=1,\dots,G_{q},$ and $h,k=1,\dots,m_{q}$.
For a given $q\in(0,1)$, the EM algorithm starts from the following complete-
data log-likelihood function:

$\displaystyle\ell_{c}(\mbox{\boldmath$\theta$}_{q})$
$\displaystyle=\sum_{i=1}^{n}\left\\{\left[\sum_{g=1}^{G_{q}}u_{i,q}(g)\log\pi_{g,q}\right]\right.$
(5)
$\displaystyle+\left.\left[\sum_{h=1}^{m_{q}}v_{i1,q}(h)\log\delta_{h,q}+\sum_{t=2}^{T_{i}}\sum_{h=1}^{m_{q}}\sum_{k=1}^{m_{q}}v_{it,q}(h,k)\log\gamma_{k\mid
h,q}\right]\right.$
$\displaystyle+\left.\left[\sum_{t=1}^{T_{i}}\sum_{g=1}^{G_{q}}\sum_{h=1}^{m_{q}}u_{i,q}(g)v_{it,q}(h)\log
f_{y\mid
b,s}(y_{it}\mid\mbox{\boldmath$b$}_{i,q}=\mbox{\boldmath$\zeta$}_{g,q},S_{it,q}=h)\right]\right\\}.$

At the $r$-th iteration, the E-step of the EM algorithm requires the
computation of the expected value of the complete-data log-likelihood in
equation (5), conditional on the observed data
$\mbox{\boldmath$y$}=(\mbox{\boldmath$y$}_{1},\dots,\mbox{\boldmath$y$}_{n})^{\prime}$
and the current parameter estimates $\mbox{\boldmath$\theta$}_{q}^{(r-1)}$.
That is, it requires the computation of

$Q(\mbox{\boldmath$\theta$}_{q}\mid\mbox{\boldmath$\theta$}_{q}^{(r-1)})=E\left[\ell_{c}(\mbox{\boldmath$\theta$}_{q})\mid\mbox{\boldmath$y$},\mbox{\boldmath$\theta$}_{q}^{(r-1)}\right].$

This means computing the posterior expectations of the indicator variables
$u_{i,q}(g)$, $v_{it,q}(h)$, and $v_{it,q}(h,k)$. As regards these latter two,
a simplification is obtained by considering the forward and backward
probabilities (Baum et al.,, 1970; Welch,, 2003) typically used in the hidden
Markov model framework; see Marino et al., (2018) for further details. In the
M-step of the algorithm, parameter estimates
$\hat{\mbox{\boldmath$\theta$}}_{q}$ are derived by maximizing
$Q(\mbox{\boldmath$\theta$}_{q}\mid\mbox{\boldmath$\theta$}_{q}^{(r-1)})$ with
respect to $\boldsymbol{\theta}_{q}$.

The E- and the M-step are alternated until convergence, which is defined as
the (relative) difference between two subsequent likelihood values being lower
than a given threshold, $\varepsilon>0.$

###  3.2 Standard errors and model selection

Following the standard procedure used in the quantile regression framework,
standard errors for model parameter estimates are derived by exploiting a
nonparametric bootstrap approach (see e.g., Buchinsky,, 1995). In detail, we
employ a block bootstrap procedure, where a resampling of the statistical
units is performed and the corresponding sequence of observed measurements is
retained to preserve within-unit dependence (Lahiri,, 1999).

Let $\hat{\mbox{\boldmath$\theta$}}^{(b)}_{q}$ denote the vector of model
parameter estimates obtained in the $b$-th bootstrap sample, $b=1,\dots,B$.
Estimates of the standard errors for the vector
$\hat{\mbox{\boldmath$\theta$}}_{q}$ correspond to the diagonal elements of
the matrix

$\hat{\textbf{V}}(\hat{\mbox{\boldmath$\theta$}}_{q})=\sqrt{\frac{1}{B}\sum_{b=1}^{B}\left(\hat{\mbox{\boldmath$\theta$}}^{(b)}_{q}-\hat{\mbox{\boldmath$\theta$}}_{q}\right)\left(\hat{\mbox{\boldmath$\theta$}}^{(b)}_{q}-\hat{\mbox{\boldmath$\theta$}}_{q}\right)^{\prime}}.$

A crucial point when dealing with finite mixtures is the choice of the number
of components and/or hidden states. For a fixed quantile level $q$, a simple
and frequently used solution is as follows: parameter estimates are computed
for varying combinations of the number of components and states,
$[G_{q},m_{q}]$, and the model with the best fit, typically measured via log-
likelihood or penalized likelihood criteria (such as AIC or BIC), is retained.
However, given that these criteria may suffer from early stopping due to lack
of progress rather than true convergence and the log-likelihood may present
multiple local maxima, the EM algorithm is initialized from different starting
points and the model corresponding to the highest likelihood value is retained
before applying the penalization.

##  4 The lqmix R package

In this section, we introduce the R package lqmix, developed to deal with
linear quantile mixture models for longitudinal data. We illustrate the main
functions for estimation and inference on the parameters of models described
in Sections 2.1-2.3 by considering the application to the labor pain benchmark
dataset (Davis,, 1991). This is described in detail in Section 4.1.

###  4.1 Labor pain data

Firstly reported by Davis, (1991) and since then analyzed by Jung, (1996) and
Geraci and Bottai, (2007) among others, labor pain data come from a randomized
clinical trial aiming at analyzing the effectiveness of a medication for
relieving labor pain in women. A total of $n=83$ women were randomized to a
treatment/placebo group, and a response variable evaluating the self-reported
amount of pain measured every $30$ minutes on a 100-$mm$ line was recorded.
Here, 0 corresponds to absence of pain, while 100 corresponds to extreme pain.
The number of available measurements per woman ranges from a minimum of 1 to a
maximum of $T_{i}=6$. That is, we are in the presence of an unbalanced
longitudinal design, for which a MAR assumption seems to be reasonable. A
total of $\sum_{i=1}^{83}T_{i}=357$ measurements are available.

Together with the outcome of interest (meas), two explanatory variables are
available: trt denotes an indicator variable identifying the group each woman
is assigned to (1 = treatment, 0 = placebo), while time identifies the
measurement occasion the response was recorded at. Data are stored in the data
frame pain included in the lqmix package, as shown below.

    > library(lqmix)
    > head(pain)
      id meas trt time
    1  1  0.0   1    1
    2  1  0.0   1    2
    3  2  0.0   1    1
    4  2  0.0   1    2
    5  2  0.0   1    3
    6  2  2.5   1    4

Data are severely skewed, and skewness changes magnitude and sign over time.
In Figure 2, we report some selected diagnostics for a linear mixed model
based on TC, Gaussian, random intercepts, using trt, time, and their
interaction as covariates.

    > library(lme4)
    > outLME = lmer(meas~trt+time+trt:time + (1|id), data=pain)
    > par(mfrow = c(1,2))
    > qqnorm(residuals(outLME), main = "")
    > qqline(residuals(outLME))
    > qqnorm(unlist(ranef(outLME)$id), main = "")
    > qqline(unlist(ranef(outLME)$id))

In particular, Figure 2 reports the Normal probability plots for model
residuals: the left panel shows unit- and time-specific residuals, while the
right one shows the empirical Bayes estimates of unit-specific random
intercepts. As it can be easily noticed, both plots indicate the presence of
potentially influential observations in the data, as well as the violation of
the Gaussian assumption for the random intercepts in the model. Therefore,
using linear quantile mixtures seems to be a reasonable choice for the
analysis of such data.

Figure 2: Labor pain data. Linear mixed effect model with Gaussian random
intercepts. Normal probability plot for unit- and time-specific residuals
(left plot) and for empirical Bayes estimates of unit-specific intercepts
(right plot).

###  4.2 The lqmix function

The main function implemented in the lqmix R package is lqmix. It allows for
the estimation of the parameters of a linear quantile mixture based on either
TC, TV or both types of discrete random coefficients. The input arguments can
be displayed on the R console as follows:

    > args(lqmix)
    function (formula, randomTC=NULL, randomTV=NULL, group, time, G=NULL,
            ΨΨ m=NULL, data, qtl=0.5, eps=10^-5, maxit=1000, se=TRUE, R=50,
           ΨΨ  start=0, parInit = list(betaf=NULL, betarTC=NULL, betarTV=
    ΨΨ  ΨΨΨΨ    Ψ  NULL, pg=NULL, delta=NULL, Gamma=NULL, scale=NULL),
    ΨΨ          verbose=TRUE, seed=NULL, parallel=FALSE)

The first mandatory argument, formula, denotes a two-side formula of the type
(resp $\sim$ Fexpr), where resp is the response variable and Fexpr is an
expression determining the fixed-coefficient vector,
$\mbox{\boldmath$x$}_{it}$. On the other side, randomTC and randomTV are two
one-side formulas of the type ($\sim$ Rexpr1) and ($\sim$ Rexpr2), where
Rexpr1 and Rexpr2 identify the columns of the TC and the TV random-coefficient
vectors, $\mbox{\boldmath$z$}_{it}$ and $\mbox{\boldmath$w$}_{it}$,
respectively. Note that both these arguments are optional, so that the user
may decide to estimate a linear quantile regression model with either TC or TV
random coefficients, or with a combination of them. Further, variables
reported in Rexpr1 and Rexpr2 must not overlap and, if any of them appears
also in the fixed effect formula (formula), the corresponding parameter is
estimated as a TC/TV random coefficient only. The arguments time and group are
strings indicating the time and the grouping variable, respectively. All such
variables are taken from the optional data frame data. The arguments G and m
are used to specify the number of mixture components and/or hidden states in
the model, respectively, while qtl allows to specify the quantile level (by
default, qtl = 0.5). The arguments se, R, and start allow to specify whether
bootstrap standard errors should be computed (by default, se = TRUE), the
number of bootstrap samples to be used for this purpose (by default, R = 50),
and the initialization rule to consider. As regards this latter, three
possible specifications are allowed: start = 0 is used for a deterministic
start of model parameters (the default option); start = 1 is used for a random
start of model parameters; start = 2 is used to consider a starting rule based
on given model parameters specified via the parInit list argument. The
arguments maxit, eps, and verbose identify the maximum number of iterations
for the EM algorithm (by default, maxit = 1000), the corresponding tolerance
level (by default eps = $\mathtt{10^{-5}}$), and whether output should be
printed (by default, verbose = TRUE). Finally, seed and parallel are devoted
to setting a seed for random number generation that may be useful when
considering a random starting rule (start=1) and for computing bootstrap
standard errors (only when se=TRUE) via parallel computing (2 cores are used),
respectively.

####  4.2.1 Estimating linear quantile mixtures with TC random coefficients

To explore features of the lqmix function, we start by considering a two-
component linear quantile mixture ($\mathtt{G_{q}=2}$) for the median
($\mathtt{qtl=0.5}$) with a TC random intercept for the analysis of the pain
data. This is estimated as follows:

    > outTC = lqmix(formula=meas~time+trt+trt:time,randomTC=~1,group="id",
    + ΨΨ ΨΨΨtime="time",G=2,data=pain)
    Model TC - qtl = 0.5
    ------------|-------------|-------------|-------------|
      iteration |      G      |      lk     |   (lk-lko)  |
    ------------|-------------|-------------|-------------|
              0 |           2 |    -1682.31 |          NA |
              8 |           2 |    -1607.11 | 4.64719e-06 |
    ------------|-------------|-------------|-------------|
    Computing standard errors: 1  ... 2  ... 3  ... 4  ... 5  ... 6  ... 7
    ... 8  ... 9  ... 10  ... 11  ... 12  ... 13  ... 14  ... 15  ... 16  ...
    17 ... 18 ... 19  ... 20  ... 21  ... 22  ... 23  ... 24  ... 25  ... 26
    ... 27  ... 28 ... 29  ... 30  ... 31  ... 32  ... 33  ... 34  ... 35  ...
    36 ... 37  ... 38 ... 39  ... 40  ... 41  ... 42  ... 43  ... 44  ... 45
     ... 46  ... 47  ... 48  ... 49  ... 50  ...

The running time for the above command is 1.091 seconds when run on an Apple
M1 architecture. This represents the architecture used for all the codes
illustrated in the paper.

In the following, we report the output of the above call to lqmix, which is an
object of class “lqmix”, obtained by using the print method of the S3 class.

    > outTC
    Model: TC random coefficients with G=2 at qtl=0.5
    ******************************************************

    ---- Observed process ----

    Fixed Coefficients:
        time      trt time.trt
     11.6669  -2.1659 -10.8335

    Time-Constant Random Coefficients:
          (Intercept)
    Comp1      1.3325
    Comp2     43.3325

    Residual scale parameter: 7.3289 - Residual standard deviation: 20.7294

    ---- Latent process ----

    Mixture probabilities:
     Comp1  Comp2
    0.6737 0.3263

    Log-likelihood at convergence: -1607.11
    Number of observations: 357 - Number of subjects: 83

Looking at the output, we may recognize two separate sections, reporting
estimates for the observed and the latent process, respectively. As regards
the observed process, the estimated fixed coefficients suggest how self-
reported pain increases as the time passes by, together with a positive effect
of the treatment under investigation (negative sign for the trt variable),
with benefits that increase with time. On the other side, the estimated random
coefficients highlight the presence of two well-separated groups of women,
reporting a low and a medium pain at the baseline, respectively. In the last
part, the estimated scale parameter and the corresponding error standard
deviation are shown. This latter corresponds to the standard deviation of an
ALD with parameters $\sigma_{q}$ and $q$:
$\frac{s^{2}(1-2q+2q^{2})}{(1-q)^{2}q^{2}}$. As regards the latent process,
estimates highlight that the 67.37% of women belong to the first mixture
component (low-pain level); the remaining 32.63% belong to the second one
(medium-pain level). All estimated parameters are stored in the object outTC
as betaf (the fixed coefficientes), betarTC (the TC random coefficients),
scale (the scale parameter), sigma.e (the conditional standard deviation of
responses derived from an ALD with parameters scale and qtl), and pg (the
component prior probabilities). Further, when leaving the statement se = TRUE
(default value), the outTC object also contains information on the estimated
standard errors of model parameters obtained via the bootstrap procedure
detailed in Section 3.2. In particular, these may be accessed by using the
prefix “se.” with model parameter names indicated above.

The following information are also stored in the the outTC object: the log-
likelihood value at convergence (lk), the number of model parameters (npar),
the values of AIC and BIC (aic and bic), the quantile level (qtl), the number
of mixture components (G), the EM threshold (eps), the total number of
subjects and observations (nsbjs and nobs), the type of missingness data are
affected by (miss), the estimated model (model), and the model call (call).

####  4.2.2 Estimating linear quantile mixtures with TV random coefficients

The same lqmix function can be used to obtain parameter estimates for a
quantile mixture with TV random coefficients. To analyze the labor pain data,
we consider a liner quantile mixture for the median ($\mathtt{qtl=0.5}$), with
a TV random intercept defined over a two-state latent space
($\mathtt{m_{q}=2}$):

    > outTV = lqmix(formula=meas~time+trt+trt:time,randomTV=~1,start=0,
    +ΨΨΨ   group="id", time="time",m=2,data=pain)
    Model TV - qtl = 0.5
    ------------|-------------|-------------|-------------|
      iteration |      m      |      lk     |   (lk-lko)  |
    ------------|-------------|-------------|-------------|
              0 |           2 |    -1688.99 |          NA |
             10 |           2 |    -1576.61 |    0.366345 |
             20 |           2 |    -1575.79 | 9.32048e-05 |
             29 |           2 |    -1575.79 | 9.43243e-06 |
    ------------|-------------|-------------|-------------|
    Computing standard errors: 1  ... 2  ... 3  ... 4  ... 5  ... 6  ... 7
    ... 8  ... 9  ... 10  ... 11  ... 12  ... 13  ... 14  ... 15  ... 16
    ...  17  ... 18  ... 19  ... 20  ... 21  ... 22  ... 23  ... 24  ...
    25  ... 26  ... 27  ... 28  ... 29  ... 30  ... 31  ... 32  ... 33
    ... 34  ... 35  ...   36  ... 37  ... 38  ... 39  ... 40  ... 41  ...
    42  ... 43  ... 44  ... 45  ... 46  ... 47  ... 48  ... 49  ... 50  ...

The running time for obtaining results is 4.933 seconds. The output of
estimation is given below.

    > outTV
    Model: TV random coefficients with m=2 at qtl=0.5
    ******************************************************

    ---- Observed process ----

    Fixed Coefficients:
        time      trt time.trt
      6.5012  -0.4923  -6.0013

    Time-Varying Random Coefficients:
        (Intercept)
    St1      0.4924
    St2     60.9926

    Residual scale parameter: 5.8337 - Residual standard deviation: 16.5003

    ---- Latent process ----

    Initial probabilities:
       St1    St2
    0.7667 0.2333

    Transition probabilities:
             toSt1  toSt2
    fromSt1 0.8883 0.1117
    fromSt2 0.0271 0.9729

    Log-likelihood at convergence: -1575.795
    Number of observations: 357 - Number of subjects: 83

Results for the observed process allow to derive similar conclusions to those
detailed in Section 4.2.1, even though the magnitude of the effects is
reduced.

As regards the latent layer of the model, the print method shows the estimates
for the initial and the transition probabilities of the hidden Markov chain.
Based on such estimates, one may conclude that 76.67% of women start in the
first hidden state which, based on the estimated time-varying random
intercepts, is characterized by a low pain level. The remaining 23.33% of
women start the study with an intermediate baseline pain level. Looking at the
estimated transition probabilities, we may notice that, regardless of the
treatment effect, labor pain tends to increase as the time passes by, with
women being in the low pain state moving towards the medium pain state with
probability equal to $0.1117$. Estimated parameters are stored in the object
outTV as betaf (the fixed coefficients), betarTV (the TV random coefficients),
scale (the scale parameter), sigma.e (the conditional standard deviation of
responses derived from an ALD with parameters scale and qtl), delta (the
initial probability vector), and Gamma (the transition probability matrix). As
detailed in the previous section, bootstrap standard errors are stored by
using the prefix “se.” to the names of parameters of interest.

####  4.2.3 Estimating linear quantile mixtures with TC and TV random
coefficients

When specifying both the randomTC and the randomTV formulas, and therefore
both G and m, the lqmix function allows for the estimation of linear quantile
mixtures based on both TC and TV, discrete, random coefficients. To analyze
labor pain data, we consider a linear quantile mixture for the median (qtl =
0.5), with a TV random intercept, with $m_{q}=2$ hidden states, and a TC
random slope for the time variable, with $G_{q}=2$ components. To estimate
such a model and consider a random start initialization (start=1), the
following command can be run:

    > outTCTV = lqmix(formula=meas~trt+time+trt:time,randomTC=~time,
    +.   randomTV=~1, group="id",time="time",m=2,G=2,data=pain, se=TRUE,
    +    start=1, seed=11)
    Model TCTV - qtl = 0.5
    ------------|-------------|-------------|-------------|-------------|
      iteration |      m      |      G      |      lk     |   (lk-lko)  |
    ------------|-------------|-------------|-------------|-------------|
              0 |           2 |           2 |    -1724.08 |          NA |
             10 |           2 |           2 |     -1591.1 |     11.8652 |
             20 |           2 |           2 |    -1554.53 |     2.06191 |
             30 |           2 |           2 |    -1545.38 |     1.47303 |
             40 |           2 |           2 |    -1540.69 |  0.00400722 |
             44 |           2 |           2 |    -1540.69 | 7.06303e-06 |
    ------------|-------------|-------------|-------------|-------------|
    Computing standard errors: 1  ... 2  ... 3  ... 4  ... 5  ... 6  ... 7
    ... 8  ... 9  ... 10  ... 11  ... 12  ... 13  ... 14  ... 15  ... 16
    ... 17 ... 18  ... 19  ... 20  ... 21  ... 22  ... 23  ... 24  ... 25
      ... 26  ... 27  ... 28  ... 29  ... 30  ... 31  ... 32  ... 33  ...
      34  ... 35  ... 36  ... 37  ... 38  ... 39  ... 40  ... 41  ... 42
      ... 43  ... 44  ... 45  ... 46  ... 47  ... 48  ... 49  ... 50  ...

This requires a running time of 8.863 seconds. The output of the above call to
lqmix, obtained through the print method of the S3 class, is reported below.

    > outTCTV
    Model: TC and TV random coefficients with m=2 and G=2 at qtl=0.5
    ******************************************************************

    ---- Observed process ----

    Fixed Coefficients:
         trt trt.time
     -8.9867  -5.1027

    Time-Constant Random Coefficients:
             time
    Comp1  5.6027
    Comp2 14.9777

    Time-Varying Random Coefficients:
        (Intercept)
    St1      8.4867
    St2     69.4867

    Residual scale parameter: 5.1196 - Residual standard deviation: 14.4803

    ---- Latent process ----

    Mixture probabilities:
     Comp1  Comp2
    0.7376 0.2624

    Initial probabilities:
       St1    St2
    0.8029 0.1971

    Transition probabilities:
             toSt1  toSt2
    fromSt1 0.9643 0.0357
    fromSt2 0.0437 0.9563

    Log-likelihood at convergence: -1540.693
    Number of observations: 357 - Number of subjects: 83

As it is clear from the output, the linear quantile mixture with both TC and
TV, discrete, random coefficients produces more refined information when
compared to those discussed in Sections 4.2.1 and 4.2.2. Also in this case, we
may recognize in the output the two sections entailing the observed and the
latent process, respectively. In the former, estimates for the fixed, the TC,
and the TV random coefficients are reported. For the latent part, the output
reports information on the mixture, the initial, and the transition
probabilities.

By looking at the results, we conclude again that self reported pain is lower
for women under medication (negative sign for the trt variable) and that
benefits increase with time. Looking at the estimated TC random coefficients
and the corresponding prior probabilities, we may distinguish two well
separated groups of women. These are characterized by a lower and a higher
linear effect of the time on pain levels: for the 73.76% of the women under
observation, pain level mildly increases in a linear manner as the time passes
by, while the remaining 26.24% experience a steeper increase of pain levels
with time. On the other side, the estimated TV random coefficients identify
two groups of women characterized by a low and a medium baseline labor pain
level, respectively. At the beginning of the observation period, the 80.29% of
women belong to the first group, while the remaining 19.71% to the second. By
looking at the estimated transition probability matrix, we do conclude that
group composition remains rather unchanged over time once controlling for the
other effects in the model ($\gamma_{hh}>0.95,h=1,2$).

As reported above, estimated parameters are stored in the object outTCTV. Now
both those related to the TC finite mixture as well as those related to the
hidden Markov chain are referenced. In detail, betaf (the fixed coefficients),
betarTC (the TC random coefficients), betarTV (the TV random coefficients,
scale (the scale parameter), sigma.e (the conditional standard deviation of
responses derived from an ALD with parameters scale and qtl), pg (the
component prior probabilities), delta (the initial probabilities), and Gamma
(the transition probabilities). Also, standard errors (if computed),as well as
additional information on the data and the estimated model, are stored in the
outTCTV object.

###  4.3 The summary method for lqmix objects

The summary method applied to objects of the class lqmix produces a summary
object that helps make inference on model parameters. As stated above,
calculation of standard errors is based on a block-bootstrap procedure,
considering by default $R=50$ replicates. The output of the summary method for
the outTCTV object described in Section 4.2 is presented below.

    > summary(outTCTV)
    Model: TC and TV random coefficients with m=2 and G=2 at qtl=0.5
    ***********************************************************************

    ---- Observed process ----

    Fixed Coefficients:
             Estimate St.Error z.value P(>|z|)
    trt       -8.9867   5.4057 -1.6624  0.0945 .
    trt.time  -5.1027   0.8940 -5.7076  <2e-16 ***

    Time-Constant Random Coefficients:
               Estimate St.Error z.value   P(>|z|)
    time_Comp1   5.6027   0.5475  10.234 < 2.2e-16 ***
    time_Comp2  14.9777   0.7552  19.832 < 2.2e-16 ***

    Time-Varying Random Coefficients:
                    Estimate St.Error z.value P(>|z|)
    (Intercept)_St1   8.4867   3.5206  2.4106  0.0156 *
    (Intercept)_St2  69.4867   5.2791 13.1626  <2e-16 ***
    ---
    Signif. codes:  0 *** 0.001 ** 0.01 * 0.05 . 0.1   1

    Residual scale parameter: 5.1196 - Residual standard deviation: 14.4803

    ---- Latent process ----

    Mixture probabilities:
          Estimate St.Err
    Comp1   0.7376 0.0678
    Comp2   0.2624 0.0678

    Initial probabilities:
        Estimate St.Err
    St1   0.8029 0.0546
    St2   0.1971 0.0546

    Transition probabilities:
                 Estimate St.Err
    fromSt1toSt1   0.9643 0.0154
    fromSt1toSt2   0.0357 0.0154
    fromSt2toSt1   0.0437 0.0511
    fromSt2toSt2   0.9563 0.0511

    Log-likelihood at convergence: -1540.693
    Number of observations: 357 - Number of subjects: 83

As before, two different sections may be distinguished. In the former,
estimates, standard errors, test statistics, and corresponding approximate
p-values for assessing significance of model parameters are reported. These
information are completed with the estimates of the scale parameter and the
corresponding conditional standard deviation of the error terms. In the latter
section of the output, estimates and standard errors for the parameters in the
latent layer of the model are reported.

###  4.4 The search_lqmix function for model selection

As described in Section 3.2, the number of mixture components $G_{q}$ and the
number of hidden states $m_{q}$ in the model are unknown quantities that need
to be estimated. Moreover, a multi-start strategy is frequently needed to
solve, at least partially, the potential multimodality of the likelihood
surface. Both issues can be addressed by means of the search_lqmix() function.
Input arguments and corresponding default values can be displayed on the R
console though the args function:

    > args(search_lqmix)
    function (formula, randomTC=NULL, randomTV=NULL, group, time, Gv=NULL,
              mv=NULL, data, method="bic", nran=0, qtl=0.5, eps=10^-5,
              maxit=1000, se=TRUE, R=50, verbose=TRUE, seed=NULL)

Most of the arguments of this function correspond to those described so far.
The remaining ones are specified as follows. Gv and mv denote vectors
identifying the range for the number of mixture components, $G_{q}$, and the
number of hidden states, $m_{q}$, to be considered, respectively, for a fixed
quantile level $q\in(0,1)$. When both arguments are specified, the search of
the optimal linear quantile mixture with TC and TV random coefficients is
performed. When only one out the two arguments is specified, a linear quantile
mixture based on either TC or TV random coefficients is considered. A linear
quantile regression with no random coefficients is also estimated either when
both mv and/or Gv include the value 1, or when Gv includes 1 and mv = NULL, or
when mv includes 1 and Gv = NULL. In this case, the function lqr() implemented
in the lqmix R package and described in the following is employed.

The argument method is used to specify the method to be used for model
selection purposes. One of three possible values is admitted: “bic” (by
default), “aic”, or “lk”. The former two identify the optimal model as that
providing the minimum value of the BIC or the AIC index, respectively. The
latter, selects the optimal model as the one corresponding to the maximum log-
likelihood value.

The argument nran specifies the number of random starts to be considered for
each value in the range identified by Gv and mv. This number is strictly
related to the type of model for which the optimal search is performed. In
detail, following the same strategy suggested by Bartolucci et al., (2017) for
the LMest R package, when a linear quantile mixture based on TC random
coefficients only is considered, the number of random initializations is equal
to $\mathtt{nran\times(G_{q}-1)}$; when only TV random coefficients are
allowed, the number of random initializations is equal to
$\mathtt{nran\times(m_{q}-1)}$; last, when both TC and TV random coefficients
are considered, the number of random initializations is equal to
$\mathtt{nran\times(G_{q}-1)\times(m_{q}-1)}$. By default, $\mathtt{nran=0}$,
so that no random initializations are considered. Last the seed argument is
used to fix a seed and ensure reproducibility of results when a multi-start
strategy based on random initializations is considered to estimate model
parameters.

We report below the results of the model selection strategy for the pain data,
when focusing on a linear quantile mixture with a TV random intercept and a TC
random slope associated to the variable time, for the quantile level $q=0.5$.
The search is done by looking for the optimal number of components ($G_{q}$)
and hidden states ($m_{q}$), both in the set $\\{1,2\\}$, by setting nran=2.
The optimal model is selected according to the BIC index.

    > sTCTV = search_lqmix(formula=meas~trt+time+trt:time,randomTC=~time,
    +   randomTV=~1,nran=2, group="id",time="time",mv=1:2,Gv=1:2,data=pain,
    +   seed=11)
    Search the optimal linear quantile mixture model
    *************************************************
    Model homogeneous - qtl = 0.5
    Random start: 0 ...
    ------------|-------------|
      iteration |      lk     |
    ------------|-------------|
              0 |    -1707.73 |
    ------------|-------------|
    Model TC - qtl = 0.5
    Random start: 0 ... 1  ... 2  ...
    ------------|-------------|-------------|-------------|
      iteration |      G      |      lk     |   (lk-lko)  |
    ------------|-------------|-------------|-------------|
              0 |           2 |     -1607.5 |          NA |
              1 |           2 |     -1607.5 | 4.03395e-08 |
    ------------|-------------|-------------|-------------|
    Model TV - qtl = 0.5
    Random start: 0 ... 1  ... 2  ...
    ------------|-------------|-------------|-------------|
      iteration |      m      |      lk     |   (lk-lko)  |
    ------------|-------------|-------------|-------------|
              0 |           2 |    -1575.79 |          NA |
              5 |           2 |    -1575.79 | 8.44501e-06 |
    ------------|-------------|-------------|-------------|
    Model TCTV - qtl = 0.5
    Random start: 0 ... 1  ... 2  ...
    ------------|-------------|-------------|-------------|-------------|
      iteration |      m      |      G      |      lk     |   (lk-lko)  |
    ------------|-------------|-------------|-------------|-------------|
              0 |           2 |           2 |    -1542.77 |          NA |
              5 |           2 |           2 |    -1542.77 | 4.09907e-06 |
    ------------|-------------|-------------|-------------|-------------|
    Computing standard errors for the optimal model: 1  ... 2  ... 3  ...
    4  ... 5  ... 6  ... 7  ... 8  ... 9  ... 10  ... 11  ... 12  ... 13
    ... 14  ... 15  ... 16  ... 17  ... 18  ... 19  ... 20  ... 21  ... 22
    ... 23  ... 24  ... 25  ... 26  ... 27  ... 28  ... 29  ... 30  ... 31
    ... 32  ... 33  ... 34  ... 35  ... 36  ... 37  ... 38  ... 39  ... 40
    ... 41  ... 42  ... 43  ... 44  ... 45  ... 46  ... 47  ... 48  ... 49
    ... 50  ...

The running time for the above command is 11.152 seconds and the output is
stored in the sTCTV object. This can be shown by means of the print method of
the class S3 as follows:

    >  sTCTV
    Opt model: TC and TV random coefficients with m=2 and G=2 at qtl=0.5
    ***************************************************************************

    ---- Observed process ----

    Fixed Coefficients:
             Estimate St.Error z.value P(>|z|)
    trt      -10.4558   5.2090 -2.0072  0.0438 *
    trt.time  -4.8088   0.8036 -5.9845  <2e-16 ***

    Time-Constant Random Coefficients:
               Estimate St.Error z.value   P(>|z|)
    time_Comp1   5.3088   0.4938  10.752 < 2.2e-16 ***
    time_Comp2  14.6814   0.6614  22.198 < 2.2e-16 ***

    Time-Varying Random Coefficients:
                    Estimate St.Error z.value P(>|z|)
    (Intercept)_St1   9.9558   3.3033  3.0139  0.0025 **
    (Intercept)_St2  70.6372   5.2557 13.4401  <2e-16 ***
    ---
    Signif. codes:  0 *** 0.001 ** 0.01 * 0.05 . 0.1   1

    Residual scale parameter: 5.1482 - Residual standard deviation: 14.5614

    ---- Latent process ----

    Mixture probabilities:
          Estimate St.Err
    Comp1   0.7421 0.0716
    Comp2   0.2579 0.0716

    Initial probabilities:
        Estimate St.Err
    St1   0.8047 0.0531
    St2   0.1953 0.0531

    Transition probabilities:
                 Estimate St.Err
    fromSt1toSt1   0.9625 0.0173
    fromSt1toSt2   0.0375 0.0173
    fromSt2toSt1   0.0417 0.0513
    fromSt2toSt2   0.9583 0.0513

    Log-likelihood at convergence: -1542.77
    Number of observations: 357 - Number of subjects: 83

As it is evident, such an output corresponds to that of the summary.lqmix
method, applied to the optimal model specification identified via the BIC
index. This, as reported in the first output line, corresponds to the model
with $G_{q}=2$ mixture components and $m_{q}=2$ hidden states. Clearly, some
differences with respect to the output reported in Section 4.3 are present,
due to the random start initialization. It is worth to highlight that, when
setting se = FALSE in search_lqmix, the print method of the S3 class is
directly used to show the the results obtained under the identified optimal
model specification.

##  5 Other methods for lqmix and search_lqmix objects

Together with the print and the summary methods introduced in the previous
sections, basic functions and methods typically available for regression model
objects are implemented in the lqmix R package. These include the logLik and
the coef functions returning the maximum log-likelihood value obtained at
convergence of the EM algorithm and the estimated fixed effect coefficients.
Further, the methods AIC and BIC may be used for model selection purposes.
When the above methods are applied to objects of class search_lqmix, these
refer to the optimal specification identified via the search_lqmix function.
Last, the plot method returns graphical representations for the outputs of the
lqmix and the search_lqmix functions. In detail, a graphical display of the
mixture component probabilities and/or of the transitions across states of the
hidden Markov chain is provided when applied to an lqmix object; the plot is
similar in spirit to those provided in the LMest R package (Bartolucci et
al.,, 2017). A plot of the values of the chosen model selection criterion for
varying number of components and/or hidden states (defined in Gv and mv,
respectively) is also provided when applied to objects of class search_lqmix.

##  6 Linear quantile regression for independent data

The package lqmix is thought for dealing with longitudinal data. However, it
also implements the function lqr() for estimating a homogeneous linear
quantile regression model via a ML approach. This is obtained by exploiting
the parallelism between the AL distribution and the asymmetric quantile loss
function. Input arguments for the lqr function and corresponding default
values are as follows:

    > args(lqr)
    function (formula, data = sys.frame(sys.parent()), qtl=0.5, se=TRUE,
              R=50, verbose=TRUE)

These arguments are identical to those described above for the functions that
are built to estimate parameters of linear quantile mixtures. As evident, also
in this case the user may specify whether standard errors need to be computed.
In this case, a block bootstrap approach is considered. By defaults, R = 50
resamples are considered.

The print and the summary method of the S3 class can be used for object of
class lqr in order to print the output and make inference on model parameters.

##  7 Conclusions

Quantile regression is a fundamental tool of analysis when (i) the response
distribution is skewed, (ii) outilying values are present, (iii) interest
entails both the center and the tails of the response distribution. When
dealing with longitudinal data, dependence between observations from the same
subject must be properly considered to avoid biased inference. Thus, linear
quantile regression models with unit-specific random coefficients may be
effectively employed.

In this paper, we present the lqmix R package to estimate mixtures of linear
quantile regression models for continuous longitudinal data, possibly affected
by random missingness. Strong and unverifiable parametric assumptions on the
random coefficients are avoided by considering a mixture specification. Either
TC or TV random coefficients, or a combination of both can be fitted, to deal
with different sources of unobserved heterogeneity affecting the longitudinal
process. That is, the package allows us to deal with unit-specific unobserved
features (omitted covariates in the model) which may vary and/or stay constant
over the observational period, while looking at the quantiles of the
conditional response variable.

Further updates of the package will include the possibility to deal with
multivariate responses, as well as with binary or ordinal responses. The
inclusion of observed covariates on the latent layer of the model represents a
further extension we aim at working on.

## Acknowledgements

The work of Ranalli has been developed under the support of “Fondo Ricerca di
Ateneo, Università degli Studi di Perugia, edition 2021, project: AIDMIX”’.
The work of Marino has been supported by the grant “Dipartimenti Eccellenti
2023-2027”, from Italian Ministry of Education, University and Research. The
work of Salvati has been supported by the grant “PRIN 2022 PNRR Prot.
P2022TB5JF: Quantification in the Context of Dataset Shift (QuaDaSh)” from
Italian Ministry of Education, University and Research.

## References

  * Abrevaya and Dahl, (2008) Abrevaya, J. and Dahl, C. M. (2008).  The effects of birth inputs on birthweight: Evidence from quantile estimation on panel data.  Journal of Business and Economic Statistics, 26:379–397.
  * Alfó et al., (2017) Alfó, M., Salvati, N., and Ranalli, M. G. (2017).  Finite mixtures of quantiles and M-quantile models.  Statistics and Computing, 27:547–570.
  * Baker, (2016) Baker, M. (2016).  QREGPD: Stata module to perform Quantile Regression for Panel Data.  Statistical Software Components, Boston College Department of Economics.
  * Bartolucci et al., (2013) Bartolucci, F., Farcomeni, A., and Pennoni, F. (2013).  Latent Markov Models for Longitudinal Data.  Chapman & Hall/CRC.
  * Bartolucci et al., (2017) Bartolucci, F., Pandolfi, S., and Pennoni, F. (2017).  LMest: An R package for latent markov models for longitudinal categorical data.  Journal of Statistical Software, 81:1–38.
  * Bates et al., (2015) Bates, D., Mächler, M., Bolker, B., and Walker, S. (2015).  Fitting linear mixed-effects models using lme4.  Journal of Statistical Software, 67:1–48.
  * Baum et al., (1970) Baum, L., Petrie, T., Soules, G., and Weiss, N. (1970).  A maximization technique occurring in the statistical analysis of probabilistic functions of Markov chains.  Annals of Mathematical Statistics, 41:164–171.
  * Buchinsky, (1995) Buchinsky, M. (1995).  Estimating the asymptotic covariance matrix for quantile regression models. A Monte Carlo study.  Journal of Econometrics, 68:303–338.
  * Davis, (1991) Davis, C. S. (1991).  Semi-parametric and non-parametric methods for the analysis of repeated measurements with applications to clinical trials.  Statistics in Medicine, 10:1959–1980.
  * Dempster et al., (1977) Dempster, A., Laird, N. M., and Rubin, D. B. (1977).  Maximum likelihood from incomplete data via the EM algorithm.  Journal of the Royal Statistical Society. Series B. Methodological, 39:1–38.
  * Diggle et al., (2002) Diggle, P. J., Heagerty, P. J., Liang, K., and Zeger, S. L. (2002).  Analysis of Longitudinal Data, volume 25 of Oxford Statistical Science Series.  Oxford University Press, second edition.
  * Einbeck et al., (2018) Einbeck, J., Darnell, R., and Hinde, J. (2018).  npmlreg: Nonparametric Maximum Likelihood Estimation for Random Effect Models.  R package version 0.46-5.
  * Farcomeni, (2012) Farcomeni, A. (2012).  Quantile regression for longitudinal data based on latent Markov subject-specific parameters.  Statistics and Computing, 22.
  * Fitzmaurice et al., (2004) Fitzmaurice, G. M., Laird, N. M., and Ware, J. H. (2004).  Applied Longitudinal Analysis.  Wiley-Interscience, Hoboken, NJ.
  * Galvao and Wang, (2015) Galvao, A. F. and Wang, L. (2015).  Efficient minimum distance estimator for quantile regression fixed effects panel data.  Journal of Multivariate Analysis, 133:1–26.
  * Geraci, (2014) Geraci, M. (2014).  Linear quantile mixed models: the lqmm package for laplace quantile regression.  Journal of Statistical Software, 57:1–29.
  * Geraci and Bottai, (2007) Geraci, M. and Bottai, M. (2007).  Quantile regression for longitudinal data using the asymmetric laplace distribution.  Biostatistics, 8:140–54.
  * Geraci and Bottai, (2014) Geraci, M. and Bottai, M. (2014).  Linear quantile mixed models.  Statistics and Computing, 24:461–479.
  * Graham et al., (2015) Graham, B. S., Hahn, J., Poirier, A., and Powell, J. L. (2015).  Quantile regression with panel data.  Technical report, National Bureau of Economic Research.
  * Hao et al., (2007) Hao, L., Naiman, D. Q., and Naiman, D. Q. (2007).  Quantile Regression.  Sage.
  * Jung, (1996) Jung, S.-H. (1996).  Quasi-likelihood for median regression models.  Journal of the American Statistical Association, 91:251–257.
  * Koenker, (2004) Koenker, R. (2004).  Quantile regression for longitudinal data.  Journal of Multivariate Analysis, 91:74–89.
  * Koenker, (2005) Koenker, R. (2005).  Quantile regression.  Cambridge University Press.
  * Koenker, (2022) Koenker, R. (2022).  quantreg: Quantile Regression.  R package version 5.94.
  * Koenker and Bache, (2014) Koenker, R. and Bache, S. H. (2014).  rqpd: Regression Quantiles for Panel Data.  R package version 0.6/r10.
  * Koenker and Bassett, (1978) Koenker, R. and Bassett, Jr., G. (1978).  Regression quantiles.  Econometrica, 46:33–50.
  * Koenker and Hallock, (2001) Koenker, R. and Hallock, K. (2001).  Quantile regression.  Journal of Economic Perspectives, 15:143–156.
  * Lahiri, (1999) Lahiri, S. (1999).  Theoretical comparisons of block bootstrap methods.  Annals of Statistics, 27:386–404.
  * Laird, (1978) Laird, N. (1978).  Nonparametric maximum likelihood estimation of a mixing distribution.  Journal of the American Statistical Association, 73:805–811.
  * (30) Lindsay, B. G. (1983a).  The geometry of mixture likelihoods: a general theory.  The Annals of Statistics, 11:86–94.
  * (31) Lindsay, B. G. (1983b).  The Geometry of Mixture Likelihoods, Part II: the Exponential Family.  The Annals of Statistics, 11:783–792.
  * Little and Rubin, (2002) Little, R. J. and Rubin, D. B. (2002).  Statistical Analysis with Missing Data.  Wiley.
  * Machado et al., (2021) Machado, J. A., Parente, P., and Santos Silva, J. (2021).  Qreg2: Stata module to perform quantile regression with robust and clustered standard errors.
  * Machado and Santos Silva, (2019) Machado, J. A. and Santos Silva, J. (2019).  Quantiles via moments.  Journal of Econometrics, 213:145–173.
  * Machado and Santos Silva, (2021) Machado, J. A. and Santos Silva, J. (2021).  Xtqreg: Stata module to compute quantile regression with fixed effects.
  * Marino et al., (2023) Marino, M. F., Alfó, M., Salvati, N., and Ranalli, M. G. (2023).  lqmix: Linear Quantile Mixture Models.  R package version 0.1.0.
  * Marino and Farcomeni, (2015) Marino, M. F. and Farcomeni, A. (2015).  Linear quantile regression models for longitudinal experiments: an overview.  METRON, 73:229–247.
  * Marino et al., (2018) Marino, M. F., Tzavidis, N., and Alfò, M. (2018).  Mixed hidden markov quantile regression models for longitudinal data with possibly incomplete sequences.  Statistical Methods in Medical Research, 27:2231–2246.
  * Pons and Melly, (2022) Pons, M. and Melly, B. (2022).  Stata commands to estimate quantile regression with panel and grouped data.  Swiss Stata Conference 2022 05, Stata Users Group.
  * R Core Team, (2019) R Core Team (2019).  R: A Language and Environment for Statistical Computing.  R Foundation for Statistical Computing, Vienna, Austria.
  * Rubin, (1976) Rubin, D. B. (1976).  Inference and missing data.  Biometrika, 63:581–592.
  * Welch, (2003) Welch, L. R. (2003).  Hidden Markov models and the Baum-Welch algorithm.  IEEE Information Theory Society Newsletter, 53:10–13.
  * Yu et al., (2003) Yu, K., Lu, Z., and Stander, J. (2003).  Quantile regression: Applications and current research areas.  Journal of the Royal Statistical Society, Series D, 52:331–350.
  * Yu and Moyeed, (2001) Yu, K. and Moyeed, R. A. (2001).  Bayesian quantile regression.  Statistics and Probability Letters, 54:437–447.
  * Zucchini and MacDonald, (2009) Zucchini, W. and MacDonald, I. (2009).  Hidden Markov models for time series, volume 110 of Monographs on Statistics and Applied Probability.  CRC Press.


# Experimental prospects for $V_{ud}$, $V_{us}$, $V_{cd}$, $V_{cs}$ and
(semi-)leptonic decays at LHCb

Adam C. S. Davis
on behalf of the LHCb Collaboration
School of Physics and Astronomy
University of Manchester
Manchester, United Kingdom

(January 15, 2019)

###### Abstract

We review the prospects of measurements of the CKM matrix elements
corresponding to the first two generations of quarks, as well as measurements
of leptonic and semileptonic charm meson decays at LHCb. We review the status
of searches for Lepton Non-Universality in the charm sector, and provide
motivation to continue such searches.
Presented at the 10th workshop on the CKM unitary triangle, Heidelberg,
Germany

##  1 Introduction

The LHCb experiment at CERN is a single arm forward spectrometer specifically
built for the reconstruction of heavy flavor hadron decays at the LHC[1].
Given the ample cross-section[2], LHCb has collected the largest sample of
charmed and beauty hadrons to date in the world. Leveraging this huge dataset
for measurements of the parameters of the CKM triangle has been one of the
major programs of the LHCb detector, especially involving bottom quarks, but
the prospects for measurements with charm and strange quarks remain relatively
unexplored. In these proceedings, we summarize the possible prospects of such
measurements at LHCb. We begin with a short review of the standard ways to
measure the elements of the CKM matrix involving only the first two
generations. Second, we highlight some of the challenges and opportunities of
reconstructing semileptonic decays at LHCb and in a hadron collider
environment. Next, we offer a few sample measurements in the charm sector. We
also review the state of lepton non-universality (LNU) measurements in the
charm sector. Finally, we present preliminary sensitivities of such
measurements by the LHCb experiment.

##  2 The landscape of $V_{qq^{\prime}}$ measurements

Most measurements of the first two generations of the CKM matrix either rely
on either fully leptonic decay, semileptonic decay or other processes
involving, for instance, nuclear decays[3]. Precision measurements
contributing to the world average of $|V_{ud}|=0.97420\pm 0.00021$ currently
come from pure vector transitions of nuclei involving $\beta$ decay, resulting
in a relative precision of 0.2 per mille, with the limiting uncertainty coming
from theory. Measurements using the decay $\pi^{+}\to\pi^{0}e^{+}\nu_{e}$ are
30 times less precise, but still in agreement with nuclear $\beta$ decays.
Only the pion decay would be possible at LHCb, but these single prong decays
are not feasible in the LHCb environment, where charged pions are effectively
stable and single kinks are not well reconstructed in the tracking system.

Measurements of $|V_{us}|$ are dominated by semileptonic kaon decays, mainly
tapping $K_{L}^{0}$ and $K^{+}$ decays, and using the ratio with pion decays
to extract $|V_{us}|$[4, 3]. The most general way to report results of
semileptonic decays are then to list $|V_{us}|\times f_{+}(q^{2}=0)$, then use
either lattice information or otherwise to extract $|V_{us}|$ for this decay.
Additional information from
$\frac{|V_{us}|^{2}}{|V_{ud}|^{2}}\frac{f_{K}^{2}}{f_{\pi}^{2}}$ and inputting
the kaon and pion decay constants from leptonic $K$ and $\pi$ decays are not
as precise, but still contribute to the world average. The relative precision
for $|V_{us}|=0.2243\pm 0.0005$[3] now stands at 2 per mille.

In the LHCb Upgrade II, given the extremely high cross-section of strange
hadrons, measurements of branching fractions at the level of $10^{-15}$ are
possible, assuming perfect reconstruction efficiencies[5]. This opens the door
to the possible first observation of the $K_{s}^{0}\to\pi^{+}\mu^{-}\nu_{\mu}$
decay, and possible contribution to the measurement of $|V_{us}|f^{+}(q^{2})$.
Additional measurements are possible utilizing the decay $\Lambda^{0}\to
p\mu^{+}\nu_{\mu}$, but in both cases, lifetime of the $K_{S}^{0}$ and the
$\Lambda$ will have significant impact on both the ability to trigger these
candidates as well as their reconstruction, given that most will decay
downstream of the vertex locator, increasing the difficulty of vertexing.

Decays of charmed mesons allow for the extraction of $|V_{cd}|$. The usual
methods involve either the purely leptonic decay $D^{+}\to\ell^{+}\nu_{\ell}$,
or semileptonic decays of a singly charmed hadron $H_{c}\to
H_{d}\ell\nu_{\ell}$. Current world averages give $|V_{cd}|=0.218\pm 0.004$,
with a relative precision of about 2%[3]. Previous measurements include
neutrino scattering, a technique which is beyond the capabilities of the LHC.
Purely leptonic single prong $D^{+}$ decays are also not favored, as in the
case of the pion decay. However, the decays
$D^{0}\to\pi^{-}\ell^{+}\nu_{\ell}$ are accessible at LHCb, especially if
originating from the decay chain $D^{*}(2010)^{+}\to D^{0}\pi_{s}^{+}$. As
these are abundantly produced, the main uncertainties will be with the lepton
chosen, as well as the reconstruction method of the missing neutrino. We will
approach how to deal with these uncertainties in the coming sections. Further
decays involving $\Lambda_{c}$ baryons are possible, though much more
challenging, as the lifetime of the $\Lambda_{c}^{+}$ baryon is a factor two
less than that of the $D^{0}$, which limits prompt production trigger
efficiency. Tagging from secondary $b$-hadron decays is possible, but suffers
from lower statistics due to the difference in cross-sections between $c-$ and
$b-$hadrons.

The extraction of $|V_{cs}|$ is often done in analogy to that of $|V_{cd}|$,
either substituting a kaon for a pion, or $D_{s}^{+}$ for the $D^{+}$ decays.
As is the case with $|V_{cd}|$, world averages of $|V_{cs}|=0.997\pm 0.017$,
places the relative precision just below 2% [3], and as with the $D^{+}$
decays, single prong $D_{s}^{+}$ decays remain out of reach of LHCb
measurements. Measurements involving on-shell $W$ decays is out of scope of
LHCb analyses. Semileptonic decays, however, are still fair game, and the
Cabibbo favored $D^{0}\to K^{-}\ell^{+}\nu_{\ell}$ decays reconstructed at
LHCb will easily have the largest statistics of any sample in the world. The
measurements then rely on the explicit cancellation of detector and
theoretical uncertainties, which motivates the measurement of
$\frac{|V_{cd}|\times f_{+}^{D^{0}\to\pi}(q^{2})}{|V_{cs}|f_{+}^{D^{0}\to
K}(q^{2})}$ using
$R=\frac{\mathcal{B}(D^{0}\to\pi\mu\nu)}{\mathcal{B}(D^{0}\to K\mu\nu)}$.
Utilizing this ratio, and measuring as a function of $q^{2}$ leaves only the
efficiencies kaons and pions to be calculated, which are well known at LHCb[1,
6].

An interesting new approach is to extract of $|V_{cd}|$ and $|V_{cs}|$ using
the reconstruction of the decays $B_{c}^{+}\to B^{0}_{(s)}\mu^{+}\nu_{\mu}$
and $B_{c}^{+}\to B^{0}_{(s)}\pi^{+}$, tapping into the rich $B_{c}$ physics
program of LHCb[7]. We propose these measurements here for the first time. For
extraction of only $|V_{cs}|$, one can also look to measure $B_{c}^{+}\to
B^{0}K^{+}$. As with the semileptonic measurements of $D^{0}\to
K^{-}\ell^{+}\nu_{\ell}$ and $D^{0}\to\pi^{-}\ell^{+}\nu_{\ell}$, taking the
ratio of different channels allows for the effective cancellation of many
uncertainties, and leaves more theoretically clean observables. Some
limitation are high statistics samples of $B_{c}$ mesons, and the knowledge of
the $B_{c}$ hadronization fractions. Both of these points will be improved by
the LHCb experiment in years to come via direct measurements, but it may also
be possible to avoid some of these constraints by instead of focusing on the
CKM elements, taking ratios with the $D$ sector to isolate the form factors
and hadronization fractions.

For the remainder of these proceedings, we will focus on the measurements
involving semileptonic charm meson decays to extract quantities of interest.

##  3 Reconstruction of semileptonic decays at LHCb

The hadron collider environment of LHCb does not allow for the conventional
techniques of reconstruction employed at $e^{+}e^{-}$ colliders; one cannot
simply rely on the beam energies and full reconstruction of a tag side of the
decay to solve for the momentum of the missing neutrino. Partial
reconstruction is already known given a displaced secondary vertex, namely the
momentum of the neutrino perpendicular to this flight direction. The remaining
component has a two fold ambiguity with respect to this flight direction. One
can either choose to not reconstruct this component or reconstruct given extra
information. If one chooses not to reconstruct the component, then can use the
so-called “corrected mass” $m_{corr}=\sqrt{p_{T}^{\prime
2}+m_{vis}^{2}}+|p_{T}^{\prime}|$; using the example of the decay $D^{0}\to
K^{-}\mu^{+}\nu_{\mu}$, the variable $p_{T}^{\prime}$ is the momentum
component of the visible system perpendicular to the $D^{0}$ flight direction,
and $m_{vis}$ is the visible invariant mass of the fully reconstructed
candidates, in this case $m(K\mu)$. This variable provides distinguishing
power, as if there is truly one massless particle missing, then it will peak
at the correct mass, and have a tail to lower values, whereas if there are
more missing particles, the distribution will be shifted to lower values and
not peak as strongly. This is show in Figure 1, generated using RapidSim[8].
Other reconstruction methods include the well-known $k-$factor approach, where
the true momentum of the $D^{0}$ candidate is given with a function of the
visible mass and a factor $k$ to solve for the visible momentum,
$p(K\mu\nu)=\frac{p(K\mu)}{k(m_{vis})}$, multivariate regression[9], and using
higher resonance decays to break the two-fold ambiguity of momentum by adding
an additional mass constraint, for instance using the cone closure method[10].
One can also use multiple methods to exploit differing correlations between
signal and background across multiple variables.

Figure 1: Corrected mass distributions for $D^{0}\to K^{-}\mu^{+}\nu_{\mu}$
(green) and $D^{0}\to K^{-}\mu^{+}\nu_{\mu}\pi^{0}$ (red) decays. Uses
RapidSim[8].

##  4 Sensitivity estimates for
$\frac{|V_{cd}|f_{+}^{D^{0}\to\pi}(q^{2})}{|V_{cs}|f_{+}^{D^{0}\to K}(q^{2})}$
at LHCb

The differential decay rate of the decay $D^{0}\to P^{-}\ell^{+}\nu_{\ell}$ at
leading order can be written as[11]

$\begin{split}\frac{d\Gamma(D^{0}\to
P^{-}\ell^{+}\nu_{\ell})}{dq^{2}}=|V_{cQ}|^{2}&\frac{G_{F}^{2}}{24\pi^{3}}\frac{(q^{2}-m_{\ell}^{2})^{2}\sqrt{E_{P}^{2}-m_{P}^{2}}}{q^{4}m_{D^{0}}^{2}}\\\
&\times\left[\left(1+\frac{m_{\ell}^{2}}{2q^{2}}\right)m_{D^{0}}(E_{P}^{2}-m_{P}^{2})|f_{+}(q^{2})|^{2}+\frac{3m_{\ell}^{2}}{8q^{2}}(m_{D^{0}}^{2}-m_{P}^{2})^{2}|f_{0}(q^{2})|^{2}\right],\end{split}$
(1)

where $Q$ represents the outgoing quark from the weak vertex, the terms
$f_{+}(q^{2})$ and $f_{0}(q^{2})$ are the vector and scalar form factors as a
function of the momentum transfer squared to the lepton-neutrino pair $q^{2}$,
respectively, used to parameterize the hadronic current. As the scalar form
factor is suppressed by the factor $m_{\ell}^{2}/m_{D^{0}}^{2}$, it is
customary to omit this form factor, leading to the simplified form

$\frac{d\Gamma(D^{0}\to
P_{x}\ell\nu)}{dq^{2}}=\frac{G_{F}}{24\pi^{3}}\left|\vec{p}_{P}\right|^{3}\left|V_{cx}\right|^{2}|f_{+}(q^{2})|^{2}.$
(2)

Therefore, taking the ratio of the differential decay rates of the decays
$D^{0}\to\pi^{-}\mu^{+}\nu_{\mu}$ and $D^{0}\to K^{-}\mu^{+}\nu_{\mu}$ results
in a measurement of $\frac{|V_{cd}|\times
f_{+}^{D\to\pi}(q^{2})}{|V_{cs}|\times f_{+}^{D\to K}(q^{2})}$. By utilizing
the cone closure technique from $D^{*+}\to D^{0}\pi^{+}$ decays, one can
accurately reconstruct the $q^{2}$ distributions of the $D^{0}$ decay.
Preliminary yields of the $D^{0}\to K^{-}\mu^{+}\nu_{\mu}$ decays from this
channel give $q^{2}$ integrated yields of roughly 5 million signal candidates
for Run 1 of the LHC[12, 13], which would yield a relative uncertainty on the
measurement of 0.2%. If the error on the form factors can at some point be
neglected, this would result in a measurement of the CKM elements at the level
of the current world average[3].

##  5 Sensitivity of measurements of Charge Parity Violation and mixing from
$D^{0}\to K\mu\nu_{\mu}$

The decay $D^{0}\to K\mu\nu$ is interesting itself from a mixing standpoint as
well: the right-sign $D^{0}\to K^{-}\mu^{+}\nu_{\mu}$ is entirely dominated by
Cabibbo-favored decay, whereas the wrong-sign decay is only accessible from
mixing. This leaves the ratio of wrong-sign to right-sign decays as

$R(t)=\frac{\mathcal{P}(D^{0}\to
K^{+}\mu^{-}\overline{\nu}_{\mu})}{\mathcal{P}(D^{0}\to
K^{-}\mu^{+}\nu_{\mu})}\propto\frac{x^{2}+y^{2}}{4}\left(\frac{t}{\tau}\right)^{2},$
(3)

where $\mathcal{P}$ represents the probability of the transition,
$x=2(m_{1}-m_{2})/(\Gamma_{1}+\Gamma_{2})$, y =
$(\Gamma_{1}-\Gamma_{2})/(\Gamma_{1}+\Gamma_{2})$, and $t$ is the $D^{0}$
decay time. Here $m_{1,2}$ is the mass of the different $D$ mass eigenstates,
and $\Gamma_{1,2}$ is the width of the different eigenstates.
Integrating over decay-time, this ratio is $R_{M}=\frac{x^{2}+y^{2}}{2}$.
Using the same 5M as a sample as a baseline and the current world averages for
$x$ and $y$, the projected number of WS decays would be around 700,
corresponding to an uncertainty on $R_{M}$ of the order of 0.01%[14]. These
projections, however, do not include form-factor uncertainties, nor shape
uncertainties in the fit. Extraction of CPV related quantities would be
possible by splitting by the charge of the pion from $D^{*+}\to D^{0}\pi^{+}$
decays.

##  6 Lepton non-universality in charm

As the flavor anomalies in the $b$-hadron sector still persist, it becomes
increasingly important to search in different regions for similar effects. One
such place is the semileptonic decays $D^{0}\to P\ell\nu_{\ell}$. As is
evident from equation 1, taking the ratio between electron and muon modes will
lead to large cancellations in form factor uncertainties, leaving a
theoretically clean ratio. This represents a measurement of lepton non-
universality in charm decays (LNU). Measurements of these types have only now
started to be completed by the BESIII collaboration[15, 16], and LHCb
measurements are expected soon.

The state of global measurements is shown in Figure 2. As pointed out in the
previous version of this plot[13], all points lie to one side of this line.
Current progress from BESIII maintains this trend, but increases the
precision.

Figure 2: Summary of measurements of lepton non-universality in $D^{0}\to
P\ell\nu$ decays. Reported is the ratio of muon ($D^{0}\to
P^{-}\mu^{+}\nu_{\mu}$) to electron ($D^{0}\to P^{-}e^{+}\nu_{e}$) decays for
$P=K^{-},\pi^{-},K^{*}(892)^{-}$. In black are the ratio of muon to electron
branching fractions as reported in the PDG 2016 edition [17], in red are the
current BESIII measurements, and in blue the estimated reach of LHCb from the
Run 1 dataset, with the central value set to the PDG 2016 value. The bottom-
most BESIII measurement is also provided in Ref. [15]. The theoretical
prediction given by S. Fajfer is shown in the red dashed line with the one
sigma uncertainty listed as the green hatched area.

Given the dataset from LHCb, and incorporating differing efficiencies for
different triggering, the Run 1 reach is estimated to increase the precision
from the 2016 average by one order of magnitude on the $q^{2}$ integrated
measurement. In the case of $D^{0}\to K^{-}\ell^{+}\nu_{\ell}$ decays,
possible new physics effects are also probed in the $q^{2}$ dependence of the
LNU ratio $R^{\mu/e}$ which are not present in the measurement of the
branching ratio alone, due mainly to form factor uncertainties[18]. The BESIII
collaboration has recently measured the $q^{2}$ dependence of this ratio, but
lacks the precision to distinguish between standard model and new physics
effects. The LHCb collaboration will provide the measurement of $R^{\mu/e}$ vs
$q^{2}$ by utilizing the $D^{*+}\to D^{0}\pi^{+}$ decay, and using the cone
closure technique. A pseudo-experiment fit is shown for illustrative purposes
of one of the two fit methods considered in Figure 3. The method is a two-
dimensional fit to the corrected mass of the $D^{0}$ meson versus the visible
mass difference $\Delta m_{visible}=m(K\ell\pi_{s})-m(K\ell)$. The different
components of peaking backgrounds and signal are shown in different colors.

Figure 3: Unofficial fit to one pseudo-experiment of $D^{0}\to
K^{-}\mu^{+}\nu_{\mu}$ (top) and $D^{0}\to K^{-}e^{+}\nu_{e}$ (bottom) decays
at LHCb. The fits show are projections onto the corrected $D^{0}$ mass and the
visible delta mass, left and right, respectively. Individual fit components
are shown in different colors.

##  7 Conclusions

LHCb has collected the largest charm sample in the world, and despite
operating in a hadron collider environment, offers strong opportunities to
increase knowledge in CKM elements in the first two generations, provide
sensitivity to mixing and CPV in the $D^{0}$ system, and unique contributions
to LNU searches in the charm sector utilizing semileptonic decays. Many
sensitivities were presented for the Run 1 dataset, with Run 2 providing an
additional $\sim 6~{}$fb-1 at $\sqrt{s}=13$ TeV. Clearly the contribution of
the LHCb experiment in this field are both motivated and needed, and
contributions will come soon. Additionally, new ideas on how to access
quantities of interest, both in rare kaon decays and leveraging the $B_{c}$
physics program will contribute to the knowledge in the sector.

## References

  * [1] LHCb collaboration, A. A. Alves Jr. et al.,  JINST 3, S08005 (2008).
  * [2] LHCb, R. Aaij et al.,  JHEP 03, 159 (2016), 1510.01707,  [Erratum: JHEP05,074(2017)].
  * [3] Particle Data Group, M. Tanabashi et al.,  Phys. Rev. D98, 030001 (2018).
  * [4] FlaviaNet Working Group on Kaon Decays, M. Antonelli et al.,  Eur. Phys. J. C69, 399 (2010), 1005.2323.
  * [5] A. A. Alves Junior et al.,  (2018), 1808.03477.
  * [6] A. Powell et al.,  PoS ICHEP2010, 020 (2010),  LHCb-PROC-2011-008.
  * [7] L. E. The LHCb Collaboration,  CERN Report No. LHCb-PUB-2018-009. CERN-LHCb-PUB-2018-009, 2018 (unpublished).
  * [8] G. Cowan, D. Craik, and M. Needham,  Computer Physics Communications 214, 239 (2017).
  * [9] G. Ciezarek, A. Lupato, M. Rotondo, and M. Vesterinen,  JHEP 02, 021 (2017), 1611.08522.
  * [10] W. E. Johns,  Measurements of the semileptonic decay of the neutral charmed meson $D^{0}\to K^{-}\mu^{+}$ muon-neutrino,  PhD thesis, Colorado U., 1995.
  * [11] S. Aoki et al.,  Eur. Phys. J. C77, 112 (2017), 1607.00299.
  * [12] LHCb Collaboration, R. Aaij et al.,  Phys. Rev. Lett. 117, 061803 (2016).
  * [13] LHCb, A. C. S. Davis,  PoS CKM2016, 025 (2017), 1703.10695.
  * [14] D. Mitzel,  Study of semileptonic $D^{0}$ decays for a measurement of charm mixing at LHCb,  PhD thesis, Physikalisches Institut Heidelberg, 2014.
  * [15] BESIII, M. Ablikim et al.,  (2018), 1810.03127.
  * [16] BESIII, M. Ablikim et al.,  Phys. Rev. Lett. 121, 171803 (2018), 1802.05492.
  * [17] Particle Data Group, C. Patrignani et al.,  Chin. Phys. C40, 100001 (2016).
  * [18] S. Fajfer, I. Nisandzic, and U. Rojec,  Phys. Rev. D91, 094009 (2015), 1502.07488.