\section{Introduction}

Recent literature on astrophysical shocks concludes with increasing conviction that a shock, such as the solar-wind termination shock, is a phenomenon in which multifluid effects play an important role. This means that, for the understanding of  global shock physics,  effects arising from different plasma species must be taken into account while all species have to fulfill the Rankine-Hugoniot conservation laws concertedly. For instance, previous studies indicate that pickup ions substantially modify the structure of the solar-wind termination shock.
In addition,  pickup ions are generally considered an important fluid component since they transport a major fraction of the entropized upstream kinetic energy in the form of thermal energy into the downstream regime of the shock \citep[see][]{decker08}. \citet{zank10} and
\citet{fahr07,fahr10,fahr11} found relations between the upstream and downstream ion distribution functions that are different for solar-wind protons and for pickup protons, respectively.
However, these refinements did not lead to a fully satisfying representation of some properties of the shocked plasma observed by Voyager 2 \citep{richardson08}.
In fact, \citet{chalov13} demonstrate that it is necessary to allow for solar-wind electrons to behave like a third, independent plasma fluid  to achieve an agreement between the model results and the plasma data presented by \citet{richardson08}.
Their parameter study shows that the best-fit results for the observed proton temperature are then achieved if the electrons are assumed to be heated to a temperature that is higher than the downstream proton temperature by a factor $\gtrsim 10$.
We refer to the works of \citet{leroy84}, \citet{tokar86}, and \citet{schwartz88} on electron heating at fast-mode shocks.
Also, plasma shock simulations, in which electrons are treated kinetically, show preferential shock heating of the plasma electrons \citep[see][]{lembege03,lembege04}. Two-stream and viscous interactions lead to a demagnetization of the electrons; i.e., they are not Lorentz-wound by the magnetic fields and attain an increase in their temperature by a factor of 50 or more on the downstream side of the shock.
In a different approach, \citet{leroy82} and \citet{goodrich84} predict that electrons carry out perpendicular electric drifts that are different from those of the ions due to the shock-electric field, thereby establishing an electric current that modifies the surface-parallel magnetic field.

With a similar level of consistency,  \citet{fahr12} describe
the conditions of the upstream and downstream plasma in bulk-frame
systems with frozen-in magnetic fields. They treat the transition from the upstream side to the downstream side of the shock as an instantaneous kinetic reaction in the velocity distribution function via the Liouville-Vlasov theorem, predicting all of the relevant downstream plasma quantities. In this model, a semikinetic representation of the multifluid termination shock with mass- and charge-specific reactions of protons and electrons to the electric shock ramp leads to excessive electron heating as described by \citet{fahr12} or \citet{fahr13}.
Based on these studies, we assume that electrons enter the downstream side as a strongly heated, massless plasma fluid that dominates the downstream plasma pressure.
A recently published multifluid termination shock reconstruction by \citet{zieger15}, which describes the electrons as a separate, independently reacting fluid, supports this notion.
This study finds that electrons are preferentially heated with respect to ions, converting overshoot velocities into thermal energy supporting the claims made by \citet{chalov13}, \citet{fahr15}, and \citet{fahr15b}. Furthermore, the study finds that the partial pressure of the heated electrons plays a dominant role in the downstream plasma flow.

In the following study, we investigate more carefully for what reason and in which way the plasma electrons  behave as a preferentially heated fluid at the shock crossing.

\section{The shock electric field and charge separation}

Strongly space-dependent electric and magnetic fields are typical for the transition region of astrophysical shocks, and none of the conventional test-particle descriptions, such as the expansion of the particle motion into different forms of electromagnetic drifts, apply under these conditions.
This is certainly the case when electric forces locally flash up and become strongly dominant over Lorentz forces so that the latter forces can be safely neglected for the first-order motion analysis. In this section, we look into the details of this special situation that arises at shocks such as the solar wind-termination shock.

We assume that the shock is characterized by a three-phase transition in space. In the first phase (the demagnetization region),  electric forces dominate the deceleration of the plasma. In the second phase, Lorentz forces dominate the deceleration over the electric forces. While in the first two phases the transition is isothermal for protons, the third phase on the downstream side represents the region in which the particle thermalization due to streaming plasma instabilities takes place. We focus on the first phase of the shock transition (the demagnetization region). We illustrate the three-phase shock in Fig.~\ref{sketch}.
\begin{figure}
  \resizebox{\hsize}{!}{\includegraphics{sketch.eps}}
  \caption{Illustration of the three-phase shock. Phase I: electric forces dominate the deceleration of the plasma (demagnetization region); Phase II: magnetic Lorentz forces dominate; Phase III: quasi-equilibrium relaxation of the plasma and thermalization through plasma instabilities.}
  \label{sketch}
\end{figure}

We first estimate the typical field conditions at a shock of this nature. In a sufficiently small surface element of the termination shock, the curvature of the shock surface is negligible. We approximate this surface element under consideration as a planar surface and, consequently, assume a one-dimensional geometry of the system. In this approximation, the shock parameters only depend on the coordinate $z$ parallel to the shock normal.
We assume a representation of the ion bulk-flow velocity $U$ as a function of  $z$ of the form
\begin{equation}\label{bulk}
U(z)=U_{1}-\frac{1}{2}\Delta U \left[ 1+\tanh \left(\frac{ z}{\alpha }\right)\right]
\end{equation}
where $\Delta U=U_{1}-U_{2}$. Indexes  1 and 2 characterize the asymptotic
upstream and downstream conditions of the flow, respectively. The quantity $\alpha $
determines the transition length at the shock. We choose the coordinate system in which $\vec B$ lies in the x-z-plane.
In addition, we assume that the $\vec B$-field is frozen into the main momentum flow  (i.e,. the ion flow). Therefore, the frozen-in condition under steady-state conditions, $\partial \vec{B}/\partial t=\nabla\times [\vec{U}\times \vec{B}]=0$,  yields
\begin{equation}
\vec{B}=\vec{B}_{x}(z)+\vec{B}_{z}
\end{equation}
in our one-dimensional coordinate system, where $B_{z}=B_{z1}=\mathrm{const.}$ and $B_{x}=B_{x1} U_{1}/U(z)$.

