Likelihood-based Inference for Skewed Responses in a Crossover Trial Setup

Savita Pareek1, Kalyan Das1, and Siuli Mukhopadhyay1, 222Corresponding author.
Email<EMAIL_ADDRESS>
1Department of Mathematics, Indian Institute of Technology Bombay,
Mumbai 400 076, India

## Abstract

This work proposes a statistical model for crossover trials with multiple
skewed responses measured in each period. A 3 $\times$ 3 crossover trial data
where different doses of a drug were administered to subjects with a history
of seasonal asthma rhinitis to grass pollen is used for motivation. In each
period, gene expression values for ten genes were measured from each subject.
It considers a linear mixed effect model with skew normally distributed random
effect or random error term to model the asymmetric responses in the crossover
trials. The paper examines cases (i) when a random effect follows a skew-
normal distribution, as well as (ii) when a random error follows a skew-normal
distribution. The EM algorithm is used in both cases to compute maximum
likelihood estimates of parameters. Simulations and crossover data from the
gene expression study illustrate the proposed approach.

Keywords. Crossover design, Mixed effect models, Skew-normal distribution, EM
algorithm.

##  1 Introduction

Crossover design is a specific type of longitudinal study in which every
subject receives different treatments in different periods. It is most
commonly used in the pharmaceutical industry and other medical fields to
investigate the safety and efficacy of new drugs or treatments. In a crossover
design, the treatment effect is compared within the subject since each subject
acts as its own control. Therefore, fewer subjects are needed than in a
parallel design in order to achieve the same level of accuracy. An in-depth
description of crossover trials can be found in the books by Senn 32, Jones
and Kenward 22.

In many clinical studies, we sometimes encounter crossover trials with
measurements on two or more response variates. For example, one may consider
the measurement of both systolic (SBP) and diastolic (DBP) blood pressure of
subjects in each period (Grender and Johnson 19) or blood sugar levels
recorded at multiple time points in each period (Putt and Chinchilli 29) or
microarray gene expression profiles of subjects measured in each period
(Leaker et al. 25). Other than multiple responses in some real-life crossover
trials, such as bioequivalence trials, we come across multi-modal or skewed
responses (Jones and Kenward 22). One solution to such cases is to apply the
log or box-cox power transformation. However, this may not be a reasonable
solution for most cases, thus causing a considerable modeling challenge for
the statistician (Feng et al. 16). The use of Bayesian methods based on
extended generalized gamma distribution and skew-t distribution has been
discussed in the literature for such bioequivalence studies with the skewed
response (de Souza et al. 12, Burger et al. 9). Using conventional analysis
techniques for such skewed multivariate data may lead to an incorrect and
biased parameter and variance estimates.

In comparison to univariate responses, crossover trials with multiple
responses measured in each period have been addressed by very few researchers,
namely (Grender and Johnson 19, Chinchilli and Esinhart 10, Putt and
Chinchilli 29, Tudor et al. 34, Johnson and Mercante 21, Pareek et al. 27).
The main reason behind the scarcity of literature may be the difficulty in
modeling the complex inter and intra-response relationships. The linear mixed
or random-effects model (Laird and Ware 24) is a widely used technique for
analyzing data from crossover studies. These models routinely assume both
random effects and random errors to be normally distributed. While the
assumption of normality is easy to execute, several authors (Zhang and
Davidian 37, Ghidey et al. 17, Drikvandi et al. 15 and Drikvandi 14) have
questioned the robustness of fitted models, specifically when data show multi-
modality and skewness. Zhang and Davidian 37 have demonstrated that in the
Framingham heart study, estimated subject-specific intercepts are not normally
distributed, and the use of normal distribution in this scenario has resulted
in less efficient inferences on intercept and treatment effect or subject
level covariate. Drikvandi et al. 15 has argued that maximum likelihood
estimates of regression parameters may be biased when the distribution of
random effects is not correctly specified. Therefore, it is of practical
interest to develop statistical models having random effects or random errors
to be skew-normal while allowing for multiple response measurements in each
period. In the last decade, a substantial amount of work has been done on
skewed responses in longitudinal studies. Arellano-Valle et al. 3 have
developed an EM-type algorithm for maximum likelihood estimation in linear
mixed effect (LME) models, assuming random effects or random errors to be
multivariate skew-normal. Lachos et al. 23 proposed the longitudinal data
modeling using linear mixed models with skew-normal independent (SNI)
distribution for random effects and normal independent (NI) distributions of
within-subject errors. SNI distributions are a sub-class of scale mixture of
skew-normal (SMSN) distribution introduced by Branco and Dey 8. Recently
Pereira and Russo 28 have developed non-linear regression models assuming
random effects to be a class of SMSN such as skew-normal, skew-t, skew-slash,
and skew-contaminated normal distributions. Schumacher et al. 31 presented the
maximum likelihood estimation using an EM-type algorithm for correlated error
terms such as the auto-regressive correlation of order p in the LME model with
SNI random effects. However, none of these authors have considered crossover
trials with multiple and skewed measurements in each period.

In this work, we address the issue of skewed responses in a multivariate
crossover setup. A mixed-effect model approach with an EM-based estimation
method is proposed. We consider the linear mixed-effects model with skew-
normally distributed random subject-specific effect or random error term. The
hierarchical representation of the model makes it feasible to use the EM-type
algorithm, which produces the closed-form expressions for E and M-steps for
special cases. We further analyze the gene expression data (Leaker et al. 25)
from a three-period three-treatment crossover design.

The rest of the article is structured as follows: after a brief discussion on
the gene expression case study in Section 2, the proposed random effects model
is presented in Section 3. In Section 4, a likelihood-based EM-type algorithm
is considered for parameter estimation. In Section 5, two simulation studies
are conducted to examine the performance of parameter estimates. The proposed
methodology’s benefits are illustrated by analyzing a gene expression dataset
in Section 6, while Section 7 provides the relevant computational specifics.
Finally, concluding remarks are provided in Section 8.

##  2 Case Study: Multivariate Crossover Trial of Oral Prednisone

We use a gene expression dataset from (Leaker et al. 25) as a case study. The
dataset is publicly available from the NCBI Gene Expression Omnibus (Clough
and Barrett 11) and can be accessed using the hyperlink, nasal mRNA data. In
the gene expression study, results from a randomized double-blind, placebo-
controlled, three-period, crossover trial are considered to evaluate the
effects of two single doses of oral prednisone (10 mg, 25 mg) on inflammatory
mediators measured in nasal exudates after nasal allergen challenge in
susceptible individuals with allergic rhinitis. All subjects have a history of
seasonal asthma rhinitis to grass pollen and a positive result from the
intraepidermal skin prick test to grass pollen extract. Seventeen subjects
were enrolled in the study and assigned to three treatment sequences/groups,
out of which five subject observations with missing values were not
considered. The main interest here is to study the effect of treatments and
genes on allergic reactions to grass pollen. The outcomes measured are a fold
change of mRNA expression levels (pg ml-1), i.e., changes in gene expression
values for ten genes recorded in the nasal allergen challenge. Subjects with
missing observations are excluded from our analysis. The study design is
described in Table 1.

Table 1: Study design: a three-way crossover trial to examine the effects of a single oral dose of prednisone (10 or 25 mg) versus placebo given before nasal allergen challenge (NAC). | Period 1 | Period 2 | Period 3
---|---|---|---
Sequence 1 (4 subjects) | 10mg | Placebo | 25 mg
Sequence 2 (4 subjects) | 25 mg | 10 mg | Placebo
Sequence 3 (4 subjects) | Placebo | 25 mg | 10 mg

As an example, sequence 1 indicates that 10 mg prednisone is given to the
subjects (subject 1 to subject 4) in period 1, followed by a washout, a
placebo is given in period 2, and again after a washout, 25 mg prednisone is
given in period 3.

Before model fitting, we ran some exploratory analysis on the gene data as
follows:

  * (i)

The density and normal Q-Q plot of the original and log-transformed responses
as given in Figure 1. We determine if the given responses are representative
of a normal population by using the Shapiro-Wilk tests (Shapiro and Wilk 33).
Table 2 interprets raw and transformed responses’ density, normal Q-Q plots,
and Shapiro-Wilk tests at 5% significance level. The observations made in
Table 2 indicate asymmetric behavior of the gene expression levels. Moreover,
the Shapiro-Wilk test on box-cox transformed responses yields a p-value
$<0.0001$, indicating that the box-cox power transformation also fails to
conform the responses to normality.

Table 2: Testing for normality using the density, Q-Q plot, and Shapiro-Wilk test for original and log transformed responses. Method | Original responses | Log-transformed responses
---|---|---
| Shapiro-Wilk test
---
(p-value)
<0.0001 | <0.0001
Density plot |  | Does not resemble
---
a bell-shaped curve.
| Does not resemble
---
a bell-shaped curve.
Normal Q-Q plot |  | Both ends of the Q-Q plot
---
deviates significantly
from the diagonal line.
| Both ends of the Q-Q plot
---
deviates significantly
from the diagonal line.
Figure 1: Density and normal Q-Q plot comparisons: observed raw responses
(Plots A and B) and log-transformed responses (Plots C and D), with true
normal density indicated by a dashed line.

  * (ii)

The objective is to evaluate the effect of different genes on the relationship
between period versus responses, treatment versus responses, and subject
versus responses, respectively. This can be achieved by creating interaction
plots. According to Figure 2, plots A and B primarily have parallel lines.
Plot C shows a slight overlap between the lines for genes 2, 9, and 10.
Nevertheless, our model does not consider subject versus gene interactions in
order to reduce model complexity and parameters, but it may be possible to
examine such interactions in the future.

Figure 2:  Interaction Plots: gene interactions by time period, treatment, and
subject. Plot A represents gene interactions based on time periods (Pi, i=1,
2, 3), Plot B represents gene interactions based on treatments (Ti, i=1, 2,
3), and Plot C represents gene interactions based on subjects (Si, i=1, 2, …,
12). Each point in Plots A, B, and C displays the average response of a given
gene (averaged over sequences).

  * (iii)

As a further step, a normal linear mixed-effects model was fitted, with
period, treatment, and gene as fixed effects and subject-specific normal
random effects. As a result of fitting the model, we constructed a Q-Q plot of
the estimated random intercepts, as well as a residual versus estimated
responses plot. From Figure 3, we see that the estimated subject-specific
intercept histogram, Q-Q plot, and the Shapiro-Wilk test (having p-value
$0.33$) depict no apparent non-normal patterns. However, the normal Q-Q plot
and the Shapiro-Wilk test (p-value $<0.0001$) of the standardized residuals
show asymmetric behavior, and the residual versus fitted values plot also
indicates non-constant variance.

Figure 3: Validation graphs for normal model fitting based on estimated random
effects and standardized residuals: Subplots A and B represent the density and
normal Q-Q Plot of the estimated random intercepts, while Subplots C and D
depict the standardized residuals vs. fitted values and normal Q-Q Plots of
the standardized residuals, with the true normal density displayed as a dashed
line.

Based on the exploratory analysis, we propose the following skewed-normal
model for the multiple gene responses.

##  3 Model and Notations

In this section, we propose a statistical model motivated by the gene
expression data to fit multiple skewed responses measured in a crossover
trial. Suppose, $y_{ijtk}$ denotes the $k^{\text{th}}$ gene expression at
$t^{\text{th}}$ time period for $j^{\text{th}}$ subject in $i^{\text{th}}$
sequence, where $i=1,2,\ldots,s$; $j=1,2,\ldots,n_{i}$; $t=1,2,\ldots,p$ and
$k=1,2,\ldots,m$. There are in total $n$ subjects, where
$n=\sum_{i=1}^{s}n_{i}$. We may write

$y_{ijtk}=\mu+\pi_{t}+\tau_{d[i,t]}+g_{k}+s_{ij}+e_{ijtk},$ (1)

where, $\mu$ is the overall mean, $\pi_{t}$ is the $t^{\text{th}}$ period
effect, $\tau_{d[i,t]}$ is the effect of treatment applied in $t^{\text{th}}$
period of $i^{\text{th}}$ sequence, $g_{k}$ is the $k^{\text{th}}$ gene
effect, $s_{ij}$ is the subject-specific random effect, and $e_{ijtk}$ is the
random error.

For ease of exposition, we rewrite eq. 1 in matrix form as

$\bm{y}_{ij}=\bm{X}_{ij}\bm{\beta}+\bm{Z}_{ij}b_{ij}+\bm{e}_{ij};~{}i=1(1)s,~{}j=1(1)n_{i}.$
(2)

For fixed values of i and j,

  * (i)

