Likelihood based missing data analysis in multivariate crossover trials

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

For gene expression data measured in a crossover trial, a multivariate mixed
effects model seems to be most appropriate. Standard statistical inference
fails to provide reliable results when some responses are missing.
Particularly for crossover studies, missingness is a serious concern as the
trial requires a small number of participants. A Monte Carlo EM (MCEM) based
technique has been adopted to deal with this situation. Along with estimation,
a MCEM likelihood ratio test (LRTs) is developed for testing the fixed effects
in such a multivariate crossover model with missing data. Intensive simulation
studies have been carried out prior to the analysis of the gene expression
data.

Keywords: Multivariate crossover trials; Monte Carlo EM algorithm; MCEM LRTs;

##  1 Introduction

In many clinical studies, we may come across crossover trials with two or more
response variates with subject dropouts and incomplete data. For example, one
may consider the measurement of both systolic (SBP)and diastolic (DBP) blood
pressure in each period (Grender and Johnson 6) or recording blood sugar
levels at multiple time points in each period (Putt and Chinchilli 14) or
microarray gene expression profiles of subjects measured in each period
(Leaker et al. 11) with incomplete data. Using conventional univariate
crossover models based on only complete cases to analyse such missing and
multivariate data, may lead to incorrect and biased parameter and variance
estimates. In a crossover trial, usually a smaller number of participants as
compared to cross-sectional studies are recruited as the same subjects are
switched between all treatment groups. Thus, ignoring subject dropout
information and performing a complete case analysis may lead to highly
underpowered studies and biased estimates. Similarly, ignoring the correlated
nature of the multiple responses recorded in each period and fitting separate
crossover models to each response may lead us to inaccurate conclusions about
the intra and inter subject variabilities. Thus, to address the issue of
missing responses in a multivariate crossover setup, we propose a mixed effect
model approach and an Monte Carlo EM (MCEM) based estimation method.

For a detailed review of univariate crossover designs, refer to the books by
(Senn 17; Jones and Kenward 9). In comparison, crossover trials with multiple
responses measured in each period have been addressed by fewer researchers,
namely (Grender and Johnson 6; Chinchilli and Esinhart 2; Putt and Chinchilli
14; Tudor et al. 19). To tackle the presence of missing/incomplete data
specifically for the $AB/BA$ design under the MAR or MNAR missing mechanism
one may refer to (Patel 13; Ho et al. 7; Matthews and Henderson 12). For more
general univariate crossover designs with incomplete data references are
namely, (Richardson and Flack 15; Chow and Shao 3; Basu and Santra 1;
Rosenkranz 16). However, none of the works cited deal with missingness in a
multivariate crossover design setup.

In this article we model and analyse multivariate crossover studies in the
presence of missing at random outcomes. This work is motivated by a $3\times
3$ crossover trial studying the effect of two single doses of oral prednisone
(10 and 25 mg) with placebo on biomarkers of mucosal inflammation and
transcriptomics (Leaker et al. 11). The gene expression changes of subjects
were measured in each period of the trial giving rise to multiple responses in
a crossover framework. A total of seventeen subjects assigned to three
treatment sequences were enrolled in the crossover study, however some missing
observations (particularly) in the third period were observed. Our task here
is to fit an appropriate statistical model to the crossover trial with
microarray expression data taking into account the multivariate structure of
the data and the missing responses. A Monte Carlo EM based mixed model
analysis is used to fit and estimate the treatment and gene effects along with
the intra and inter subject variability in the presence of missing at random
outcomes. The mixed model assumed, includes the fixed effects of treatment,
period and genes, and their interactions while the random effects are subject
and subject gene interaction effect. We used maximum likelihood estimation
coupled with the MCEM algorithm to determine the estimated model parameters
and variance components. For finding the asymptotic covariance matrix,
multiple imputation method of (Els and Ryan 5) was considered. Likelihood
ratio tests (LRTs) based on the MCEM estimators were formulated. We used
detailed simulations to study the properties of the estimators and power of
the LRTs proposed. A real data example based on the gene expression levels
measured in a $3\times 3$ crossover study is used to illustrate the
statistical model and its estimation.

To the best of our knowledge this is the first work involving multivariate
crossover studies with missing responses. Though originally developed for
agricultural sciences, crossover trials are nowadays frequently used in
clinical trials and biological studies to evaluate treatment effects. These
trials are most useful for comparing treatments for chronic or long term
diseases, such as asthma, rheumatism, hypertension etc. Note that one useful
advantage of crossover studies is that they require a smaller number of
participants. However, due to smaller number of recruits, reliable statistical
inference is harder to make when some responses are randomly missing.

##  2 Case study: multivariate crossover trial of oral prednisone

We use a gene expression dataset from (Leaker et al. 11) as a case study. The
dataset is publicly available from the NCBI Gene Expression Omnibus (Clough
and Barrett 4). 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 also a positive
result from the intraepidermal skin prick test to grass pollen extract. There
are a total of seventeen subjects enrolled in the study, assigned to three
treatment sequences/groups. Some missing observations (particularly) in the
third period are also recorded. The main interest here is to study the effect
of treatments and genes on allergic reaction to grass pollen.

The crossover trial consists of 3 distinct sequences of treatments assigned to
a total of 17 subjects. The sequences of treatments are: Sequence 1 (10 mg,
placebo, 25 mg), Sequence 2 (25 mg, 10 mg, placebo), Sequence 3 (placebo, 25
mg, 10 mg). Subject 1 to subject 6 are in sequence 1, subject 7 to subject 12
are in sequence 2 and subject 13 to subject 17 are in sequence 3. For each
subject, 10 gene expressions are measured by giving treatments (10 mg
prednisone, 25 mg prednisone and placebo) in three periods with a wash out
period of 4 weeks. As an example, Sequence 1 indicates that 10 mg prednisone
is given to the subjects (subject 1 to subject 6) at period 1 followed by a
washout, period 2 placebo is given and again after a washout, 25 mg prednisone
is given at period 3. The study comprises of missing responses in period 3 of
various sequences. In sequence 1, response values of subjects 2, 5 are
missing; in sequence 2, response values of subjects 8, 10 are missing and in
sequence 3, response values from subject 14 are missing. The overall
percentage of missing responses in the data is $9\%$, with $11\%$ missing in
sequences 1 and 2, and $6\%$ in sequence 3.

Prior to model fitting we ran some exploratory analysis of the gene data using
complete cases. The following facts were revealed:

  * (i)

Gene by treatment and gene by period interactions were seen to be almost
absent for each treatment sequence, see fig. 1.

  * (ii)

From fig. 2 it is fair to assume no subject by gene interaction effect exists
in our model.

Figure 1: Period by gene and treatment by gene interaction plots; Plots on the
left (A, B, C) is period by gene interaction plots for sequences 1, 2 and 3
respectively. Plots on the right (D, E, F) is treatment by gene interaction
plots for sequences 1, 2 and 3 respectively. Figure 2: Subject by gene
interaction plot; Plots A, B, C is subject by gene interaction plots for
sequence 1, 2 and 3 respectively

Considering these two points, in the following section, we frame a specific
model that pertains to the above microarray expression data.

