# A Global Model Approach to Robust Few-Shot SAR Automatic Target Recognition

Nathan Inkawhich
Air Force Research Laboratory

###### Abstract

In real-world scenarios, it may not always be possible to collect hundreds of
labeled samples per class for training deep learning-based SAR Automatic
Target Recognition (ATR) models. This work specifically tackles the few-shot
SAR ATR problem, where only a handful of labeled samples may be available to
support the task of interest. Our approach is composed of two stages. In the
first, a global representation model is trained via self-supervised learning
on a large pool of diverse and unlabeled SAR data. In the second stage, the
global model is used as a fixed feature extractor and a classifier is trained
to partition the feature space given the few-shot support samples, while
simultaneously being calibrated to detect anomalous inputs. Unlike competing
approaches which require a pristine labeled dataset for pretraining via meta-
learning, our approach learns highly transferable features from unlabeled data
that have little-to-no relation to the downstream task. We evaluate our method
in standard and extended MSTAR operating conditions and find it to achieve
high accuracy and robust out-of-distribution detection in many different few-
shot settings. Our results are particularly significant because they show the
merit of a global model approach to SAR ATR, which makes minimal assumptions,
and provides many axes for extendability.

###### Index Terms:

Few Shot Learning, Automatic Target Recognition, Out-of-Distribution
Detection, Deep Learning.

##  I Introduction

Most modern Synthetic Aperture Radar (SAR) Automatic Target Recognition (ATR)
algorithms use Deep Neural Networks (DNNs) to perform the underlying
recognition function [1]. While DNNs have shown to be very effective in ideal
laboratory-like conditions, they present several challenges when being
considered for deployment in “real-world” scenarios. One challenge is that
DNNs are data-hungry, and often utilize hundreds to thousands of training
samples per class to achieve state-of-the-art accuracy. Another key challenge
is that DNNs are notoriously susceptible to producing erroneous predictions on
anomalous/out-of-distribution (OOD) inputs (i.e., inputs that fall outside of
the categories in the training distribution). These challenges spark two main
concerns. First, what if we, as the model practitioners, cannot obtain
hundreds of labeled samples per class for training, and instead can only come
up with $\sim$5 (i.e., a few shots)? And second, in these few-shot settings,
does OOD detection become harder because the models are primarily seeking
generalization and have not been given enough information to learn highly
detailed/nuanced representations of the in-distribution (ID) classes?

Figure 1: Overview of our two-stage approach. In Stage 1 we train the feature
extractor on a pool of unlabeled data and in Stage 2 we learn a classifier on
top of the feature extractor for a given few shot learning problem.

Our goal in this paper is to address the stated concerns by developing an
accurate and robust SAR target classifier given only a few labeled samples per
class. We also strive to keep a minimal assumption set to prioritize
applicability and extendability. Our high level approach involves two stages
and is shown in Fig. 1. In the first stage, a global representation model is
trained on a diverse and unlabeled pool of SAR data such that it is capable of
extracting quality features from (nearly) any SAR chip. In Stage 2, for an
arbitrary downstream (M-way, N-shot) few-shot learning (FSL) task we train a
light-weight classifier using the global model as a fixed feature extractor.

A key part of our methodology is how we train the global model. Importantly,
we do not assume to have labeled data in this step, which eliminates the
option for a meta-learning approach [2, 3]. To learn quality representations
we instead aggregate publicly available data from recent papers and leverage a
modern unsupervised learning paradigm called Self-Supervised Learning (SSL).
Specifically, we use the SimCLR [4] method to represent the SSL class of
algorithms in this work. SSL has shown great promise for transfer and few-shot
learning in the natural imagery domain [5], but has not been studied widely in
the context of SAR ATR. Our hypothesis is that these modern algorithms can
learn a highly flexible / expressive / useful / transferable feature extractor
from unlabeled SAR data sourced from different sensors, imaging modes,
polarizations, resolutions, and target types. Further, the unlabeled data does
not necessarily have to represent the anticipated downstream conditions. This
is contrary to several lines of concurrent research whose goal is to generate
high-fidelity synthetic imagery as close to the downstream task conditions as
possible for supervised pre-training [6].

The second key part of our methodology is how we improve robustness to OOD
inputs. This is accomplished in the classifier training stage and given our
stated operating environment it is nearly a “free-lunch.” Practically, we can
use the SSL model’s pretraining dataset as an OOD dataset for outlier exposure
(OE) training [7, 8]. Intuitively, the classifier is taught to be accurate and
confident on labeled task data while being minimally confident on the OOD
data. This objective creates a confidence calibration effect which can be
exploited to detect novel OOD samples during deployment [9].

Using MSTAR [10] as our primary few-shot test environment, we find that our
methodology can generate highly accuracy models in both standard and extended
operating conditions. We also find that by OE training in Stage 2, our models
can reliably detect a spectrum of fine- to coarse-grained OOD types unseen
during training. Finally, we uncover and discuss an important trade-off
between OOD generalization and detection. Overall, we make the following
contributions:

  * •

We show a proof-of-concept that a highly transferable global SAR feature
extractor can be trained without labels on diverse data and used to great
effect in downstream FSL tasks;

  * •

We develop a robust FSL classification scheme that can reliably detect and
reject OOD inputs at test time without adding overhead or assumptions;

  * •

We uncover a critical trade-off between OOD detection and generalization that
directly motivates future work.

##  II Methodology

###  II-A Global Model Training

####  II-A1 Data Collection

Our global model approach is highly data driven. In order to learn a
generalized and transferable SAR feature extractor, intuition says we must
train it on data from different sensors, imaging modes, polarizations,
resolutions, target/scene types, etc. Thus, the first step in our methodology
is to aggregate data from the following public sources: SAR-Ships [11]; HRSID
[12]; FUSAR [13]; SSDD [14]; LS-SSDD [15]; Dual-Pol Ships [16]; SRSDD [17];
and CVDome [18]. After the necessary pre-processing (e.g., chipping from a
full frame), we are left with $\sim$100k unlabeled SAR chips which we refer to
as $\mathcal{D}_{pretrain}$. For reference, Fig. 2 shows a handful of samples
from this set. An important note is that $90\%$ of $\mathcal{D}_{pretrain}$
chips contain ship-like targets and none of the chips contain MSTAR land-
vehicle targets. Given that MSTAR will be our downstream task for FSL
evaluation, this sets up for a true test of generalization. However, in
practice if data resembling the anticipated downstream task were available it
could easily be leveraged here.

####  II-A2 Representation Learning

With $\mathcal{D}_{pretrain}$ curated we can execute Stage 1, which is to
train the global model feature extractor $f_{\phi}$. We use the powerful
SimCLR [4] contrastive learning algorithm, which trains by enforcing that two
“views” of the same image lie near each-other in a feature space, while views
of different images lie far away. Each view is created via aggressive
stochastic augmentations that preserve the underlying semantic features but
force the model to learn an informative feature set via minimizing the
contrastive objective. We specifically leverage the normalized temperature-
scaled cross-entropy loss (NT-Xent) from [4], and refer the reader to that
work for a more technical description of the algorithm and [19] for a PyTorch
implementation.

Perhaps the most unique part about the application of SimCLR (and many other
SSL algorithms) to the SAR ATR setting is the composition of augmentations
needed to create informative views. Since the SimCLR algorithm was developed
and tested in the natural imagery domain, the standard augmentation pipeline
contains transforms that are relevant to RGB images. However, several of the
key augmentations such as color jittering and random gray-scaling [19] do not
fit with our magnitude-only gray-scale SAR modality. Thus, we create a new
augmentation pipeline made of seven image-level transformations that work to
create translation invariance, and to manipulate the signal-to-noise ratio and
frequency content of the chips. Included are random resized cropping, random
horizontal flipping, element-wise power-scaling, Gaussian noising and
filtering, and linear re-scaling. See Supplemental A for more information
about our transformation pipeline. Lastly, we mention that we do not believe
our design choices made in this section are “optimal.” Important areas of
future work are to experiment with other SSL algorithms (see Supplemental C)
and to develop more domain-relevant augmentations that may lead to higher
quality learned features.

Figure 2: Samples from $\mathcal{D}_{pretrain}$, sourced from: SAR-Ships [11];
HRSID [12]; FUSAR [13]; SSDD [14]; LS-SSDD [15]; Dual-Pol Ships [16]; SRSDD
[17]; and CVDome [18]. TABLE I: MSTAR Standard Operating Condition Accuracy
Results (with 95% confidence interval)

| 2-way | 5-way | 10-way
---|---|---|---
Method | 1-shot | 5-shot | 10-shot | 20-shot | 25-shot | 1-shot | 5-shot | 10-shot | 20-shot | 25-shot | 1-shot | 5-shot | 10-shot | 20-shot | 25-shot
Scratch-AConv | 59.6$\pm$1.2 | 73.6$\pm$1.2 | 84.3$\pm$0.9 | 92.1$\pm$0.7 | 93.8$\pm$0.6 | 31.6$\pm$0.6 | 55.9$\pm$0.8 | 72.3$\pm$0.8 | 85.7$\pm$0.5 | 88.5$\pm$0.4 | 20.9$\pm$0.3 | 47.4$\pm$0.4 | 65.0$\pm$0.4 | 80.7$\pm$0.2 | 84.5$\pm$0.2
Scratch-RN18 | 55.8$\pm$1.2 | 74.9$\pm$1.4 | 84.1$\pm$1.1 | 89.6$\pm$0.8 | 92.1$\pm$0.6 | 30.6$\pm$0.7 | 53.7$\pm$0.8 | 67.0$\pm$0.7 | 79.3$\pm$0.6 | 83.2$\pm$0.5 | 20.7$\pm$0.3 | 43.4$\pm$0.4 | 58.9$\pm$0.4 | 74.0$\pm$0.4 | 78.5$\pm$0.3
1-16 |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
SimCLR+basic | 64.5$\pm$1.4 | 83.0$\pm$1.1 | 91.4$\pm$0.7 | 96.1$\pm$0.4 | 97.4$\pm$0.2 | 38.8$\pm$0.7 | 68.4$\pm$0.7 | 82.6$\pm$0.5 | 91.8$\pm$0.3 |  94.0$\pm$0.2 | 28.3$\pm$0.3 | 59.2$\pm$0.3 | 75.2$\pm$0.2 |  87.7$\pm$0.1 |  90.6$\pm$0.1
SimCLR+OE |  65.4$\pm$1.5 |  86.4$\pm$1.0 |  92.2$\pm$0.8 |  96.5$\pm$0.4 |  97.5$\pm$0.2 |  39.9$\pm$0.7 |  69.7$\pm$0.6 |  83.3$\pm$0.5 |  92.3$\pm$0.3 | 93.8$\pm$0.2 |  28.8$\pm$0.3 |  60.5$\pm$0.3 |  75.6$\pm$0.2 | 87.0$\pm$0.1 | 89.3$\pm$0.1