$\bm{y}_{ij}=(y_{ij11},\ldots,y_{ij1m},y_{ij21},\ldots,y_{ij2m},\cdots,y_{ijp1},\ldots,y_{ijpm})^{T}$
is a response vector of length $pm$.

  * (ii)

$\bm{X}_{ij}$ is the design matrix of order $pm\times(p+t+m-2)$ corresponding
to the fixed effects, where,

$\bm{X}_{ij}=(\bm{1}_{pm},(\bm{0}_{m\times(p-1)},(\bm{I}_{p-1}\otimes\bm{1}_{m})^{T})^{T}_{pm\times(p-1)},\bm{T}_{pm\times(t-1)},\bm{1}_{p}\otimes(\bm{0}_{m-1},\bm{I}_{m-1})^{T}).$

The ${pm\times(t-1)}$ matrix of treatment effects is represented as
$\bm{T}=(\bm{T}_{1}^{T},\ldots\bm{T}_{p}^{T})^{T}$, where,
$\bm{T}_{u}=(\bm{1}_{m}\otimes(a_{2},\ldots,a_{t}))_{m\times(t-1)}$ for
$u=1(1)p$, and $a_{l}$ for $l=2(1)p$ is an indicator variable which takes
value $1$ if the $l^{\text{th}}$ treatment is assigned to the $u^{\text{th}}$
period and $0$ otherwise.

  * (iii)

$\bm{\beta}$ is the parameter vector of length $(p+t+m-2)$ corresponding to
the fixed effects,
$\bm{\beta}=(\mu,\pi_{2},\dots\pi_{p},\tau_{2},\ldots,\tau_{t},g_{2},\ldots,g_{m}).$

  * (iv)

$\bm{Z}_{ij}=\bm{1}_{pm}$ is the design matrix corresponding to the random
effects. Consider, for example, two responses being measured for each subject
in every period of a $2\times 2$ crossover design with treatment sequences AB
and BA having five subjects in each. Then for $i=1,2;j=1,2,\ldots,5$, the
design matrices $\bm{X}_{ij},\bm{Z}_{ij}$ and parameter vector $\bm{\beta}$
will be

$\bm{X}_{1j}=\begin{bmatrix}1&0&0&0\\\ 1&0&0&1\\\ 1&1&1&0\\\
1&1&1&1\end{bmatrix},\>\quad\bm{X}_{2j}=\begin{bmatrix}1&0&1&0\\\ 1&0&1&1\\\
1&1&0&0\\\ 1&1&0&1\end{bmatrix},\>\quad\bm{Z}_{ij}=\begin{bmatrix}1\\\ 1\\\
1\\\ 1\end{bmatrix},$ $\bm{\beta}=(\mu,\pi_{2},\tau_{2},g_{2}).$

  * (v)

Moreover, ${b}_{ij}=s_{ij}$ is the subject-specific random effect, and
$\bm{e}_{ij}=(e_{ij11},\ldots,e_{ijpm})^{T}$, is the random error vector. We
assume an independent error structure to reduce the number of parameters and
for ease of calculation in the maximum likelihood estimation. In our
computations, we consider two scenarios:

    * Scenario 1: error is skew-normal

$\bm{e}_{ij}\overset{ind}{\sim}SN_{pm}(\bm{0},\sigma_{e}^{2}\bm{I}_{pm},(\lambda_{e},\ldots,0));\;\;{b}_{ij}\overset{ind}{\sim}N(0,\sigma_{s}^{2}).$

    * Scenario 2: random effect is skew-normal

$\bm{e}_{ij}\overset{ind}{\sim}N_{pm}(\bm{0},\sigma_{e}^{2}\bm{I}_{pm});\;\;{b}_{ij}\overset{ind}{\sim}SN(0,\sigma_{s}^{2},\lambda_{s}).$

  * (vi)

The parameters of interest are,
$\bm{\theta}=(\bm{\beta},\sigma_{e}^{2},\sigma_{s}^{2},\lambda_{e},\lambda_{s}).$

For ease of the readers, we briefly define the univariate and multivariate
skew normal distribution and their parameters used in the following sections
(more details are given in Section 9.1).

  * (a)

Univariate skew-normal (SN) variate (Azzalini 4): If a random variable $W$ has
the density function

$f(w;\lambda)=2\phi(w)\Phi(\lambda w);\;\;-\infty<w<\infty,$

where $\phi$, $\Phi$ are the standard normal density and distribution
function, respectively, then we say $W$ is a skew-normal variate with
parameter $\lambda$, or, $W\sim SN(0,1,\lambda)$.

  * (b)

Multivariate skew-normal (SN) variate (Arellano-Valle et al. 3): An $n$
dimensional random vector $\bm{W}_{1}$ follows a skew-normal distribution with
location vector $\bm{\mu}\in\mathbb{R}^{n}$, a dispersion matrix $\bm{\Sigma}$
(a $n\times n$ positive definite matrix) and skewness vector
$\bm{\lambda}\in\mathbb{R}^{n}$, if its pdf is given by,

$f_{\bm{W}_{1}}(\bm{w}_{1})=2\phi_{n}(\bm{w}_{1}|\bm{\mu},\bm{\Sigma})\Phi_{1}(\bm{\lambda}^{T}\bm{\Sigma}^{-1/2}(\bm{w}_{1}-\bm{\mu})),\;\;\bm{w}_{1}\in\mathbb{R}^{n}.$

We denote it by $\bm{W}_{1}\sim SN_{n}(\bm{\mu},\bm{\Sigma},\bm{\lambda})$,
further $\bm{W}_{1}$ can be written as
$\bm{W}_{1}{=}\bm{\mu}+\bm{\Sigma}^{1/2}\bm{W}$, where, $\bm{W}$ is a
standardized multivariate skew-normal vector. Throughout this article, skew-
normal distributions are denoted by the letter SN, whereas normal
distributions are identified by the letter N.

##  4 Maximum Likelihood Estimation

As discussed in Verbeke and Molenberghs 35, parameter estimation in the models
as specified in eq. 2 is based on the marginal distribution of the response
$\bm{y}_{ij}$ unless they are analyzed in the Bayesian framework. However, to
compute the marginal distribution of $\bm{y}_{ij}$ we need to use complex
Monte Carlo integration since it involves the skew-normal distribution. There
is also no explicit method for directly maximizing the marginal likelihood,
and one has to resort to numerical maximization (Arellano-Valle et al. 3).

In this work, we use the EM algorithm (Dempster et al. 13, Arellano-Valle et
al. 3, 23), a generic iterative approach for maximum likelihood estimation in
models with random effects and incomplete data. Every iteration of the EM
algorithm has two steps: the Expectation Step (E-step) and the Maximization
Step (M-step), which increases the likelihood function and typically converges
to the local or global maximum of the likelihood function (Wu 36).

In this section, we discuss parameter estimation using the EM algorithm for
two specific cases, (i) when the error is SN, and (ii) when the random effect
is SN. The necessary results for the execution of the expectation step of the
EM algorithm can be found in Section 9.2. The corresponding maximization step
is solved using the first-order Newton-Raphson equation.

###  4.1 ML Estimation when Errors are SN

In the case of skew-normal errors, in eq. 2, we assume

$\bm{e}_{ij}\sim
SN_{pm}(\textbf{0},\bm{V}_{1ij},\bm{\lambda}_{er}),\;\text{where,}\;\bm{V}_{1ij}=\sigma_{e}^{2}\bm{I}_{pm},\;\text{and}\;\bm{\lambda}_{er}=(\lambda_{e},0,0,\ldots,0).$

Using Lemma 9.1 from the Section 9.2, we have

$\bm{e}_{ij}\overset{d}{=}\bm{V}_{1ij}^{1/2}\bm{\delta}_{e}|U_{0ij}|+\bm{V}_{1ij}^{1/2}(\textbf{I}_{pm}-\bm{\delta}_{e}\bm{\delta}_{e}^{T})^{1/2}\textbf{U}_{1ij},\;\;$

where,
$\bm{\delta}_{e}=\frac{\bm{\lambda}_{er}}{\sqrt{1+\bm{\lambda}_{er}^{T}\bm{\lambda}_{er}}}$,
$U_{0ij}\sim N(0,1)$ independently of $\bm{U}_{1ij}\sim
N_{n}(\bm{0},\bm{I}_{n})$. Thus, eq. 2 can be written as

$\bm{y}_{ij}=\bm{X}_{ij}\bm{\beta}+\bm{d}_{ij}t_{ij}+\bm{r}_{ij},$ (3)

where,

$\bm{d}_{ij}=\bm{V}_{1ij}^{1/2}\bm{\delta}_{e},\;\;t_{ij}=|U_{0ij}|,\;t_{ij}\overset{iid}{\sim}HN_{1}(0,1),$

and,
$\bm{r}_{ij}=\bm{Z}_{ij}{b}_{ij}+\bm{V}_{1ij}^{1/2}(\textbf{I}_{pm}-\bm{\delta}_{e}\bm{\delta}_{e}^{T})^{1/2}\textbf{U}_{1ij},$
follows:

$\bm{r}_{ij}\overset{ind}{\sim}N_{pm}(\textbf{0},\bm{V}_{ij}),\;\text{for}\;\bm{V}_{ij}=\sigma_{s}^{2}\bm{Z}_{ij}\bm{Z}_{ij}^{T}+\sigma_{e}^{2}\bm{R},\;\;\bm{R}=\textbf{I}_{pm}-\bm{\delta}_{e}\bm{\delta}_{e}^{T}.$

Moreover, for $i=1(1)s$ and $j=1(1)n_{i}$, $t_{ij}\perp\bm{r}_{ij}$. Since
$\bm{r}_{ij}$ has a zero mean, it can be used in residual analyses to assess
the adequacy of the model. The right side of eq. 3 has a mean
$\sigma_{e}\bm{\delta}_{e}\sqrt{2/\pi}$, which is used to correct the
intercept in the fixed effects. Under the above setup, the conditional model
where conditioning refers to $t_{ij}$ takes the form

$\bm{y}_{ij}|t_{ij}\sim
N_{pm}(\bm{\mu}_{ij}^{*},\bm{V}_{ij}),\;\bm{\mu}_{ij}^{*}=\bm{X}_{ij}\bm{\beta}+\bm{d}_{ij}t_{ij}.$

The complete data log-likelihood can then be expressed as

$\displaystyle l_{c}(\bm{\theta}_{c})$
$\displaystyle=\log(\prod_{i=1}^{s}\prod_{j=1}^{n_{i}}f(\bm{y}_{ij},t_{ij}))=\sum_{i=1}^{s}\sum_{j=1}^{n_{i}}\left(\log
f(\bm{y}_{ij}|t_{ij})+\log f(t_{ij})\right)$
$\displaystyle\propto-\frac{1}{2}\sum_{i=1}^{s}\sum_{j=1}^{n_{i}}\left(\log|\bm{V}_{ij}|+(\bm{y}_{ij}-\bm{\mu}_{ij}^{*})^{T}{\bm{V}_{ij}}^{-1}(\bm{y}_{ij}-\bm{\mu}_{ij}^{*})+t_{ij}^{2}\right).$
(4)

The E-step consists of calculating the expected value of the complete data
log-likelihood given the observed data and current parameter estimates.

For fixed i, j, we have

$\displaystyle Q_{ij}(\bm{\theta}_{c}|\bm{\theta}_{c}^{(r)})$
$\displaystyle=E_{t_{ij}|\bm{y}_{ij},\bm{\theta}_{c}^{(r)}}[l_{c}(\bm{\theta}_{c})]=\int
l_{c}(\bm{\theta}_{c})f(t_{ij}|\bm{y}_{ij},\bm{\theta}_{c}^{(r)})dt_{ij}.$ (5)

Using $l_{c}(\bm{\theta}_{c})$ from eq. 4 in eq. 5, the expression for
$Q_{ij}$ is as follows,

$\displaystyle Q_{ij}(\bm{\theta}_{c}|\bm{\theta}_{c}^{(r)})$
$\displaystyle=-\frac{1}{2}[\log|\bm{V}_{ij}|+E_{t_{ij}|\bm{y}_{ij},\bm{\theta}_{c}^{(r)}}(t_{ij}^{2})$
$\displaystyle+E_{t_{ij}|\bm{y}_{ij},\bm{\theta}_{c}^{(r)}}(\bm{y}_{ij}-\bm{X}_{ij}\bm{\beta}-\bm{d}_{ij}t_{ij})^{T}{\bm{V}_{ij}}^{-1}(\bm{y}_{ij}-\bm{X}_{ij}\bm{\beta}-\bm{d}_{ij}t_{ij})]$
$\displaystyle=-\frac{1}{2}[\log|\bm{V}_{ij}|+T_{ij}^{02}+\bm{d}_{ij}^{T}{\bm{V}_{ij}}^{-1}\bm{d}_{ij}T_{ij}^{02}$
$\displaystyle+(\bm{y}_{ij}-\bm{X}_{ij}\bm{\beta})^{T}{\bm{V}_{ij}}^{-1}(\bm{y}_{ij}-\bm{X}_{ij}\bm{\beta}-2\bm{d}_{ij}T_{ij}^{01})],$