##  3 Model and notation

Suppose $y_{ijkl}$ denotes measurement of the $l^{\text{th}}$ gene expression
from the $k^{\text{th}}$ subject in the $j^{\text{th}}$ period assigned to the
$i^{\text{th}}$ treatment sequence, where $i=1,2,\ldots,s$; $j=1,2,\ldots,p$;
$k=1,2,\ldots,n_{i}$ and $l=1,2,\ldots,m$. So there are in total $n$ subjects,
where $n=\sum_{i=1}^{s}n_{i}$. Assuming a normal random effects model we may
write,

$y_{ijkl}=\mu+\pi_{j}+\tau_{d[i,j]}+g_{l}+s_{ik}+e_{ijkl}$ (1)

where $\mu$ is an intercept, $\pi_{j}$ is the fixed effect of the $j$th
period, $\tau_{d[i,j]},$ $(d[i,j]=1,2,\ldots,t)$ is the fixed effect
associated with the treatment applied in period $j$ of the sequence $i$,
$g_{l}$ is the fixed effect of the $l^{\text{th}}$ gene, $s_{ik}$ is a random
effect associated with the $k^{\text{th}}$ subject in sequence $i$, and
$e_{ijkl}$ is a random error. We assume that,
$e_{ijkl}\overset{ind}{\sim}{N}(0,\sigma_{e}^{2})\;\;\forall\,i,j,k,l$, and
$s_{ik}\overset{ind}{\sim}{N}(0,\sigma_{s}^{2})\;\;\forall\,i,k$. Also,
$s_{ik},e_{ijkl}$ are independent. Note, we have not considered an interaction
effect of gene with period and treatment in the above model, in accordance
with the exploratory analysis results for the case study data. Based on the
above assumptions, we have
$\text{var}(y_{ijkl})=\sigma_{e}^{2}+\sigma_{s}^{2}$, and
$\text{cov}(y_{ijkl},y_{ij^{{}^{\prime}}kl^{{}^{\prime}}})=\sigma_{s}^{2}$ for
all $j\neq j^{{}^{\prime}}$, $l\neq l^{{}^{\prime}}$. Thus, variance-
covariance matrix corresponding to the $k^{\text{th}}$ subject administered
the $i^{\text{th}}$ treatment sequence is given by,

$var(\bm{y}_{ik})=(\bm{\Sigma}_{ik})_{pm\times
pm}=\begin{bmatrix}\sigma_{e}^{2}\bm{R}&\sigma_{s}^{2}\bm{1}_{m}\bm{1}_{m}^{T}&.&.&.&\sigma_{s}^{2}\bm{1}_{m}\bm{1}_{m}^{T}\\\
\sigma_{s}^{2}\bm{1}_{m}\bm{1}_{m}^{T}&\sigma_{e}^{2}\bm{R}&.&.&.&\sigma_{s}^{2}\bm{1}_{m}\bm{1}_{m}^{T}\\\
\vdots\\\
\sigma_{s}^{2}\bm{1}_{m}\bm{1}_{m}^{T}&\sigma_{s}^{2}\bm{1}_{m}\bm{1}_{m}^{T}&.&.&.&\sigma_{e}^{2}\bm{R}\end{bmatrix}$

where,
$\bm{R}=\bm{I}_{m}+\frac{\sigma_{s}^{2}}{\sigma_{e}^{2}}\bm{1}_{m}\bm{1}_{m}^{T}$
and $\bm{1}_{m}^{T}=(1,1,\ldots,1)_{1\times m}$.

Using matrix notations we rewrite equation 1 as

$\bm{y}_{i}=\bm{X}_{i}\bm{\beta}+\bm{Z}_{i}\bm{b}_{i}+\bm{e}_{i},~{}i=1,2,\ldots,s$
(2)

where,

  * (i)

$\bm{y}_{i}=(\bm{y}_{i1}^{T},\bm{y}_{i2}^{T},\ldots,\bm{y}_{iq}^{T},\dots,\bm{y}_{in_{i}}^{T})^{T}$
is a response vector of length $pmn_{i}$, and
$\bm{y}_{iq}=(y_{i1q1},\ldots,y_{ipqm})_{pm\times 1}^{T}$.

  * (ii)

$\bm{X}_{i}=(\bm{X}_{i1}^{T},\bm{X}_{i2}^{T},\ldots,\bm{X}_{iq}^{T},\ldots,\bm{X}_{in_{i}}^{T})^{T}_{pmn_{i}\times(p+t+m-2)}$
is the design matrix corresponding to the fixed effects, and

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

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

  * (iii)

$\bm{Z}_{i}=\bm{I}_{n_{i}}\otimes\bm{1}_{pm}$ is the design matrix
corresponding to the random effects.

  * (iv)

Moreover, $\bm{b}_{i}=(s_{i1},s_{i2},\ldots,s_{in_{i}})^{T}$ is the
${n_{i}\times 1}$ vector of random effects, and
$\bm{e}_{i}=(\bm{e}_{i1}^{T},\bm{e}_{i2}^{T},\ldots,\bm{e}_{iq}^{T},\ldots,\bm{e}_{in_{i}}^{T})^{T}$
is the ${pmn_{i}\times 1}$ random error vector, where
$\bm{e}_{iq}=(e_{i1q1},\ldots,e_{ipqm})^{T}$, it is assumed

$\bm{e}_{i}\sim{N}(\bm{0},\sigma_{e}^{2}\bm{I}_{pmn_{i}});\;\;\bm{b}_{i}\sim{N}_{n_{i}}(\bm{0},\bm{D}_{i}),\;\bm{D}_{i}=\sigma_{s}^{2}\bm{I}_{n_{i}}.$

  * (v)

The parameters of interest are,
$\bm{\beta}=(\mu,\pi_{1},\dots\pi_{p-1},\tau_{1},\ldots,\tau_{t-1},g_{1},\ldots,g_{m-1})$
and $\sigma^{2}=(\sigma_{e}^{2},\sigma_{s}^{2}).$

Under the above assumptions, the conditional distribution of $\bm{y}_{i}$
given the random effects $\bm{b}_{i}$ is of the form

$\bm{y}_{i}|\bm{\beta},\sigma_{e}^{2},\bm{b}_{i}\sim{N}_{pmn_{i}}(\bm{X}_{i}\bm{\beta}+\bm{Z}_{i}\bm{b}_{i},\sigma_{e}^{2}\bm{I}_{pmn_{i}}),$

while the marginal distribution of $\bm{y}_{i}$ is

$\bm{y}_{i}|\bm{\beta},\sigma_{e}^{2},\bm{D}_{i}\sim{N}_{pmn_{i}}(\bm{X}_{i}\bm{\beta},\bm{\Sigma}_{i}),\;\bm{\Sigma}_{i}=\bm{Z}_{i}\bm{D}_{i}\bm{Z}_{i}^{T}+\sigma_{e}^{2}\bm{I}_{pmn_{i}}$

##  4 Estimation with missing at random responses

