Multivariate skew generalized t distributions

Multivariate skew generalized t (SGT) distributions are special cases of multivariate skew scale mixture of normal (SSMN) distributions [28] (pages 102-103) which we first introduce. A random variable Z is said to follow a p−variate SSMN distribution with location μ, scale Ω, and shape λ, if it can be represented as [29] (page 20, Eq 3.12):

Lemma 1 (see S1 Appendix for a proof) The free parameters (μ, δ, Ω and ν) of a SSMN distribution with the representation in Eq (2) are identifiable if and only if i) the scale mixing distribution F_U(⋅|ν) is identifiable and ii) F_U(⋅|ν) does not satisfy $$ F_U\left(u\right|\nu )=H\left(\frac{u}{{\nu }_k}\right|{\nu }_{-k})$$ for any element ν_k of ν and any distribution function H(⋅|ν_−k) where ν_−k is the vector ν without the element ν_k. If U has a probability density function (pdf) f_U(u|ν) for all u > 0, then the condition ii) is equivalent to f_U(⋅|ν) does not satisfy $$ f_U\left(u\right|\nu )=\frac{1}{{\nu }_k}{h}_U\left(\frac{u}{{\nu }_k}\right|{\nu }_{-k})$$ for any pdf h_U(⋅|ν_−k).

On setting $$ c=\sqrt{2/\pi }$$ and defining the expectations $$ {\tilde{U}}_t={\text{E}}_U\left\{U^{-t/2}\right\}$$, and assuming that $$ {\tilde{U}}_t<\infty $$ for the required expectations, the first two central moments of a SSMN vector Z are given by [28] (pages 109-110):

We notice from the expressions for skewness Eq (5) and kutosis Eq (6) indices that the parameter λ controls the shape of the distribution only through the working shape parameter δ = (1 + λ^⊤ λ)^−1/2 Ω^1/2 λ. This quantity is invariant under marginalization, i.e. by the stochastic representation in Eq (2), for any arbitrary subset of Z, the working shape parameter is the corresponding subset of δ. It is worth noticing however that the quantity δ cannot be specified unrestrictedly, independently of Ω. Indeed, we observe that δ = (1 + λ^⊤ λ)^−1/2 Ω^1/2 λ implies that λ = (1 + λ^⊤ λ)^1/2 Ω^−1/2 δ. This in turn gives λ^⊤ λ = (1 + λ^⊤ λ)δ^⊤ Ω⁻¹ δ which yields λ^⊤ λ(1−δ^⊤ Ω⁻¹ δ) = δ^⊤ Ω⁻¹ δ. We then get $$ {\mathbf{\lambda }}^{\top }\mathbf{\lambda }=\frac{{\mathit{\delta }}^{\top }{\mathbf{\Omega }}^{-1}\mathit{\delta }}{1-{\mathit{\delta }}^{\top }{\mathbf{\Omega }}^{-1}\mathit{\delta }}$$ so that $$ 1+{\mathbf{\lambda }}^{\top }\mathbf{\lambda }=\frac{1}{(1-{\mathit{\delta }}^{\top }{\mathbf{\Omega }}^{-1}\mathit{\delta })}$$, i.e. 1 + λ^⊤ λ = (1 − δ^⊤ Ω⁻¹ δ)⁻¹ provided δ^⊤ Ω⁻¹ δ ≠ 1. Therefore, λ is recorvered from δ and Ω under the constraint δ^⊤ Ω⁻¹ δ < 1 as:

Lemma 2 (see S2 Appendix for a proof) Let us consider the SGT distribution $$ {\mathcal{\text{SGT}}}_p(\mathit{\mu },\mathbf{\Omega },\mathbf{\lambda },\mathit{\nu })$$ with ν = (ν, ν₀)^⊤ and pdf in Eq (10). Set $$ {\mathbf{\Omega }}^{*}=\frac{{\nu }_0}{\nu }\mathbf{\Omega }$$. The following statements hold:

SGt_p(z|μ, Ω, λ, ν) = St_p(z|μ, Ω*, λ, ν);

If $$ \mathbf{Z}\sim {\mathcal{\text{SGT}}}_p(\mathit{\mu },\mathbf{\Omega },\mathbf{\lambda },\nu )$$ then $$ \mathbf{Z}\sim {\mathcal{\text{ST}}}_p(\mathit{\mu },{\mathbf{\Omega }}^{*},\mathbf{\lambda },\nu )$$.

Lemma 2 indicates that any SGT vector is a rescaled version of a ST vector. However, in the frame of binary data models, the use of a SGT distribution instead of a simple ST distribution as link function allows to control the scale of the link function through the parameter ν₀ [9]. Specifically, a skew generalized t-link allows to define a skewed and heavy-tailled binary mixed model where fixed effects have the same scale and interpretation as in the mixed probit model in Eq (1). Interestingly, the popular logit and probit links for binary data can be recast as special cases of the cdf of the SGT class of distributions. Indeed, the normal distribution is a limiting case of SGT distributions when ν₀ = ν → ∞ and λ = 0. Moreover, the logistic distribution is well appoximated by a rescaled Student t distribution with appropriate degrees of freedom [8] (page 228). These constatations make the SGT distributions appropriate for extending the traditional probit and logistic GLMMs, accounting for skewness and heavy tail behaviors. To this end, we present in the next section some results on truncated multivariate SGT distributions since binary data can reflect truncation of latent continuous variables.

Truncated multivariate skew generalized t distributions

As seen for the mixed probit model in Eq (1), models for binary data can be defined by truncating latent variables following continuous distributions. We define in this section a class of truncated multivariate skew generalized t distributions which are useful for a latent variable representation of skew generalized t-link binary data models. We also give expressions to evaluate some joint moments of a truncated multivariate skew generalized t distribution and a gamma distribution, as they prove useful in the implementation of the EM algorithm for the skew generalized t-link model.

Let $$ {\mathcal{\text{TSGT}}}_p(\mathit{\mu },\mathbf{\Omega },\mathbf{\lambda },\nu ,\mathbb{A})$$ represent a p−variate skew generalized t (SGT) vector restricted to a p-dimensional hyperplane $$ \mathbb{A}$$; with μ a p−vector (location), Ω a p × p positive definite matrix (scale), λ a p−vector (shape) and ν = (ν, ν₀)^⊤ a vector of positive scalars (degrees of freedom). The pdf of $$ \mathbf{Z}\sim {\mathcal{\text{TSGT}}}_p(\mathit{\mu },\mathbf{\Omega },\mathbf{\lambda },\mathit{\nu },\mathbb{A})$$ is:

In the frame of correlated binary data models, the truncation region $$ \mathbb{A}$$ typically has the form $$ \mathbb{A}={\mathbb{A}}_1\times {\mathbb{A}}_2\times \dots \times {\mathbb{A}}_p$$ where $$ {\mathbb{A}}_k$$ are real intervals of the form $$ {\mathbb{A}}_k=(-\infty ,a_k]$$ or $$ {\mathbb{A}}_k=(a_k,\infty )$$, for $$ a_k\in \mathbb{R}$$ (k = 1, 2, ⋯, p). Let us consider for instance a vector Y of three binary outcomes obtained by truncating the elements of a 3−variate SGT vector $$ \mathbf{Z}\sim {\mathcal{\text{SGT}}}_3(\mathit{\mu },\mathbf{\Omega },\mathbf{\lambda },\nu )$$: Y_k = 0 if Z_k ≤ 0 and Y_k = 1 if Z_k > 0. In practice, however, only the binary outcomes (Y) are observable whereas the latent outcome Z is unobservable. Suppose one observes the binary outcomes y = (1, 0, 1)^⊤. This implies that the corresponding value z of the latent vector Z satisfies z₁ > 0, z₂ ≤ 0, and z₃ > 0. The conditional distribution of Z given Y = y (required for maximum likelihood estimation using the EM algorithm) is thus $$ {\mathcal{\text{SGT}}}_3(\mathit{\mu },\mathbf{\Omega },\mathbf{\lambda },\mathit{\nu })$$ truncated to the region $$ {\mathbb{A}}_{e.g.}=(0,\infty )\times (-\infty ,0]\times (0,\infty )$$, i.e. $$ \mathbf{Z}|\mathbf{Y}=\mathbf{y}\sim {\mathcal{\text{TSGT}}}_p(\mathit{\mu },\mathbf{\Omega },\mathbf{\lambda },\mathit{\nu },{\mathbb{A}}_{e.g.})$$ as defined in Eq (13).

We shall use the simplified notation $$ {\mathcal{\text{TSGT}}}_p(\mathit{\mu },\mathbf{\Omega },\mathbf{\lambda },\mathit{\nu },\mathbf{a})$$ with $$ \mathbf{a}\in {\mathbb{R}}^p$$ to denote a truncated SGT distribution $$ {\mathcal{\text{TSGT}}}_p(\mathit{\mu },\mathbf{\Omega },\mathbf{\lambda },\nu ,\mathbb{A})$$ when $$ \mathbb{A}$$ is the right truncated hyperplane $$ \mathbb{A}=\{\mathbf{z}\in {\mathbb{R}}^p|z_1\le a_1,\dots ,z_p\le a_p\}$$. In this case, α_st = SGT_p(a|μ, Ω, λ, ν) with SGT_p(⋅|μ, Ω, λ, ν) the cdf of the p−variate ST distribution. This corresponds for instance to the situation where all binary outcomes are zeros. When λ = 0, the right truncated SGT distribution $$ {\mathcal{\text{TSGT}}}_p(\mathit{\mu },\mathbf{\Omega },0,\mathit{\nu },\mathbf{a})$$ is a right truncated GT distribution denoted $$ {\mathcal{\text{TGT}}}_p(\mathit{\mu },\mathbf{\Omega },\mathit{\nu },\mathbf{a})$$ whose pdf and cdf are respectively denoted TGt_p(⋅|μ, Ω, ν, a) and TGT_p(⋅|μ, Ω, ν, a). When ν₀ = ν, the right truncated SGT distribution is a right truncated ST distribution denoted $$ {\mathcal{\text{TST}}}_p(\mathit{\mu },\mathbf{\Omega },\mathbf{\lambda },\mathit{\nu },\mathbf{a})$$ with pdf TSt_p(⋅|μ, Ω, λ, ν, a) and cdf TST_p(⋅|μ, Ω, λ, ν, a). If both λ = 0 and ν₀ = ν, the distribution is reduced to a right truncated t distribution denoted $$ {\mathcal{\text{TT}}}_p(\mathit{\mu },\mathbf{\Omega },\nu ,a)$$ with pdf Tt_p(⋅|μ, Ω, ν, a) and cdf TT_p(⋅|μ, Ω, ν, a). In the above example, if y = (0, 0, 0)^⊤, then the truncation region becomes $$ {\mathbb{A}}_{e.g.}=(-\infty ,0]\times (-\infty ,0]\times (-\infty ,0]$$. Since all truncation points are zeros, we shall write in this case $$ \mathbf{Z}|\mathbf{Y}=\mathbf{y}\sim {\mathcal{\text{TSGT}}}_p(\mathit{\mu },\mathbf{\Omega },\mathbf{\lambda },\mathit{\nu },\mathbf{a})$$ with a = (0, 0, 0)^⊤ using the above simplified notation.

The implementation of an EM algorithm for a SGT distribution based binary data model requires joint moments of the form $$ \overline{u_r}=\text{E}\left\{U^{r/2}\right\}$$, $$ \overline{u_r{\mathbf{z}}_s}=\text{E}\left\{U^{r/2}{\mathbf{Z}}^{\left(s\right)}\right\}$$, $$ \overline{{\tau }_r}=\text{E}\left\{U^{r/2}{\zeta }_1\left(U^{1/2}\alpha \right)\right\}$$ and $$ \overline{{\tau }_r{\mathbf{z}}_s}=\text{E}\left\{U^{r/2}{\zeta }_1\left(U^{1/2}\alpha \right){\mathbf{Z}}^{\left(s\right)}\right\}$$ for s ∈ {1, 2}, Z⁽¹⁾ = Z and Z⁽²⁾ = Z Z^⊤, $$ U\sim \mathcal{G}amma(\nu /2,{\nu }_0/2)$$, $$ \mathbf{Z}\sim {\mathcal{\text{TSGT}}}_p(\mathit{\mu },\mathbf{\Omega },\mathbf{\lambda },\nu ,\mathbb{A})$$, α = λ^⊤ Ω^−1/2(Z − μ), ζ₁(x) = ϕ(x)/Φ(x) with ϕ(⋅) the pdf of the standard normal distribution, and $$ \mathbb{A}$$ is an hyperplane of the form $$ \mathbb{A}={\mathbb{A}}_1\times \dots \times {\mathbb{A}}_p$$ with $$ {\mathbb{A}}_k\in \{(-\infty ,a_k],(a_k,\infty \left)\right\}$$, (a_k, ∞)} for a = (a₁, ⋯, a_p)^⊤. Proposition 1 hereafter will be useful for the derivation of $$ \overline{u_r}$$, $$ \overline{u_r{\mathbf{z}}_s}$$, $$ \overline{{\tau }_r}$$ and $$ \overline{{\tau }_r{\mathbf{z}}_s}$$.

