# Incompletely observed nonparametric factorial designs with repeated
measurements: A wild bootstrap approach

Lubna Amro1∗ 111
1Mathematical Statistics and Applications in Industry, Faculty of Statistics,
Technical University of Dortmund, Germany
2Institute of Biometry and Clinical Epidemiology, Charité—Universitätsmedizin
Berlin, Berlin, Germany
3 Berlin Institute of Health (BIH), Berlin, Germany Frank Konietschke2,3 222
Corresponding author:
Lubna Amro, Mathematical Statistics and Applications in Industry, Faculty of
Statistics, Technical University of Dortmund, Germany.
Email<EMAIL_ADDRESS>Markus Pauly1 333

###### Abstract

In many life science experiments or medical studies, subjects are repeatedly
observed and measurements are collected in factorial designs with multivariate
data. The analysis of such multivariate data is typically based on
multivariate analysis of variance (MANOVA) or mixed models, requiring complete
data, and certain assumption on the underlying parametric distribution such as
continuity or a specific covariance structure, e.g., compound symmetry.
However, these methods are  usually not applicable when discrete data or even
ordered categorical data are present. In such cases, nonparametric rank-based
methods that do not require stringent distributional assumptions are the
preferred choice. However, in the multivariate case, most rank-based
approaches have only been developed for complete observations. It is the aim
of this work is to develop asymptotic correct procedures that are capable of
handling missing values, allowing for singular covariance matrices and are
applicable for ordinal or ordered categorical data. This is achieved by
applying a wild bootstrap procedure in combination with quadratic form-type
test statistics. Beyond proving their asymptotic correctness, extensive
simulation studies validate their applicability for small samples. Finally,
two real data examples are analyzed.

###### keywords:

Missing Values; Wild bootstrap; Rank tests; Repeated measures; Ordered
categorical data; Nonparametric hypotheses.

##  1 Introduction

Factorial designs have a long history in several scientific fields such as
biology, ecology or medicine (Traux and Carkhuff, 1965; Nesnow et al., 1998;
Wildsmith et al., 2001; Biederman et al., 2008; Lekberg et al., 2018). The
reasons are simple: They are an efficient way to study main and interaction
effects between different factors. Typically, such designs are inferred by
parametric mean-based procedures such as linear mixed models, or multivariate
analysis of variance (MANOVA). These procedures,  however, rely on restrictive
distributional assumptions, such as multivariate normality, or special
dependencies (Lawley, 1938; Bartlett, 1939; Davis, 2002; Johnson and Wichern,
2007; Fitzmaurice et al., 2012).
However, assessing multivariate normality or a specific type of covariance
matrix is difficult  in practice (Micceri, 1989; Qian and Huang, 2005;
Bourlier et al., 2008; Xu and Cui, 2008), especially when the sample sizes are
small. In particular, Erceg-Hurn and Mirosevich (2008) pointed out that
”Researchers relying on statistical tests (e.g., Levene’s test) to identify
assumption violations may frequently fail to detect deviations from normality
and homoscedasticity that are large enough to seriously affect the type-I
error rate and power of classic parametric tests”(p. 600). In fact, classical
MANOVA procedures are usually non robust to deviations and may result in
inaccurate decisions caused by possibly conservative or inflated type-I error
rates (Arnau et al., 2012; Pesarin and Salmaso, 2012; Brombin et al., 2013;
Konietschke et al., 2015; Pauly et al., 2015b; Friedrich et al., 2017a;
Friedrich and Pauly, 2018). Moreover, if ordinal, or ordered categorical data
are present, mean-based approaches are not applicable. A tempting alternative
are nonparametric rank-based methods, which are applicable for non-normal
data, in particular discrete data or even ordered categorical data. Their key
features are their robustness and their invariance under monotone
transformations of the data. Consequently, Akritas and Arnold (1994); Akritas
and Brunner (1997); Munzel and Brunner (2000); Brunner and Puri (2001);
Akritas (2011); Brunner et al. (2017), and Friedrich et al. (2017b) proposed
nonparametric ranking methods for all kinds of factorial designs. Thereby,
hypotheses are no longer formulated in terms of means (which do not exist for
ordinal data), distribution functions are used instead. However, all these
methods are only applicable for completely observed factorial designs data and
cannot be used to analyze multivariate data with missing values.

In contrast, there are only a few methods which are applicable in case of
missing values, require no parametric assumptions and also lead to valid
inferences in case of arbitrary covariance structures, skewed distributions or
unbalanced experimental designs.  Examples in the nonparametric framework for
the matched pairs designs include the ranking methods proposed by Akritas et
al. (2002); Konietschke et al. (2012); Fong et al. (2018). A promising
approach for factorial designs with repeated measurements is given by the
proposals of Brunner et al. (1999) and Domhof et al. (2002) who recommended
two types of quadratic forms for testing nonparametric hypotheses in terms of
distribution functions: rank based Wald and ANOVA-type statistics. The Wald-
type test is an asymptotically valid test which usually needs large sample
sizes to obtain accurate test decisions, see Brunner (2001); Friedrich et al.
(2017b) for the case of complete observations and our simulation study below.
Apart from that, the ANOVA-type test is based upon an approximation of its
distribution with scaled $\chi^{2}$-distributions. As the latter does not
coincide with the ANOVA-type test limiting distribution under $H_{0}$, the
test is in general not asymptotically correct and also exhibits a more or less
conservative behaviour under small sample sizes, see Brunner (2001) for the
case of complete observations. Another applicable technique is the
nonparametric imputation method proposed by Gao (2007), which possesses a good
type-I error control but a lower power behavior than the methods by Brunner et
al. (1999) and Domhof et al. (2002), see the simulation study of Gao (2007).

Thus, the aims of the present paper are (i) to provide statistical tests  for
hypotheses formulated in distribution functions that are capable of treating
missing values in factorial designs, (ii) work without parametric assumptions
such as continuity of the distribution functions, or nonsingular covariance
matrices, (iii) are asymptotically correct while (iv) showing a satisfactorily
type-I error control and good power properties. To accomplish this, we
propose three different quadratic-form type test statistics and equip them
with a nonparametric Wild bootstrap procedure for calculating critical values
to enhance their small sample performance.
The paper is organized as follows: First we introduce the statistical model
and the hypotheses of interest in the next section. In Section 3, we introduce
the test statistics and analyze their asymptotic behaviour. In Section 4, the
proposed wild bootstrap technique is explained. Section 5 displays the results
from our extensive simulation study and two real data examples from a
fluvoxamine study and a skin disorder clinical trial are analyzed in Section
6. All proofs as well as additional simulation results can be found in the
supplement.

To facilitate the presentation, we introduce the following notations: Let
$\boldsymbol{I}_{d}$ denote the d-dimensional unit matrix,
$\boldsymbol{J}_{d}$ denote the $d\times d$ matrix of $\boldsymbol{1}$’s i.e.
$\boldsymbol{J}_{d}=\boldsymbol{1}_{d}\boldsymbol{1}_{d}^{T}$, where
$\boldsymbol{1}_{d}=(1,...,1)^{T}_{d\times 1}$ denotes the d-dimensional
column vector and
$\boldsymbol{P}_{d}=\boldsymbol{I}_{d}-\frac{1}{d}\boldsymbol{J}_{d}$ is the
so-called d-dimensional centering matrix. Finally, by
$\boldsymbol{A}\otimes\boldsymbol{B}$ we denote the Kronecker product of the
matrices $\boldsymbol{A}$ and $\boldsymbol{B}$.

##  2 Statistical model and nonparametric hypotheses

We consider a nonparametric repeated measures model with $a$ independent and
potentially unbalanced treatment groups and $d$ different time points given by
independent random vectors

$\displaystyle\boldsymbol{X}_{ik}=(X_{i1k},...,X_{idk})^{T},\quad
i=1,...,a,\quad k=1,...,n_{i},$ (2.1)

where $X_{ijk}\sim
F_{ij}(x)=\frac{1}{2}[F^{+}_{ij}(x)+F^{-}_{ij}(x)],i=1,...,a,j=1,...,d,k=1,...,n_{i}$.
Here, $F^{+}_{ij}(x)=P(X_{ijk}\leq x)$ is the right continuous version and
$F^{-}_{ij}(x)=P(X_{ijk}<x)$ is the left continuous version of the
distribution function. Using the normalized version $F_{ij}(x)$ is useful for
handling ties and including continuous as well as discontinuous distribution
functions (Brunner et al., 2018). To include the case of missing values, we
follow the notation of Brunner et al. (1999) and let

$\displaystyle\lambda_{ijk}=\begin{cases}1,\quad\text{if $X_{ijk}$ is
observed}\\\ 0,\quad\text{if $X_{ijk}$ is non-
observed}\end{cases}i=1,...,a,j=1,...,d,k=1,...,n_{i}.$ (2.2)

Moreover, let $n=\sum_{i=1}^{a}n_{i}$ denote the total number of subjects and
let
$N=\sum_{i=1}^{a}\sum_{j=1}^{d}\sum_{k=1}^{n_{i}}\lambda_{ijk}$ denote the
total number of observations.