In this section we illustrate how an EM type algorithm can be used for maximum
likelihood estimation (MLE) in multivariate crossover model in a missing set
up. We have used the MCEM algorithm of (Wei and Tanner 20) for parameter
estimation. Note that for ignorable missing (MAR) assumption in the data, it
is not necessary to incorporate the parametric model for the missing data
mechanism into the complete data log likelihood. Therefore the complete data
log likelihood for the $i$th sequence is given by $f(\bm{y}_{i},\bm{b}_{i})$,
where $\bm{y}_{i}$ may contain missing values. We assume the covariates to be
fully observed.

In the E-step we calculate the expected value of the complete data log
likelihood given the observed data and current parameter estimates. Suppose,
$\bm{y}_{i}=(\bm{y}_{mis,i}^{T},\bm{y}_{obs,i}^{T})^{T}$ where
$\bm{y}_{mis,i}$ is the $m_{i}\times 1$ vector of missing components of
$\bm{y}_{i}$. Note, both $\bm{b}_{i}$ and $\bm{y}_{mis,i}$ are unobserved, and
must be integrated over. To execute the EM algorithm the steps listed were
followed from (Ibrahim and Molenberghs 8). Thus the E-step for $i$-th sequence
at $(t+1)$th iteration.

$\displaystyle Q_{i}(\bm{\gamma}|\bm{\gamma}^{(t)})$
$\displaystyle=E\left(l(\bm{\gamma};\bm{y}_{i},\bm{b}_{i}|\bm{y}_{obs,i},\bm{\gamma}^{(t)})\right)$
$\displaystyle=\int\int\log\left[f(\bm{y}_{i}|\bm{\beta},\sigma_{e}^{2},\bm{b}_{i})\right]f(\bm{y}_{mis,i},\bm{b}_{i}|\bm{y}_{obs,i},\bm{\gamma}^{(t)})d\bm{b}_{i}d\bm{y}_{mis,i}$
$\displaystyle+\int\int\log\left[f(\bm{b}_{i}|D_{i})\right]f(\bm{y}_{mis,i},\bm{b}_{i}|\bm{y}_{obs,i},\bm{\gamma}^{(t)})d\bm{b}_{i}d\bm{y}_{mis,i}$
$\displaystyle\equiv I_{1}+I_{2},$ (3)

where
$\bm{\gamma}^{(t)}=(\bm{\beta}^{(t)},\sigma_{e}^{2(t)},\sigma_{s}^{2(t)})$.

To integrate out $\bm{b}_{i}$ from $I_{1}$ and $I_{2}$ we write

$\displaystyle f(\bm{y}_{mis,i},\bm{b}_{i}|\bm{y}_{obs,i},\bm{\gamma}^{(t)})$
$\displaystyle=f(\bm{b}_{i}|\bm{y}_{i},\bm{\gamma}^{(t)})f(\bm{y}_{mis,i}|\bm{y}_{obs,i},\bm{\gamma}^{(t)})$
$\text{where,}\;\;\;\bm{b}_{i}|\bm{y}_{i},\bm{\gamma}^{(t)}\sim
N(\bm{b}_{0i}^{(t)},\bm{\Sigma}_{0i}^{(t)})\;\;\;\text{and,}$
$\displaystyle\bm{b}_{0i}^{(t)}$
$\displaystyle=\bm{D}_{i}^{(t)}\bm{Z}_{i}^{{}^{\prime}}(\bm{Z}_{i}\bm{D}_{i}^{(t)}\bm{Z}_{i}^{{}^{\prime}}+\sigma_{e}^{2}\bm{I}_{pmn_{i}})^{-1}(\bm{y}_{i}-\bm{X}_{i}\bm{\beta}^{(t)})$
$\displaystyle=\bm{\Sigma}_{0i}^{(t)}\bm{Z}_{i}^{{}^{\prime}}(\bm{y}_{i}-\bm{X}_{i}\bm{\beta}^{(t)})/\sigma_{e}^{2(t)}$
$\displaystyle\bm{\Sigma}_{0i}^{(t)}$
$\displaystyle=\bm{D}_{i}^{(t)}-\bm{D}_{i}^{(t)}\bm{Z}_{i}^{{}^{\prime}}(\bm{Z}_{i}\bm{D}_{i}^{(t)}\bm{Z}_{i}^{{}^{\prime}}+\sigma_{e}^{2}\bm{I}_{pmn_{i}})^{-1}\bm{Z}_{i}\bm{D}_{i}^{(t)}$
$\displaystyle=\left[\sigma_{e}^{-2(t)}\bm{Z}_{i}^{{}^{\prime}}\bm{Z}_{i}+(\bm{D}_{i}^{(t)})^{-1}\right]^{-1}.$

Thus,

$\displaystyle I_{1}$ $\displaystyle=\frac{-pmn_{i}}{2}\log
2\pi-\frac{pmn_{i}}{2}\log\sigma_{e}^{2}-\frac{1}{2\sigma_{e}^{2}}(Tr(\bm{Z}_{i}^{T}\bm{Z}_{i}\bm{\Sigma}_{0i}^{(t)})$
$\displaystyle+\int(\bm{y}_{i}-\bm{X}_{i}\bm{\beta}-\bm{Z}_{i}\bm{b}_{0i}^{(t)})^{T}(\bm{y}_{i}-\bm{X}_{i}\bm{\beta}-\bm{Z}_{i}\bm{b}_{0i}^{(t)})f(\bm{y}_{mis,i}|\bm{y}_{obs,i},\bm{\gamma}^{(t)})d\bm{y}_{mis,i}).$
(4)

and,

$\displaystyle I_{2}$ $\displaystyle=\frac{-(m+1)n_{i}}{2}\log
2\pi-\frac{1}{2}\log(\det\bm{D}_{i})-\frac{1}{2}Tr(\bm{D}_{i}^{-1}\bm{\Sigma}_{0i}^{(t)})$
$\displaystyle-\frac{1}{2}\int({\bm{b}_{0i}^{(t)}}^{T}\bm{D}_{i}^{-1}\bm{b}_{0i}^{(t)})f(\bm{y}_{mis,i}|\bm{y}_{obs,i},\bm{\gamma}^{(t)})d\bm{y}_{mis,i}.$
(5)

We sample $v_{i1},v_{i2},\ldots v_{ic_{i}}$ from
$[\bm{y}_{mis,i}|\bm{y}_{obs,i},\bm{\gamma}^{(t)}]$, which has the density of
the form,

$\bm{y}_{mis,i}|\bm{y}_{obs,i},\bm{\gamma}^{(t)}\sim{N}_{m_{i}}(\bm{\mu}_{1}+\bm{\Sigma}_{12}\bm{\Sigma}_{22}^{-1}(\bm{y}_{obs,i}-\bm{\mu}_{2}),\bm{\Sigma}_{11}-\bm{S}_{12}\bm{\Sigma}_{22}^{-1}\bm{\Sigma}_{21}).$
(6)

where,

