\section{Introduction}

Research on the concurrent and post-flight effects of in-flight exposures to altitude is an understudied, but burgeoning area of work in environmental health. Studies on the impact of altitude on heart health have found that arrhythmic events can occur during transitions and appear to have age-related differences \citep{Behn2014}. These cardiac health effects persist in flight conditions as \cite{Silva2016} demonstrate, showing that altered cardiovascular functioning in participants within 24 hours after commercial flight. \cite{Chatterjee2021} examine cardiac arrests in air travels both in-flight and post-flight finding that rapid intervention is essential. In a controlled setting, \cite{Meyer2019} describe the results of a controlled crossover experiment where participants were exposed to flight conditions in a hypobaric chamber. Those authors found significant, concurrent increases in heart rate and heart rate variability over the duration of the flight when compared to control. Despite these existing studies, little work has been done to examine the lagged effects of exposure to altitude on post-flight heart health.

Lagged, exposure-related effects in flight may take on different forms. One such effect is blood oxygen saturation levels will naturally fluctuate over time but are also impacted by flight conditions. Specifically, desaturation will occur when at altitude \citep{Meyer2019}. This exposure-related effect, related to being exposed to altitude, may have a deleterious impact on post-flight cardiac health as measured by heart rate and heart rate variability. In-flight blood oxygen saturation was collected in the study described by \cite{Meyer2019} every five-minutes over the duration of a six hour simulated flight, but not analyzed. Its relationship to post-flight heart and heart rate variability also not explored. A natural choice of model for analysis is the distributed lag model which relates lagged effects of some exposure variable on an outcome taken after the exposure (or at least concurrently with the last lagged effect).


Distributed lag models or DLMs are commonly used methods for analyzing lagged effects of exposure on various outcomes and a large body of literature exists regarding their formulation and estimation. We briefly review a limited selection of this literature. DLMs may be unconstrained, where each lag has its own effect in a model, or constrained, where some transformation has been performed to reduce the parameter space \citep{Zanobetti2000}. For example, with regards to the penalization of correlated lagged effects, \cite{Welty2009} consider penalized spline. \cite{Chen2018} and \cite{Chen2019} employ ridge-type shrinkage to respectively compare between unconstrained DLMs to constrained and consider interactions of multiple pollutants. \cite{Liu2018} uses a novel group and fused lasso penalty and kernel regression for multi-pollutant mixtures. \cite{Wilson2017} use principal component bases to constrain lagged effects for identifying windows of susceptibility that may interact with different covariates.

Other authors consider various types of clustered effects in DLMs. For instance, \cite{Chang2015} clusters hazards by week in a hierarchical time-to-event DLM. \cite{Liu2017} examine spatial-temporal models with lagged effects of pollution exposure. For binary outcomes, \cite{WarrenLuben2020} develop a model that estimates subject-specific main effects that vary by both spatial location and time. From the epidemiology literature, \cite{Baek2016a} explore hierarchical DLMs to analyze associations between the built environment and health with spatial effects. In the context of a case-control design, \cite{WarrenKong2020} build spatial models to identify critical windows of susceptibility to exposure. While these clustered and hierarchical models often incorporate random intercepts (or their equivalence), none examine the case where lags are allowed to randomly vary by subject. However, in the case of chamber study, the lagged in-flight exposures may vary by subject.

In this manuscript, we propose a unified modeling framework for the use of longitudinal distributed lag models or LDLMs. Our models can accommodate subject-specific, or random, effects that can potentially arise from more complex structures when lagged exposures are longitudinally sampled, including under the crossover experimental design. We employ penalized splines to estimate lagged effects and to model random lags. As the models may be computationally intensive, we conduct estimation using variational Bayesian inference. This necessitates the development of model selection tools for variational LDLMs as it is common in longitudinal analyses to select the random effect structure before considering the mean model \citep{Fitz2011}. Thus, we derive variational AICs (VAICs) for our LDLMs and determine decision rules for model selection. We show in simulation that the VAICs performs well in identifying the correctly specified model and that our models have good properties in terms of bias and mean integrated squared error (MISE). We also investigate two interval construction procedures, point-wise intervals and simultaneous bands, for the crossover design case. Finally, we show that the variational algorithm admits a form of the estimated mean that is well-suited for use in a Frequentist test for the equivalence of two curves, applied to the difference resulting from crossover design.

The motivation for this work comes from \cite{Meyer2019}, a hypobaric chamber study examining the impact of exposure to altitude on cardiac health in a crossover design. While the concurrent, in-flight effects have previously be studied \citep{Meyer2019}, the post-flight effects of in-flight exposures have not. To our knowledge, no previous studies have assessed the impact of lagged blood oxygen saturation/desaturation on post-flight heart health. We investigate these effects using the full LDLM framework for analysis, conducting model selection and inference for two outcomes of interest: post-flight heart rate and post-flight root Mean Squared Standard Deviations of N-N intervals (rMSSD), a measure of heart rate variability. We examine the impact of exposure to flight both in the full data and in a subgroup analysis of the pre-specified groups.

The remainder of this manuscript is organized as follows: Section~\ref{s:data} describes out motivating data in detail. In Section~\ref{s:stat}, we establish our framework for estimating LDLMs, conducting model selection, and performing inference. Section~\ref{s:sim} details our simulation study where we assess the discrimination of the VAIC and Frequentist characteristics of our inference procedures. In Section~\ref{s:chamber}, we conduct our analysis of the chamber study data examining the impact of lagged in-flight blood oxygen saturation on post-flight heart rate and rMSSD. Finally, we provide a discussion of our work in Section~\ref{s:disc}.

\section{Motivating Data}
\label{s:data}

\cite{Meyer2019} present a crossover design study examining the impact of exposure to altitude in-flight on concurrent heart metrics. This study serves as the motivation for our LDLM. As part of the study, participants entered a hypobaric pressure chamber outfitted with airplane seats on two separate days. On one day, the chamber's air pressure was depressurized to the equivalent of 7,000 ft altitude to simulate a commercial flight. This flight day can be considered the active treatment. On the other day, the chamber was turned on but the air pressure was not depressurized. This day is the control treatment. The order in which participants received the treatment, flight or control, was block randomized.

During the days, participant blood oxygen saturation was monitored every five-minutes for the duration of the simulated flight (or control) for approximately six hours. Post-flight measurements of heart health were also monitored including heart rate and heart rate variability as measured by rMSSD. In the original study, \cite{Meyer2019} found the greatest impact of concurrent exposure to altitude on heart rate and rMSSD. We aim to see if these effects persist in the first five minute measurement taken immediately after flight (or control). While \cite{Meyer2019} use flight day and control day as their proxy for the biological effects of exposure to altitude, we use the five-minute blood oxygen saturation as our lagged exposure. This exposure will vary by day, flight vs control, with each participant serving as their own control. Additional variables of interest to the analysis including an experiment effect to control for the order in which the treatment occurred, the time of day the participant entered the chamber which can impact autonomic responses, and...