The $z$-component of the ion momentum equation is written as
\begin{equation}\label{momentum}
m_{\mathrm p}U\frac{\mathrm dU}{\mathrm dz}=-\frac{1}{n_{\mathrm p}}\frac{\mathrm dP_{\mathrm p}}{\mathrm dz} +eE,
\end{equation}
where $m_{\mathrm p}$ is the proton mass $e$ is the elementary charge, $E$ is the electric field that is parallel to the flow direction $\hat{\vec e}_z$, $P_{\mathrm p}$ is the proton pressure, and $n_{\mathrm p}$ is the proton density  \citep[see also][]{verscharen08,fahr13}.
In our isothermal approximation ($T_{\mathrm p}=\mathrm{const.}$) for the first phase of the shock transition, we find for the pressure-gradient term in the momentum equation by applying the flux-conservation requirement  $\Phi _{\mathrm e}=\Phi _{\mathrm p}=n_{\mathrm p}U=\mathrm{const}$ that
\begin{equation}
-\frac{1}{n_{\mathrm p}}\frac{\mathrm dP_{\mathrm p}}{\mathrm dz}=\frac{k_{\mathrm B}T_{\mathrm p}}{U}\frac{\mathrm dU}{\mathrm dz},
\end{equation}
where $k_{\mathrm B}$ is the Boltzmann constant.
This leads to the following expression for the electric field
\begin{equation}\label{Ecompl}
E=\frac{m_{\mathrm p}}{e}U\frac{\mathrm dU}{\mathrm dz}\left(1-\frac{1}{2}\frac{c_{\mathrm p}^2}{U^2}\right),
\end{equation}
where $c_{\mathrm p}=\sqrt{2k_{\mathrm B}T_{\mathrm p}/m_{\mathrm p}}$ is the ion thermal speed.
The correction term $c_{\mathrm p}^2/U^2$ on the right-hand side is small compared to 1 for high-Mach-number shocks in the first phase of the shock transition, which justifies neglecting the pressure gradient.
With Eq.~(\ref{bulk}), the electric field as a function of $z$ is given by
\begin{equation}
E=-\frac{m_{\mathrm p}}{2e\alpha}U\,\Delta U\mathrm{sech}^2\left(\frac{z}{\alpha}\right).
\end{equation}

The electric field $E$ is self-consistently generated by the charge separation between protons and electrons according to Gauss' law:
\begin{equation}
\frac{\mathrm dE}{\mathrm dz}=4\pi e\left(n_{\mathrm p}-n_{\mathrm e}\right).
\end{equation}

\section{Dominance of electric forces over Lorentz forces}\label{electric}

We can assume that, as long as the electric force is much stronger than
the Lorentz force on the electrons, the motion of a single electron with  electron velocity $v_{\mathrm ez}$ corresponds to a permanent linear acceleration according to
\begin{equation}
m_{\mathrm e}v_{\mathrm ez}\frac{\mathrm d v_{\mathrm ez}}{\mathrm dz}=-eE.
\end{equation}
Integrating leads to
\begin{equation}
\left. \frac{1}{2}m_{\mathrm e}v_{\mathrm ez}^{2}\right| _{-\Delta }^{+\Delta
}=-\int\limits _{-\Delta }^{+\Delta }eE\,\mathrm dz.
\end{equation}
We assume that the initial energy of an electron with a velocity equal to the bulk-flow velocity on the upstream side, $\epsilon _{\mathrm e0}=m_{\mathrm e}U_{1}^{2}/2$, is negligible compared to the resulting downstream energy. We note that protons and electrons have equal upstream bulk velocities $U_1$ as well as equal upstream densities $n_{1}$. With this assumption and Eq.~(\ref{Ecompl}) in the limit $c_{\mathrm p}^2/U^2\ll 1$, we find for the $z$-dependent kinetic energy of the electron that
\begin{multline}\label{DSE}
\frac{1}{2}m_{\mathrm e}v_{\mathrm ez}^{2}=-\int\limits _{-\infty }^{z}eE(z^{\prime})\mathrm dz^{\prime}=\frac{m_{\mathrm p}}{2\alpha}\Delta U\int\limits_{-\infty}^{z}U(z^{\prime})\mathrm{sech}^2\left(\frac{z^{\prime}}{\alpha}\right)\mathrm dz^{\prime}\\
=\frac{1}{2}m_{\mathrm p}\Delta U\left\{\left(U_1-\frac{\Delta U}{2}\right)\left[1+\tanh\left(\frac{z}{\alpha}\right)\right]+\frac{\Delta U}{4}\mathrm{sech}^2\left(\frac{z}{\alpha}\right)\right\}.
\end{multline}
For the downstream electron kinetic energy, we obtain with Eq.~(\ref{DSE}) that
\begin{multline}
\frac{1}{2}m_{\mathrm e}v_{\mathrm ez}^{2}=\frac{m_{\mathrm p}}{2\alpha}\Delta U\int\limits_{-\infty}^{+\infty}U(z)\mathrm{sech}^2\left(\frac{z}{\alpha}\right)\mathrm dz\\
= \frac{1}{2}m_{\mathrm p} U_1^2\left(\frac{s^2-1}{s^2}\right)=0.84\frac{1}{2}m_{\mathrm p}U_{1}^2,
\end{multline}
with the observed compression ratio $s=U_{1}/U_{2}\approx 2.5$
\citep[see][]{richardson08}.
Therefore, the shock electrons may consume about 84\% of the upstream proton kinetic energy, which corresponds to about 0.7~keV or $v_{\mathrm ez2}\approx  10^{9}\,\mathrm{cm/s}$, respectively.

We choose an electron with $v_{\mathrm ez1}=U_1$ and follow its trajectory through the shock. This electron fulfills $v_{\mathrm ez}(z)=U_{\mathrm e}(z)$ throughout its transit.
With this condition and flux conservation, $n_{\mathrm e}U_{\mathrm e}=\mathrm{const.}$, we directly determine the upstream-to-downstream density ratio for the electrons as
\begin{equation}
\frac{n_{\mathrm e2}}{n_1}=\frac{U_1}{U_{\mathrm e2}}=\sqrt{\frac{m_{\mathrm e}}{m_{\mathrm p}}}\sqrt{\frac{s^2}{s^2-1}}=0.025,
\end{equation}
which translates to $U_{\mathrm e2}=40U_{1}$.