$\begin{pmatrix}\bm{y}_{mis,i}\\\
\bm{y}_{obs,i}\end{pmatrix}=\begin{pmatrix}\bm{X}_{i_{m_{i}\times 1}}^{(1)}\\\
\bm{X}_{i_{l_{i}\times
1}}^{(2)}\end{pmatrix}\bm{\beta}+\begin{pmatrix}\bm{Z}_{i_{m_{i}\times
1}}^{(1)}\\\ \bm{Z}_{i_{l_{i}\times
1}}^{(2)}\end{pmatrix}\bm{b_{i}}+\begin{pmatrix}\bm{e}_{m_{i}}\\\
\bm{e}_{l_{i}}\end{pmatrix};\;\;m_{i}+l_{i}=pmn_{i}$
$\bm{y}_{i}=\begin{pmatrix}\bm{y}_{mis,i}\\\
\bm{y}_{obs,i}\end{pmatrix}\sim{N}(\bm{\mu},\bm{\Sigma}),$
$\bm{\mu}=\begin{pmatrix}\bm{\mu}_{1}\\\
\bm{\mu}_{2}\end{pmatrix}=\begin{pmatrix}\bm{X}_{i}^{(1)}\bm{\beta}\\\
\bm{X}_{i}^{(2)}\bm{\beta}\end{pmatrix},$
$\bm{\Sigma}=\begin{pmatrix}\bm{\Sigma}_{11}&\bm{\Sigma}_{12}\\\
\bm{\Sigma}_{21}&\bm{\Sigma}_{22}\end{pmatrix}=\begin{pmatrix}\bm{Z}_{i}^{(1)}\bm{D}_{i}{\bm{Z}_{i}^{(1)}}^{T}&\bm{Z}_{i}^{(1)}\bm{D}_{i}{\bm{Z}_{i}^{(2)}}^{T}\\\
\bm{Z}_{i}^{(2)}\bm{D}_{i}{\bm{Z}_{i}^{(1)}}^{T}&\bm{Z}_{i}^{(2)}\bm{D}_{i}{\bm{Z}_{i}^{(2)}}^{T}\end{pmatrix}+\sigma_{e}^{2}\begin{pmatrix}\bm{I}_{m_{i}}&0\\\
0&\bm{I}_{l_{i}}\end{pmatrix}.$

Using the data $\bm{y}_{i}^{(k)}=(\bm{v}_{ik}^{T},\bm{y}_{obs,i}^{T})^{T}$ and
$\bm{b}_{0i}^{(tk)}=\bm{\Sigma}_{0i}^{(t)}\bm{Z}_{i}^{T}(\bm{y}_{i}^{(k)}-\bm{X}_{i}\bm{\beta}^{(t)})\sigma_{e}^{2(t)}$,
eq. 4, eq. 5 and eq. 6 the $E$-step for the $i^{\text{th}}$ sequence (eq. 3)
in the $(t+1)^{\text{th}}$ iteration takes the form

$\displaystyle Q_{i}(\bm{\gamma}|\bm{\gamma}^{(t)})$
$\displaystyle=\frac{-pmn_{i}}{2}\log
2\pi-\frac{pmn_{i}}{2}\log\sigma_{e}^{2}$
$\displaystyle-\frac{1}{2\sigma_{e}^{2}}\bigg{(}Tr(\bm{Z}_{i}^{T}\bm{Z}_{i}\bm{\Sigma}_{0i}^{(t)})+\frac{1}{c_{i}}\displaystyle\sum_{k=1}^{c_{i}}(\bm{y}_{i}^{(k)}-\bm{X}_{i}\bm{\beta}-\bm{Z}_{i}\bm{b}_{0i}^{(tk)})^{T}(\bm{y}_{i}^{(k)}-\bm{X}_{i}\bm{\beta}-\bm{Z}_{i}\bm{b}_{0i}^{(tk)})\bigg{)}$
$\displaystyle-\frac{(m+1)n_{i}}{2}\log
2\pi-\frac{1}{2}\log(\det\bm{D}_{i})-\frac{1}{2}Tr(\bm{D}_{i}^{-1}\bm{\Sigma}_{0i}^{(t)})$
$\displaystyle-\frac{1}{2c_{i}}\displaystyle\sum_{k=1}^{c_{i}}{\bm{b}_{0i}^{(tk)}}^{T}\bm{D}_{i}^{-1}\bm{b}_{0i}^{(tk)}.$

Obviously, the E-step for all the sequences is given by

$Q(\bm{\gamma}|\bm{\gamma}^{(t)})=\displaystyle\sum_{i=1}^{s}Q_{i}(\bm{\gamma}|\bm{\gamma}^{(t)}).$

In the M-step, we maximize $Q(\bm{\gamma}|\bm{\gamma}^{(t)})$. The closed
forms are obtained for $(\bm{\beta},\sigma_{e}^{2},\sigma_{s}^{2})$ by
computing

$\frac{\partial Q}{\partial\bm{\beta}}=0,~{}~{}~{}\frac{\partial
Q}{\partial\sigma_{e}^{2}}=0,~{}~{}~{}\frac{\partial
Q}{\partial\sigma_{s}^{2}}=0$ $\displaystyle\bm{\beta}^{(t+1)}$
$\displaystyle=\bigg{(}\displaystyle\sum_{i=1}^{s}\bm{X}_{i}^{T}\bm{X}_{i}\bigg{)}^{-1}\displaystyle\sum_{i=1}^{s}\bm{X}_{i}^{T}\frac{1}{c_{i}}\displaystyle\sum_{k=1}^{c_{i}}(\bm{y}_{i}^{(k)}-\bm{Z}_{i}\bm{b}_{0i}^{(tk)})$
$\displaystyle\sigma_{e}^{2(t+1)}$
$\displaystyle=\frac{1}{pmn}\displaystyle\sum_{i=1}^{s}\bigg{(}\frac{1}{c_{i}}\displaystyle\sum_{k=1}^{c_{i}}(\bm{y}_{i}^{(k)}-\bm{X}_{i}\bm{\beta}^{(t+1)}-\bm{Z}_{i}\bm{b}_{0i}^{(tk)})^{T}(\bm{y}_{i}^{(k)}-\bm{X}_{i}\bm{\beta}^{(t+1)}-\bm{Z}_{i}\bm{b}_{0i}^{(tk)})$
$\displaystyle+Tr(\bm{Z}_{i}^{T}\bm{Z}_{i}\bm{\Sigma}_{0i}^{(t)})\bigg{)}$
$\displaystyle\sigma_{s}^{2(t+1)}$
$\displaystyle=\frac{1}{n}\displaystyle\sum_{i=1}^{s}\bigg{(}Tr(\bm{\Sigma}_{0i}^{(t)})+\frac{1}{c_{i}}\displaystyle\sum_{k=1}^{c_{i}}{\bm{b}_{0i}^{(tk)}}^{T}\bm{I}_{n_{i}}\bm{b}_{0i}^{(tk)}\bigg{)}$

The above $E$ and $M-$ steps are repeated until convergence is achieved.

To find the asymptotic covariance matrix, the multiple imputation technique
has been considered. For imputing a missing value we sampled from
$[\bm{b}_{i},\bm{y}_{mis,i}|\hat{\bm{\gamma}},\bm{y}_{obs,i}]$ following the
steps below, note
$\hat{\bm{\gamma}}=(\hat{\bm{\beta}},\hat{\sigma_{e}}^{2},\hat{\sigma_{s}}^{2})$
are the MCEM estimates,

  1. (i)