Proposition 1 (see S3 Appendix for a proof) Let $$ \mathbf{Z}\sim {\mathcal{\text{TSGT}}}_p(\mathit{\mu },\mathbf{\Omega },\mathbf{\lambda },\mathit{\nu },\mathbb{A})$$ with ν = (ν, ν₀)^⊤, $$ U\sim \mathcal{G}amma(\nu /2,{\nu }_0/2)$$ and set α = λ^⊤ Ω^−1/2(Z − μ). Then, for any real r > − ν and an integrable function g(⋅) of Z:

By the mean of a simple linear transformation, we obtain from Proposition 1 the joint expectations $$ \overline{u_r}$$, $$ \overline{u_r{\mathbf{z}}_s}$$, $$ \overline{{\tau }_r}$$ and $$ \overline{{\tau }_r{\mathbf{z}}_s}$$ in terms of moments of a truncated multivariate skew t distribution.

Corollary 1 (see S4 Appendix for a proof) Let $$ \mathbf{Z}\sim {\mathcal{\text{TSGT}}}_p(\mathit{\mu },\mathbf{\Omega },\mathbf{\lambda },\mathit{\nu },\mathbb{A})$$ with ν = (ν, ν₀)^⊤, $$ U\sim \mathcal{G}amma(\nu /2,{\nu }_0/2)$$, $$ \mathbf{a}\in {\mathbb{R}}^p$$ and $$ \mathbb{A}={\mathbb{A}}_1\times \dots \times {\mathbb{A}}_p$$ with $$ {\mathbb{A}}_k=(-\infty ,a_k]$$ or $$ {\mathbb{A}}_k=(a_k,\infty )$$. Then, on setting $$ \mathit{\delta }={\left(\frac{{\nu }_0}{\nu }\right)}^{1/2}{(1+{\mathbf{\lambda }}^{\top }\mathbf{\lambda })}^{-1/2}{\mathbf{\Omega }}^{1/2}\mathbf{\lambda }$$, $$ \overline{\mathbf{\Omega }}=\frac{{\nu }_0}{\nu }\mathbf{\Omega }-\mathit{\delta }{\mathit{\delta }}^{\top }$$, A = diag(A₁, ⋯, A_p) with A_k = 1 if $$ {\mathbb{A}}_k=(-\infty ,a_k]$$ and A_k = −1 if $$ {\mathbb{A}}_k=(a_k,\infty )$$, a* = Aa, μ* = A μ, $$ {\mathbf{\Omega }}^{*}=\frac{{\nu }_0}{\nu }\mathbf{A}\mathbf{\Omega }\mathbf{A}$$, $$ {\overline{\mathbf{\Omega }}}^{*}=\mathbf{A}\overline{\mathbf{\Omega }}\mathbf{A}$$, λ* = Aλ and α_st = ST_p(a*|μ*, Ω*, λ*, ν),

For a practical use of Corollary 1, the cumulative multivariate skew t distribution is required. To this end, the function pmst of the package sn [35] in R freeware [36] is an appropriate routine.

Moments of truncated multivariate skew generalized t distributions

The evaluation of expectations $$ \text{E}\left\{{\mathbf{X}}_{u,r}^{\left(s\right)}\right\}$$ involved in Corollary 1 calls for general expressions for the first and second order moments of truncated multivariate SGT distributions. These moments are required in the implementation of an EM algorithm for a SGT distribution based binary data model. The moments have been derived for truncated multivariate t distributions by [34] and were used by [24] in their implementation of the EM algorithm for a t-link GLMM. We present in this section the expressions for the first two moments of the multivariate SGT distributions, relying on the Theorem 1 of [37] and the moments of truncated multivariate t distributions available from [34] (see also [38]).

Let $$ \mathbf{Z}\sim {\mathcal{\text{TSGT}}}_p(\mathit{\mu },\mathbf{\Omega },\mathbf{\lambda },\mathit{\nu },\mathbf{a})$$ with ν = (ν, ν₀)^⊤ and $$ \mathbf{a}\in {\mathbb{R}}^p$$, i.e. a p−variate SGT vector restricted to the right truncated hyperplane $$ \mathbb{A}=\{\mathbf{x}\in {\mathbb{R}}^p|x_1\le a_1,\dots ,x_p\le a_p\}$$. The pdf of Z is:

Proposition 2 (see S5 Appendix for a proof) Let $$ \mathbf{Z}\sim {\mathcal{\text{TST}}}_p(0,\mathbf{R},\mathbf{\lambda },\mathit{\nu },\mathbf{a})$$ with R a correlation matrix. Then,

The following corollary gives the first two moments of a general right truncated SGT vector $$ \mathbf{Z}\sim {\mathcal{\text{TSGT}}}_p(\mathit{\mu },\mathbf{\Omega },\mathbf{\lambda },\mathit{\nu },\mathbf{a})$$ with ν = (ν, ν₀)^⊤.

Corollary 2 Let $$ \mathbf{Z}\sim {\mathcal{\text{TSGT}}}_p(\mathit{\mu },\mathbf{\Omega },\mathbf{\lambda },\mathit{\nu },\mathbf{a})$$ with ν = (ν, ν₀)^⊤. Then,

When ν → ∞, the truncated multivariate SGT family has the truncated multivariate skew normal family as a limiting case (see S5 Appendix for a definition and formulas for computing the first two moments).

Model specification and marginal log-likelihood

The skew generalized t-link GLMM (SGTLM) considered in this work is defined as:

In the SGTLM, the distribution of a single latent outcome Z_ij is $$ \mathcal{\text{SGT}}({\mathit{\eta }}_{ij}-c{\tilde{U}}_1{\upsilon }_0{\mathit{\delta }}_{\epsilon },{\omega }_{\epsilon }^2,{\mathit{\delta }}_{\epsilon },\nu )$$ where $$ {\omega }_{\epsilon }^2=1+{\mathit{\delta }}_{\epsilon }^2$$ and $$ \mathcal{\text{SGT}}(\mathit{\mu },{\omega }^2,\mathbf{\lambda },\nu )$$ denotes a univariate SGT distribution with location μ, scale ω², shape λ and degrees of freedom ν. Therefore, on denoting SGT(⋅|μ, ω², δ, ν) the cdf of a scalar SGT distribution $$ \mathcal{\text{SGT}}(\mathit{\mu },{\omega }^2,\mathit{\delta },\mathit{\nu })$$, the success probability of an outcome Y_ij given the random effects b_i is $$ P(Y_{ij}=1|{\mathbf{b}}_i)=SGT\left(0\right|{\mathit{\eta }}_{ij}-c{\tilde{U}}_1{\upsilon }_0{\mathit{\delta }}_{\epsilon },{\omega }^2,{\mathit{\delta }}_{\epsilon },\mathit{\nu })$$. Unlike in the common probit model (PM) (see Eq (1)), for a given cluster i, the n_i latent outcomes Z_ij are not independent given the random effects b_i. Indeed, on using Eq (4) and setting $$ {\tilde{U}}_2=\frac{\nu }{\nu -2}$$, the variance-covariance matrix of Z_i given b_i is

The positive constant υ₀ in the SGTLM controls the scale of the latent variable Z_ij and thus the scale of the model link function. Indeed, from Eq (28), the conditional variance of Z_ij is $$ \text{V}ar\left\{Z_{ij}\right|{\mathbf{b}}_i\}={\upsilon }_0^2[{\tilde{U}}_2+({\tilde{U}}_2-c^2{\tilde{U}}_1^2\left){\mathit{\delta }}_{\epsilon }^2\right]$$. Setting υ₀ = 1 would yield a skew t-link model (i.e. ν = (ν, ν)^⊤). However, to make fixed effects in the SGTLM comparable with fixed effects in the common probit model (PM) characterized by a link function with a unit scale (i.e. Var{Z_ij|b_i} = 1), we have set

By Eq (2), the SGTLM has the stochastic representation

Proposition 3 (see S6 Appendix for a proof). Let us consider the latent vector Z_i and the binary variable Y_ij in Eq (27) and define $$ {\mathit{\mu }}_i={\mathbf{X}}_i\beta -c{\tilde{U}}_1{\mathbf{\Delta }}_i$$ with $$ {\mathbf{\Omega }}_i={\overline{\mathbf{\Omega }}}_i+{\mathbf{\Delta }}_i{\mathbf{\Delta }}_i^{\top }$$, $$ {\overline{\mathbf{\Omega }}}_i={\upsilon }_0^2{\mathbf{I}}_{n_i}+{\mathbf{W}}_i\overline{\mathbf{D}}{\mathbf{W}}_i^{\top }$$, $$ {\mathbf{\Delta }}_i={\upsilon }_0{\delta }_{\epsilon }{\mathbf{J}}_{n_i}+{\mathbf{W}}_i\mathit{\delta }$$, and the related shape parameter $$ {\mathbf{\lambda }}_i={\mathbf{\Omega }}_i^{-1/2}{\mathbf{\Delta }}_i{(1+{\mathbf{\Delta }}_i^{\top }{\overline{\mathbf{\Omega }}}_i^{-1}{\mathbf{\Delta }}_i)}^{1/2}$$. Then $$ {\mathbf{Z}}_i\sim {\mathcal{\text{ST}}}_{n_i}({\mathit{\mu }}_i,{\mathbf{\Omega }}_i,{\mathbf{\lambda }}_i,\nu )$$. Furthermore, the vector of binary outcomes Y_i has a multivariate Bernoulli distribution with joint probability mass function,

Eq (32) conveniently expresses for a value y_i of Y_i, the likelihood as a cumulative probability of a ST distribution whose location, scale and shape parameters depend on y_i, using the identities P(Z_ij > z_ij) = P(−Z_ij < −z_ij) and sign(Z_ij) = 2Y_ij − 1 where sign(⋅) returns the sign of its argument. On using Eq (4) on the distribution of Z_i given in Proposition 3, the variance-covariance matrix of the outcomes at the latent scale is

The parameters of the SGTLM include β, δ_ε, δ, $$ vech\left(\overline{\mathbf{D}}\right)$$ and ν where the vech(⋅) operator returns the lower triangle elements of its matrix argument. In order to avoid non-regular likelihood problems occurring in models based on the Student distribution and its extensions (in particular when ν is close to zero) [40], we follow some recent related works [24, 41] and first consider ν as known, focusing on $$ \mathit{\theta }={({\mathit{\beta }}^{\top },{\delta }_{\epsilon },{\mathit{\delta }}^{\top },{vech\left(\overline{\mathbf{D}}\right)}^{\top })}^{\top }$$. Classical inference on θ is based on the marginal likelihood of the observed data y. Using Proposition 3, the joint marginal log-likelihood of n independent clusters y_i (i = 1, 2, ⋯, n) is:

Model identifiability

Estimations in the skew generalized t-link model (SGTLM) may produce inconsistent results which would induce unreliable and misleading conclusions, if the model is not identifiable. It is thus of great importance to check whether different points in the parameter space can be distinguished from observations y_i. We analyse in this section the identifiability of the SGTLM and indicate when it is necessary to restrict the parameter space to ensure reliable inference from observed data. We restrict attention to the case υ₀ = 1 (ignoring Eq (29)) since υ₀ is an artificial device only included to ensure a unit variance in the conditional link function as in the traditional probit mixed model.

Although not sufficient, the identifiability of the SGTLM requires both the marginal random effects distribution and the conditional model given random effects to be identifiable. The identifiability of the random effects distribution follows from the identifiability of multivariate skew t distributions. We survey the identifiability of the conditional model before turning to the marginal model.

• Conditional identifiability