To formulate the null hypothesis in this nonparametric setup, let
$\boldsymbol{F}=(F_{11},...,F_{ad})^{T}$ denote the vector of the distribution
functions $F_{ij},i=1,...,a,j=1,...,d.$ and $\boldsymbol{C}$ denote a contrast
matrix , i.e., $\boldsymbol{C1}=\boldsymbol{0}$ where
$\boldsymbol{1}=(1,...,1)^{T}$ and $\boldsymbol{0}=(0,...,0)^{T}$. Then, the
null hypotheses are formulated by
$H_{0}:\\{\boldsymbol{C}\boldsymbol{F}=\boldsymbol{0}\\}$. This framework
covers different factorial repeated measures designs. For example, the
hypothesis of no treatment group effect, i.e.
$H_{0}^{G}:\\{\bar{F}_{1.}=...=\bar{F}_{a.}\\}$ , is equivalently written in
matrix notation as
$H_{0}^{G}:\\{(\boldsymbol{P}_{a}\otimes\frac{1}{d}\boldsymbol{1}_{d}^{T})\boldsymbol{F}=\boldsymbol{0}\\}$.
Similarly, the hypothesis of no time effect, i.e.
$H_{0}^{T}:\\{\bar{F}_{.1}=...=\bar{F}_{.d}\\}$ is equivalently written as
$H_{0}^{T}:\\{(\frac{1}{a}\boldsymbol{1}_{a}^{T}\otimes\boldsymbol{P}_{d})\boldsymbol{F}=\boldsymbol{0}\\}$,
and the hypothesis of no interaction effect between treatment and time is
written as
$H_{0}^{GT}:\\{(\boldsymbol{P}_{a}\otimes\boldsymbol{P}_{d})\boldsymbol{F}=\boldsymbol{0}\\}$.

To entail a parameter for describing differences between distributions,
Brunner et al. (1999) and Domhof et al. (2002) considered the relative
marginal effects

$\displaystyle p_{ij}=\int{H(x)dF_{ij}(x)},$ (2.3)

where
$H(x)=N^{-1}\sum_{i=1}^{a}\sum_{j=1}^{d}\sum_{k=1}^{n_{i}}\lambda_{ijk}F_{ij}(x)$
is the weighted average of all distribution functions in the experiment.
Estimators thereof are given by plugging-in the empirical versions of
$F_{ij}(x)$ and $H(x)$

$\displaystyle\hat{F}_{ij}(x)$
$\displaystyle=\frac{1}{\lambda_{ij.}}\sum_{k=1}^{n_{i}}\lambda_{ijk}\hat{F}_{ijk}(x)$
(2.4)
$\displaystyle=\frac{1}{\lambda_{ij.}}\sum_{k=1}^{n_{i}}\lambda_{ijk}c(x-X_{ijk}),\quad\lambda_{ij.}=\sum_{k=1}^{n_{i}}\lambda_{ijk},$
(2.5) $\displaystyle\hat{H}(x)$
$\displaystyle=\frac{1}{N}\sum_{i=1}^{a}\sum_{j=1}^{d}\sum_{k=1}^{n_{i}}c(x-X_{ijk}),$
(2.6)

where, $c(u)$ is the normalized version of the counting function, i.e.
$c(u)=0,1/2$ or $1$ according as $u<0,u=0$ or $u>0$ and $c(x-X_{ijk})$ is
assumed to equal $0$ if the observation $X_{ijk}$ is missing. Note that
$\hat{F}_{ij}(x)=0$ in case of $\lambda_{ij.}=0$. Thus, the relative marginal
effects $p_{ij}$ are estimated by

$\displaystyle\hat{p}_{ij}$
$\displaystyle=\int\hat{H}d\hat{F}_{ij}=\frac{1}{\lambda_{ij.}}\sum_{k=1}^{n_{i}}\lambda_{ijk}\hat{H}(X_{ijk})$
(2.7)
$\displaystyle=\frac{1}{\lambda_{ij.}}\sum_{k=1}^{n_{i}}\lambda_{ijk}\frac{1}{N}(R_{ijk}-\frac{1}{2}),$
(2.8)

where $R_{ijk}$ is the mid-rank of $X_{ijk}$ among all $N$ dependent and
independent observations. For convenience, the relative marginal effect
estimators are collected in the vector
$\hat{\boldsymbol{p}}=(\hat{p}_{11},...,\hat{p}_{ad})^{T}$.

In the sequel, we derive asymptotic theory under the following sample size
assumption and missing values:

######  Assumption 1.

$\frac{\lambda_{ij.}}{n}\rightarrow\kappa_{i}\in(0,1)\quad i=1,...,a,\quad
j=1,...,d$, as $n_{0}\equiv min\\{\lambda_{ij.}\\}\rightarrow\infty$.

This ensures that there are ’enough’ subjects without missing values in each
group.

Asymptotic distribution. Brunner et al. (1999) derived the asymptotic
distribution of $\sqrt{n}\boldsymbol{C}\hat{\boldsymbol{p}}$ under the null
hypothesis $H_{0}:\boldsymbol{C}\boldsymbol{F}=\boldsymbol{0}$. Setting
$\hat{\boldsymbol{Y}}_{ijk}=\hat{H}(X_{ijk})$ and
$\boldsymbol{Y}_{ijk}=H(X_{ijk})$, they proved that under $H_{0}$,
$\sqrt{n}\boldsymbol{C}\hat{\boldsymbol{p}}$ and
$\sqrt{n}\boldsymbol{C}\bar{\boldsymbol{Y}}_{.}$ are asymptotically
equivalent. Accordingly, they showed that
$\sqrt{n}\boldsymbol{C}\hat{\boldsymbol{p}}$ follows under $H_{0}$,
asymptotically, a multivariate normal distribution with expectation
$\boldsymbol{0}$ and covariance matrix
$\boldsymbol{C\boldsymbol{V}_{n}C^{T}}$. Here

$\displaystyle\boldsymbol{V}_{n}=\bigoplus\limits_{i=1}^{d}\frac{n}{n_{i}}\boldsymbol{V}_{i}=\bigoplus\limits_{i=1}^{d}\frac{n}{n_{i}^{2}}\sum_{k=1}^{n_{i}}\boldsymbol{\Lambda}_{ik}\boldsymbol{V}_{ik}\boldsymbol{\Lambda}_{ik},$
(2.9)

where
$\boldsymbol{\Lambda}_{ik}=n_{i}diag\\{\frac{\lambda_{i1k}}{\lambda_{i1.}},...,\frac{\lambda_{idk}}{\lambda_{id.}}\\}$,
and $\boldsymbol{V}_{ik}=Cov(\boldsymbol{Y}_{ik})$ denotes the covariance
matrix of the random vectors
$\boldsymbol{Y}_{ik}=(H(X_{i1k}),...,H(X_{idk}))^{T}$.

Brunner et al. (1999) propose to estimate the covariance matrix
$\boldsymbol{V}_{n}$ by

$\displaystyle\boldsymbol{\hat{V}}_{n}=\bigoplus\limits_{i=1}^{r}\frac{n}{n_{i}}\boldsymbol{\hat{V}}_{i},$
(2.10)

where $\boldsymbol{\hat{V}}_{i}=[\hat{v}_{i}(j,j^{\prime})]$ with diagonal and
off-diagonal elements $\hat{v}_{i}(j,j)$ and $\hat{v}_{i}(j,j^{\prime})$,
respectively, defined as

$\displaystyle\hat{v}_{i}(j,j)=\frac{n_{i}\sum_{k=1}^{n_{i}}\lambda_{ijk}[R_{ijk}-\bar{R}_{ij.}]^{2}}{(N^{2})\lambda_{ij.}(\lambda_{ij.}-1)},$
(2.11)

$\displaystyle\hat{v}_{i}(j,j^{\prime})=\frac{n_{i}\sum_{k=1}^{n_{i}}\lambda_{ijk}\lambda_{ij^{\prime}k}[(R_{ijk}-\bar{R}_{ij.})(R_{ij^{\prime}k}-\bar{R}_{ij^{\prime}.})]}{(N^{2})((\lambda_{ij.}-1)(\lambda_{ij^{\prime}.}-1)+\Delta_{i,jj^{\prime}}-1)}.$
(2.12)

Here,
$\bar{R}_{ij.}=\frac{1}{\lambda_{ij.}}\sum_{k=1}^{n_{i}}\lambda_{ijk}(R_{ijk})$,
and
$\Delta_{i,jj^{\prime}}=\sum_{k=1}^{n_{i}}\lambda_{ijk}\lambda_{ij^{\prime}k}$.
Under Assumption 1, Brunner et al. (1999) have shown that
$\boldsymbol{\hat{V}}_{i}$ is a consistent estimator of $\boldsymbol{V}_{i}$.

##  3 Statistics and Asymptotics

In this section, we propose three different quadratic forms for testing the
null hypothesis: a Wald-type statistic (WTS), an ANOVA-type statistic (ATS) as
suggested in Brunner et al. (1999) and Domhof et al. (2002) as well as a
modified ANOVA-type statistic (MATS). For mean-based analyses, the latter was
proposed by Friedrich and Pauly (2018).

The rank version of the WTS is defined as

$\displaystyle
T_{W}=n\boldsymbol{\hat{p}}^{T}\boldsymbol{C}^{T}[\boldsymbol{C}\boldsymbol{\hat{V}}_{n}\boldsymbol{C}^{T}]^{+}\boldsymbol{C}\boldsymbol{\hat{p}},$
(3.1)

where $[\boldsymbol{B}]^{+}$ denotes the Moore-penrose inverse of a matrix
$\boldsymbol{B}$. Its asymptotic null distribution is summarized below.

######  Theorem 3.1.

Under Assumption $(1)$ and $\boldsymbol{V}_{n}\rightarrow\boldsymbol{V}>0$ as
$n_{0}\rightarrow\infty$, the statistic $T_{W}$ has under the null hypothesis
$H_{0}:\boldsymbol{C}\boldsymbol{F}=\boldsymbol{0}$, asymptotically, as
$n_{0}\rightarrow\infty$, a central $\chi^{2}_{f}$-distribution with
$f=rank(\boldsymbol{C})$ degrees of freedom.

