11institutetext: Neurobiology Research Unit, Rigshospitalet 22institutetext:
Section of Biostatistics, University of Copenhagen 33institutetext: Molecular
Imaging Branch, National Institute of Mental Health (NIMH) 44institutetext:
Department of Computer Science, University of Copenhagen
44email<EMAIL_ADDRESS>
# A sensitivity analysis to quantify the impact of neuroimaging preprocessing
strategies on subsequent statistical analyses

Brice Ozenne 1122 0000-0001-9694-2956 Martin Nørgaard 3344 1111-2222-3333-4444
Cyril Pernet 11 2222–3333-4444-5555 Melanie Ganz 1144 2222–3333-4444-5555

###### Abstract

Even though novel imaging techniques have been successful in studying brain
structure and function, the measured biological signals are often contaminated
by multiple sources of noise, arising due to e.g. head movements of the
individual being scanned, limited spatial/temporal resolution, or other issues
specific to each imaging technology. Data preprocessing (e.g. denoising) is
therefore critical. Preprocessing pipelines have become increasingly complex
over the years, but also more flexible, and this flexibility can have a
significant impact on the final results and conclusions of a given study. This
large parameter space is often referred to as multiverse analyses. Here, we
provide conceptual and practical tools for statistical analyses that can
aggregate multiple pipeline results along with a new sensitivity analysis
testing for hypotheses across pipelines such as “no effect across all
pipelines” or “at least one pipeline with no effect”. The proposed framework
is generic and can be applied to any multiverse scenario, but we illustrate
its use based on positron emission tomography data.

###### Keywords:

neuroimaging preprocessing multiverse analyses

##  1 Introduction

Modern neuroimaging techniques have provided unique opportunities to measure
complex signaling pathways in the living human brain with the goal of
identifying reliable biomarkers of disease states and treatment outcomes [10].
Data arising from state-of-the-art neuroimaging techniques are, however, often
contaminated with noise confounds such as motion-related artefacts, affecting
both the spatial and temporal correlation structure of the data [9]. Carefully
designed preprocessing steps have been developed to remove unwanted noise
sources, but in the absence of a "ground truth" it remains a major challenge
to evaluate the impact of preprocessing choices on subsequent statistical
analyses and results. Over time, preprocessing pipelines (i.e. a set of
preprocessing steps) have become more complex and flexible, and this increase
in researcher degrees of freedom (termed multiverse analyses) has consistently
been shown to affect the outcomes of neuroimaging studies [3, 8, 2]. The most
common approach in the neuroimaging field is, to date, to use a single
pipeline and ignore the heterogeneity of preprocessing choices. This approach
not only makes abstraction of the multitude of possible results but is likely
also sub-optimal because the best pipeline is more often than not, study,
population or even subject dependent [8, 4]. More concerning, neuroscientists
might be tempted to “tune” the pipeline in order to obtain the most satisfying
results. This will generally lead to spurious and non-reproducible results
since the variability induced by the choice of pipeline is not independent
from the results. However, since it is neither realistic nor optimal to move
toward a single unified preprocessing pipeline, there is an urgent need for a
statistical framework allowing to explore results among many preprocessing
pipelines in a principled way.

The aim of this work is thus to provide a statistical framework that can
aggregate the evidence from multiverse analyses to produce conclusions robust
to the choice of the pipeline. More specifically, the present paper proposes a
statistical sensitivity analysis providing:

  1. (i)

visualizations of the heterogeneity of several preprocessing pipelines

  2. (ii)

estimation of a global effect across all preprocessing pipelines

  3. (iii)

quantification of the proportion of pipelines with evidence for an effect

  4. (iv)

a statistical framework for testing hypotheses across pipelines such as “no
effect across all pipelines”, “at least one pipeline with no effect”

The corresponding software can be found at *anonymous*, including all code to
reproduce all simulations and figures.

##  2 Materials and experimental settings

###  2.1 Data

We use two different data sources in our analyses: in sillico and real data.
For the in sillico data, different noise structures were chosen to reflect
different configurations of pipelines and the sample size was varied to
encompass small to larger scale clinical studies. In the real data analysis,
we mimic how real neuroimaging studies compare an intervention to a reference
measurement. Pipelines were selected independently of the intervention data
using the healthy/placebo arm of [6] using results from [8]. The proposed
sensitivity analysis was illustrated on the intervention arm of the study
where the follow-up value was compared to a reference value taken from a
normative serotonergic atlas [1].

To simulate in sillico data, we consider the simple case of a single brain
measurement ($R=1$), with a single binary exposure ($P=1$) following a
Bernoulli distribution with parameter $\pi=0.5$ (i.e. two balanced groups) and
no covariates ($C=0$). Latent $Y$ values for the brain measurement are
simulated using a normal distribution with variance 1 and mean $\beta$ times
the exposure $X$ values where $\beta=0$ (null hypothesis) or $\beta=0.5$
(alternative hypothesis). The observed $Y$ is simulated for each pipeline
adding pipeline specific noise to the latent $Y$. This noise is simulated
using a multivariate normal distribution with mean 0 and variance
$\Sigma_{\text{scenario 1}}$, $\Sigma_{\text{scenario 2}}$, and
$\Sigma_{\text{scenario 3}}$ depending on the scenario. In scenario 1, we
simulated many pipelines ($J=20$) with correlated homoscedastic noise, in
scenario 2 a few pipelines ($J=6$) with uncorrelated heteroscedastic noise,
and in scenario 3 many pipelines ($J=20$) with correlated heteroscedastic
noise. Scenario 2 and 3 included one pipeline with high signal to noise ratio
(SNR), i.e. low variance, and many pipelines with low SNR, i.e., high
variance. The sample size was varied from $n=10$ to $n=500$ in each group,
such that the smallest sample size was well below the number of parameters
($2\times J$ mean parameters, $J^{2}$ variance-covariance parameters) in
scenario 1 and 3 and the asymptotic regime was reached for the largest sample
size. We generate 10,000 datasets per scenario and sample size - this provides
sufficient precision about the mean, standard deviation, and rejection rate to
neglect the Monte Carlo uncertainty. This data will be used to assess the
large sample size properties of the procedure (bias, relative efficacy, type 1
error control) in finite samples.

To illustrate the use of the proposed sensitivity analysis on real data, we
utilize neuroimaging results from a placebo-controlled, double-blinded,
clinical study [6]. The study was registered and approved by the ethics
committee for the capital region of Copenhagen (protocol-ID: H-2-2010-108) and
registered as a clinical trial: www.clinicaltrials.gov under the trial ID
NCT02661789. All subjects provided written informed consent prior to
participation, in accordance with The Declaration of Helsinki II. The aim of
the study was to assess the association between the emergence of depressive
symptoms and change in cerebral serotonin transporter (SERT) availability
following a hormonal treatment ($p=1$). Data is available from the CIMBI
database (www.cimbi.dk) upon request. It consists of structural Magnetic
Resonance Imaging (MRI) and Positron Emission Tomography (PET) imaging data
for $60$ healthy females who underwent a baseline scan, received either
Placebo ($n=30$) or a GnRHa implant intervention ($n=30$), and participated in
a follow-up scan. SERT availability estimates were extracted for each subject
for 28 subcortical and cortical regions, and averaged across hemispheres,
producing a final sample of $R=14$ regions per subject and pipeline. These
regions (amygdala, thalamus, putamen, caudate, anterior cingulate cortex,
hippocampus, orbital frontal cortex, superior frontal cortex, occipital
cortex, superior temporal gyrus, insula, inferior temporal gyrus, parietal
cortex, and entorhinal cortex) were chosen because they cover the entire
brain, and many are target regions in published serotonin transporter (SERT)
PET studies. No covariates where consider and ($C=0$).

###  2.2 Data preprocessing

Five preprocessing steps were used to curate the data and estimate the SERT
availability (outcome measure). These steps include, motion correction
(with/without), co-registration (4 options), delineation of volumes of
interest (3 options), partial volume correction (4 options), and kinetic
modeling for quantification of SERT availability (MRTM, SRTM, Non-invasive
Logan and MRTM2). More information about the preprocessing choices can be
found in [8]. The combination of individual preprocessing steps leads to a
number of $J=2\times 3\times 4^{3}=384$ possible combinations.

###  2.3 Notation and assumptions

In the following section, we use generic notations to describe the proposed
sensitivity analysis. We are interested in relating $R$ brain measurements
($\boldsymbol{Y}=\left(Y_{1},\ldots,Y_{R}\right)$) to $P$ exposures or
treatments ($\boldsymbol{X}=\left(X_{1},\ldots,X_{P}\right)$) accounting for
$C$ covariates ($\boldsymbol{W}=\left(W_{1},\ldots,W_{C}\right)$). We consider
a set of $J$ pipelines used to preprocess the neuroimaging data. For a given
pipeline $j\in\\{1,\ldots,J\\}$ we fit a statistical model with parameters
$\theta_{j}$ that relates $\boldsymbol{Y}$ processed by pipeline $j$,
$\boldsymbol{X}$, and $\boldsymbol{W}$. We then obtain from this model an
estimate $\widehat{\psi}_{k}$ of the effect of interest (denoted $\psi$). In
our real life example we use, for each pipeline, a paired t-test to compare
the observed change in SERT availability to an atlas value so for
$j\in\\{1,\ldots,J\\}$,
$\widehat{\theta}_{j}=(\widehat{\psi}_{j},\widehat{\sigma}^{2}_{j})$ where
$\widehat{\psi}_{j}$ is the empirical mean and $\sigma^{2}_{j}$ the empirical
variance of the change in SERT availability (processed with pipeline $j$)
between baseline of follow-up.

We make the following working assumptions: first the observed data
$\left(\mathcal{O}_{i}\right)_{i\in\\{1,\ldots\\}}=\left(\boldsymbol{y}_{i},\boldsymbol{x}_{i},\boldsymbol{w}_{i}\right)_{i\in\\{1,\ldots,n\\}}$
correspond to independent and identically distributed replicates of
$\left(\boldsymbol{Y},\boldsymbol{X},\boldsymbol{W}\right)$. Second we have
chosen a set of reasonable pipelines, meaning that the estimated effects
$\widehat{\psi}_{1},\ldots,\widehat{\psi}_{J}$ found in the follow-up
statistical analysis will converge to the right value $\psi$ as the sample
size increases. This set can include pipelines distorting the signal
$\boldsymbol{Y}$ (e.g. adding a fixed value) if that has no consequence,
asymptotically, on the estimated effect (the mean change is not biased as the
added value cancels out when substracting the baseline value to the follow-up
value). Finally, when considering asymptotic results we will consider a fixed
number of pipelines and let the sample size $n$ increase to infinity.

##  3 Proposed sensitivity analysis

To be able to draw conclusions across pipelines, we not only need the result
of each pipeline but also some information about how they relate. If all
pipelines were equally reliable and equally different, we would weight each
pipeline equally. If there exists one independent pipeline and a block of
correlated pipelines, all equally reliable, then we would weight the
independent pipeline more compared with other pipelines. By treating pipelines
as black boxes, we can investigate their relation in terms of each
observations influence on the estimated effects across pipelines,
$\widehat{\boldsymbol{\psi}}=\left(\widehat{\psi}_{1},\ldots,\widehat{\psi}_{J}\right)$.
This relation is fully characterized by the joint distribution of the effects.
Once estimated, we can extract summaries of this distribution, e.g. an average
value, and carry out statistical tests, e.g. testing the compatibility between
the observations and the joint distribution that would have been observed
under a specific hypothesis.

###  3.1 Estimating the joint distribution across pipelines

The joint distribution could be obtained using a multivariate model, e.g.
modelling data from all pipelines at once using a mixed model. Because of the
complexity of the dependency among pipelines, this is, however, rarely
feasible with the available sample size. Instead, and this matches common
practice, we perform the same analysis separately for each pipeline and obtain
a vector of estimated associations $\widehat{\boldsymbol{\psi}}$ with their
standard errors
$\boldsymbol{\sigma}_{\widehat{\psi}}=\left(\sigma_{\widehat{\psi}_{1}},\ldots,\sigma_{\widehat{\psi}_{J}}\right)$.
Using tools from the semi-parametric theory (see [7] and [12] for more
details), we can approximate the influence of each observation on the estimate
by a random variable called the influence function, denoted
$\varphi_{\widehat{\psi}_{j}}$ for pipeline $j$, and satisfying:

$\displaystyle\sqrt{n}\left(\widehat{\psi}_{j}-\psi\right)=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\varphi_{\widehat{\psi}_{j}}(\mathcal{O}_{i})+o_{p}(1)$

where $o_{p}(1)$ denotes a residual term that convergences toward zero in
probability as the sample size tends to infinity. As shown in appendix A,
$\varphi_{\widehat{\psi}_{j}}$ has a simple expression when
$\widehat{\psi}_{j}$ is the empirical mean or an element of a maximum
likelihood (ML) estimator. Since this decomposition applies to all pipelines,
we get from the multivariate central limit theorem that the joint distribution
of the estimates is asymptotically multivariate normal. It has mean $\psi$ and
its variance-covariance, denoted $\Sigma_{\widehat{\boldsymbol{\psi}}}$, is
the same as the one of
$\boldsymbol{\varphi}_{\widehat{\boldsymbol{\psi}}}=\left(\varphi_{\widehat{\psi}_{1}},\ldots,\varphi_{\widehat{\psi}_{J}}\right)$
divided by $n$. Note that with limited number of observations, typically when
$n<J$, the estimated variance-covariance
$\widehat{\Sigma}_{\widehat{\boldsymbol{\psi}}}$ based on the estimated
influence function is not guaranteed to be positive definite.

###  3.2 Testing hypotheses across pipelines

The global null hypothesis “no effect across all pipelines” can be tested
using a max-test approach: more extreme realizations would correspond to
larger values of the maximum statistic
$\widehat{\boldsymbol{t}}_{\max}=\max(|\widehat{t}_{1}|,..,|\widehat{t}_{J}|)$
where $|.|$ denotes the absolute value and
$\widehat{t}_{j}=\frac{\widehat{\psi}_{j}}{\sigma_{\widehat{\psi}_{j}}}$. A
p-value may therefore be computed by integrating the joint density under the
null hypothesis outside of the domain
$\mathcal{D}(\widehat{\boldsymbol{t}}_{\max})=[-\widehat{\boldsymbol{t}}_{max},\widehat{\boldsymbol{t}}_{max}]^{\otimes
J}$ (see figure A in appendix B). Here, we use the notation $[a,b]^{\otimes
J}=\prod_{j=1}^{J}[a,b]$ that represents the Cartesian product between $J$
intervals $[a,b]$. The value $t_{c}$, such that the integral outside
$\mathcal{D}(t_{c})$ equals $\alpha$, provides a critical threshold for the
estimated test statistics
$\left(|\widehat{t}_{1}|,..,|\widehat{t}_{J}|\right)$. This threshold can also
be used to derive confidence intervals.

The null hypothesis “at least one pipeline with no effect” is an intersection
union test. As such it can be rejected if and only if all the un-adjusted
p-values relative to each pipeline are below $\alpha$ or equivalently if the
largest un-adjusted p-value is below $\alpha$.

The proportion of pipelines where there is evidence for an effect $\eta$ can
be estimated as
$\tilde{\eta}=\frac{1}{J}\sum_{j=1}^{J}\mathds{1}_{\left|\widehat{t}_{j}\right|<t_{c}}$
where $\mathds{1}_{.}$ denotes the indicator function. One drawback with this
non-parametric estimator is that it is a non-smooth function of
$\widehat{t}_{j}$, making the associated uncertainty difficult to evaluate.
Instead, one can use a parametric approach, assuming normally distributed test
statistics:

$\displaystyle\widehat{\eta}=\frac{1}{J}\sum_{j=1}^{J}1-\mathbb{P}\left[-t_{c}<\widehat{t}_{j}<t_{c}\right]=1-\frac{1}{J}\sum_{j=1}^{J}\Phi\left(t_{c}-\widehat{t}_{j}\right)-\Phi\left(-t_{c}-\widehat{t}_{j}\right)$

which is a smooth (but complex) function of the model parameters. Here $\Phi$
refers to the cumulative distribution function of a standard normal
distribution. The uncertainty about the estimator can therefore be derived
using a non-parametric bootstrap or a delta method where
$Var(\widehat{\eta})=\frac{\partial\widehat{\eta}}{\partial\Theta}\Sigma_{\widehat{\Theta}}{\frac{\partial\widehat{\eta}}{\partial\Theta}}^{\intercal}$
where $\Theta=\left(\theta_{j}\right)_{j\in\\{1,\ldots,J\\}}$ is the set of
parameters of the statistical model across pipelines.

###  3.3 Visualizing the heterogeneity across pipelines

We have derived that the estimated associations are, asymptotically, normally
distributed. They can therefore be summarized by their expectation and
variance-covariance matrix (i.e. standard errors and correlation matrix). We
suggest two graphical displays to visualize the heterogeneity of the results
across pipelines:

  * •

A heatmap of the estimated correlation among estimates, obtained by converting
$\widehat{\Sigma}_{\widehat{\boldsymbol{\psi}}}$ into a correlation matrix.
111Re-ordering the pipelines may be useful to better visualize blocks of
pipelines that are especially correlated. This can for instance be performed
by converting R to a dissimilarity matrix and then use hierarchical
clustering.

  * •

A forest plot displaying
$\left(\widehat{\psi}_{1},\ldots,\widehat{\psi}_{J}\right)$ and
$\left(\widehat{\sigma}_{\widehat{\psi}_{1}},\ldots,\widehat{\sigma}_{\widehat{\psi}_{J}}\right)$
through the confidence intervals, possibly using the previously established
re-ordering.

###  3.4 Estimating a global effect across pipelines

Several methods can be used for estimating a global effect across pipelines. A
naive method would be to compute the mean of the estimated associations:

$\displaystyle\widehat{\Psi}_{\text{average}}=\frac{1}{J}\sum_{j=1}^{J}\widehat{\psi}_{j}$

This estimator will, however, not be efficient if some pipelines lead to more
precise estimates, i.e.
$\left(\sigma_{\widehat{\psi}_{1}},\ldots,\sigma_{\widehat{\psi}_{J}}\right)$
are not equal. Intuitively, we would like to pool the estimates with weights
inversely proportional to the standard errors such that we put more weight on
precise estimates:

$\displaystyle\widehat{\Psi}_{\text{pool-
se}}=\sum_{j=1}^{J}w^{se}_{j}\widehat{\psi}_{j}\text{ where
}w^{se}_{j}=\frac{1/\sigma^{2}_{\widehat{\psi}_{j}}}{\sum_{j=1}^{J}1/\sigma^{2}_{\widehat{\psi}_{j}}}$

However, we also need to take into account the correlation between estimates.
Indeed, perfectly correlated estimates should weight as if there was only one
estimate. To do so, we can use the following GLS estimator of the global
effect:

$\displaystyle\widehat{\Psi}_{\text{GLS}}=\left({\mathbf{1}}^{\intercal}\widehat{\Sigma}_{\widehat{\Psi}}^{-1}\mathbf{1}\right)^{-1}{\mathbf{1}}^{\intercal}\widehat{\Sigma}_{\widehat{\Psi}}^{-1}\widehat{\boldsymbol{\psi}}=\sum_{j=1}^{J}w^{GLS}_{j}\widehat{\psi}_{j}$
(1)

where $\mathbf{1}$ is a column vector filled with ones. Equation 1 can be
shown to be equivalent to performing a spectral decomposition of
$\widehat{\Sigma}_{\widehat{\boldsymbol{\psi}}}$, and use the eigenvectors to
combine estimates into independent components that can be pooled according to
weights proportional to the eigenvalues (appendix C.1). This is used when
$\Sigma_{\widehat{\boldsymbol{\psi}}}$ is singular to compute
$\widehat{\Psi}_{\text{GLS}}$, by restricting the spectral decomposition to
the eigenvalues above a given threshold ($\epsilon=10^{-10}$). In the simple
case where $R=1$, $p=1$, $C=0$, brain measurements are jointly normally
distributed, $X$ is binary, and there is no missing data, the GLS estimator
can be shown to be asymptotically efficient (appendix C.2).

Cox and colleagues [5] studied a similar estimator when $J=2$ and found that,
under unequal variance, the global estimate can be outside of the range of the
(pipeline specific) estimates. As a remedy to this unpleasant behavior, we
propose a constrained GLS estimator, denoted
$\widehat{\Psi}_{\text{constrained GLS}}$, which constrains the weight of each
estimate to be at most 1 in absolute value and ensures that the sum of the
weights is 1:

$\displaystyle w_{j}^{\text{constrained
GLS}}=\frac{w_{j}^{GLS}}{\kappa+\max_{j\in\\{1,\ldots,J\\}}|w_{j}^{GLS}|}+\frac{1}{J}\left(1-\frac{1}{\kappa+\max_{j\in\\{1,\ldots,J\\}}|w_{j}^{GLS}|}\right)$

where $\kappa$ is chosen to satisfy the constraints. Note that if some of the
pipelines may induce some bias, the previous approaches will propagate this
bias and therefore be unsatisfactory. Systematic differences between pipelines
can be investigated by comparing the estimates between pipelines, e.g.
$\widehat{\psi}_{j}-\widehat{\psi}_{j^{\prime}}$ and using
$\widehat{\Sigma}_{\widehat{\boldsymbol{\psi}}}$ to obtain the corresponding
uncertainty
$\mathbb{V}ar\left[\widehat{\psi}_{j}-\widehat{\psi}_{j^{\prime}}\right]=\mathbb{V}ar\left[\widehat{\psi}_{j}\right]+Var\left[\widehat{\psi}_{j^{\prime}}\right]-2\mathbb{C}ov\left(\widehat{\psi}_{j},\widehat{\psi}_{j^{\prime}}\right)$.

##  4 Results

###  4.1 Simulation results - pooled estimator

We compare the performance of the proposed estimators when using unbiased
pipelines on simulated data. For each dataset $\Sigma_{\text{scenario 1}}$,
$\Sigma_{\text{scenario 2}}$, and $\Sigma_{\text{scenario 3}}$, we computed
the four previously described estimators of the global effect,
$\widehat{\Psi}_{\text{average}}$, $\widehat{\Psi}_{\text{pool-se}}$,
$\widehat{\Psi}_{\text{GLS}}$, and $\widehat{\Psi}_{\text{constrained GLS}}$.
P-values for $\widehat{\Psi}_{\text{pool-se}}$, $\widehat{\Psi}_{\text{GLS}}$,
and $\widehat{\Psi}_{\text{constrained GLS}}$ were computed neglecting the
uncertainty of the weights ($w_{j}^{\text{se}}$, $w_{j}^{\text{GLS}}$,
$w_{j}^{\text{constrained GLS}}$).

Weights: Based on $\Sigma_{\text{scenario 1}}$, $\Sigma_{\text{scenario 2}}$,
and $\Sigma_{\text{scenario 3}}$, we can compute how
$\widehat{\Psi}_{\text{average}}$, $\widehat{\Psi}_{\text{pool-se}}$, and
$\widehat{\Psi}_{\text{GLS}}$ (or $\widehat{\Psi}_{\text{constrained GLS}}$)
would weight the results from each pipeline if the variance-covariance matrix
of the pipeline estimates was known (Figure 1). In scenario 1,
$\widehat{\Psi}_{\text{average}}$ and $\widehat{\Psi}_{\text{pool-se}}$ would
provide equal weight to all pipelines with a weight of 5%, while
$\widehat{\Psi}_{\text{GLS}}$ would weight by 2% the correlated pipelines and
by 14% the uncorrelated pipelines (in this paragraph all weights are rounded
for readability). In scenario 2, $\widehat{\Psi}_{\text{average}}$ equally
weights all pipelines by 17% while $\widehat{\Psi}_{\text{pool-se}}$ and
$\widehat{\Psi}_{\text{GLS}}$ would use the following weights: 8%, 82%, 4%,
3%, 2%, 1%, favoring the high SNR pipeline. In scenario 3,
$\widehat{\Psi}_{\text{average}}$ would equally weight all pipelines by 5%,
$\widehat{\Psi}_{\text{pool-se}}$ would weight each correlated pipeline by
3.8% and the remaining pipelines by 38% (high SNR pipeline), 2%, 1.3%, 1%,
0.6%, while $\widehat{\Psi}_{\text{GLS}}$ would weight by 0.6% the correlated
pipelines and by 81% (high SNR pipeline), 4%, 3%, 2%, 1% the remaining
pipelines.