Conditional on the random effects b_i and for fixed degrees of freedom parameter ν, the identifiability of the SGTLM reduces to the identifiability of the fixed effects skew-probit model. This follows because the skew t inverse link function is an average of the skew-probit inverse link function with respect to the gamma mixing distribution. The identifiability of the skew-probit model with one covariate has been recently shown to depend on the nature (binary/continuous) of the covariate in the model [42]. Indeed, the fixed effect skew-probit model is not identifiable in the absence of any covariate (i.e. each X_i is a column of n_i ones) [42] (page 1624, Proposition 2.1) or in the presence of a binary covariate [42] (page 1626, Proposition 2.2). On the other hand, the fixed effect skew-probit model is identifiable when the covariate is continuous [42] (page 1627, Proposition 2.3). Extension to the case of multiple covariates is straightforwardly obtained by requiring the covariate matrix $$ \mathbf{X}={({\mathbf{X}}_1^{\top },{\mathbf{X}}_2^{\top },\dots ,{\mathbf{X}}_n^{\top })}^{\top }$$ to be of full column rank as in the classical linear regression model context. Whenever binary covariates are considered or no covariate is considered, we advocate to set δ_ϵ = 0 so that the conditional model reduces to a classical probit model.

• Marginal identifiability

From Proposition 3, it appears that when υ₀ = 1 the paramaters δ_ε and δ enter the marginal distribution of Y_i only through the marginal working shape $$ {\mathbf{\Delta }}_i={\delta }_{\epsilon }{\mathbf{J}}_{n_i}+{\mathbf{W}}_i\mathit{\delta }$$ whose j^th element can be written $$ {\mathbf{\Delta }}_{ij}={\delta }_{ϵ}+{\mathbf{W}}_{ij}^{\top }\mathit{\delta }$$. As a result, caution is required when learning the model parameters from some realizations y_i of Y_i. Indeed, if the model includes a random intercept term, i.e.W_ij has the form $$ {\mathbf{W}}_{ij}={(1,{\mathbf{W}}_{1ij}^{\top })}^{\top }$$ where $$ {\mathbf{W}}_{1ij}\in {\mathbb{R}}^{q-1}$$, we can partition the random effects working shape as $$ \mathit{\delta }={({\delta }_0,{\mathit{\delta }}_1^{\top })}^{\top }$$ with $$ {\mathit{\delta }}_1\in {\mathbb{R}}^{q-1}$$ so that the j^th element of the marginal working shape reads $$ {\mathbf{\Delta }}_{ij}={\delta }_{ϵ}+{\delta }_0+{\mathbf{W}}_{1ij}^{\top }{\mathit{\delta }}_1$$. Therefore, only the sum δ_ϵ + δ₀ could be estimated and it would not be possible to distinguish in the observed skewness the part due to the random intercept from the part due to the conditional link function. This confounding issues may actually be avoided by considering more complex models based for instance on the fundamental skew distributions [43]. Recall that skewness is introduced in the SGTLM through a hidden standard half normal variable, namely V_i. As opposed to the unique standard half normal variable V_i used for both Z_i and b_i in the SGTLM, the use of two different standard half normal hidden variables for Z_i and b_i [44] (page 667 eq 5-6) or two different standard half multivariate normal hidden vectors [45] (page 420 eq 2.2) remove the confounding problem.

Fortunately, the non identifiability of the skewness of conditional link function and random intercept does not affect the success probability of the response since this only depends on Δ_i, but not on the individual values of δ_ϵ and δ₀. However, since the conditional link scale depends on δ_ϵ through υ₀, the confounding problem affects the scale and thus the interpretation of fixed effects. Moreover, inference on the random intercept is affected since the random intercept variance and skewness depend on δ₀. For example, δ_ϵ + δ₀ = 0 only indicates null marginal skewness, and in no way absence of link function and random intercept skewness which could be equally strong but of opposite signs. Thus, only a lower bound can be given to the random intercept variance: $$ {\tilde{U}}_2{\overline{D}}_{11}+({\tilde{U}}_2-c^2{\tilde{U}}_1^2){\mathit{\delta }}_0^2\ge {\tilde{U}}_2{\overline{D}}_{11}$$ where $$ {\overline{D}}_{11}$$ is the first diagonal element of $$ \overline{\mathbf{D}}$$. To rule out this peculiar situation where the model is not marginally fully identifiable, some previous works on skew normal/skew t distributions have considered the restriction δ_ϵ = 0 (regardless of the presence of a random intercep term) for instance in the context of linear mixed effects models [30, 46, 47] (page 1492 eq 2, page 4100 eq 4, and page 309 eq 3.2 respectively), multivariate measurement error models [48] (page 35, Eq 4.11) and non linear mixed effect models [49] (page 7 eq 10), but no argument was given to support this choice.

In the very common situation where the mixed model includes a random intercept term, prior information on the shape of the link function and/or the random intercept is required to place a meaningful restriction on the parameter space by setting for instance δ_ϵ = 0 or δ₀ = 0 or δ_ϵ = δ₀. In the absence of such information, we advocate to consider the restriction δ₀ = 0 because the success probability of a response may exibit skewness, irrespective of the presence of random effects. This restriction thus allows to recover a fixed effects skew generalized t-link model when no random effect is considered. Overall, the restriction δ₀ = 0 simply expresses the unability of the SGTLM to capture any additional skewness structure from the data through the inclusion of only random intercept. For completeness, we develop in the next section an estimation procedure for the full model in Eq (27), since the introduction of any equality restriction on δ_ϵ and δ₀ can be straightforwardly reduced to δ₀ = 0.

The two restrictions discussed in this section are related to the structure of the quantity Δ_i and are required only for some specific data structures (models including a binary covariate or models with a random intercept only). However, even when a restriction is required on δ_ε or the first element of δ, the quantity Δ_i itself can remain unrestricted. When a restriction is required, it forces the skewness from a data to be summarized either by δ_ε or elements of δ. Overall, the SGT-link function is always allowed to be skewed (unconditional link). But some designs do not allow to distinguish skewness in the conditional link function from skewness in the distribution of random effects.

Maximum likelihood inference

Estimation via the EM algorithm