Wald-type statistics of similar form are used in many diferent situations,
e.g., in heteroscedatic mean-based analyses (Krishnamoorthy and Lu, 2010; Xu
et al., 2013; Konietschke et al., 2015; Friedrich and Pauly, 2018; Amro et
al., 2019b) and even more complex regression models in survival analyses
(Martinussen and Scheike, 2007; Dobler et al., 2019b). However, the
convergence of the WTS to its limiting $\chi^{2}$-distribution is usually slow
and large sample sizes are required to obtain adequate results (Vallejo et
al., 2010; Konietschke et al., 2015; Pauly et al., 2015a; Smaga, 2017). Thus,
Brunner et al. (1999) and Domhof et al. (2002) proposed an alternative
quadratic form by deleting the variance $\boldsymbol{\hat{V}}_{n}$ involved in
the computation of the WTS, resulting in the ANOVA-type statistic (ATS)
defined as

$\displaystyle
T_{A}=\frac{1}{tr(\boldsymbol{T}\boldsymbol{\hat{V}}_{n})}n\boldsymbol{\hat{p}}^{T}\boldsymbol{T}\boldsymbol{\hat{p}},$
(3.2)

where
$\boldsymbol{T}=\boldsymbol{C}^{T}[\boldsymbol{C}\boldsymbol{C}^{T}]^{+}\boldsymbol{C}$
is a projection matrix. Note that
$H_{0}:\boldsymbol{T}\boldsymbol{F}=\boldsymbol{0}\Leftrightarrow\boldsymbol{C}\boldsymbol{F}=\boldsymbol{0}$
because $\boldsymbol{C}^{T}[\boldsymbol{C}\boldsymbol{C}^{T}]^{+}$ is a
generalized inverse of $\boldsymbol{C}$.
It is also worth to note, that different to $T_{W}$, the asymptotic
distribution of $T_{A}$ can also be derived if $\boldsymbol{V}$ is singular.

######  Theorem 3.2.

Under Assumption $(1)$ and under the null hypothesis
$H_{0}:\boldsymbol{C}\boldsymbol{F}=\boldsymbol{0}$, the test statistic
$T_{A}$ has asymptotically, as $n_{0}\rightarrow\infty$, the same distribution
as the random variable

$\displaystyle\boldsymbol{A}=\sum_{i=1}^{a}\sum_{j=1}^{d}\zeta_{ij}B_{ij}/tr(\boldsymbol{T}\boldsymbol{V}_{n}),$
(3.3)

where $B_{ij}\overset{\text{i.i.d}}{\sim}\chi^{2}_{1}$ and the weights
$\zeta_{ij}$ are the eigenvalues of $\boldsymbol{T}\boldsymbol{V}_{n}$.

The limiting distribution is approximated by a scaled $g\chi^{2}_{f}$
distribution. Brunner et al. (1999) proposed a Box (1954)-type approximation
such that the first two moments of $T_{A}$ and $g\chi^{2}_{f}$ approximately
coincide where $f$ can be estimated by
$\hat{f}=\frac{(tr(\boldsymbol{T}\boldsymbol{\hat{V}}_{n}))^{2}}{tr(\boldsymbol{T}\boldsymbol{\hat{V}}_{n}\boldsymbol{T}\boldsymbol{\hat{V}}_{n})}$.
Thus, the distribution of $T_{A}$ is approximated by a central
$F(\hat{f},\infty)$-distribution. The corresponding ANOVA-type test
$\phi_{A}=\mathds{1}\\{T_{A}>F_{\alpha}(\hat{f},\infty)\\}$, where
$F_{\alpha}(\hat{f},\infty)$ denotes the $(1-\alpha)$-quantile of the
$F(\hat{f},\infty)$-distribution.

Despite its advantage of being applicable in case of singular covariance
matrices, the ATS has the drawback of being an approximative test and thus,
even its asymptotic exactness cannot be guaranteed.

Another possible test statistic is the modified version of the ANOVA type-
statistic (MATS) that was developed by Friedrich and Pauly (2018) for mean-
based MANOVA models. Here, we transfer it to the ranked-based set-up, where it
is given by

$\displaystyle
T_{M}=n\boldsymbol{\hat{p}}^{T}\boldsymbol{C}^{T}[\boldsymbol{C}\hat{\boldsymbol{D}}_{n}\boldsymbol{C}^{T}]^{+}\boldsymbol{C}\boldsymbol{\hat{p}},$
(3.4)

where $\hat{\boldsymbol{D}}_{n}=diag(\boldsymbol{\hat{V}}_{n})_{ii}$. In this
way it is a compromise between the WTS and the ATS as it only uses the
diagonal entries of $\hat{\boldsymbol{V}}_{n}$ for the multivariate
studentization. In fact,  the nonsingularity assumption of $\boldsymbol{V}$ is
not needed to derive its asymptotics. It is replaced by the weaker requirement
that $\mathbf{D}=diag(\boldsymbol{V})_{ii}>0$, which is fulfilled when all
it’s diagonal elements $v_{i}(j)$ are positive. This is the same as
$Var(H(X_{ij1}))>0$ for all $i\in\\{1,...,a\\},j\in\\{1,...,d\\}$. It is
supposed to be met in a wide range of applicable settings, except only cases
where, for example, at least one time point of any of the treatment groups is
a discrete variable with very few distinct values.

######  Theorem 3.3.

Under Assumption $(1)$ and  assuming that $v_{i}(j,j)>0$ for all
$i\in\\{1,...,a\\},j\in\\{1,...,d\\}$, the test statistic $T_{M}$, has under
the null hypothesis $H_{0}:\boldsymbol{C}\boldsymbol{F}=\boldsymbol{0}$,
asymptotically, as $n_{0}\to\infty$, the same distribution as the random
variable

$\boldsymbol{M}=\sum_{i=1}^{a}\sum_{j=1}^{d}\tilde{\zeta}_{ij}\tilde{B}_{ij},$

where $\tilde{B}_{ij}\overset{\text{i.i.d}}{\sim}\chi_{1}^{2}$ and the weights
$\tilde{\zeta}_{ij}$ are the eigenvalues of
$\boldsymbol{V}_{n}^{1/2}\boldsymbol{C}^{T}[\boldsymbol{C}\boldsymbol{D}\boldsymbol{C}^{T}]^{+}\mathbf{C}\boldsymbol{V}_{n}^{1/2}$
and $\mathbf{D}=diag(\boldsymbol{V})_{ii}$.

Since, the limit distributions of both, the ATS and the MATS depend on unknown
weights $\zeta_{i}$ and $\tilde{\zeta}_{i}$, we cannot directly calculate
critical values. In addition, the $\chi^{2}_{f}-$approximation to $T_{W}$ is
rather slow. To this end, we develop asymptotically correct testing procedures
based on bootstrap versions of $T_{W},T_{A},$ and $T_{M}$ in the subsequent
section.

##  4 Wild bootstrap approach

We consider a wild bootstrap approach to derive new asymptotically valid
testing procedures with good finite sample properties. To this end, let
$\boldsymbol{Z}_{ik}=(\boldsymbol{R}_{ik}-\bar{\boldsymbol{R}}_{i.})$ denote
the centered rank vectors, where
$\boldsymbol{R}_{ik}=(R_{i1k},...,R_{idk})^{T}$,
$\bar{\boldsymbol{R}}_{i.}=(\bar{R}_{i1.},...,\bar{R}_{id.})^{T}$ and
$\bar{R}_{ij.}=\frac{1}{\lambda_{ij.}}\sum_{k=1}^{n_{i}}\lambda_{ijk}(R_{ijk})$.
Moreover, let $W_{ik}$ denote independent and identically distributed random
weights with $E(W_{ik})=0$ and $Var(W_{ik})=1$. Although, there are different
possible choices for these random weights (Mammen, 1992; Davidson and
Flachaire, 2008), some particular choices have become popular. Following the
investigations in Friedrich et al. (2017b) for the corresponding complete case
scenario, we use Rademacher random variables, which are defined by
$P(W_{ik}=-1)=P(W_{ik}=1)=1/2$. Then, a wild bootstrap sample is defined as

$\displaystyle\boldsymbol{Z}_{ik}^{*}=W_{ik}.\boldsymbol{Z}_{ik},\quad
i=1,...,a,\quad k=1,...,n_{i}.$ (4.1)

In other words, $\boldsymbol{Z}_{ik}^{*}$ is a symmetrization of the rank
vector $\boldsymbol{Z}_{ik}$. Now, we can define the bootstrap version of the
relative effect estimator $\hat{p}_{ij}$

$\displaystyle\hat{p}_{ij}^{*}=\frac{1}{\lambda_{ij.}}\sum_{k=1}^{n_{i}}\lambda_{ijk}\frac{1}{N}(Z_{ijk}^{*})=\frac{1}{\lambda_{ij.}}\sum_{k=1}^{n_{i}}\lambda_{ijk}\frac{1}{N}(W_{ik}(R_{ijk}-\bar{R}_{ij.})).$
(4.2)

We pool them into the vector
$\hat{\boldsymbol{p}}^{*}=(\hat{p}_{11}^{*},...,\hat{p}_{ad}^{*})^{T}$ to
obtain the bootstrap counterpart of $\hat{\boldsymbol{p}}$.
In the same way, the bootstrap covariance matrix estimator
$\boldsymbol{\hat{V}}_{i}^{*}=[\hat{v}_{i}^{*}(j,j^{\prime})]$ with diagonal
and off-diagonal elements $\hat{v}_{i}^{*}(j,j)$ and
$\hat{v}_{i}^{*}(j,j^{\prime})$, respectively, is given by