The electric acceleration of
the electron, however, only dominates compared to the acceleration due to Lorentz forces as long as
\begin{equation}\label{Edom}
\left|eE\right|\geq \left|\frac{e}{c}v_{\mathrm ez}B_{x}\right|.
\end{equation}
This condition limits the resulting velocity to
\begin{equation}
\left|v_{\mathrm ez}\right|\leq \left| \frac{c}{B_{x}}E\right|=\left| \frac{1}{2\alpha} \frac{\Delta U}{U_1}\frac{m_{\mathrm p}c}{eB_{x1}}U^2\mathrm{sech}^2\left(\frac{z}{\alpha}\right)\right|,
\end{equation}
which we write in normalized units as
\begin{equation}\label{thr}
\left|\frac{v_{\mathrm ez}}{U_1}\right|\le\left|\frac{1}{2}\frac{U_1}{\Omega_{\mathrm px1}\alpha}\left(\frac{U}{U_1}\right)^2\frac{\Delta U}{U_1}\mathrm{sech}^2\left(\frac{z}{\alpha}\right)\right|,
\end{equation}
where $\Omega_{\mathrm px1}=eB_{x1}/m_{\mathrm p}c$. We define the normalized threshold velocity $v_{\mathrm{thr}}/U_1$ as the right-hand side of Eq.~(\ref{thr}).
Beyond this limit, Lorentz forces begin to compete with electric forces and need to be taken into account. We show this condition in Fig.~\ref{lorentz_cond}.
\begin{figure}
  \resizebox{\hsize}{!}{\includegraphics{lorentz_cond.eps}}
  \caption{Electron bulk velocity $U_{\mathrm e}$ from Eq.~(\ref{DSE}), threshold velocity $v_{\mathrm{thr}}$ from Eq.~(\ref{thr}), and $U_{\mathrm e}+c_{\mathrm e 1}$, where $c_{\mathrm e1}=\sqrt{2k_{\mathrm B}T_{\mathrm e1}/m_{\mathrm e}}$ is the upstream electron thermal speed, as functions of the spatial coordinate $z$. We use the following parameters: $s=2.5$, $U_1=4\times 10^7\,\mathrm{cm/s}$, $B_{x1}=2\times 10^{-7}\,\mathrm G$, $\alpha=10^8\,\mathrm{cm}$, and $T_{\mathrm e1}=2\times 10^4\,\mathrm K$.}
  \label{lorentz_cond}
\end{figure}
The black solid line shows the threshold velocity $v_{\mathrm {thr}}$ in units of the upstream bulk velocity $U_1$.  Considering again an electron with $v_{\mathrm ez}(z)\approx U_{\mathrm e}(z)$, we determine the profile of $U_{\mathrm e}$ as a function of $z$ from Eq.~(\ref{DSE}). We show this profile as the red dashed line in Fig~\ref{lorentz_cond}.  The bulk velocity fulfills the condition in Eq.~(\ref{thr}) in the range between $z\approx-1.5\alpha$ and $z\approx -0.1\alpha$. We also plot the sum of the electron bulk velocity and the upstream electron thermal speed as the green dashed line. This line indicates the region in which the majority of the solar-wind electrons fulfill Eq.~(\ref{thr}). This condition is fulfilled between $z\approx -1.2\alpha$ and $z\approx -0.1\alpha$. As these considerations show, the demagnetization region is restricted to the region $z<0$.

In the limit of a quasi-parallel shock, $\Omega_{\mathrm px1}\rightarrow 0$ which leads to $v_{\mathrm{thr}}\rightarrow\infty$. Therefore, the quasi-parallel shock is solely a result of the electrostatic interaction while Lorentz forces are negligible. On the other hand, a quasi-perpendicular shock still has a finite $v_{\mathrm{thr}}$, which can be $>U_{\mathrm e}$ depending on the finite shock parameters $U_1$, $\Omega_{\mathrm px1}$, and $\alpha$.

\section{Motion under first-order Lorentz forces}\label{motion}

We assume that, within a finite region, the electric forces are strongly dominant according to Eq.~(\ref{Edom}).
The zeroth-order motion of an electron is then determined by
\begin{equation}
\frac{\mathrm d\vec v_{\mathrm e}^{0}}{\mathrm dt}=-\frac{e}{m_{\mathrm e}} \vec E,
\end{equation}
while the first-order Lorentz force is given by
\begin{equation}
\vec{K}_{\mathrm L}^{1}=-\frac{e}{c}\left[\vec{v}_{\mathrm e}^{0}\times \vec{B}\right].
\end{equation}

Under the assumption that the electric field is purely oriented along the $z$-direction, we obtain
\begin{equation}
{v_{\mathrm e}^{0}}^2(z)-{v_{\mathrm e1}^{0}}^{2}=-2\int\limits_{-\infty}^{z}\frac{e}{m_{\mathrm e}}E(z^{\prime})\,\mathrm dz^{\prime}.
\end{equation}
We furthermore assume that the upstream electron velocity $v_{\mathrm e1}^{0}$
is negligible.
Now we allow for a first-order velocity component in the $y$-direction (i.e., the $\hat{\vec e}_z\times \hat{\vec{b}}$-direction, where $\hat{\vec b}=\vec B/|\vec B|$) due to the first-order Lorentz force $K_{\mathrm L}^1$:
\begin{equation}\label{KL}
v_{\mathrm ey}^{1}=\int\limits_{-\infty}^z \frac{1}{m_{\mathrm e}}K_{\mathrm L}^{1}(z^{\prime})\frac{\mathrm dt}{\mathrm dz^{\prime}}\,\mathrm dz^{\prime},
\end{equation}
where $K_{\mathrm L}^1$ is given by
\begin{equation}
K_{\mathrm L}^{1}(s)=-\frac{e}{c}\left|\left[\vec{v}_{\mathrm e}^{0}(z)\times \vec{B}(z)\right]\right|=-\frac{e}{c}B_x(z)\sqrt{-\frac{2e}{m_{\mathrm e}}\int\limits_{-\infty}^{z}E(z^{\prime})\,\mathrm dz^{\prime}}.
\end{equation}
With Eq.~(\ref{KL}), we find
\begin{equation}
v_{\mathrm ey}^{1}=-\int\limits_{-\infty}^z \frac{e}{m_{\mathrm e}c}B_x(z^{\prime})\sqrt{-\frac{2e}{m_{\mathrm e}}\int\limits_{-\infty}^{z^{\prime}}E(z^{\prime\prime})\,\mathrm dz^{\prime\prime}}%
\frac{\mathrm dt}{\mathrm dz^{\prime}}\,\mathrm dz^{\prime}.
\end{equation}
Using $\mathrm dt/\mathrm dz=1/v_{\mathrm e}^{0}(z)$, we obtain the solution
\begin{equation}
v_{\mathrm ey}^{1}=-\int\limits_{-\infty}^{z}\frac{e}{m_{\mathrm e}c}B_x(z^{\prime})\,\mathrm dz^{\prime}=-\int\limits_{-\infty}^{z}\Omega_{\mathrm e x}(z^{\prime})\,\mathrm dz^{\prime}=-z\left\langle \Omega_{\mathrm ex}\right\rangle _{z},
\end{equation}
where $\Omega_{\mathrm ex}(z)$ denotes the local electron gyro-frequency based on $B_x$, and $\left\langle
\Omega_{\mathrm ex}\right\rangle _{z}$ is the average gyro-frequency in the integration domain.