Initial estimates $({\bm{b}_{i}}_{0},{\bm{y}_{mis,i}}_{0})$ were taken from
the respective marginal distributions i.e. $\bm{b}_{i}\sim
N(\bm{0},\sigma_{s}^{2}\bm{I}_{n_{i}})$ and
$\bm{y}_{mis,i}$ $\sim
N(\bm{X}_{i}^{(1)}\bm{\beta},\bm{Z}_{i}^{(1)}\bm{D}_{i}{\bm{Z}_{i}^{(1)}}^{T}+\sigma_{e}^{2}\bm{I}_{m_{i}})$.

  2. (ii)

At the $k(\geq 1)$th iteration $({\bm{b}_{i}}_{k},{\bm{y}_{mis,i}}_{k})$ were
sampled using ${\bm{b}_{i}}_{k}$ from
$\bm{b}_{i}|\bm{y}_{obs,i},{\bm{y}_{mis,i}}_{k-1}$ and ${\bm{y}_{mis,i}}_{k}$
from ${\bm{y}_{mis,i}}|\bm{y}_{obs,i},{\bm{b}_{i}}_{k}$.

We generated $b$ such samples of $[\bm{b}_{i},\bm{y}_{mis,i}]$. The parameter
estimates and their variances were computed based on the $b$ sets of imputed
values. Repeating the procedure $m_{0}$ times the final variance estimates
were taken to be: (mean of the imputed variances)+$(1+1/m_{0})$(empirical
variances of the imputed point estimates).

##  5 Hypothesis testing and Power compuation

A MCEM-LRT test is developed for testing the fixed effects in a multivariate
crossover model with missing data. Here our interest is in testing the fixed
effects, however the LRT procedure may also be used to test the variance
components. Specifically, for testing the treatment effect we frame our
hypotheses as

$H_{0}:\tau_{1}=\tau_{2}=\ldots=\tau_{t}\;\;\text{vs}\;\;H_{1}:\tau_{1},\tau_{2},\ldots,\tau_{t}\text{
are not all equal }$

The LRT statistic is $\Lambda=2\log L\big{(}\hat{\theta}_{full}\big{)}-2\log
L\big{(}\hat{\theta}_{\text{reduced}}\big{)},$ where $\log
L\big{(}\hat{\theta}_{\text{full}}\big{)}$ is the likelihood function for the
full model containing all the parameters, and $\log
L\big{(}\hat{\theta}_{\text{reduced}}\big{)}$ is the likelihood function value
for the reduced model under $H_{0}$. Under the normality assumption and
certain other regularity conditions (Shao 18), when $H_{0}$ is true, $S^{2}$
follows a central $\chi^{2}_{df}$ where $\mathrm{df}$ is the difference
between the number of parameters in the full model and the reduced model. The
null hypothesis, $H_{0}$, is rejected if the observed value of $S^{2}$ exceeds
the $(1-\alpha)^{\text{th}}$ quantile of $\chi^{2}_{df}$.

For testing the period effect or the fixed effect of the response variates we
may use the respective hypotheses,

$H_{0}:\pi_{1}=\pi_{2}=\ldots=\pi_{p}\;\;\text{vs}\;\;H_{1}:\pi_{1},\pi_{2},\ldots,\pi_{p}\text{
are not all equal},$
$H_{0}:g_{1}=g_{2}=\ldots=g_{m}\;\;\text{vs}\;\;H_{1}:g_{1},g_{2},\ldots,g_{m}\text{
are not all equal }.$

Note for EM-LRT we have assumed that the central limit theorem holds for the
MCEM estimator. Thus, in the next section we carry out extensive simulations
to assess the EM-LRT test and its power for various combinations of sample
size and $\%$ of missingness.

##  6 Simulation study

We discuss a detailed simulation study to asses the performance of the MCEM
estimators and the LRT test based on these estimators. For data generation, we
assumed a crossover trial with two treatment sequences $\\{AB,BA\\}$ in two
periods. In each period four response variates were measured. The model is
represented as,

$\displaystyle y_{ijkl}$
$\displaystyle=\beta_{0}+\beta_{\text{per}}\text{Period}_{1}+\beta_{\tau}\text{Trt}_{1}+\beta_{r_{1}}\text{Res}_{1}+\beta_{r_{2}}\text{Res}_{2}$
$\displaystyle+\beta_{r_{3}}\text{Res}_{3}+s_{ik}+e_{ijkl};\quad i,j=1,2,\quad
l=1,2,3,4,\quad k=1,2,\ldots 50$

where $y_{ijkl}$ denotes the $l^{\text{th}}$ response value from the
$k^{\text{th}}$ subject in the $j^{\text{th}}$ period of the $i^{\text{th}}$
sequence; $\beta_{0}$ is the intercept; Period1 is an indicator variable
corresponding to the first period effect; $\operatorname{Trt}_{1}$ is an
indicator variable corresponding to the treatment effect $A$ and Res1, Res2,
Res3 are indicator variables corresponding to the response variates, where
Res${}_{i}=1$ for the $i^{\text{th}}$ response variate and $0$ otherwise for
$i=1,2,3$. Also, as defined before $s_{ik}$ is the subject specific random
effect and $e_{ijkl}$ is the error term.

$s_{ik}\sim N\left(0,\sigma_{s}^{2}\right),\quad e_{ijkl}\sim
N\left(0,\sigma_{e}^{2}\right);\;\;s_{ik}\perp e_{ijkl}$

In matrix notations,

$\bm{y}_{i}=\bm{X}_{i}\bm{\beta}+\bm{Z}_{i}\bm{b}_{i}+\bm{e}_{i},$

where, for $i=1,2$, $\bm{e}_{i}\sim N_{400}(\bm{0},1.44\bm{I})$,
$\bm{b}_{i}\sim N_{50}(\bm{0},0.49\bm{I})$, $\bm{X}_{i}$ and $\bm{Z}_{i}$
matrices are constructed similarly as in eq. 2. The true values of the
components of $\bm{\beta}$ are given in the first column of Table 1.

We considered different probabilities of missingness $(15\%,\,25\%,\,35\%)$ in
each sequence and the missing mechanism to be MAR. Since the non-response
model is assumed to be ignorable (Laird 10) we use a Bernoulli random variable
with probability $0.15,\,0.25,\,0.35$, respectively to generate the missing
data. The MCEM algorithm discussed in section 4 is used for parameter
estimation and SEs (standard errors) of the estimators are obtained via the
imputation method. At every iteration of the MCEM algorithm, one thousand
samples were generated for missing observations.

The simulations were repeated 500 times. Table 1 shows the average simulation
results in terms of parameter estimates, SEs, bias, and empirical coverages
for different missing percentages. Absolute average bias was computed as
$|\frac{\sum_{w=1}^{500}(\hat{\beta}_{uw}-\beta_{u})}{500}|$, where
$\hat{\beta}_{uw}$ is the $u^{\text{th}}$ component of $\hat{\bm{\beta}}$ for
the $w^{\text{th}}$ simulation and $\beta_{u}$ is the true value. Whereas, the
empirical coverage probability (ECP) was taken to be the proportion of the
times the $95\%$ asymptotic confidence interval of $\hat{\beta}_{uw}$
contained the true parameter value $\beta_{u}$ in $500$ simulations. To
compute the SEs, one hundred imputations with a burn-in of $1000$ were used.