$\displaystyle\hat{v}_{i}^{*}(j,j)=\frac{n_{i}\sum_{k=1}^{n_{i}}\lambda_{ijk}[Z^{*}_{ijk}-\bar{Z}^{*}_{ij.}]^{2}}{(N^{2})\lambda_{ij.}(\lambda_{ij.}-1)},$
(4.3)

$\displaystyle\hat{v}_{i}^{*}(j,j^{\prime})=\frac{n_{i}\sum_{k=1}^{n_{i}}\lambda_{ijk}\lambda_{ij^{\prime}k}[(Z^{*}_{ijk}-\bar{Z}^{*}_{ij.})(Z^{*}_{ij^{\prime}k}-\bar{Z}^{*}_{ij^{\prime}.})]}{(N^{2})((\lambda_{ij.}-1)(\lambda_{ij^{\prime}.}-1)+\Delta_{i,jj^{\prime}}-1)}.$
(4.4)

From this, the bootstrapped versions of the quadratic forms, i.e. the Wald-
type statistic $T_{W}^{*}$, the ANOVA-type statistic $T_{A}^{*}$ and the
modified ANOVA-type statistic $T_{M}^{*}$ are computed as

$\displaystyle
T_{W}^{*}=n\boldsymbol{\hat{p}}^{*T}\boldsymbol{C}^{T}[\boldsymbol{C}\boldsymbol{\hat{V}}_{n}^{*}\boldsymbol{C}^{T}]^{+}\boldsymbol{C}\boldsymbol{\hat{p}}^{*},$
(4.5) $\displaystyle
T_{A}^{*}=\frac{1}{tr(\boldsymbol{T}\boldsymbol{\hat{V}}^{*}_{n})}n\boldsymbol{\hat{p}}^{*T}\boldsymbol{T}\boldsymbol{\hat{p}}^{*},$
(4.6) $\displaystyle
T_{M}^{*}=n\boldsymbol{\hat{p}}^{*T}\boldsymbol{C}^{T}[\boldsymbol{C}\hat{\boldsymbol{D}}^{*}_{n}\boldsymbol{C}^{T}]^{+}\boldsymbol{C}\boldsymbol{\hat{p}}^{*},$
(4.7)

where
$\boldsymbol{\hat{V}}_{n}^{*}=\bigoplus\limits_{i=1}^{a}\frac{n}{n_{i}}\boldsymbol{\hat{V}}_{i}^{*}$
and $\hat{\boldsymbol{D}}^{*}_{n}=diag(\boldsymbol{\hat{V}}_{n}^{*})_{ii}$. To
get an asymptotically valid bootstrap test, we have to assure that the
conditional distribution of the Wald, ANOVA-type, and MATS-type bootstrap
statistics $T_{W}^{*}$, $T_{A}^{*}$, and $T_{M}^{*}$ approximate the null
distribution of $T_{W}$, $T_{A}$, and $T_{M}$, respectively.

######  Theorem 4.1.

Under Assumption $(1)$, the following results hold:

  1. 1.

For $i=1,...,a$, the conditional distribution of
$\sqrt{n}\boldsymbol{\hat{p}}^{*}_{i}$, given the data $\boldsymbol{X}$,
converges weakly to the multivariate
$N(0,\kappa_{i}^{-1}\boldsymbol{V}_{i})$-distribution in probability.

  2. 2.

The conditional distribution of $\sqrt{n}\boldsymbol{\hat{p}}^{*}$, given the
data $\boldsymbol{X}$, converges weakly to the multivariate
$N(0,\bigoplus\limits_{i=1}^{a}\kappa_{i}^{-1}\boldsymbol{V}_{i})$-distribution
in probability.

Thus, the distributions of $\sqrt{n}\boldsymbol{C}\boldsymbol{\hat{p}}^{*}$
and $\sqrt{n}\boldsymbol{C}(\boldsymbol{\hat{p}}-\boldsymbol{p})$
asymptotically coincide under the null hypothesis $H_{0}$.

######  Theorem 4.2.

Under Assumption $(1)$, for any choice $[-]\in\\{A,M,W\\}$, the conditional
distribution of $T_{[-]}^{*}$ converges weakly to the null distribution of
$T_{[-]}$ in probability for any choice of
$\boldsymbol{p}\in\mbox{$\mathbb{R}$}^{ad}$. In particular we have that

$\displaystyle\sup_{x\in\mbox{$\mathbb{R}$}}|P_{\boldsymbol{p}}(T_{[-]}^{*}\leq
x|\boldsymbol{X})-P_{H_{0}}(T_{[-]}\leq x)|\xrightarrow{\text{p}}0$

holds, provided that $\boldsymbol{V}>0$ or $v_{i}(j,j)>0$ for all
$i\in\\{1,...,a\\},j\in\\{1,...,d\\}$ is fulfilled in case of $T_{W}$ or
$T_{M}$, respectively.

Therefore, the corresponding wild bootstrap tests are given by
$\phi^{*}_{W}=\mathds{1}\\{T_{W}>c^{*}_{w}\\}$,
$\phi^{*}_{A}=\mathds{1}\\{T_{A}>c^{*}_{A}\\}$, and
$\phi^{*}_{M}=\mathds{1}\\{T_{M}>c^{*}_{M}\\}$, where $c^{*}_{w}$,
$c^{*}_{A}$, and $c^{*}_{M}$ denoting the $(1-\alpha)$ quantiles of the
conditional bootstrap distributions of $T_{W}^{*}$, $T_{A}^{*}$, and
$T_{M}^{*}$ given the data, respectively.

Theorem 4.2 implies that the wild bootstrap tests are of asymptotic level
$\alpha$ under the null hypothesis and consistent for any fixed alternative
$H_{1}:\boldsymbol{C}\boldsymbol{F}\neq\boldsymbol{0}$, i.e., they have
asymptotically power $1$. In addition, it follows from Jansen et al. (2003)
that they have the same local power under contiguous alternatives as their
original tests.

Figure 1: Type-I error of the asymptotic ANOVA-type test $T_{A}$ and the
bootstrapped tests $T_{W}^{*}$, $T_{A}^{*}$, and $T_{M}^{*}$ based on several
hypotheses of interest for varying time points $d\in\\{4,8\\}$. Each boxplot
summarizes the type-I error results from 96 different simulation scenarios for
this hypothesis. For all individual simulations, see the tables in the
supplement. Figure 2: Type-I error of the asymptotic ANOVA-type test $T_{A}$
and the bootstrapped tests $T_{W}^{*}$, $T_{A}^{*}$, and $T_{M}^{*}$ based on
missing rates $r=\\{10,30\\}$ for varying time points $d\in\\{4,8\\}$. Each
boxplot summarizes the type-I error results from 144 different simulation
scenarios for this missing rate. For all individual simulations, see the
tables in the supplement.

(a) Normal (b) Chisquare

Figure 3:  Type-I error simulation results ($\alpha=0.05$) of the tests
$T_{W}$ $({\color[rgb]{0,0,0}\textendash\textendash\textendash})$, $T_{A}$
$({\color[rgb]{0,0,1}\textendash\textendash\quad\textendash\textendash}),$
$T_{W}^{*}$ $({\color[rgb]{.5,0,.5}-\boldsymbol{\cdot}-})$, $T_{A}^{*}$
$({\color[rgb]{.75,.5,.25}\boldsymbol{\cdots}})$, and $T_{M}^{*}$
$({\color[rgb]{0,0,0}\textbf{-- - --}})$ under different covariance structures
with sample sizes $(n_{1},n_{2})=(15,15)$ and $d=4$ for varying percentages of
MCAR data $r\in\\{10\%,20\%,30\%,40\%,50\%,60\%\\}$ with observations
generated from (a) a Normal and (b) a $\chi^{2}_{15}$ distribution,
respectively.

##  5 Monte-Carlo simulations

The above results are valid for large sample sizes. To analyze the finite
sample behavior of the asymptotic quadratic tests and their wild bootstrap
counterparts described in Sections 3 and 4, we conduct extensive simulations.
As an assessment criteria, all procedures were studied with respect to their

  1. (i)

type-I error rate control at level $\alpha=5\%$ and their

  2. (ii)

power to detect deviations from the null hypothesis.

All simulations were operated by means of the R computing environment, version
3.2.3 R Core Team (2013) and each setting was based on $10,000$ simulation
runs and $B=999$ bootstrap runs. The algorithm for the computation of the
$p$-value of the wild bootstrap tests in any of the test statistics
$T\in\\{T_{A},T_{W},T_{M}\\}$ is as follows:

  1. 1.

For the given incomplete multi-factor data, calculate the observed test
statistic, say $T$.

  2. 2.

Compute the rank values $\boldsymbol{Z}_{ik}$.

  3. 3.

Using i.i.d Rademacher weights ${W}_{i1},...,{W}_{in_{i}}$, generate
bootstrapped rank values $\boldsymbol{Z}_{ik}^{*}=W_{ik}.\boldsymbol{Z}_{ik}$.

  4. 4.

Calculate the value of the test statistic for the bootstrapped sample $T^{*}$.

  5. 5.

Repeat the Steps $3$ and $4$ independently $B=999$ times and collect the
observed test statistic values in ${T^{*}_{b},b=1,.....,B}$.

  6. 6.

Finally, estimate the wild bootstrap $p$-value as
$P-value=\frac{\sum_{b=1}^{B}I(T^{*}_{b}>=T)}{B}$.