Figure 1: Large sample weights used by $\widehat{\Psi}_{\text{average}}$,
$\widehat{\Psi}_{\text{pool-se}}$, and the GLS estimators
($\widehat{\Psi}_{\text{GLS}}$ or $\widehat{\Psi}_{\text{constrained GLS}}$)
to combine the pipeline specific estimates in scenario 3. Weights relative to
the correlated pipeline are shown in shades of gray (first fifteen blocks);
weights relative to the independent pipeline are shown using rainbow colors
(last five blocks). In the legend the variance for a given pipeline effect is
indicated in parenthesis.

Estimate and bias: No significant bias was found. Under the null, the
proportion of pipelines with evidence for an effect ranged between 5% and 7%
and was relatively stable with sample size. Under the alternative, the
proportion of pipelines increased from 6.3% in scenario 1 to 98% in scenario 1
when the sample size increased from $n=20$ to $n=500$. A similar behavior was
observed in the other scenarios (up to 63% in scenario 2 and up to 79% in
scenario 3).
TODO: add random intercept model

Efficiency: Similar results were obtained for $\beta=0$ and $\beta=0.5$, so we
will only discuss the former case (Figure 2). In scenario 1,
$\widehat{\Psi}_{\text{average}}$ and $\widehat{\Psi}_{\text{pool-se}}$ showed
similar empirical variance, whereas $\widehat{\Psi}_{\text{GLS}}$ showed
higher empirical variance in very small samples (+159% for $n=10$, +13% for
$n=25$) and lower empirical variance at larger sample sizes (e.g. -12% for
$n=500$). In scenario 2, $\widehat{\Psi}_{\text{pool-se}}$ and
$\widehat{\Psi}_{\text{GLS}}$ showed similar empirical variance (slightly
higher for $\widehat{\Psi}_{\text{GLS}}$ in low sample sizes and slightly
lower afterwards), both smaller compared to $\widehat{\Psi}_{\text{average}}$,
e.g. 10% for $n=10$ and -24% for $n=500$ for the GLS estimator. Results in
scenario 3 were similar to scenario 1 up to a larger decrease in variance in
large samples for GLS estimators (-28% for n=500).
$\widehat{\Psi}_{\text{constrained GLS}}$ was similar to
$\widehat{\Psi}_{\text{GLS}}$ but with better small sample properties (at most
+4% in standard error compared to $\widehat{\Psi}_{\text{average}}$).

Type 1 error: Across all scenarios, the type 1 errors for
$\widehat{\Psi}_{\text{average}}$ and $\widehat{\Psi}_{\text{pool-se}}$ were
well controlled except in very small samples where small deviations from the
nominal level were observed (maximum of 7% for
$\widehat{\Psi}_{\text{average}}$ and 8% for $\widehat{\Psi}_{\text{pool-
se}}$). When neglecting the uncertainty about the weights, the type 1 error
control of $\widehat{\Psi}_{\text{GLS}}$ and
$\widehat{\Psi}_{\text{constrained GLS}}$ were only controlled for large
samples, i.e. for $n=500$; large type 1 error rates were found in very small
samples (>10%).

Figure 2: Upper panel: empirical standard deviation of the estimated common
effect for each estimator, scenario, and sample size under the null
hypothesis. For large samples the green line (GLS) is covered by the blue line
(constrained GLS) in all scenarios. In scenario 1 and 3, the red line (pool-
average) is covered by the yellow line (pool se). In scenario 1, the GLS
estimators have the same or higher standard deviation as the average for small
samples and lower standard deviation (approx. 12% lower) for large samples.
Fig. X Lower panel: rejection rate under the null hypothesis (i.e. type 1
error). Shaded area represents the Monte Carlo uncertainty.

###  4.2 Real-world application

We illustrate the statistical sensitivity analysis with data described in
section 2.1 in which the SERT availability was assessed after a drug
intervention and compared to normative values. We start by studying the
behavior of the four statistical estimators (pooled GLS, pooled constrained
GLS, pooled average and pooled SE) to estimate a common effect across
pipelines for the given null hypothesis, reported as a forest plot in Figure 3
(left panel). The dashed vertical line in Figure 3 represents the normative
value, and all horizontal error bars represent the estimated effect (mean and
95% CI) for a given pipeline and estimator. Across a reasonable set of
preprocessing pipelines, 3 of the 8 selected reject the null hypothesis (as
indicated by the non-overlapping CI with the normative value) with estimated
percent differences between groups ranging between -9% (pipeline 6) and 2.5%
(pipeline 1). The pooled constrained GLS, pooled SE and pooled average, all
fail to reject a common effect across pipelines, whereas the GLS estimator
rejects the null hypothesis. However, when inspecting the pipeline weights for
the GLS estimator for this data, it assigned a very high weight to four
pipelines (i.e. weight above 1 in absolute values), leading to an unreliable
estimate as illustrated by the large standard deviation found in the
simulation study for low sample sizes (Figure 2). The constrained GLS
estimator did not exhibit this problem and had weights between -0.79 and 1.

Figure 3: Left panel: forest plot of the estimated SERT availability in the
amygdala for the intervention group (point and full line) vs. the normative
values (dashed line) for each pipeline and the four proposed pooled
estimators.
Fig. X Right panel: correlation of the estimated SERT availability between
pipelines.
Pipeline 1: with motion correction (MC), boundary based registration (BBR)
using the time-weighted average PET image (twa), and MRTM2 as kinetic modeling
choice.
Pipeline 2: MC, normalized mutual information registration (NMI) using twa,
MRTM2.
Pipeline 3: MC, BBR using the average PET image, and MRTM2.
Pipeline 4: MC, NMI using the average PET image, and MRTM2.
Pipeline 5: MC, BBR_twa, and MRTM.
Pipeline 6: no motion correction (nMC), BBR_twa, and MRTM.
Pipeline 7: nMC, BBR_twa, and MRTM2.
Pipeline 8: MC, BBR_twa, and SRTM.

Figure 3 (right panel) shows a heatmap for the estimated correlations across
preprocessing pipelines, ranging from 0.68 to 1. Pipelines 1-4 show a very
high correlation with each other (only varying the registration choices),
whereas pipelines 5, 6 and 8 are equally correlated with each other (different
kinetic models), and pipeline 7 is somewhat in between (no motion correction).
The heatmap captures important differences in the correlation structure
between pipelines, suggesting that not all pipelines perform similarly with
moderate levels of unexplained variance. Pipelines 4 to 8 exhibit numerically
smaller correlation compared to pipelines 1 to 4 (and similar variance), which
explains why GLS estimators produce estimates with lower numerical values
compared to pooled estimators ignoring the correlation, as GLS estimators
assign more weight to the last four pipelines.

Notably, no correction for multiple comparisons was carried out across regions
and pipelines at this point. The rationale for not including this (as should
otherwise always be carried out) is that we wanted to make our analysis as
comparable as possible to the Neuroimaging Analysis Replication and Prediction
Study [2], where each participating institution analyzed the data using their
own established pipeline and tested only a single region in a hypothesis-
driven fashion.

##  5 Discussion

Looking first at the simulated data results, we observe that, asymptotically,
$\widehat{\Psi}_{\text{GLS}}$ is the best estimator, leading to more precise
estimates than $\widehat{\Psi}_{\text{average}}$ or
$\widehat{\Psi}_{\text{pool-se}}$. It is, however, not suited for small sample
sizes or large numbers of pipelines, as it exhibits large variance. We
proposed an alternative estimator $\widehat{\Psi}_{\text{constrained GLS}}$
which can be seen as regularizing the GLS estimator toward the empirical
average in small samples. Even though this considerably improves the small
sample performance of the estimator, other regularization approaches may lead
to further gain. An alternative approach could be to regularize the estimated
variance-covariance matrix between pipeline specific estimates, e.g. using
graphical lasso. However, this is challenging since the correlation structure
among pipelines is typically complex and non-sparse. The proposed method to
quantify the uncertainty of $\widehat{\Psi}_{\text{constrained GLS}}$ is
numerically fast, but unreliable in small samples or with a large number of
pipelines. However, these latter points are beyond the scope of this paper,
and is left for future work. In the real-world application, we observed that
all estimators could be readily applied, and three of them performed as
expected based on the results from individual pipelines. Only the
$\widehat{\Psi}_{\text{GLS}}$ estimator performed differently and was the only
one rejecting the null hypothesis hinting at an effect across pipelines. This
is though, due to the estimator not being able to fit the weights properly in
small sample sizes. Since PET neuroimaging studies are rarely beyond sample
sizes of n>50, other estimators should be used. We recommend using
$\widehat{\Psi}_{\text{pool-se}}$ or $\widehat{\Psi}_{\text{constrained
GLS}}$: the latter when pipelines are not similarly related and the sample
size is moderate to large, otherwise the former.

The framework that we are proposing is not without limitations. The underlying
assumption of combining results across pipelines in our analysis is that all
pipelines are unbiased. This can however not always be guaranteed. Alternative
approaches (e.g. STAPLE [14]) could be used to reduce the bias of the pooled
estimators by assuming a majority of unbiased pipelines or identifying
clusters of pipelines and pooling pipeline-specific estimates within clusters.

##  6 Conclusion

In this work, we have developed a statistical sensitivity analysis that can
quantify the impact of different preprocessing choices on subsequent
statistical analyses. As has been reported in previous studies, we observe
that the influence of preprocessing pipelines on subsequent statistical
analysis can be quite large. Hence, we provide tools for statistical analyses
that can aggregate multiple analyses of the same data. We introduce four
statistical estimators, $\widehat{\Psi}_{\text{average}}$,
$\widehat{\Psi}_{\text{pool-se}}$, $\widehat{\Psi}_{\text{GLS}}$, and
$\widehat{\Psi}_{\text{constrained GLS}}$ to combine the pipeline specific
estimates. This enables testing hypotheses across pipelines, such as “no
effect across all pipelines” or “at least one pipeline with no effect”. The
proposed framework is generic and can be applied to any imaging modality.

## References

  * [1] Beliveau, V., et al.: A high-resolution in vivo atlas of the human brain’s serotonin system. Journal of Neuroscience 37(1), 120–128 (2017)
  * [2] Botvinik-Nezer, R., et al.: Variability in the analysis of a single neuroimaging dataset by many teams. Nature 582(7810), 84–88 (2020)
  * [3] Carp, J.: On the plurality of (methodological) worlds: estimating the analytic flexibility of fMRI experiments. Frontiers in neuroscience 6, 149 (2012)
  * [4] Churchill, N.W., et al.: An automated, adaptive framework for optimizing preprocessing pipelines in task-based functional MRI. PLoS ONE 10(7), 1–25 (2015)
  * [5] Cox, M., Eiø, C., Mana, G., Pennecchi, F.: The generalized weighted mean of correlated quantities. Metrologia 43(4), S268 (2006)
  * [6] Frokjaer, V.G., et al.: Role of serotonin transporter changes in depressive responses to sex-steroid hormone manipulation: a positron emission tomography study. Biological psychiatry 78(8), 534–543 (2015)
  * [7] Kennedy, E.H.: Semiparametric theory. arXiv:1709.06418 (2017)
  * [8] Nørgaard, M., , et al.: Optimization of preprocessing strategies in positron emission tomography (PET) neuroimaging: a [11C] DASB PET study. Neuroimage 199, 466–479 (2019)
  * [9] Nørgaard, M., Ozenne, B., Svarer, C., Frokjaer, V.G., Schain, M., Strother, S.C., Ganz, M.: Preprocessing, prediction and significance: Framework and application to brain imaging. International Conference on Medical Image Computing and Computer-Assisted Intervention pp. 196–204 (2019)
  * [10] Poldrack, R.A., et al.: Establishment of Best Practices for Evidence for Prediction: A Review. JAMA Psychiatry 77(5), 534–540 (2020)
  * [11] Stenbæk, D.S., Fisher, P.M., Ozenne, B., Andersen, E., Hjordt, L.V., McMahon, B., Hasselbalch, S.G., Frokjaer, V.G., Knudsen, G.M.: Brain serotonin 4 receptor binding is inversely associated with verbal memory recall. Brain and behavior 7(4), e00674 (2017)
  * [12] Tsiatis, A.A.: Semiparametric theory and missing data. Springer (2006)
  * [13] Van der Vaart, A.W.: Asymptotic statistics, vol. 3. Cambridge university press (2000)
  * [14] Warfield, S.K., Zou, K.H., Wells, W.M.: Simultaneous truth and performance level estimation (STAPLE): An algorithm for the validation of image segmentation. IEEE Transactions on Medical Imaging 23(7), 903–921 (2004). https://doi.org/10.1109/TMI.2004.828354