From Table 1, one notes that absolute average biases for $15\%$ missing
responses, ranges between 2$\%$ \- 5$\%$. As we increase the $\%$ of
missingness to $35\%$, the range of the average bias change to 4$\%$ \-
11$\%$. The average MSEs of the estimates range from $1\%-3\%$ for all three
missing percentages, showing that the corresponding standard errors are also
quite low. These simulation results imply that the MCEM estimators work well
for finite samples.

Table 1: Simulation results based on five hundred samples Parameter | MLE | SE | Bias | MSE | ECP
---|---|---|---|---|---
| 15% missing responses
Intercept (4.50) | 4.53 | 0.11 | 0.03 | 0.03 | 0.91
Period1 (0.20) | 0.15 | 0.08 | 0.05 | 0.01 | 0.92
Trt1 (1.06) | 1.01 | 0.08 | 0.05 | 0.01 | 0.94
Res1 (0.46) | 0.44 | 0.11 | 0.02 | 0.03 | 0.94
Res2 (1.09) | 1.07 | 0.11 | 0.02 | 0.03 | 0.95
Res3 (0.50) | 0.47 | 0.11 | 0.03 | 0.03 | 0.96
$\sigma_{e}^{2}$ (1.44) | 1.42 | 0.07 | 0.02 | 0.01 | 0.92
$\sigma_{s}^{2}$ (0.49) | 0.52 | 0.08 | 0.03 | 0.02 | 0.93
| 25% missing responses
Intercept (4.50) | 4.54 | 0.10 | 0.04 | 0.02 | 0.90
Period1 (0.20) | 0.13 | 0.08 | 0.07 | 0.02 | 0.88
Trt1 (1.06) | 0.98 | 0.08 | 0.08 | 0.02 | 0.89
Res1 (0.46) | 0.42 | 0.11 | 0.04 | 0.02 | 0.93
Res2 (1.09) | 1.05 | 0.11 | 0.04 | 0.02 | 0.94
Res3 (0.50) | 0.46 | 0.11 | 0.04 | 0.02 | 0.95
$\sigma_{e}^{2}$ (1.44) | 1.41 | 0.06 | 0.03 | 0.01 | 0.92
$\sigma_{s}^{2}$ (0.49) | 0.54 | 0.08 | 0.05 | 0.02 | 0.92
| 35% missing responses
Intercept (4.50) | 4.55 | 0.10 | 0.05 | 0.02 | 0.89
Period1 (0.20) | 0.09 | 0.08 | 0.11 | 0.02 | 0.75
Trt1 (1.06) | 0.96 | 0.08 | 0.10 | 0.02 | 0.77
Res1 (0.46) | 0.40 | 0.10 | 0.06 | 0.02 | 0.91
Res2 (1.09) | 1.02 | 0.10 | 0.07 | 0.02 | 0.91
Res3 (0.50) | 0.43 | 0.10 | 0.07 | 0.02 | 0.93
$\sigma_{e}^{2}$ (1.44) | 1.40 | 0.06 | 0.04 | 0.01 | 0.89
$\sigma_{s}^{2}$ (0.49) | 0.58 | 0.09 | 0.09 | 0.02 | 0.92

In order to check the asymptotic normality of the MCEM estimators, empirical
probability density curve of each of the standardized estimator was compared
with the standard normal density curve. Figure 3 displays empirical PDF of the
standardized MCEM estimates with varying missing percentages and number of
simulations. From the fig. 3, note as the number of simulations increase, the
MCEM estimators (standardized one) are closer to asymptotic normality. The R
programs are available on request from the first author.

Figure 3: Empirical PDF of standardized MCEM estimates; Plots (A, C)
corresponds to cases with 100 simulated data sets when missingness is 15% and
25% respectively. Plots (B, D) corresponds to cases with 500 simulated data
sets when missingness is 15% and 25% respectively.

Further, we tested the following hypotheses using the MCEM-LRT procedure
detailed in section 5 and computed the power of the following tests:

  * (i)

Tests for response variate effects:
$H_{0}:\beta_{r_{1}}=\beta_{r_{2}}=\beta_{r_{3}}=0$ versus $H_{1}$: At least
one $\beta_{r_{i}}\neq 0$ for $i=1,2,3$.

  * (ii)

We also performed the pairwise tests of response variates,

  * (iii)

To test for the treatment effect we used, $H_{0}:\beta_{\tau}=0$ vs
$H_{1}:\beta_{\tau}\neq 0$

For empirical power computations we list the steps:

  * (a)

For data simulation we used the true values of the parameters from Table 1.
The proposed MCEM algorithm was used to find the estimates of the unknown
parameters and variance components, $\theta$, and denoted as
$\hat{\theta}_{full}$.

  * (b)

To obtain the restricted estimates we maximised the Q- function under $H_{0}$,
and computed $\hat{\theta}_{reduced}$.

  * (c)

Using $\hat{\theta}_{full}$ and $\hat{\theta}_{reduced}$, the test statistic
$T=2Q(\hat{\theta})-2Q(\hat{\theta}_{0})$ was determined and compared with
$C=\chi^{2}_{(1-\alpha,df)}$, where df is the difference between the number of
parameters in full and reduced model.

  * (d)

The steps (a) – (d) were repeated $1000$ times.

  * (e)

The empirical power was computed as the proportion of times $H_{0}$ is
rejected.

Figures 4 and 5 show the variation in the power function with respect to
sample size and $\%$ of missingness $(15\%,25\%\;\text{and}\;35\%)$. From
these plots we note that power increases with sample size (or number of total
subjects) and decreases with $\%$ of missingness. Note, power values in Figure
4 plot C are low due to the small effect size (difference between
$\beta_{r_{1}}$ and $\beta_{r_{3}}$ is small). Similarly, in Figure 5 we see
the power values decrease as we decrease the effect size (or as alternative
value of $\beta_{\tau}$ decreases).

Figure 4: Power values of MCEM LRTs; The points plotted indicate the empirical
proportion of test (by use of a nominal level $\alpha$ = 0.05) that rejected
the $H_{0}$ among 1000 simulated data sets. The plot on the left (A)
correspond to the case with testing hypothesis on significance of response
variates. Plots on the right (B, C, D) correspond to the case with testing
hypothesis on pairwise comparison of response variates. Figure 5: Power values
of MCEM LRTs; The points plotted indicate the empirical proportion of test (by
use of a nominal level $\alpha$ = 0.05) that rejected the $H_{0}$ among 1000
simulated data sets. Plots A, B, C, D correspond to the case with testing
hypothesis on treatment when alternative value of treatment effect is 1.06,
-0.7, 0.5, 0.2, respectively.

##  7 Gene case study results

Based on the exploratory analysis discussed in Section 2, we fitted the
following normal cross over model:

$\displaystyle y_{ijkl}$
$\displaystyle=\beta_{0}+\beta_{1}\text{Period}_{2}+\beta_{2}\text{Period}_{3}+\beta_{3}\mathrm{Trt}_{2}+\beta_{4}\mathrm{Trt}_{3}+\beta_{5}\text{Gene}_{2}$
$\displaystyle+\beta_{6}\text{Gene}_{3}+\beta_{7}\text{Gene}_{4}+\beta_{8}\text{Gene}_{5}+\beta_{9}\text{Gene}_{6}+\beta_{10}\text{Gene}_{7}$
$\displaystyle+\beta_{11}\text{Gene}_{8}+\beta_{12}\text{Gene}_{9}+\beta_{13}\text{Gene}_{10}+s_{ik}+e_{ijkl}$
(7)

where for the real data $i,j=1,2,3;\,l=1(1)10$ and $k=1(1)n_{i}$, where
$n_{1}=n_{2}=6$ and $n_{3}=5$. The random effects $s_{ik}$ and $e_{ijkl}$ were
assumed to be independent normal variates with variances $\sigma_{s}^{2}$ and
$\sigma_{e}^{2}$, respectively. The period, treatment and gene effects were
taken to be fixed and were represented by corresponding indicator variables.
We assumed that the missing data mechanism is MAR.

Table 2 reports the maximum likelihood estimates of the parameters
$(\bm{\beta},\sigma_{e}^{2},\sigma_{s}^{2})$ from the normal crossover model
presented in eq. 7. For the MCEM algorithm, ANOVA estimates were used as
initial values of the parameters. Standard errors of the estimates and their
corresponding p-values obtained by the MCEM-LRT test procedure are also
presented in section 4. One hundred imputations, each with burn-in of one
thousand, were used for standard error estimation. In order to check the model
adequacy computed the root mean square error (RMSE) of marginal residuals
along with AIC and BIC values. The missing data estimates were compared with
complete case analysis results.

From table 2, we observe that gene2, gene3, gene4, gene5 gene6, gene7, gene8,
gene9, gene10 were significant at level of significance $0.10$ in both
complete case and missing case results. From the negative coefficients of the
gene effects, we saw that the presence of that particular gene has a
decreasing effect on the allergic reaction. Also, we note the AIC and BIC
values of MAR responses model were much smaller as compared to complete case
analysis. However, if instead of the responses we worked with log(responses)
then the AIC, BIC values reduced further. From the RMSE of the residuals, we
saw that the best fit was obtained by fitting the normal crossover model (eq.
7) to log(responses).

We also artificially increased the number of missing responses from $9\%$ to
$21\%$ and $24\%$. The corresponding results are given in table 3. We used a
Bernoulli random variable for creating the artificially missing responses. For
the $21\%$ case, we assumed that if there is a missing response for any gene
in a specific period then all other responses are also missing. This mimicked
the property of the actual real data set. However, in the $24\%$ missingness
case we assumed that in a specific period only some gene responses (not
necessarily all) were missing. Under the higher missing percentage also, the
gene effects were seen to be highly significant. The AIC, BIC and RMSE values
for responses and log(responses) were also computed.

Table 2: Results based on the Gene data Variable | MLE | SE | p-value | MLE | SE | p-value | MLE | SE | p-value
---|---|---|---|---|---|---|---|---|---
| Complete Cases | MAR Responses | log MAR Responses
Intercept | 3.49 | 0.02 |  $<$0.0001 | 3.49 | 0.02 |  $<$0.0001 | 1.25 | 0.01 |  $<$0.0001
Period2 | 0.001 | 0.01 | 0.96 | -0.002 | 0.01 | 0.73 | -0.003 | 0.01 | 0.58
Period3 | 0.01 | 0.01 | 0.40 | 0.01 | 0.01 | 0.39 | 0.004 | 0.01 | 0.35
Trt2 | -0.01 | 0.01 | 0.70 | -0.01 | 0.01 | 0.58 | -0.003 | 0.01 | 0.57
Trt3 | 0.01 | 0.01 | 0.57 | 0.01 | 0.01 | 0.58 | 0.002 | 0.01 | 0.61
Gene2 | -0.93 | 0.03 |  $<$0.0001 | -0.94 | 0.02 |  $<$0.0001 | -0.32 | 0.01 |  $<$0.0001
Gene3 | -1.59 | 0.03 |  $<$0.0001 | -1.59 | 0.02 |  $<$0.0001 | -0.61 | 0.01 |  $<$0.0001
Gene4 | -0.33 | 0.03 |  $<$0.0001 | -0.33 | 0.02 |  $<$0.0001 | -0.10 | 0.01 |  $<$0.0001
Gene5 | -0.04 | 0.03 | 0.14 | -0.04 | 0.02 | 0.08 | -0.02 | 0.01 | 0.28
Gene6 | -1.27 | 0.03 |  $<$0.0001 | -1.27 | 0.02 |  $<$0.0001 | -0.45 | 0.01 |  $<$0.0001
Gene7 | -1.48 | 0.03 |  $<$0.0001 | -1.47 | 0.02 |  $<$0.0001 | -0.55 | 0.01 |  $<$0.0001
Gene8 | -2.03 | 0.03 |  $<$0.0001 | -2.02 | 0.02 |  $<$0.0001 | -0.87 | 0.01 |  $<$0.0001
Gene9 | -0.82 | 0.03 |  $<$0.0001 | -0.82 | 0.02 |  $<$0.0001 | -0.27 | 0.01 |  $<$0.0001
Gene10 | -0.65 | 0.03 |  $<$0.0001 | -0.63 | 0.02 |  $<$0.0001 | -0.20 | 0.01 |  $<$0.0001
$\sigma_{e}^{2}$ | 0.01 | 0.001 | - | 0.01 | 0.001 | - | 0.002 | 0.0001 | -
$\sigma_{s}^{2}$ | 0.001 | 0.0004 | - | 0.001 | 0.0004 | - | 0.0001 | 0.0001 | -
AIC | -545.5 |  |  | -1853.84 |  |  | -2686.15 |  |
BIC | -483.3 |  |  | -1786.10 |  |  | -2618.40 |  |
RMSE | 0.11 |  |  | 0.11 |  |  | 0.05 |  |
Table 3: Results based on the Gene data when we increase missing values Variable | MLE | SE | p-value | MLE | SE | p-value
---|---|---|---|---|---|---
| 21% missing |
| original responses | log responses
Intercept | 3.49 | 0.02 |  $<$0.0001 | 1.25 | 0.01 |  $<$0.0001
Period2 | -0.003 | 0.01 | 0.48 | -0.003 | 0.004 | 0.39
Period3 | 0.01 | 0.01 | 0.28 | 0.005 | 0.004 | 0.21
Treatment2 | -0.001 | 0.01 | 0.53 | -0.0004 | 0.004 | 0.57
Treatment3 | 0.003 | 0.01 | 0.74 | 0.001 | 0.004 | 0.53
Gene2 | -0.93 | 0.02 |  $<$0.0001 | -0.31 | 0.01 |  $<$0.0001
Gene3 | -1.60 | 0.02 |  $<$0.0001 | -0.61 | 0.01 |  $<$0.0001
Gene4 | -0.33 | 0.02 |  $<$0.0001 | -0.09 | 0.01 |  $<$0.0001
Gene5 | -0.04 | 0.02 | 0.06 | -0.01 | 0.01 | 0.19
Gene6 | -1.27 | 0.02 |  $<$0.0001 | -0.45 | 0.01 |  $<$0.0001
Gene7 | -1.47 | 0.02 |  $<$0.0001 | -0.55 | 0.01 |  $<$0.0001
Gene8 | -2.02 | 0.02 |  $<$0.0001 | -0.86 | 0.01 |  $<$0.0001
Gene9 | -0.82 | 0.02 |  $<$0.0001 | -0.27 | 0.01 |  $<$0.0001
Gene10 | -0.63 | 0.02 |  $<$0.0001 | -0.20 | 0.01 |  $<$0.0001
$\sigma_{e}^{2}$ | 0.01 | 0.001 | - | 0.002 | 0.0001 | -
$\sigma_{s}^{2}$ | 0.001 | 0.0003 | - | 0.0001 | 0.0001 | -
AIC | -1907.4 |  |  | -2758.3 |  |
BIC | -1839.6 |  |  | -2690.5 |  |
RMSE | 0.09 |  |  | 0.04 |  |
| 24% missing |
| original responses | log responses
Intercept | 3.49 | 0.02 |  $<$0.0001 | 1.25 | 0.01 |  $<$0.0001
Period2 | -0.01 | 0.01 | 0.24 | -0.005 | 0.005 | 0.18
Period3 | 0.004 | 0.01 | 0.45 | 0.002 | 0.005 | 0.38
Treatment2 | -0.004 | 0.01 | 0.45 | -0.002 | 0.005 | 0.42
Treatment3 | 0.01 | 0.01 | 0.29 | 0.005 | 0.005 | 0.25
Gene2 | -0.93 | 0.02 |  $<$0.0001 | -0.31 | 0.01 |  $<$0.0001
Gene3 | -1.58 | 0.02 |  $<$0.0001 | -0.60 | 0.01 |  $<$0.0001
Gene4 | -0.33 | 0.02 |  $<$0.0001 | -0.10 | 0.01 |  $<$0.0001
Gene5 | -0.03 | 0.02 | 0.08 | -0.01 | 0.01 | 0.21
Gene6 | -1.26 | 0.02 |  $<$0.0001 | -0.45 | 0.01 |  $<$0.0001
Gene7 | -1.47 | 0.02 |  $<$0.0001 | -0.55 | 0.01 |  $<$0.0001
Gene8 | -2.02 | 0.02 |  $<$0.0001 | -0.87 | 0.01 |  $<$0.0001
Gene9 | -0.82 | 0.02 |  $<$0.0001 | -0.27 | 0.01 |  $<$0.0001
Gene10 | -0.64 | 0.02 |  $<$0.0001 | -0.20 | 0.01 |  $<$0.0001
$\sigma_{e}^{2}$ | 0.01 | 0.001 | - | 0.003 | 0.0001 | -
$\sigma_{s}^{2}$ | 0.001 | 0.0003 | - | 0.0001 | 0.00004 | -
AIC | -1830.9 |  |  | -2652.6 |  |
BIC | -1763.2 |  |  | -2584.8 |  |
RMSE | 0.10 |  |  | 0.04 |  |