It is noteworthy that the electron motion in the $y$-direction does not depend on the electric-field configuration in the transition region as long as the Lorentz force can be considered to be of higher
order than the electric force.
Therefore, the curling of the particle trajectory around the frozen-in field only happens after a time $\tau_{\mathrm{Le}}$ when the action of the Lorentz force changed the velocity by about its magnitude,
$(e/c)v_{\mathrm ez}B_{x} \tau _{\mathrm{Le}}\approx m_{\mathrm e}v_{\mathrm ez}$, yielding $\tau _{\mathrm{Le}}\approx m_{\mathrm e}c/eB_x=1/\Omega _{\mathrm ex}$. During this period, the electron moves in the $z$-direction by the amount $\Delta z=\tau_{\mathrm{Le}}v_{\mathrm ez}=v_{\mathrm ez}/\Omega _{\mathrm ex}$.
With our results for $v_{\mathrm ez}$ from Sect.~\ref{electric}, we find that the permitted extent of the overall transition region can roughly be
estimated as $\Delta z\approx v_{\mathrm ez}/\Omega _{\mathrm ex}\approx  10^{8}\,\mathrm{cm}$.

\section{Discussion and conclusions}

In this paper, we study the action of shock-electric fields on
electrons and ions entering from the upstream regime of the shock into the shock
transition region. We split the latter into three consecutive regions:
the first is demagnetized and operates like an electric double layer
where the action of the upflashing shock-electric field strongly dominates
over Lorentz forces. In the second region, Lorentz forces due to the piled-up
magnetic fields compete with or even dominate compared to electric forces. In the third phase, the overshooting energy is transferred into proton and electron heating by the action of kinetic plasma instabilities.

In Sect.~\ref{electric}, we determine the electron velocity profile. The downstream electron velocity before entering phase three of the shock is about 40 times the upstream velocity, while the downstream electron density is about 0.025 the upstream density. The shocked plasma opens up a new thermodynamic degree
of freedom allowing that a substantial fraction (about 84 percent) of the
upstream ion kinetic energy gets stored in the overshoot velocities of the
electrons which, in the second region of the shock transition, is deposited into a shell in velocity space
centered around the downstream ion bulk velocity. In that case, a major fraction of the upstream ion bulk-flow energy is available to be transferred into electron thermal energy.
We determine the threshold velocity $v_{\mathrm{thr}}$, which defines the separation between phases one and two of the shock transition. The electron motion is dominated by electric forces for all electrons with a velocity $<v_{\mathrm{thr}}$. For typical termination-shock parameters, we find that the region in which this condition is fulfilled for electrons with the electron bulk velocity extends from $z=-1.5\alpha$ to $z=-0.1\alpha$.

In Sect.~\ref{motion}, we treat Lorentz-forces as a first-order correction to the plasma deceleration. We show that the first-order shock-perpendicular component of the electron velocity is determined by the average electron gyro-frequency in the shock layer. With this result, we determine the size of the demagnetization region (phase one) as $\Delta z\approx 10^8\,\mathrm{cm}$ for typical solar-wind parameters.

In earlier work  \citep[see][]{chalov13,fahr13,fahr15}, we emphasized the fact that the creation of energetic electrons at the plasma passage over the shock is
essential to fulfill the thermodynamic entropy requirements \citep{fahr15} and to arrive at downstream plasma properties that nicely fit the Voyager-2 measurements of \citet{richardson08}. \citet{zieger15} recently published a multifluid study of the solar-wind termination shock, which strongly supports our claim for the occurrence of energetic downstream
electrons by showing that the Voyager-2 measurements allow for a reasonably good theoretical fit with the model results only if the appearance of
energetic downstream electrons is taken into account.

The remaining question as to why these predicted energetic
electrons were not detected by the Voyager plasma analyzers,
can most probably be answered along the argumentation developed in a recent
paper by \citet{fahr15b}. These authors show that the energetic electrons create
a strongly increased electric charge-up of the spacecraft detectors and
thereby, owing to electric screening,  impede the detectors in measuring
countable fluxes of these electrons.

We intend to continue the studies presented in this letter by kinetically analyzing the physics relevant for phases II and III of the three-phase shock.

\begin{acknowledgements}

We appreciate helpful discussions with Marty Lee. This work was supported in part by NSF/SHINE grant AGS-1460190.

\end{acknowledgements}

\bibliographystyle{aa}


\section{Introduction}

Functional data analysis (FDA) has become a popular and useful framework for applications in which the unit of measurement is a function, curve or image. Conceptually, FDA leverages the underlying data structure, often temporal or spatial, to improve understanding of patterns and variation. A wide array of tools have been developed for the functional data setting, for example, functional principal component analysis (FPCA) and regression models using functional responses \citep{ramsay2005, morris2015, sorensen2013}. The basic unit of observation is the curve $Y_i(t)$ for subjects $i \in \ldots, I$ in the cross-sectional setting and $Y_{ij}(t)$ for subject $i$ at visit $j \in \ldots, J_i$ for the multilevel or longitudinal structure. Methods for functional data are typically presented in terms of continuous functions, but in practice data are observed on a discrete grid that may be sparse or dense at the subject level and that may be the same across subjects or irregular.

Many methods for FDA have standard visualization approaches that clarify the results of analyses; examples include scree plots for FPCA and coefficient function plots for function on scalar regression. Clear visualizations aid in exploratory analysis and help to communicate results to non-statistical collaborators. However, creating useful plots is often time consuming and must be repeated each time a model is changed, and no software currently exists to facilitate this process.

The \texttt{refund.shiny} package \citep{refund.shiny} creates interactive visualizations for functional data analyses, allowing researchers to create common graphics for standard analyses with just a few lines of code. Currently, \texttt{refund.shiny} builds plots for functional principal component analysis (FPCA), multilevel FPCA (MFPCA), time-varying FPCA (TV-FPCA), and function-on-scalar regression (FoSR). The workflow separates analysis and visualization steps: analyses are performed by functions in the \texttt{refund} package \citep{refund} and interactive visualizations are generated by the \texttt{plot\_shiny()} function in the \texttt{refund.shiny} package. Changes to the analysis -- increasing the number of retained principal components, for example, or augmenting a regression model with new predictors -- are easily incorporated into the graphical interface. User interaction with the displayed graphics facilitates comparisons and streamlines navigation between visualizations.