The choice of the value of υ₀ in Eq (29) is in line with one of our purposes: rescale fixed effects so that they have the same interpretation as in the mixed probit model. There is no need to define υ₀ depending on situations, because in our proposal, υ₀ is fixed. However, since υ₀ is simply a scaling factor, it may be given any positive value during estimation, as long as the estimates are rescaled after the convergence of the estimation procedure so that υ₀ is finally given by Eq (29). Indeed, because routines are basically written for the skew t distributions, we used υ₀ = 1 in the EM algorithm and rescaled the estimates at the end of the procedure. Let us consider the complete data $$ {\mathbf{y}}_{co{m}_i}=\{{\mathbf{y}}_i,{\mathbf{z}}_i,{\mathbf{b}}_i,u_i,v_i\}$$. Because y_i only retains the signs of elements of z_i, the joint density of y_i and z_i is $$ f_{\mathbf{Y}\mathbf{Z}}({\mathbf{y}}_i,{\mathbf{z}}_i)=f_{\mathbf{Z}}\left({\mathbf{z}}_i\right)\times I_{{\mathbb{A}}_i}\left({\mathbf{z}}_i\right)$$ where $$ {\mathbb{A}}_i={\mathbb{A}}_{i1}\times \dots \times {\mathbb{A}}_{i{n}_i}$$ with $$ {\mathbb{A}}_{ij}=(-\infty ,0]$$ if y_ij = 0 and $$ {\mathbb{A}}_{ij}=(0,\infty )$$ if y_ij = 1. The density of $$ {\mathbf{y}}_{co{m}_i}$$ is thus $$ f_{{\mathbf{y}}_{com}}\left({\mathbf{y}}_{co{m}_i}\right)=f_{\mathbf{Z},\mathbf{b},U,V}({\mathbf{z}}_i,{\mathbf{b}}_i,u_i,v_i)\times I_{{\mathbb{A}}_i}\left({\mathbf{z}}_i\right)$$. Hence by Bayes’s rule and in light of Eq (27) with υ₀ = 1, the density of $$ {\mathbf{y}}_{co{m}_i}$$ is:

where f_V(v_i) = 2ϕ(v_i) I_(0,∞)(v_i) and f_G(⋅|ν/2, ν/2) is given in Eq (9). By Eq (37) and on setting $$ N={\sum }_{i=1}^n{n}_i$$ and $$ C_{\ell }=-\frac{N+n(q+1)}{2}\text{log}\left(2\pi \right)+n\text{log}2$$, the complete data log-likelihood is:

where $$ {\gamma }_i={\mathbf{z}}_i-{\mathbf{W}}_i{\mathbf{b}}_i-{\mathit{\delta }}_{\epsilon }(v_i{u}_i^{-1/2}-c{\tilde{U}}_1){\mathbf{J}}_{n_i}$$, and tr(⋅) is the trace operator. Let $$ {\widehat{\mathit{\theta }}}^{\left(k\right)}$$ the estimate of θ at the k^th EM iteration. The E-step of the (k + 1)^th iteration finds the expectation $$ Q(·|{\widehat{\mathit{\theta }}}^{\left(k\right)})$$ of ℓ(⋅|y_com) given the observed data y and the current parameter estimate $$ {\widehat{\mathit{\theta }}}^{\left(k\right)}={({\widehat{\mathit{\beta }}}^{\left(k\right)\top },{\widehat{\delta }}_{\epsilon }^{\left(k\right)},{\widehat{\mathit{\delta }}}^{\left(k\right)\top },{vech\left({\widehat{\overline{\mathbf{D}}}}^{\left(k\right)}\right)}^{\top })}^{\top }$$:

where $$ {\widehat{u_2\gamma }}_i^{\left(k\right)}={\widehat{u_2\mathbf{z}}}_i^{\left(k\right)}-{\mathbf{W}}_i{\widehat{u_2\mathbf{b}}}_i^{\left(k\right)}-{\delta }_{\epsilon }({\widehat{vu}}_i^{\left(k\right)}-c{\tilde{U}}_1{\widehat{u_2}}_i^{\left(k\right)}){\mathbf{J}}_{n_i}$$ and

so that we have

Note that the conditional expectation of $$ \text{log}{I}_{{\mathbb{A}}_i}\left({\mathbf{z}}_i\right)$$ is 0 since given y_i, $$ I_{{\mathbb{A}}_i}\left({\mathbf{z}}_i\right|{\mathbf{y}}_i)=1$$. The E-step thus reduces to the computation of the conditional expectations $$ {\widehat{u_2\mathbf{b}}}_i^{\left(k\right)}=\text{E}\{U_i{\mathbf{b}}_i,{\widehat{\mathit{\theta }}}^{\left(k\right)}\}$$, $$ {\widehat{vu\mathbf{b}}}_i^{\left(k\right)}=\text{E}\left\{V_i{U}_i^{1/2}{\mathbf{b}}_i\right|{\mathbf{y}}_i,{\widehat{\mathit{\theta }}}^{\left(k\right)}\}$$, $$ {\widehat{u_2\mathbf{b}\mathbf{z}}}_i^{\left(k\right)}=\text{E}\left\{U_i{\mathbf{b}}_i{\mathbf{Z}}_i^{\top }\right|{\mathbf{y}}_i,{\widehat{\mathit{\theta }}}^{\left(k\right)}\}$$, $$ {\widehat{u_2{\mathbf{b}}_2}}_i^{\left(k\right)}=\text{E}\left\{U_i{\mathbf{b}}_i{\mathbf{b}}_i^{\top }\right|{\mathbf{y}}_i,{\widehat{\mathit{\theta }}}^{\left(k\right)}\}$$, $$ {\widehat{vu}}_i^{\left(k\right)}=\text{E}\left\{V_i{U}_i^{1/2}\right|{\mathbf{y}}_i,{\widehat{\mathit{\theta }}}^{\left(k\right)}\}$$, $$ {\widehat{v_2}}_i^{\left(k\right)}=\text{E}\left\{V_i^2\right|{\mathbf{y}}_i,{\widehat{\mathit{\theta }}}^{\left(k\right)}\}$$, $$ {\widehat{u_2}}_i^{\left(k\right)}=\text{E}\left\{U_i\right|{\mathbf{y}}_i,{\widehat{\mathit{\theta }}}^{\left(k\right)}\}$$, $$ {\widehat{u_2\mathbf{z}}}_i^{\left(k\right)}=\text{E}\left\{U_i{\mathbf{Z}}_i\right|{\mathbf{y}}_i,{\widehat{\mathit{\theta }}}^{\left(k\right)}\}$$, $$ {\widehat{vu\mathbf{z}}}_i^{\left(k\right)}=\text{E}\left\{V_i{U}_i^{1/2}{\mathbf{Z}}_i\right|{\mathbf{y}}_i,{\widehat{\mathit{\theta }}}^{\left(k\right)}\}$$, $$ {\widehat{u_2{\mathbf{z}}_2}}_i^{\left(k\right)}=\text{E}\left\{U_i{\mathbf{Z}}_i{\mathbf{Z}}_i^{\top }\right|{\mathbf{y}}_i,{\widehat{\mathit{\theta }}}^{\left(k\right)}\}$$ and $$ {\widehat{\text{log}{u}_2}}_i^{\left(k\right)}=\text{E}\left\{\text{log}{U}_i\right|{\mathbf{y}}_i,{\widehat{\mathit{\theta }}}^{\left(k\right)}\}$$. The expressions for these expectations (except $$ {\widehat{\text{log}{u}_2}}_i^{\left(k\right)}$$) are given in the following result where we have dropped the supraindex (k) for simplicity.