##  Appendix A Influence function

In the univariate case (i.e. $R=P=1$) without covariate ($C=0$), the effect of
interest $\psi$ may be the Pearson’s correlation coefficient:

$\displaystyle\psi=\frac{\mathbb{E}[XY]-\mathbb{E}[X]\mathbb{E}[Y]}{\sqrt{\mathbb{E}[X^{2}]-E[X]^{2}}\sqrt{\mathbb{E}[Y^{2}]-\mathbb{E}[Y]^{2}}}$

when X is continuous or the mean difference when X is binary.

$\displaystyle\psi=\mathbb{E}[Y|X=1]-\mathbb{E}[Y|X=0]$

Here $\mathbb{E}[.]$ denotes the expectation and $\mathbb{E}[.|.]$ the
conditional expectation. The later case leads to simple expression:
$\widehat{\psi}_{j}$ is the empirical mean difference between the groups

$\displaystyle\widehat{\psi}_{j}=\frac{1}{n}\sum_{i=1}^{n}\left(\frac{y_{ij}x_{i}}{\pi}-\frac{y_{ij}(1-x_{i})}{1-\pi}\right)$

where $\pi$ denotes the proportion of observations with $X=1$, $y_{ij}$ the
brain signal for individual $i$ processed by pipeline $j$, and $x_{i}$ the
exposure value for individual $i$. The previous expression is equivalent to:

$\displaystyle\sqrt{n}\left(\widehat{\psi}_{j}-\psi_{j}\right)=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\varphi_{\widehat{\psi}_{j}}(\mathcal{O}_{i})$

where
$\varphi_{\widehat{\psi}_{j}}(\mathcal{O}_{i})=\frac{y_{ij}x_{i}}{\pi}-\frac{y_{ij}(1-x_{i})}{1-\pi}-\psi_{j}$,
$\mathcal{O}_{j}=\left(y_{i1},\ldots,y_{iJ},x_{i}\right)$, and $\psi_{j}$ is
the large sample value of $\widehat{\psi}_{j}$. Denoting $Y_{j}$ the random
variable representing the brain signal processed by pipeline $j$, the
estimator can be seen as the empirical average of independent realizations of
a new random variable $\frac{Y_{j}X}{\pi}-\frac{Y_{j}(1-X)}{1-\pi}-\psi_{j}$,
called the influence function of $\widehat{\psi}_{j}$. Thus from the
multivariate central limit theorem we get that the joint distribution of the
estimates is asymptotically multivariate normal. It has mean $\psi$ (under our
assumption that $\forall j\in\\{1,\ldots,J\\},\psi=\psi_{j}$) and its
variance-covariance, denoted $\Sigma_{\widehat{\boldsymbol{\psi}}}$, is the
same as the one of
$\boldsymbol{\varphi}_{\widehat{\boldsymbol{\psi}}}=\left(\varphi_{\widehat{\psi}_{1}},\ldots,\varphi_{\widehat{\psi}_{J}}\right)$
divided by $n$. Because $\boldsymbol{\varphi}_{\widehat{\boldsymbol{\psi}}}$
involves some unknown parameters like $\pi$ and $\psi_{j}$ we do not observe
it and cannot directly estimate $\Sigma_{\widehat{\boldsymbol{\psi}}}$.
However, by plugging our estimates of these unknown parameters we can
approximate $\varphi_{\widehat{\psi}_{j}}$ as
$\widehat{\varphi}_{\widehat{\psi}_{j}}(\mathcal{O}_{ij})=\frac{y_{ij}x_{i}}{\frac{1}{n}\sum_{i=1}^{n}x_{i}}-\frac{y_{ij}(1-x_{i})}{1-\frac{1}{n}\sum_{i=1}^{n}x_{i}}-\frac{1}{n}\sum_{i=1}^{n}\left(\frac{y_{ij}x_{i}}{\pi}-\frac{y_{ij}(1-x_{i})}{1-\pi}\right)$
and approximate $\Sigma_{\widehat{\boldsymbol{\psi}}}$.

In a more general case, we would define a statistical model
$\mathcal{M}(\Theta)$ relating $\boldsymbol{X}$ and $\boldsymbol{Y}$ via a
parameter $\psi$. $\psi$ may be an element of $\Theta$, the set of model
parameters, or a function of elements of $\Theta$. For instance one could use
a latent variable model (LVM) with two latent variables, one summarizing the
brain measurements and another summarizing the exposure variables. $\psi$ is
then the coefficient relating the two latent variables. See Figure 1 of [11]
for a graphical representation of a LVM - in this example the latent variable
“LVu” represents the PET measurement and the latent variable “LVpos” the
memory relative to positive word. In this more general case the previous
decomposition does not hold exactly but up to a residual term $o_{p}(1)$ which
converges to 0 as the sample approaches infinity:

$\displaystyle\sqrt{n}\left(\widehat{\psi}_{j}-\psi_{j}\right)=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\varphi_{\widehat{\psi}_{j}}(\mathcal{O}_{i})+o_{p}(1)$