###  II-B Classifier Training

In Stage 2 (ref. Fig. 1) we use the global model $f_{\phi}$ as a fixed feature
extractor and learn a classifier $c_{\theta}$ for the given (M-way, N-shot)
FSL task, described by
$\mathcal{D}_{support}=\\{(x_{i},y_{i})\\}_{i=1\ldots(M\times N)}$. To be
clear, “ways” refers to the number of categories in the label-space and
“shots” is number of labeled training chips per category. Unlike related works
in the natural imagery domain that use simple linear regression or nearest
neighbor-style classifiers [5], we find it beneficial to utilize a 2-layer
neural network which has the flexibility to learn non-linear decision
boundaries.

We specifically consider two ways to train $c_{\theta}$. The first, called the
basic method, is to train under a vanilla supervised learning objective which
minimizes empirical risk on $\mathcal{D}_{support}$ samples only.
Mathematically, the basic training method is described as

$\min_{\theta}\mathop{\mathbb{E}}_{(x,y)\sim\mathcal{D}_{support}}\big{[}L(c_{\theta}(f_{\phi}(x)),y)\big{]},$
(1)

where $L$ is the cross-entropy loss between the classifier’s prediction and
the ground truth label $y$. The intuition for this training method is to learn
a classifier that produces accurate and confident predictions on the ID data.