Proposition 4 (see S7 Appendix for a proof). Consider the random variables Y_ij, Z_i, b_i, U_i, and V_i as defined in Eq (27) with υ₀ = 1, and an update $$ \widehat{\mathit{\theta }}={({\widehat{\mathit{\beta }}}^{\top },{\widehat{\delta }}_{\epsilon },{\widehat{\mathit{\delta }}}^{\top },{vech\left(\widehat{\overline{\mathbf{D}}}\right)}^{\top })}^{\top }$$ of the model parameter θ. Let $$ {\widehat{\overline{\mathbf{\Omega }}}}_i={\mathbf{I}}_{n_i}+{\mathbf{W}}_i\widehat{\overline{\mathbf{D}}}{\mathbf{W}}_i^{\top }$$, $$ {\widehat{\mathbf{r}}}_i=\widehat{\overline{\mathbf{D}}}{\mathbf{W}}_i^{\top }{\widehat{\overline{\mathbf{\Omega }}}}_i^{-1}$$, $$ {\widehat{\Lambda }}_i=({\mathbf{I}}_q-{\widehat{\mathbf{r}}}_i{\mathbf{W}}_i)\widehat{\overline{\mathbf{D}}}$$, $$ {\widehat{\mathbf{\Delta }}}_i={\widehat{\delta }}_{\epsilon }{\mathbf{J}}_{n_i}+{\mathbf{W}}_i\widehat{\mathit{\delta }}$$, $$ {\widehat{\mathbf{s}}}_i=\widehat{\mathit{\delta }}-{\widehat{\mathbf{r}}}_i{\widehat{\mathbf{\Delta }}}_i$$, $$ {\widehat{\mathit{\mu }}}_i={\mathbf{X}}_i\widehat{\mathit{\beta }}-c{\tilde{U}}_1{\widehat{\mathbf{\Delta }}}_i$$, $$ {\widehat{M}}_i={(1+{\widehat{\mathbf{\Delta }}}_i^{\top }{\widehat{\overline{\mathbf{\Omega }}}}_i^{-1}{\widehat{\mathbf{\Delta }}}_i)}^{-1/2}$$, $$ {\widehat{\mathbf{\Omega }}}_i={\widehat{\overline{\mathbf{\Omega }}}}_i+{\widehat{\mathbf{\Delta }}}_i{\widehat{\mathbf{\Delta }}}_i^{\top }$$, $$ {\widehat{\mathbf{\lambda }}}_i={\widehat{M}}_i^{-1}{\widehat{\mathbf{\Omega }}}_i^{-1/2}{\widehat{\mathbf{\Delta }}}_i$$, $$ {{\mathbf{Z}}_0}_i={\widehat{\mathbf{\Omega }}}_i^{-1/2}({\mathbf{Z}}_i-{\widehat{\mathit{\mu }}}_i)$$. Then:

where $$ {\widehat{u_2\alpha }}_i={\widehat{M}}_i{\widehat{\mathbf{\Delta }}}_i^{\top }{\widehat{\overline{\mathbf{\Omega }}}}_i^{-1}({\widehat{u_2\mathbf{z}}}_i-{\widehat{u_2}}_i{\widehat{\mathit{\mu }}}_i)$$, $$ {\widehat{\tau \alpha }}_i={\widehat{M}}_i{\widehat{\mathbf{\Delta }}}_i^{\top }{\widehat{\overline{\mathbf{\Omega }}}}_i^{-1}({\widehat{\tau \mathbf{z}}}_i-{\widehat{\tau }}_i{\widehat{\mathit{\mu }}}_i)$$, $$ {\widehat{u_2\alpha \mathbf{z}}}_i={\widehat{M}}_i({\widehat{u_2{\mathbf{z}}_2}}_i-{\widehat{u_2\mathbf{z}}}_i{\widehat{\mathit{\mu }}}_i^{\top }){\widehat{\overline{\mathbf{\Omega }}}}_i^{-1}{\widehat{\mathbf{\Delta }}}_i$$, $$ {\widehat{u_2{\alpha }_2}}_i={\widehat{M}}_i^2[{\widehat{u_2}}_i{\left({\widehat{\mathbf{\Delta }}}_i^{\top }{\widehat{\overline{\mathbf{\Omega }}}}_i^{-1}{\widehat{\mathit{\mu }}}_i\right)}^2+{\widehat{\mathbf{\Delta }}}_i^{\top }{\widehat{\overline{\mathbf{\Omega }}}}_i^{-1}({\widehat{u_2{\mathbf{z}}_2}}_i-2{\widehat{u_2\mathbf{z}}}_i{\widehat{\mathit{\mu }}}_i^{\top }\left){\widehat{\overline{\mathbf{\Omega }}}}_i^{-1}{\widehat{\mathbf{\Delta }}}_i\right]$$, and the expectations $$ {\widehat{u_2\mathbf{z}}}_i=\text{E}\left\{U_i{\mathbf{Z}}_i\right|{\mathbf{y}}_i,\widehat{\mathit{\theta }}\}$$, $$ {\widehat{u_2{\mathbf{z}}_2}}_i=\text{E}\left\{U_i{\mathbf{Z}}_i{\mathbf{Z}}_i^{\top }\right|{\mathbf{y}}_i,\widehat{\mathit{\theta }}\}$$, $$ {\widehat{\tau }}_i=\text{E}\left\{U_i^{1/2}{\zeta }_1\left(U_i^{1/2}{\mathbf{\lambda }}^{\top }{\mathbf{Z}}_{0i}\right)\right|{\mathbf{y}}_i,\widehat{\mathit{\theta }}\}$$, $$ {\widehat{\tau \mathbf{z}}}_i=\text{E}\left\{U_i^{1/2}{\zeta }_1\left(U_i^{1/2}{\mathbf{\lambda }}^{\top }{\mathbf{Z}}_{0i}\right){\mathbf{Z}}_i\right|{\mathbf{y}}_i,\widehat{\mathit{\theta }}\}$$ and $$ {\widehat{u_2}}_i=\text{E}\left\{U_i\right|{\mathbf{y}}_i,\widehat{\mathit{\theta }}\}$$ are to be evaluated directly using Corollary 1 applied to the conditional latent vector $$ {\mathbf{Z}}_i|{\mathbf{Y}}_i={\mathbf{y}}_i,\mathit{\theta }=\widehat{\mathit{\theta }}\sim {\mathcal{\text{TST}}}_{n_i}({\widehat{\mathit{\mu }}}_i,{\widehat{\mathbf{\Omega }}}_i,{\widehat{\mathbf{\lambda }}}_i,\widehat{\nu },{\mathbb{A}}_i)$$, $$ {\mathbb{A}}_i={\mathbb{A}}_{i1}\times \dots \times {\mathbb{A}}_{i{n}_i}$$ with $$ {\mathbb{A}}_{ij}=(-\infty ,0]$$ if y_ij = 0 and $$ {\mathbb{A}}_{ij}=(0,\infty )$$ if y_ij = 1.