where, $T_{ij}^{01}$, and $T_{ij}^{02}$ are obtained as follows using Lemmas
9.2, and 9.3 from the Section 9.2;

$T_{ij}^{01}=E[t_{ij}|\bm{y}_{ij},\bm{\theta}_{c}^{(r)}]=\eta_{ij}+\frac{\phi(\frac{\eta_{ij}}{\zeta_{ij}})}{\Phi_{1}(\frac{\eta_{ij}}{\zeta_{ij}})}\zeta_{ij}.$
$T_{ij}^{02}=E[t_{ij}^{2}|\bm{y}_{ij},\bm{\theta}_{c}^{(r)}]=\eta_{ij}^{2}+\zeta_{ij}^{2}+\frac{\phi(\frac{\eta_{ij}}{\zeta_{ij}})}{\Phi_{1}(\frac{\eta_{ij}}{\zeta_{ij}})}\eta_{ij}\zeta_{ij},$
$\text{here,}\;\eta_{ij}=\frac{\bm{d}_{ij}^{T}\bm{V}_{ij}^{-1}(\bm{y}_{ij}-\bm{X}_{ij}\bm{\beta})}{1+\bm{d}_{ij}^{T}\bm{V}_{ij}^{-1}\bm{d}_{ij}}\;\;\text{and}\;\zeta_{ij}^{2}=\frac{1}{1+\bm{d}_{ij}^{T}\bm{V}_{ij}^{-1}\bm{d}_{ij}}.$

Considering subjects are independent, the E-step yields,

$Q(\bm{\theta}_{c}|\bm{\theta}_{c}^{(r)})=\sum_{i=1}^{s}\sum_{j=1}^{n_{i}}Q_{ij}(\bm{\theta}_{c}|\bm{\theta}_{c}^{(r)}).$
(6)

In the second step of the algorithm, i.e., the M-step, our task is to maximize
$Q(\bm{\theta}_{c}|\bm{\theta}_{c}^{(r)})$, we do this in two parts. Here the
maximization of $Q-$function corresponding to the fixed effects is
straightforward and closed-form expressions are available. However, for
variance components, maximization of $Q-$function has to be done numerically
using the first-order Newton-Raphson equation. In the M-step as a first step,
we find the $(r+1)^{th}$ estimate of ${\bm{\beta}}$, $\bm{\beta}^{(r+1)}$, by
maximizing

$\displaystyle Q(\bm{\beta})$
$\displaystyle=\sum_{i=1}^{s}\sum_{j=1}^{n_{i}}-\frac{1}{2}[(\bm{y}_{ij}-\bm{X}_{ij}\bm{\beta})^{T}{\bm{V}_{ij}}^{-1}(\bm{y}_{ij}-\bm{X}_{ij}\bm{\beta}-2\bm{d}_{ij}T_{ij}^{01})],$

which yields

$\displaystyle\bm{\beta}^{(r+1)}$
$\displaystyle=(\sum_{i=1}^{s}\sum_{j=1}^{n_{i}}\bm{X}_{ij}^{T}{\bm{V}_{ij}}^{-1(r)}\bm{X}_{ij})^{-1}\sum_{i=1}^{s}\sum_{j=1}^{n_{i}}\bm{X}_{ij}^{T}{\bm{V}_{ij}}^{-1(r)}(\bm{y}_{ij}-\bm{d}_{ij}T_{ij}^{01(r)}).$

We use these estimated $\bm{\beta}^{(r+1)}$ in the second step to obtain the
variance components. The three components
$(\sigma_{e}^{2(r+1)},\sigma_{s}^{2(r+1)},\lambda_{e}^{(r+1)})$ are found
numerically by solving the first order derivatives of the $Q-$ function (eq.
6), using the Newton-Raphson method. Suppose
${\bm{\xi}}=({\xi}_{[1]},{\xi}_{[2]},{\xi}_{[3]})=(\sigma_{e}^{2},\sigma_{s}^{2},\lambda_{e})$,
then,

$\displaystyle\frac{\partial Q}{\partial{\xi}_{[v_{0}]}}$
$\displaystyle=\sum_{i=1}^{s}\sum_{j=1}^{n_{i}}-\frac{1}{2}[\text{tr}(\bm{V}_{ij}^{-1}\bm{V}_{ij{\xi}_{[v_{0}]}})+\bm{d}_{ij}^{T}{\bm{V}^{*}_{ij{\xi}_{[v_{0}]}}}\bm{d}_{ij}T_{ij}^{02}+2\bm{d}_{ij}^{T}\bm{V}_{ij}^{-1}{\bm{d}_{ij{\xi}_{[v_{0}]}}}T_{ij}^{02}$
$\displaystyle+(\bm{y}_{ij}-\bm{X}_{ij}\bm{\beta}^{(r+1)})^{T}{\bm{V}^{*}_{ij{\xi}_{[v_{0}]}}}(\bm{y}_{ij}-\bm{X}_{ij}\bm{\beta}^{(r+1)}-2\bm{d}_{ij}T_{ij}^{01})$
$\displaystyle-2(\bm{y}_{ij}-\bm{X}_{ij}\bm{\beta}^{(r+1)})^{T}\bm{V}_{ij}^{-1}{\bm{d}_{ij{\xi}_{[v_{0}]}}}T_{ij}^{01}],\;\;v_{0}=1,2,3,$

where,

$\displaystyle\bm{V}_{ij\bm{\xi}}$
$\displaystyle=(\frac{\partial\bm{V}_{ij}}{\partial\sigma_{e}^{2}},\frac{\partial\bm{V}_{ij}}{\partial\sigma_{s}^{2}},\frac{\partial\bm{V}_{ij}}{\partial\lambda_{e}})=(\bm{R},\bm{Z}_{ij}\bm{Z}_{ij}^{T},\sigma_{e}^{2}\bm{R}_{l}),$
$\displaystyle\bm{d}_{ij\bm{\xi}}$
$\displaystyle=(\frac{\partial\bm{d}_{ij}}{\partial\sigma_{e}^{2}},\frac{\partial\bm{d}_{ij}}{\partial\sigma_{s}^{2}},\frac{\partial\bm{d}_{ij}}{\partial\lambda_{e}})=(\frac{\bm{\delta}_{e}}{2\sqrt{\sigma_{e}^{2}}},\bm{0},\sigma_{e}\bm{\delta}_{el}),$
$\displaystyle\bm{V}^{*}_{ij\bm{\xi}}$
$\displaystyle=\frac{\partial\bm{V}^{-1}_{ij}}{\partial\bm{\xi}}=-\bm{V}_{ij}^{-1}\bm{V}_{ij{\bm{\xi}}}\bm{V}_{ij}^{-1},$
$\displaystyle\text{and}\;\bm{\delta}_{el}$
$\displaystyle=\frac{\partial\bm{\delta}_{e}}{\partial\lambda_{e}}=\frac{(1,0,\ldots,0)}{(1+\lambda_{e}^{2})^{3/2}},\;\;\bm{R}_{l}=\frac{\partial\bm{R}}{\partial\lambda_{e}}=-\bm{\delta}_{e}\bm{\delta}_{el}^{T}-\bm{\delta}_{el}\bm{\delta}_{e}^{T}.$

For the parameter vector
$\bm{\xi}=(\sigma_{e}^{2},\sigma_{s}^{2},\lambda_{e})$, the updated estimates
at $(r+1)^{\text{th}}$ iteration are then given by

$\bm{\xi}^{(r+1)}=\bm{\xi}^{(r)}-\bm{H}^{-1(r)}\frac{\partial
Q}{\partial\bm{\xi}}\Bigr{|}_{\bm{\xi}=\bm{\xi}^{(r)}},\;\text{where,}\;\bm{H}=\frac{\partial^{2}Q}{\partial\bm{\xi}\partial\bm{\xi}^{{}^{\prime}}}.$

Detailed computations of the $\bm{H}$ matrix are provided in Section 9.3. We
define the convergence criterion as the difference between estimated values at
$(r+1)^{\text{th}}$ iteration and $(r)^{\text{th}}$ iteration being less than
$5\times 10^{-3}$. The E and M-steps are iterated until convergence.

###  4.2 ML Estimation when Random Effects are SN

For skew-normal random effects, i.e., when ${b}_{ij}\sim
SN(0,\sigma_{s}^{2},\lambda_{s})$, applying Lemma 9.1 of Section 9.2, we can
express ${b}_{ij}$ in eq. 2 as

${b}_{ij}\overset{d}{=}\sigma_{s}{\delta}_{b}|U_{0ij}|+\sigma_{s}(1-\delta_{b}^{2})^{1/2}{U}_{1ij},\;\;$

where, ${\delta}_{b}=\frac{{\lambda_{s}}}{\sqrt{1+\lambda_{s}^{2}}}$,
$U_{0ij},\;\text{and}\;U_{1ij}$ are independent standard normals. As a result,
eq. 2 becomes

$\bm{y}_{ij}=\bm{X}_{ij}\bm{\beta}+\bm{d}_{ij}t_{ij}+\bm{r}_{ij},$ (7)

where,

$\bm{d}_{ij}=\bm{Z}_{ij}\sigma_{s}\delta_{b};\;t_{ij}=|U_{0ij}|,\
\text{and}\;t_{ij}\overset{iid}{\sim}HN_{1}(0,1)$
$\text{For}\;\bm{r}_{ij}=\bm{Z}_{ij}\sigma_{s}(1-\delta_{b}^{2})^{1/2}\bm{U}_{1ij}+\bm{e}_{ij},$

we see that,

$\bm{r}_{ij}\overset{ind}{\sim}N_{pm}(\textbf{0},\bm{V}_{ij}),\;\text{where}\;\bm{V}_{ij}=\sigma_{s}^{2}\bm{Z}_{ij}R\bm{Z}_{ij}^{T}+\sigma_{e}^{2}\bm{I}_{pm},\;\text{and}\;R=1-\delta_{b}^{2}.$

Further, $t_{ij}\perp\bm{r}_{ij}$; $i=1(1)s$ and $j=1(1)n_{i}$. Equation 7 has
a mean $\bm{Z}_{ij}\sigma_{s}\delta_{b}$ that is used to correct the intercept
in the fixed effects. Under the above setting, the conditional model where
conditioning refers to $t_{ij}$ takes the form,

$\bm{y}_{ij}|t_{ij}\sim
N_{pm}(\bm{\mu}_{ij}^{*},\bm{V}_{ij}),\;\bm{\mu}_{ij}^{*}=\bm{X}_{ij}\bm{\beta}+\bm{d}_{ij}t_{ij}.$

As in Section 4.1, the complete data log-likelihood can be expressed as
follows:

$\displaystyle l_{c}(\bm{\theta}_{c})$
$\displaystyle\propto-\frac{1}{2}\sum_{i=1}^{s}\sum_{j=1}^{n_{i}}\left(\log|\bm{V}_{ij}|+(\bm{y}_{ij}-\bm{\mu}_{ij}^{*})^{T}{\bm{V}_{ij}}^{-1}(\bm{y}_{ij}-\bm{\mu}_{ij}^{*})+t_{ij}^{2}\right).$

The E-step, as discussed for the skewed random error case, consists of
calculating the expected value of complete data log-likelihood given the
observed data and current parameter estimates.

$\displaystyle Q_{ij}(\bm{\theta}_{c}|\bm{\theta}_{c}^{(r)})$
$\displaystyle=-\frac{1}{2}[\log|\bm{V}_{ij}|+T_{ij}^{02}+\bm{d}_{ij}^{T}{\bm{V}_{ij}}^{-1}\bm{d}_{ij}T_{ij}^{02}$
$\displaystyle+(\bm{y}_{ij}-\bm{X}_{ij}\bm{\beta})^{T}{\bm{V}_{ij}}^{-1}(\bm{y}_{ij}-\bm{X}_{ij}\bm{\beta}-2\bm{d}_{ij}T_{ij}^{01})]$