While the basic method can achieve high ID accuracy, it gives no guidance to
the model for how to behave when an OOD input is encountered. Without
violating or adding any assumptions, we can do something clever to greatly
improve the classifier’s ability to handle OOD inputs. Specifically, we can
re-purpose the $\mathcal{D}_{pretrain}$ dataset from Stage 1 as an outlier
exposure (OE) set in Stage 2. The intuition for OE [7] is to teach the model a
confidence calibration on ID and OOD inputs – it should make accurate and
confident predictions on ID data and minimally confident predictions on OOD
data. Functionally, the minimum confidence state is when the predicted
probability is $\nicefrac{{1}}{{\\#classes}}$, so to achieve the desired
effect we set the training targets of the OE data to be a Uniform distribution
over the classes $\mathcal{U}_{\mathcal{C}}$, while the targets of the ID data
remain 1-hots. The complete OE objective is

$\min_{\theta}\mathop{\mathbb{E}}_{(x,y)\sim\mathcal{D}_{support}}\big{[}L(c_{\theta}(f_{\phi}(x)),y)\big{]}\\\
+\lambda\mathop{\mathbb{E}}_{\tilde{x}\sim\mathcal{D}_{pretrain}}\big{[}L(c_{\theta}(f_{\phi}(\tilde{x})),\mathcal{U}_{\mathcal{C}})\big{]},$
(2)

where $\lambda$ is a weighting factor between the ID- and OOD-focused terms
that we set to 0.5 here [7].

###  II-C OOD Detection

The final detail of our methodology is how we detect OOD samples during
deployment. To conceptualize the process, think about each test input $x$
being assigned a real-valued score $\mathcal{S}_{ID}(x)$ as a measure of its
ID-ness. During operation, if $\mathcal{S}_{ID}(x)\geq\beta_{thresh}$ (an
application dependent threshold), then the input would be considered ID and
the classifier would release its prediction over the set of known classes.
Else, if $\mathcal{S}_{ID}(x)<\beta_{thresh}$, the input is considered OOD and
the system abstains from releasing the prediction. In this work we use a
temperature-scaled Maximum Softmax Probability detector [20] to produce
$\mathcal{S}_{ID}(x)$ scores, described as

$\mathcal{S}_{ID}(x)=\max_{i\in\mathcal{C}}\frac{exp\big{(}c_{\theta}^{(i)}(f(x))/\tau\big{)}}{\sum_{j}^{\mathcal{|C|}}exp\big{(}c_{\theta}^{(j)}(f(x))/\tau\big{)}}.$
(3)

Here, $\tau$ is a temperature hyperparameter which we set to 100 via
validation. The intuition for this detection scheme is straightforward: ID
samples should be predicted by the classifier with higher confidence than OOD
samples. This fits naturally with the OE training goal, making it a logical
choice in our system. We note that OOD detection is currently a popular
research topic and future works may investigate alternative scoring
mechanisms.

##  III Experiments

To show the merit of our approach to the robust few-shot SAR ATR problem we
perform three primary experiments and one analysis. The first experiment (Sec.
III-A) shows results in MSTAR’s Standard Operating Condition (SOC) and the
second (Sec. III-B) investigates generalization performance in an Extended
Operating Condition (EOC). The third experiment (Sec. III-C) measures OOD
detection performance against a spectrum of granularities and the analysis
(Sec. III-D) discusses a critical trade-off between OOD detection and
generalization.

All of the experiments follow a similar setup. In Stage 1, we initialize the
feature extractor as a ResNet-18 (RN18) [21] backbone with output dimension
512, and the Projector network as a 3-layer neural net [19] with output
dimension 128. These components are trained with an ADAM optimizer for 200
epochs, with 1024 batch size, weight decay of 1e-4, and a learning rate of
3e-4 following a cosine decay schedule. For the SimCLR NT-Xent loss we use a
temperature of 0.01. In Stage 2, we discard the Projector network, fix the
weights of the RN18, and initialize our 2-layer neural net classifier. We
train the classifier for 500 iterations, with label smoothing, using an ADAM
optimizer with cosine decayed learning rate starting at 1e-3. On the data
side, we train with $64\times 64$ px crops and use Gaussian noise and random
flipping augmentations. Finally, the results in this document are averages
over 250 randomized runs.

###  III-A Standard Operating Conditions

To measure SOC performance on MSTAR we train the models on 17∘ elevation
imagery and test the accuracy of classifying 15∘ data [22]. As mentioned,
because we do not assume to have labeled pretraining data, meta-learning-based
FSL methods are inappropriate for comparison (also see Supplemental B). For
baselines we instead train supervised models from scratch on the
$\mathcal{D}_{support}$ samples only. We examine two architectures of varying
complexities: A-ConvNet [22] and RN18 (displayed as Scratch-{AConv,RN18}). The
MSTAR SOC results are shown in Table I for #-ways={2,5,10} and
#-shots={1,5,10,20,25}, covering a gamut of potential scenarios.

To start, we find several intuitive takeaways. First, accuracy increases with
the number of shots; and second, as the number of ways increases the models
need more shots to achieve the same performance. In all cases, our SimCLR-
based models significantly outperform the models trained from scratch, and in
most cases the OE classifier has a slight edge over the basic classifier. When
only given 10-shots, our average margin of improvement over the best baseline
is +9.8%. On full 10-class MSTAR, we reach 90% average accuracy at
$\sim$25-shots, which is a 6% improvement over the best baseline. Lastly, we
want to emphasize the implications of these results. The self-supervised
SimCLR backbone, which has been trained on mostly ships, is able to separate
the MSTAR data by class even though it has never been trained on such data!

###  III-B Extended Operating Conditions

For the MSTAR EOC test we train the classification models on 17∘ elevation
imagery and test on 30∘ data [22]. The large difference in collection geometry
causes a sizable change in the target signatures and offers a more challenging
test of generalization. Our results are reported in Table II. We see several
similar trends to the SOC results w.r.t. the #-shots and #-ways. In all cases,
the SimCLR-based models outperform the train-from-scratch baselines, while the
basic and OE classifiers perform very similarly to each other. Interestingly,
the largest margins over the baselines are at the lowest number of shots. For
example, at 5-shots our average improvement over the best baseline is +9%.

TABLE II: MSTAR Extended Operating Condition Accuracy Results

| 2-way | 4-way
---|---|---
Method |  5-shot |  10-shot |  20-shot |  25-shot |  5-shot |  10-shot |  20-shot |  25-shot
Scratch-AConv | 68.9 | 78.4 | 87.0 | 89.3 | 53.8 | 66.2 | 74.3 | 76.3
Scratch-RN18 | 72.1 | 82.1 | 86.6 | 86.9 | 56.0 | 65.8 | 73.2 | 74.5
1-9 |  |  |  |  |  |  |  |
SimCLR+basic | 78.8 | 85.4 | 90.3 | 90.7 | 65.2 | 74.2 | 80.5 | 82.0
SimCLR+OE | 80.9 | 85.7 | 89.6 | 89.8 | 65.3 | 74.1 | 79.7 | 81.0

###  III-C Out-of-Distribution Detection

Our goal in the third experiment is to examine OOD detection performance
across a spectrum of difficulties/granularities. To do this we leverage OOD
data from various sources. For a particularly hard OOD test we use a holdout
scheme. Starting with a 7-way classifier, where in each test iteration 7 of
the 10 MSTAR classes are randomly selected for use in the ID label space, we
use data from the remaining 3 classes as Holdout OOD samples. For medium
difficulty OOD, we use data from SARSIM-Roads, SARSIM-Medium, and SARSIM-Grass
[6]. This is synthetic SAR data which has been generated to specifically match
the collection conditions of MSTAR with different backgrounds (note, we remove
the Bulldozer and Tank classes to avoid potential ID/OOD ambiguities). For
coarse-grained OOD we generate FakeData by creating random Uniform noise chips
and use the MNIST test set which contains hand-written digits. Both are
clearly OOD w.r.t. any SAR classifier. Fig. 3 shows examples from each of
these datasets. Importantly, keep in mind that none of the test OOD data is in
$\mathcal{D}_{pretrain}$ (the OE set), meaning the model has never seen it
before.

Table III shows the average OOD detection performance (as measured by % AUROC
[20]) of a 7-way MSTAR SOC classifier at different #-shots. Note, the AUROC
value can be interpreted as the probability that an ID test input would have a
greater $\mathcal{S}_{ID}$ score than an OOD input [23], where 100% is a
perfect detector. At each setting we examine the performance of the basic
versus OE classifier, which we know both achieve high ID accuracy from Sec.
III-A. The most important takeaway from Table III is that the OE classifier is
better than the basic classifier in all cases by a large margin.
Interestingly, for the basic model the coarse-grained OOD sets are the hardest
to detect, while the OE classifier detects these almost perfectly. This
surprising behavior is something we plan to study in a future work. Another
interesting finding is that detection of Holdout samples improves
significantly as #-shots increases, while detection of the other OOD types
stays relatively constant w.r.t. #-shots. Finally, a somewhat expected result
is that for OE models, the Holdout samples are the most difficult to detect,
while the medium- and coarse-grained OOD inputs are less challenging. Given
that these results still show room for improvement, we highly encourage future
system designers (whether in few-shot scenarios or not) to use OOD-aware
classifier training schemes such as OE to build robustness to the many forms
of OOD inputs.

Figure 3: Sample test chips from each dataset used in the OOD experiment.
TABLE III: MSTAR Out-of-Distribution Detection Results for 7-way classifier (%
AUROC)

| Notation = basic / OE
---|---
OOD Set | 1-shot | 5-shot | 10-shot | 20-shot | 25-shot
Holdout | 54.1 / 56.0 | 64.1 / 67.2 | 70.2 / 76.2 | 78.0 / 81.9 | 81.7 / 84.0
SARSIM-R | 59.6 / 93.2 | 56.3 / 95.9 | 58.8 / 96.6 | 66.0 / 97.2 | 68.3 / 97.3
SARSIM-M | 55.5 / 88.4 | 50.8 / 90.5 | 53.9 / 92.4 | 60.3 / 93.5 | 61.2 / 93.6
SARSIM-G | 44.8 / 87.9 | 40.3 / 87.0 | 44.2 / 88.1 | 50.6 / 87.9 | 52.5 / 87.2
FakeData | 42.1 / 98.5 | 27.0 / 99.5 | 21.7 / 99.6 | 21.6 / 99.5 | 21.3 / 99.5
MNIST | 15.1 / 94.1 | 15.5 / 95.8 | 16.6 / 96.3 | 21.1 / 95.5 | 22.2 / 95.8

###  III-D Generalization vs. Detection Trade-off

Our final experiment is an analysis that deliberates the question: can we have
it all? Specifically, can we have high accuracy on ID data in SOC and high
accuracy on ID classes in EOC and reliable OOD detection? To answer this
question we create Fig. 4, which shows density plots of $\mathcal{S}_{ID}(x)$
scores measured on a SimCLR+OE model for ID test data in SOC (SOC-ID), ID test
data in EOC (EOC-ID), and Holdout OOD data in SOC (Holdout-OOD).

Firstly, we observe that the SOC-ID and Holdout-OOD are relatively separable,
which matches the findings in Table III. However, the EOC-ID and Holdout-OOD
are not very separable, as both are predicted with relatively low scores.
Thus, with our method if you require good Holdout-OOD detection then most of
the EOC-ID would also get flagged as OOD. Conversely, if you require high
recall on EOC-ID then many Holdout-OOD would pass through the detector. These
results indicate that we cannot have it all with this design, and highlight a
serious trade-off to be considered in future work. Finally, while this trade-
off may not be completely unique to the FSL setting, our previous observation
that less shots makes Holdout-OOD detection harder (Sec. III-C) means this
trade-off may be more dire in the few-shot setting.

Figure 4: Density plots of $\mathcal{S}_{ID}$ scores for ID data in SOC, ID
data in EOC, and Holdout OOD data. Also shown are hypothetical thresholds to
achieve 80% True Positive Rates on ID data in SOC and EOC.

##  IV Conclusion

We now confirm our initial hypothesis (Sec. I) that modern SSL techniques like
SimCLR are capable of learning highly generalizable SAR feature extractors
from large pools of diverse and unlabeled SAR data. As evidence, our
experiments show that a model trained mostly on unlabeled SAR ships is able to
provide an informative-enough feature space to do few-shot classification of
MSTAR targets. This result provides merit to the under-studied global model
approach to SAR ATR and confirms that SAR features can be highly transferable,
even across sensors, imaging modes, target types, etc. At a more concrete
level, we provide a method for performing robust few-shot SAR ATR in a limited
environment where there is no labeled data for pretraining. We also provide a
methodology that boosts OOD detection performance of the classifier without
adding significant assumptions or overhead. We hope that this work motivates
further study on few-shot SAR ATR and representation learning for SAR, with
specific considerations for robustness to OOD inputs and the trade-off between
OOD detection and generalization.

Lastly, we provide some suggestions for future work. We believe there are
potential gains in developing improved domain-relevant augmentations for both
Stage 1 and Stage 2 training. We also believe that adding synthetically
generated data crafted to be relevant to the downstream task may improve the
quality of learned features. Next, one may investigate the use of different
SSL algorithms and backbone model architectures. Finally, an interesting study
would be to evaluate when it becomes beneficial to fine-tune the feature
extractor instead of keeping it fixed during Stage 2 training.

Public Release Number: AFRL-2022-3418

## References

  * [1] U. K. Majumder, E. P. Blasch, and D. A. Garren, _Deep Learning for Radar and Communications Automatic Target Recognition_. Artech House, 2020.
  * [2] J. Snell, K. Swersky, and R. Zemel, “Prototypical networks for few-shot learning,” in _Advances in Neural Information Processing Systems_ , I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, Eds., vol. 30, 2017.
  * [3] L. Zhang, X. Leng, S. Feng, X. Ma, K. Ji, G. Kuang, and L. Liu, “Domain knowledge powered two-stream deep network for few-shot sar vehicle recognition,” _IEEE Transactions on Geoscience and Remote Sensing_ , vol. 60, pp. 1–15, 2022.
  * [4] T. Chen, S. Kornblith, M. Norouzi, and G. E. Hinton, “A simple framework for contrastive learning of visual representations,” in _International Conference on Machine Learning_ , 2020, pp. 1597–1607.
  * [5] Y. Tian, Y. Wang, D. Krishnan, J. B. Tenenbaum, and P. Isola, “Rethinking few-shot image classification: A good embedding is all you need?” in _European Conference on Computer Vision_ , vol. 12359, 2020, pp. 266–282.
  * [6] D. Malmgren-Hansen, A. Kusk, J. Dall, A. A. Nielsen, R. Engholm, and H. Skriver, “Improving sar automatic target recognition models with transfer learning from simulated data,” _IEEE Geoscience and Remote Sensing Letters_ , vol. 14, no. 9, pp. 1484–1488, 2017.
  * [7] D. Hendrycks, M. Mazeika, and T. G. Dietterich, “Deep anomaly detection with outlier exposure,” in _International Conference on Learning Representations_ , 2019.
  * [8] N. Inkawhich, E. Davis, M. Inkawhich, U. K. Majumder, and Y. Chen, “Training sar-atr models for reliable operation in open-world environments,” _IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing_ , vol. 14, pp. 3954–3966, 2021.
  * [9] N. Inkawhich, J. Zhang, E. K. Davis, R. Luley, and Y. Chen, “Improving out-of-distribution detection by learning from the deployment environment,” _IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing_ , vol. 15, pp. 2070–2086, 2022.
  * [10] T. Ross, S. Worrell, V. Velten, J. Mossing, and M. Bryant, “Standard sar atr evaluation experiments using the mstar public release data set,” in _SPIE Conference on Algorithms for Synthetic Aperture Radar Imagery V_ , 1998\.
  * [11] Y. Wang, C. Wang, H. Zhang, Y. Dong, and S. Wei, “A sar dataset of ship detection for deep learning under complex backgrounds,” _Remote Sensing_ , vol. 11, no. 7, p. 765, Mar 2019.
  * [12] S. Wei, X. Zeng, Q. Qu, M. Wang, H. Su, and J. Shi, “Hrsid: A high-resolution sar images dataset for ship detection and instance segmentation,” _IEEE Access_ , vol. 8, pp. 120 234–120 254, 2020.
  * [13] X. Hou, W. Ao, Q. Song, J. Lai, H. Wang, and F. Xu, “Fusar-ship: building a high-resolution SAR-AIS matchup dataset of gaofen-3 for ship detection and recognition,” _Sci. China Inf. Sci._ , vol. 63, no. 4, 2020.
  * [14] T. Zhang, X. Zhang, J. Li, X. Xu, B. Wang, X. Zhan, Y. Xu, X. Ke, T. Zeng, H. Su, I. Ahmad, D. Pan, C. Liu, Y. Zhou, J. Shi, and S. Wei, “Sar ship detection dataset (ssdd): Official release and comprehensive data analysis,” _Remote Sensing_ , vol. 13, no. 18, 2021.
  * [15] T. Zhang, X. Zhang, X. Ke, X. Zhan, J. Shi, S. Wei, D. Pan, J. Li, H. Su, Y. Zhou, and D. Kumar, “Ls-ssdd-v1.0: A deep learning dataset dedicated to small ship detection from large-scale sentinel-1 sar images,” _Remote Sensing_ , vol. 12, no. 18, 2020.
  * [16] Y. Hu, Y. Li, and Z. Pan, “A dual-polarimetric sar ship detection dataset and a memory-augmented autoencoder-based detection method,” _Sensors_ , vol. 21, no. 24, 2021.
  * [17] S. Lei, D. Lu, X. Qiu, and C. Ding, “Srsdd-v1.0: A high-resolution sar rotation ship detection dataset,” _Remote Sensing_ , vol. 13, no. 24, 2021\.
  * [18] K. E. Dungan, C. Austin, J. Nehrbass, and L. C. Potter, “Civilian vehicle radar data domes,” in _SPIE Algorithms for Synthetic Aperture Radar Imagery XVII_ , 2010.
  * [19] T. Silva, “Simclr,” https://github.com/sthalles/SimCLR, 2022.
  * [20] S. Liang, Y. Li, and R. Srikant, “Enhancing the reliability of out-of-distribution image detection in neural networks,” in _International Conference on Learning Representations_ , 2018.
  * [21] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” in _IEEE Conference on Computer Vision and Pattern Recognition_ , 2016.
  * [22] S. Chen, H. Wang, F. Xu, and Y. Jin, “Target classification using the deep convolutional networks for sar images,” _IEEE Transactions on Geoscience and Remote Sensing_ , vol. 54, no. 8, pp. 4806–4817, 2016.
  * [23] D. Hendrycks and K. Gimpel, “A baseline for detecting misclassified and out-of-distribution examples in neural networks,” in _International Conference on Learning Representations_ , 2017.

## Supplemental Materials

### A. SimCLR Augmentations for SAR Data

In this section we will discuss the augmentation pipeline used to train the
SimCLR model in more detail. Fig. 6 shows the PyTorch transform block that
defines the augmentations. Notice, several of the operations are from
torchvision.transforms and we leave description of those to the official docs
(https://pytorch.org/vision/0.8/transforms.html). The remaining transforms are
defined as follows:

  * •

ClipAndScale(a,b) \- linearly re-scale the range of pixel values in the image
given a randomly chosen maximum value between a and b. Pseudo:
$\text{max\\_val}=random.uniform(a,b);x_{new}=clip(x_{orig},0,\text{max\\_val})/\text{max\\_val}$.

  * •

PowScale(a,b) \- pixel-wise exponentiation using a randomly selected value
between a and b. Pseudo:
$\text{val}=random.choice([random.uniform(a,1),random.uniform(1,b)]);x_{new}=(x_{orig})^{\text{val}}$.

  * •

SpeckleNoise(a,b) \- randomly replace a subset of pixel values with values
chosen from a Uniform(0,1) distribution. The size of the subset in % total
pixels is decided with $random.uniform(a,b)$.

  * •

GaussianNoise(a,b) \- apply pixel-wise additive Gaussian noise using
$stdev=random.uniform(a,b)$.

  * •

GaussianBlur \- convolve the chip with a Gaussian filter.

Finally, Fig. 6 shows several examples of augmented pairs that are generated
with this transform block. Recall, the SimCLR training objective encourages
the two different views in the pair to lie near each-other in projection
space, while views from different pairs should lie far away.

Figure 5: PyTorch-like code describing transforms used to train SimCLR model.

Figure 6: Examples of augmented pairs.

### B. Comparison With A Powerful Recent Method for Few-Shot MSTAR

In the main document we do not compare our method to meta-learning-based FSL
algorithms. As mentioned several times, this is because meta-learning
approaches assume the existence of a labeled dataset for pretraining the
feature extractor, which is a violation of our basic assumption set. However,
for the sake of completeness we now present a detailed comparison with a
state-of-the-art meta-learning approach to few-shot MSTAR called Domain
Knowledge Powered Two Stream Deep Network (DKTS-N) [3]. Below are some notable
points in no particular order:

  * •

DKTS-N uses the labeled SARSIM dataset [6] as the source of pretraining data
for meta-learning. We argue that this is essentially an ideal scenario because
SARSIM is specifically/conveniently synthesized to match the conditions of
MSTAR (i.e., the downstream FSL task). We believe that this implies a more
significant assumption than any we make in our work. Also, it is unclear if
the method would be nearly as effective if the pretraining dataset was not so
closely related to the downstream task. Finally, although our method is not
reliant on such a related dataset, in practice we could leverage more data by
adding it to $\mathcal{D}_{pretrain}$.

  * •

DKTS-N reports results using a preprocessing function that performs perfect
azimuth angle normalization on both the training and test data. We believe
this is done using meta-data provided in the MSTAR dataset. They also show
that as the normalization function incurs errors, performance can degrade
significantly. It is unclear how one would achieve perfect azimuth angle
normalization in the “wild,” and thus we do not use it here. However, if such
a function were developed, we believe it could complement our model in a
similar way to theirs.

  * •

DKTS-N does not consider OOD detection in the few-shot classifier and thus is
susceptible to producing erroneous predictions on OOD data at test time.

  * •

DKTS-N involves a time-intensive inference procedure which scales poorly as
#-shots increases (refer Table VII of [3]). This is because for each test
input it performs an iterative optimization procedure w.r.t each support
sample. Our inference procedure is not nearly as complex, and is a single
forward pass through a DNN (e.g., a RN18).

  * •

While DKTS-N shows very high performance in MSTAR SOC, our method actually
outperforms it in several of the EOC settings. Specifically, in the (4-way,
10-shot) and (4-way, 25-shot) cases our method is over 3% more accurate.

### C. Additional results using a different self-supervised pretraining
algorithm

In the main manuscript we use the SimCLR algorithm to pretrain the feature
extractor in Stage 1 of the pipeline. SimCLR is often used as a standard
baseline for comparison in Self-Supervised Learning (SSL) literature and our
intention is to use it as a representative of the SSL-class of algorithms. To
be clear, our framework is not tied to the SimCLR algorithm and we believe
that any modern SSL algorithm may be used for pretraining. To show the
potential impact of a different SSL algorithm we run experiments with the
Bootstrap Your Own Latent† (BYOL) method. BYOL uses a non-contrastive
distillation scheme to perform the representation learning. This is thought to
offer some advantages over SimCLR because the learning signal is not reliant
on many negative samples to contrast the positive pair with. Because the field
of SSL is advancing so quickly it is unclear which SSL paradigm will
ultimately emerge as the “best,” but for now these two methodologies make for
an interesting comparison. As mentioned in the paper, an important future work
is still to consider using different SSL algorithms for pretraining within our
framework.

Table IV shows the extended few-shot classification results in MSTAR Standard
Operation Conditions, where the Scratch and SimCLR numbers are copy-pasted
from Table I in the main paper. The BYOL training parameters are nearly
identical to the SimCLR parameters, including the extractor architecture,
optimization procedure, augmentation scheme, and training schedule. The only
method specific parameter is the EMA momentum which we set to $\tau=0.9995$.
In the table we see that BYOL actually outperforms the Scratch and SimCLR
models in all scenarios. Interestingly, the benefits of BYOL over SimCLR are
most clear in the very low-shot scenarios (i.e., 1- and 5-shots) where the
average margin of improvement is $\sim$5%. Focusing on the 10-way scenario,
BYOL outperforms the best scratch model by at least 12% in the 1-10 shot
range. Finally, we note that BYOL-based models trained with 10-shots perform
on par with Scratch models trained with 20-shots! The ease at which BYOL is
swapped into our framework is a key advantage and speaks to the flexibility
and extendability we were striving for.

TABLE IV: Extended MSTAR Standard Operating Condition Accuracy Results (with
95% confidence interval)

| 2-way | 5-way | 10-way
---|---|---|---
Method | 1-shot | 5-shot | 10-shot | 20-shot | 25-shot | 1-shot | 5-shot | 10-shot | 20-shot | 25-shot | 1-shot | 5-shot | 10-shot | 20-shot | 25-shot
Scratch-AConv | 59.6$\pm$1.2 | 73.6$\pm$1.2 | 84.3$\pm$0.9 | 92.1$\pm$0.7 | 93.8$\pm$0.6 | 31.6$\pm$0.6 | 55.9$\pm$0.8 | 72.3$\pm$0.8 | 85.7$\pm$0.5 | 88.5$\pm$0.4 | 20.9$\pm$0.3 | 47.4$\pm$0.4 | 65.0$\pm$0.4 | 80.7$\pm$0.2 | 84.5$\pm$0.2
Scratch-RN18 | 55.8$\pm$1.2 | 74.9$\pm$1.4 | 84.1$\pm$1.1 | 89.6$\pm$0.8 | 92.1$\pm$0.6 | 30.6$\pm$0.7 | 53.7$\pm$0.8 | 67.0$\pm$0.7 | 79.3$\pm$0.6 | 83.2$\pm$0.5 | 20.7$\pm$0.3 | 43.4$\pm$0.4 | 58.9$\pm$0.4 | 74.0$\pm$0.4 | 78.5$\pm$0.3
1-16 |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
SimCLR+basic | 64.5$\pm$1.4 | 83.0$\pm$1.1 | 91.4$\pm$0.7 | 96.1$\pm$0.4 | 97.4$\pm$0.2 | 38.8$\pm$0.7 | 68.4$\pm$0.7 | 82.6$\pm$0.5 | 91.8$\pm$0.3 | 94.0$\pm$0.2 | 28.3$\pm$0.3 | 59.2$\pm$0.3 | 75.2$\pm$0.2 | 87.7$\pm$0.1 | 90.6$\pm$0.1
SimCLR+OE | 65.4$\pm$1.5 | 86.4$\pm$1.0 | 92.2$\pm$0.8 | 96.5$\pm$0.4 | 97.5$\pm$0.2 | 39.9$\pm$0.7 | 69.7$\pm$0.6 | 83.3$\pm$0.5 | 92.3$\pm$0.3 | 93.8$\pm$0.2 | 28.8$\pm$0.3 | 60.5$\pm$0.3 | 75.6$\pm$0.2 | 87.0$\pm$0.1 | 89.3$\pm$0.1
1-16 |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
BYOL+basic | 72.4$\pm$1.9 | 86.7$\pm$1.1 | 93.8$\pm$0.6 | 97.1$\pm$0.3 | 98.0$\pm$0.2 | 44.7$\pm$0.9 | 74.5$\pm$0.6 | 86.0$\pm$0.5 | 93.9$\pm$0.2 | 94.8$\pm$0.2 | 32.8$\pm$0.4 | 65.1$\pm$0.3 | 79.5$\pm$0.2 | 89.8$\pm$0.1 | 92.1$\pm$0.1
BYOL+OE | 70.1$\pm$1.8 | 87.6$\pm$1.1 | 94.7$\pm$0.5 | 97.8$\pm$0.2 | 98.1$\pm$0.2 | 45.9$\pm$0.8 | 75.9$\pm$0.7 | 87.4$\pm$0.4 | 93.8$\pm$0.2 | 94.9$\pm$0.2 | 33.7$\pm$0.4 | 66.2$\pm$0.4 | 79.8$\pm$0.2 | 88.7$\pm$0.1 | 90.7$\pm$0.1

[$\dagger$] Jean-Bastien Grill, Florian Strub, Florent Altché, Corentin
Tallec, Pierre H. Richemond, Elena Buchatskaya, Carl Doersch, Bernardo Avila
Pires, Zhaohan Daniel Guo, Mohammad Gheshlaghi Azar, Bilal Piot, Koray
Kavukcuoglu, Rémi Munos, and Michal Valko. “Bootstrap your own latent: A new
approach to self-supervised learning.” In Advances in Neural Information
Processing Systems, 2020.


A QUANTILE SHIFT APPROACH TO MAIN EFFECTS AND INTERACTIONS IN A 2-BY-2 DESIGN

Rand R. Wilcox1∗ & Guillaume A. Rousselet2∗
1\. Dept. of Psychology, University of Southern California, Los Angeles, CA
90089-1061, USA orcid.org/0000-0002-2524-2976
2\. School of Psychology and Neuroscience, College of Medical, Veterinary and
Life
Sciences, University of Glasgow, 62 Hillhead Street, G12 8QB, Glasgow, UK
∗ Corresponding authors<EMAIL_ADDRESS><EMAIL_ADDRESS>

###### Abstract

When comparing two independent groups, shift functions are basically
techniques that compare multiple quantiles rather than a single measure of
location, the goal being to get a more detailed understanding of how the
distributions differ. Various versions have been proposed and studied. This
paper deals with extensions of these methods to main effects and interactions
in a between-by-between, 2-by-2 design. Two approaches are studied, one that
compares the deciles of the distributions, and one that has a certain
connection to the Wilcoxon–Mann–Whitney method. For both methods, we propose
an implementation using the Harrell–Davis quantile estimator, used in
conjunction with a percentile bootstrap approach. We report results of
simulations of false and true positive rates.

Keywords: shift function, deciles, factorial design, interaction, quantile
estimation

##  1 Introduction

When comparing two distributions, certainly the most common approach is to
focus on a single measure of location, typically the mean or median. An
alternative approach is to compare multiple quantiles with the goal of getting
a more detailed understanding of where groups differ and by how much
(Rousselet et al., 2017). Several methods have been derived for dealing with
this issue (e.g., Doksum & Sievers, 1976; Goldman & Kaplan, 2018; Lombard,
2005; Wilcox, 1995; Wilcox et al., 2014). The goal in this paper is to suggest
extensions of these methods when dealing with a between-by-between, 2-by-2
design.

Here, two distinct approaches are studied. The first defines interactions and
main effects with means replaced by a collection of quantiles with an emphasis
on the deciles. For example, if $\theta_{jk}$ denotes the .2 quantile
corresponding to level $j$ of the first factor and level $k$ of second factor,
one goal is to test

$H_{0}:\theta_{11}-\theta_{12}=\theta_{21}-\theta_{22},$ (1)

which mimics the usual notion of an interaction in an obvious way. Here,
however, the goal is to use multiple quantiles and to assess how well the
familywise error rate (the probability of one or more Type I errors) is
controlled. Of course, a related issue is computing a reasonably accurate
confidence interval for $\theta_{11}-\theta_{12}-\theta_{21}+\theta_{22}$. As
is evident, main effects can be addressed as well. For example, one can test

$H_{0}:\theta_{11}+\theta_{12}=\theta_{21}+\theta_{22},$ (2)

for a collection of quantiles.

The second approach, when dealing with an interaction, has a certain
connection to a rank-based method proposed by Patel and Hoel (1973) that in
turn has a connection with the classic Wilcoxon–Mann–Whitney test. To describe
the Patel–Hoel approach, let $X_{jk}$ denote four independent random variables
where $X_{11}$ and $X_{12}$ correspond to the first level of the first factor
in a 2-by-2 design, while $X_{21}$ and $X_{22}$ correspond to the second level
of the first factor. Let $p_{1}=P(X_{11}<X_{12})$ and
$p_{2}=P(X_{21}<X_{22})$. The hypothesis of no interaction is

$H_{0}:p_{1}=p_{2}.$ (3)

And there is the issue of computing a $1-\alpha$ confidence interval for
$p_{1}-p_{2}$. Wilcox (2022, section 7.9.2) describes a method for making
inferences about this measure of effect size that performs well in
simulations. For some related rank-based methods, see Gao and Alvo (2005), as
well as De Neve and Thas (2017).

Let $X_{ijk}$ denote a random sample from the $j$th level of the first factor
and the $k$th level of the second factor ($i=1,\ldots,n_{jk}$; $j=1$, 2;
$k=1$, 2). For convenience momentarily focus on $p_{1}$ and let
$D_{ih}=X_{i11}-X_{h12}$ ($i=1,\ldots,n_{11}$; $h=1,\ldots,n_{12}$). An
estimate of $p_{1}$ is simply

$\hat{p}_{1}=\frac{1}{n_{1}n_{2}}\sum\sum I(D_{ih}),$ (4)

where the indicator function $I(D_{ih})=1$ if $D_{ih}<0$; otherwise
$I(D_{ih})=0$. The estimator $\hat{p}_{1}$ is the estimator used by the
classic Wilcoxon–Mann–Whitney test.

Let $\theta_{1}$ denote the median of $D$ and notice that

$H_{0}:p_{1}=.5$ (5)

is the same as

$H_{0}:\theta_{1}=0.$ (6)

The parameter $\theta_{1}$ is defined based on level 1 of the first factor.
Let $\theta_{2}$ denote the analog of $\theta_{1}$ when dealing with level 2
of the first factor. Then an analog of the Patel–Hoel interaction is
$\theta_{1}-\theta_{2}$. Here, however, the goal is to consider the broader
issue of comparing the deciles of these two distributions. More formally, let
$q_{1}$ and $q_{2}$ denote the $q$th quantile of the distribution $D$ for
level 1 of the first factor and level 2 of first factor, respectively. The
goal is to test

$H_{0}:q_{1}=q_{2}$ (7)

and to compute a $1-\alpha$ confidence interval for $q_{1}-q_{2}$ for
$q=.1,.2,\ldots,.9$. A second goal is to control the FWE rate in a reasonably
accurate manner.

Notice that when testing (1), there is no distinction between
$\theta_{11}-\theta_{12}=\theta_{21}-\theta_{22}$ and
$\theta_{11}-\theta_{22}=\theta_{12}-\theta_{22}$. That is, interchanging the
rows and columns does not alter the estimated effect size. However, when
dealing with (7), interchanging the rows and columns can yield different
results.

To illustrate this point, consider, for example, a situation where for the
first level of the first factor, both levels of the second factor have
standard normal normal distributions, while for the second level of the first
factor, the two levels of the second factor have lognormal distributions. Data
were generated as just described based on sample sizes of 50 and the 50th
quantile was estimated. The estimate of $q_{1}-q_{2}$ was .027. But
interchanging the rows and columns, now the estimate was $-.088$. The code to
reproduce this example can be found in the R notebook `apd_ex.Rmd`, which is
available as part of the companion reproducibility package for this article
(Wilcox & Rousselet, 2023) and on GitHub (Rousselet, 2023).

##  2 The Proposed Methods

First consider testing (1). The first issue is choosing a reasonable quantile
estimator from among the many estimators that might be used. The focus here is
on the estimator derived by Harrell and Davis (1982) with the understanding
that no single estimator dominates in terms of efficiency. There are several
quantile estimators that are based on a weighted average of just two order
statistics (Hyndman & Fan, 1996). In contrast, the Harrell–Davis estimator
uses a weighted average of all the order statistics. A possible appeal of the
Harrell–Davis estimator is that it has better efficiency under normality
compared to using a weighted average of only two order statistics. Moreover,
it is known that extant methods for making inferences about quantiles, when
using only two order statistics, can be unsatisfactory when there are tied
values (Wilcox, 2022). An exception, when comparing two independent groups, is
the usual sample median when coupled with a percentile bootstrap method.

There are several quantile estimators, in addition to the Harrell–Davis
estimator, that use all of the order statistics (e.g., Liu et al., 2022;
Navruz, G. & Özdemir, 2020). A possible criticism is that these estimators
have a breakdown point of only $1/n$. That is, the minimum number of order
statistics that must be altered to be make the estimate arbitrarily large is
one. In practice this issue might not be a serious concern because the extreme
order statistics get a relatively small weight. In a situation where the
breakdown point is an issue, one possibility is to use the trimmed
Harrell–Davis estimator derived by Akinshin (2022).

Let $U$ be a random variable having a beta distribution with parameters
$a=(n+1)q$ and $b=(n+1)(1-q)$. Let

$W_{i}=P\left(\frac{i-1}{n}\leq U\leq\frac{i}{n}\right).$

The Harrell–Davis estimate of the $q$th quantile is

$\hat{\theta}_{q}=\sum_{i=1}^{n}W_{i}X_{(i)},$ (8)

where $X_{(1)}\leq\cdots\leq X_{(n)}$ are the values written in ascending
order. The beta weights used to calculate the deciles for a sample size of
$n=50$ are illustrated in Figure 1.

Figure 1:  Beta weights used to calculate the Harrell-Davis estimates of the
deciles with a sample size of 50.

As an alternative to the Harrell–Davis estimator, we considered the default
quantile estimator in R, called using the command `quantile(x,
probs=seq(0.1,0.9,0.1), type=7)` to compute the deciles. This estimator, in
addition to being widely used, relies on two order statistics (Hyndman & Fan,
1996), and could thus be more robust to outliers than the Harrell–Davis
estimator in some situations.

To test (1), we considered a percentile bootstrap method combined with the
Harrell–Davis estimator and the quantile(type=7) estimator. The percentile
bootstrap method begins by generating bootstrap samples by sampling with
replacement $n_{jk}$ values from the data associated with the $j$th level of
the first factor and the $k$th level of second factor yielding $X^{*}_{ijk}$
($i=1,\ldots,n_{jk}$). Based on these bootstrap samples, compute the $q$th
quantile using the Harrell–Davis estimator, or the quantile(type=7) estimator,
yielding $\theta^{*}_{jk}$ followed by

$\Psi^{*}=\theta^{*}_{11}-\theta^{*}_{12}-\theta^{*}_{21}+\theta^{*}_{22}.$
(9)

Repeat this process $B$ times yielding $\Psi^{*}_{1},\ldots,\Psi^{*}_{B}$. Let
$A$ denote the number of $\Psi^{*}$ values that are less than zero and let $D$
denote the number of $\Psi^{*}$ values that are equal to zero. Let

$P=\frac{A}{B}+.5\frac{D}{B}.$ (10)

A p-value, when testing (1), is $2\min\\{P,1-P\\}$. The term $D/B$ is
important when dealing with tied values (e.g., Wilcox, 2022). To compute a
$1-\alpha$ confidence interval, put $\Psi^{*}_{1},\ldots,\Psi^{*}_{B}$ in
ascending order yielding $\Psi^{*}_{(1)}\leq\cdots\leq\Psi^{*}_{(B)}$. Let
$\ell=\alpha B/2$, rounded to the nearest integer. Let $u=n-\ell$. Then a
$1-\alpha$ confidence interval for $\Psi$ is

$(\Psi^{*}_{(\ell+1)},\Psi^{*}_{(u)}).$ (11)

Here, $B=2000$ is used.

Note that the same bootstrap samples were used for each of the quantiles being
compared. An alternative approach is to use separate bootstrap samples for
each test to be performed (Wilcox, 1995). Both approaches were considered and
there is no indication that separate bootstrap samples offer a practical
advantage in terms of controlling the Type I error probability. Results of the
simulations comparing the two bootstrap approaches are available in the R
notebooks `sim_fp_b1b9.Rmd` and `sim_fp_apd_b1b9.Rmd` (Wilcox & Rousselet,
2023). Using the same bootstrap samples for all of the tests performed
considerably reduces execution time, which is why it is assumed henceforth.

When performing $C$ tests, there is the issue of controlling the FWE
(familywise error) rate, meaning the probability of one or more Type I errors.
Two approaches are considered here. The first is Hochberg’s (1988) improvement
on the Bonferroni method. Let $p_{1},\dots,p_{C}$ be the p-values associated
with the $C$ tests. Put these p-values in descending order, and label the
results $p_{[1]}\geq p_{[2]}\geq\cdots\geq p_{[C]}$. Set $k=0$ and proceed as
follows:

  1. 1.

Increment $k$ by 1. If

$p_{[k]}\leq\frac{\alpha}{k},$

stop and reject all hypotheses having a p-value less than or equal to
$p_{[k]}$.

  2. 2.

If $p_{[k]}>\alpha/k$, repeat step 1.

  3. 3.

Repeat steps 1 and 2 until a significant result is obtained or all $C$
hypotheses have been tested.

The second method for controlling the FWE rate was Benjamini and Hochberg
(1995), which is aimed at controlling the false discovery rate (FDR). That is,
the goal is to control the expected proportion of Type I errors among the null
hypotheses that are correct. It is known that there are situations where the
Benjamini–Hochberg method ensures that the false discovery rate is less than
or equal to $\alpha$, but it does not necessarily control the FWE rate
(Hommel, 1988). However, the simulations described in section 3, below,
indicated that when using Hochberg’s method, the FWE rate is always below the
nominal level. Consequently, there was interest in how well the
Benjamini–Hochberg method performs. It is readily verified that the
Benjamini–Hochberg method always has as much or more power than Hochberg’s
method. Consequently, provided the Benjamini–Hochberg method controls the FWE
rate for the situation at hand, it has a practical advantage over Hochberg’s
method.

The Benjamini–Hochberg method is applied simply by replacing
$p_{[k]}\leq\alpha/k$ in step 1 of Hochberg’s method with

$p_{[k]}\leq\frac{(C-k+1)\alpha}{C}$ (12)

As for testing (7), again both the Harrell–Davis and the quantile(type=7)
estimators were considered, in conjunction with a percentile bootstrap method.
Bootstrap samples are generated as before yielding estimates of $q_{1}$ and
$q_{2}$, which are labelled $q^{*}_{1}$ and $q^{*}_{2}$. This process is
repeated $B$ times and the results are used to compute a $1-\alpha$ confidence
interval for $q_{1}-q_{2}$, as well as a p-value, in essentially the same
manner as done when testing (1).

##  3 Simulation Results

Simulations were used to check the FWE rate when testing (1) and (7), as well
as how the power of these methods compare to the classic ANOVA F test as well
as a method for comparing 20% trimmed mean, which is described in Wilcox,
2022, section 7.4.3). Data were generated from four continuous distributions
as well as three discrete distributions. The discrete distributions were a
Poisson distribution having mean 9, and two beta-binomial distributions, one
with parameter $r=1$, the other with $r=9$, and the other parameters set to
$s=9$, and 10 bins. The three discrete distributions were included as a
partial check on the impact of tied values. The four continuous distributions
were standard normal, mixed normal, lognormal and mixed lognormal. The
distribution of the mixed normal is

$H(x)=.9\Phi(x)+.1\Phi(x),$ (13)

where $\Phi(x)$ is the standard normal distribution. The mixed normal is a
symmetric distribution with heavy tails, roughly meaning that outliers are
likely to occur. The mixed lognormal distribution is given by (13), but with
$\Phi$ replaced by the lognormal distribution.

Based on over 1,500 estimates of skewness and kurtosis reported in various
journal articles, Cain et al. (2017) report that 99% of the estimates were
less than the skewness and kurtosis of a lognormal distribution. This suggests
that if a method performs reasonably when dealing with a lognormal
distribution, it is highly likely to perform reasonably well in practice.
However, a possible concern is that when dealing with heavy-tailed
distributions, the standard error of the usual kurtosis estimator can be quite
high even when the sample size is fairly large. Moreover, the usual estimate
of kurtosis can grossly under estimate the true value.

Consider, for example, a lognormal distribution, which has kurtosis 113.9.
Based on a sample of 100,000, the kurtosis of the lognormal distribution was
estimated and this process was repeated 1000 times. It was found that 79% of
the estimates were less than the true value. The median estimate was 82. This
process was repeated using the contaminated lognormal distribution which is
skewed and very heavy-tailed. The kurtosis of this distribution was estimated
based on a sample size of one million. This process was repeated 1000 times
yielding estimates ranging between 242 and 16400. The median estimate was 429.
Consequently, the contaminated lognormal distribution was used here as an
additional check on how well the methods perform. The code for these
simulations is available in the notebook `kurtosis_estimation.Rmd` (Wilcox &
Rousselet, 2023).

Simulations were performed for sample sizes 20, 30, 40, …, 100, using 10,000
iterations. This was done using both types of quantile estimators and for the
two main effects and the interaction. A p-value was computed for each of the
nine quantiles to be compared. In terms of controlling the FWE rate, both
Hochberg’s method and the Benjamini–Hochberg method were considered. And as
previously indicated, this was done for seven distributions. So in terms of
Type I errors there is a total of 3402 results (7 distributions * 2 quantile
estimators * 3 contrasts * 9 deciles * 9 sample sizes).

Finally, simulations were performed dealing with power. Various situations
were considered, including shifting the four continuous distributions used to
estimate Type I errors, varying the skewness of _g-and-h_ distributions
(Hoaglin, 1985), and the skewness of Poisson and beta-binomial distributions.
Complete details, including the code that was used, are reported in files
available on figshare (Wilcox & Rousselet, 2023) and on GitHub (Rousselet,
2023). The main R packages used to perform the simulations and illustrate the
results were `Rcpp` (Eddelbuettel & Francois, 2011), `ggplot2` (Wickham, 2016)
and `cowplot` (Wilke, 2017).

First we consider the results for testing (1), then for testing (7).

###  3.1  Results when testing (1)

The results regarding the FWE rate were highly consistent over the conditions
considered, making it a simple matter to briefly summarize the relative merits
of the methods being considered. All indications are that the actual FWE rate
is less than the nominal .05 level using the Hochberg method as well as the
Benjamini–Hochberg method. Consequently, the Benjamini–Hochberg method is
recommended because it always has as much or more power than Hochberg’s
method.

####  3.1.1 Type I errors

To illustrate some of the type I error simulation results, Figure 2 contains
results for the normal and lognormal populations, separately for the
interaction and the two main effects. Results with Hochberg’s method are
omitted but detailed figures are available in the notebook `sim_fp.Rmd`. The
FWER without correction for multiple comparison (grey lines) is lower than
expected if the 9 tests were independent, which in theory would be
$1-(1-\alpha)^{9}=0.37$. The results in 2 and 3 are about half of that
expected value, supporting the use of the FDR correction over Hochberg’s.

As a sanity check, results for the two types of ANOVAs are included,
confirming false positives very close to the nominal level when sampling from
a normal population. Although the percentile bootstrap led to conservative
estimates for both quantile estimators, the Harrell–Davis estimator clearly
outperforms quantile(type=7).

Figure 2:  False positive results for normal and lognormal populations. FWER
(across quantiles) are plotted as a function of sample size for the
interaction and the two main effects A and B. ANOVA.M = ANOVA using means.
ANOVA.TM = ANOVA using 20% trimmed means. FDR (Benjamini–Hochberg method) and
No corr. (no correction) refer to the percentile bootstrap method in
conjunction with the Harrell–Davis estimator (HD) and the quantile(type=7)
estimator (QT7). Horizontal ribbons: dark grey indicates Bradley’s (1978)
satisfactory range [0.025; 0.075]; light grey indicates Bradley’s (1978) ideal
range [0.045; 0.055].

The gap between the two estimators increases when we consider samples from
discrete populations, because the performance of quantile(type=7) deteriorates
while that of Harrell–Davis remains stable (Figure 3).

Bradley (1978) suggested, as a general guide, that when testing at the .05
level, the actual level should be between .025 and .075. This criterion was
met in most situations for the Harrell–Davis estimator but not for the
quantile(type=7) estimator.

Figure 3:  False positive results for Poisson and one of the two beta-binomial
populations. See details in the Figure 2 caption.

The higher performance of the Harrell-Davis estimator compared to the
quantile(type=7) estimator can be observed at individual deciles as well
(Figure 4). However, for the first and last deciles, the type I error rate is
higher than the nominal level when using the Harrell–Davis estimator and n=20,
a result that confirms earlier observations (Wilcox et al., 2014).

Figure 4:  False positive results at individual deciles for normal and beta-
binomial populations. Type I error rates are plotted at each decile,
separately for the different sample sizes (n) and for the Harrell–Davis
estimator (HD) and the quantile(type=7) estimator (QT7). Horizontal ribbons:
dark grey indicates Bradley’s (1978) satisfactory range [0.025; 0.075]; light
grey indicates Bradley’s (1978) ideal range [0.045; 0.055].

####  3.1.2 Power

Not surprisingly, there are situations where inferences based on means or a
20% trimmed have more power. But there are situations where comparing deciles
provides more power: no single method dominates. To illustrate, Figure 5
presents results for normal and lognormal populations. In each case, data were
generated for the four groups by sampling from the standard normal and
lognormal distributions, before shifting the groups by different amounts. When
sampling from normal populations, the ANOVA on means dominates other methods.
To compare methods, we report familywise power for the decile methods
(bootstrap p values with and without FDR correction), the probability of at
least one rejection among the 9 tests. The ANOVA on 20% trimmed means is a bit
less powerful, followed by the bootstrap method using the Harrell–Davis
estimator, and last the bootstrap method combined with the quantile(type=7)
estimator. Switching to lognormal populations, the power of both ANOVA tests
is dramatically lower than the bootstrap approach. This figure and the next
one were generated using the notebook `sim_tp.Rmd`.

Figure 5:  Power results for normal and lognormal populations. ANOVA.M = ANOVA
using means. ANOVA.TM = ANOVA using 20% trimmed means. FDR (Benjamini–Hochberg
method) and No corr. (no correction) refer to the percentile bootstrap method
in conjunction with the Harrell–Davis estimator (HD) and the quantile(type=7)
estimator (QT7). The black horizontal line marks 0.05.

Figure 6 presents results from populations in which tied values were common.
In both cases, the ANOVA on means dominates other methods. When using the
bootstrap approach, the gap between the Harrell–Davis estimator and the
quantile(type=7) estimator is larger than what was observed when sampling from
continuous populations.

Figure 6:  Power results for Poisson and beta-binomial populations. See
details in Figure 5 caption.

###  3.2 Compare deciles of distributions of all pairwise differences – test
(7)

To assess the bootstrap method aimed at testing (7), we used the same approach
employed in the previous section. Now only the interaction is considered. The
simulations and illustrations of the type I error rates can be found in the
notebook `sim_fp_apd.Rmd`. For power, see notebook `sim_tp_apd.Rmd`.

####  3.2.1 Type I errors

Again we observed FWE rates much lower than expected if the deciles were
independent (Figure 7). For continuous distributions (panels A-D in Figure 7),
results were very similar for the Harrell–Davis and quantile(type=7)
estimators. Both methods were conservative. When sampling from distributions
in which tied values are likely, now the Harrell–Davis estimator outperforms
the quantile(type=7) method (panels E and F in Figure 7).

Figure 7:  FWER results for the comparison of the deciles of distributions of
all pairwise differences.

For continuous distributions, the similarity in performance between the two
quantile methods is evident at the level of individual deciles as well (Figure
8). Under normality, all deciles were associated with type I error rates close
to the nominal level, irrespective of the sample size (panel A in Figure 8).
In the worst situation tested, when sampling from a contaminated lognormal,
results are a bit conservative, especially for the extreme deciles. In both
situation, there is very little separating the two quantile methods.

Figure 8:  Type I error rates for the comparison of the deciles of
distributions of all pairwise differences: continuous distributions.

When sampling from a Poisson distribution (mean = 9 in all groups), the type I
error rates for the Harrell–Davis estimator remain near 0.05, irrespective of
sample size (panel A in Figure 9). However, the quantile(type=7) is
conservative and the situation deteriorates with increasing sample size. In
the most extreme situation considered, when sampling from a beta-binomial
distribution with $r=1,s=9,nbin=10$, the type I error rates were lower for
both quantile estimators relative to the Poisson case (panel B in Figure 9),
or when sampling from a beta-binomial distribution with $r=9$ (not illustrated
here, but see notebook `sim_fp_apd.Rmd`). Although the situation got worse
with increasing sample size for both estimators, Harrell–Davis outperformed
quantile(type=7) in all situations.

Figure 9:  Type I error rates for the comparison of the deciles of
distributions of all pairwise differences: discrete distributions.

####  3.2.2 Power

Under normality, the ANOVA on means performed best, followed by ANOVA on
trimmed means and finally the bootstrap method (panel A in Figure 10). When
sampling from lognormal distributions, power was low for the ANOVAs relative
to the bootstrap method, and much more so when making inferences about means
(panel B in Figure 10). For the contaminated normal distributions, again the
ANOVA on means performed poorly, but now the ANOVA on trimmed means dominates
the bootstrap approach (panel C in Figure 10). If sampling from contaminated
lognormal distributions, now the bootstrap method is the most powerful (panel
D in Figure 10). In all these situations, the Harrell–Davis and
quantile(type=7) estimators gave very similar results. Finally, in the
presence of tied values, the ANOVA on means dominated the other approaches,
and the Harrell–Davis estimator led to higher power than the quantile(type=7)
estimator (panels E and F in Figure 10).

Figure 10:  Power results for the comparison of the deciles of distributions
of all pairwise differences. See details in Figure 5 caption.

##  4 An Illustration

Both methods are illustrated using data dealing with perceived health (PH)
among older adults (Clark, et al., 2011). The first factor consists of two
educational groups: those who did not complete high school and those who have
some college or technical training. The other two groups are based on a
measure of depressive symptoms (CESD). One group corresponds to participants
with a CESD score greater than 15, which is often taken to indicate mild
depression or worse. The other level consists of participants with CESD scores
less than or equal to 15. The four groups are defined like this:

  * •

$A_{1}B_{1}$ = lower education, lower CESD,

  * •

$A_{1}B_{2}$ = lower education, higher CESD,

  * •

$A_{2}B_{1}$ = higher education, lower CESD,

  * •

$A_{2}B_{2}$ = higher education, higher CESD.

Perceived health results are illustrated for the four groups in Figure 11A.
The figure was generated using the notebook `examples.Rmd`.

Figure 11:  Example: comparison of the deciles of perceived health. (A)
Perceived health scores in the 4 groups, with superimposed deciles indicated
by horizontal lines. Medians appear as thicker lines. (B-C) Comparisons
between B1 and B2 (lower/higher CESD) at each level of A. (D-E) Main effects.
(F) Interaction. Grey disks: decile differences plotted as a function of the
deciles in one group. Thick vertical lines mark 50% confidence intervals. Thin
vertical lines mark 95% confidence intervals. Vertical dashed lines mark
medians.

A 2x2 ANOVA on means returns these p values: main effect of A (education) =
0.001; main effect of B (depression) $<$ 0.0001; interaction = 0.09. Should we
conclude that we have failed to obtain sufficient evidence about the presence
of an interaction? This conclusion would be appropriate if the populations
were symmetric and differed only in central tendency. However, the plot of
marginals suggests differences in skewness and spread (Figure 11A). In keeping
with this observation, considering the deciles reveals a more complex picture,
as illustrated in panels B-F of Figure 11, with patterns of non-uniform group
differences. Panels B and C illustrate the shift functions comparing B1 and B2
at each level of A: the decile differences between two groups are plotted as a
function of the deciles in one group. Panels D and E illustrate the main
effects. The values along the x axis (A1 and B1) correspond to the deciles of
observations pooled across groups (for instance A1 = (A1B1, A1B2). Computing
the average of the deciles leads to very similar graphs. Finally, panel F
illustrates the interaction, which is highly non-linear, growing from the
first decile to the median, and then decreasing and reversing sign.

The output of the `decinter` R function that tests (1) is shown in Table 1. As
can be seen, the unadjusted p-values suggest that there is an interaction for
the .2-.6 quantiles. Shown are the adjusted p-values based on Hochberg’s
method (using the R function `p.adjust`) in order to underscore the practical
advantage of the Benjamini–Hochberg method. As indicated, no significant
difference is found at the .05 level using Hochberg’s method. But using the
Benjamini–Hochberg correction when testing (1), the adjusted p-values, when
dealing with the .3, .4 and .5 quantiles are all equal to .030.

Table 1: Result for perceived health when testing (1)

         Quant Est.Lev 1 Est.Lev 2        Dif     ci.low      ci.up p-value p.adj
     [1,]   0.1  2.137558  5.295364  -3.157806  -9.898611  3.9567186   0.389 0.775
     [2,]   0.2  2.334238  9.091485  -6.757246 -13.234037 -0.3524981   0.039 0.195
     [3,]   0.3  2.865681 11.685416  -8.819735 -15.196697 -2.2297824   0.010 0.070
     [4,]   0.4  2.670866 12.699564 -10.028698 -16.657294 -2.9629910   0.008 0.070
     [5,]   0.5  3.098901 13.724232 -10.625331 -17.541216 -3.2211028   0.009 0.070
     [6,]   0.6  4.539717 14.729066 -10.189350 -16.841309 -0.8789187   0.036 0.195
     [7,]   0.7  6.823277 11.374696  -4.551418 -14.407102  5.1928292   0.348 0.775
     [8,]   0.8  7.708862  5.545448   2.163415  -8.824631  9.6404128   0.775 0.775
     [9,]   0.9  6.201298  2.407586   3.793712  -4.813915 10.4323243   0.362 0.775

When testing (7), plots of the difference scores help provide perspective.
Figure 12A illustrates this point. For each level of education (A1 = lower
level; A2 = higher level), every participant with a lower CESD score was
compared to every participant with a higher CESD score. As previously noted,
when two distributions are identical, the distribution of $D$ is symmetric
about zero. Figure 12A suggests that for both groups the distribution of the
difference scores is shifted to the right, but with a stronger shift for group
2 (completed high school – lower panel in Figure 12A). A positive difference
indicates higher perceived health in not depressed participants relative to
depressed participants. For the second group, testing the hypothesis that the
median of the difference scores is zero, the estimate is 10.5, with a [4.2,
14.7] 95% confidence interval, $p=0$. In the first group the median is 4.21
[0, 8.4], $p=0.0765$.

Figure 12:  Example: comparison of the deciles of all pairwise differences of
perceived health. (A) Distributions of all pairwise differences between
participants with lower and higher CESD scores. Participants with lower levels
of education are shown at the top, participants with higher levels at the
bottom. The vertical lines mark the deciles, with a thicker line for the
median. (B). Interaction plot; see details in Figure 11 caption.

Table 2 summarizes the results when testing (7). Here, the .1, .25, .5, .75
and .9 quantiles are used, which is the default for the R function `iband`
that was used. Now the unadjusted p-values indicate an interaction in the
upper tails of the two distributions. For example, the estimates of the 0.75
quantile indicate that for the first level of the first factor (did not
complete high school), when comparing the not depressed group to the depressed
group, there is 25% chance of getting a difference between perceived health
values greater than 14.75, while for the second group there is 25% chance of
getting a difference greater than 22. But using the Hochberg adjustment shown
here, no significant difference is found at the .05 level and this remains the
case when using the Benjamini–Hochberg method.

Table 2: Result for Perceived health When Testing (7)

         Quant  Est.Lev 1  Est.Lev 2       Dif    ci.low      ci.up p-value p.adj
    [1,]  0.10 -15.380695 -14.229966 -1.150729  -8.58284  4.6709548   0.607 0.607
    [2,]  0.25  -6.342360  -2.168466 -4.173895 -10.51543  2.0383410   0.165 0.330
    [3,]  0.50   4.206226  10.466616 -6.260390 -12.44074  0.1001971   0.058 0.174
    [4,]  0.75  14.746649  22.041140 -7.294491 -12.55662 -1.5587962   0.013 0.065
    [5,]  0.90  24.529027  29.460025 -4.930998 -10.51350 -0.1045500   0.039 0.156

##  5 Concluding Remarks

The two new methods presented in this article help gain a deeper understanding
of where and by how much groups differ in a 2x2 factorial design. These
methods are implemented in the functions `decinter` to test (1) and `iband` to
test (7). The R functions `decinter` and `iband` can be found in the file
`Rallfun-v41.txt`, downloadable from https://osf.io/nvd59/quickfiles. These
functions, the R code for the simulation and the figures are available in the
reproducibility package for the article (Wilcox & Rousselet, 2023) and on
GitHub (Rousselet, 2023). The reproducibility package also contains Rcpp code
that is much faster to execute than the base R version (see notebook
`examples.Rmd`). By default, these functions use 2000 bootstrap samples and
correct for multiple comparison using the FDR correction from Benjamini &
Hochberg (1995), which outperformed the FWER correction of Hochberg (1988) in
our simulations. Even with the FDR correction, the two methods presented here
remain conservative, so it would be worthwhile to explore other correction
strategies. We have already considered several alternative methods, but they
do not improve matters (Benjamini & Yekutieli, 2001; Benjamini et al., 2006;
Blanchard & Roquain, 2008). Preliminary investigations of yet other approaches
to control the FWER, such as using a maximum statistic distribution (Nichols &
Holmes, 2002), have not revealed any method that could significantly improve
power in all situations, making recommendations difficult. While this issue
requires further work, for the moment we recommend to use the FDR correction
from Benjamini & Hochberg (1995) by default.

More generally, the obvious concern with comparing multiple quantiles is that
power might be negatively impacted due to controlling the FWE rate. But power
might also be negatively impacted when focusing on a single quantile simply
because other quantiles are being ignored. The example in the previous section
demonstrated the risk of drawing conclusions from a single measure of central
tendency or quantile. It is also worth keeping in mind that there is no free
lunch in inference: methods that can reveal more complex patterns in the data,
such as those proposed here, necessarily require larger sample sizes to reveal
where and by how much distributions differ. For applications to the deciles
considered here, given the increased uncertainty associated with the
estimation of the 1st and 9th deciles, we recommend sample sizes of at least
30, echoing earlier recommendations (Wilcox et al., 2014).

Several strategies are worth considering to boost power, starting with testing
more specific hypotheses involving only a subset of quantiles. For instance,
one could imagine a cross-validation approach in which a large dataset is
split between a discovery set and a testing set. Instead of testing each
quantile individually, one could also look for specific patterns across
quantiles, such as stochastic dominance, which is characterised by all
quantile differences having the same sign, or differences in spread, which are
characterised by monotonic trends across quantile differences (Rousselet et
al. 2017). Other obvious strategies are to consider different quantile
estimators and bootstrap methods. Preliminary investigations suggest that
using the quantile(type=8) estimator recommended by Hyndman & Fan (1996)
improves type I error rates and power relative to the quantile(type=7)
estimator for instance, but is still outperformed by the Harrell–Davis
estimator in all situations. As for parametric methods, Goldman & Kaplan
(2018) have proposed a fast and powerful extension of the Kolmogorov–Smirnov
test to compare all the quantiles of two independent groups (Kolmogorov, 1992;
Stephens, 1992). However, their approach assumes no tied values, and it is
unclear how it could be generalised to deal with interactions.

A bootstrap approach provides enough flexibility to deal with a variety of
experimental designs. As such, this work could be extended to mixed designs
with within- and between-subject factors. A covariate could be handled using
different strategies, including non-parametric methods with smoothers (Wilcox,
1997; Wilcox, 2021, chapter 12).

##  6 References

Akinshin, A. (2022). Trimmed Harrell-Davis quantile estimator based on the
highest density interval of the given width. Communications in Statistics -
Simulation and Computation. Online. DOI: 10.1080/03610918.2022.2050396

Benjamini, Y. & Hochberg, Y. (1995) Controlling the False Discovery Rate: A
Practical and Powerful Approach to Multiple Testing. Journal of the Royal
Statistical Society: Series B (Methodological), 57, 289–300.

Benjamini, Y., Krieger, A.M., & Yekutieli, D. (2006) Adaptive Linear Step-up
Procedures That Control the False Discovery Rate. Biometrika, 93, 491–507.

Benjamini, Y. & Yekutieli, D. (2001) The control of the false discovery rate
in multiple testing under dependency. The Annals of Statistics, 29, 1165–1188.

Blanchard, G. & Roquain, E. (2008) Two simple sufficient conditions for FDR
control. Electronic Journal of Statistics, 2, 963–992.

Bradley, J. V. (1978) Robustness? British Journal of Mathematical and
Statistical Psychology, 31, 144–152.
doi.org/10.1111/j.2044-8317.1978.tb00581.x

Cain, M., Zhang, Z. & Yuan, K.-H. (2017). Univariate and multivariate skewness
and kurtosis for measuring nonnormality: Prevalence, influence and estimation.
Behavioral Research, 49, 1716–1735. DOI 10.3758/s13428-016-0814-1.

Clark, F., Jackson, J., Carlson, M., Chou, C.-P., Cherry, B. J., Jordan-Marsh
M., Knight, B. G., Mandel, D., Blanchard, J., Granger, D. A., Wilcox, R. R.,
Lai, M. Y., White, B., Hay, J., Lam, C., Marterella, A. & Azen, S. P. (2011).
Effectiveness of a lifestyle intervention in promoting the well-being of
independently living older people: results of the Well Elderly 2 Randomised
Controlled Trial. Journal of Epidemiology and Community Health, 66, 782–790.
doi:10.1136/jech.2009.099754.

De Neve, J. & Thas, O. (2017). A Mann–Whitney type effect measure of
interaction for factorial designs. Communications in Statistics - Theory and
Methods, 46, 11243–11260. DOI: 10.1080/03610926.2016.1263739.

Doksum, K. A. & Sievers, G. L. (1976). Plotting with confidence: graphical
comparisons of two populations. Biometrika, 63, 421–434.
doi.org/10.2307/2335720.

Eddelbuettel, D. & Francois, R. (2011) Rcpp: Seamless R and C++ Integration.
J. Stat. Soft., 40, 1–18.

Gao, X. & Alvo, M. (2005). A nonparametric test for interaction in two-way
layouts. Canadian Journal of Statistics, 33, 529–543.

Goldman, M. & Kaplan, D. M. (2018). Comparing distributions by multiple
testing across quantiles or CDF values. Journal of Econometrics, 206, 143–166.

Harrell, F. E. & Davis, C. E. (1982). A new distribution-free quantile
estimator. Biometrika, 69, 635–640. https://doi.org/10.1093/biomet/69.3.635.

Hoaglin, D.C. (1985) Summarizing Shape Numerically: The g-and-h Distributions.
In Exploring Data Tables, Trends, and Shapes. John Wiley & Sons, Ltd, pp.
461–513.

Hochberg, Y. (1988) A sharper Bonferroni procedure for multiple tests of
significance. Biometrika, 75, 800–802.

Hommel, G. (1988). A stagewise rejective multiple test procedure based on a
modified Bonferroni test. Biometrika,75 , 383–386.

Hyndman, R. J. & Fan, Y. (1996) Sample quantiles in statistical packages,
American Statistician 50, 361–365. doi: 10.2307/2684934.

Kolmogorov, A. (1992) On the Empirical Determination of a Distribution
Function. In Kotz, S. & Johnson, N.L. (eds), Breakthroughs in Statistics:
Methodology and Distribution, Springer Series in Statistics. Springer, New
York, NY, pp. 106–113.

Liu, X., Song, Y. & Zhang, K. (2022). An exact bootstrap-based bandwidth
selection rule for kernel quantile estimators.Communications in Statistics -
Simulation and Computation, online. DOI: 10.1080/03610918.2022.2110595

Lombard, F. (2005). Nonparametric confidence bands for a quantile comparison
function. Technometrics, 47, 364–369.

Nichols, T.E. & Holmes, A.P. (2002) Nonparametric permutation tests for
functional neuroimaging: A primer with examples. Human Brain Mapping, 15,
1–25.

Navruz, G. & Özdemir, A. F. (2020). A new quantile estimator with weights
based on a subsampling approach. British Journal of Mathematical and
Statistical Psychology, 73, 506–521. https://doi.org/10.1111/bmsp.12198

Patel, K. M. & Hoel, D. G. (1973). A nonparametric test for interaction in
factorial experiments. Journal of the American Statistical Association, 68,
615–620. https://doi.org/10.2307/2284788

Rousselet, G. A. (2023) Code for article A Quantile Shift Approach To Main
Effects And Interactions In A 2-By-2 Design. GitHub.
https://github.com/GRousselet/decinter

Rousselet, G.A., Pernet, C.R., & Wilcox, R.R. (2017) Beyond differences in
means: robust graphical methods to compare two groups in neuroscience.
European Journal of Neuroscience, 46, 1738–1748.

Stephens, M.A. (1992) Introduction to Kolmogorov (1933) On the Empirical
Determination of a Distribution. In Kotz, S. & Johnson, N.L. (eds),
Breakthroughs in Statistics: Methodology and Distribution, Springer Series in
Statistics. Springer, New York, NY, pp. 93–105.

Wickham, H. (2016) Ggplot2: Elegant Graphics for Data Analysis, 2nd edn, Use
R! Springer International Publishing.

Wilcox, R.R. (1995) Comparing Two Independent Groups Via Multiple Quantiles.
Journal of the Royal Statistical Society. Series D (The Statistician), 44,
91–99.

Wilcox, R.R. (1997) ANCOVA based on comparing a robust measure of location at
empirically determined design points. British Journal of Mathematical and
Statistical Psychology, 50, 93–103.

Wilcox, R. (2022). Introduction to Robust Estimation and Hypothesis Testing
5th Edition. San Diego, CA: Academic Press

Wilcox, R. R., Erceg-Hurn, D., Clark, F. & Carlson, M. (2014). Comparing two
independent groups via the lower and upper quantiles. Journal of Statistical
Computation and Simulation, 84, 1543–1551. DOI: 10.1080/00949655.2012.754026

Wilcox, R. R. & Rousselet, G. A. (2023). A quantile shift approach to main
effects and interactions in a 2-by-2 design reproducibility package. figshare.
Online resource. https://doi.org/10.6084/m9.figshare.22927520.v1

Wilke, C.O. (2017) cowplot: Streamlined Plot Theme and Plot Annotations for
“ggplot2.”