as indicated in the main text of this article. This decomposition exists for
any estimator $\widehat{\psi}$ derived from an M-estimator [13] (section 5.3),
including likelihood-based estimators. Denote by $\widehat{\theta}_{j}$ the ML
estimator and $\widehat{\Psi}_{j}={c}^{\intercal}\widehat{\theta}_{j}$ the
parameter of interest ($c$ may be a vector starting by 1 and followed by 0’s,
i.e. selects the first element of $\widehat{\theta}_{j}$. The corresponding
influence function only involves the first two derivatives of the log-
likelihood (formula 3.6 in [12]):

$\displaystyle\varphi_{\widehat{\psi}_{j}}(\mathcal{O}_{i})=-{c}^{\intercal}\mathbb{E}\left[\frac{\partial\mathcal{S}_{j}(\mathcal{O}_{i},\theta_{j})}{\partial\theta_{j}}\right]^{-1}\mathcal{S}_{j}(\mathcal{O}_{i},\theta_{j})$

with $\mathcal{S}_{j}(\mathcal{O}_{i},\theta_{j})$ being the score for
individual $i$ when considering pipeline $j$, i.e. vector containing the first
derivatives of the log-likelihood contribution of individual $i$. Once the
influence function has been estimated for each individual, we can use it to
obtain a consistent estimator of $\Sigma_{\widehat{\boldsymbol{\psi}}}$:

$\displaystyle\widehat{\Sigma}_{\widehat{\boldsymbol{\psi}}}=\frac{1}{n}\sum_{i=1}^{n}{\boldsymbol{\varphi}_{\widehat{\boldsymbol{\psi}}}(\mathcal{O}_{i})}^{\intercal}\boldsymbol{\varphi}_{\widehat{\boldsymbol{\psi}}}(\mathcal{O}_{i})$

where
${\boldsymbol{\varphi}_{\widehat{\boldsymbol{\psi}}}(\mathcal{O}_{i})}^{\intercal}$
denotes the transpose of the $J$-dimensional vector of influence functions
relative to individual $i$. Practically speaking the above allows us assess
the variance-covariance matrix of effects across pipelines, e.g., if pipelines
are completely independent in terms of their estimated effect, then
$\widehat{\Sigma}_{\widehat{\boldsymbol{\psi}}}$ would be a diagonal matrix
with the uncertainty of the estimated effect per pipeline in the diagonal.
However, if the estimated effects across pipelines are correlated, the matrix
would not be sparse.

##  Appendix B Integration of the Gaussian density - bivariate case

Figure A: Density of the two dimensional standard normal distribution (colored
surface). The black line delimits the domain
$\mathcal{D}(\widehat{\boldsymbol{t}}_{\max})$ where $J=2$,
$\widehat{t}_{1}=1.5$, $\widehat{t}_{2}=3$, and
$\Sigma_{\widehat{\boldsymbol{\psi}}}$ is the identity matrix with two rows
and two columns. The area under the blue surface external to the black line
corresponds to the p-value relative to the first test, adjusted for two tests.

##  Appendix C Reformulation of the GLS estimator

###  C.1 Via a spectral decomposition

$\widehat{\Sigma}_{\widehat{\psi}}$ being a symmetric semi-positive definite
matrix, it admits the decomposition
$\widehat{\Sigma}_{\widehat{\psi}}=QD{Q}^{\intercal}$ where $Q$ is an
orthogonal matrix (i.e. $Q{Q}^{\intercal}=I_{J}$, the identity matrix) and $D$
is a diagonal matrix with non negative values
$\left(\lambda_{1},\ldots,\lambda_{J}\right)$. Thus:

$\displaystyle\widehat{\boldsymbol{\psi}}_{\text{GLS}}=\left({\mathbf{1}}^{\intercal}\widehat{\Sigma}_{\widehat{\Psi}}^{-1}\mathbf{1}\right)^{-1}{\mathbf{1}}^{\intercal}\widehat{\Sigma}_{\widehat{\Psi}}^{-1}\widehat{\boldsymbol{\psi}}=\left({\mathbf{1}}^{\intercal}QD^{-1}{Q}^{\intercal}\mathbf{1}\right)^{-1}{\mathbf{1}}^{\intercal}QD^{-1}{Q}^{\intercal}\widehat{\boldsymbol{\psi}}$

We first note that $\bar{Q}={\mathbf{1}}^{\intercal}Q$ is a vector with
element $q$ equal to the sum (column-wise) of the eigenvectors. Then:

$\displaystyle{\mathbf{1}}^{\intercal}QD^{-1}{Q}^{\intercal}\mathbf{1}=\sum_{j=1}^{J}\bar{q}_{j}\lambda^{-1}\bar{q}_{j}=\sum_{j=1}^{J}w_{j}$

Where $w_{j}=\left.\bar{q}_{j}\right.^{-2}/\lambda_{j}$. Moreover
${\mathbf{1}}^{\intercal}QD^{-1}{Q}^{\intercal}$ is a vector with $k$-th
elements
$\sum_{j=1}^{J}\bar{q}_{j}\lambda_{j}^{-1}q_{kj}=\sum_{j=1}^{J}w_{j}q^{*}_{kj}$
where $q_{kj}^{*}=q_{kj}/\bar{q}_{j}$

Therefore

$\displaystyle\widehat{\Psi}_{\text{GLS}}=\frac{1}{\sum_{j=1}^{J}w_{j}}\sum_{j=1}^{J}w_{j}{q_{j}^{*}}^{\intercal}\widehat{\boldsymbol{\psi}}$

###  C.2 Via joint modeling

Consider the simple case of a single continuous brain measurement ($R=1$), a
single binary exposure ($P=1$) with equal probability of being $0$ and $1$, no
covariate ($C=0$), and no missing value. We can use a joint linear model:

$\displaystyle Y_{ij}=\alpha_{j}+\beta X_{i}+\varepsilon_{ij}\text{ where
}\left(\varepsilon_{i1},\ldots,\varepsilon_{iJ}\right)\sim\mathcal{N}\left(0,\Sigma_{\varepsilon}\right)$

Denote by
$\widehat{\boldsymbol{\mu}}_{g}=\left(\widehat{\mu}_{g1},\ldots,\widehat{\mu}_{gJ}\right)$
the vector empirical mean in each group (i.e. one relative to $X=1$ another to
$X=0$) of sample size $\frac{n}{2}$. Since the mean and variance are
sufficient statistics in Gaussian models, the joint linear model is equivalent
to:

$\displaystyle\widehat{\mu}_{gj}=\alpha_{j}+\beta X_{g}+e_{gj}\text{ where
}\left(e_{g1},\ldots,e_{gJ}\right)\sim\mathcal{N}\left(0,2\Sigma_{\varepsilon}/n\right)$

Denote by
$\Delta\widehat{\boldsymbol{\mu}}=\left(\Delta\widehat{\mu}_{1},\ldots,\Delta\widehat{\mu}_{J}\right)$
the vector of difference in mean, the previous model implies:

$\displaystyle\Delta\widehat{\mu}_{j}=\beta+\epsilon_{j}\text{ where
}\left(\epsilon_{1},\ldots,\epsilon_{J}\right)\sim\mathcal{N}\left(0,4\Sigma_{\varepsilon}/n\right)$

whose maximum likelihood solution is
$\widehat{\beta}=\left({\mathbf{1}}^{\intercal}\widehat{\Sigma}_{\epsilon}^{-1}\mathbf{1}\right)^{-1}{\mathbf{1}}^{\intercal}\widehat{\Sigma}_{\epsilon}^{-1}\Delta\widehat{\boldsymbol{\mu}}$.

Now consider the GLS estimator of the common exposure effect
$\widehat{\boldsymbol{\psi}}_{\text{GLS}}=\left({\mathbf{1}}^{\intercal}\widehat{\Sigma}_{\widehat{\Psi}}^{-1}\mathbf{1}\right)^{-1}{\mathbf{1}}^{\intercal}\widehat{\Sigma}_{\widehat{\Psi}}^{-1}\widehat{\boldsymbol{\psi}}$
based on the pipeline specific Ordinary Least Squares (OLS) estimators of the
exposure effect
$\widehat{\boldsymbol{\psi}}=\left(\widehat{\psi}_{1},\ldots,\widehat{\psi}_{J}\right)$.
Denote by $\Theta_{j}=\left(\alpha_{j},\psi_{j}\right)$ the mean parameters of
each pipeline specific model and $\sigma^{2}_{j}$ the residual variance
parameter. We have:

  * •

$\widehat{\boldsymbol{\psi}}=\left({X}^{\intercal}X\right)^{-1}{X}^{\intercal}Y_{.j}$
where $Y_{.j}$ denotes the brain measurements across individuals relative to
the $j$-th pipeline.

  * •

$\widehat{\Sigma}_{\widehat{\Psi}}$ has elements
$\sum_{i=1}^{n}{\left(Y_{ij}-X_{i}\Theta_{j}\right)}^{\intercal}\left(Y_{ij}-X_{i}\Theta_{j}\right)\left(X_{i}\left({X}^{\intercal}X\right)^{-1}\right)^{2}{c}^{\intercal}$
where $c=(0,1)$. This follows from the fact the the score relative to
$\widehat{\psi}_{j}$ is
$\frac{1}{\sigma_{j}^{2}}{X}^{\intercal}\left(Y_{j}-X\Theta_{j}\right)c$ and
the variance covariance $\sigma_{j}^{2}\left({X}^{\intercal}X\right)^{-1}$.

With a single binary covariate (and an intercept),
${X}^{\intercal}X=\begin{bmatrix}n&n/2\\\ n/2&n\end{bmatrix}$ whose inverse is
$\begin{bmatrix}2/n&-2/n\\\ -2/n&4/n\end{bmatrix}$ Hence the second element of
$X_{i}({X}^{\intercal}X)^{-1}$ is either $-2/n$ or $2/n$ so
$\left(X_{i}({X}^{\intercal}X)^{-1}\right){c}^{\intercal}=4/n$. Therefore,
assuming that the observations are sorted by group:

$\displaystyle\widehat{\psi}_{j}$
$\displaystyle=\frac{n}{2}\sum_{i=1}^{n/2}Y_{ij}-\frac{n}{2}\sum_{i=n/2+1}^{n}Y_{ij}=\Delta\widehat{\mu}_{j}$
$\displaystyle\widehat{\Sigma}_{\widehat{\psi}}$
$\displaystyle=\frac{4}{n^{2}}\sum_{i=1}^{n}{\left(Y_{i.}-X\Theta\right)}^{\intercal}\left(Y_{i.}-X\Theta\right)=\frac{4}{n}\widehat{\Sigma}_{\varepsilon}$

So
$\widehat{\boldsymbol{\psi}}_{\text{GLS}}=\left({\mathbf{1}}^{\intercal}\widehat{\Sigma}_{\varepsilon}^{-1}\mathbf{1}\right)^{-1}{\mathbf{1}}^{\intercal}\widehat{\Sigma}_{\varepsilon}^{-1}\Delta\widehat{\boldsymbol{\mu}}$.
If the common effect models holds (i.e.
$\widehat{\psi}_{1}=\ldots=\widehat{\psi}_{J}=\beta$) then
$\widehat{\Sigma}_{\varepsilon}$ is an unbiased estimate of
$\Sigma_{\varepsilon}$ so $\widehat{\boldsymbol{\psi}}_{\text{GLS}}$ is
asymptotically equivalent to the ML estimator $\widehat{\beta}$ and thus
asymptotically efficient.


# A Speed Bump: SN 2021aefx Shows that Doppler Shift Alone can Explain Early-
Excess Blue Flux in Some Type Ia Supernovae

C. Ashall Institute for Astronomy, University of Hawaii, 2680 Woodlawn Drive,
Honolulu, HI 96822,USA  J. Lu Department of Physics, Florida State University,
77 Chieftan Way, Tallahassee, FL 32306, USA  B. J. Shappee Institute for
Astronomy, University of Hawaii, 2680 Woodlawn Drive, Honolulu, HI 96822,USA
C. R. Burns The Observatories of the Carnegie Institution for Science, 813
Santa Barbara Street, Pasadena, CA 91101, USA  E. Y. Hsiao Department of
Physics, Florida State University, 77 Chieftan Way, Tallahassee, FL 32306, USA
S. Kumar Department of Physics, Florida State University, 77 Chieftan Way,
Tallahassee, FL 32306, USA  N. Morrell Carnegie Observatories, Las Campanas
Observatory, Casilla 601, La Serena, Chile  M. M. Phillips Carnegie
Observatories, Las Campanas Observatory, Casilla 601, La Serena, Chile  M.
Shahbandeh Department of Physics, Florida State University, 77 Chieftan Way,
Tallahassee, FL 32306, USA  E. Baron Homer L. Dodge Department of Physics and
Astronomy, University of Oklahoma, 440 W. Brooks, Rm 100, Norman, OK
73019-2061, USA  Hamburger Sternwarte, Gojenbergsweg 112, D-21029 Hamburg,
Germany  Department of Physics, The George Washington University, Corcoran
Hall, 725 21st Street NW, Washington, DC 20052, USA  K. Boutsia Carnegie
Observatories, Las Campanas Observatory, Casilla 601, La Serena, Chile  P. J.
Brown George P. and Cynthia Woods Mitchell Institute for Fundamental Physics
and Astronomy, Department of Physics and Astronomy, Texas A&M University,
College Station, TX 77843, USA  J. M. DerKacy Homer L. Dodge Department of
Physics and Astronomy, University of Oklahoma, 440 W. Brooks, Rm 100, Norman,
OK 73019-2061, USA  L. Galbany Institute of Space Sciences (ICE, CSIC), Campus
UAB, Carrer de Can Magrans, s/n, E-08193 Barcelona, Spain  Institut d’Estudis
Espacials de Catalunya (IEEC), E-08034 Barcelona, Spain  P. Hoeflich
Department of Physics, Florida State University, 77 Chieftan Way, Tallahassee,
FL 32306, USA  K. Krisciunas George P. and Cynthia Woods Mitchell Institute
for Fundamental Physics and Astronomy, Department of Physics and Astronomy,
Texas A&M University, College Station, TX 77843, USA  P. Mazzali  Astrophysics
Research Institute, Liverpool John Moores University, IC2, Liverpool Science
Park, 146 Brownlow Hill, Liverpool L3 5RF, UK  A. L. Piro The Observatories of
the Carnegie Institution for Science, 813 Santa Barbara Street, Pasadena, CA
91101, USA  M. D. Stritzinger Department of Physics and Astronomy, Aarhus
University, Ny Munkegade 120, DK-8000 Aarhus C, Denmark  N. B. Suntzeff George
P. and Cynthia Woods Mitchell Institute for Fundamental Physics and Astronomy,
Department of Physics and Astronomy, Texas A&M University, College Station, TX
77843, USA

(Received 11th April 2022 ; Revised 12th May 2022; Accepted 22nd May 2022)

###### Abstract

We present early-time photometric and spectroscopic observations of the Type
Ia Supernova (SN Ia) 2021aefx. The early time $u$-band light curve shows an
excess flux when compared to normal SNe Ia. We suggest that the early-excess
blue flux may be due to a rapid change in spectral velocity in the first few
days post explosion, produced by the emission of the Ca ii H&K feature passing
from the $u$ to the $B$ bands on the time scale of a few days. This effect
could be dominant for all SNe Ia which have broad absorption features and
early-time velocities over 25,000 km s-1. It is likely to be one of the main
causes of early-excess u-band flux in SNe Ia which have early-time high-
velocities. This effect may also be dominant in the UV filters, as well as in
places where the SN spectral energy distribution is quickly rising to longer
wavelengths. The rapid change in velocity can only produce a monotonic change
(in flux-space) in the $u$-band. For objects which explode at lower
velocities, and have a more structured shape in the early-excess emission,
there must also be an additional parameter producing the early-time diversity.
More early time observations, in particular early spectra, are required to
determine how prominent this effect is within SNe Ia.

Supernova

††journal: ApJL

##  1 Introduction

Type Ia Supernovae (SNe Ia) are thought to come from the thermonuclear
explosion of at least one Carbon-Oxygen White Dwarf (WD) in a binary system
(Whelan & Iben, 1973; Iben & Tutukov, 1984). Their peak luminosities are
commonly used to measure luminosity distances to define the Hubble flow (e.g.,
Phillips, 1993; Riess et al., 1998; Perlmutter et al., 1999). Yet, to date,
determining the exact progenitor scenario (e.g. single degenerate vs double
degenerate) or explosion mechanism (Chandrasekhar mass vs sub-Chandrasekhar
mass) has eluded the community (e.g., Maoz et al., 2014; Blondin et al., 2017;
Hoeflich et al., 2017; Ashall et al., 2018; Galbany et al., 2019; Jha et al.,
2019).

In recent years high-cadence transient surveys such as the All-Sky Automated
Survey for Supernovae (ASAS-SN; Shappee et al. 2014; Kochanek et al. 2017),
the Asteroidal Terrestrial-impact Last Alert System (ATLAS; Tonry et al.
2018), the Distance less than 40
(DLT40)111http://dark.physics.ucdavis.edu/dlt40/DLT40, and the Zwicky
Transient Facility (ZTF; Bellm et al. 2019), have dramatically increased the
number of SNe Ia caught within hours to days of explosion. The radiative
transfer through the outer layers which can only be seen up to a few days
after explosion allows us to model the explosion and measure the chemical
properties of the progenitor system.

There have only been a handful of SNe Ia that have multi-band photometry and
were caught within hours of first light. A subset of these objects are thought
to have early-excess blue emission, two-component rising light curves, or blue
bumps in their light curves. Some of these SNe are: SN 2012cg (Marion et al.,
2016; Shappee et al., 2018), SN 2017cbv (Hosseinzadeh et al., 2017),
ASASSN-18bt/2018oh (Shappee et al., 2019; Dimitriadis et al., 2019), SN
2019yvq (Miller et al., 2020; Tucker et al., 2021), and SN 2020hvf (Jiang et
al., 2021). For a summary of SNe Ia with early-excess light-curve emission,
see Jiang et al. (2018). Another subset of SNe Ia have smooth rising light
curves such as SN 2011fe (Nugent et al., 2011), ASASSN-14lp (Shappee et al.,
2016), and SN 2015F (Cartier et al., 2017). Furthermore, the early-time color
curves of SNe Ia have shown diversity (Stritzinger et al., 2018; Bulla et al.,
2020).

The cause of the early-time light curve diversity is unknown, but various
effects have been proposed. These include: 1) the collision of the ejecta with
a companion star in a single degenerate system (Marietta et al., 2000; Kasen,
2010); 2) the presence of 56Ni in the outer layer, which can be mixed out from
the center of the explosion or produced in situ (Piro & Morozova, 2016; Ni et
al., 2022); or 3) the interaction of the ejecta with circumstellar material
(Gerardy et al., 2007; Dragulin & Hoeflich, 2016; Piro & Morozova, 2016; Jiang
et al., 2021).

However, to date, no work has been carried out on determining if rapidly
changing spectral features can produce the observed diversity in the rising
light curves. In this work, we quantify the effect that a rapidly changing
spectral feature has on the fixed-filter early light curves. We use
observations of SN 2021aefx as a test case. In section 2 the data are
presented, in section 3 we discuss how the data are calibrated, in section 4
the results are shown, and the discussion is presented in section 5.

##  2 Data

###  2.1 Discovery