We illustrate the tools  in \texttt{refund.shiny} using a single dataset, which we describe briefly here. The diffusion tensor imaging (\texttt{DTI}) dataset available in the \texttt{refund} package includes cerebral white matter tracts for multiple sclerosis patients and healthy controls. White matter tracts are collections of axons, projections of neurons that transmit electrical signals and are coated by a fatty substance called myelin \citep{greven2010, goldsmith2011a, staicu2012modeling}. DTI is a magnetic resonance imaging modality that measures diffusion of water in the brain; because water movement is restricted in white matter fibers, DTI allows the quantification of white matter tract integrity. The \texttt{DTI} dataset contains tract profiles -- continuous summaries of tract properties along their major axis -- for 142 subjects across multiple visits, with a median of 4 scans per subject. The dataset includes tract profiles for several tracts, the PASAT score (a continuous variable that indicates brain reactivity and attention span), subject sex, subject ID, visit number, and time of visit \citep{strauss2006}. Because we observe tract profiles for each subject over time, the DTI dataset is a functional dataset with longitudinal structure; in order to use the same dataset across examples we sometimes neglect this structure or subset the data. The following code can be used to install \texttt{refund} and \texttt{refund.shiny} and load the \texttt{DTI} data:

\setstretch{1.1}
\begin{verbatim}
> install.packages("refund.shiny")
> library(refund.shiny)
> library(refund)
> data(DTI)
\end{verbatim}
\setstretch{1.8}

Sections \ref{sec:fpca}, \ref{sec:mfpca}, \ref{sec:tfpca}, and \ref{sec:fosr} each provide a brief methodological overview of an analysis technique for FDA and describe the corresponding interactive visualization tools in the \texttt{refund.shiny} package. Section \ref{sec:structure} details the structure of the \texttt{refund.shiny} package. We close in section \ref{sec:conc} with a discussion.

\section{Functional Principal Component Analysis}
\label{sec:fpca}

We start with FPCA, one of the most common exploratory tools for functional datasets.

\subsection{FPCA Model}
\label{subsec:fpca_mod}

FPCA characterizes modes of variability by decomposing functional observations into population level basis functions and subject-specific scores \citep{ramsay2005}. The basis functions have a clear interpretation, analogous to that of PCA: the first basis function explains the largest direction of variation, and each subsequent basis function describes less. The FPCA model is typically written
\beqn
\label{eq:fpca_model}
	Y_i(t) = \mu(t) + \sum_{k=1}^{K} c_{ik}\psi_{k}(t) + \epsilon_i(t)