where, $T_{ij}^{01}=E[t_{ij}|\bm{y}_{ij},\bm{\theta}_{c}^{(r)}]$,
$T_{ij}^{02}=E[t_{ij}^{2}|\bm{y}_{ij},\bm{\theta}_{c}^{(r)}].$ The E-step for
all the subjects yields,

$Q(\bm{\theta}_{c}|\bm{\theta}_{c}^{(r)})=\sum_{i=1}^{s}\sum_{j=1}^{n_{i}}Q_{ij}(\bm{\theta}_{c}|\bm{\theta}_{c}^{(r)}).$
(8)

Similar to M-step in Section 4.1, we maximize
$Q(\bm{\theta}_{c}|\bm{\theta}_{c}^{(r)})$ in two parts. In the first part we
find $\bm{\beta}^{(r+1)}$ by maximising eq. 8 with respect to $\bm{\beta}$,

$\displaystyle\bm{\beta}^{(r+1)}$
$\displaystyle=(\sum_{i=1}^{s}\sum_{j=1}^{n_{i}}\bm{X}_{ij}^{T}{\bm{V}_{ij}}^{-1(r)}\bm{X}_{ij})^{-1}\sum_{i=1}^{s}\sum_{j=1}^{n_{i}}\bm{X}_{ij}^{T}{\bm{V}_{ij}}^{-1(r)}(\bm{y}_{ij}-\bm{d}_{ij}T_{ij}^{01(r)}).$

While in the second part, we numerically solve for
$(\sigma_{e}^{2(r+1)},\sigma_{s}^{2(r+1)},\lambda_{s}^{(r+1)})$ using the
Newton-Raphson method. Based on Section 4.1,
${\bm{\xi}}=({\xi}_{[1]},{\xi}_{[2]},{\xi}_{[3]})=(\sigma_{e}^{2},\sigma_{s}^{2},\lambda_{s})$,
then

$\displaystyle\frac{\partial Q}{\partial{\xi}_{[v_{0}]}}$
$\displaystyle=\sum_{i=1}^{s}\sum_{j=1}^{n_{i}}-\frac{1}{2}[\text{tr}(\bm{V}_{ij}^{-1}\bm{V}_{ij{\xi}_{[v_{0}]}})+\bm{d}_{ij}^{T}{\bm{V}^{*}_{ij{\xi}_{[v_{0}]}}}\bm{d}_{ij}T_{ij}^{02}+2\bm{d}_{ij}^{T}\bm{V}_{ij}^{-1}{\bm{d}_{ij{\xi}_{[v_{0}]}}}T_{ij}^{02}$
$\displaystyle+(\bm{y}_{ij}-\bm{X}_{ij}\bm{\beta})^{T}{\bm{V}^{*}_{ij{\xi}_{[v_{0}]}}}(\bm{y}_{ij}-\bm{X}_{ij}\bm{\beta}-2\bm{d}_{ij}T_{ij}^{01})$
$\displaystyle-2(\bm{y}_{ij}-\bm{X}_{ij}\bm{\beta})^{T}\bm{V}_{ij}^{-1}{\bm{d}_{ij{\xi}_{[v_{0}]}}}T_{ij}^{01}],v_{0}=1,2,3,$

where,

$\displaystyle\bm{V}_{ij\bm{\xi}}$
$\displaystyle=(\frac{\partial\bm{V}_{ij}}{\partial\sigma_{e}^{2}},\frac{\partial\bm{V}_{ij}}{\partial\sigma_{s}^{2}},\frac{\partial\bm{V}_{ij}}{\partial\lambda_{s}})=(\bm{I}_{pm},R\bm{Z}_{ij}\bm{Z}_{ij}^{T},\sigma_{s}^{2}R_{l}\bm{Z}_{ij}\bm{Z}_{ij}^{T}),$
$\displaystyle\bm{d}_{ij\bm{\xi}}$
$\displaystyle=(\frac{\partial\bm{d}_{ij}}{\partial\sigma_{e}^{2}},\frac{\partial\bm{d}_{ij}}{\partial\sigma_{s}^{2}},\frac{\partial\bm{d}_{ij}}{\partial\lambda_{s}})=(\bm{0},\frac{\bm{Z}_{ij}\delta_{b}}{2\sqrt{\sigma_{s}^{2}}},\bm{Z}_{ij}\sigma_{s}\delta_{bl}),$
$\displaystyle\bm{V}^{*}_{ij\bm{\xi}}$
$\displaystyle=\frac{\partial\bm{V}^{-1}_{ij}}{\partial\bm{{\xi}}}=-\bm{V}_{ij}^{-1}\bm{V}_{ij{\bm{\xi}}}\bm{V}_{ij}^{-1},$
$\displaystyle\delta_{bl}$
$\displaystyle=\frac{\partial\delta_{b}}{\partial\lambda_{s}}=\frac{1}{(1+\lambda_{s}^{2})^{3/2}},\;\;\text{and}\;R_{l}=\frac{\partial
R}{\partial\lambda_{s}}=-2\delta_{b}\delta_{bl}.$

For the parameter vector
$\bm{\xi}=(\sigma_{e}^{2},\sigma_{s}^{2},\lambda_{s})$, the updated estimates
at the $(r+1)^{\text{th}}$ iteration are given by

$\bm{\xi}^{(r+1)}=\bm{\xi}^{(r)}-\bm{H}^{-1(r)}\frac{\partial
Q}{\partial\bm{\xi}}\Bigr{|}_{\bm{\xi}=\bm{\xi}^{(r)}},\bm{H}=\frac{\partial^{2}Q}{\partial\bm{\xi}\partial\bm{\xi}^{{}^{\prime}}}.$

The elements of the $\bm{H}$ matrix are given in Section 9.4. The E and
M-steps are iterated till convergence is achieved.

##  5 Simulation Studies

To assess the performance of the proposed estimators, we present simulation
studies for the two cases, (i) errors are skew-normal and (ii) random effects
are skew-normal. For data generation in both cases, we assume a crossover
trial with three treatment sequences $\\{ABC,BCA,CAB\\}$ in three periods. Two
simulation scenarios with 30 and 50 subjects assigned, respectively to each
treatment sequence are considered. In each period, four response variates are
measured. The model is represented as,

$\displaystyle y_{ijtk}$
$\displaystyle=\beta_{0}+\beta_{p_{2}}\text{Per}_{2}+\beta_{p_{3}}\text{Per}_{3}+\beta_{\tau_{2}}\text{Trt}_{2}+\beta_{\tau_{3}}\text{Trt}_{3}+\beta_{g_{2}}\text{Gene}_{2}+\beta_{g_{3}}\text{Gene}_{3}$
$\displaystyle+\beta_{g_{4}}\text{Gene}_{4}+\beta_{1}w_{ij}+s_{ij}+e_{ijtk};\;i,t=1,2,3;\;j=1(1)n_{i};\;k=1(1)4,$

where $y_{ijtk}$ denotes the $k^{\text{th}}$ response value from the
$j^{\text{th}}$ subject in the $t^{\text{th}}$ period of the $i^{\text{th}}$
sequence; $\beta_{0}$: intercept; Per2, Per3: indicator variables
corresponding to time/period effects, Trt2, Trt3: indicator variables
corresponding to the treatment effects, and Gene2, Gene3, Gene4: indicator
variables corresponding to a subject’s gene expression level. Following are
the indicated variables with respect to the period, treatment, and gene
effects for $u=2,3$ and $v=2,3,4$,

$\text{Per}_{u}/\text{Trt}_{u}=\begin{cases}1,&\text{if $u^{\text{th}}$
period/treatment}\\\ 0,&\text{otherwise}\end{cases}$
$\text{Gene}_{v}=\begin{cases}1,&\text{if $v^{\text{th}}$ gene expression
level}\\\ 0,&\text{otherwise}\end{cases}$

$w_{ij}$ represents an individual-level covariate taking values in
$\\{0,1,2\\}$. In scenario 1, when there are 30 subjects in each sequence,
$w_{ij}$ takes the value 0 for subjects 1 to 10, 1 for subjects 11 to 20, and
2 for subjects 21 to 30. Considering scenario 2, when $n_{i}=50$ for each
sequence, $w_{ij}$ is 0 for subjects 1 to 18, 1 for subjects 19 to 34, and 2
for subjects 35 to subject 50. Also, $s_{ij}$ and $e_{ijtk}$ are the subject-
specific random effect and the random error terms, respectively. In matrix
notations,

$\bm{y}_{ij}=\bm{X}_{ij}\bm{\beta}+\bm{Z}_{ij}{b}_{ij}+\bm{e}_{ij};\;\;i=1,2,3,\;\;j=1,2,\ldots,n_{i},$

where, $\bm{X}_{ij}$ and $\bm{Z}_{ij}$ matrices are as described in eq. 2. The
true values of the components of $\bm{\beta}$ are given in the first column of
Tables 3 and 4, respectively. For our simulations, we take the following true
values for the variance components:

  * (i)

For the case where the errors are skew-normal,

$b_{ij}\sim N(0,0.64),\bm{e}_{ij}\sim
SN_{pm}(\bm{0},2\bm{I}_{pm},(3,0,\ldots,0)).$

Hence, $\sigma_{e}^{2}=2,\sigma_{s}^{2}=0.64,$ and $\lambda_{e}=3.$

  * (ii)

For the skewed random effect case,

$b_{ij}\sim SN(0,3,4),\bm{e}_{ij}\sim N_{pm}(\bm{0},0.72\bm{I}_{pm}).$

Thus, $\sigma_{e}^{2}=0.72,\sigma_{s}^{2}=3,$ and $\lambda_{s}=4.$

Two hundred Monte Carlo data sets were generated for each of the above
parameter settings. The model fitting results under parameter settings (i) and
(ii) are compared with the case where we assume both random effects and random
errors to be normally distributed. The Akaike Information Criterion (AIC)
(Glosup and Axelrod 18) was applied to select the best-fitting model.

After the estimation of all of the model parameters, the fitted distribution
of $\bm{r}_{ij}$ given in eq. 3, eq. 7 can be plotted along with the data
scatter in order to assess the model’s adequacy. In an alternative approach to
evaluating model fitting, Mahalanobis-type distances may be used, which were
first introduced by Healy 20 for multivariate normal distribution distances.
Let $\bm{y}=(\bm{y}_{1},\ldots,\bm{y}_{m})$ where the $i^{\text{th}}$
component, $\bm{y}_{i}$, is sampled from
$SN_{n}(\bm{\mu},\bm{\Sigma},\bm{\lambda})$, the Mahalanobis-type distances
are defined as,

$d_{i}=(\bm{y}_{i}-\bm{\hat{\mu}})^{T}\bm{\hat{\Sigma}}^{-1}(\bm{y}_{i}-\bm{\hat{\mu}}),\;\;i=1,2,\ldots,m,$