SN 2021aefx was discovered by the DLT40 transient survey on
2021-11-11T12:30:43.78 (MJD 59529.52) at 17.24 mag, with a last non-detection
on 2021-11-06T07:52:17 (MJD 59524.33) at 19.35 mag. Both the discovery and
last non-detection were in a clear filter (Valenti et al., 2021). SN 2021aefx
was classified as a high-velocity SN Ia using a spectrum obtained $\sim$8
hours (2021-11-11T20:55:27; MJD 59529.87) after discovery (Bostroem et al.,
2021). ASAS-SN (Shappee et al. 2014; Kochanek et al. 2017) places tighter
constraints on the last non-detection with a limiting 3-sigma observation of
17.70 mag in the $g$-band on 2021-11-09T08:37:39 (MJD 59527.36). The host of
SN 2021aefx (NGC 1566) has a heliocentric redshift of 0.0050 (Allison et al.,
2014), which corresponds to a distance modulus of 31.43$\pm$0.15 mag when
corrected for the local peculiar motion of the Virgo, Great Attractor, and
Shapley clusters (Mould et al., 2000). Alternatively, using the Tip of the Red
Giant Branch (TRGB) produces a distance modulus of 31.27$\pm$0.49 mag (Sabbi
et al., 2018).

Figure 1: (Top:) Multi-band photometric observations of SN 2021aefx. The ASAS-
SN last non-detection and DLT40 discovery are indicated with open triangle and
filled star, respectively. The dotted lines show the interpolation of the
early time light curves that is used for color correcting the spectra. The
diamonds represent the synthetic magnitudes of integrating the mangled spectra
described in Section 3. The vertical blue lines on the bottom mark the epochs
of the observed spectra of SN 2021aefx. (Middle:) Power-law fits to the rising
$u$-band light. The dash-dotted line is the f$\propto$t2 fit, the dashed line
allows for the power-law index to vary, and the solid line is the two
component fit. (Bottom:) The early-time $B-V$ color curves of SN 2021aefx and
SN 2018aoz. Both SNe have an initial period where they turn red before
becoming blue. The error bars are plotted but are smaller than the points. The
plotted light- and color-curves are corrected for Milky-Way extinction.

###  2.2 Photometry

Multi-band ($uBVgri$) follow-up observations were initiated by the Precision
Observations for Infant Supernovae Explosions (POISE, Burns et al. 2021) a day
after discovery. The photometry was obtained on the 1m Swope telescope at Las
Campanas Observatory. Observations were acquired twice per night in order to
look for small scale fluctuations in the light curve. The photometry was
reduced and analyzed using the same methods outlined by Krisciunas et al.
(2017) and Phillips et al. (2019). The uncertainties shown in Fig. 1 are
derived from the usual propagation of errors assuming Poisson statistics, as
well as the uncertainty in the nightly zero-points derived from observations
of standard stars (see Krisciunas et al. 2017). A further source of potential
error which is not included is the contamination of host galaxy light. While
visual inspection indicates this is likely very small, final photometry and
errors will require observations of the host galaxy after the SN has faded. A
log of the photometric observations can be found in Table 1. The early time
light curve is presented in the top panel of Fig. 1. The middle panel in Fig.
1 demonstrates that the rise is not well characterized by a simple constant-
temperature fireball expansion f $\propto$ t2, which gives an invalid time of
explosion that is after the discovery. However, there is no good reason why
the f $\propto$ t2 should provide an exact match to the data (Nugent et al.,
2011). A power-law that allows the index to vary (f $\propto$ t3.34) fits the
$u$-band rise better, but still underestimates the flux of the first three
data points. In fact, the $u$-band light curve may be best described as having
a two-component rise, which may be identified as early-excess emission in the
two days following discovery. For the two component fit and the fireball model
the 13 data points shown in the middle panel of Fig. 1 used in the fitting
process. The two-component fit has six free parameters and a reduced
$\chi^{2}$ of 28. Whereas the fireball model had two free parameters and a
reduced $\chi^{2}$ of 664, and the one-component power-law fit has three free
parameters and reduced $\chi^{2}$ is 276. We note that the two component rise
in SN 2021aefx is different from the one seen in the Kepler filter for
ASASSN-18bt (Shappee et al., 2019). This is because the Kepler filter covers a
much larger wavelength range ($\sim$4000-9000Å). Therefore, in SN 2021aefx the
effect likely comes from a different physical process.

A full detailed analysis of the later epochs of SN 2021aefx will be presented
in a future paper. However, here we provide some important information about
the SN for context. SN 2021aefx peaked at an apparent $B$-band magnitude of
11.99$\pm$0.01 mag, and has a decline rate of $\Delta{\rm m_{15}(B)}$=1.01
$\pm$ 0.06 mag.222$\Delta{\rm m_{15}(B)}$ is the change in $B$-band magnitude
between maximum light and 15 days past maximum. This corresponds to an
absolute $B$-band magnitude of –19.28$\pm$0.49 mag, using the TRGB distance of
$\mu$=31.27 mag, and a foreground Galactic extinction of $E(B-V)_{MW}$=0.01
mag (Schlafly & Finkbeiner, 2011). Throughout this work the data have been
correct for foreground Galactic extinction. This places it in the normal part
of the luminosity width relation, and in a similar location to SN 2009ig,
which has $M_{B}$= –19.46$\pm$0.50 mag and $\Delta{\rm m_{15}(B)}$ =
0.90$\pm$0.07 mag (Marion et al., 2013).

Furthermore, the $B-V$ color curve of SN 2021aefx has the unusual behavior
where it starts red, and then turns blue within the first two days, see Fig 1.
This is unlike any SN in the $B-V$ sample of Stritzinger et al. (2018), but is
similar to SN 2018aoz. However, the evolution of SN 2018aoz is more extreme
compared to SN 2021aefx. For SN 2018aoz this effect has been interpreted as
enhanced line blanketing in the outer layers of the ejecta caused by iron
group elements (Ni et al., 2022).

Table 1: Swope e2v optical photometry of SN 2021aefx.aaMagnitudes in the Swope+e2v natural system. MJD | PhasebbRelative to the $B$-band maximum of 2021aefx at MJD=59547.29. | $u$ | $g$ | $B$ | $V$ | $r$ | $i$
---|---|---|---|---|---|---|---
59530.2 |  $-$16.9 | 17.39(0.02) | 16.48(0.01) | 16.24(0.01) | 16.04(0.01) | 16.05(0.01) | 16.57(0.01)
59530.4 |  $-$16.7 | 17.34(0.04) | 16.43(0.01) | 16.15(0.01) | 15.93(0.01) | 15.94(0.01) | 16.47(0.01)
59531.2 |  $-$15.9 | 17.26(0.02) | 16.06(0.01) | 15.71(0.01) | 15.38(0.01) | 15.41(0.01) | 16.04(0.01)
59532.1 |  $-$15.0 | 16.77(0.02) | 15.46(0.01) | 15.15(0.01) | 14.86(0.01) | 14.91(0.01) | 15.54(0.01)
59532.3 |  $-$14.8 | 16.60(0.02) | 15.28(0.01) | 14.98(0.01) | 14.74(0.01) | 14.77(0.01) | 15.40(0.02)
59533.1 |  $-$14.0 | 16.03(0.02) | 14.74(0.01) | 14.53(0.01) | 14.34(0.01) | 14.37(0.01) | 14.99(0.01)
59533.3 |  $-$13.8 | 15.88(0.02) | 14.56(0.01) | 14.39(0.01) | 14.24(0.01) | $\cdots$ | 14.90(0.01)
59534.1 |  $-$13.0 | 15.34(0.02) | 14.09(0.01) | 13.98(0.01) | 13.90(0.01) | 13.95(0.01) | 14.53(0.01)
59534.3 |  $-$12.8 | 15.16(0.02) | 13.94(0.01) | $\cdots$ | 13.8(0.01) | 13.85(0.01) | 14.41(0.01)
59535.1 |  $-$12.0 | 14.70(0.02) | 13.61(0.01) | 13.56(0.01) | 13.54(0.01) | 13.57(0.01) | 14.12(0.01)
59535.3 |  $-$11.8 | $\cdots$ | 13.50(0.01) | 13.46(0.01) | 13.48(0.01) | 13.49(0.01) | 14.03(0.01)
59536.1 |  $-$11.0 | 14.18(0.02) | 13.24(0.01) | 13.20(0.01) | 13.22(0.01) | 13.25(0.01) | 13.75(0.01)
59536.3 |  $-$10.8 | 14.06(0.02) | 13.14(0.01) | 13.12(0.01) | 13.18(0.01) | 13.18(0.01) | 13.66(0.01)
59537 |  $-$10.1 | $\cdots$ | 12.96(0.01) | $\cdots$ | 12.97(0.01) | 12.97(0.01) | 13.42(0.01)
59537.1 |  $-$10.0 | 13.80(0.04) | $\cdots$ | 12.90(0.01) | $\cdots$ | $\cdots$ | $\cdots$
59537.3 |  $-$9.9 | 13.66(0.02) | 12.87(0.01) | 12.89(0.01) | 12.92(0.01) | $\cdots$ | 13.40(0.01)
59538.1 |  $-$9.1 | 13.40(0.02) | 12.71(0.01) | 12.69(0.01) | 12.75(0.01) | 12.77(0.01) | 13.19(0.01)
59538.3 |  $-$8.9 | 13.36(0.02) | 12.66(0.01) | 12.66(0.01) | 12.72(0.01) | 12.74(0.01) | 13.15(0.01)
59539.2 |  $-$8.0 | 13.17(0.02) | 12.52(0.01) | 12.52(0.01) | 12.60(0.01) | $\cdots$ | 13.04(0.01)
59539.3 |  $-$7.9 | $\cdots$ | 12.50(0.01) | $\cdots$ | 12.58(0.01) | 12.58(0.01) | $\cdots$
59539.4 |  $-$7.8 | 13.12(0.02) | $\cdots$ | 12.50(0.01) | $\cdots$ | $\cdots$ | 13.01(0.01)
59540.1 |  $-$7.1 | 12.98(0.01) | 12.40(0.01) | 12.38(0.01) | 12.46(0.01) | 12.49(0.01) | 12.90(0.01)
59540.3 |  $-$6.9 | 12.92(0.02) | 12.38(0.01) | 12.35(0.01) | 12.44(0.01) | 12.46(0.01) | 12.86(0.01)
59541.1 |  $-$6.1 | 12.82(0.05) | 12.27(0.01) | 12.26(0.01) | 12.36(0.01) | 12.40(0.01) | 12.81(0.01)

###  2.3 Spectroscopy

We present seven spectra of SN 2021aefx. Four were taken from the Transient
Name Server database,333https://www.wis-tns.org/ one was acquired using MIKE
on the 6.5 m Magellan Clay Telescope, and two were obtained using IMACS on the
6.5 m Magellan Baade Telescope at Las Campanas Observatory. The spectra were
reduced using standard iraf444The Image Reduction and Analysis Facility (iraf)
is distributed by the National Optical Astronomy Observatory, which is
operated by the Association of Universities for Research in Astronomy, Inc.,
under cooperative agreement with the National Science Foundation. packages
following the method as described by Hamuy et al. (2006). A log of the spectra
can be found in Table 2. The spectra are presented in Fig. 2.

Table 2: Journal of Spectroscopic Observations. UT Date | MJD | PhaseaaRelative to the $B$-band maximum of 2021aefx at MJD=59547.29. | Telescope/Instrument | Source
---|---|---|---|---
|  | [day] |  |
2021-11-11 | 59529.87 | -17.3 | SALT + RSS | TNS
2021-11-13 | 59531.27 | -15.9 | NTT + EFOSC2 | ePESSTO+
2021-11-14 | 59532.08 | -15.1 | Alpy600 + 35cm R-C | TNS
2021-11-21 | 59539.29 | -7.9 | Clay + MIKE | POISE
2021-12-03 | 59551.21 | 3.9 | Alpy600 + 35cm R-C | TNS
2021-12-09 | 59557.35 | 10.0 | Baade + IMACS | POISE
2021-12-16 | 59564.27 | 16.9 | Baade + IMACS | POISE
Figure 2: Spectra comparison of SN 2021aefx and a variety of other SNe Ia. SN
2021aefx is most similar to SNe 2009ig and 2012fr. Interestingly, the
comparison also shows that 2021aefx is not the same as SN 2012cg or SN 2017cbv
despite them having an early-excess blue flux. Throughout this work, we use
the spectra of SN 2009ig to fill in gaps in wavelength and phase.  Figure 3:
(Left:) The six color-adjusted spectra chosen for this analysis, including
three early spectra of SN 2021aefx and three spectra of SN 2009ig. The spectra
are color-matched to the observed photometry of SN 2021aefx. (Top right:) An
example of the merging process between SN 2021aefx and SN 2009ig. Before the
merging process, the spectra of SN 2009ig are shifted in velocity to match the
velocity of SN 2021aefx at the same phase. The gray shade marks the
overlapping region used to match the flux scale. (Bottom right:) Velocity
profiles of the minimum of the Si ii $\lambda$6355 feature in SN 2011fe, SN
2009ig, and SN 2021aefx.