\eeqn
\noindent
where $\mu(t)$ is the population mean,  $\psi_{k}(t)$ are a set of orthonormal population-level basis functions, $c_{ik}$ are subject-specific scores with mean zero and variance $\lambda_k$, and $\epsilon_i(t)$ are residual curves. Estimated basis functions $\widehat{\psi}_1(t), \widehat{\psi}_2(t), \ldots, \widehat{\psi}_{K}(t)$ and corresponding variances $\widehat{\lambda}_1 \geq \widehat{\lambda}_2 \geq \ldots \geq \widehat{\lambda}_K$ are obtained from a truncated Karhunen-Lo\`eve decomposition of the sample covariance $\widehat{\Sigma}(s,t) = \widehat{\mbox{Cov}}\left(Y_i(s), Y_i(t)\right)$. In practice, the covariance $\widehat{\Sigma}(s,t)$ is often smoothed using a bivariate smoother that omits entries on the main diagonal to avoid a ``nugget effect" attributable to measurement error, and scores are estimated in a mixed model framework \citep{yao2005, goldsmith2013b}. The truncation lag $K$ is often chosen so that the resulting approximation accounts for at least 95\% of observed variance.

\subsection{Graphics for FPCA}
\label{subsec:fpca_graph}

Our example uses the \texttt{fpca.sc()} function from the \texttt{refund} package. Several other implementations of FPCA are available in \texttt{refund}, including \texttt{fpca.face()}, \texttt{fpca.ssvd()}, and \texttt{fpca2s()}, all of which are compatible with \texttt{refund.shiny}. The number of functional principal components (FPCs) is chosen by percent variance explained, with the default set to 0.99. See \texttt{?plot\_shiny} for examples.

Graphics for FPCA are implemented by the code below:

\setstretch{1.1}
\begin{verbatim}
> fit.fpca = fpca.sc(Y = DTI$cca)
> plot_shiny(obj = fit.fpca)
\end{verbatim}
\setstretch{1.8}
\noindent
Executing this code produces a user interface with five tabs. The first tab shows  $\widehat{\mu}(t) \pm \sqrt{\widehat{\lambda}_k} \widehat{\psi}_k(t)$, and includes a drop-down menu through which the user can select $k$ (an example for a similar tab, based on multilevel data, is shown in Section \ref{sec:mfpca}). The second tab presents static scree plots of the eigenvalues $\widehat{\lambda}_k$ and the percent variance explained by each eigenvalue. The third tab shows $\widehat{\mu}(t) + \sum_{k=1}^{K} c_{k} \widehat{\psi}_k(t)$, and includes slider bars through which the values of $c_{k}$ can be set; adjusting the sliders allows the user to see a fitted curve for a hypothetical subject with the selected combination of scores. The fourth tab allows users to assess quality-of-fit by plotting fitted and observed values for any subject in the dataset.

\begin{figure}[h]
  \centering
     \begin{tabular}{cc}
       \includegraphics[width=1\textwidth]{tab1_fpca.png}
     \end{tabular}
     \caption{Screenshot showing tab 5 of the interactive graphics for FPCA. A scatterplot of FPC loadings $\widehat{c}_{ik}$ against $\widehat{c}_{ik'}$ is shown in the upper plot, and $k$ and $k'$ are selected using drop-down menus at the left. The lower plot shows fitted curves for all subjects. In the scatterplot, a subset of estimated loadings can be selected by clicking-and-dragging to create a blue box; blue curves in the plot of fitted values correspond to selected points in upper plot.}
    \label{fig:fpca}
\end{figure}

The fifth tab for the interactive graphic produced by the code above is shown as a static plot in Figure~\ref{fig:fpca}. A scatterplot of estimated FPC loadings $\widehat{c}_{ik}$ against $\widehat{c}_{ik'}$ is shown in the upper plot, and $k$ and $k'$ are selected using drop-down menus at the left. The lower plot shows fitted curves for all subjects. In the scatterplot, a subset of FPC loadings can be selected by clicking-and-dragging to create a blue box; blue curves in the plot of fitted values correspond to selected subjects in upper plot. In Figure~\ref{fig:fpca} the first and second FPCs are selected for the $x$ and $y$ axes of the score plot, respectively, and several subjects that have negative values for FPC 1 are highlighted. Fitted values for these subjects are clustered at the top of the $y$-axis, indicating that the first FPC largely represents a vertical shift from the mean. A working example of \texttt{refund.shiny} for FPCA on a different dataset is available at \href{https://jeff-goldsmith.shinyapps.io/FPCA}{https://jeff-goldsmith.shinyapps.io/FPCA}.

\section{Multilevel Functional Principal Components Analysis}
\label{sec:mfpca}

Multilevel functional principal component analysis (MFPCA) extends the ideas of FPCA to functional data with a multilevel structure.

\subsection{MFPCA Model}
\label{subsec:mfpca_mod}

Multilevel functional data are increasingly common in practice; in the case of our DTI example, this structure arises from multiple clinical visits made by each subject. MFPCA models the within-subject correlation induced by repeated measures as well as the between-subject correlation modeled by classic FPCA. This leads to a two-level FPC decomposition, where level 1 concerns subject-specific effects and level 2 concerns visit-specific effects. Population-level basis functions and subject-specific scores are calculated for both levels \citep{di2009, di2014}. The MFPCA model is:
\beqn
\label{eq:mpca_model}
	Y_{ij}(t) = \mu(t) + \eta_j(t) + \sum_{k_1=1}^{K_1} 	c_{ik}^{(1)}\psi_{k}^{(1)}(t) + \sum_{k_2=1}^{K_2}c^{(2)}_{ijk}\psi_{k}^{(2)}(t) + \epsilon_{ij}(t)
\eeqn

\noindent
where $\mu(t)$ is the population mean,  $\eta_j(t)$ is the visit-specific shift from the overall mean, $\psi_{k}^{(1)}(t)$ and $\psi_{k}^{(2)}(t)$ are the eigenfunctions for levels 1 and 2, respectively, and $c_{ik}^{(1)}$ and $c^{(2)}_{ijk}$ are the subject-specific and subject-visit-specific scores. Often, visit-specific means $\eta_j(t)$ are not of interest and can be omitted from the model. Estimation for MFPCA extends the approach for FPCA: estimated between- and within-covariances $\widehat{\Sigma}^{(1)}(s,t) = \widehat{\mbox{Cov}}(Y_{ij}(s), Y_{ij'}(t))$ for $j \neq j'$ and $\widehat{\Sigma}^{(2)}(s,t) = \widehat{\mbox{Cov}}(Y_{ij}(s), Y_{ij}(t))$ are derived from the observed data, smoothed, and decomposed to obtain eigenfunctions and values. Given these objects, scores are estimated in a mixed-model framework.

\subsection{Graphics for MFPCA}
\label{subsec:mfpca_graph}

MFPCA is implemented in the \texttt{mfpca.sc()} function from the \texttt{refund} package. By default, \texttt{mfpca.sc()} does not calculate visit-means, but they can be calculated by specifying the \texttt{mfpca.sc()} argument \texttt{twoway = TRUE}.

Graphics for MFPCA are implemented by the code below:

\setstretch{1.1}
\begin{verbatim}
 > Y = DTI$cca
 > id = DTI$ID
 > fit.mfpca = mfpca.sc(Y = Y, id = id, twoway = FALSE)
 > plot_shiny(fit.mfpca)
\end{verbatim}
\setstretch{1.8}

\noindent
This code produces an interface with five tabs, which is similar to the interface for FPCA but includes features unique to multilevel analyses. Tabs 1, 2, 3, and 5 for MFPCA are $\widehat{\mu}(t) \pm \sqrt{\widehat{\lambda}_{k_L}^{(L)}} \, \widehat{\psi}_{k_L}^{(L)}(t)$, static scree plots of the estimated eigenvalues $\widehat{\lambda}_{k_L}^{(L)}$, $\widehat{\mu}(t) + \sum_{k_L=1}^{K_L} c_{k_L} ^{(L)} \widehat{\psi}_{k_L}^{(L)}(t)$, and scatterplots of FPC scores (similar to Figure \ref{fig:fpca}), respectively. These mirror the tabs for FPCA and include inset sub-tabs to toggle between level, $L$, to display results for level 1 or level 2. The fourth tab plots fitted and observed values for any user-selected subject in the dataset; the user can display all visits for the selected subject or choose a subset of visits. The first tab for the interactive visualization produced by the code above is displayed in Figure~\ref{fig:mfpca}, and shows $\widehat{\mu}(t) \pm \sqrt{\widehat{\lambda}_{2}^{(1)}} \, \widehat{\psi}_{2}^{(1)}(t)$.

\begin{figure}[h]
  \centering
     \begin{tabular}{cc}
       \includegraphics[width=1\textwidth]{mfpccTab1.png}
     \end{tabular}
     \caption{Screenshot showing tab 1 of the interactive graphic for MFPCA. The plot at right shows $\widehat{\mu}(t) \pm \sqrt{\widehat{\lambda}_{k_L}^{(L)}} \widehat{\psi}_{k_L}^{(L)}(t)$; $k_L$ is chosen by the drop-down menu in at left, and the user can switch between level $L$ by clicking \textit{Level 1} or \textit{Level 2} inset tabs at the top left.}
    \label{fig:mfpca}
\end{figure}

\section{Time-varying Functional Principal Component Analysis} \label{sec:tfpca}

Time-varying functional principal component analysis (TV-FPCA) extends the ideas of FPCA to model functional data that are observed repeatedly in a longitudinal framework. In contrast to MFPCA, TV-FPCA accounts for the actual time of visit $T_{ij}$ at which the functional object $Y_{ij}(\cdot)$ is recorded; this allows us to study the time-varying behavior of the underlying true process and make predictions of full trajectory at an unobserved visit time \citep{park2015lfda}. Other modeling methods for longitudinal functional data that incorporate the actual visit times $T_{ij}$ include \citet{greven2010} and \citet{chen2013repeated}.

\subsection{TV-FPCA Model} \label{sec: TFPCA model}

TV-FPCA \citep{park2015lfda} model for $Y_{ij}(t) = Y_i(t, T_{ij})$ is given as follows:
\begin{equation}
Y_{ij}(t) = \mu(t, T_{ij})  + \sum_{k = 1}^{K} c_{ik}(T_{ij}) \psi_{k}(t) + \epsilon_{ij}(t),
\end{equation}
where $\mu(t, T_{ij})$ is the population mean that is assumed to vary smoothly over $t$ and visit time $T_{ij}$, $\psi_k(t)$ are orthogonal basis functions,
$c_{ik}(T_{ij})$ are corresponding loadings that vary over $T_{ij}$ with mean zero and variance $\lambda_{k}$, and $\epsilon_{ij}(t)$ are residual curves. The time-varying scores $c_{ik}(t_{ij})$ are uncorrelated over $i$, but correlated over $j$. Estimation of the TV-FPCA  model components entails: 1) estimation of the population mean by using bi-variate smoothing, 2) estimation of the marginal covariance $ \Sigma(s,t) = \int {\mbox{Cov}} \{Y_i(s, T), Y_i(t,T)\} g(T) \,dT$, where $g(T)$ is the density of the $T_{ij}$'s using the observed data, smoothing and decomposing it to get the eigenfunctions/eigenvalues $\widehat \psi_k(t)$ and $\widehat \lambda_k$; 3) estimation of the $k$th component covariance $\widehat G_k (T, T') = {\mbox{Cov}} \{c_{ik}(T) c_{ik}(T')\}$. The last step is carried out using either linear random effects, implying $c_{ik}(T) = b^{(k)}_{0i} + b^{(k)}_{1i}T$ or FPCA implying $c_{ik}(T)  = b_{ik1}\phi_{k1}(T)+\ldots + b_{ikL_k}\phi_{kL_k}(T)
$.
By modeling these longitudinal dynamics, the time-varying coefficient function $c_{ik}(\cdot)$ can be used to predict scores at any longitudinal time $T$ and, as a result, to predict the full response trajectory $Y_i(\cdot, T)$.

\subsection{Graphics for TV-FPCA}

TV-FPCA is implemented in the \verb|fpca.lfda()| function in the \verb|refund| package. In Section \ref{sec: TFPCA model}, we have used $t$ to denote the functional argument for consistency with the rest of the paper; however to maintain consistency with the notations used in \citet{park2015lfda}, the \verb|plot_shiny()| function for TV-FPCA uses $s$ to denote the functional argument and $T$ to denote the longitudinal time.

Graphics for TV-FPCA are implemented by the code below:

\linespread{1}
\begin{verbatim}
> MS <- subset(DTI, case ==1)
> index.na <- which(is.na(MS$cca)); Y <- MS$cca; Y[index.na] <- fpca.sc(Y)$Yhat[index.na]
> id <- MS$ID
> visit.index <- MS$visit
> visit.time <- MS$visit.time/max(MS$visit.time)
> fit.tfpca <- fpca.lfda(Y = Y, subject.index = id,
+                        visit.index = visit.index, obsT = visit.time,
+                        LongiModel.method = `lme')
> plot_shiny(fit.tfpca)
\end{verbatim}
\linespread{1.6}
The code produces an interface with two tabs. Tab 1 shows exploratory plots and includes three inset sub-tabs. The first sub-tab, shown in Figure \ref{fig:tfpca}, plots the observed curves for any user-selected subject, and includes options to display the observed curves of all subjects in the background and to display the estimated pointwise mean curve, denoted by $m(t)$. The second sub-tab allows the user to see the longitudinal changes of the observed curves for a user-selected subject $i$; a slider bar animates the subject's visit times and highlights the corresponding observed curve in the plot. The last sub-tab shows two plots of the actual visit times $T_{ij}$: the bottom plot presents static histogram of visit times of all subjects, while the top plot presents all of observed visit times on a horizontal line to help visualize the sparsity of the longitudinal sampling.