whose approximate reference distribution is $\chi^{2}_{n}$ (Azzalini and
Capitanio 6. From these $d_{i}$’s, we obtain QQ-plots. Plotting nominal
probability values against the theoretical cumulative probabilities of the
observed Mahalanobis distances. These plots are also called a Healy-type plot
(Schumacher et al. 31). An appropriately fitted model should produce a
straight line with a unit slope through the origin in a Healy-type plot.

Tables 3 and 4 show the average simulation results in terms of parameter
estimates, standard errors (SEs), and absolute bias. The SEs are estimated
using the Hessian matrix, and the average absolute bias for the true value
$\beta_{u}$ is computed as
$\frac{\sum_{w=1}^{200}|\hat{\beta}_{uw}-\beta_{u}|}{200}$, where
$\hat{\beta}_{uw}$ is the $u^{\text{th}}$ component of $\hat{\bm{\beta}}$ for
the $w^{\text{th}}$ simulation.

Table 3: Simulation study: results based on 200 Monte Carlo samples with a different number of subjects. Maximum likelihood estimates, SEs, and absolute bias for the cases $b_{ij}\sim N(0,0.6),\;\bm{e}_{ij}\sim SN_{pm}(\bm{0},2\bm{I}_{pm},(3,0,\ldots,0))$ and when both $\bm{e}_{ij}$, $b_{ij}$ are N. True parameter values are given in parentheses.  |  If $\bm{e}_{ij}$ is SN, $b_{ij}$ is N |  Both $\bm{e}_{ij}$ , $b_{ij}$ are N
---|---|---
Parameter | Estimate | SE | $\left|\text{Bias}\right|$ | Estimate | SE | $\left|\text{Bias}\right|$
$n_{i}$=30 for each sequence
$\beta_{0}$(2.1) | 2.0769 | 0.0505 | 0.1459 | 1.5482 | 0.1846 | 0.5225
$\beta_{p_{2}}$(2.4) | 2.3925 | 0.0767 | 0.0811 | 2.1271 | 0.1036 | 0.2729
$\beta_{p_{3}}$(1.1) | 1.0826 | 0.0797 | 0.0807 | 0.8172 | 0.1036 | 0.2828
$\beta_{\tau_{2}}$(0.9) | 0.8999 | 0.0793 | 0.0770 | 0.9026 | 0.1036 | 0.0810
$\beta_{\tau_{3}}$(2.1) | 2.1118 | 0.0533 | 0.0737 | 2.1130 | 0.1036 | 0.0793
$\beta_{g_{2}}$(1.5) | 1.4750 | 0.0820 | 0.0882 | 1.1211 | 0.1197 | 0.3789
$\beta_{g_{3}}$(2.0) | 1.9880 | 0.0928 | 0.0864 | 1.6341 | 0.1197 | 0.3659
$\beta_{g_{4}}$(3.4) | 3.3853 | 0.0661 | 0.0814 | 3.0314 | 0.1197 | 0.3686
$\beta_{1}$(1.8) | 1.7959 | 0.0501 | 0.1049 | 1.7966 | 0.1136 | 0.1045
$\sigma_{e}^{2}$(2.0) | 1.9774 | 0.1182 | 0.0746 | 1.9343 | 1.3904 | 0.0914
$\sigma_{s}^{2}$(0.6) | 0.6216 | 0.0461 | 0.0922 | 0.6166 | 0.7821 | 0.0939
$\lambda_{e}$(3.0) | 3.2824 | 0.2384 | 0.8926 | - | - | -
$n_{i}$=50 for each sequence
$\beta_{0}$(2.1) | 2.0690 | 0.0322 | 0.1120 | 1.5314 | 0.1410 | 0.5391
$\beta_{p_{2}}$(2.4) | 2.4009 | 0.0603 | 0.0642 | 2.1345 | 0.0806 | 0.2655
$\beta_{p_{3}}$(1.1) | 1.0919 | 0.0634 | 0.0624 | 0.8254 | 0.0806 | 0.2746
$\beta_{\tau_{2}}$(0.9) | 0.9039 | 0.0619 | 0.0670 | 0.9058 | 0.0806 | 0.0666
$\beta_{\tau_{3}}$(2.1) | 2.1062 | 0.0413 | 0.0621 | 2.1058 | 0.0806 | 0.0648
$\beta_{g_{2}}$(1.5) | 1.4886 | 0.0651 | 0.0668 | 1.1333 | 0.0930 | 0.3667
$\beta_{g_{3}}$(2.0) | 2.0075 | 0.0728 | 0.0802 | 1.6522 | 0.0930 | 0.3478
$\beta_{g_{4}}$(3.4) | 3.4049 | 0.0507 | 0.0726 | 3.0497 | 0.0930 | 0.3503
$\beta_{1}$(1.8) | 1.8011 | 0.0414 | 0.0708 | 1.8015 | 0.0877 | 0.0695
$\sigma_{e}^{2}$(2.0) | 1.9906 | 0.0867 | 0.0580 | 1.9471 | 1.3952 | 0.0720
$\sigma_{s}^{2}$(0.6) | 0.6260 | 0.0388 | 0.0726 | 0.6221 | 0.7867 | 0.0740
$\lambda_{e}$(3.0) | 3.4737 | 0.1884 | 0.8793 | - | - | -
Table 4: Simulation study: results based on 200 Monte Carlo samples with a different number of subjects. Maximum likelihood estimates, SEs, and absolute bias for the cases $b_{ij}\sim SN(0,3,4),\;\bm{e}_{ij}\sim N_{pm}(\bm{0},0.7\bm{I}_{pm})$ and when both $\bm{e}_{ij}$, $b_{ij}$ are N. True parameter values are given in parentheses. |  If $\bm{e}_{ij}$ is N, $b_{ij}$ is SN |  Both $\bm{e}_{ij}$ , $b_{ij}$ are N
---|---|---
Parameter | Estimate | SE | $\left|\text{Bias}\right|$ | Estimate | SE | $\left|\text{Bias}\right|$
$n_{i}$=30 for each sequence
$\beta_{0}$(3.3) | 3.3548 | 0.1023 | 0.1535 | 3.3597 | 0.1983 | 0.1584
$\beta_{p_{2}}$(2.4) | 2.4004 | 0.0515 | 0.0497 | 2.4004 | 0.0632 | 0.0497
$\beta_{p_{3}}$(1.1) | 1.1023 | 0.0610 | 0.0509 | 1.1023 | 0.0632 | 0.0509
$\beta_{\tau_{2}}$(0.9) | 0.8993 | 0.0617 | 0.0544 | 0.8993 | 0.0632 | 0.0544
$\beta_{\tau_{3}}$(2.1) | 2.0979 | 0.0545 | 0.0562 | 2.0979 | 0.0632 | 0.0562
$\beta_{g_{2}}$(1.5) | 1.5017 | 0.0693 | 0.0598 | 1.5017 | 0.0730 | 0.0598
$\beta_{g_{3}}$(2.0) | 2.0008 | 0.0663 | 0.0598 | 2.0008 | 0.0730 | 0.0598
$\beta_{g_{4}}$(3.4) | 3.3996 | 0.0513 | 0.0623 | 3.3996 | 0.0730 | 0.0623
$\beta_{1}$(1.8) | 1.7865 | 0.0580 | 0.1021 | 1.7846 | 0.1442 | 0.1132
$\sigma_{e}^{2}$(0.7) | 0.7205 | 0.0110 | 0.0288 | 0.7205 | 0.8485 | 0.0288
$\sigma_{s}^{2}$(3.0) | 2.9948 | 0.1923 | 0.5095 | 1.1948 | 1.0894 | 1.8052
$\lambda_{s}$(4.0) | 4.2827 | 0.1771 | 1.2641 | - | - | -
$n_{i}$=50 for each sequence
$\beta_{0}$(3.3) | 3.3617 | 0.0759 | 0.1124 | 3.3668 | 0.1495 | 0.1239
$\beta_{p_{2}}$(2.4) | 2.3965 | 0.0400 | 0.0380 | 2.3965 | 0.0491 | 0.0380
$\beta_{p_{3}}$(1.1) | 1.1003 | 0.0473 | 0.0378 | 1.1003 | 0.0491 | 0.0378
$\beta_{\tau_{2}}$(0.9) | 0.8993 | 0.0479 | 0.0433 | 0.8993 | 0.0491 | 0.0433
$\beta_{\tau_{3}}$(2.1) | 2.0973 | 0.0423 | 0.0414 | 2.0973 | 0.0491 | 0.0414
$\beta_{g_{2}}$(1.5) | 1.5051 | 0.0538 | 0.0463 | 1.5051 | 0.0567 | 0.0463
$\beta_{g_{3}}$(2.0) | 1.9985 | 0.0515 | 0.0482 | 1.9985 | 0.0567 | 0.0482
$\beta_{g_{4}}$(3.4) | 3.3974 | 0.0398 | 0.0418 | 3.3974 | 0.0567 | 0.0418
$\beta_{1}$(1.8) | 1.7868 | 0.0441 | 0.0750 | 1.7811 | 0.1105 | 0.0901
$\sigma_{e}^{2}$(0.7) | 0.7230 | 0.0086 | 0.0190 | 0.7230 | 0.8502 | 0.0190
$\sigma_{s}^{2}$(3.0) | 2.9621 | 0.1472 | 0.4180 | 1.1872 | 1.0867 | 1.8128
$\lambda_{s}$(4.0) | 4.2618 | 0.1381 | 1.1209 | - | - | -

The results in Table 3 demonstrate that the estimates for bias and standard
errors are lower when using the SN assumption, compared to assuming both the
random error and random effect are normal. The decrease in bias and standard
errors is more pronounced as the number of subjects per sequence increases to
50. Additionally, the proposed skew-normal model was deemed the most
appropriate fit by 89% of the AIC values when the actual random error was
determined to be skew-normal, suggesting that the proposed algorithm is
successful in handling skewed crossover data.

The results shown in Table 4 indicate that the parameter estimates for period,
treatment, and gene effects in the skew-normal model are comparable to those
in the normal model, however, their standard errors differ. Specifically, the
standard error for $\beta_{g_{4}}$ and $\beta_{1}$ in the skew-normal case is
approximately half of that in the normal case. As previously reported by other
researchers (Lachos et al. 23), the estimate of the covariate $w_{ij}$’s
parameter $\beta_{1}$ is impacted by both standard error and bias.
Additionally, when the random effect is actually skew-normal, 83% of AIC
values suggest that the skew-normal fit is the best model, demonstrating that
fitting standard model results in less accurate estimates for skewed crossover
data.

The accuracy of the proposed SN model is evaluated by analyzing the chi-square
Q-Q plot of the Mahalanobis distances ($\bm{d}_{ij}$) and the normal Q-Q plot
of the standardized residual vectors ($\bm{r}_{ij}$). The top panels of
Figures 4 and 5 show that the quantiles of the proposed model are closely
aligned with the diagonal reference line. Furthermore, the empirical density
plot of the standardized maximum likelihood estimate is created, and the
bottom panels of Figures 4 and 5 indicate that they conform to the normal
distribution curve. This leads us to the conclusion that the proposed model is
appropriate for the data at hand.

Figure 4: Assesing model adequacy: a comparison of SN model ($\bm{e}_{ij}$ is
SN, and $b_{ij}$ is N) with the normal model for 30 subject sequences. The top
panel compares the Mahalanobis distances to the theoretical $\chi^{2}_{12}$
distribution and the standardized residual vectors to a standard normal
distribution. The ‘N fit’ refers to a model fitted using normal random errors
and random effects, while the ‘SN fit’ refers to the proposed model. The
bottom panel displays the estimated fixed effects through empirical density
plots. Figure 5: Assesing model adequacy: a comparison of SN model
($\bm{e}_{ij}$ is N and $b_{ij}$ is SN) with the normal model for 30 subject
sequences. The top panel compares the Mahalanobis distances to the theoretical
$\chi^{2}_{12}$ distribution and the standardized residual vectors to a
standard normal distribution. The ‘N fit’ refers to a model fitted using
normal random errors and random effects, while the ‘SN fit’ refers to the
proposed model. The bottom panel displays the estimated fixed effects through
empirical density plots.

##  6 Gene Case Study Results

In our preliminary investigation, we discovered that both the raw and log-
transformed responses were skewed. We hypothesized that the period and
treatment effects would be uniform across the ten gene expression levels and
that specific interaction terms could be disregarded. Additionally, utilizing
a model with a non-normal error distribution may result in a more favorable
outcome for fitting the gene data. Thus, we fit the following statistical
model with the underlying regression equation as given in eq. 9:

$\displaystyle y_{ijtk}$
$\displaystyle=\beta_{0}+\beta_{p_{2}}\text{Per}_{2}+\beta_{p_{3}}\text{Per}_{3}+\beta_{\tau_{2}}\mathrm{Trt}_{2}+\beta_{\tau_{3}}\mathrm{Trt}_{3}+\beta_{g_{2}}\text{Gene}_{2}$
$\displaystyle+\beta_{g_{3}}\text{Gene}_{3}+\beta_{g_{4}}\text{Gene}_{4}+\beta_{g_{5}}\text{Gene}_{5}+\beta_{g_{6}}\text{Gene}_{6}+\beta_{g_{7}}\text{Gene}_{7}$
$\displaystyle+\beta_{g_{8}}\text{Gene}_{8}+\beta_{g_{9}}\text{Gene}_{9}+\beta_{g_{10}}\text{Gene}_{10}+s_{ij}+e_{ijtk},$
(9)

where for the gene data $i,t=1,2,3;\,k=1(1)10$ and $j=1(1)n_{i}$, and
$n_{1}=n_{2}=n_{3}=4$. Scenarios with different assumptions on $s_{ij}$ and
$e_{ijtk}$ were used. The period, treatment, and gene effects were assumed to
be fixed and represented by corresponding indicator variables. The AIC and
Bayesian information criterion (BIC) were used to compare the various cases.
These cases are,

  * Case 1: A model with normally distributed random error and random effects, i.e., when both $b_{ij}$ and $\bm{e}_{ij}$ are N.

  * Case 2: A model with independent multivariate normal distribution for the random error and a univariate skew-normal distribution for random effects, i.e., when $b_{ij}$ is SN and $\bm{e}_{ij}$ is N.

  * Case 3: A model with independent multivariate skew-normal distribution for random error and univariate normal distribution for random effects, i.e., when $b_{ij}$ is N and $\bm{e}_{ij}$ is SN.

Table 5 reports the maximum likelihood estimates and estimated asymptotic
standard errors of the parameters
$(\bm{\beta},\sigma_{e}^{2},\sigma_{s}^{2},\lambda)$ for the three cases. As
an initial value of mean and variance parameters, estimates obtained by
fitting a normal linear mixed effect model were used in the EM algorithm. From
Table 5, we observe that parameter estimates and standard errors of Case 2 are
close to the estimates obtained under the normality assumption (Case 1),
implying that the asymmetry is not detected in random effects. However, Case 3
suggests the best fit based on AIC and BIC values, supporting the argument of
departure from the normality of the residuals.

The adequacy of the selected model in Case 3 is assessed by examining the
Mahalanobis distances $\bm{d_{ij}}$. A Kolmogorov–Smirnov (KS) test (Massey
26) is most commonly used to test the goodness of fit of data to a theoretical
distribution. KS test of Mahalanobis distances $\bm{d_{ij}}$ gives a p-value
of 0.06, which indicates the calculated distances are generated from the chi-
square distribution at the 5% significance level. Furthermore, we have
constructed chi-square Q-Q plots for the Mahalanobis distances $\bm{d_{ij}}$
by which we can compare the observed and expected values of the Mahalanobis
distances. Figure 6 shows that the observed values are close to the expected
values, indicating that the model defined in Case 3 is appropriately fitted.

Figure 6: Model fit comparison for Case 3: Mahalanobis Distances vs. Theoretical $\chi^{2}_{30}$ Distribution. Table 5: Analysis of gene data: Maximum likelihood estimates, SEs, based on fitting cases 1 to 3. Comparing cases 1, 2, and 3 using AIC and BIC values.  |  Case 1 (both $b_{ij}$ , $\bm{e}_{ij}$ is N) |  Case 2 ($b_{ij}$ is SN, $\bm{e}_{ij}$ is N) |  Case 3 ($b_{ij}$ is N, $\bm{e}_{ij}$ is SN)
---|---|---|---
Parameter | Estimate | SE | Estimate | SE | Estimate | SE
$\beta_{0}$ | 3.4896 | 0.022 | 3.4896 | 0.0187 | 3.4486 | 0.0069
$\beta_{p_{2}}$ | 0.0007 | 0.0138 | 0.0007 | 0.0138 | 0.0117 | 0.0128
$\beta_{p_{3}}$ | 0.0118 | 0.0138 | 0.0118 | 0.0138 | 0.0228 | 0.0127
$\beta_{\tau_{2}}$ | -0.0055 | 0.0138 | -0.0055 | 0.0138 | -0.0088 | 0.0117
$\beta_{\tau_{3}}$ | 0.0081 | 0.0138 | 0.0081 | 0.0138 | 0.0016 | 0.0118
$\beta_{g_{2}}$ | -0.9284 | 0.0252 | -0.9284 | 0.0239 | -0.8915 | 0.0188
$\beta_{g_{3}}$ | -1.5865 | 0.0252 | -1.5865 | 0.021 | -1.5496 | 0.0173
$\beta_{g_{4}}$ | -0.3275 | 0.0252 | -0.3275 | 0.0251 | -0.2906 | 0.0182
$\beta_{g_{5}}$ | -0.0381 | 0.0252 | -0.0381 | 0.0252 | -0.0012 | 0.0171
$\beta_{g_{6}}$ | -1.2646 | 0.0252 | -1.2646 | 0.0226 | -1.2278 | 0.0183
$\beta_{g_{7}}$ | -1.4755 | 0.0252 | -1.4755 | 0.0216 | -1.4386 | 0.0177
$\beta_{g_{8}}$ | -2.0328 | 0.0252 | -2.0328 | 0.0178 | -1.9959 | 0.0146
$\beta_{g_{9}}$ | -0.8211 | 0.0252 | -0.8211 | 0.0242 | -0.7842 | 0.0188
$\beta_{g_{10}}$ | -0.6482 | 0.0252 | -0.6482 | 0.0246 | -0.6113 | 0.0187
$\sigma_{e}^{2}$ | 0.0115 | 0.0219 | 0.0115 | 0.0001 | 0.0117 | 0.0001
$\sigma_{s}^{2}$ | 0.0005 | 0.107 | 0.0005 | 0.0004 | 0.0003 | 0.0002
$\lambda$ | - | - | 0.001 | 0.3856 | 3.8849 | 0.0496
AIC | -545.50 | -528.08 | -553.64
BIC | -483.32 | -465.90 | -491.46

##  7 Computational Details

For simulation and real data analyses, R programming (R Core Team 30) has been
used. We used R version 4.1.0 under Windows 10 (64-bit), with an Intel core I5
processor and 4GB of RAM.

The simulation specifications outlined in section Section 6 were utilized to
create 100 simulated data sets, each containing 30 subjects per sequence. It
took approximately 1.8 hours to run the simulation when $b_{ij}$ was
designated as SN and $\bm{e}_{ij}$ was N. The EM algorithm converged in an
average of 45-50 steps in these data sets. In the scenario where $b_{ij}$ was
designated as N and $\bm{e}_{ij}$ was SN, the simulation took approximately 45
minutes to run and the EM algorithm converged in an average of 15-20 steps.
The normal linear mixed model was then fit to the data using the ‘lme’
function from the nlme library in R.

The gene data exploratory analysis (as shown in Section 2) and model adequacy
plots (as shown in Section 5) were depicted through Q-Q plots, and empirical
density plots using the ggplot2, ggpubr, and tidyverse libraries. In gene data
analysis, the EM algorithm converged in 2 steps for Case 2, where $b_{ij}$ is
SN and $\bm{e}_{ij}$ is N, whereas 14 steps were required for convergence in
Case 3, where $b_{ij}$ is N and $\bm{e}_{ij}$ is SN. Using the R platform with
similar settings, a gene data model can be fitted in approximately 4-5
minutes. The R programs used are available at the provided hyperlink,
rprograms.

##  8 Concluding Remarks

We have proposed an EM algorithm-based estimation for multiple and skewed
crossover data by moderating the assumptions on random effect and model error
densities. Closed-form expression are obtained for fixed effects, and variance
and skewness components are obtained using the first-order NR equation. A
small simulation study is carried out to highlight the potential gain in the
efficiency of some parameters when the normality assumption does not hold,
with some extra computing cost.

The approach used in this paper can be used in treating other multivariate
models, such as having different treatment, period, and gene effects for each
response variate. We believe this idea is applicable by assuming both the
random effect and random error to be skew normally distributed and having the
general variance-covariance structure. It is also possible to consider
dependence structures in random errors, such as damped exponential
correlations and serially autoregressive correlations of order $p$. As a
consequence, the estimation problem becomes more complex. One can also apply
the proposed technique to study other classes of asymmetric distributions such
as skew-t, skew-slash, and skew-contaminated normal.

Acknowledgement We thank Dr. Atanu Bhattacharjee, Tata Memorial Center,
Mumbai, India, for his assistance in obtaining the gene data set.

## References

  * Arellano-Valle and Azzalini [2009] R.B. Arellano-Valle and A. Azzalini.  The centred parametrization for the multivariate skew-normal distribution.  _Journal of Multivariate Analysis_ , 100(4):816, 2009.
  * Arellano-Valle and Genton [2005] R.B. Arellano-Valle and Marc G. Genton.  On fundamental skew distributions.  _Journal of Multivariate Analysis_ , 96(1):93–116, 2005.
  * Arellano-Valle et al. [2005] R.B. Arellano-Valle, H. Bolfarine, and V.H. Lachos.  Skew-normal Linear Mixed Models.  _Journal of Data Science_ , 3:415–438, 2005.
  * Azzalini [1985] A. Azzalini.  A Class of Distributions Which Includes the Normal Ones.  _Scandinavian Journal of Statistics_ , 12(2):171–178, 1985.
  * Azzalini and Capitanio [1999] A. Azzalini and A. Capitanio.  Statistical applications of the multivariate skew normal distribution.  _Journal of the Royal Statistical Society. Series B: Statistical Methodology_ , 61(3):579–602, 1999.
  * Azzalini and Capitanio [2014] A. Azzalini and A. Capitanio.  _The Skew-Normal and Related Families_.  Cambridge, 2014.
  * Azzalini and Dalla Valle [1996] A. Azzalini and A. Dalla Valle.  The multivariate skew-normal distribution.  _Biometrika_ , 83(4):715–726, 1996.
  * Branco and Dey [2001] Márcia D. Branco and Dipak K. Dey.  A general class of multivariate skew-elliptical distributions.  _Journal of Multivariate Analysis_ , 79(1):99–113, 2001.
  * Burger et al. [2021] Divan Aristo Burger, Robert Schall, and Sean van der Merwe.  A robust method for the assessment of average bioequivalence in the presence of outliers and skewness.  _Pharmaceutical Research_ , 38(10):1697–1709, 2021.
  * Chinchilli and Esinhart [1996] Vernon M. Chinchilli and James D. Esinhart.  Design and analysis of intra-subject variability in cross-over experiments.  _Statistics in Medicine_ , 15:1619–1634, 1996.
  * Clough and Barrett [2016] Emily Clough and Tanya Barrett.  The Gene Expression Omnibus database.  _Methods in Molecular Biology_ , 1418:93–110, 2016.
  * de Souza et al. [2016] Roberto Molina de Souza, Jorge Alberto Achcar, Edson Zangiacomi Martinez, and Josmar Mazucheli.  The use of asymmetric distributions in average bioequivalence.  _Statistics in Medicine_ , 35(15):2525–2542, 2016\.
  * Dempster et al. [1977] A.P. Dempster, N.M. Laird, and D.B. Rubin.  Maximum likelihood from incomplete data via the EM algorithm.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 39:1–38, 1977.
  * Drikvandi [2019] Reza Drikvandi.  Nonlinear mixed-effects models with misspecified random-effects distribution.  _Pharmaceutical Statistics_ , 19(3):187–201, 2019\.
  * Drikvandi et al. [2016] Reza Drikvandi, Geert Verbeke, and Geert Molenberghs.  Diagnosing mis- specification of the random-effects distribution in mixed models.  _Biometrics_ , 73(1):63–71, 2016.
  * Feng et al. [2014] Changyong Feng, Hongyue Wang, Naiji Lu, Tian Chen, Hua He, Ying Lu, and Xin M. Tu.  Log-transformation and its implications for data analysis.  _Shanghai Archives of Psychiatry_ , 26(2):105–109, 2014.
  * Ghidey et al. [2004] Wendimagegn Ghidey, Emmanuel Lesaffre, and Paul Eilers.  Smooth Random Effects Distribution in a Linear Mixed Model.  _Biometrics_ , 60(4):945–953, 2004.
  * Glosup and Axelrod [1994] J.G. Glosup and M.C. Axelrod.  Use of the AIC with the EM Algorithm: A Demonstration of a Probability Model Selection Technique.  In _Joint Statistical Meeting_ , 1994.
  * Grender and Johnson [1993] Julie M. Grender and William D. Johnson.  Analysis of crossover designs with multivariate response.  _Statistics in Medicine_ , 12(1):69–89, 1993\.
  * Healy [1968] M.J.R. Healy.  Multivariate Normal Plotting.  _Journal of the Royal Statistical Society. Series C (Applied Statistics)_ , 17(2):157–161, 1968.
  * Johnson and Mercante [1996] W.D. Johnson and D.E. Mercante.  Analyzing multivariate data in crossover designs using permutation tests.  _Journal of Biopharmaceutical Statistics_ , 6(3):327–342, 1996.
  * Jones and Kenward [2003] Byron Jones and Michael G. Kenward.  _Design and Analysis of Cross-Over Trials_.  Chapman & Hall/CRC, second edition, 2003.
  * Lachos et al. [2010] V.H. Lachos, Pulak Ghosh, and R.B. Arellano-Valle.  Likelihood based inference for skew-normal independent linear mixed models.  _Statistica Sinica_ , 20(1):303–322, 2010.
  * Laird and Ware [1982] Nan M. Laird and James H. Ware.  Random-Effects Models for Longitudinal Data.  _Biometrics_ , 38(4):963–974, 1982.
  * Leaker et al. [2016] B.R. Leaker, V.A. Malkov, R. Mogg, M.K. Ruddy, G.C. Nicholson, A.J. Tan, C. Tribouley, and G. Chen.  The nasal mucosal late allergic reaction to grass pollen involves type 2 inflammation ( IL-5 and IL-13 ), the inflammasome ( IL-1 b ), and complement.  _Nature_ , 10(2):408–420, 2016.
  * Massey [1951] Frank J. Massey.  The Kolmogorov-Smirnov Test for Goodness of Fit.  _Journal of the American Statistical Association_ , 46(253):68–78, 1951.
  * Pareek et al. [2021] Savita Pareek, Kalyan Das, and Siuli Mukhopadhyay.  Likelihood-based missing data analysis in multivariate crossover trials.  _arXiv pre-print_ , 2021.
  * Pereira and Russo [2019] Marcos Antonio Alves Pereira and Cibele Maria Russo.  Nonlinear mixed-effects models with scale mixture of skew-normal distributions.  _Journal of Applied Statistics_ , 46(9):1602–1620, 2019.
  * Putt and Chinchilli [1999] Mary Putt and Vernon M. Chinchilli.  A mixed effects model for the analysis of repeated measures cross-over studies.  _Statistics in Medicine_ , 18(22):3037–3058, 1999\.
  * R Core Team [2022] R Core Team.  _R: A Language and Environment for Statistical Computing_.  R Foundation for Statistical Computing, Vienna, Austria, 2022.  URL https://www.R-project.org/.
  * Schumacher et al. [2021] Fernanda L. Schumacher, V.H. Lachos, and Larissa A. Matos.  Scale mixture of skew-normal linear mixed models with within-subject serial dependence.  _Statistics in Medicine_ , 40(7):1790–1810, 2021\.
  * Senn [2002] Stephen Senn.  _Cross-Over Trials in Clinical Research_.  John Wiley & Sons, Ltd., 2002.
  * Shapiro and Wilk [1965] S.S. Shapiro and M.B. Wilk.  An Analysis of Variance Test for Normality (Complete Samples).  _Biometrika_ , 52(3):591–611, 1965.
  * Tudor et al. [2000] Gail E. Tudor, Gary G. Koch, and Diane Catellier.  Statistical methods for crossover designs in bioenvironmental and public health studies.  _Handbook of Statistics_ , 18:571–614, 2000.
  * Verbeke and Molenberghs [2000] Geert Verbeke and Geert Molenberghs.  _Linear Mixed Models for Longitudinal Data_.  Springer-Verlag New York, Inc., 2000.
  * Wu [1983] C.F. Jeff Wu.  On the convergence properties of the EM algorithm.  _Annals of Statistics_ , 11(1):95–103, 1983.
  * Zhang and Davidian [2001] Daowen Zhang and Marie Davidian.  Linear Mixed Models with Flexible Distributions of Random Effects for Longitudinal Data.  _Biometrics_ , 57:795–802, 2001.

##  9 Appendix

###  9.1 Skew-Normal Distribution

The following is a brief overview of the skew-normal distribution and the
terminology that we have used in our analysis.

  * (i)

Univariate skew-normal variate (Azzalini 4): If a random variable $W$ has the
density function

$f(w;\lambda)=2\phi(w)\Phi(\lambda w),\;\;-\infty<w<\infty.$

where $\phi$, $\Phi$ are the standard normal density and distribution
function, respectively, then we say $W$ is a skew-normal variate with
parameter $\lambda$, or, $W\sim SN(0,1,\lambda)$. The random variable $W$ has
additive representation in terms of normal and half normal distribution, i.e.,
if $U_{0},U_{1}$ are independent standard normals, $\delta\in(-1,1)$ and
‘$\overset{d}{=}$’ meaning is ‘distributed as’ then

$W\overset{d}{=}\sqrt{1-\delta^{2}}U_{0}+\delta|U_{1}|,\;\;\delta=\frac{\lambda}{\sqrt{(1+\lambda^{2})}}.$

The parameters $(0,1,\lambda)$ are not the true mean, variance, and skewness
of the random variate $W$. The moments of skew-normal can be found from its
additive representation as

$\mu_{w}=\delta\sqrt{2/\pi},\;\;\sigma^{2}_{w}=1-\delta^{2}(2/\pi),\;\;\gamma_{1}(W)=\frac{4-\pi}{2}\frac{\mu_{w}^{3}}{\sigma_{w}^{3}}.$
(10)

For applied work, we need to introduce location and scale parameters. If
$W_{1}$ is a continuous random variable with location and scale parameters
$\xi,\omega$, respectively then the variable $W_{1}=\xi+\omega W$ will be
skew-normal with parameters $\xi,\omega^{2},\lambda$, with density function

$f_{W_{1}}(w_{1})=2\phi(w_{1}|\xi,\omega^{2})\Phi\left(\lambda\frac{w_{1}-\xi}{\omega}\right),\;-\infty<w_{1}<\infty.$

We use the notation $W_{1}\sim SN(\xi,\omega^{2},\lambda)$. When $\lambda=0$,
it reduces to the normal distribution. Skewed-normal distributions, which are
parametric families governed by moments up to third order, may need a larger
sample size than the traditional rule of thumb of n=30 (Azzalini and Capitanio
6, Chapter 3).

The parameters $(\xi,\omega^{2},\lambda)$ are called direct parameters (DP) as
they appear in the density function and are used to regulate the mean,
variance, and skewness parameters. However, the expected and observed Fisher
information matrix is singular when $\lambda$ is close to 0 (Azzalini and
Capitanio 6). Azzalini 4 proposed that MLE inference be based on centered
parameters (CP) to overcome the singularity problem. CP parameters
$(\mu,\sigma^{2},\gamma)$ represent the actual mean, variance, and skewness
parameters that are derived from an additive representation (eq. 10). CP
parameters $(\mu,\sigma^{2},\gamma)$ are more familiar than DP parameters
$(\xi,\omega^{2},\lambda)$, which makes their interpretation easy. Because CP
parameters are interpretable and close to asymptotic normal distributions,
they are preferred over DP parameters. As a result, CP parameters are suitable
for constructing confidence intervals and other inference methods (Azzalini
and Capitanio 5). More properties of this distribution can be found in
Arellano-Valle and Genton 2.

  * (ii)

Multivariate skew-normal variate (Arellano-Valle et al. 3): An $n$ dimensional
random vector $\bm{W}_{1}$ follows a skew-normal distribution with location
vector $\bm{\mu}\in\mathbb{R}^{n}$, a dispersion matrix $\bm{\Sigma}$ (a
$n\times n$ positive definite matrix) and skewness vector
$\bm{\lambda}\in\mathbb{R}^{n}$, if its pdf is given by

$f_{\bm{W}_{1}}(\bm{w_{1}})=2\phi_{n}(\bm{w}_{1}|\bm{\mu},\bm{\Sigma})\Phi_{1}(\bm{\lambda}^{T}\bm{\Sigma}^{-1/2}(\bm{w}_{1}-\bm{\mu})),\;\;\bm{w}_{1}\in\mathbb{R}^{n}.$

We denote it by $\bm{W}_{1}\sim SN_{n}(\bm{\mu},\bm{\Sigma},\bm{\lambda})$,
further $\bm{W}_{1}$ can be written as
$\bm{W}_{1}{=}\bm{\mu}+\bm{\Sigma}^{1/2}\bm{W}$, where, $\bm{W}$ is a
standardized skew-normal vector. This definition is derived from the
fundamental skew-normal distribution introduced by Azzalini and Dalla Valle 7,
Azzalini and Capitanio 5. The standardised skew-normal, $\bm{W}\sim
SN(\bm{0},\bm{I},\bm{\lambda})$ can also be represented as

$\bm{W}\overset{d}{=}\bm{\delta}|U_{0}|+(\bm{I}_{n}-\bm{\delta}\bm{\delta}^{T})^{1/2}\bm{U}_{1},\;\;\bm{\delta}=\frac{\bm{\lambda}}{\sqrt{1+\bm{\lambda}^{T}\bm{\lambda}}},$

where, $U_{0}\sim N(0,1)$ independent of $\bm{U}_{1}\sim
N_{n}(\bm{0},\bm{I}_{n})$. One important aspect of the multivariate skew
normal is that the joint independence of a random vector holds only if, at
most, one of them is marginally skew-normal (Azzalini and Capitanio 6). In
Arellano-Valle and Azzalini 1, a centered parameterization (CP) is proposed
for multivariate SN distributions.

###  9.2 Some Results for Maximising the Likelihood

We will use the following results from Arellano-Valle et al. 3 to implement
the two steps of EM for maximizing the likelihood.

######  Lemma 9.1.

Let $\textbf{W}\sim{SN}_{n}({\bm{\lambda}})$. Then

$\textbf{W}\overset{d}{=}\bm{\delta}|X_{0}|+(\textbf{I}_{n}-\bm{\delta}\bm{\delta}^{T})^{1/2}\textbf{X}_{1},\;\;\text{where}\;\bm{\delta}=\frac{\bm{\lambda}}{\sqrt{1+\bm{\lambda^{T}}\bm{\lambda}}},$

$X_{0}\sim N_{1}(0,1)$ independent of $\textbf{X}_{1}\sim
N_{n}(\textbf{0},\textbf{I}_{n})$ and ‘$\overset{d}{=}$’ meaning ‘distributed
as’.

######  Lemma 9.2.

Suppose that $\bm{Y}|T=t\sim N_{n}(\bm{\mu}+\bm{d}t,\bm{\Psi})$ and $T\sim
HN_{1}(0,1)$ (the standardized half normal distribution). Let
$\bm{\Sigma}=\bm{\Psi}+\bm{d}\bm{d}^{T}$. Then the joint distribution of
$(\bm{Y}^{t},T)^{T}$ can be written as

$f_{\bm{Y},T}(\bm{y},t|\bm{\theta,\lambda})=2\phi_{n}(\bm{y|\mu,\Sigma})\phi_{1}(t|\eta,\zeta^{2})I\\{t>0\\}.$

The marginal distribution of $\bm{Y}$ after integrating out $t$ is given by

$f_{\bm{Y}}(\bm{y}|\bm{\theta,\lambda})=2\phi_{n}(\bm{y|\mu,\Sigma})\Phi_{1}(\eta|\zeta),$
$\text{where,}\;\eta=\frac{\bm{d}^{T}\bm{\Psi}^{-1}(\bm{y}-\bm{\mu})}{1+\bm{d}^{T}\bm{\Psi}^{-1}\bm{d}}\;\;\text{and}\;\zeta^{2}=\frac{1}{1+\bm{d}^{T}\bm{\Psi}^{-1}\bm{d}}.$

######  Lemma 9.3.

Under the condition in Lemma 9.2,

$E[T^{k}|\bm{y}]=E[X^{k}|X>0],$

where $X\sim N_{1}(\eta,\zeta^{2})$ with $\eta$ and $\zeta^{2}$ given above.
Particularly,

$E[T|\bm{y}]=\eta+\frac{\phi(\frac{\eta}{\zeta})}{\Phi_{1}(\frac{\eta}{\zeta})}\zeta,$
$\text{and}\;E[T^{2}|\bm{y}]=\eta^{2}+\zeta^{2}+\frac{\phi(\frac{\eta}{\zeta})}{\Phi_{1}(\frac{\eta}{\zeta})}\eta\zeta.$

###  9.3 Elements of Hessian Matrix for SN Errors

Following the notation defined in Section 4.1, we proceed as follows to find
the second-order derivative of the Q-function. For $a=1(1)3,b=a(1)3$, we have,

$\displaystyle\frac{\partial^{2}Q}{\partial{\xi}_{[a]}\partial{\xi}_{[b]}}$
$\displaystyle=\sum_{i=1}^{s}\sum_{j=1}^{n_{i}}-\frac{1}{2}[\text{tr}(\bm{V}_{ij{\xi}_{[b]}}^{*}\bm{V}_{ij{\xi}_{[a]}}+\bm{V}_{ij}^{-1}\bm{V}_{ij{\xi}{\xi}_{[a,b]}})$
$\displaystyle+\bm{d}_{ij}^{T}{\bm{V}^{*}_{ij{\xi}\xi_{[a,b]}}}\bm{d}_{ij}T_{ij}^{02}+2\bm{d}_{ij}^{T}{\bm{V}^{*}_{ij{\xi}_{[a]}}}{\bm{d}_{ij{\xi}_{[b]}}}T_{ij}^{02}$
$\displaystyle+2(\bm{d}_{ij}^{T}\bm{V}_{ij}^{-1}{\bm{d}_{ij{\xi}{\xi}_{[a,b]}}}+{\bm{d}_{ij{\xi}_{[b]}}}\bm{V}_{ij}^{-1}{\bm{d}_{ij{\xi}_{[a]}}}+\bm{d}_{ij}^{T}\bm{V}_{ij{\xi}_{[b]}}^{*}{\bm{d}_{ij{\xi}_{[a]}}})T_{ij}^{02}$
$\displaystyle+(\bm{y}_{ij}-\bm{X}_{ij}\bm{\beta})^{T}{\bm{V}^{*}_{ij{\xi}\xi_{[a,b]}}}(\bm{y}_{ij}-\bm{X}_{ij}\bm{\beta}-2\bm{d}_{ij}T_{ij}^{01})$
$\displaystyle-2(\bm{y}_{ij}-\bm{X}_{ij}\bm{\beta})^{T}{\bm{V}^{*}_{ij{\xi}_{[a]}}}{\bm{d}_{ij{\xi}_{[b]}}}T_{ij}^{01}-2(\bm{y}_{ij}-\bm{X}_{ij}\bm{\beta})^{T}{\bm{V}^{*}_{ij{\xi}_{[b]}}}{\bm{d}_{ij{\xi}_{[a]}}}T_{ij}^{01}$
$\displaystyle-2(\bm{y}_{ij}-\bm{X}_{ij}\bm{\beta})^{T}{\bm{V}^{-1}_{ij}}{\bm{d}_{ij{\xi}{\xi}_{[a,b]}}}T_{ij}^{01}],$

where,

$\displaystyle\bm{V}_{ij{\xi}{\xi}}$
$\displaystyle=\frac{\partial\bm{V}_{ij\bm{\xi}}}{\partial\bm{{\xi}}}=(\bm{V}_{ij{\xi}{\xi}_{[a,b]}})_{3\times
3}=\begin{bmatrix}\frac{\partial^{2}\bm{V}_{ij}}{\partial\sigma_{e}^{2}\partial\sigma_{e}^{2}}&\frac{\partial^{2}\bm{V}_{ij}}{\partial\sigma_{e}^{2}\partial\sigma_{s}^{2}}&\frac{\partial^{2}\bm{V}_{ij}}{\partial\sigma_{e}^{2}\partial\lambda_{e}}\\\
&\frac{\partial^{2}\bm{V}_{ij}}{\partial\sigma_{s}^{2}\partial\sigma_{s}^{2}}&\frac{\partial^{2}\bm{V}_{ij}}{\partial\sigma_{s}^{2}\partial\lambda_{e}}\\\
&&\frac{\partial^{2}\bm{V}_{ij}}{\partial\lambda_{e}\partial\lambda_{e}}\end{bmatrix},$
$\displaystyle\frac{\partial^{2}\bm{V}_{ij}}{\partial\sigma_{e}^{2}\partial\lambda_{e}}$
$\displaystyle=\bm{R}_{l},\;\;\frac{\partial^{2}\bm{V}_{ij}}{\partial\lambda_{e}\partial\lambda_{e}}=\sigma_{e}^{2}\bm{R}_{ll},$
$\displaystyle\bm{R}_{ll}$
$\displaystyle=\frac{\partial\bm{R}_{l}}{\partial\lambda}=-2(\bm{\delta}_{e}\bm{\delta}_{ell}^{T}+\bm{\delta}_{el}\bm{\delta}_{el}^{T}),$
$\displaystyle{\bm{V}^{*}_{ij{\xi}\xi_{[a,b]}}}$
$\displaystyle=\frac{\partial^{2}\bm{V}_{ij}^{-1}}{\partial{\xi}_{[a]}\partial{\xi}_{[b]}}=-\bm{V}_{ij}^{-1}\bm{V}_{ij{{\xi}{\xi}}_{[a,b]}}\bm{V}_{ij}^{-1}-2\bm{V}_{ij{\xi}_{[b]}}^{*}\bm{V}_{ij{{\xi}_{[a]}}}\bm{V}_{ij}^{-1}.$

Remaining entries of the matrix $\bm{V}_{ij{\xi}{\xi}}$ are zero.

$\displaystyle\bm{d}_{ij{\xi}{\xi}}$
$\displaystyle=\frac{\partial\bm{d}_{ij\bm{\xi}}}{\partial\bm{{\xi}}}=(\bm{d}_{ij{\xi}{\xi}_{[a,b]}})_{3\times
3}=\begin{bmatrix}\frac{\partial^{2}\bm{d}_{ij}}{\partial\sigma_{e}^{2}\partial\sigma_{e}^{2}}&0&\frac{\partial^{2}\bm{d}_{ij}}{\partial\sigma_{e}^{2}\partial\lambda_{e}}\\\
&0&0\\\
&&\frac{\partial^{2}\bm{d}_{ij}}{\partial\lambda_{e}\partial\lambda_{e}}\end{bmatrix},$
$\displaystyle\frac{\partial^{2}\bm{d}_{ij}}{\partial\sigma_{e}^{2}\partial\sigma_{e}^{2}}$
$\displaystyle=\frac{-\bm{\delta}_{e}}{4{\sigma_{e}}^{3}},\;\;\frac{\partial^{2}\bm{d}_{ij}}{\partial\sigma_{e}^{2}\partial\lambda_{e}}=\frac{\bm{\delta}_{el}}{2\sqrt{\sigma_{e}^{2}}},$
$\displaystyle\frac{\partial^{2}\bm{d}_{ij}}{\partial\lambda_{e}\partial\lambda_{e}}$
$\displaystyle=\sigma_{e}\bm{\delta}_{ell},\;\;\bm{\delta}_{ell}=\frac{\partial\bm{\delta}_{el}}{\partial\lambda_{e}}=\frac{-3\lambda_{e}(1,0,\ldots,0)}{(1+\lambda_{e}^{2})^{5/2}}.$

###  9.4 Elements of Hessian Matrix for SN Random Effect

Based on the notation defined in Section 4.2, we proceed as follows to find
the second order derivative of the Q-function. For $a=1(1)3,b=a(1)3$, we have,

$\displaystyle\frac{\partial^{2}Q}{\partial{\xi}_{[a]}\partial{\xi}_{[b]}}$
$\displaystyle=\sum_{i=1}^{s}\sum_{j=1}^{n_{i}}-\frac{1}{2}[\text{tr}(\bm{V}_{ij{\xi}_{[b]}}^{*}\bm{V}_{ij{\xi}_{[a]}}+\bm{V}_{ij}^{-1}\bm{V}_{ij{\xi}{\xi}_{[a,b]}})$
$\displaystyle+\bm{d}_{ij}^{T}{\bm{V}^{*}_{ij{\xi}\xi_{[a,b]}}}\bm{d}_{ij}T_{ij}^{02}+2\bm{d}_{ij}^{T}{\bm{V}^{*}_{ij{\xi}_{[a]}}}{\bm{d}_{ij{\xi}_{[b]}}}T_{ij}^{02}$
$\displaystyle+2(\bm{d}_{ij}^{T}\bm{V}_{ij}^{-1}{\bm{d}_{ij{\xi}{\xi}_{[a,b]}}}+{\bm{d}_{ij{\xi}_{[b]}}}\bm{V}_{ij}^{-1}{\bm{d}_{ij{\xi}_{[a]}}}+\bm{d}_{ij}^{T}\bm{V}_{ij{\xi}_{{[b]}}}^{*}{\bm{d}_{ij{\xi}_{[a]}}})T_{ij}^{02}$
$\displaystyle+(\bm{y}_{ij}-\bm{X}_{ij}\bm{\beta})^{T}{\bm{V}^{*}_{ij{\xi}\xi_{[a,b]}}}(\bm{y}_{ij}-\bm{X}_{ij}\bm{\beta}-2\bm{d}_{ij}T_{ij}^{01})$
$\displaystyle-2(\bm{y}_{ij}-\bm{X}_{ij}\bm{\beta})^{T}{\bm{V}^{*}_{ij{\xi}_{[a]}}}{\bm{d}_{ij{\xi}_{[b]}}}T_{ij}^{01}-2(\bm{y}_{ij}-\bm{X}_{ij}\bm{\beta})^{T}{\bm{V}^{*}_{ij{\xi}_{[b]}}}{\bm{d}_{ij{\xi}_{[a]}}}T_{ij}^{01}$
$\displaystyle-2(\bm{y}_{ij}-\bm{X}_{ij}\bm{\beta})^{T}{\bm{V}^{-1}_{ij}}{\bm{d}_{ij{\xi}{\xi}_{[a,b]}}}T_{ij}^{01}],$

where,

$\displaystyle\bm{V}_{ij{\xi}{\xi}}$
$\displaystyle=\frac{\partial\bm{V}_{ij\bm{\xi}}}{\partial\bm{{\xi}}}=(\bm{V}_{ij{\xi}{\xi}_{[a,b]}})_{3\times
3}=\begin{bmatrix}\frac{\partial^{2}\bm{V}_{ij}}{\partial\sigma_{e}^{2}\partial\sigma_{e}^{2}}&\frac{\partial^{2}\bm{V}_{ij}}{\partial\sigma_{e}^{2}\partial\sigma_{s}^{2}}&\frac{\partial^{2}\bm{V}_{ij}}{\partial\sigma_{e}^{2}\partial\lambda_{s}}\\\
&\frac{\partial^{2}\bm{V}_{ij}}{\partial\sigma_{s}^{2}\partial\sigma_{s}^{2}}&\frac{\partial^{2}\bm{V}_{ij}}{\partial\sigma_{s}^{2}\partial\lambda_{s}}\\\
&&\frac{\partial^{2}\bm{V}_{ij}}{\partial\lambda_{s}\partial\lambda_{s}}\end{bmatrix},$
$\displaystyle\frac{\partial^{2}\bm{V}_{ij}}{\partial\sigma_{s}^{2}\partial\lambda_{s}}$
$\displaystyle=R_{l}\bm{Z}_{ij}\bm{Z}_{ij}^{T},\;\;\frac{\partial^{2}\bm{V}_{ij}}{\partial\lambda_{s}\partial\lambda_{s}}=\sigma_{s}^{2}R_{ll}\bm{Z}_{ij}\bm{Z}_{ij}^{T},$
$\displaystyle R_{ll}$ $\displaystyle=\frac{\partial
R_{l}}{\partial\lambda_{s}}=-2(\delta_{b}\delta_{bll}+\delta_{bl}^{2}),$
$\displaystyle{\bm{V}^{*}_{ij{\xi}\xi_{[a,b]}}}$
$\displaystyle=\frac{\partial^{2}\bm{V}_{ij}^{-1}}{\partial{\xi}_{[a]}\partial{\xi}_{[b]}}=-\bm{V}_{ij}^{-1}\bm{V}_{ij{{\xi}{\xi}}_{[a,b]}}\bm{V}_{ij}^{-1}-2\bm{V}_{ij{\xi}_{[b]}}^{*}\bm{V}_{ij{{\xi}_{[a]}}}\bm{V}_{ij}^{-1}.$

Remaining entries of the matrix $\bm{V}_{ij{\xi}{\xi}}$ are zero.

$\displaystyle\bm{d}_{ij{\xi}{\xi}}$
$\displaystyle=\frac{\partial\bm{d}_{ij\bm{\xi}}}{\partial\bm{{\xi}}}=(\bm{d}_{ij{\xi}{\xi}_{[a,b]}})_{3\times
3}=\begin{bmatrix}0&0&0\\\
&\frac{\partial^{2}\bm{d}_{ij}}{\partial\sigma_{s}^{2}\partial\sigma_{s}^{2}}&\frac{\partial^{2}\bm{d}_{ij}}{\partial\sigma_{s}^{2}\partial\lambda_{s}}\\\
&&\frac{\partial^{2}\bm{d}_{ij}}{\partial\lambda_{s}\partial\lambda_{s}}\end{bmatrix},$
$\displaystyle\frac{\partial^{2}\bm{d}_{ij}}{\partial\sigma_{s}^{2}\partial\sigma_{s}^{2}}$
$\displaystyle=\frac{-\bm{Z}_{ij}\delta_{b}}{4{\sigma_{s}}^{3}},\;\;\frac{\partial^{2}\bm{d}_{ij}}{\partial\sigma_{s}^{2}\partial\lambda_{s}}=\frac{\bm{Z}_{ij}\delta_{bl}}{2\sqrt{\sigma_{s}^{2}}},$
$\displaystyle\frac{\partial^{2}\bm{d}_{ij}}{\partial\lambda_{s}\partial\lambda_{s}}$
$\displaystyle=\bm{Z}_{ij}\sigma_{s}\delta_{bll},\;\;\delta_{bll}=\frac{\partial\delta_{bl}}{\partial\lambda_{s}}=\frac{-3\lambda_{s}}{(1+\lambda_{s}^{2})^{5/2}}.$


See pages 1-last of paper-icprs-2017.pdf