The M-step jointly maximizes $$ Q(·|{\widehat{\mathit{\theta }}}^{\left(k\right)})$$ over θ. This yields the following updating expressions for θ.

Proposition 5 (see S8 Appendix for a proof) Consider an identifiable SGTLM as defined in Eq (27) with υ₀ = 1; and an estimate $$ {\widehat{\mathit{\theta }}}^{\left(k\right)}$$ of θ. Set $$ {\widehat{{\mathbf{S}}_1}}_i^{\left(k\right)}={\widehat{u_2\mathbf{z}}}_i^{\left(k\right)}-{\mathbf{W}}_i{\widehat{u_2\mathbf{b}}}_i^{\left(k\right)}$$, $$ {\widehat{S_2}}_i^{\left(k\right)}={\widehat{vu}}_i^{\left(k\right)}-c{\tilde{U}}_1{\widehat{u_2}}_i^{\left(k\right)}$$, $$ {\widehat{S_3}}_i^{\left(k\right)}={\widehat{v_2}}_i^{\left(k\right)}-2c{\tilde{U}}_1{\widehat{vu}}_i^{\left(k\right)}+c^2{\tilde{U}}_1^{2\left(k\right)}{\widehat{u_2}}_i^{\left(k\right)}$$, $$ {\widehat{{\mathbf{S}}_4}}_i^{\left(k\right)}={\widehat{vu\mathbf{z}}}_i^{\left(k\right)}-{\mathbf{W}}_i{\widehat{vu\mathbf{b}}}_i^{\left(k\right)}$$, $$ {\widehat{\mathbf{T}}}_1^{\left(k\right)}={\left[{\sum }_{i=1}^n{\widehat{u_2}}_i^{\left(k\right)}{\mathbf{X}}_i^{\top }{\mathbf{X}}_i\right]}^{-1}$$ and $$ {\widehat{\mathbf{T}}}_2^{\left(k\right)}={\sum }_{i=1}^n{\widehat{S_2}}_i^{\left(k\right)}{\mathbf{J}}_{n_i}^{\top }{\mathbf{X}}_i$$. At EM iteration (k + 1), the updates of β, δ_ε, δ and $$ \overline{\mathbf{D}}$$ are given by:

At convergence of the EM algorithm, we obtain the estimate $$ \tilde{\mathit{\theta }}\left(\mathit{\nu }\right)$$ of θ. The corresponding estimate of the variance-covariance matrix of random effects is $$ {\widehat{\Sigma }}_b={\tilde{U}}_2\widehat{\overline{\mathbf{D}}}+({\tilde{U}}_2-c^2{\tilde{U}}_1^2)\widehat{\mathit{\delta }}{\widehat{\mathit{\delta }}}^{\top }$$. More generally, when ν is not actually known, the M-step of the EM algorithm can be extended to include a profiled marginal log-likelihood maximization step. Indeed, at EM iteration k, we notice that the estimate $$ {\tilde{\mathit{\theta }}}^{\left(k\right)}\left(\nu \right)$$ of θ depends on ν only through $$ {\tilde{U}}_1$$ and the profiled marginal log-likelihood L_ν(⋅) for ν can be obtained by simply substituting $$ {\tilde{\mathit{\theta }}}^{\left(k\right)}\left(\nu \right)$$ for θ in Eq (36). We can thus find $$ {\widehat{\nu }}^{(k+1)}$$ using a one dimensional optimization routine (e.g. optimize in R) to maximize L_ν(⋅). Then, the update of the parameter θ becomes $$ {\widehat{\mathit{\theta }}}^{(k+1)}={\tilde{\mathit{\theta }}}^{\left(k\right)}\left({\widehat{\nu }}^{(k+1)}\right)$$. The use of the profiled marginal log-likelihood instead of a profiled version of $$ Q(·|{\widehat{\mathit{\theta }}}^{\left(k\right)})$$ can provide substantial gain of efficiency [50] and mostly helps bypass the calculation of $$ {\widehat{\text{log}{u}_2}}_i^{\left(k\right)}$$ which does not have any known closed form.

It is worthwhile noticing however that the inclusion of a profiled marginal log-likelihood maximization step would prevent the convergence of the whole estimation procedure if the marginal log-likelihood in Eq (36) as a function of ν is unbounded for a particular dataset. This issue especially when ν is close to zero. Another challenge associated to the estimation of ν is time. The use of this strategy requires a very fast routine to compute cumulative probabilities of the skew t distributions. As an alternative route to estimate ν, we point out the model selection approach of [51] (page 893). It consists in setting a grid of feasible values of ν and obtaining a sequence of estimates $$ \tilde{\mathit{\theta }}\left(\nu \right)$$ of θ. Then, the couple ν and $$ \tilde{\mathit{\theta }}\left(\nu \right)$$ maximizing the marginal log-like-lihood in Eq (36) is taken as the estimates of ν and θ.