Tab 2 shows estimated model components and predictions, and includes 8 inset sub-tabs. Sub-tabs 1 and 2 present static images of the estimated mean surface $\widehat{\mu}(t, T)$ and estimated marginal covariance $\widehat{\Sigma}(s,t)$. Sub-tabs 3, 4, and 5 illustrate the first step of estimation, and plot estimates of eigenfunctions $\widehat{\psi}_k(t)$, $m(t) \pm 2\sqrt{\widehat{\lambda}_k}\widehat{\psi}_k(t)$, and static scree plots of the estimated eigenvalues $\widehat{\lambda}_k$, respectively. Sub-tab 6 shows the estimated covariance of the time-varying loadings $c_{ik}(\cdot)$ for user-specified $k$. Sub-tab 7 shows the prediction of $c_{ik}(T)$ for any user-selected subject $i$ and component $k$; it also has an option of displaying predicted values of $c_{ik}(T)$ for all subjects in the background. Lastly, sub-tab 8 shows the prediction of a full response trajectory $Y_i(\cdot, T)$ for user-selected subject $i$ in animation with change of values across $21$ equi-spaced grid of points of $T$ in the range of observed visit times of all subjects.

\begin{figure}[h]
  \centering
     \begin{tabular}{cc}
       \includegraphics[width=1\textwidth]{tfpca.png}
     \end{tabular}
     \caption{Screenshot showing Tab 1 of the interactive graphic for TV-FPCA. The plot shows observed data of the selected subject.}
    \label{fig:tfpca}
\end{figure}

\section{Function-on-Scalar Regression}
\label{sec:fosr}

In many cases, a length $p$ vector of scalar covariates $\bx_i= [x_{i1}, \ldots, x_{ip}]$ is observed in addition to the function $Y_i(t)$. In these situations, it is often of interest to model the conditional expectation of the functional response as it depends on the scalar predictors; indeed, this problem has been the focus of a large literature \citep{brumback1998, guo2002, morris2003, morris2006, reiss2010b, scheipl2015, goldsmith2015a, goldsmith2015b}.

\subsection{FoSR Model}
\label{subsec:fosr_mod}

The most common function-on-scalar regression model is
\beqn
\label{eq:fosr_model}
		Y_{i}(t)  =  \beta_{0}(t) + \sum_{k = 1}^{p} x_{ik} \beta_{k}(t) + \epsilon_{i}(t)