We considered a two-way layout design with $a=2$ independent groups and two
different time points $d\in\\{4,8\\}$ underlying discrete and continuous
distributions.  As in Konietschke et al. (2015) and Bathke et al. (2018) for
the case of complete observations, we investigated balanced situations with
sample size vectors $(n_{1},n_{2})\in\\{(5,5),(10,10),(20,20)\\}$, and
unbalanced situation with sample size vector $(n_{1},n_{2})=(10,20)$. In
particular, we studied three different kinds of hypotheses: the hypothesis of
no group effect
$H_{0}^{G}:(\boldsymbol{P}_{a}\otimes\frac{1}{d}\boldsymbol{1}_{d}^{T})$, no
time effect
$H_{0}^{T}:(\frac{1}{a}\boldsymbol{1}_{a}^{T}\otimes\boldsymbol{P}_{d})$, as
well as no interaction effect
$H_{0}^{GT}:(\boldsymbol{P}_{a}\otimes\boldsymbol{P}_{d})$.

For each scenario, we generated missingness within MCAR as well as MAR
frameworks as described below: For the MCAR mechanism, we simulated
$\boldsymbol{X}_{ik}=(\delta_{1ik}X_{i1k},...,\delta_{idk}X_{idk}),\quad
k=1,...,n_{i}$, for independent Bernoulli distributed $\delta_{i1k}\sim B(r)$
and a zero entry was interpreted as a missing observation. The missing
probability was chosen from $r\in\\{10\%,30\%\\}$.

For the MAR mechanism, we considered several pairs of features
$\\{X_{obs},X_{miss}\\}$. For each pair, there is a determining feature
$X_{obs}$ that determines the missing pattern of its corresponding $X_{miss}$
(Santos et al., 2019). Thus, for $d=4$, we define the pairs
$\\{X_{i1k},X_{i2k}\\}$ and $\\{X_{i3k},X_{i4k}\\}$. Whereas, for $d=8$, we
define the following pairs
$\\{X_{i1k},X_{i2k}\\},\\{X_{i1k},X_{i3k}\\},\\{X_{i6k},X_{i7k}\\},$ and
$\\{X_{i6k},X_{i8k}\\}$. Two different MAR scenarios are investigated - MAR1
and MAR2-. In MAR1, for each pair, we divided $\mathbf{X}$ into three groups
based on their $X_{obs}$ values corresponding to a $2\sigma-$rule: the first
group is given by
$\\{\boldsymbol{X}_{ik}:X_{obs}\in(-\infty,-2\sigma_{1}),k=1,..,n_{i}\\}$, the
second by
$\\{\boldsymbol{X}_{ik}:X_{obs}\in(-2\sigma_{1},2\sigma_{1}),k=1,..,n_{i}\\}$
and the last group by
$\\{\boldsymbol{X}_{ik}:X_{obs}\in(2\sigma_{1},\infty),k=1,..,n_{i}\\}$, where
$\sigma_{1}$ is the variance of $X_{obs}$. Then, we randomly inserted missing
values on $X_{miss}$ based on the following missing percentages: $15\%$ for
group one and three and $30\%$ for the second group.

For generating MAR2 scenario, we considered the median of each $X_{obs}$ to
define the missing pattern of $X_{miss}$ (Zhu et al., 2012; Pan et al., 2015).
For each pair from above, two groups were defined: the first one is
$\\{\boldsymbol{X}_{ik}:X_{obs}\in(-\infty,median(X_{obs})],k=1,..,n_{i}\\}$,
and the second is
$\\{\boldsymbol{X}_{ik}:X_{obs}\in(median(X_{obs}),\infty),k=1,..,n_{i}\\}$.
Then, missing values were inserted on $X_{miss}$ as follows: $10\%$ for group
one and $30\%$ for the second group.

###  5.1 Continuous data

To investigate the type-I error control of the suggested methods, data samples
were generated from the model

$\boldsymbol{X}_{ik}\sim\boldsymbol{F}_{i}(\boldsymbol{\mu}_{0},\boldsymbol{\Sigma}_{i}),\quad
i=1,2,\quad k=1,...,n_{i},$

where, $\boldsymbol{F}_{i}(\boldsymbol{\mu}_{0},\boldsymbol{\Sigma}_{i})$
represents a multivariate distribution with expectation vector
$\boldsymbol{\mu}_{0}$ and covariance matrix $\boldsymbol{\Sigma}_{i}$.
Marginal data distributions were generated from different symmetric
distributions (normal, double exponential) as well as skewed distributions
(lognormal, $\chi^{2}_{15}$) (similar to Pauly et al. (2015a)). We used normal
copulas to generate the dependency structures of the repeated measurements
using the R package copula (Yan et al., 2007). For the covariance matrix
$\boldsymbol{\Sigma}_{i}$, we investigated the three following covariance
structures

(AR) Setting 1:

$\boldsymbol{\Sigma}_{i}=(\rho^{|l-j|})_{l,j\leq d},\rho=0.6$ for $i=1,2$,

(CS) Setting 2:

$\boldsymbol{\Sigma}_{i}=\boldsymbol{I}_{t}$, for $i=1,2,$

(TP) Setting 3:

$\boldsymbol{\Sigma}_{i}=(d-|l-j|)_{l,j\leq d},$ for $i=1,2$,

representing an autoregressive structure (Setting $1$), compound symmetry
pattern (Setting $2$) and a linear Toeplitz covariance structure (Setting
$3$). These covariance settings were inspired by the simulations studies in
Konietschke et al. (2015); Friedrich et al. (2017b) and Umlauft et al. (2017).

Type-I Error Results. The type-I error simulation results of the studied
procedures for testing the hypotheses of no time effect $H_{0}^{T}$ and no
group $\times$ time interaction $H_{0}^{GT}$ under the MCAR and MAR frameworks
are shown in Tables S.1 - S.8 (MCAR framework) and Tables S.9 - S.16 (MAR
framework) in the supplement. It can be readily seen that the suggested
bootstrap approaches based upon $T_{W}^{*},T_{A}^{*}$ and $T_{M}^{*}$ tend to
result in quite accurate type-I error rate control for most hypotheses under
symmetric as well as skewed distributions and under MCAR and MAR mechanisms.
The type-I error control is surprisingly not affected by less stringent
missing mechanisms and the bootstrap tests are robust under fairly large
amounts of missing observations. In addition, the data dependency structures
do not affect the quality of the approximations. Moreover, there was no
significant impact of the sample sizes balanced or unbalancedness in the
control of type-I error. On the other hand, the asymptotic ANOVA-type test
$T_{A}$ also shows a quite accurate type-I error control for large sample
sizes. However, under small sample sizes, $T_{A}$ tends to be sensitive to the
missing rates. In particular, it exhibits a liberal behavior for larger
missing rates. In contrast, the asymptotic Wald test $T_{W}$ shows an
extremely liberal behavior in all considered situations and under all
investigated missing mechanisms. A closer look at the type-I error simulation
results for all considered settings under the MCAR framework is provided.
Compact boxplots factorized in terms of hypotheses of interest or missing
rates are given in Figure 1 and 2, respectively.  They are based upon
different sample size settings (4) x different distributions (4) x different
missing rates (2) x different covariance structures (3) x different time
points (2) x different hypotheses (3). Due to the extreme liberal behavior of
the asymptotic Wald test, it has been removed from those plots in order to get
a better view of the remaining tests behavior. It can be clearly seen that the
asymptotic ANOVA-type test has the less accurate control among all other
considered tests. The hypothesis of interest and number of time points
obviously impact the quality of the approximations of the asymptotic ANOVA-
type test. Furthermore, the type-I error control behavior of the bootstrapped
MATS depends on the hypothesis of interest. Consequently, the bootstrapped
Wald and ANOVA-type tests are recommended over all considered methods.

In order to cover the effect of increasing missing rates, we additionally
studied type-I error control for $a=2$ groups, $d=4$ time points,
$(n_{1}=15,n_{2}=15)$ sample sizes with
$r\in\\{10\%,20\%,30\%,40\%,50\%,60\%\\}$ covering missingness in observations
ranging from $10\%$ to $60\%$. Figure 3 and Figures S.1 - S.2 in the
supplement summarize type-I error rate control for these settings under
symmetric and asymmetric distributions. The results indicate that the
asymptotic Wald test $T_{W}$ tends to be liberal in all considered situations.
In particular, It is extremely liberal when testing the hypotheses of no time
effect ($H_{0}^{T}$) and no interaction effect ($H_{0}^{GT}$). In contrast,
the asymptotic ANOVA-type test $T_{A}$ tends to be sensitive to missing rates
and hypothesis type. In particular, it exhibits an accurate or liberal
behavior for small or large missing rates, respectively. Moreover, it shows a
quite constant liberal behavior when testing the hypothesis of no group effect
($H_{0}^{G}$). Its behaviour appears to be independent of the covariance
pattern. In contrast, the suggested bootstrap approaches tend to control
type-I error rate more accurate over the range of missing rates $r$ for almost
all settings.

Figure 4:  Type-I error simulation results ($\alpha=0.05$) of the tests
$T_{W}$ $({\color[rgb]{0,0,0}\textendash\textendash\textendash})$, $T_{A}$
$({\color[rgb]{0,0,1}\textendash\textendash\quad\textendash\textendash}),$
$T_{W}^{*}$ $({\color[rgb]{.5,0,.5}-\boldsymbol{\cdot}-})$, $T_{A}^{*}$
$({\color[rgb]{.75,.5,.25}\boldsymbol{\cdots}})$, and $T_{M}^{*}$
$({\color[rgb]{0,0,0}\textbf{-- - --}})$ for ordinal data under MCAR framework
with sample sizes $(n_{1},n_{2})=\\{(10,10),(10,20),(20,10),(20,20)\\}$ and
$d=4$.

###  5.2 Ordinal data

In order to address all the goals outlined above, we simulated ordinal data.
The observations were simulated similar to Brunner and Langer (2000) as
follows:

$\displaystyle\boldsymbol{X}_{ijk}=int\big{(}4\frac{cZ_{ik}+Y_{ijk}}{c+1}\big{)}+1,\quad
i=1,...,a,\quad k=1,...,n_{i},\quad j=1,...,d,$

where $Z_{ik}$ and $Y_{ijk}$ are independently uniformly distributed in the
interval $[0,1]$, $c>0$ is a constant and $int(x)$ indicates the integer part
of $x$. The elements of $\boldsymbol{X}_{ijk}$ take values between $1$ and
$4$. The correlation between $\boldsymbol{X}_{ij^{`}k}$ and
$\boldsymbol{X}_{ijk}$ is determined by the choice of the constant $c$. In our
simulation study, we considered $c=1$ assuring a compound symmetric covariance
structure. We considered $a=2$ groups and $d\in\\{4,8\\}$ dimensions, and the
same sample sizes as in the continuous data settings.

The type-I error results of the considered methods under MCAR and MAR
frameworks are summarized in Figure 4 and Figures S.3 - S.5 in the supplement.
It can be seen that the asymptotic Wald test is too liberal in most situations
and under all considered hypotheses. The ANOVA-type test of Brunner et al.
(1999) shows much better behavior. In contrast, the type-I error control for
the bootstrap-based procedures is the best, particularly for the bootstrapped
Wald and ANOVA-type tests which are also less effected by an increased missing
rate or strict MAR assumptions.

Figure 5: Power simulation results of the tests $T_{W}$
$({\color[rgb]{0,0,0}\textbf{\textendash\textendash\textendash}})$, $T_{A}$
$({\color[rgb]{0,0,1}\textbf{
\textendash\textendash\quad\textendash\textendash}})$ , $T_{W}^{*}$
$({\color[rgb]{.5,0,.5}\textbf{--}\boldsymbol{\cdot}\textbf{--}})$,
$T_{A}^{*}$ $({\color[rgb]{.75,.5,.25}\boldsymbol{\cdots}})$, and $T_{M}^{*}$
$({\color[rgb]{0,0,0}\textbf{-- - --}})$ under different covariance structures
with sample size $n=15$ and $d=4$ under alternative $1$, under the MCAR
framework with missing rate $r=30\%$ with observations generated from a Normal
(upper row) and a Double exponential (bottom row) distribution, respectively.
Figure 6: Power simulation results of the tests $T_{W}$
$({\color[rgb]{0,0,0}\textbf{\textendash\textendash\textendash}})$, $T_{A}$
$({\color[rgb]{0,0,1}\textbf{
\textendash\textendash\quad\textendash\textendash}})$ , $T_{W}^{*}$
$({\color[rgb]{.5,0,.5}\textbf{--}\boldsymbol{\cdot}\textbf{--}})$,
$T_{A}^{*}$ $({\color[rgb]{.75,.5,.25}\boldsymbol{\cdots}})$, and $T_{M}^{*}$
$({\color[rgb]{0,0,0}\textbf{-- - --}})$ under different covariance structures
with sample size $n=15$ and $d=4$ under alternative $1$, under the MCAR
framework with missing rate $r=30\%$ with observations generated from a
$\chi^{2}_{15}$ (upper row) and a Lognormal (bottom row) distribution,
respectively.

###  5.3 Power

In order to assess the empirical power of all studied methods, we considered a
one-sample repeated measures design with $d=4$ repeated measures, sample size
$n=15$, and covariance structures as given in Settings $1-3$ for various
distributions. Data was generated by

$\displaystyle\boldsymbol{X}_{1k}\sim\boldsymbol{F}_{1}(\boldsymbol{\mu}_{0},\boldsymbol{\Sigma}_{1})+\boldsymbol{\mu}_{1},\quad
k=1,...,15,$

where, we were interested in detecting two specific alternatives

  * •

Alternative 1: $\boldsymbol{\mu}_{1}=(0,0,\zeta,\zeta)$,

  * •

Alternative 2: $\boldsymbol{\mu}_{1}=(0,0,0,\zeta)$,

for varying shift parameter $\zeta\in\\{0,0.5,1,1.5,2,2.5,3\\}$.

The power analysis results of the considered methods under the MCAR framework
for several distributions, involving various covariance settings for detecting
Alternative $1$ are summarized in Figures 5 \- 6 (missing rate $r=30\%$) and
Figures S.6 - S.7 in the supplement ($r=10\%$). The simulation results for
investigating Alternative $2$ are displayed in Figures S.8 - S.11 in the
supplement. The power analysis results of the considered methods under the
respective MAR framework are summarized in Figures S.12 - S.15 (MAR1 scenario)
and Figures S.16 - S.19 (MAR2 scenario).

As the asymptotic Wald test is too liberal compared to the other studied
methods, its power function is larger. Moreover, the bootstrapped Wald test
exhibits the lowest power behavior, while the bootstrapped MATS and ANOVA-type
have a quite similar power behavior as the Brunner et al. (1999) ANOVA-type
test. The differences between the procedures is less pronounced under the MAR
framework.

To sum up, we recommend the bootstrap ANOVA-type statistic. It exhibits the
overall best type-I-error control combined with a good power behaviour and
needs the less stringent assumptions for application.

(a) Side effect (b) Therapeutic effect

Figure 7: Frequencies of the side effect and therapeutic effect scores
observed in the fluvoxamine trial.

##  6 Application to empirical data

###  6.1 The fluvoxamine trial

In this section, we re-examine a clinical fluvoxamine study. It has already
been analyzed by Molenberghs and Lesaffre (1994); Molenberghs et al. (1997);
Van Steen et al. (2001); Jansen et al. (2003); Molenberghs and Verbeke (2004),
and Molenberghs and Kenward (2007). The study has been performed to establish
the profile of fluvoxamine in ambulatory clinical psychiatric operations.
Hereby, 315 patients who suffer from depression, panic disorder, and/or
obsessive-compulsive disorder were scored every two weeks over six weeks of
treatment (d=3). At each clinical visit, scores for both side effect and
therapeutic effect scales which are based upon about 20 psychiatric symptoms
were recorded. The side effect scale ranges from $1$ to $4$. The lower the
score, the better the clinical record. For example, score $1$ stands for ”no
side effect” while score $4$ indicates that ”the side effect surpasses the
therapeutic effect”. Similarly, the therapy effect is a $4$ category ordinal
scale: (1) ”no improvement or worsening”; (2) ”minimal improvement, not
changing functionality”; (3) ”moderate improvement, partial disappearance of
symptoms”; and (4) ”important improvement, almost disappearance of symptoms”.
The higher the score, the better the patients health.

Several patients missed the recording of their measurements in some sessions
which leaded to a large amount of missing values. A closer look to our data
shows that, from the total of 315 initially recruited patients, 14 patients
dropped off, 31 patients were scored on the first session only, and 44
patients on Session 1 and Session 3 only. Thus, only 224 patients have
complete observations. Two patients were excluded from the analysis due to
their non-monotone missing pattern. This leaves us with 299 patients. Waffle
plots representing the distributions of the side effect and therapeutic effect
among the three sessions are shown in Figure 7. We aim to test the hypotheses
whether side effect or therapeutic effect scores are significantly different
between the three sessions for patients with psychiatric disorder. To this
end, we applied all considered testing methods; asymptotic Wald and ANOVA-type
tests ($T_{W},T_{A}$), and the bootstrap procedures ($T_{W}^{*}$, $T_{A}^{*}$,
$T_{M}^{*}$) to detect the null hypothesis
$H_{0}^{T}:\\{\boldsymbol{C}\boldsymbol{F}=\boldsymbol{0}\\}$. The results are
summarized in Table 1. It can be seen that all tests indicate a significant
difference between the therapeutic effect scores as well as the side effect
scores of the three sessions (two-sided p-value $<0.00001$). Moreover, we
conclude that the clinical outcome of the patients significantly improves
after three sessions. These findings coincide with that in Molenberghs and
Kenward (2007).

Table 1: P-values of the fluvoxamine study. Effect | $T_{W}$ | $T_{A}$ | $T_{W}^{*}$ | $T_{A}^{*}$ | $T_{M}^{*}$
---|---|---|---|---|---
Therapeutic effect | $2.3e^{-79}$ | $2.2e^{-86}$ | 0 | 0 | 0
Side effect | $2.3e^{-10}$ | $4.1e^{-12}$ | 0 | 0 | 0

###  6.2 The skin disorder trial

Here, we study data from a randomized, multi-center, parallel group study for
treating a skin condition. Treating skin conditions can be tough, thus the
goal of the study was to assess the severe rate of the skin condition over
time and to compare the efficiency and safety of two continuous therapy
treatments drug and placebo. Patients were randomly assigned to drug or
placebo therapy treatment. Prior to treatment, patients were assessed to
determine the initial severity of the skin condition (moderate or severe). At
three follow-up visits, the treatment outcome was measured according to a
five-point ordinal response scale that assess the extent of improvement (1 =
rapidly improving, 2 = slowly improving, 3 = stable, 4 = slowly worsening, 5 =
rapidly worsening). The study consists of 88 and 84 subjects allocated to the
active treatment group and the placebo group, respectively. And, the
proportion of missing observations is around $30\%$. The distribution of
patients improvement across the treatment groups and the follow-up visits is
displayed in Figure 8. The study is described in full detail in Landis et al.
(1988) and is published in Davis (2002).

Figure 8: Frequencies of observed patients treatment outcome in the skin
disorder trial.

Similar to the fluvoxamine study above, we applied all asymptotic and
bootstrap procedures to infer the following null hypotheses: ”no group
effect”, ”no time effect ” and ”no group $\times$ time effect”. The results
are summarized in Table 2: All approaches reject the null hypothesis of no
group effect, the null hypothesis of no time effect and the null hypothesis of
no treatment group $\times$ time interaction. This implies that the clinical
outcome of the patients significantly improves with time and this progression
is significantly different between the two treatment groups drug and placebo.
Therefore, the data is further analyzed and split by the factor initial
severity and the analysis is replicated separately for each baseline severity
level (moderate or severe). The results are provided in Table 3.

It can be seen from Table 3 that all approaches detect a significant group
effect as well as a significant time effect in both moderate and severe
groups. In contrast, a significant group $\times$ time interaction effect
arises only in the moderate severity group.This indicates that the change in
clinical outcomes of patients of moderate severity varies over time depending
on treatment group membership. All five approaches share the previous
findings.

Table 2: P-values of the skin disorder trial. Hypothesis | $T_{W}$ | $T_{A}$ | $T_{W}^{*}$ | $T_{A}^{*}$ | $T_{M}^{*}$
---|---|---|---|---|---
Group | 0 | 0 | 0 | 0 | 0
Time | 0 | 0 | 0 | 0 | 0
Group $\times$ Time | $0.032$ | $0.016$ | $0.026$ | $0.01$ | $0.009$
Table 3: P-values of the split skin disorder trial based on initial severity. Hypothesis | Initial severity = Moderate |  | Initial severity = Severe
---|---|---|---
| $T_{W}$ | $T_{A}$ | $T_{W}^{*}$ | $T_{A}^{*}$ | $T_{M}^{*}$ |  | $T_{W}$ | $T_{A}$ | $T_{W}^{*}$ | $T_{A}^{*}$ | $T_{M}^{*}$
Group | 0 | 0 | 0 | 0 | 0 |  | 0 | 0 | 0 | 0 | 0
Time | $0.003$ | $0.004$ | $0.003$ | $0.002$ | $0.002$ |  | 0 | 0 | 0 | 0 | 0
Group $\times$ Time | $0.041$ | $0.032$ | $0.043$ | $0.027$ | $0.029$ |  | $0.423$ | $0.318$ | $0.447$ | $0.338$ | $0.335$

###  6.3 Simulations for data with similar attributes to the trials data

It was also interesting to discover the type-I error rate control of the tests
under similar attributes to the data sets of the fluvoxamine and skin disorder
clinical trials. The data sets reflect large sample sizes and a small or
moderate amount of missing values from either a one sample repeated
measurements design or a two way layout design. Simulation results for the
type-I error rate of the studied procedures for $(n=299,d=3)$ and
($n_{1}=88,n_{2}=84,d=3$) sample sizes are presented in Table S.21 and Table
S.22 in the supplement, respectively. It can be seen that most tests are
robust under both settings and control type-I error rate accurately. Only the
Wald type tests exhibit a liberal behaviour for some of the lognormal
settings.

##  7 Summary

Multi-group repeated measures design with ordinal or skewed observations are
quite common. If the observations are additionally subject to missing values
existing methods for testing null hypotheses in terms of distribution function
may be either liberal (Wald-type statistics) or run (asymptotically) on a
wrong type-I-error level (ANOVA-type statistic, (Domhof et al., 2002)). To
this end we investigated three alternatives based upon resampling. We proved
their asymptotic validity and analyzed their small sample behaviour regarding
type-I-error control and power in extensive simulations. Under all of the $5$
considered methods, an ANOVA-type statistic with critical values calculated by
means of a Wild bootstrap approach exhibits the best behaviour and is
recommended.

In the future, we will include the present methodology into the R package
MissPair (Amro et al., 2019a). Moreover, we plan to extend our investigations
to general MANOVA settings, e.g., extending the results of Dobler et al.
(2019a) to the situation with missing values.

## Acknowledgement

Lubna Amro and Markus Pauly acknowledge the support of the German Research
Foundation (DFG) (No. PA 2409/3-2). Frank Konietschke’s work was also
supported by the German Research Foundation (DFG) (No. KO 4680/3-2).

## References

  * Akritas [2011] Michael G. Akritas.  _Nonparametric Models for ANOVA and ANCOVA Designs_ , pages 964–968.  Springer Berlin Heidelberg, Berlin, Heidelberg, 2011.  ISBN 978-3-642-04898-2.  10.1007/978-3-642-04898-2_415.
  * Akritas and Arnold [1994] Michael G Akritas and Steven F Arnold.  Fully nonparametric hypotheses for factorial designs i: Multivariate repeated measures designs.  _Journal of the American Statistical Association_ , 89(425):336–343, 1994.
  * Akritas and Brunner [1997] Michael G Akritas and Edgar Brunner.  A unified approach to rank tests for mixed models.  _Journal of Statistical Planning and Inference_ , 61(2):249–277, 1997.
  * Akritas et al. [2002] Michael G Akritas, Jouni Kuha, and D Wayne Osgood.  A nonparametric approach to matched pairs with missing data.  _Sociological methods & research_, 30(3):425–454, 2002.
  * Amro et al. [2019a] Lubna Amro, Frank Konietschke, and Markus Pauly.  R package misspair.  https://github.com/lubnaamro/MissPair, 2019a.
  * Amro et al. [2019b] Lubna Amro, Markus Pauly, and Burim Ramosaj.  Asymptotic based bootstrap approach for matched pairs with missingness in a single-arm.  _arXiv preprint arXiv:1912.04902_ , 2019b.
  * Arnau et al. [2012] Jaume Arnau, Roser Bono, María J Blanca, and Rebecca Bendayan.  Using the linear mixed model to analyze nonnormal data distributions in longitudinal designs.  _Behavior research methods_ , 44(4):1224–1238, 2012.
  * Bartlett [1939] Maurice S Bartlett.  A note on tests of significance in multivariate analysis.  In _Mathematical Proceedings of the Cambridge Philosophical Society_ , volume 35, pages 180–185. Cambridge University Press, 1939.
  * Bathke et al. [2018] Arne C Bathke, Sarah Friedrich, Markus Pauly, Frank Konietschke, Wolfgang Staffen, Nicolas Strobl, and Yvonne Höller.  Testing mean differences among groups: multivariate and repeated measures analysis with minimal assumptions.  _Multivariate behavioral research_ , 53(3):348–359, 2018.
  * Biederman et al. [2008] Lori A Biederman, Thomas W Boutton, and Steven G Whisenant.  Nematode community development early in ecological restoration: The role of organic amendments.  _Soil Biology and Biochemistry_ , 40(9):2366–2374, 2008.
  * Bourlier et al. [2008] Virginie Bourlier, Alexia Zakaroff-Girard, Alexandra Miranville, Sandra De Barros, Marie Maumus, Coralie Sengenès, Jean Galitzky, Max Lafontan, Fredrik Karpe, KN Frayn, et al.  Remodeling phenotype of human subcutaneous adipose tissue macrophages.  _Circulation_ , 117(6):806, 2008.
  * Box [1954] George EP Box.  Some theorems on quadratic forms applied in the study of analysis of variance problems, i. effect of inequality of variance in the one-way classification.  _The annals of mathematical statistics_ , 25(2):290–302, 1954.
  * Brombin et al. [2013] Chiara Brombin, Edoardo Midena, and Luigi Salmaso.  Robust non-parametric tests for complex-repeated measures problems in ophthalmology.  _Statistical methods in medical research_ , 22(6):643–660, 2013.
  * Brunner [2001] Edgar Brunner.  Asymptotic and approximate analysis of repeated measures designs under heteroscedasticity.  _Mathematical Statistics with Applications in Biometry_ , page 2, 2001\.
  * Brunner and Langer [2000] Edgar Brunner and Frank Langer.  Nonparametric analysis of ordered categorical data in designs with longitudinal observations and small sample sizes.  _Biometrical Journal: Journal of Mathematical Methods in Biosciences_ , 42(6):663–675, 2000.
  * Brunner and Puri [2001] Edgar Brunner and Madan L Puri.  Nonparametric methods in factorial designs.  _Statistical papers_ , 42(1):1–52, 2001.
  * Brunner et al. [1999] Edgar Brunner, Ulrich Munzel, Madan L Puri, et al.  Rank-score tests in factorial designs with repeated measures.  _Journal of Multivariate Analysis_ , 70(2):286–317, 1999.
  * Brunner et al. [2017] Edgar Brunner, Frank Konietschke, Markus Pauly, and Madan L Puri.  Rank-based procedures in factorial designs: Hypotheses about non-parametric treatment effects.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 79(5):1463–1485, 2017.
  * Brunner et al. [2018] Edgar Brunner, Arne C Bathke, and Frank Konietschke.  _Rank and Pseudo-Rank Procedures for Independent Observations in Factorial Designs_.  Springer, 2018.
  * Davidson and Flachaire [2008] Russell Davidson and Emmanuel Flachaire.  The wild bootstrap, tamed at last.  _Journal of Econometrics_ , 146(1):162–169, 2008\.
  * Davis [2002] Charles S Davis.  _Statistical methods for the analysis of repeated measurements_.  Springer Science & Business Media, 2002.
  * Dobler et al. [2019a] Dennis Dobler, Sarah Friedrich, and Markus Pauly.  Nonparametric manova in meaningful effects.  _Annals of the Institute of Statistical Mathematics_ , pages 1–26, 2019a.
  * Dobler et al. [2019b] Dennis Dobler, Markus Pauly, and ThomasH Scheike.  Confidence bands for multiplicative hazards models: Flexible resampling approaches.  _Biometrics_ , 75(3):906–916, 2019b.
  * Domhof et al. [2002] Sebastian Domhof, Edgar Brunner, and D Wayne Osgood.  Rank procedures for repeated measures with missing values.  _Sociological methods & research_, 30(3):367–393, 2002.
  * Erceg-Hurn and Mirosevich [2008] David M Erceg-Hurn and Vikki M Mirosevich.  Modern robust statistical methods: an easy way to maximize the accuracy and power of your research.  _American Psychologist_ , 63(7):591, 2008.
  * Fitzmaurice et al. [2012] Garrett M Fitzmaurice, Nan M Laird, and James H Ware.  _Applied longitudinal analysis_ , volume 998.  John Wiley & Sons, 2012.
  * Fong et al. [2018] Youyi Fong, Ying Huang, Maria P Lemos, and M Juliana Mcelrath.  Rank-based two-sample tests for paired data with missing values.  _Biostatistics_ , 19(3):281–294, 2018.
  * Friedrich and Pauly [2018] Sarah Friedrich and Markus Pauly.  Mats: Inference for potentially singular and heteroscedastic manova.  _Journal of Multivariate Analysis_ , 165:166–179, 2018\.
  * Friedrich et al. [2017a] Sarah Friedrich, Edgar Brunner, and Markus Pauly.  Permuting longitudinal data in spite of the dependencies.  _Journal of Multivariate Analysis_ , 153:255–265, 2017a.
  * Friedrich et al. [2017b] Sarah Friedrich, Frank Konietschke, and Markus Pauly.  A wild bootstrap approach for nonparametric repeated measurements.  _Computational Statistics & Data Analysis_, 113:38–52, 2017b.
  * Gao [2007] Xin Gao.  A nonparametric procedure for the two-factor mixed model with missing data.  _Biometrical journal_ , 49(5):774–788, 2007.
  * Jansen et al. [2003] Ivy Jansen, Geert Molenberghs, Marc Aerts, Herbert Thijs, and Kristel Van Steen.  A local influence approach applied to binary data from a psychiatric study.  _Biometrics_ , 59(2):410–419, 2003.
  * Johnson and Wichern [2007] Richard Arnold Johnson and Dean W Wichern.  _Applied multivariate statistical analysis_ , volume 6.  Prentice hall Upper Saddle River, NJ, 2007.
  * Konietschke et al. [2012] Frank Konietschke, Solomon W Harrar, Katharina Lange, and Edgar Brunner.  Ranking procedures for matched pairs with missing data—asymptotic theory and a small sample approximation.  _Computational Statistics & Data Analysis_, 56(5):1090–1102, 2012.
  * Konietschke et al. [2015] Frank Konietschke, Arne C Bathke, Solomon W Harrar, and Markus Pauly.  Parametric and nonparametric bootstrap methods for general manova.  _Journal of Multivariate Analysis_ , 140:291–301, 2015\.
  * Krishnamoorthy and Lu [2010] K Krishnamoorthy and Fei Lu.  A parametric bootstrap solution to the manova under heteroscedasticity.  _Journal of Statistical Computation and Simulation_ , 80(8):873–887, 2010.
  * Landis et al. [1988] J Richard Landis, Michael E Miller, Charles S Davis, and Gary G Koch.  Some general methods for the analysis of categorical data in longitudinal studies.  _Statistics in medicine_ , 7(1-2):109–137, 1988\.
  * Lawley [1938] DN Lawley.  A generalization of fisher’s z test.  _Biometrika_ , 30(1/2):180–187, 1938.
  * Lekberg et al. [2018] Ylva Lekberg, James D Bever, Rebecca A Bunn, Ragan M Callaway, Miranda M Hart, Stephanie N Kivlin, John Klironomos, Beau G Larkin, John L Maron, Kurt O Reinhart, et al.  Relative importance of competition and plant–soil feedback, their synergy, context dependency and implications for coexistence.  _Ecology letters_ , 21(8):1268–1281, 2018.
  * Mammen [1992] Enno Mammen.  Bootstrap, wild bootstrap, and asymptotic normality.  _Probability Theory and Related Fields_ , 93(4):439–455, 1992.
  * Martinussen and Scheike [2007] Torben Martinussen and Thomas H Scheike.  _Dynamic regression models for survival data_.  Springer Science & Business Media, 2007.
  * Micceri [1989] Theodore Micceri.  The unicorn, the normal curve, and other improbable creatures.  _Psychological bulletin_ , 105(1):156, 1989.
  * Molenberghs and Kenward [2007] Geert Molenberghs and Michael Kenward.  _Missing data in clinical studies_ , volume 61.  John Wiley & Sons, 2007.
  * Molenberghs and Lesaffre [1994] Geert Molenberghs and Emmanuel Lesaffre.  Marginal modeling of correlated ordinal data using a multivariate plackett distribution.  _Journal of the American Statistical Association_ , 89(426):633–644, 1994.
  * Molenberghs and Verbeke [2004] Geert Molenberghs and Geert Verbeke.  Meaningful statistical model formulations for repeated measures.  _Statistica Sinica_ , pages 989–1020, 2004.
  * Molenberghs et al. [1997] Geert Molenberghs, Michael G Kenward, and Emmanuel Lesaffre.  The analysis of longitudinal ordinal data with nonrandom drop-out.  _Biometrika_ , 84(1):33–44, 1997.
  * Munzel and Brunner [2000] Ullrich Munzel and Edgar Brunner.  Nonparametric methods in multivariate factorial designs.  _Journal of Statistical Planning and Inference_ , 88(1):117–132, 2000.
  * Nesnow et al. [1998] Stephen Nesnow, Marc J Mass, Jeffrey A Ross, Anthony J Galati, Guy R Lambert, Chris Gennings, Walter H Carter Jr, and Gary D Stoner.  Lung tumorigenic interactions in strain a/j mice of five environmental polycyclic aromatic hydrocarbons.  _Environmental Health Perspectives_ , 106(suppl 6):1337–1346, 1998.
  * Pan et al. [2015] Ruilin Pan, Tingsheng Yang, Jianhua Cao, Ke Lu, and Zhanchao Zhang.  Missing data imputation by k nearest neighbours based on grey relational structure and mutual information.  _Applied Intelligence_ , 43(3):614–632, 2015\.
  * Pauly et al. [2015a] Markus Pauly, Edgar Brunner, and Frank Konietschke.  Asymptotic permutation tests in general factorial designs.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 77(2):461–473, 2015a.
  * Pauly et al. [2015b] Markus Pauly, David Ellenberger, and Edgar Brunner.  Analysis of high-dimensional one group repeated measures designs.  _Statistics_ , 49(6):1243–1261, 2015b.
  * Pesarin and Salmaso [2012] Fortunato Pesarin and Luigi Salmaso.  A review and some new results on permutation testing for multivariate problems.  _Statistics and Computing_ , 22(2):639–646, 2012\.
  * Qian and Huang [2005] Hui-Rong Qian and Shuguang Huang.  Comparison of false discovery rate methods in identifying genes with differential expression.  _Genomics_ , 86(4):495–503, 2005.
  * R Core Team [2013] R Core Team.  _R: A Language and Environment for Statistical Computing_.  R Foundation for Statistical Computing, Vienna, Austria, 2013.  URL http://www.R-project.org/.
  * Santos et al. [2019] Miriam Seoane Santos, Ricardo Cardoso Pereira, Adriana Fonseca Costa, Jastin Pompeu Soares, João Santos, and Pedro Henriques Abreu.  Generating synthetic missing data: A review by missing mechanism.  _IEEE Access_ , 7:11651–11667, 2019.
  * Smaga [2017] Łukasz Smaga.  Bootstrap methods for multivariate hypothesis testing.  _Communications in Statistics-Simulation and Computation_ , 46(10):7654–7667, 2017.
  * Traux and Carkhuff [1965] CB Traux and Robert R Carkhuff.  Personality change in hospitalized mental patients during group psychotherapy as a function of the use of alternate sessions and vicarious therapy pretraining.  _Journal of Clinical Psychology_ , 1965.
  * Umlauft et al. [2017] Maria Umlauft, Marius Placzek, Frank Konietschke, and Markus Pauly.  Wild bootstrapping rank-based procedures: Multiple testing in nonparametric split-plot designs.  _arXiv preprint arXiv:1710.04532_ , 2017.
  * Vallejo et al. [2010] G Vallejo, MP Fernández, and Pablo Esteban Livacic-Rojas.  Analysis of unbalanced factorial designs with heteroscedastic data.  _Journal of Statistical Computation and Simulation_ , 80(1):75–88, 2010.
  * Van Steen et al. [2001] Kristel Van Steen, Geert Molenberghs, Geert Verbeke, and Herbert Thijs.  A local influence approach to sensitivity analysis of incomplete longitudinal ordinal data.  _Statistical Modelling_ , 1(2):125–142, 2001\.
  * Wildsmith et al. [2001] SE Wildsmith, GE Archer, AJ Winkley, PW Lane, and PJ Bugelski.  Maximization of signal derived from cdna microarrays.  _Biotechniques_ , 30(1):202–208, 2001.
  * Xu and Cui [2008] Jin Xu and Xinping Cui.  Robustified manova with applications in detecting differentially expressed genes from oligonucleotide arrays.  _Bioinformatics_ , 24(8):1056–1062, 2008.
  * Xu et al. [2013] Li-Wen Xu, Fang-Qin Yang, Shuang Qin, et al.  A parametric bootstrap approach for two-way anova in presence of possible interactions with unequal variances.  _Journal of Multivariate Analysis_ , 115:172–180, 2013\.
  * Yan et al. [2007] Jun Yan et al.  Enjoy the joy of copulas: with a package copula.  _Journal of Statistical Software_ , 21(4):1–21, 2007.
  * Zhu et al. [2012] Bing Zhu, Changzheng He, and Panos Liatsis.  A robust missing value imputation method for noisy data.  _Applied Intelligence_ , 36(1):61–74, 2012.