The earliest spectrum of SN 2021aefx (at –17.3 d555Throughout this work, all
phases (in days) are provided relative to rest frame $B$-band maximum.) was
obtained during the early-excess emission phase in the $u$-band light curve.
This early-time spectrum shows high-velocity and broad spectral features.
Whereas, by –7.9 d the spectrum resembles a normal SN Ia. At this phase it is
located between the Core normal and Broad Line SNe in the Branch diagram
(Branch et al., 2009). At –17.3 d, the velocity of the absorption minimum of
the Si ii $\lambda$6355 feature reaches $\sim$30,000 km s-1(see bottom right
panel of Fig. 3). Within 1.4 d, this velocity drops to $\sim$21,000 km s-1.
This deceleration is larger than other SNe Ia caught in these early stages,
for example SN 2011fe declined from 16,000 km s-1 to 13,000 km s-1 over
similar epochs. In fact, at early times all of the spectral lines in SN
2021aefx are higher in velocity compared to SN 2011fe, as is seen in Fig. 2.

Although it is apparent that there is diversity in the rising light curves of
SNe Ia, few spectra have been obtained or analyzed during these early phases.
SN 2021aefx is the first public example of a SN with an early excess $u$-band
emission which also has high velocity spectra in the early phases. Given that
the velocity of the spectral lines rapidly decreases over the course of a few
days, it is natural to ask the question: Could a rapid change in velocity of
the ejecta produce the early excess u-band emission? For example, it could be
the case that at early times the spectral features are so shifted to the blue
that they affect the broad-band photometry. The blue shift may add flux to the
$u$-band. We investigate this further below. In fact, the sensitivity of
broad-band photometry to spectral features near the edges of the photometric
bandpasses was noted as far back as Suntzeff et al. (1988) when comparing CTIO
and SAAO $I$-band photometry. which could differ by 0.3 mag.

##  3 Constructing spectrophotometric data

While the spectra cover most of the $u$ band, they do not cover the entire
bandpass so careful consideration is warranted. Therefore, assumptions must be
made about the nature of the SED in this region. In order to establish which
objects have the same ionization state and line ratios, SN 2021aefx is
compared to a range of SNe Ia which have early-time spectra. This is done to
ensure that any extrapolation to the blue is made with a SN which has similar
ejecta properties. At early times, 2021aefx does not resemble SN 2011fe, SN
2012cg, or SN 2017cbv, but is remarkably similar to SN 2012fr (Contreras et
al., 2018) and SN 2009ig (Foley et al., 2012; Marion et al., 2013). Figure 2
shows these comparisons. The latter objects are SNe with high-velocities and
broad features at early times. For example, at $\sim$–15 d the pseudo
equivalent width (pEW) of the Si ii 6355 feature for SN 2021aefx, SN 2009ig
and SN 2012fr is larger than 150 $\AA$, whereas the pEW of SN 2011fe is
$\sim$100 $\AA$. Conveniently, some early-time spectra of SN 2009ig extend
into the blue and cover the whole wavelength range of the $u$-band.

In order to extend the spectra of SN 2021aefx to cover the whole of the
$u$-band region, we merge the blue part of the spectra of SN 2009ig onto the
spectra of SN 2021aefx. However, before doing this, the spectra of SN 2009ig
are matched to the closest phase of SN 2021aefx, and then shifted in velocity
to match the optical features of SN 2021aefx. The assumption made here is that
if SNe Ia are very similar in the $BVgri$ bands, they should also be similar
in the $u$-band. We do note that SNe which are optical twins might not always
have similar mid and far-UV flux. For example, SNe 2011fe and 2011by were
spectroscopic twins in the optical but differed in the mid and far UV (Foley &
Kirshner, 2013; Graham et al., 2015). However, when the spectra are normalized
in the 4000-5500Å region at –10 d SNe 2011fe and 2011by only differ in
$u$-band flux by 6%.

An example of the merged spectrum is presented in the top right panel of
Figure 3. SN 2009ig does not have $u$-band photometry early enough to
determine if it has an early-excess blue flux. However, at the earliest time
96% of the $u$-band flux in the merged spectrum region comes from the data of
SN 2021aefx, and only 4% comes from the added/merged data. To determine if the
choice of SN for the merging process affects the results, we also ran the
trial using SN 2011fe as the merged spectrum. Using SN 2011fe only 9% of the
flux in the $u$-band comes from SN 2011fe. Therefore, the merging process does
not significantly affect our results.

The merged spectra are then color matched (also referred to as mangled) with
the $uBVgri$ bands to the observed photometry of SN 2021aefx (see the left
panel of Fig. 3). The color matching is done using the mangle function in
SuperNovae in object oriented Python (SNooPy) program (Burns et al., 2011).
Cubic splines were fit with control points located at the effective
wavelengths of each filter. In addition, two "stabilizing" filters are
introduced 100 Å to the blue of the bluest filter and to the red of the
reddest filter. These extra control points are given the same flux as through
the bluest and reddest filters, forcing the spline to have zero slope at the
ends. Integrating the spectra through the $uBVgri$ filters at each epoch
correctly produces the observed light curve, as is seen in Fig. 1.

##  4 The effects of rapidly changing velocity profiles

Having established that the dataset used in the analysis is
spectrophotometric, we can then begin to explore if different velocity
evolutions of the ejecta would produce changes in the early time behavior of
the light curve. The first approach is to study the effect of velocity on the
$u$\- and $B$-band photometry at a single epoch. In the top left panel of Fig.
4, the –17.3 d spectra of SN 2021aefx are presented at various velocities,
shifted by increments of 1000 km s-1. Note that where we discuss the velocity
shift, it is relative to the minimum of the Si ii $\lambda$6355 feature. At
each velocity-shift the synthetic $u$\- and $B$-band magnitudes are
calculated.666It is important to note, that we are moving the spectra without
moving the band-pass. In this case the zero point does not need recalculating.
However, if the band-pass was moved the zero point would need recalculating.
The assumption made here is that the underlying chemical distribution and
state of the spectra stay the same at every velocity (i.e. the ionization
state and line ratios do not change), and only a wavelength shift is applied.
The inlay in the figure demonstrates how the synthetic magnitude changes as a
function of the velocity shift. Across a velocity range of 14,000 km s-1,
starting at the observations (30,000 km s-1) and decreasing to 16,000 km s-1,
the synthetic $B$-band stays constant. However, the synthetic $u$-band drops
by $\sim$1 mag. This demonstrates that, at a given time, spectra with higher
velocities are likely to have brighter $u$-band magnitudes. The drop in flux
is caused by the steepness of the spectrum across the $u$ passband,
specifically the emission component from the Ca ii H&K feature is being
redshifted from the $u$\- to the $B$-band as the velocity of the spectrum
decreases. In the middle left panel of Fig. 4 we also show the same velocity
test carried out with the $-14.2$ d spectrum of SN 2009ig. The difference in
in the $u$-band magnitude between the low velocity and high velocity spectrum
is 0.87 mag.

The next stage of the analysis is to determine if various velocity evolution
pathways of the spectral time series significantly affect the rising flux of
the SN. To do this, a range of velocity profiles was created between the
evolution of SN 2021aefx and SN 2011fe; see the right panel of Fig. 4. There
is no significant difference between the low-velocity and high-velocity
profiles in the $B$ band, however there is a clear difference in the early
$u$-band behavior. The shape of the rising light curve in $u$ changes
dramatically within the first few days between the high-velocity and low-
velocity time series. The high-velocity time series produces an apparent
early-excess flux within the $u$-band. If the different synthetic $u$-bands
are fit with a power law, the high velocity time series requires a two-
component power law fit, whereas the low-velocity time series can be fit by a
power law of f $\propto$ t7.7.

In the $u$-band the low-velocity time series is 0.96 mag fainter than the
high-velocity time series (i.e. the observations) at –17.3 d. The low-velocity
time series rises smoothly and is similar to that of a ‘normal’ SN Ia light
curve. However, the high-velocity time series has an early-excess flux in the
$u$ band. These results suggest that in this case, this excess flux is caused
by line profiles shifting to shorter wavelengths. In particular, at early
times the Ca ii H&K emission is located within the $u$ band for SNe that have
high velocities, but not for SNe Ia that have lower velocities. An example of
how different velocity gradients affects the underlying opacity can be seen in
Fig. 1 & 2 of Hoeflich et al. (1993).

Figure 4: (Top left:) Illustration of how changing the velocity of the –17.3
day spectrum can affect the $u$-band and $B$-band synthetic magnitudes. The
scale on the right denotes the velocity shift, where 30,000 km s-1, is the
velocity of the observed spectra. (Middle left:) Similar to the top left panel
but with the spectrum of SN 2009ig at -14.2 d (Foley et al., 2012). In this
panel we also present the magnitude change in the swift $U$-band. (Bottom
left:) Similar to the top left panel but with the spectrum of SN 2011fe at
-13.3 d (Mazzali et al., 2014) taken with the Hubble Space Telescope.
Interestingly, the drop in magnitude effect is larger in the $UVM2$ band than
the $u$-band, on top of this the drop in $UVW1$ is similar in magnitude to
that of the $u$-band. (Right:) The synthetic light curves (bottom panel) of
the spectrophotometric time series with different velocity evolution. The
upper panel shows the chosen velocity evolution for this analysis. These are
velocity profiles between SN 2021aefx and SN 2011fe. The $u$-band light curves
show significant changes in magnitudes and shape for various velocity
profiles, whereas the $B$-band remains relatively constant. At –17.3 d the $u$
band changes by 0.96 mag between the high-velocity and low-velocity models.
Conversely, the $B$-band only changes by 0.03 mag at the same epoch.

Many SNe Ia that have early-excess emission in the $u$-band also have excess-
emission at mid-UV wavelengths. SN 2021aefx does not have early time UV
observations available, therefore we cannot test the effects directly for this
object. However, as a test case we use the -13.3 d spectrum of SN 2011fe to
determine if the rapid drop in spectral flux can also produce early-excess
flux in the UV (see the bottom left panel of Fig. 4). For SN 2011fe the Swift
$UVM2$ and $UVW1$ magnitudes drops by 0.74 mag and 0.41 mag respectively over
a velocity change from 30,000 km s-1 to 13,000 km s-1, but stays constant in
the $B$-band. We note that the HST spectra does not go far enough into the
blue to test that $UVw2$ filter. While the spectrum of SN 2011fe is taken at a
different time compared to the phase where the blue early-excess is seen, it
acts as a proof of concept. Generally, if there is a large slope in the flux
across a passband then a rapid change in velocity of the profiles can cause
additional structure within the rising light curves. However, we note that
there is another effect of velocities at these mid-UV wavelength ranges, since
the opacity in this region is almost entirely dominated by line blanketing due
to iron group elements (DerKacy et al., 2020). Increasing velocity tends to
depress the flux in the UV since the increased velocity leads to an increased
effective opacity in this region (Wang et al., 2012), which our assumptions
neglects. This effect is clearly seen in SNe 2017erp and 2021fxy (J.M. DerKacy
et al., in preparation).

Finally, we tested how the $u$-band flux changes with spectra of SNe Ia which
do not have broad spectral features (see the bottom left panel of Fig. 4). By
using the spectrum of SN 2011fe at –13.3 d the minimum of the features were
artificially increased in velocity to 30,000 km s-1. The change in $u$-band
magnitude between the high-velocity spectrum and low-velocity spectrum is 0.5
mag demonstrates that velocity has an important effect. However, the effect of
the velocity shift is smaller in SN 2011fe than in SN 2021aefx because the
features in SN 2011fe are narrower. For SN 2021aefx with its broader features,
the short wavelength edge of the absorption of Ca ii H&K will eat into the
emission region of the adjacent shorter wavelength features, producing excess
line blanketing and further reducing the flux in the region. Therefore, we
conclude the effect seen in SN 2021aefx is produced by both the high-velocity
and the broadness of the spectral features.