##  8 Concluding remarks

In this paper we studied statistical model to analyse incomplete data in
multivariate crossover setup. We presented the MCEM-LRT to test the fixed
effects in the proposed set up, and exhibited the power analysis. This idea
can be extended to non-ignorable missing (NMAR) case, where the parametric
model for missing or the the missing data mechanism needs to be incorporated
into the complete data log likelihood. As a consequence the estimation problem
becomes more complex.

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

Conflict of Interest The authors have declared no conflict of interest.

## References

  * Basu and Santra [2010] Sanjib Basu and Sourav Santra.  A joint model for incomplete data in crossover trials.  _Journal of Statistical Planning and Inference_ , 140(10):2839–2845, 2010.
  * 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.
  * Chow and Shao [1997] Shein-chung Chow and Jun Shao.  Statistical methods for two-sequence three-period cross-over designs with incomplete data.  _Statistics in Medicine_ , 16:1031–1039, 1997.
  * Clough and Barrett [2016] Emily Clough and Tanya Barrett.  The Gene Expression Omnibus database.  _Methods in Molecular Biology_ , 1418:93–110, 2016.
  * Els and Ryan [2000] Goetghebeur Els and Louise Ryan.  Semiparametric regression analysis of interval-censored data.  _Biometrics_ , 56(4):1139–1144, 2000.
  * 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\.
  * Ho et al. [2012] Weang K. Ho, John N.S. Matthews, Robin Henderson, Daniel Farewell, and Lauren R. Rodgers.  Dropouts in the AB/BA crossover design.  _Statistics in Medicine_ , 31(16):1675–1687, 2012\.
  * Ibrahim and Molenberghs [2009] Joseph G. Ibrahim and Geert Molenberghs.  Missing data methods in longitudinal studies: A review.  _Test_ , 18(1):1–43, 2009.
  * Jones and Kenward [2003] Byron Jones and Michael G. Kenward.  _Design and Analysis of Cross-Over Trials_.  Chapman & Hall/CRC, second edition, 2003.
  * Laird [1988] Nan M. Laird.  Missing data in longitudinal studies.  _Statistics in Medicine_ , 7:305–315, 1988.
  * Leaker et al. [2016] BR Leaker, VA Malkov, R Mogg, MK Ruddy, GC Nicholson, AJ 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.
  * Matthews and Henderson [2013] John N.S. Matthews and Robin Henderson.  Two-period, two-treatment crossover designs subject to non-ignorable missing data.  _Biostatistics_ , 14(4):626–638, 2013.
  * Patel [1985] H. I. Patel.  Analysis of incomplete data in a two-period crossover design with reference to clinical trials.  _Biometrika_ , 72(2):411–418, 1985.
  * 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\.
  * Richardson and Flack [1992] Barbra A. Richardson and Virginia F. Flack.  The analysis of incomplete data in the three-period two-treatment crossover design for clinical trials.  _Controlled Clinical Trials_ , 13(5):381, 1992\.
  * Rosenkranz [2015] Gerd K. Rosenkranz.  Analysis of cross-over studies with missing data.  _Statistical Methods in Medical Research_ , 24(4):420–433, 2015.
  * Senn [2002] Stephen Senn.  _Cross-Over Trials in Clinical Research._ John Wiley & Sons, Ltd., 2002.
  * Shao [2003] Jun Shao.  _Mathematical Statistics_.  Springer International Publishing, second edition, 2003.
  * Tudor et al. [2000] Gail E. Tudor, Gary G. Koch, and Diane Catellier.  20 Statistical methods for crossover designs in bioenvironmental and public health studies.  _Handbook of Statistics_ , 18(1980):571–614, 2000\.
  * Wei and Tanner [1990] Greg C.G. Wei and Martin A. Tanner.  A Monte Carlo Implementation of the EM Algorithm and the Poor Man’s Data Augmentation Algorithm.  _Journal of the American Statistical Association_ , 85(411):699–704, 1990.