In total, there were 26 participants available for analysis. Three subgroups were identified a priori, related to the relative health of the participants. All participants were over the age of 55 and categorized as being healthy (11 participants), smokers (6), or those with diagnosed cardiovascular disease (9). These three subgroups are the focus of a subgroup analysis to see if the impact of lagged blood oxygen saturation varies by group. All participants provided informed written consent. Protocols for the study were approved by the Human Subjects Committees at the Harvard School of Public Health (Institutional Review Board Protocol \#P15170-101), the Civilian Aerospace Medical Institute, and the University of Oklahoma Medical Center.

\section{Statistical Methods}
\label{s:stat}

\subsection{Variational Bayesian Inference}


	Let ${\mbox{\boldmath$\theta$}}$ be a vector of parameters and $\textbf{Y}$ denote a vector of observed data. The posterior distribution for ${\mbox{\boldmath$\theta$}}$ given $\textbf{Y}$ is $p({\mbox{\boldmath$\theta$}} | \textbf{Y}) = p(\textbf{Y}, {\mbox{\boldmath$\theta$}})\big/p(\textbf{Y})$. Typically, $p(\textbf{Y})$ is intractable and posterior estimates need to be obtained algorithmically. For the arbitrary density $q$, $p(\textbf{Y})$ is bounded below by $\underline{p}(\textbf{Y}; {\mbox{\boldmath$\theta$}})$ where
	\begin{align*}
		\underline{p}(\textbf{Y}; q) = \exp\left[ \int q({\mbox{\boldmath$\theta$}}) \log\left\{ \frac{p(\textbf{Y}, {\mbox{\boldmath$\theta$}})}{q({\mbox{\boldmath$\theta$}})} \right\} d{\mbox{\boldmath$\theta$}} \right].
	\end{align*}
	Variational algorithms work by maximizing $\underline{p}(\textbf{Y}; q)$ over a class of densities, $q$, that are tractable. This in turn minimizes the Kullback-Leibler Divergence between $q({\mbox{\boldmath$\theta$}})$, and the posterior itself, $p({\mbox{\boldmath$\theta$}} | \textbf{Y})$.
	
	Given a partition of ${\mbox{\boldmath$\theta$}}$ in to $M$ subcomponents, $\{{\mbox{\boldmath$\theta$}}_1, \ldots, {\mbox{\boldmath$\theta$}}_M\}$, we use the mean field variational Bayesian approximation to construct $q$, specifically $q({\mbox{\boldmath$\theta$}}) = \prod_{m = 1}^M q_m({\mbox{\boldmath$\theta$}}_m)$. This approach approximates the posterior as the product of these $q$-densities which are analogous to conditional posterior densities obtained for a Gibbs sampler \citep{Gelman2013}. Optimal densities are obtained iteratively with convergence, and therefore minimization, occurring when changes in the variation lower bound, $\underline{p}(\textbf{y}; q)$, become negligible. For a thorough introduction to variational Bayesian techniques with examples, see \cite{OrmerodWand2010} and \cite{Blei2017}. 	

\subsection{Variational Longitudinal DLMs}

Let $Y_{ij}$ denote the $j$th response of the $i$th subject. We define the following models:
	\begin{align}
		Y_{ij} &= \textbf{x}_{ij}'{\mbox{\boldmath$\beta$}} + \textbf{l}_{ij}'{\mbox{\boldmath$\gamma$}}_j +  \textbf{l}_{ij}'\textbf{g}_i + \varepsilon_{ij} \text{ (random lag) and} \label{eq:dlmrl}\\
		Y_{ij} &= \textbf{x}_{ij}'{\mbox{\boldmath$\beta$}} + \textbf{l}_{ij}'{\mbox{\boldmath$\gamma$}}_j + b_i  + \varepsilon_{ij}\text{ (random intercept)}\label{eq:dlmri}
	\end{align}
	where $\textbf{x}_{ij}$ and $\textbf{l}_{ij}$ are $1\times P$ and $1\times\ell$ vectors, respectively, of fixed and lagged effects with corresponding vectors of coefficients ${\mbox{\boldmath$\beta$}}$ and ${\mbox{\boldmath$\gamma$}}_j$. We allow the lagged effects to vary by $j$ for the case where each subject receives both treatments under a crossover design. Thus, each treatment has its own estimated effect. The target of inference under the crossover design is the difference curve, ${\mbox{\boldmath$\delta$}} = {\mbox{\boldmath$\gamma$}}_1 - {\mbox{\boldmath$\gamma$}}_0$ where $j = 0$ could denote the lagged effects under a control and $j = 1$ could denote an active treatment. For the longitudinal case, we drop the $j$ index and ${\mbox{\boldmath$\gamma$}}_j = {\mbox{\boldmath$\gamma$}}$. We assume $\varepsilon_{ij} \stackrel{iid}{\sim} N(0, \sigma^2)$ for both models.
	
	To reduce the number of parameters that need to be estimated and to model the functional form of the distributed lag, we perform a basis expansion of the lagged effects. Suppose ${\mbox{\boldmath$\Theta$}}$ is the $\ell\times K$ matrix of known basis functions. The basis expanded coefficients are ${\mbox{\boldmath$\gamma$}}_j = {\mbox{\boldmath$\Theta$}}{\mbox{\boldmath$\gamma$}}_j^S$ (for the crossover model), ${\mbox{\boldmath$\gamma$}} = {\mbox{\boldmath$\Theta$}}{\mbox{\boldmath$\gamma$}}^S$ (longitudinal model), and $\textbf{g}_i = {\mbox{\boldmath$\Theta$}}\textbf{g}_i^S$ for some selection of knots, $K$---the basis expansions do not necessarily have to have the same number of knots. Under the crossover design, the spline-based versions of equations~\eqref{eq:dlmrl} and~\eqref{eq:dlmri} are
	\begin{align}
		Y_{ij} &= \textbf{x}_{ij}'{\mbox{\boldmath$\beta$}} + \textbf{l}_{ij}'{\mbox{\boldmath$\Theta$}}{\mbox{\boldmath$\gamma$}}_j^S +  \textbf{l}_{ij}'{\mbox{\boldmath$\Theta$}}\textbf{g}_i^S + \varepsilon_{ij} \text{ (random lag) and} \label{eq:sprl}\\
		Y_{ij} &= \textbf{x}_{ij}'{\mbox{\boldmath$\beta$}} + \textbf{l}_{ij}'{\mbox{\boldmath$\Theta$}}{\mbox{\boldmath$\gamma$}}_j^S + b_i + \varepsilon_{ij}\text{ (random intercept)}. \label{eq:spri}
	\end{align}
	When the data arises from a longitudinal design, ${\mbox{\boldmath$\gamma$}}_j^S = {\mbox{\boldmath$\gamma$}}^S$.
	
		For both the fixed and random lag components, we implement penalized B-splines or P-splines using the mixed model representation. Thus, ${\mbox{\boldmath$\Theta$}}$ is a matrix of B-spline basis functions of size dependent upon the number of knots selected and the number of lags. The corresponding penalty matrix, $\mathcal{P} = \xi D_0 + (1-\xi)D_2$, is weighted between the zeroth derivative matrix ($D_0$) and the second derivative matrix ($D_2$). The parameter $\xi$ controls the desired tradeoff between shrinkage and smoothness, with values near 0 favoring shrinkage. Depending on the model, the penalty priors on ${\mbox{\boldmath$\gamma$}}_j^S$ and ${\mbox{\boldmath$\gamma$}}^S$ are ${\mbox{\boldmath$\gamma$}}_j^S \sim N(\textbf{0}, \lambda_{j} \mathcal{P}_j^{-1})$ and ${\mbox{\boldmath$\gamma$}}^S \sim N(\textbf{0}, \lambda_{\gamma} \mathcal{P}_{\gamma}^{-1})$, where $\lambda_{j}$ and $\lambda_{\gamma}$ are tuning-parameters and the penalty matrix is allowed to vary in size under the crossover design. We select a small number of knots, setting $K_j = 8$ for $j = 0, 1$ or $K_{\gamma} = 8$ under either design.
				
		The prior on the subject-specific lags, $\textbf{g}_i^S$, is $\textbf{g}_i^S \stackrel{iid}{\sim} N\left[\textbf{0}, \lambda_{g} \mathcal{P}_g^{-1}\right]$ where $\lambda_{g}$ is a tuning-parameter. The penalty matrix, $\mathcal{P}_g$, has the same general form as $\mathcal{P}$ with dimension varying on the choice of $K_g$. We take the number of knots for the random lagged effect to be $K_g = n$, where $n$ is the total number of subjects. This choice corresponds to allowing one knot for each subject which effectively reduces the random lag dimension, in the transformed space, to a single subject-specific parameter. Subject-specific lags are still estimated after transformation but as our model generates estimates in the transformed space, this choice gives us a large dimension reduction while maintaining separate parameters for each subject.

Regardless of model, we place weakly informative priors on the components ${\mbox{\boldmath$\beta$}}$, $\beta_p \stackrel{iid}{\sim} N(0, \sigma_b^2)$ with $\sigma_b^2$ fixed at something large. For the random intercept model, the vector $\textbf{b}$ consists of the $b_i$ which are independent normals, $b_i \stackrel{iid}{\sim} N(0, \sigma_u^2)$. We place inverse-gamma priors (and hyper-priors) on all variance components: $\sigma^2 \sim IG(a_e, b_e)$,  $\lambda_{j} \sim IG(a_{j}, b_{j})$ or $\lambda_{\gamma} \sim IG(a_{\gamma}, b_{\gamma})$, and $\sigma_u^2 \sim IG(a_u, b_u)$ or $\lambda_{g} \sim IG(a_{g}, b_{g})$, depending on the model. The hyper-parameters for each variance prior are set to something small.

Using a mean field variational Bayesian approximation, we obtain approximation densities or $q$-densities for each component of the spline-based models in equations~\eqref{eq:sprl} and~\eqref{eq:spri}. We let ${\mbox{\boldmath$\theta$}}$ generically denote the coefficients which vary by model and design. The random lag model has ${\mbox{\boldmath$\theta$}} = [\begin{array}{ccc} {\mbox{\boldmath$\beta$}} & {\mbox{\boldmath$\gamma$}}^S &  \textbf{g}^S \end{array}]$ (longitudinal) or ${\mbox{\boldmath$\theta$}} = [\begin{array}{cccc} {\mbox{\boldmath$\beta$}} & {\mbox{\boldmath$\gamma$}}_0^S & {\mbox{\boldmath$\gamma$}}_1^S & \textbf{g}^S \end{array}]$ (crossover) for all transformed subject-specific lag effects $\textbf{g}^S$. In general, the $q$-densities are then
\begin{align*}
	q({\mbox{\boldmath$\theta$}}) &\sim N\left[{\mbox{\boldmath$\mu$}}_{q({\mbox{\boldmath$\theta$}})}, {\mbox{\boldmath$\Sigma$}}_{q({\mbox{\boldmath$\theta$}})}\right],\\
	q(\sigma^2) &\sim IG\left[ a_e + \frac{N}{2}, B_{q(\sigma^2)} \right],\\
	q(\lambda_j) &\sim IG\left[ a_j + \frac{K_j}{2}, B_{q(\lambda_j)}\right], \text{ and}\\
	q(\lambda_g) &\sim IG\left[ a_g + \frac{K_g}{2}, B_{q(\lambda_g)}\right],
\end{align*}
	where the $j$ can take on 0 or 1 (indicating crossover treatments) or $\gamma$ (longitudinal). The subscript $q(\cdot)$ notation indicates the parameter to which the quantity belongs under the mean field approximation. For the random intercept model, $\textbf{g}^S$ is replaced with $\textbf{b}$, $a_g$ with $a_u$, $b_g$ with $b_u$, $K_g$ with $K_u$, and $\lambda_g$ with $\sigma^2_u$. The specific values for the quantities ${\mbox{\boldmath$\mu$}}_{q({\mbox{\boldmath$\theta$}})}$,  ${\mbox{\boldmath$\Sigma$}}_{q({\mbox{\boldmath$\theta$}})}$, $B_{q(\sigma^2)}$, $B_{q(\lambda_0)}$, $B_{q(\lambda_1)}$, and $B_{q(\lambda_g)}$ are in Table~\ref{t:algo} along with the algorithm for the crossover design for the design matrix $\textbf{C}$ which contains matrix versions of the data vectors from equations~\eqref{eq:sprl} and~\eqref{eq:spri}.
	
	Let $\textbf{Y}$ be a $\left(\sum_{i=1}^n m_i\right)\times 1$ vector where $m_i$ denotes the number of measurements per subject $i$, $\textbf{X}$ be a $\left(\sum_{i=1}^n m_i\right)\times P$ matrix for $P$ fixed effects, and $\textbf{U}$ be a $\sum_{i=1}^n m_i \times n$ block diagonal matrix. For the main distributed lag effects, $\textbf{L}$ is a block diagonal $\left(\sum_{i=1}^n m_i\right)\times 2\ell$ matrix under the crossover design and $\left(\sum_{i=1}^n m_i\right)\times \ell$ under the longitudinal design. The subject-specific distributed lag design matrix $\textbf{W} = diag\{\textbf{L}_i\}$ is a block diagonal matrix with blocks consisting of the subject-specific matrices $\textbf{L}_i$ formed by stacking the $\textbf{l}_{ij}$. Thus, $\textbf{W}$ has dimension $\left(\sum_{i=1}^n m_i\right)\times \ell n$.  In the case of the crossover design, as in the chamber study, $m_i = 2\ \forall\ i$ and the full design matrix is $\textbf{C} = [\textbf{X}\ \textbf{L}(\textbf{I}_2\otimes\Theta)\ \textbf{W}(\textbf{I}_n \otimes\Theta)]$ for the random distributed lag and $\textbf{C} = [\textbf{X}\ \textbf{L}(\textbf{I}_2\otimes{\mbox{\boldmath$\Theta$}})\ \textbf{U}]$ for the random intercept. The design matrices simplify under the longitudinal design to $\textbf{C} = [\textbf{X}\ \textbf{L}{\mbox{\boldmath$\Theta$}}\ \textbf{W}(\textbf{I}_n \otimes{\mbox{\boldmath$\Theta$}})]$ (random lag) or $\textbf{C} = [\textbf{X}\ \textbf{L}{\mbox{\boldmath$\Theta$}}\ \textbf{U}]$ (random intercept). The variational lower bound for the crossover design random lag model is
	\begin{align*}
		\log[ \underline{p}(\textbf{y}; q) ] &= \frac{1}{2}(P + K_0 + K_1 + K_g) - \frac{N}{2}\log(2\pi) - \frac{P}{2}\log(\sigma_b^2) \\
			&+ \frac{1}{2}\log\left(|{\mbox{\boldmath$\Sigma$}}_{q({\mbox{\boldmath$\theta$}})}|\right) - \frac{1}{2\sigma_b^2}\left[ {{\mbox{\boldmath$\mu$}}_{q({\mbox{\boldmath$\beta$}})}}'{\mbox{\boldmath$\mu$}}_{q({\mbox{\boldmath$\beta$}})} + \text{tr}\left\{{\mbox{\boldmath$\Sigma$}}_{q({\mbox{\boldmath$\beta$}})}\right\} \right] \\
			&- a_e\log(b_e) - \left(a_e + \frac{N}{2}\right)\log(B_{q(\sigma^2)}) + \log\left(\Gamma\left(a_e + \frac{N}{2}\right)\right) - \log\left(\Gamma(a_e)\right) \\
			&+ a_0\log(b_0) - \left(a_0 + \frac{K_0}{2}\right)\log(B_{q(\lambda_0)}) + \log\left(\Gamma\left(a_0 + \frac{K_0}{2}\right)\right) - \log(\Gamma(a_0))	\\
			&+ a_1\log(b_1) - \left(a_1 + \frac{K_1}{2}\right)\log(B_{q(\lambda_1)}) + \log\left(\Gamma\left(a_1 + \frac{K_1}{2}\right)\right) - \log(\Gamma(a_1))	\\
			&+ a_g\log(b_g) - \left(a_g + \frac{K_g}{2}\right)\log(B_{q(\lambda_g)}) + \log\left(\Gamma\left(a_g + \frac{K_g}{2}\right)\right) - \log(\Gamma(a_g)),
	\end{align*}
	where ${\mbox{\boldmath$\mu$}}_{q({\mbox{\boldmath$\beta$}})}$ is ${\mbox{\boldmath$\mu$}}_{q({\mbox{\boldmath$\theta$}})}$ subset to the quantities corresponding to ${\mbox{\boldmath$\beta$}}$, likewise for ${\mbox{\boldmath$\Sigma$}}_{q({\mbox{\boldmath$\beta$}})}$. For the random intercept model, the last line replaces with $a_g$ with $a_u$, $b_g$ with $b_u$, $K_g$ with $K_u$, and $\lambda_g$ with $\sigma^2_u$.

Code to implement all algorithms in \texttt{R} is available at \url{https://github.com/markjmeyer/LDLM}.

	\begin{table}[b]
		\centering
		\caption{Algorithm for the random lag model under the crossover design. The algorithm for the random intercept model replaces $B_{q(\lambda_g)}$ with $B_{q(\sigma_u^2)}$ and sets $\mathcal{P}_g = I_{N/2\times N/2}$. \label{t:algo}}
		\begin{tabular}{lll}
			\hline
			\hline
			Set: & $B_{q(\sigma^2)}, B_{q(\lambda_0)}, B_{q(\lambda_1)}, B_{q(\lambda_{\gamma})}, B_{q(\lambda_g)} > 0$ & \\
			Iterate: &  & \\
			 & Crossover: \\
			 & \footnotesize \ \ \ \ \  ${\mbox{\boldmath$\Sigma$}}_{q({\mbox{\boldmath$\theta$}})} = \left[ \frac{a_e + \frac{N}{2}}{B_{q(\sigma^2)}}\textbf{C}'\textbf{C} + \text{blockdiag}\left\{ (\sigma_b^2)^{-1} I_{p\times p}, \frac{a_0 + \frac{1}{2} K_0}{B_{q(\lambda_0)}} \mathcal{P}_0, \frac{a_1 + \frac{1}{2} K_1}{B_{q(\lambda_1)}} \mathcal{P}_1, \frac{a_g + \frac{1}{2} K_g}{B_{q(\lambda_g)}} \mathcal{P}_g  \right\}\right]^{-1}$ & \\
			 & Longitudinal: \\
			 & \footnotesize \ \ \ \ \  ${\mbox{\boldmath$\Sigma$}}_{q({\mbox{\boldmath$\theta$}})} = \left[ \frac{a_e + \frac{N}{2}}{B_{q(\sigma^2)}}\textbf{C}'\textbf{C} + \text{blockdiag}\left\{ (\sigma_b^2)^{-1} I_{p\times p}, \frac{a_{\gamma} + \frac{1}{2} K_{\gamma}}{B_{q(\lambda_{\gamma})}} \mathcal{P}_{\gamma}, \frac{a_g + \frac{1}{2} K_g}{B_{q(\lambda_g)}} \mathcal{P}_g  \right\}\right]^{-1}$ & \\
			 &  & \\
			 & \footnotesize ${\mbox{\boldmath$\mu$}}_{q({\mbox{\boldmath$\theta$}})} = \left( \frac{a_e + \frac{N}{2}}{B_{q(\sigma^2)}}\right){\mbox{\boldmath$\Sigma$}}_{q{\mbox{\boldmath$\theta$}})}\textbf{C}'\textbf{Y} $ & \\
			 &  & \\
			 & \footnotesize $B_{q(\sigma^2)} = b_e + \frac{1}{2}\left[ \left\{\textbf{Y} - \textbf{C}{\mbox{\boldmath$\mu$}}_{q({\mbox{\boldmath$\theta$}})} \right\}'\left\{\textbf{Y} - \textbf{C}{\mbox{\boldmath$\mu$}}_{q({\mbox{\boldmath$\theta$}})} \right\} + \text{tr}\left\{ \textbf{C}'\textbf{C} \right\}{\mbox{\boldmath$\Sigma$}}_{q({\mbox{\boldmath$\theta$}})} \right]$ & \\
			 &  & \\
			 & Crossover: \\
			 & \footnotesize \ \ \ \ \   $B_{q(\lambda_0)} = b_ 0 + \frac{1}{2}\left[ {{\mbox{\boldmath$\mu$}}_{q({\mbox{\boldmath$\gamma$}}_0^S)}}'{\mbox{\boldmath$\mu$}}_{q({\mbox{\boldmath$\gamma$}}_0^S)} + \text{tr}\left\{\frac{a_0 + \frac{1}{2} K_0}{B_{q(\sigma_0^2)}} \mathcal{P}_0 \right\} \right]$ & \\
			 & \footnotesize \ \ \ \ \   $B_{q(\lambda_1)} = b_ 1 + \frac{1}{2}\left[ {{\mbox{\boldmath$\mu$}}_{q({\mbox{\boldmath$\gamma$}}_1^S)}}'{\mbox{\boldmath$\mu$}}_{q({\mbox{\boldmath$\gamma$}}_1^S)} + \text{tr}\left\{\frac{a_1 + \frac{1}{2} K_1}{B_{q(\sigma_1^2)}} \mathcal{P}_1 \right\} \right]$ & \\
			 & Longitudinal: \\
			 & \footnotesize \ \ \ \ \   $B_{q(\lambda_{\gamma})} = b_ {\gamma} + \frac{1}{2}\left[ {{\mbox{\boldmath$\mu$}}_{q({\mbox{\boldmath$\gamma$}}_{\gamma}^S)}}'{\mbox{\boldmath$\mu$}}_{q({\mbox{\boldmath$\gamma$}}_{\gamma}^S)} + \text{tr}\left\{\frac{a_{\gamma} + \frac{1}{2} K_{\gamma}}{B_{q(\sigma_{\gamma}^2)}} \mathcal{P}_{\gamma} \right\} \right]$ & \\
			 &  & \\
			 & \footnotesize $B_{q(\lambda_g)} = b_ g + \frac{1}{2}\left[ {{\mbox{\boldmath$\mu$}}_{q({\mbox{\boldmath$\gamma$}}_g^S)}}'{\mbox{\boldmath$\mu$}}_{q({\mbox{\boldmath$\gamma$}}_g^S)} + \text{tr}\left\{\frac{a_g + \frac{1}{2} K_g}{B_{q(\sigma_g^2)}}  \mathcal{P}_g \right\} \right]$ & \\
			 &  & \\
			\multicolumn{2}{l}{until changes in $\log[ \underline{p}(\textbf{y}; q) ]$ are negligible.}  &\\
			\hline
			\hline
		\end{tabular}
	\end{table}


\subsection{Model Selection Criterion}


 \cite{You2014} propose a variational AIC (VAIC) for multiple linear regression. In that context, they show that the VAIC has good asymptotic properties and converges in probability to the standard AIC. The general form of their VAIC is
	\begin{align*}
		VAIC \equiv - 2\log{p(\textbf{Y} | {\mbox{\boldmath$\theta$}}^*)} + 2P^*_D
	\end{align*}
where ${\mbox{\boldmath$\theta$}}^* = E_q({\mbox{\boldmath$\theta$}})$ and $P^*_D = 2\log{p(\textbf{Y} | {\mbox{\boldmath$\theta$}}^*)} - 2 E_q[\log{p(\textbf{Y}|{\mbox{\boldmath$\theta$}})}]$. The expectation is taken with respect to the variational approximation densities, i.e. the $q$-densities. We derive $\log{p(\textbf{Y} | {\mbox{\boldmath$\theta$}}^*)}$ and $E_q[\log{p(\textbf{Y}|{\mbox{\boldmath$\theta$}})}]$ for the LDLM to obtain
	\begin{align*}
		\log{p(\textbf{Y} | {\mbox{\boldmath$\theta$}}^*)} = &-\frac{N}{2} \log(2\pi) -\frac{N}{2}\left[ \log\left\{B_{q(\sigma^2)}\right\} - \log\left(a_e +\frac{N}{2} - 1\right) \right] \\
		&- \frac{1}{2}\frac{a_e +\frac{N}{2} - 1}{ B_{q(\sigma^2)}}\left\{\textbf{Y} - \textbf{C}{\mbox{\boldmath$\mu$}}_{q({\mbox{\boldmath$\theta$}})}\right\}'\left\{\textbf{Y} - \textbf{C}{\mbox{\boldmath$\mu$}}_{q({\mbox{\boldmath$\theta$}})}\right\} \text{ and}\\
		E_q[ \log{p(y | \theta)} ] = &-\frac{N}{2} \log(2\pi) + \frac{N}{2} \left[ \psi\left(a_e +\frac{N}{2}\right) - \log\left\{ B_{q(\sigma^2)} \right\}   \right] \\
		&- \frac{1}{2}  \frac{a_e + \frac{N}{2}}{ B_{q(\sigma^2)} } \left[ \text{tr}\left\{ \textbf{C}{\mbox{\boldmath$\Sigma$}}_{q({\mbox{\boldmath$\theta$}})} \textbf{C}' \right\} + \left\{\textbf{Y} - \textbf{C}{\mbox{\boldmath$\mu$}}_{q({\mbox{\boldmath$\theta$}})} \right\}'\left\{ \textbf{Y} - \textbf{C}{\mbox{\boldmath$\mu$}}_{q({\mbox{\boldmath$\theta$}})} \right\} \right],
	\end{align*}
	where $\psi(\cdot)$ denotes the digamma function and ${\mbox{\boldmath$\theta$}}$ depends on model specification, as above. Combining, the VAIC for the LDLM is then
	\begin{align*}
		VAIC &=  N\log\left(a_e +\frac{N}{2} - 1\right) + N \log(2\pi) - 2N \psi\left(a_e +\frac{N}{2}\right)  + N\log\left\{ B_{q(\sigma^2)} \right\}  \\
		&+ 2  \frac{a_e + \frac{N}{2}}{ B_{q(\sigma^2)} } \text{tr}\left\{ \textbf{C}{\mbox{\boldmath$\Sigma$}}_{q({\mbox{\boldmath$\theta$}})} \textbf{C}' \right\} +   \frac{a_e + \frac{N}{2} + 1}{ B_{q(\sigma^2)} } \left\{\textbf{Y} - \textbf{C} {\mbox{\boldmath$\mu$}}_{q({\mbox{\boldmath$\theta$}})} \right\}'\left\{\textbf{Y} - \textbf{C} {\mbox{\boldmath$\mu$}}_{q({\mbox{\boldmath$\theta$}})} \right\}.
	\end{align*}
	This formulation applies to both the random intercept and random lag models where ${\mbox{\boldmath$\Sigma$}}_{q({\mbox{\boldmath$\theta$}})}$ and ${\mbox{\boldmath$\mu$}}_{q({\mbox{\boldmath$\theta$}})}$ are the values from the algorithm upon convergence and will depend on the model. Because the theoretical results in \cite{You2014} were derived for non-longitudinal models, we empirically examine the properties of the VAIC we derive for use in selecting the structure of the random effects in LDLMs.

As with the standard AIC, the smallest VAIC often suggests the best model fit. We refer to this as the minimum decision rule. Parsimonious models are, however, also desirable when model fits are close. Thus other rules might impose a threshold on the absolute difference in VAICs between candidate models. When the difference is under the threshold, the more parsimonious model, i.e. the random intercept model, would be selected. We examine the minimum decision rule along with absolute difference thresholds of 2, 5 and 10 in our simulation study. 

\subsection{Interval Estimation for Crossover Design}

The tendency of variational algorithms to underestimate the variance of mean model components, even when the means themselves are unbiased and consistent, is well documented \citep{Gelman2013}. Any inferential developments must take this into consideration and thoroughly investigate the properties of intervals derived from variational Bayesian estimation. Our motivating data comes from a crossover design study of the impact of in-flight exposures on post-flight cardiac health related outcomes. Accordingly, we consider the construction of intervals for only ${\mbox{\boldmath$\delta$}} = {\mbox{\boldmath$\gamma$}}_1 - {\mbox{\boldmath$\gamma$}}_0$, the difference in treatment-specific curves.

We examine two different interval constructions for ${\mbox{\boldmath$\delta$}}$: a point-wise interval using the variational estimates and a simultaneous band, as described by \cite{Ruppert2003}. For the point-wise intervals, we construct variance for each element of ${\mbox{\boldmath$\delta$}}$, $\delta(t)$, in the usual fashion: $Var[\delta(t)] = Var[\gamma_1(t)] + Var[\gamma_0(t)] - 2Cov[\gamma_1(t), \gamma_0(t)]$ with variances and covariances estimated using the corresponding components of ${\mbox{\boldmath$\Sigma$}}_{q({\mbox{\boldmath$\theta$}})}$, upon convergence. Since the $q$-densities for both ${\mbox{\boldmath$\gamma$}}_0$ and ${\mbox{\boldmath$\gamma$}}_1$ are Gaussian, the distribution of their difference is also Gaussian. Using this, and the variances for each $\delta(t)$, we construct an interval at each lag, separately. Combining all lag-specific intervals gives a point-wise interval for ${\mbox{\boldmath$\delta$}}$.

Because the point-wise intervals rely on the variance estimates from the variational estimation, they may not achieve nominal coverage. Thus, we also construct a joint interval for ${\mbox{\boldmath$\delta$}}$ using simultaneous bands for semiparametric curves. Similar to \cite{Ruppert2003}, we note that
\begin{align}
	\left[ \begin{array}{c}
		\hat{{\mbox{\boldmath$\gamma$}}}_0^S - {\mbox{\boldmath$\gamma$}}_0^S \\
		\hat{{\mbox{\boldmath$\gamma$}}}_1^S - {\mbox{\boldmath$\gamma$}}_1^S \end{array}\right] \stackrel{\cdot}{\sim} N\left[ \textbf{0}, {\mbox{\boldmath$\Sigma$}}_{q({\mbox{\boldmath$\gamma$}}_0^S, {\mbox{\boldmath$\gamma$}}_1^S)} \right], \label{eq:gam}
\end{align}
where $\stackrel{\cdot}{\sim}$ denotes approximately follows. The $100(1-\alpha)$\% simultaneous band for ${\mbox{\boldmath$\delta$}}$ is $\hat{{\mbox{\boldmath$\delta$}}} \pm m_{1-\alpha}\widehat{\text{SD}}(\hat{{\mbox{\boldmath$\delta$}}})$, where $\widehat{\text{SD}}(\hat{{\mbox{\boldmath$\delta$}}})$ is calculated in the fashion as in the point-wise interval construction. The key difference is the critical value, $m_{1-\alpha}$, which is the $(1-\alpha)$ quantile of the random variable
\begin{align}
	\sup_t \left| \frac{\hat{\delta}(t) - \delta}{\widehat{\text{SD}}\{\hat{\delta}(t)\}} \right| \approx \max_t \left| \frac{\hat{\delta}(t) - \delta}{\widehat{\text{SD}}\{\hat{\delta}(t)\}} \right|. \label{eq:max}
\end{align}
To obtain $m_{1-\alpha}$, we generate 10,000 samples from~\eqref{eq:gam}, construct~\eqref{eq:max} for $\delta = 0$, and find the $(1-\alpha)$ quantile.

\subsection{Global Test of Differences}

\cite{Zhang2000} and \cite{Zhang2003} discuss hypothesis testing for comparing two or more semiparametric curves. The test statistic, which we refer to as the ZLS, requires the existence of a vector, $\textbf{c}(t)$, such that the estimate of a sermiparametric function, $f(t)$, can be written as $\hat{f}(t) = \textbf{c}(t)'Y$ at time $t$. We apply this test in the crossover design setting to construct a global test of the differences, ${\mbox{\boldmath$\delta$}}$, being equal to zero. Our variational algorithm in Table~\ref{t:algo} admits an explicit form of the $\textbf{c}(t)$ vector for estimating each element of ${\mbox{\boldmath$\delta$}}$, $\delta(t)$, which can be derived from the mean component, ${\mbox{\boldmath$\mu$}}_{q({\mbox{\boldmath$\theta$}})}$.

Let ${\mbox{\boldmath$\Sigma$}}_{q({\mbox{\boldmath$\gamma$}}_0^S)}$ and ${\mbox{\boldmath$\Sigma$}}_{q({\mbox{\boldmath$\gamma$}}_1^S)}$ denote the covariance matrices that correspond to the vectors of treatment-specific effects ${\mbox{\boldmath$\gamma$}}_0^S$ and ${\mbox{\boldmath$\gamma$}}_1^S$, respectively. Further, let $\textbf{C}_0$ and $\textbf{C}_1$ be matrices whose columns consist of the columns of the matrix $\textbf{C}$ that correspond to these treatment-specific effects. Denote the following vectors ${\mbox{\boldmath$\eta$}}_0(t) = \left( \frac{a_e + \frac{N}{2}}{B_{q(\sigma^2)}}\right){\mbox{\boldmath$\Sigma$}}_{q({\mbox{\boldmath$\gamma$}}_0^S)}\textbf{C}_{0, t}$ and ${\mbox{\boldmath$\eta$}}_1(t) =  \left( \frac{a_e + \frac{N}{2}}{B_{q(\sigma^2)}}\right){\mbox{\boldmath$\Sigma$}}_{q({\mbox{\boldmath$\gamma$}}_1^S)}\textbf{C}_{1, t}$, for the $t$th columns of $\textbf{C}_0$ and $\textbf{C}_1$. Then, $\textbf{c}_0(t) = {\mbox{\boldmath$\eta$}}_0(t) \Theta'$ and $\textbf{c}_1(t) = {\mbox{\boldmath$\eta$}}_1(t)\Theta'$ are vectors such that $\hat{\gamma}_0(t) = \textbf{c}_0(t)'Y$ and $\hat{\gamma}_1(t) = \textbf{c}_1(t)'Y$. This allows us to make the observation that, upon convergence of the algorithm, $\textbf{c}(t) = \textbf{c}_1(t) - \textbf{c}_0(t)$ is a vector such that $\hat{\delta}(t) = \textbf{c}(t)'Y$. This vector forms the foundation of our global test of $H_0 : \delta(t) = 0\ \forall\ t$.

The statistic for the ZLS test applied the crossover LDLM has the form
\begin{align}
	G(Y) = \int_0^{\ell} Y'\textbf{c}(t)\textbf{c}(t)'Y dt = Y'\textbf{S} Y, \label{eq:zls}
\end{align}
for $\textbf{S} = \int_0^{\ell} \textbf{c}(t)\textbf{c}(t)' dt$. As in \cite{Zhang2000} and \cite{Zhang2003}, we approximate the distribution of $G(Y)$ using \cite{Satterthwaite1943} as a scaled chi-squared. Thus, for the mean, $\textbf{E}$, and variance, $\textbf{V}$, of $Y$ under $H_0$, the mean and variance of $G(Y)$ are $e = \textbf{E}'\textbf{S}\textbf{E} + tr(\textbf{S}\textbf{V})$ and $\psi = 2tr\left\{ (\textbf{S}\textbf{V})^2 \right\} + 4\textbf{E}'\textbf{S}\textbf{V}\textbf{S}\textbf{E}$. Under the null, $\textbf{S}\textbf{E}$ is negligible, so these quantities are frequently approximated with $e \approx tr(\textbf{S}\textbf{V})$ and $\psi \approx 2tr\left\{ (\textbf{S}\textbf{V})^2 \right\}$. Setting the approximate versions of $e$ and $\psi$ equal to the mean and variance of a scaled chi-squared, $\kappa \chi^2_{\nu}$, results in the scaling factor being $\kappa = \psi/(2e)$ and the degrees of freedom being $\nu = 2e^2/\psi$. The $p$-value for the test can be approximated by finding $Pr\left[\chi^2_{\nu} > G(Y)/\kappa\right]$.


\section{Simulation Study}
\label{s:sim}

We design our simulation to evaluate several aspects of the LDLM framework with to goal of identifying decision rules of when the VAIC selects the ``correct'' model and how misspecification impacts the model's operating characteristics. Our simulation aims to address this using varying sample sizes of $n = 25, 50, 100$, a sampling density on the lag that can be $\ell = 60$ or $120$, and three true effects:  a peak effect, cyclical effect, and sigmoidal effect. For each combination of settings, we generate 500 datasets under a ``true'' random intercept model and 500 datasets under a ``true'' random lag model. We apply both LDLMs to each dataset to observe how misspecification impacts absolute bias, mean integrated squared error (MISE), interval coverage, size, and power. Descriptions of the functions for each of the true effects are in the Supplementary Material. We vary the structure of the design as well considering a crossover design and longitudinal designs with two and three measurements. Under the crossover design, we include a scaling factor to differentiate the ``true'' effects between the control and treatment conditions. This factor also allows us to examine the power of the ZLS test under the crossover design. The number of scaling factors varies by sample size with $s$ being a sequence from 1 to 5.5, by 0.1 when $n = 25$, a sequence from 1 to 4 by 0.1 when $n = 50$, and 1 to 2.5 by 0.1 when $n = 100$. Larger sample sizes require smaller scales to achieve high power, hence the dependence of the number of scaling factors on $n$. Thus, the number of datasets used to generate Tables~\ref{t:correct} to~\ref{t:cover} vary by sample size with 23,000 when $n = 25$, 15, 500 when $n = 50$, and 8,000 when $n = 100$.


\subsection{Model Selection Discrimination and Estimation}

The goal of the VAIC is to select a ``best'' model, but ``best'' model can vary depending on the model selection approach taken. We examine the accuracy of four different decision rules for selecting between the random intercept and random lag models: a minimum decision rule and absolute differences $< 2$, $<5$, and $<10$. For the latter decision rules, the random intercept model is selected when the absolute differences falls under the threshold. Table~\ref{t:correct} contains the percent of correctly specified models, by decision rule, for the peak effect under a crossover design when $\ell = 60$. Even under the strictest decision rule, the VAIC will still select the ``correct'' model at least 85\% of the time, depending on sample size, effect, and true underlying model.

When using the absolute difference based decision rules, an approach that favors parsimonious models, the accuracy when the true model is the random intercept model increases as the threshold increases. This does not simultaneously impact accuracy when the true model is the random lag model which maintains a high percent of correctly identified models regardless of the effect, sample size, or decision rule. When the threshold is set to 5, the VAIC correctly identifies the model in excess of at least 93\% of the time. 

\begin{table}[b]
\centering
\caption{Percent correctly identified, by decision rule, for the crossover design with $\ell = 60$. Int. is short for intercept.\label{t:correct}}
\begin{tabular}{llrcccc}
  \hline
\multirow{2}{*}{Effect} & \multirow{2}{*}{True Model} & \multirow{2}{*}{$n$}  & \multicolumn{4}{c}{Decision Rule}  \\
\cline{4-7}
  &  &  & $\min$ & $< 2$ & $< 5$ & $< 10$ \\ 
  \hline
Peak & Int. & 25 & 93.52\% & 95.45\% & 97.43\% & 98.99\% \\  
      &  & 50 & 94.40\% & 95.71\% & 96.59\% & 98.10\% \\
      &  & 100 & 97.60\% & 98.14\% & 98.88\% & 99.34\% \\
   \cline{3-7}
 & Lag  & 25 & 99.99\% & 99.98\% & 99.98\% & 99.98\%  \\  
     &  & 50 & 100.0\% & 100.0\% & 100.0\% & 100.0\% \\
     &  & 100 & 100.0\% & 100.0\% & 100.0\% & 100.0\% \\
   \cline{2-7}
Cyclical  &  Int.  & 25 & 86.93\% & 89.80\% & 93.48\% & 97.45\%  \\  
       &  & 50 & 90.83\% & 92.78\% & 95.14\% & 97.77\% \\
       &  & 100 & 96.41\% & 97.00\% & 97.80\% & 98.79\%  \\
   \cline{3-7}
  & Lag  & 25 & 100.0\% & 100.0\% & 100.0\% & 100.0\%   \\  
      &  & 50 & 100.0\% & 100.0\% & 100.0\% & 100.0\%  \\
      &  & 100 & 100.0\% & 100.0\% & 100.0\% & 100.0\%  \\
   \cline{2-7}
Sigmoidal & Int.  & 25 & 85.79\% & 89.60\% & 93.64\% & 97.70\% \\   
      &  & 50 & 90.81\% & 92.90\% & 95.35\% & 97.95\% \\ 
      &  & 100 & 96.92\% & 97.00\% & 97.65\% & 98.67\% \\ 
   \cline{3-7}
 & Lag  & 25 & 100.0\% & 100.0\% & 100.0\% & 100.0\%  \\  
     &  & 50 & 100.0\% &  100.0\% & 100.0\% & 100.0\%  \\
     &  & 100 & 100.0\% & 100.0\% & 100.0\% & 100.0\%  \\
   \hline
\end{tabular}
\end{table}

	Table~\ref{t:bias} presents the absolute bias and MISE for the crossover design when $\ell = 60$, under correct and incorrect model specifications. Overall, absolute bias and MISE are lower when the true model is the random intercept model and all metrics shrink toward zero as $n$ increases. When the model is correctly specified, the absolute bias and MISE are lower than when the model is incorrectly specified. However, in some cases the trade-off is minor such as when the sample size is small or when the true model is the random intercept model. For larger sample sizes when the true model is the random lag model, the bias and MISE can be much larger when the model is incorrectly specified. Similar simulation results for the percent correctly identified, absolute bias, and MISE when $\ell = 120$, as well as under the longitudinal design, are in the Supplementary Material. 

\begin{table}[b]
\centering
\caption{LDLM average absolute bias and MISE for the crossover design with $\ell = 60$. Int. is short for intercept.\label{t:bias}}
\begin{tabular}{llrccccc}
  \hline
  & & & \multicolumn{5}{c}{Specified Model} \\
\multirow{2}{*}{Effect} & \multirow{2}{*}{True Model} & \multirow{2}{*}{$n$}  & \multicolumn{2}{c}{Bias} & & \multicolumn{2}{c}{MISE} \\
\cline{4-5} \cline{7-8}
   &  &  & Lag & Int.  &  & Lag & Int. \\ 
  \hline
Peak  & Int.  & 25 & 0.1620 & 0.1435 & & 0.0643 & 0.0469 \\  
      &  & 50 & 0.0910 & 0.0865 & & 0.0155 & 0.0160 \\
      &  & 100 & 0.0564 & 0.0505 & & 0.0057 & 0.0046 \\
   \cline{3-8}
 & Lag  & 25 & 0.4913 & 0.6224 & & 0.4868 & 0.7901 \\   
     &  & 50 & 0.2149 & 0.4205 & & 0.0768 & 0.3033  \\
     &  & 100 & 0.1330 & 0.2738 & & 0.0303 & 0.1301 \\
   \cline{2-8}
Cyclical & Int.  & 25 & 0.1101 & 0.1027 &  & 0.0231 & 0.0193 \\  
      &  & 50 & 0.0653 & 0.0594 &  & 0.0073 & 0.0061 \\ 
      &  & 100 & 0.0427 & 0.0396 &  & 0.0031 & 0.0027 \\
   \cline{3-8}
 & Lag  & 25 & 0.4753 & 0.6139 &  & 0.4527 & 0.7656 \\  
     &  & 50 & 0.2034 & 0.4145 &  & 0.0685 & 0.2915 \\ 
     &  & 100 & 0.1267 & 0.2717 &  & 0.0276 & 0.1279 \\
      \cline{2-8}
Sigmoidal &  Int.  & 25 & 0.1107 & 0.1032 &  & 0.0232 & 0.0194 \\   
      &  & 50 & 0.0654 & 0.0595 &  & 0.0074 & 0.0062 \\ 
      &  & 100 & 0.0428 & 0.0396 &  & 0.0031 & 0.0027 \\ 
   \cline{3-8}
 & Lag  & 25 & 0.4741 & 0.6130 &  & 0.4503 & 0.7636 \\     
     &  & 50 & 0.2033 & 0.4140 &  & 0.0684 & 0.2908 \\
     &  & 100 & 0.1267 & 0.2715 &  & 0.0276 & 0.1277 \\
   \hline
\end{tabular}
\end{table}


\subsection{Interval Coverage and Global Test Characteristics}

Frequentist evaluation of variational results is increasingly more common \citep{You2014, Wang2019, Zhang2020}, with variational models frequently producing ideal Frequentist properties. We now examine the point-wise coverage of both the point-wise and simultaneous intervals for ${\mbox{\boldmath$\delta$}}$. To calculate coverage, we determine whether or not the ``true'' effect at each time point $t$ is contained the interval and then average over all $\ell$ lags and simulated datasets. The intervals used to generate the coverage in Table~\ref{t:cover} were all constructed at the 95\% level. The most notable result is that, regardless of sample size, true model, or effect, the simultaneous intervals always cover the truth. However, the point-wise intervals are often more efficient as they are closer to nominal for a number of combinations of settings. In general, when the true model is the random lag model and for sample sizes 50 and 100, the point-wise intervals achieve nominal to near-nominal coverage. When the true model is the random lag, the point-wise interval coverage is closer to nominal for correctly specified models and sample sizes of $n = 50$ or 100. At its worst, in smaller sample sizes the point-wise interval under-covers by as much as 9.1\%.

\begin{table}[b]
\centering
\caption{Coverage by interval type. Int. is short for intercept. All intervals are constructed at the 95\% level.\label{t:cover}}
\begin{tabular}{llrccccc}
  \hline
  & & & \multicolumn{5}{c}{Specified Model} \\
\multirow{2}{*}{Effect} & \multirow{2}{*}{True Model} & \multirow{2}{*}{$n$}  & \multicolumn{2}{c}{Point-wise} & & \multicolumn{2}{c}{Simultaneous} \\
\cline{4-5} \cline{7-8}
   &  &  & Lag & Int.  &  & Lag & Int. \\ 
  \hline
Peak  & Int.  & 25 & 0.909 & 0.905 & & 1.000 & 1.000 \\ 
      &  & 50 & 0.909 & 0.893 & & 1.000 & 1.000 \\
      &  & 100 & 0.903 & 0.902 & & 1.000 & 1.000 \\ 
   \cline{3-8}
 & Lag  & 25 & 0.859 & 0.953 & & 1.000 & 1.000 \\ 
     &  & 50 & 0.925 & 0.907 & & 1.000 & 1.000 \\ 
     &  & 100 & 0.931 & 0.902 & & 1.000 & 1.000 \\ 
   \cline{2-8}
Cyclical & Int.  & 25 & 0.936 & 0.932 & & 1.000 & 1.000 \\   
      &  & 50 & 0.939 & 0.935 & & 1.000 & 1.000 \\ 
      &  & 100 & 0.953 & 0.947 & & 1.000 & 1.000 \\
   \cline{3-8}
 & Lag  & 25 & 0.861 & 0.956 & & 1.000 & 1.000 \\   
     &  & 50 & 0.930 & 0.910 & & 1.000 & 1.000 \\ 
     &  & 100 & 0.939 & 0.903 & & 1.000 & 1.000 \\ 
      \cline{2-8}
Sigmoidal &  Int.  & 25 & 0.933 & 0.931 & & 1.000 & 1.000 \\   
      &  & 50 & 0.938 & 0.935 & & 1.000 & 1.000 \\
      &  & 100 & 0.953 & 0.947&  & 1.000 & 1.000 \\
   \cline{3-8}
 & Lag  & 25 & 0.861 & 0.955 & & 1.000 & 1.000 \\    
     &  & 50 & 0.930 & 0.910 & & 1.000 & 1.000 \\ 
     &  & 100 & 0.940 & 0.903 & & 1.000 & 1.000 \\ 
   \hline
\end{tabular}
\end{table}

Figures~\ref{f:power} and~\ref{f:size} contain the power and size, respectively, for the ZLS test under the crossover design. The rows correspond to different effects while the columns correspond to different ``true'' models. The solid green lines are the power or size when the specified model is a random lag while dashed blue lines are the power or size when the specified model is the random intercept. Both Figures present only the results from when $n = 25$. The graphs in Figure~\ref{f:power} vary by the scaling factor. Overall, as the scaling factor increases the power of the ZLS test increases as well. Power is always higher when the model is correctly specified. The test is the most powerful when the random  intercept model is the true model. We calculate power for Figure~\ref{f:power} using an $\alpha$ of 0.05.

The size of the ZLS test varies both by the level of the test what the true underlying model is, and whether or not the model was correctly specified. The graph in Figure~\ref{f:size} vary by the level, $\alpha$, and include a reference line to denote nominal (dotted black line). When the random intercept model is the true model, the test from the misspecified model is quite conservative. Even the correctly specified model has size below nominal for most $\alpha$. When the true model is the random lag model, however, the size of the test is closer to nominal particularly when the model has been correctly specified. The ZLS test is closest to nominal when $\alpha$ is below 0.04, although even at an $\alpha$ of 0.05 the size when the random lag model is the true model is between 0.035 and 0.04. Similar simulation results for coverage, size, and power when $\ell = 120$, along with size and power for lager sample sizes, are in the Supplementary Material. 


\begin{figure}[b]
	\centering
	\includegraphics[width = 6in]{power25}
	\caption{Power curves under all three effects. Rows denote true effects, columns correspond to the true model in the titles. Solid green lines indicate a specified random lag model, dashed blue lines indicate a specified random intercept model. The x-axis, s, indicates the scaling factor dictating the difference in the curves. \label{f:power}}
\end{figure}

\begin{figure}[b]
	\centering
	\includegraphics[width = 6in]{size25}
	\caption{Size curves under all three effects. Rows denote true effects, columns correspond to the true model in the titles. Solid green lines indicate a specified random lag model, dashed blue lines indicate a specified random intercept model. The x-axis indicates the $\alpha$ under which the test was conducted. For reference, we add a dotted black line to denote `nominal' under each $\alpha$. \label{f:size}}
\end{figure}

\section{Post-flight Effects of In-flight Exposure}
\label{s:chamber}

\subsection{Primary Analyses}



\subsection{Subgroup Analyses}

\section{Discussion}
\label{s:disc}


\bigskip
\begin{center}
{\large\bf SUPPLEMENTARY MATERIAL}
\end{center}

Additional results from the simulation study and analysis are in the supplementary materials. We make the \texttt{R} code used to produce the simulations and all analyses described in this paper, along with the relevant datasets, available at \url{https://github.com/markjmeyer}.


\bibliographystyle{plainnat}