\eeqn
where the $\beta_{k}(t)$ are fixed effects associated with scalar covariates and the $\epsilon_{i}(t)$ are residual curves. The coefficients $\beta_k(t)$ are interpreted analogously to coefficients in a (non-functional) multiple linear regression -- as the expected change in response for each one unit change in the predictor -- with the exception that they, like the outcome, are defined over $t$. Many estimation and inferential strategies are available for model~(\ref{eq:fosr_model}); a popular approach is to expand coefficients $\beta_{k}(t)$ using a spline basis, which allows one to recast (\ref{eq:fosr_model}) as a traditional linear regression model and focus estimation on a vector of unknown spline coefficients. Our example uses the \texttt{bayes\_fosr()} function in the \texttt{refund} package, which uses a rich cubic B-spline basis and estimates spline coefficients in a Bayesian framework with priors specified to enforce smoothness in the resulting coefficient functions. Both a Gibbs sampler and a computationally efficient variational approximation are available in \texttt{refund}.

\subsection{Graphics for FoSR}
\label{subsec:fosr_graph}

Graphics for FoSR are implemented by the code below:

\setstretch{1.1}
\begin{verbatim}
> DTI = DTI[complete.cases(DTI),]
> fit.fosr = bayes_fosr(cca ~ pasat + sex, data = DTI)
> plot_shiny(fit.fosr)
\end{verbatim}
\setstretch{1.8}

\noindent
This code produces a interface with four tabs, each showing plots associated with model~\ref{eq:fosr_model}. The first tab is a plot of the observed data with the option to color curves by a user-selected covariate; this builds intuition analogously to scatterplots for non-functional regression. The second tab shows $\widehat{\beta}_0(t)+\sum_{k=1}^p x_{k} \widehat{\beta}_k(t)$, where values of $x_k$ can bet set by slider bars for continuous covariates or drop-down menus for categorical covariates; adjusting the sliders or drop-down menus shows the estimated conditional expectation for a specified predictor vector. The third tab, illustrated in Figure \ref{fig:fosr}, shows estimated coefficient functions $\widehat{\beta}_k(t)$ with pointwise confidence intervals for the covariate $x_k$ selected in a drop-down menu. The fourth tab is a plot of the residual curves $\widehat{\epsilon}_{i}(t)$ and allows for identification of median and outlying curves by band depth \citep{lopezpintado2009, sun2011, sun2012}; the user can also choose to 'rainbowize by depth', which colors the curves from the median outward based on depth.

\begin{figure}[h]
  \centering
     \begin{tabular}{cc}
       \includegraphics[width=1\textwidth]{FoSR_tab3.png}
     \end{tabular}
     \caption{Screenshot showing tab 3 of the interactive graphic for FoSR. The plot shows the estimated coefficient function $\widehat{\beta}_k(t)$ for the selected covariate $x_k$ with pointwise confidence intervals.}
    \label{fig:fosr}
\end{figure}

\section{Code Structure of the \texttt{refund.shiny} Package}
\label{sec:structure}

We now briefly describe the code infrastructure used to create the \texttt{refund.shiny} package.

As indicated in the introduction, the workflow separates visualization from analysis in the following way. First, one analyzes a dataset using a function in the \texttt{refund} package. The functions in \texttt{refund} take discretely observed functional data as input, perform an analysis, and return an object whose class corresponds to the method used. For example, the \texttt{fpca.sc()} function return as object of class \texttt{fpca} and the \texttt{bayes.fosr()} function returns an object of class \texttt{fosr}. The primary function in \texttt{refund.shiny}, \texttt{plot\_shiny()}, is a generic function whose behavior depends on the class of the object passed as an argument. Because of this structure, the user experience is uniform across a variety of analyses; this also suggests a development strategy for the addition of interactive graphics as new analysis techniques become available. Lastly, by separating the analysis and visualization steps, it is possible for analysis functions developed outside of the \texttt{refund} package to return objects of a defined class and thereby take advantage of the plotting capabilities we describe.

The interactive graphics in the \texttt{refund.shiny} are built on RStudio's \texttt{R} package \texttt{shiny} \citep{shiny}, which significantly reduces the barriers to producing webpage-style representations of analysis results in \texttt{R}. Other examples of interactive graphics that utilize the \texttt{shiny} framework are \texttt{shinyMethyl} \citep{fortin2014} for visualization of high-dimensional genomic data and \texttt{shinystan} \citep{shinystan} for exploring Bayesian models fit using Markov Chain Monte Carlo. In \texttt{refund.shiny} the plots within tabs are produced using \texttt{ggplot2} \citep{ggplot2}; it is possible to export each plot as a PDF or to save the corresponding \texttt{ggplot} object to the user's \texttt{R} workspace for further manipulation.

\section{Concluding Remarks}
\label{sec:conc}

Visualization has long been acknowledged as a central tool in data analysis. For functional datasets, the need for useful graphics is compounded: data are inherently complex, high-dimensional and structured. Although a robust literature for functional data exists and many methods have standard graphical representations, the creation of these graphics is often time consuming. The \texttt{refund.shiny} package was developed to ease this process by producing a visualization framework for several common functional data analyses. By leveraging new tools for interactivity, \texttt{refund.shiny} responds to user input and actions and, in so doing, can build intuition for analyses in both statisticians and practitioners. The interfaces produced by \texttt{refund.shiny} using the \texttt{shiny} framework are web applications, rendered locally by a web browser. These applications can be hosted publicly and may, in the spirit of  ``visuanimations`` \citep{genton2015}, be included as important parts of scientific papers and reports.

We use an analytic workflow that separates modeling from visualization. Doing so allows several methods and implementations to take advantage of the same visualization software; as an example, \texttt{fpca.sc()}, \texttt{fpca.face()}, \texttt{fpca.ssvd()}, and \texttt{fpca2s()} implement different methods for FPCA but are all compatible with \texttt{plot\_shiny()}. This produces an intuitive user experience and leaves open the possibility for future approaches to FPCA or FoSR to use the \texttt{refund.shiny} package for visualization with minimal effort. Similarly, this workflow is amenable to the development of interactive visualizations for additional functional data analyses in future iterations of the package.

\section{Acknowledgments}
\label{sec:ack}

The third author's research was supported partially by National Science Foundation DMS 0454942 and National Institutes of Health grants R01 NS085211 and R01 MH086633. The last author's research was supported in part by Award R01HL123407 from the National Heart, Lung, and Blood Institute and by Award R21EB018917 from the National Institute of Biomedical Imaging and Bioengineering. The MRI/DTI data were collected at Johns Hopkins University and the Kennedy-Krieger Institute.

\setstretch{1.2}

\bibliographystyle{jasa}