After the submission of this work Hosseinzadeh et al. (2022) presented early
time observations of SN 2021aefx. They showed that early-time open-filter
light curve of SN 2021aefx had structure which may be interpreted as an early-
excess flux. However, the early-time shoulder in the DLT40 open filter is
different from the early-excess blue flux see in the $u$ and $UV$ filters, but
is similar in shape our $V$-band light curves. In this work we concentrate on
the nature of the $u$ and $UV$ early flux as this shows the largest early-
excess flux. Our work here demonstrates that for this particular SN, the
Doppler effect accounts for most of the excess in the $UV$ and $u$-bands
leaving little room for additional effects. Therefore, any model fitting of
early SNe Ia must consider the effect of Doppler shift on the $u$ and $UV$
light curves.

##  5 Discussion

We provide a cautionary tale that the early time diversity seen in $u$-band
light curves of SNe Ia can be caused, in part, by the rapid change in the
velocity of the Ca II H+K feature. This prediction is only true for SNe Ia
which have early-time high-velocity and broad spectral profiles (i.e. objects
which have expansion velocities higher than 25,000km s-1 in the earliest
phases).

The effect is produced by the Ca ii H&K emission rapidly passing from the $u$
band and into the $B$ band causing a change in flux and change in slope of the
$u$-band light curve. The rapid change in velocity can produce various
$u$-band light curves which monotonically increase, but it cannot explain all
of the early time diversity in SNe Ia light curves, such as the more complex
bump-like shapes seen in objects like SN 2017cbv (Hosseinzadeh et al., 2017).
In fact, given the fact that SN 2017cbv has ejecta velocities lower than
$\eqsim$25,000 km s-1 we would not expect the effect discussed here to be the
dominant cause. Furthermore, in normal velocity SNe Ia such as SN 2011fe this
effect is also not expected to be large. However, for a light curve of an
object such as SN 2012cg (Marion et al., 2016) the effect discussed here may
be able to explain the early blue excess flux. Only spectra taken in this
earliest phase would reveal if the lines are of high enough velocity or broad
enough for the effect to be dominant.

The work here should also provide a word of warning that if one is to
interpret the $u$-band flux between observations at different telescopes of
the same SN, a precise knowledge of the transmission function is required
otherwise variation between $u$-band flux may be due to whether the blue edge
of the Ca ii H&K feature is in the bandpass or not.

To fully understand the effects that a rapid change in velocity has on the
light curve requires high-cadence spectral data which cover the whole of the
rest-frame $u$ band, and are obtained within hours after the explosion. Future
follow-up observations of transients should ensure that the whole wavelength
range is observed.

Finally, although this work may not explain all of the diversity within early
SNe Ia light curves, our results should be carefully considered when
interpreting physical meaning from SNe Ia which have early-time high ejecta
velocities and early-excess UV and $u$-band light curves.

##  6 Acknowledgments

We thank the referee for their helpful and constructive comments. C.A. and
B.J.S. are supported by NSF grants AST-1907570, AST-1908952, AST-1920392, and
AST-1911074. M.D.S. of funded in part by an Experiment grant (# 28021) from
the Villum FONDEN, and by a project 1 grant (#8021-00170B) from the
Independent Research Fund Denmark (IRFD) PH acknowledge support by National
Science Foundation (NSF) grant AST- 1715133\. E.B. and J.D. are supported in
part by NASA grant 80NSSC20K0538. This work has been generously supported by
the National Science Foundation under grants AST-1008343, AST-1613426,
AST-1613455, and AST1613472. This paper includes data gathered with the 6.5
meter Magellan Telescopes located at Las Campanas Observatory, Chile. We would
like to thank the technical staff for constant support for observations on the
Swope telescope. The early time spectrum which was critical for this analysis
came from SALT through Rutgers University time via program 2021-1-MLT-007 (PI:
Jha). This paper includes data gathered with the 6.5 meter Magellan Telescopes
located at Las Campanas Observatory, Chile L.G. acknowledges financial support
from the Spanish Ministerio de Ciencia e Innovación (MCIN), the Agencia
Estatal de Investigación (AEI) 10.13039/501100011033, and the European Social
Fund (ESF) "Investing in your future" under the 2019 Ramón y Cajal program
RYC2019-027683-I and the PID2020-115253GA-I00 HOSTFLOWS project, from Centro
Superior de Investigaciones Científicas (CSIC) under the PIE project
20215AT016, and the program Unidad de Excelencia María de Maeztu
CEX2020-001058-M.

## References

  * Allison et al. (2014) Allison, J. R., Sadler, E. M., & Meekin, A. M. 2014, MNRAS, 440, 696
  * Ashall et al. (2018) Ashall, C., Mazzali, P. A., Stritzinger, M. D., et al. 2018, MNRAS, 477, 153
  * Bellm et al. (2019) Bellm, E. C., Kulkarni, S. R., Graham, M. J., et al. 2019, PASP, 131, 018002
  * Blondin et al. (2017) Blondin, S., Dessart, L., Hillier, D. J., & Khokhlov, A. M. 2017, MNRAS, 470, 157
  * Bostroem et al. (2021) Bostroem, K. A., Jha, S. W., Randriamampandry, S., et al. 2021, Transient Name Server Classification Report, 2021-3888, 1
  * Branch et al. (2009) Branch, D., Chau Dang, L., & Baron, E. 2009, PASP, 121, 238
  * Bulla et al. (2020) Bulla, M., Miller, A. A., Yao, Y., et al. 2020, ApJ, 902, 48
  * Burns et al. (2021) Burns, C., Hsiao, E., Suntzeff, N., et al. 2021, The Astronomer’s Telegram, 14441, 1
  * Burns et al. (2011) Burns, C. R., Stritzinger, M., Phillips, M. M., et al. 2011, AJ, 141, 19
  * Cartier et al. (2017) Cartier, R., Sullivan, M., Firth, R. E., et al. 2017, MNRAS, 464, 4476
  * Contreras et al. (2018) Contreras, C., Phillips, M. M., Burns, C. R., et al. 2018, ApJ, 859, 24
  * DerKacy et al. (2020) DerKacy, J. M., Baron, E., Branch, D., et al. 2020, ApJ, 901, 86
  * Dimitriadis et al. (2019) Dimitriadis, G., Foley, R. J., Rest, A., et al. 2019, ApJ, 870, L1
  * Dragulin & Hoeflich (2016) Dragulin, P., & Hoeflich, P. 2016, ApJ, 818, 26
  * Foley & Kirshner (2013) Foley, R. J., & Kirshner, R. P. 2013, ApJ, 769, L1
  * Foley et al. (2012) Foley, R. J., Challis, P. J., Filippenko, A. V., et al. 2012, ApJ, 744, 38
  * Galbany et al. (2019) Galbany, L., Ashall, C., Höflich, P., et al. 2019, A&A, 630, A76
  * Gerardy et al. (2007) Gerardy, C. L., Meikle, W. P. S., Kotak, R., et al. 2007, ApJ, 661, 995
  * Graham et al. (2015) Graham, M. L., Foley, R. J., Zheng, W., et al. 2015, MNRAS, 446, 2073
  * Hamuy et al. (2006) Hamuy, M., Folatelli, G., Morrell, N. I., et al. 2006, PASP, 118, 2
  * Hoeflich et al. (1993) Hoeflich, P., Mueller, E., & Khokhlov, A. 1993, A&A, 268, 570
  * Hoeflich et al. (2017) Hoeflich, P., Hsiao, E. Y., Ashall, C., et al. 2017, ApJ, 846, 58
  * Hosseinzadeh et al. (2017) Hosseinzadeh, G., Sand, D. J., Valenti, S., et al. 2017, ApJ, 845, L11
  * Hosseinzadeh et al. (2022) Hosseinzadeh, G., Sand, D. J., Lundqvist, P., et al. 2022, arXiv e-prints, arXiv:2205.02236
  * Iben & Tutukov (1984) Iben, I., J., & Tutukov, A. V. 1984, ApJS, 54, 335
  * Jha et al. (2019) Jha, S. W., Maguire, K., & Sullivan, M. 2019, Nature Astronomy, 3, 706
  * Jiang et al. (2018) Jiang, J.-a., Doi, M., Maeda, K., & Shigeyama, T. 2018, ApJ, 865, 149
  * Jiang et al. (2021) Jiang, J.-a., Maeda, K., Kawabata, M., et al. 2021, ApJ, 923, L8
  * Kasen (2010) Kasen, D. 2010, ApJ, 708, 1025
  * Kochanek et al. (2017) Kochanek, C. S., Shappee, B. J., Stanek, K. Z., et al. 2017, PASP, 129, 104502
  * Krisciunas et al. (2017) Krisciunas, K., Contreras, C., Burns, C. R., et al. 2017, AJ, 154, 211
  * Maoz et al. (2014) Maoz, D., Mannucci, F., & Nelemans, G. 2014, ARA&A, 52, 107
  * Marietta et al. (2000) Marietta, E., Burrows, A., & Fryxell, B. 2000, ApJS, 128, 615
  * Marion et al. (2013) Marion, G. H., Vinko, J., Wheeler, J. C., et al. 2013, ApJ, 777, 40
  * Marion et al. (2016) Marion, G. H., Brown, P. J., Vinkó, J., et al. 2016, ApJ, 820, 92
  * Mazzali et al. (2014) Mazzali, P. A., Sullivan, M., Hachinger, S., et al. 2014, MNRAS, 439, 1959
  * Miller et al. (2020) Miller, A. A., Magee, M. R., Polin, A., et al. 2020, ApJ, 898, 56
  * Mould et al. (2000) Mould, J. R., Huchra, J. P., Freedman, W. L., et al. 2000, ApJ, 529, 786
  * Ni et al. (2022) Ni, Y. Q., Moon, D.-S., Drout, M. R., et al. 2022, Nature Astronomy, arXiv:2202.08889
  * Nugent et al. (2011) Nugent, P. E., Sullivan, M., Cenko, S. B., et al. 2011, Nature, 480, 344
  * Perlmutter et al. (1999) Perlmutter, S., Aldering, G., Goldhaber, G., et al. 1999, ApJ, 517, 565
  * Phillips (1993) Phillips, M. M. 1993, ApJ, 413, L105
  * Phillips et al. (2019) Phillips, M. M., Contreras, C., Hsiao, E. Y., et al. 2019, PASP, 131, 014001
  * Piro & Morozova (2016) Piro, A. L., & Morozova, V. S. 2016, ApJ, 826, 96
  * Riess et al. (1998) Riess, A. G., Filippenko, A. V., Challis, P., et al. 1998, AJ, 116, 1009
  * Sabbi et al. (2018) Sabbi, E., Calzetti, D., Ubeda, L., et al. 2018, ApJS, 235, 23
  * Schlafly & Finkbeiner (2011) Schlafly, E. F., & Finkbeiner, D. P. 2011, ApJ, 737, 103
  * Shappee et al. (2018) Shappee, B. J., Piro, A. L., Stanek, K. Z., et al. 2018, ApJ, 855, 6
  * Shappee et al. (2014) Shappee, B. J., Prieto, J. L., Grupe, D., et al. 2014, ApJ, 788, 48
  * Shappee et al. (2016) Shappee, B. J., Piro, A. L., Holoien, T. W. S., et al. 2016, ApJ, 826, 144
  * Shappee et al. (2019) Shappee, B. J., Holoien, T. W. S., Drout, M. R., et al. 2019, ApJ, 870, 13
  * Stritzinger et al. (2018) Stritzinger, M. D., Shappee, B. J., Piro, A. L., et al. 2018, ApJ, 864, L35
  * Suntzeff et al. (1988) Suntzeff, N. B., Hamuy, M., Martin, G., Gomez, A., & Gonzalez, R. 1988, AJ, 96, 1864
  * Tonry et al. (2018) Tonry, J. L., Denneau, L., Heinze, A. N., et al. 2018, PASP, 130, 064505
  * Tucker et al. (2021) Tucker, M. A., Ashall, C., Shappee, B. J., et al. 2021, ApJ, 914, 50
  * Valenti et al. (2021) Valenti, S., Sand, D. J., Wyatt, S., et al. 2021, Transient Name Server Discovery Report, 2021-3864, 1
  * Wang et al. (2012) Wang, X., Wang, L., Filippenko, A. V., et al. 2012, ApJ, 749, 126
  * Whelan & Iben (1973) Whelan, J., & Iben, Icko, J. 1973, ApJ, 186, 1007