\section{Introduction}
Stroke (i.e., cerebrovascular disease) is the fourth leading cause of death in the general population and {third} --- behind heart disease and cancer --- among those aged 85 and older, with rates increasing exponentially with age \citep{deaths}. Geographic trends in stroke mortality have been studied as far back as the 1960s, when \citet{strokebelt} identified a region of the southeastern United States (US) stretching from Mississippi to North Carolina which had the highest rates of stroke mortality --- a region which would become known as the ``stroke belt.'' Later work by \citet{shiftingbelt} noticed an apparent shift in the stroke belt, observing that regions of the Mississippi River Valley appeared in the highest decile of mortality rates in the early 1990s where they had previously not. More recently, \citet{linda} have studied geographic trends in stroke hospitalizations from 1995 to 2006, noting that, among Medicare beneficiaries ages 65 and older, this shift in the stroke belt has persisted, stretching further into parts of Texas and Oklahoma.
In addition to changing geographic patterns, numerous studies have observed the overall declines in stroke mortality \citep[e.g.,][]{regards:declines,gillum:2011:stroke}.

While many of these studies have age-adjusted their data,
accounting for the variation in age distributions among counties and disparities in stroke mortality across age ranges, this
precludes inference within individual age groups.
{On the other hand, this aggregation and data standardization step also helps mitigate the issue of small population sizes and low number of stroke deaths in many US counties, an issue that can lead to unreliable rate estimates
and is only exacerbated when the data are stratified by a factor such as age group.
In this work, however, our goal is to investigate spatiotemporal trends in stroke mortality by jointly modeling data from three age-based subpopulations, permitting more reliable inference at the county level for each age group while preserving the ability to compute age-adjusted rates.}
More specifically, we look to build a complex multivariate space-time model which
borrows strength across space, time, and age group.
{We also propose a new tool for measuring declines in mortality which accounts for temporal changes in mortality rates.}

To achieve this,
a natural starting point is the family of models based on the conditionally autoregressive (CAR) model proposed by \citet{besag}.  Since its extension to the fully Bayesian setting in \citet{bym}, CAR models have sparked a wealth of research
in the disease mapping context for both spatial \citep[e.g.,][]{besag:poisson,besag:higdon} and spatiotemporal applications.
While these early examples were all based on the standard univariate CAR model, \citet{gelfand:mcar}
developed methods for general multivariate CAR (MCAR) models, inspiring novel approaches for both multiple and spatiotemporal disease mapping {\citep[e.g.,][]{JBC,QBC,hcar,m-b:2013,b-r:2015}.}
More recently, \citet{quick:waller} proposed a special case of the MCAR of \citet{gelfand:mcar} --- referred to as the multivariate space-time CAR ($\mbox{$\text{MSTCAR}$}$) model --- for the purpose of analyzing county-level heart disease mortality rates in the US over time for various race/gender groups.
For a more complete coverage of the recent advances in spatial and space-time modeling, see \citet{BCG}.

{
In any discussion of the spatiotemporal modeling literature, we would be remiss not to mention the subject of \emph{separability} --- i.e., models where the spatiotemporal covariance can be decomposed into the product of a purely spatial covariance and a purely temporal covariance.  While separable covariance structures offer computational benefits, \citet{stein2005} highlights some drawbacks of the separability assumption, and concerns over the utility of separable models have motivated the development of classes of nonseparable covariance functions in the univariate continuous space, continuous time setting \citep[e.g.,][]{cressie:huang, gneiting2002}.
In cases where both space and time are discrete, as encountered in this study,
\citet{knorr-held:2000} has discussed a variety of possible space-time interactions, but here again the focus has been on a single outcome, with extensions to the multivariate space-time setting being more recent developments.  In addition to the nonseparable $\mbox{$\text{MSTCAR}$}$ proposed by \citet{quick:waller}, \citet{jon} implement a reduced rank multivariate spatiotemporal mixed effects model which is designed to analyze high dimensional data efficiently.  That said, both of these approaches restrict their attention to the case of Gaussian data. {Other methods \citep[e.g.,][]{JBC,m-b:2013} allow for varying spatial structures by utilizing \emph{proper} MCAR models --- while these approaches are feasible when the number of spatial regions in the spatial domain
is small, the large number of US counties
prohibits the use of proper MCAR models.}
}

Here we extend the nonseparable $\mbox{$\text{MSTCAR}$}$ model
to the generalized linear model setting to analyze spatiotemporal trends in the dataset comprised of stroke mortality counts among those aged 65--74, 75--84, and 85+ described in Section~\ref{sec:data}.
Specifically, due to the rarity of stroke deaths,
the Gaussian assumption of \citet{quick:waller} may not be appropriate.
As such, in Section~\ref{sec:methods} we detail our approach for embedding the $\mbox{$\text{MSTCAR}$}$ model into a Poisson likelihood akin to \citet{bym}, {in addition to presenting the saved person-years (SPY) tool for measuring patterns in declines}.
We then analyze the stroke mortality data in Section~\ref{sec:anal}, where we {discover different spatiotemporal trends in mortality rates across age groups and observe evidence that a separable model would be inappropriate for these data.
Finally, we summarize our findings and offer some concluding remarks in Section~\ref{sec:disc}.}

\section{Data Description}\label{sec:data}
The study population for this analysis includes all US residents aged 65 or older. In order to assess differences across the high-risk age ranges, we separate our data into $N_g =3$ groups: those aged 65--74, those 75--84, and those 85$+$.
Annual counts of stroke-related deaths per county per age group were obtained from the National Vital Statistics System (NVSS) of the National Center for Health Statistics (NCHS).
Due to inconsistencies in the manner in which death records were recorded prior to 1973,
we restrict the analysis to data from 1973--2013 ($N_t = 41$ years) to ensure valid comparisons across time.
{Deaths from stroke were defined as those for which the underlying cause of death was cerebrovascular disease according to the 8th, 9th, and 10th revisions of the International Classification of Diseases
(ICD; {ICD--8: 430--438}; ICD--9: 430--438; ICD-10: I60--69).
{Based on the {comparability ratios reported by \citet{icd8:icd9} and \citet{icd9:icd10}, which indicate a high degree of similarity between the three revisions of the ICD}, we assume that this definition is consistent over the 41 year study period.}
The geographic unit used in this analysis was the county (or county equivalent).  Given changes in county definitions during the study period affecting ten counties (e.g., the merging/splitting of counties), a single set of $N_s$ = 3,099 regions (henceforth referred to as counties) from the contiguous lower 48 states (including the District of Columbia) was used for the entire study period. Annual population counts were based on the bridged-race intercensal estimates provided by {\citet{census}}.
{While the number of individuals in each age bracket has risen considerably ({by 89\%, 101\%, and 267\%, respectively}), the number of stroke-related deaths for individuals 65--84 has decreased nearly 60\% since 1973.
From a public health perspective, this reduction in stroke mortality is a great achievement.  From a statistical perspective, however,
this can lead to concerns over the reliability of county-level mortality rate estimates based on so few events. Figures illustrating the national population and death trends for each age group can be found in Web Appendix~B.
}

\begin{comment}
\begin{figure}[t]
    \centering
        \includegraphics[width=.6\textwidth]{time_trends.png}

    \caption{Stroke-related death rates over time for each of the age groups compared to the stroke-related death rate for individuals over 65. {\color{red}I feel like this figure is uninteresting(?) because all three groups are so different (and thus it's hard to really see anything).}}
    \label{fig:trend}
\end{figure}
\end{comment}

\section{Methods}\label{sec:methods}
\subsection{Statistical model}\label{sec:mstcar}
Letting $Y_i$ and $n_i$ denote the incidence of disease and the population at risk in county $i$, \citet{bym} proposed a model of the form
\begin{equation}
Y_{i} \sim \mbox{$\text{Pois}$}\left(n_i \exp\left[ {\bf x} _i \ensuremath{\boldsymbol{\beta}} + Z_i + \phi_i\right]\right), \text{for $i=1,\ldots,N_s$} \label{eq:bym}
\end{equation}
where $ {\bf x} _i$ denotes a $p$-vector of covariates with corresponding regression coefficients, $ \ensuremath{\boldsymbol{\beta}}$, $Z_{i}$ is a spatial random effect, and $\phi_i  \mbox{$\stackrel{\text{ind}}{\sim}$} N\left(0,\tau^2\right)$ is an exchangeable random effect.  In their work, \citet{bym} modeled $ {\bf Z} =\left(Z_1,\ldots,Z_{N_s}\right)'$ as arising from an intrinsic conditional autoregressive (CAR) model, which has the conditional distribution
\begin{equation}
Z_i\,\vert\,  {\bf Z} _{(i)}, \ensuremath{\sigma}^2 \sim N\left(\sum_{j=1}^{N_s} w_{ij} Z_j \slash \sum_{j=1}^{N_s} w_{ij},  \ensuremath{\sigma}^2\slash \sum_{j=1}^{N_s} w_{ij}\right) \label{eq:car}
\end{equation}
where $ {\bf Z} _{(i)}$ denotes the vector $ {\bf Z} $ with the $i$th element removed and $w_{ij}=1$ if $i$ and $j$ are neighbors {(denoted $i\sim j$)} and 0 otherwise. Recommendations for prior distributions for $ \ensuremath{\sigma}^2$ and $\tau^2$ are offered by \citet{bernardinelli} and \citet{waller:carlin}.

Extending~\eqref{eq:bym} and \eqref{eq:car} to a setting consisting of multiple spatial surfaces is straightforward.  Letting $Y_{ikt}$ denote the number of deaths in county $i$ during year $t$ for age group $k$, we model
\begin{equation}
Y_{ikt} \sim \mbox{$\text{Pois}$}\left(n_{ikt} \exp\left[ {\bf x} _{ikt} \ensuremath{\boldsymbol{\beta}}_{kt} + Z_{ikt} + \phi_{ikt}\right]\right),  \label{eq:bym_k}
\end{equation}
for $i=1,\ldots,N_s$, $k=1,\ldots,N_g$, and $t=1,\ldots, N_t$ where $\phi_{ikt}  \mbox{$\stackrel{\text{ind}}{\sim}$} N\left(0,\tau_{k}^2\right)$. To account for the multivariate spatiotemporal association in the data,
we follow the MSTCAR model of \citet{quick:waller} --- itself a special case of the $\mbox{$\text{MCAR}$}$ of \citet{gelfand:mcar} --- and let {$ {\bf Z}  = \left( {\bf Z} _{1\cdot\cdot}',\ldots, {\bf Z} _{N_s\cdot\cdot}'\right)' \sim\mbox{$\text{MCAR}$}\left(1, \mbox{$\Sigma$} _{\eta}\right)$},
\begin{align*}
\pi\left( {\bf Z} \,\vert\,  \mbox{$\Sigma$} _{\eta}\right) &\propto \vert  \mbox{$\Sigma$} _{\eta}\vert^{-(N_s-1)\slash 2} \exp\left[-\frac{1}{2} {\bf Z} '\left\{(D-W)\otimes \mbox{$\Sigma$} _{\eta}^{-1}\right\} {\bf Z} \right]\\
\text{and}\;\; {\bf Z} _{i\cdot\cdot}\,\vert\,  {\bf Z} _{(i)\cdot\cdot}, \mbox{$\Sigma$} _{\eta} &\sim N\left(\sum_{j\sim i}  {\bf Z} _{j\cdot\cdot} \slash m_i, \frac{1}{m_i} \mbox{$\Sigma$} _{\eta}\right),
\end{align*}
where $ {\bf Z} _{i\cdot\cdot} = \left( {\bf Z} _{i1\cdot},\ldots, {\bf Z} _{iN_g\cdot}\right)'$, $ {\bf Z} _{ik\cdot} = \left(Z_{ik1},\ldots,Z_{ikN_t}\right)'$, $W$ is an adjacency matrix with elements $w_{ij}$, $D$ is an $N_s\times N_s$ diagonal matrix with elements $m_i = \sum_{j=1}^{N_s} w_{ij}$, $ \mbox{$\Sigma$} _{\eta}$ denotes the $N_tN_g \times N_tN_g$ covariance structure for our $N_t$ years and $N_g$ age groups and $\otimes$ denotes the Kronecker product.
The $\eta$ subscript is a reference to the construction of $ {\bf Z} $ from \citet{quick:waller} where the authors began by defining $ {\bf v} _{\iota\cdot t} \mbox{$\stackrel{\text{iid}}{\sim}$} N(\textbf{0},G_t)$, where $\iota = 1,\ldots, N_s-1$, to be a collection of independent $N_g$-dimensional random variables with covariance $G_t$ for $\iota=1,\ldots,(N_s-1)$ and $t=1,\ldots, N_t$.  Using $R_k = R\left(\rho_k\right)$ to denote an age group-specific temporal correlation matrix based on an autoregressive order 1 (denoted AR(1)) model and letting $\widetilde{R}_k$ to be the Cholesky decomposition of $R_k$ such that $\widetilde{R}_k \widetilde{R}_k' = R_k$, the authors define $ \mbox{\boldmath $ \eta $} _{\iota k\cdot} = \widetilde{R}_k  {\bf v} _{\iota k\cdot}$ where $ {\bf v} _{\iota k\cdot} = \left(v_{\iota k1},\ldots,v_{\iota kN_t}\right)'$.  We then find that $ \mbox{\boldmath $ \eta $} _{\iota\cdot\cdot} \sim N\left(\textbf{0}, \mbox{$\Sigma$} _{\eta}\right)$, where
\begin{align}\label{eq:Sig_eta}
 \mbox{$\Sigma$} _{\eta} =
\begin{bmatrix}
\widetilde{R}_{1,1}^* & \textbf{0} & \textbf{0}\\
\vdots & \ddots & \textbf{0} \\
\widetilde{R}_{N_t,1}^* & \cdots & \widetilde{R}_{N_t,N_t}^*
\end{bmatrix}
\begin{bmatrix}
G_1 & \textbf{0} & \textbf{0}\\
\textbf{0} & \ddots & \textbf{0} \\
\textbf{0} & \textbf{0} & G_{N_t}
\end{bmatrix}
\begin{bmatrix}
\widetilde{R}_{1,1}^* & \cdots & \widetilde{R}_{N_t,1}^*\\
\textbf{0} & \ddots & \vdots \\
\textbf{0} & \textbf{0} & \widetilde{R}_{N_t,N_t}^*
\end{bmatrix}
\end{align}
and $\widetilde{R}_{t,t'}^*$ denotes the $N_g\times N_g$ diagonal matrix with elements $\left\{\widetilde{R}_k\right\}_{t,t'}$ for $k=1,\ldots,N_g$; $ {\bf Z} $ is then constructed from the $ \mbox{\boldmath $ \eta $} _{\iota\cdot\cdot}$ using the eigenvalues and eigenvectors of the matrix $D-W$ \citep[see][]{rue:held}. This structure is then denoted $ {\bf Z}  \sim \mbox{$\text{MSTCAR}$}\left(G_{1},\ldots,G_{N_t}, \mbox{\boldmath $\rho$} \right)$.

\subsection{Hierarchical model and computational details}
While a Poisson model like~\eqref{eq:bym_k} is a straightforward extension of~\eqref{eq:bym}, such models can also pose computational challenges, particularly for large dimensions.  For instance, the full conditional of $ {\bf Z} _{i\cdot\cdot}$, given by
\begin{align}
\pi\left( {\bf Z} _{i\cdot\cdot}\,\vert\,  {\bf Y} , {\bf Z} _{(i)\cdot\cdot}, \ensuremath{\boldsymbol{\beta}}, \mbox{\boldmath $\phi$}, \mbox{$\Sigma$} _{\eta}\right) \propto& \prod_{k=1}^{N_g} \prod_{t=1}^{N_t} \mbox{$\text{Pois}$}\left(Y_{ikt}\,\vert\, n_{ikt} \exp\left[ {\bf x} _{ikt} \ensuremath{\boldsymbol{\beta}}_{kt} + Z_{ikt} + \phi_{ikt}\right]\right) \notag\\
&\times \pi\left( {\bf Z} _{i\cdot\cdot}\,\vert\,  {\bf Z} _{(i)\cdot\cdot}, \mbox{$\Sigma$} _{\eta}\right)
\end{align}
is \emph{not} a known distribution.  That is, if we use a Markov chain Monte Carlo (MCMC) algorithm to estimate the posterior distribution of our model parameters, this model may require the use of large multivariate Metropolis updates within our Gibbs sampler. \citet{besag:poisson} and \citet{knorr-held:rue} suggest a reparameterization of~\eqref{eq:bym_k} which involves integrating $\phi_{ikt}$ out of the model, yielding a Gaussian full conditional for $ {\bf Z} _{i\cdot\cdot}$ and requiring Metropolis updates for $\theta_{ikt} =  {\bf x} _{ikt} \ensuremath{\boldsymbol{\beta}}_{kt} + Z_{ikt} + \phi_{ikt}$.  Fortunately, the $\theta_{ikt}$ are independent of one another given $ {\bf Y} $ and the other model parameters, so these Metropolis updates can be conducted independently and in parallel.

We complete our hierarchical model by specifying the following prior distributions for our other model parameters:
a vague
prior for $ \ensuremath{\boldsymbol{\beta}}$, a weakly informative inverse Gamma prior for each $\tau_k^2$, a beta prior for each $\rho_k$, and an inverse Wishart prior for each $G_t$ with hyperparameter $G$, itself modeled using a Wishart prior.
While this
structure on the covariance matrices
is likely unnecessary given the number of spatial regions in the data \citep[see the discussion of prior sensitivity in spatial models by][]{bernardinelli}, this
comes at little-to-no computational cost (see {Web Appendix~A.6}) and offers a convenient means for specifying proper priors.
Putting these pieces together, our full hierarchical model is as follows:
\begin{align}\label{eq:hier}
\pi\left( \ensuremath{\boldsymbol{\beta}}, {\bf Z} ,G,G_1,\ldots,G_t, \mbox{\boldmath $\rho$} ,\tau_k^2, \mbox{\boldmath $ \theta $} \,\vert\, {\bf Y} \right) \propto&
\prod_{i,k,t} \mbox{$\text{Pois}$}\left(Y_{ikt}\,\vert\, n_{ikt}\exp\left[\theta_{ikt}\right]\right)
\times N\left( \mbox{\boldmath $ \theta $} \,\vert\, X \ensuremath{\boldsymbol{\beta}}+ {\bf Z} , \mbox{$\Sigma$} _{\theta}\right)\notag\\
&\times \mbox{$\text{MSTCAR}$}\left( {\bf Z} \,\vert\, G_{1},\ldots,G_{N_t}, \mbox{\boldmath $\rho$} \right) \notag\\
&\times \prod_{t=1}^{N_t} \mbox{$\text{InvWish}$}\left(G_t\,\vert\, G,\nu\right) \times \mbox{$\text{Wish}$}\left(G\,\vert\, G_0,\nu_0\right)\notag\\
&\times \prod_{k=1}^{N_g} \left[\mbox{$\text{Beta}$}\left(\rho_k\,\vert\, a_{\rho},b_{\rho}\right) \times \mbox{$\text{IG}$}\left(\tau_k^2\,\vert\, a_{\tau},b_{\tau}\right)\right] ,
\end{align}
where $ \mbox{$\Sigma$} _{\theta}$ is a diagonal matrix of size $N_sN_gN_t$ with elements $\tau_{k}^2$ and $X$ is the $(N_sN_gN_t\times p)$ matrix of covariates.

{While full details for implementing this model in an MCMC framework are provided in {Web Appendix~A}, we would be remiss to not discuss the computational burden associated with fitting a nonseparable model as opposed to a {separable} model; i.e., letting $\rho_k = \rho$ for $k=1,\ldots,N_g$ and $G_t = G$ for $t=1,\ldots,N_t$ corresponds to fitting a separable multivariate space-time model with $ \mbox{$\Sigma$} _{\eta} = R\left(\rho\right) \otimes G$.  First note that by using an AR(1) model for time, we can compute the $\widetilde{R}_{t,t'}^*$ elements of $ \mbox{$\Sigma$} _{\eta}$ in closed-form, reducing the burden of computing $ \mbox{$\Sigma$} _{\eta}^{-1}$ from an $N_tN_g \times N_tN_g$ matrix inversion to a series of $N_g\times N_g$ matrix inversions. Furthermore, while the nonseparable $\mbox{$\text{MSTCAR}$}$ model contains more parameters than its separable counterpart, the additional computational burden associated with its implementation in an MCMC framework is negligible.  Specifically, the computations necessary to construct the full conditional distributions for each $G_t$ are simply a partitioning of those necessary for constructing the full conditional distribution of $G$ in a separable model.}

When implementing this model in an MCMC framework, we have found that proper specification of the initial values can be crucial to facilitating convergence in a timely fashion. In particular, the parameters which require Metropolis updates --- $\rho_k$ and $\theta_{ikt}$ --- should be treated with care. For instance, we recommend initializing $\rho_k$ to be large (say, 0.90) if a high degree of temporal correlation is expected.
More importantly,
we recommend initializing
\begin{align}
\theta_{ikt}^{(0)} =
\begin{cases}
\log \left(Y_{ikt}\slash n_{ikt}\right), &\text{if $Y_{ikt} > \epsilon$}\\
\log \left(\sum_{i} Y_{ikt} \slash \sum_{i} n_{ikt}\right), &\text{if $Y_{ikt}\le \epsilon$}
\end{cases},
\end{align}
where $\epsilon\ge 0$ is some small nonnegative integer, as this will allow the model to learn about parameters such as $ \ensuremath{\boldsymbol{\beta}}$ and $ {\bf Z} $ early on in the process; in practice, letting $\epsilon=0$ has been sufficient.  When $\theta_{ikt}$ is poorly initialized, however, the MCMC algorithm may take a large number of iterations to recover,
resulting in a chain which is slow to converge.

{
\subsection{Assessment of reliability}\label{sec:reliable}
The primary motivation for this work is to achieve more reliable age-specific mortality rate estimates for these data.  In order to assess the reliability of our estimates, we will begin by generating \emph{synthetic} death counts from the posterior predictive distribution for $Y_{ikt}^{*}$,
\begin{align*}
Y_{ikt}^{*(\ell)}\,\vert\, \theta_{ikt}^{(\ell)} \sim \mbox{$\text{Pois}$}\left(n_{ikt}\exp\left[\theta_{ikt}^{(\ell)}\right]\right), \text{for $i=1,\ldots,N_s$, $j=1,\ldots,N_t$, $k=1,\ldots,N_g$},
\end{align*}
and for $\ell=1,\ldots,L$, where $\theta_{ikt}^{(\ell)}$ denotes the $\ell$th (post-burn-in) sample from the posterior distribution of $\theta_{ikt}$. We can then compute the 95\% CI for $Y_{ikt}^*$ from these posterior predictions and determine the proportion of each county's $N_tN_g$ 95\% CIs that contain the observed $Y_{ikt}$ to estimate the coverage probability. As a baseline for comparison purposes, we will compare our results to those generated from an empirical Bayesian Poisson-gamma model of the form
\begin{align}
 \mbox{\boldmath $\gamma$} \,\vert\,  {\bf Y} ,a_{kt},b_{kt} \propto \prod_{ikt} \mbox{$\text{Pois}$}\left(Y_{ikt}\,\vert\, n_{ikt}\gamma_{ikt}\right) \times \mbox{$\text{Gamma}$}\left(\gamma_{ikt} \,\vert\, a_{kt},b_{kt}\right)\label{eq:poisgam}
\end{align}
where $a_{kt} = \sum_{i} Y_{ikt} \slash \sum_i n_{ikt} \times b_{kt}$ and $b_{kt} = 1000$, indicating a prior distribution equivalent to 1,000 additional persons with an event rate equal to the national average. We can then generate synthetic death counts from the resulting posterior predictive distribution, $Y_{ikt}^{\dagger (\ell)}\,\vert\, \gamma_{ikt}^{(\ell)} \sim \mbox{$\text{Pois}$}\left(n_{ikt}\gamma_{ikt}^{(\ell)}\right)$ as before.  The method which yields coverage probabilities near 0.95 more consistently will be deemed more reliable.
}

\subsection{Tools for measuring temporal changes in mortality}
When studying temporarily-varying mortality rates, it is often of interest to measure the decline from the beginning of the study period to the end.  Letting $\lambda_{ikt} = \exp\left[\theta_{ikt}\right]$ denote the county-specific mortality rate for group $k$ at time $t$ and letting $\Delta_{ik}(t,t') = \left(\lambda_{ikt} - \lambda_{ikt'}\right)\slash\lambda_{ikt}$ denote the percent change from time $t$ to $t'> t$ for group $k$ in county $i$, an obvious choice would be to compute $\Delta_{ik}(1,N_t)$ for each county and each group.  Similarly, one could define $\Delta_{i}(t,t') = \left(\sum_{k} n_{ikt}\lambda_{ikt} - \sum_{k} n_{ikt'} \lambda_{ikt'}\right)\slash \sum_{k} n_{ikt}\lambda_{ikt}$ and compute $\Delta_{i}(1,N_t)$ as an estimate of the county's percent decline.

The drawback of these quantities is that they only account for the rates at the beginning and end of the study period, ignoring the intervening periods.  Here, we propose a measure we refer to as the
``saved person-years'' --- or SPY ---
which can be computed as
\begin{align}\label{eq:lsaa}
\mbox{$\text{SPY}$}_i &= \frac{1}{N_t-1} \sum_{t=2}^{N_t} \frac{\sum_k n_{ikt} \left(\lambda_{ik1} \left[1-\Delta_{\cdot k}(1,t)\right] - \lambda_{ikt}\right)}{\sum_{k} n_{ikt}} \times \text{100,000},
\end{align}
where $\Delta_{\cdot k}(1,t) = \left(\lambda_{\cdot k1} - \lambda_{\cdot k N_t}\right)\slash \lambda_{\cdot k 1}$ denotes the nationwide average decline from time 1 to time $t$ with $\lambda_{\cdot kt} = \sum_i n_{ikt} \lambda_{ikt} \slash \sum_i n_{ikt}$.
Note that $\mbox{$\text{SPY}$}_i$ is essentially a measure of the deviation between the expected rate for county $i$ at time $t$ if the county had declined at the rate of the national average --- i.e., $\lambda_{ik1} \left[1-\Delta_{\cdot k}(1,t)\right]$ ---
and the rate we estimate from the model. While this quantity should not replace the crude measure of the percent decline, it is a simple and easily interpretable investigative tool that
tells a more thorough story
of a region's trajectory over the study period, as we will demonstrate in our analysis.

In practice, we estimate $\theta_{ikt}$ by obtaining samples from its posterior distribution --- say, $\theta_{ikt}^{(\ell)}$ for iteration $\ell$ of our MCMC algorithm.  As such, both $\Delta_{i}(t,t')$ and $\mbox{$\text{SPY}$}_i$ can be computed for each iteration of our MCMC algorithm, resulting in posterior distributions for each of these quantities.  From this point forward, however, we restrict our attention to the posterior medians of $\Delta_{i}(t,t')$ and $\mbox{$\text{SPY}$}_i$ unless otherwise stated.

\section{Analysis of the Stroke Mortality Data}\label{sec:anal}
In the absence of covariate information, we fit the hierarchical model in~(\ref{eq:hier}) to the stroke mortality data described in Section~\ref{sec:data} using
an intercept term for each combination of year and age group \citep[as required when using an improper CAR model, per][]{besag95}, forcing the random effects to account for a substantial amount of the spatio-temporal variability in the data.
We place an informative beta prior on each $\rho_k$ to encourage higher temporal correlations in the model,
and we use a vague inverse Wishart prior for each of the $G_t$.
When running the MCMC algorithm,
we thinned our posterior samples for $\theta_{ikt}$ by removing 9 out of 10 samples --- while this is not theoretically necessary, it {reduced the burden of storing excess samples for our nearly 400,000 random effects}. Estimates provided are based on posterior medians, and 95\% credible intervals (95\% CI) were obtained by taking the 2.5- and 97.5-percentiles from the thinned post-burn-in samples. Additional figures, including animations displaying temporal evolutions in the geographic trends, can be found in {Web Appendix~B}.

Before delving into the epidemiologic findings, we evaluate some of the numerous variance parameters permitted by the use of a nonseparable model. While the $\mbox{$\text{MSTCAR}$}$ model in Section~\ref{sec:mstcar} is derived using temporally-varying covariance matrices, $G_t$, these parameters are not necessarily of direct interest as they are the variance parameters for $ {\bf v} _{\ell\cdot t}$, and thus they are \emph{not} directly interpretable with respect to the mortality rates. Instead, we need to use the posterior samples of $G_t$ and $\rho_k$ to construct $ \mbox{$\Sigma$} _{\eta}$ from~\eqref{eq:Sig_eta}.  These values coincide with the conditional covariance matrix of $ {\bf Z} _{i\cdot\cdot}$ (when scaled by the number of neighbors, $m_i$), and thus \emph{are} {interpretable with respect to the log mortality rates}.
{Furthermore, patterns in these parameters can be easily interpreted, as well.
For instance, Figure~\ref{fig:var} displays the posterior estimates for the diagonal elements of $ \mbox{$\Sigma$} _{\eta}$ corresponding to each age group.
Here, the declines in Figures~\ref{fig:var2} and~\ref{fig:var3} suggest a higher degree of spatial smoothing in later years than at the beginning of the study period, as the corresponding $Z_{ikt}$ become less free to deviate from their neighbors over time. As we will see, this may be in part due to declines in the mortality rates, themselves, as lower rates may also imply a smaller range of rates.}

\begin{figure}[t]
    \begin{center}
        \subfigure[Ages 65--74]{\includegraphics[width=.32\textwidth]{var_figs_6.png}\label{fig:var1}}
        \subfigure[Ages 75--84]{\includegraphics[width=.32\textwidth]{var_figs_7.png}\label{fig:var2}}
        \subfigure[Ages 85+]{\includegraphics[width=.32\textwidth]{var_figs_8.png}\label{fig:var3}}
    \end{center}
    \caption{Temporal evolution of the diagonal elements of $ \mbox{$\Sigma$} _{\eta}$.  While the scale of these parameters is not directly interpretable, declines in these variance parameters suggest an increase in spatial smoothing over time (particularly for those 75+).}
    \label{fig:var}
\end{figure}

{
To assess the reliability of our estimates, we follow the approach set forth in Section~\ref{sec:reliable}. We begin by generating $L=1000$ replicates for each $Y_{ikt}$ from the posterior predictive distribution corresponding to both the $\mbox{$\text{MSTCAR}$}$ model fit in~\eqref{eq:hier}, as well as the empirical Bayesian Poisson-gamma model from~\eqref{eq:poisgam}. After computing the 95\% CI of the replicates from both models and comparing these intervals to the observed $Y_{ikt}$, we find that the mean county-specific coverage from the $\mbox{$\text{MSTCAR}$}$ is 97.7\%, whereas the Poisson-gamma yields a coverage of 99.6\%.  What is immediately clear from this result is that while the intervals from the $\mbox{$\text{MSTCAR}$}$ are consistently more narrow than those from the Poisson-gamma model, this increase in precision due to the $\mbox{$\text{MSTCAR}$}$ model is valid.  {Figures related to this reliability assessment are included in Web Appendix~B}.
}

Having demonstrated the necessity and utility of the $\mbox{$\text{MSTCAR}$}$ model for these data, we now shift our attention to the rate estimates themselves.
For the sake of illustration of the temporal trends and the $\mbox{$\text{SPY}$}$ tool, we highlight two counties: one from the heart of the stroke belt (Jefferson County, AL), and one from the opposite side of the country (King County, WA).
Figure~\ref{fig:gt} displays time trends of the estimated rates for these two counties
for each of the age groups, along with the national averages.  Here, we observe that while these counties exhibit similar trends for the 65--74 age group --- where King County consistently outperforms both Jefferson County and the nation as a whole --- the trends for the remaining age groups are quite different.  This is particularly true for the eldest age group in our study, as Jefferson County experienced such a sharp decline from 1973 to 1990 so as to not only pass King County, but also the national average.  This period of consistent declines was followed by stagnant rates through the early 1990s and an \emph{increase} in rates until the early 2000s, leading to the county ending the study period among the worst in the nation {(albeit with much lower rates than in the 1970s)}.

\begin{figure}[t]
    \begin{center}
        \subfigure[Ages 65--74]{\includegraphics[width=.32\textwidth]{grouptrend_1.png}\label{fig:gt1}}
        \subfigure[Ages 75--84]{\includegraphics[width=.32\textwidth]{grouptrend_2.png}\label{fig:gt2}}
        \subfigure[Ages 85+]{\includegraphics[width=.32\textwidth]{grouptrend_3.png}\label{fig:gt3}}
    \end{center}
    \caption{Comparison of the rate trends over time for two selected counties (solid lines), along with the national average (dashed line).  Gray bands denote 95\% CI. Also displayed for each county is the expected rate had the county declined at the same rate as the national average (dotted lines).
    }
    \label{fig:gt}
\end{figure}

One aspect highlighted by these figures, however, is that computing a simple percent decline does not tell the whole story.  For instance, while Jefferson County experienced some degree of increasing rates in each of the three age groups during the late 1990s, simply measuring the age-standardized decline from 1973 to 2013 would overlook the strides the county took during the first half of the study period, when Jefferson County declined at a rate much faster than the national average.  The same cannot be said, however, for King County,
which underperformed during the first 30 years of the study, despite experiencing an overall decline that outperformed the national average.
This discrepancy can be observed by investigating the $\mbox{$\text{SPY}$}_i$ tool for each county.
Here, we find that Jefferson County saved 136.6 (98.8, 174.9) person-years per 100,000, while the $\mbox{$\text{SPY}$}_i$ for King County was -34.0 (-61.7, -5.3) person-years per 100,000, reinforcing our claim that
King County's declines lagged behind the national average.
A map of the $\mbox{$\text{SPY}$}_i$ values for all 3,099 counties can be found in Figure~\ref{fig:lsaa}, where we find that parts of the Deep South outperformed much of the nation, and a comparison
to the percent decline can be found in Figure~B.4 of the Web Appendix.}

\begin{figure}[t]
    \centering

        \includegraphics[width=.5\textwidth]{LSAA_noleg.png}
        \includegraphics[width=.10\textwidth]{SPY_leg.png}

    \caption{Map of the saved person-years (SPY)

    measure from~\eqref{eq:lsaa}, which measures the average difference between the model-estimated mortality rate of a county (per 100,000) and the expected mortality rate assuming a rate of decline equivalent to the national average.}
    \label{fig:lsaa}
\end{figure}

Turning our attention to the geographic patterns in stroke death rates presented in Figure~\ref{fig:rates}, we find substantial differences between age groups.
For the youngest subpopulation (ages 65--74), the clear geographic pattern shown in Figure~\ref{fig:6a} prominently highlights the so-called ``stroke belt'' in the rates from 1973, and the map of the percent declines in Figure~\ref{fig:6d} --- with large declines along the East Coast and smaller declines in the region stretching from Texas to the Dakotas --- seems to indicate the ``shift'' in the stroke belt identified by \citet{shiftingbelt}.
These patterns, however, are much less evident among the older age groups, especially for the eldest population.  Here, the rates in Figure~\ref{fig:8a} exhibit far less spatial clustering, while Figure~\ref{fig:8d} suggests that the rate of declines in mortality for those 85 and older was generally slower nationwide compared to those observed for those aged 65--74 and 75--84.

\begin{comment}
\begin{figure}[t]
    \begin{center}
        \subfigure[Ages 65--74: 1973]{\includegraphics[width=.45\textwidth]{Age/longquin_6_1.png}\label{fig:6a}}

        \subfigure[Ages 65--74: Declines]{\includegraphics[width=.45\textwidth]{declineleg_6.png}\label{fig:6d}}

    \end{center}
    \caption{Maps of the stroke mortality rates at the beginning of the study period (1973) and {the percent declines over the course of the study period}. Note that estimates for counties with fewer than 100 people in an age bracket in 1973 are suppressed.}
    \label{fig:rates}
\end{figure}
\end{comment}

\begin{sidewaysfigure}[t]
    \begin{center}
        \subfigure[Ages 65--74: 1973]{\includegraphics[width=.32\textwidth]{longquin_6_1.png}\label{fig:6a}}
        \subfigure[Ages 75--84: 1973]{\includegraphics[width=.32\textwidth]{longquin_7_1.png}\label{fig:7a}}
        \subfigure[Ages 85+: 1973]{\includegraphics[width=.32\textwidth]{longquin_8_1.png}\label{fig:8a}}\\
        \subfigure[Ages 65--74: Declines]{\includegraphics[width=.32\textwidth]{declineleg_6.png}\label{fig:6d}}
        \subfigure[Ages 75--84: Declines]{\includegraphics[width=.32\textwidth]{declineleg_7.png}\label{fig:7d}}
        \subfigure[Ages 85+: Declines]{\includegraphics[width=.32\textwidth]{declineleg_8.png}\label{fig:8d}}
    \end{center}
    \caption{Maps of the stroke mortality rates at the beginning of the study period (1973) and {the percent declines over the course of the study period}. Note that estimates for counties with fewer than 100 people in an age bracket in 1973 are suppressed.}
    \label{fig:rates}
\end{sidewaysfigure}

\begin{comment}
\begin{figure}[t]
    \begin{center}
        \subfigure[Ages 65--74]{\includegraphics[width=.32\textwidth]{bhamb_1.png}\label{fig:bham1}}
        \subfigure[Ages 75--84]{\includegraphics[width=.32\textwidth]{bhamb_2.png}\label{fig:bham2}}
        \subfigure[Ages 85+]{\includegraphics[width=.32\textwidth]{bhamb_3.png}\label{fig:bham3}}
    \end{center}
    \caption{Number of deaths (per 100,000) for Jefferson County, AL. {\color{red}need to add ``(per 100,000)'' to figures}.  Need to decide if this is worth putting into the paper.}
    \label{fig:bham}
\end{figure}
\end{comment}

\section{Discussion}\label{sec:disc}
This paper has extended the $\mbox{$\text{MSTCAR}$}$ model of \citet{quick:waller} to a generalized linear model for the purposes of analyzing a dataset comprised of county-level counts of stroke mortality, allowing for a nonseparable multivariate spatiotemporal dependence structure. Our analysis revealed spatiotemporal trends in stroke mortality that varied by age group, in addition to the nationwide reduction in
rates previously noted in the literature \citep[e.g.,][]{gillum:2011:stroke,linda}.  We also observed differing aspects of the western shift in
parts of the South
for each age group, as
identified in the total population by \citet{shiftingbelt}, and
explored the impact of non-linear trends in stroke mortality
via the SPY tool.

When modeling event rates for a rare event such as stroke mortality, it is important to leverage as much information as possible to achieve reliable estimates.
In addition to incorporating spatial structure into the model --- allowing for regional patterns such as the stroke belt to lend support to less populated counties --- the $\mbox{$\text{MSTCAR}$}$ accounts for temporal correlation between observations in consecutive years and multivariate dependencies, such as those between observations in different age-brackets.
These additional sources of information can be invaluable when dealing with outlying counts, a problem which particularly plagues counties with small population sizes where an increase/decrease of a single death can dramatically change the observed mortality rate.
Overcoming this challenge will be paramount in our future work, where we wish to use the $\mbox{$\text{MSTCAR}$}$ model to analyze
racial and gender disparities in heart disease and stroke mortality across various age groups.
While there is computational burden associated with implementing the $\mbox{$\text{MSTCAR}$}$ for such large multivariate datasets, these analyses will provide incredible insight into these various disparities.

\begin{comment}
Another area for future research is to assess how the effect of known risk factors for stroke mortality differ across the various age groups.  Based on a comparison of Figures~\ref{fig:6a}--\ref{fig:8a}, as well as Figures~B.4--B.6 in Web Appendix~B, it appears as though the degree of spatial clustering in the mortality rates declines from the 65--74 age group to the 85+ age group.
This could suggest that while stroke mortality amongst younger populations may be heavily influenced by environmental factors --- such as {diet, smoking, and the socioeconomic factors discussed by \citet{linda}} --- many of the risk factors for stroke mortality among those ages 85 and older may be less geographically structured.
\end{comment}

{Another area for future research is understanding the factors that contribute to differential geographic patterns by age group in both the baseline 1973 stroke mortality rates as well as the patterns of declining stroke mortality rates.
While it is well known that the risk for stroke increases with age, the spatiotemporal patterns of stroke mortality by age group have not been documented previously.  Hypotheses for understanding the observed differential spatio-temporal patterns in stroke mortality by age group include, but are not limited to, the following categories: 1) spatio-temporal differences in the relative contributions of decreasing case fatality rates and incidence rates by age group {\citep[e.g.,][]{case_fatality,geo_variation}}; 2) differential influence of living conditions (e.g. socioeconomic resources, access to quality health care, access to healthy food and recreational environments, etc.) or changes in those living conditions, on stroke mortality by age group {\citep[e.g.,][]{tassone,neighborhood}}; or 3) differences in the accuracy of death certificate reporting by age group due to more co-morbidities and competing conditions of death in the older ages.
}

Lastly, a topic worth discussing is that publicly available data for rare events such as deaths from stroke can be difficult to find due to data
privacy issues.
While data such as the number of stroke-related deaths are publicly available at the county-level from NCHS (via CDC Wonder), subsets of data with less than 10 events in a geographic region are suppressed beginning in 1989 \citep{cdc:sharing}.
This results in over {80}\% of the data points for those aged 65--74 and nearly 70\% of the over 380,000 data points used in this analysis being suppressed to the public.
To overcome this
privacy issue while still preserving utility, \citet{quick:zero} have explored the
risks associated with the generation of {synthetic} data
for rare event data for a single population for a single year of data.
As an area of active research,
we hope to use the methodology proposed here to generate reliable synthetic public-use data for small areas
which respects the spatial-, temporal-, and multivariate structures in the true data, thereby providing greater access to complete, high quality data.

\bibliographystyle{jasa}


\section{Introduction}
\label{Intro}

Criticality safety assessments require a prediction of the effective neutron multiplication factor (k$_{\rm eff}$) below a sufficient safety margin. This predicted value is derived using a validated calculation method with validated computer codes, e.g. so called criticality codes to calculate the k$_{\rm eff}$ of an application case. The validation of a criticality code can be achieved by recalculations of suitable critical experiments performed in laboratories and documented and evaluated e.g. in \cite{ICSBEP}. In recent years, several authors discussed the fact that depending on the application case and the choice of experiments, the effect of correlated experimental data on the determination of the bias, its uncertainty, and the resulting safety margins has to be considered \cite{Ivanova,NCSDCorr, NCSDBias, Hoefer2015, ICNCCorr, ICNCBias, Sobes, BJICNC, MONK}. The questions arising in the field of determination and handling of integral experimental covariance matrices in the process of code validation are also discussed in the Expert Group on Uncertainty Analysis for Criticality Safety Assessment (UACSA), a sub-group of the Working Party on Nuclear Criticality Safety (WPNCS) of the Nuclear Energy Agency (NEA) within the Organization for Economic Co-operation and Development (OECD). Actual questions which arose recently are: How to treat given sets of similar experimental data without knowing all exact statistical dependencies; and further, what are the implications on modeling these experiments in a code validation procedure regarding the consideration of the complete integral experimental correlation or covariance matrices?

In this manuscript we address these questions by following parts of the group’s proposal for a benchmark called Role of Integral Experiment Covariance Data for Criticality Safety Validation \cite{UACSABMIV}. In contrast to the benchmark proposal we focus on a reduced number of experiments but a total of nine different modeling approaches.

With the following analysis we add a new perspective to the discussion and show the effect of different modeling approaches for the same set of experimental data on the resulting integral covariance or correlation matrices.

Correlated data can arise if different experiments share parts of the experimental setup, measurement systems, or other relevant parameters. Some experiments described in the ICSBEP are not performed as single experiments, but slight variations of a setup were repeatedly investigated and published as a series of the same experiment. This is e.g. the case for LEU-COMP-THERM- 039 (LCT-39), where the number and location of empty positions  in a fuel rod grid were varied. In the following work we focus on the experimental data from experiments numbers 6, 7, and 8 from this series described in detail in \cite{ICSBEP, UACSABMIV, ICNCCorr} and references therein. The critical experiments consist of water moderated low enriched uranium fuel rods with a thermal neutron spectrum. The experimental setups are 22$\times$22 arrays consisting of 363 (459, 448) fuel rods for experiment 6 (7, 8) and 121 (25, 36) empty spots, respectively. For further details we refer to \cite{ICSBEP, ICNCCorr}. Clearly these experiments share certain components, and treating them as individual statistical independent data sets in the process of validation probably would not be appropriate. Hence, the determination of the integral covariance or correlation matrix of the experiments is a crucial step on the way to determine a bias of the calculated application case k$_{\rm eff}$.
\section{Methods and Parameters}
\label{Methods}

For the determination of the integral covariance matrices of k$_{\rm eff}$ and the corresponding correlation matrices we use a Monte Carlo Sampling approach, and SUnCISTT \cite{SUnCISTT} to steer and evaluate the numerous SCALE 6.1.2 \cite{SCALE} calculations. For two sets A and B of n sampled neutron multiplication factors k$_{\rm eff}$, the covariance ${\rm cov}^{\rm AB}$ is defined as
\begin{equation}
{\rm cov}^{\rm AB}=\frac{1}{n-1}\sum_{i=1}^n{ \left({k_{\rm eff}}^{{\rm A},i}-\overline{{k_{\rm eff}}}^{{\rm A},i}\right)
 \left({k_{\rm eff}}^{{\rm B},i}-\overline{{k_{\rm eff}}}^{{\rm B},i}\right) }
\label{Cov}
\end{equation}
with $\overline{{k_{\rm eff}}}$ symbolizing the expectation value of $k_{\rm eff}$, in our case the sample mean. The covariance can be interpreted as a measure of how much the $k_{\rm eff}$ of the two sets change simultaneously. A positive covariance indicates the following  monotonic connection between the two sets: large (or low) values in A correspond to large (or low) values in B. A negative covariance indicates the opposite behavior: large values in A correspond to low values in B. Due to its linearity the covariance gives only a tendency of the connection of two sets of random variables. To get comparable statements for more than two sets the covariance can be normalized with the standard deviation $\sigma$ to get the correlation coefficient ${\rm cor}$:
\begin{equation}
{{\rm cor}}^{\rm AB}=\frac{1}{\sigma_A \sigma_B}{\rm cov^{{\rm AB}}}
\label{ck}
\end{equation}
The correlation coefficient is a dimensionless measure of the linear dependence of two sets of random variables and takes values between +1 (complete positive linear connection) and -1 (complete negative linear connection). The confidence interval around the calculated ${\rm cor}$ is due to the value limitation of ${\rm cor}$ non-symmetric and must be determined by using transformations to so called Fisher’s distribution $z$ \cite{Fisher}. The latter is almost normally distributed and depends on the sample correlation coefficient:
\begin{equation}
z({{\rm cor}})=\frac{1}{2} {\rm ln}\left(\frac{1+{\rm cor}}{1-{\rm cor}}\right)
\label{zck}
\end{equation}
The corresponding tolerance intervals are calculated via $z\pm{\rm CL}\times\delta$ with the confidence level CL (e.g. 1.96 for the 95$\%$ confidence interval) and the sample standard deviation $\delta=\sqrt{(n-3)^{-1}}$.
The resulting z values are then transformed back to ${\rm cor}$ values. All given confidence intervals in the following work are 95$\%$ intervals.

Following a Monte Carlo Sampling approach, each value describing the experiment has to be interpreted as a distribution function. This means in turn, that the definition and interpretation of the experimental parameters and their uncertainties is essential. It strongly depends on the quality of the experimental data and availability of precise uncertainty specifications. To circumvent the problem of determining suitable distribution function for each parameter, we apply the ones given in \cite{UACSABMIV} which are also listed in \ref{Tab.1}. All experimental parameters are supposed to follow a uniform ${\rm U}$(a,b) or normal distribution N($\mu$,$\sigma$). Assuming the three experiments LCT-39 6, 7, and 8 to be statistical independent gives a correlation coefficient close to zero. Results for this assumption are shown for the correlation of k$_{\rm eff}$ values calculated by KENO V.a using the parameters given in table \ref{Tab.2} for 250 Monte Carlo samples for each experiment. The underlying model assumptions for the results of Figure 1 are very simple and straight forward: It is assumed, that the fuel rods are all identical in composition and position within its unit cell. In consequence, the modeling of one experiment consists basically of a $22\times22$ array of identical unit cells for the fuel rods and the empty positions respectively.
\begin{figure}[ht]
	\centering
		\includegraphics[width=.75\textwidth]{ScA_woCorr_Cbar_CI.eps}
	\caption{ Correlation coefficients of the experiments LCT-39 6, 7, and 8 assuming all fuel rods to be identical and no statistical dependence between each experimental setup.}
	\label{fig:Pic1}
\end{figure}

\begin{table}[hp]
	\centering
		\begin{tabular*}{.99\textwidth}{@{\extracolsep{\fill}}|p{5.2cm}|p{2cm}|r|}
			\hline
{\bf Model parameters} 	& {\bf Type of variation} 		&	{\bf Distribution functions} \\ \hline\hline
			Fuel diameter [cm]		 	& depends on scenario	& N(0.7892, 0.0017) \\ \hline
			Fuel lengths [cm]			&	depends on scenario	& N(89.7, 0.3) \\ \hline
			Fuel density [g/cm$^3$]			& depends on scenario	& N(10.38, 0.0133 ) \\ \hline
			Fuel content 234U	[At.-$\%$] &depends on scenario	& N(0.0307, 0.0005 )  \\ \hline
			Fuel content 235U [At.-$\%$]	&depends on scenario	& N(4.79525, 0.002)  \\ \hline
			Fuel content 236U	[At.-$\%$]&depends on scenario	& N(0.1373, 0.0005)  \\ \hline
			Boron concentration [atom/barn$\times$cm$\times10^{-8}$]&	depends on scenario	& N($6.9037$ ,
$0.8)$ \\ \hline
			Critical water height [cm] &	individual	& N($\mu,\sigma$) dep. on experiment\\ \hline
			Angle of fuel rod	& individual	&U( 0, 2$\pi$) \\ \hline
			Offset of grid hole x [cm]	&individual	& N(0, 0.00742) \\ \hline
			Offset of grid hole y	[cm] &individual	& N(0, 0.00742)  \\ \hline
			Hole diameter	[cm] &depends on scenario	& N(0.0105, 0.0085)  \\ \hline
			Inner cladding diameter	 [cm] &depends on scenario	&U (0.81, 0.83)\\ \hline
			Cladding thickness [cm] &	depends on scenario	&U (0.055, 0.065)    \\ \hline
		\end{tabular*}
	\caption{All model parameters and their distribution characteristics, following the suggestions of the benchmark proposal \cite{UACSABMIV}.}
	\label{Tab.1}
\end{table}

\begin{table}[htp]
	\centering
		\begin{tabular*}{.99\textwidth}{@{\extracolsep{\fill}}|l|r|p{4cm}|}
			\hline
			{\bf Code} & {\bf Parameter} & {\bf Value}\\ \hline
			{\multirow{4}{*} {KENO V.a}} & Nuclear data library & ENDF/B-VII (ce)	\\
															& Neutrons per generation & 10,000 \\
															& Skipped generations & 500 \\
															& $\sigma_{\rm MC}$ & 5$\times10^{-4}$(Sc. A to D); 1$\times10^{-4}$(Sc. E to H)\\ \hline
			SUnCISTT   & Number of samples & 250 \\ \hline

		\end{tabular*}
	\caption{Used codes and cornerstones of calculations. KENO V.a is taken from the CSAS5 sequence of SCALE 6.1.2}
	\label{Tab.2}
\end{table}
\section{Modeling Assumptions}
\label{Modeling}

Having determined all relevant parameters and their distribution functions, a calculation model is built to calculate the neutron transport equations and determine the neutron multiplication factor. Obviously the model should be as close as possible to the experimental setup to get reasonable results. However, in the statistical interpretation of experimental series like the one investigated in this article, the available data might leave some freedom of choice. The results shown in figure \ref{fig:Pic1} represent a model simplification of the experimental setup by assuming all fuel rods in one sample to be identical. However, it appears to be more reasonable that due to manufacturing tolerances of the experimental equipment individual fuel rods may vary in both, their individual composition and position within the unit cell. The position of the fuel rod in the unit cell is then limited by the grid hole. For some simplicity we assume the fuel rod to be always vertical, meaning a 90 degrees angle to the horizontal plane.
\begin{figure}
	\centering
		\includegraphics{FR_Model.eps}
	\caption{Modeling approach for the displacement of grid hole and fuel rod in a unit cell (not to scale)}
	\label{fig:FR_Model}
\end{figure}
The modeling approach for the fuel rod displacement is depicted in figure \ref{fig:FR_Model}. For the modeling of the experiments in KENO V-a, this implies each fuel rod to be simulated within its own unit cell, which we assume to have fixed dimension for all fuel rods. According to figure \ref{fig:FR_Model}, position of the grid hole might be displaced from the center of the unit cell by $\delta_x$   and $\delta_y$ in x and y direction. The center of the fuel rod itself might again be displaced in x- and y-direction, denoted by the radial displacement  R and angle $\theta$. In our modeling approach R is indirectly defined by the assumption that the fuel rod is in contact with the grid hole.

We chose eight different modeling approaches, scenarios A to H, depending on assumptions on the fuel similarity and position of each single fuel rod. Scenarios A to E assume the fuel to be identical for all fuel rods in all experiments. One can argue that this might be a reasonable approximation, based on the assumption of a very accurate fuel fabrication process with only tiny tolerances. However, scenarios F and G assume a set of the maximal needed fuel rods (484 for the 22$\times$22 grid array), each statistical independent. These fuel rods are placed for all experiments in a fixed position for scenario F or randomly for each experiment in scenario G, see figure \ref{fig:Scetch_ScF+G}. Finally, scenario H assumes all fuel rods in every experiment to be statistical independent. The statistical dependence of the fuel between two experiments decreases from scenario E to H. The results shown in figure \ref{fig:Pic1} assume the same modeling assumptions as scenario A but any correlations were neglected, meaning no statistical dependence between each experimental setup. This scenario is named ‘NoCor’.
\begin{table}[hp]
	\centering
		\begin{tabular*}{.99\textwidth}{@{\extracolsep{\fill}}|l|p{2cm}|p{2cm}|p{1.8cm}|p{1.8cm}|p{1.8cm}|}
			\hline
			{ Scenario} & { Grid hole displacement } &{ Grid hole diameter } &{
			Inner cladding diameter} &{ Cladding thickness } &{ Fuel variation } \\ \hline \hline
			NoCor & no; rod centered & shared  & shared & shared  & shared \\ \hline
			A			&	no; rod centered  & shared  & shared & shared  & shared  \\ \hline
			B			& $\delta_x$,$\delta_y$ & shared & shared	& shared &	shared \\ \hline
			C			&$\delta_x$,$\delta_y$	& individual &	shared &	shared &	shared \\ \hline
			D			&$\delta_x$,$\delta_y$  &	individual &	individual	&	shared	& shared \\ \hline
			E			&$\delta_x$,$\delta_y$  &	individual & 	individual	&	individual 	&	shared \\ \hline
			F			&$\delta_x$,$\delta_y$	&	individual	&	individual	&	individual	&	484 FR, fixed pos. \\ \hline
			G			&$\delta_x$,$\delta_y$	&	individual	&	individual	&	individual	&	484 FR, random pos.\\ \hline
			H			&$\delta_x$,$\delta_y$	&	individual	& individual	&	individual	&	individual\\ \hline
		\end{tabular*}
	\caption{Modeling assumptions for fuel rod geometries and compositions of one sample for each different scenario. The variation of fuel in the last column means the variation of the diameter, length, density and enrichment of the fuel as well as the boron impurity. Scenario NoCor (scenario H) is identical to A (scenario G), except for neglecting statistical dependencies between experiments.}
	\label{Tab.3}
\end{table}

\begin{figure}[hp]
	\centering
		\includegraphics[width=0.75\textwidth]{Scetch_ScF+G.eps}
	\caption{Sketch of modeling approaches for scenarios F (left) , G and H (right). The blue boxes represent the fixed library of individual generated fuel rods. The black and white squares represent the three experimental setups LCT-39 6, 7, and 8 from top to bottom. A black dot represents a fuel rod, a white one an empty spot. The simulation corresponding to the left part of the picture assumes a fixed position for each fuel rod in each experimental setup. E.g. FR1 is always in the top left spot for every experiment. The right part of the figure depicts the assumption of each fuel rod being randomly placed in the grid for each experiment.}
	\label{fig:Scetch_ScF+G}
\end{figure}

\section{Results}
\label{Results}

The analysis allover required a total of 6.750 SCALE inputs with up to 20.000 lines per input file. The calculations were performed using a total of 55.000 CPU-h and 882 TByte-h. The results then were processed and statistically analyzed using SUnCISTT.

The resulting $k_{\rm eff}$ values for each experimental data set and modeling scenario are shown in figure \ref{fig:keffs_StdDev}. We found a good agreement within the 2-$\sigma$ range of the experimental data ${k_{\rm eff}}^{\rm exp}=1.0$ ($\pm0.0012$ for exp. 7,8) ($\pm0.0009$ for exp. 6) given in \cite{ICSBEP}) and our results. The SCALE calculations with the applied continuous energy library ce\_v7\_endf (based on ENDF/B-VII) in the CSAS5 sequence systematically underestimates $k_{\rm eff}$, which is a known effect for low enriched uranium setups \cite{SCALEVal}. The larger error bars of the Monte Carlo approach of Scenario A in comparison to the error propagation approach done in the ICSBEP Handbook are not attributed to a general difference between the two methods. They rather arise from a different interpretation of the system parameter uncertainties. In the original experiment description \cite{Valduc} the uncertainty of the inner cladding diameter ($\pm$ 0.01 cm) and the cladding thickness ($\pm$ 0.005 cm) are reported to be independent. The uncertainty of the outer cladding diameter is obtained by error propagation. In the ICSBEP evaluation, the uncertainty of the cladding thickness is split equally between inner and outer diameter. This results in an uncertainty of the outer cladding diameter of $\pm$ 0.0025 cm, which reduces its impact on the uncertainty of $k_{\rm eff}$ significantly. The original evaluation assumes further a Gaussian distribution by dividing the half tolerance by $\sqrt{3}$. The resulting distributions for the outer cladding diameter and their impact on $k_{\rm eff}$ are different for both considerations. However, the values based on the original literature were used and for further details we refer to \cite{ICNCCorr}. Scenarios E to H in which the individual variations of the parameters partly cancel out each other, have significant lower error bars.

\begin{figure}[hp]
	\centering
		\includegraphics[width=0.75\textwidth]{keffs_StdDev.eps}
	\caption{Resulting $k_{\rm eff}$ values for 250 Monte Carlo samples for each experiment and scenario. The error bars indicate the standard deviation. While the nominal values remain fairly constant for all scenarios, the standard deviation decreases significantly from scenario D to F.}
	\label{fig:keffs_StdDev}
\end{figure}

In figure \ref{fig:All_9_Sc_CI} nine colored plots are shown for the ${\rm cor}$ values of the experiments LCT-39 6, 7, and 8 as well as the ${\rm cor}$ value and the 95$\%$ confidence interval.

The results show correlation coefficients around 0 for the scenarios NoCor and H, as expected, since there are no relevant parameters with shared values between the individual experiments. Note that the difference between the NoCor and H scenario is the variation of the fuel rods: In contrast to NoCor, in scenario H each fuel rod in each experiment is simulated individually and statistically independent. This difference is mapped in the sensitivity plots in figure \ref{fig:Fig7} and \ref{fig:Fig8} which show the correlation coefficients of each parameter with the resulting $k_{\rm eff}$. Note, that the performed sensitivity analysis shows the impact of the actual variation of each parameter on the $k_{\rm eff}$ uncertainty. However, we do not perform a sensitivity analysis by varying only one parameter at a time. This means, our sensitivities depend on the chosen distribution functions and their characteristics. Changing these assumptions in our approach might lead to a different sensitivity profile. This approach was chosen since we are interested in determining the contribution of each varied input parameter on the uncertainty of $k_{\rm eff}$ for given modeling assumptions.

While for the scenario NoCor the most relevant parameters are the cladding inner radius and thickness, the only important parameter for scenario H is the critical water height. It is notable, that in this case the different interpretations of the given experimental data lead to comparable ${\rm cor}$ values but totally different sensitivity profiles.  The highest correlation coefficients for scenario NoCor are the ones for the cladding inner diameter and thickness, and for the radius of the fuel. The only dominant parameter for scenario H is the critical water height.
Scenarios A, B, C and D show all ${\rm cor}$ values close to 1 with only little deviations between the different ${\rm cor}$ values of the scenarios. The ${\rm cor}$ values for scenario A to C are even the same within the 95$\%$ confidence interval. Their corresponding sensitivity profiles show huge similarities: The three largest ${\rm cor}$ values are the cladding inner radius and thickness and the radius of the fuel. For scenario B and C the U-235 weight-$\%$ plays a more prominent role. The sensitivity profiles for scenario D show a different behavior since the inner cladding diameter here is varied individually for each fuel rod. The leading contribution to the sensitivity profile now solely results from the cladding thickness.
\begin{figure}[hp]
	\centering
		\includegraphics[width=0.95\textwidth]{All_9_Sc_CI.eps}
	\caption{Matrices of the correlation coefficients for the 9 scenarios. The upper part of each matrix shows the ck values and the ranges for the 95$\%$ confidence interval. The lower part shows the nominal ${\rm cor}$ values color coded. The scale varies from dark red for ${\rm cor}$  = 1 to white for ${\rm cor}$ = -0.2. The resulting correlation coefficients for scenario NoCor and H are 0 within the 95$\%$ confidence interval. Scenarios A to D show ${\rm cor}$ values close to 1, while scenarios E to G show a slight decreasing of ${\rm cor}$ from approximately 0.75 to 0.5.}
	\label{fig:All_9_Sc_CI}
\end{figure}

Scenario E shows correlation coefficients of approximately 0.75 between the experiments 6 and 7 and 0.7 between experiments 7 and 8. The difference between the ${\rm cor}$ values is due to lower number of fuel rods in experiment 6 compared to 7 and 8. Thus, the individual variation of the cladding inner radius and thickness for each fuel rod affects the correlation coefficient of the experiments 7 and 8 more. The corresponding sensitivity profile shows the fuel radius as the leading parameter. A mild impact is shown by the fuel density and critical water height (additional the weight-$\%$ for U-235 and the fuel height for LCT-39 7).

Scenarios F and G show significantly smaller correlation coefficients between the experiments. The difference to Scenario E is that now also the fuel content of each fuel rod is varied individually. This can be seen in the sensitivity profile of both scenarios in figure \ref{fig:Fig8}, where the dominant parameter is the critical water height.

\begin{figure}[htbp]
	\centering
		\includegraphics[width=0.75\textwidth]{ScVsCk.eps}
	\caption{Correlation coefficients for the pairs of experiment 6 and 7 (black), 6 and 8 (red), and 7 and 8 (blue) of Series LCT-39 for the 9 scenarios. The error bars represent the 95$\%$ confidence interval. Scenarios A, B, C, and D show high correlation coefficients, all close to 1. For scenarios E to G the coefficients decrease to approximately 0.7, 0.6, and 0.5 respectively. The correlation coefficients for scenario H are 0 within the 95$\%$ confidence interval, like the ones for the NoCor scenario.}
	\label{fig:ScVsCk}
\end{figure}
\begin{table}[hp]
	\centering
	\begin{tabular*}{.99\textwidth}{@{\extracolsep{\fill}}|r|l|r|l|}
			\hline
					Key word 	&	Parameter	&	Key word	&	Parameter \\ \hline\hline
				rad\texttt{\_}CladIn	&	cladding inner raduis	&	w$\%$\texttt{\_}U235	& weight-$\%$ U-235 \\ \hline
				{thick\texttt{\_}Clad}	&	cladding thickness	&		w$\%$\texttt{\_}U236	&	weight-$\%$ U-236 \\ \hline
				{rad\texttt{\_}Fuel}		&	fuel radius	&	{height\texttt{\_}Water}	&	water height \\ \hline
				{height\texttt{\_}Fuel}	&	fuel height	&	{rad\texttt{\_}Hole}	&	hole radius \\ \hline
				{dens\texttt{\_}Fuel} 	&	fuel density	&	{delta\texttt{\_}Hole\_X}	& $\delta_x$ \\ \hline
				{dens\texttt{\_}B10}		&	B-10 density	&	{delta\texttt{\_}Hole\_Y}	&	$\delta_y$ \\ \hline
				w$\%$\texttt{\_}U234	&	weight-$\%$ U-234	&	{angle\texttt{\_}Rod}	&	$\theta$ \\ \hline

		\end{tabular*}
	\caption{Parameters and key words used in figure \ref{fig:Fig7} and \ref{fig:Fig8}}
	\label{Tab.4}
\end{table}

\begin{figure}[hp]
	\centering
		\includegraphics[width=0.95\textwidth]{Fig7.eps}
	\caption{ Correlation coefficients of each individual model parameter listed in table 4 and the resulting 250 $k_{\rm eff}$ values for scenarios NoCor and A. The profiles are similar due to identical modeling assumptions. The only difference between the two scenarios is the assumption of correlations in scenario A. Note the different scale of the color coded representation of ${\rm cor}$ w.r.t. figure \ref{fig:All_9_Sc_CI}.}
	\label{fig:Fig7}
\end{figure}

\begin{figure}[hp]
	\centering
		\includegraphics[width=0.95\textwidth]{Fig8.eps}
	\caption{Correlation coefficients for each model parameter and the resulting 250 $k_{\rm eff}$ values. From top to bottom the scenarios B to H are shown. The increased number of parameters compared to the scenarios NoCor and A is due to the variations of single fuel rods. The color coded representation of the ${\rm cor}$ shows values from -1 to 1.}
	\label{fig:Fig8}
\end{figure}

\section{Discussion and Conclusions}
\label{Conclusions}

In this work we discussed nine different modeling approaches for a given set of experimental data, leading to different correlation coefficients and sensitivity profiles. Using a Monte Carlo approach, we calculated 250 samples for each experimental setup and scenario to obtain the resulting $k_{\rm eff}$ values (figure \ref{fig:keffs_StdDev}). The steering and analyzing of the SCALE6.1.2 CSAS5 sequences where done using SUnCISTT.

Within each scenario we calculated for each pair of experiments the corresponding correlation coefficient and the 95$\%$ confidence intervals of this coefficient (figures \ref{fig:All_9_Sc_CI}, \ref{fig:ScVsCk}). We showed for each experiment and scenario the impact of the variation of each input parameter on the resulting $k_{\rm eff}$ by calculating the corresponding correlation coefficients (figure \ref{fig:Fig7}, \ref{fig:Fig8}).

For the combination of water moderated, low enriched Uranium rods modeled with the criticality code KENO V-a we found that the correlation coefficients between the $k_{\rm eff}$'s of the experiments LCT-39 6, 7, 8 varied between 0 and 1 within the 95$\%$ confidence interval.

The modeling assumption leading to scenario A and NoCor are identical, except that NoCor neglects correlations completely. The same holds for scenarios G and H, which are identical, but H neglects the correlations between experiments due to fuel similarities.

Varying all geometrical parameters affecting the outer cladding radius for each fuel rod separately leads to a significant decrease of the resulting correlation coefficient compared to the results derived from the assumption of all fuel rods having identical geometrical parameters (however, scenarios B, C and D still comprise high correlations). We found a significant drop of the ${\rm cor}$ value from scenario D to E (\ref{fig:ScVsCk}) as well as a significant drop of the $k_{\rm eff}$ uncertainty (figure \ref{fig:keffs_StdDev}). The main contribution to the $k_{\rm eff}$ uncertainty in scenario D stems from the cladding thickness, which in scenario E plays no role since it is varied for each fuel rod individually (figure \ref{fig:Fig8}).

Scenarios F and G show a further drop of the correlation coefficient, but within the 95$\%$ confidence interval the ${\rm cor}$ values of the two scenarios overlap (figure \ref{fig:ScVsCk}). The difference of the two assumptions, knowing the exact position of each fuel rod for Scenario F or randomize their position in the grid for scenario G has a comparable smaller effect on ${\rm cor}$ than the assumption of a finite number of fuel rods. It is notable, that the sensitivity analysis shows the sole dependence of the $k_{\rm eff}$ uncertainty on the critical water height (figure \ref{fig:Fig8}).

The different modeling assumptions might all be justified based on expert judgment. However, the sensitivity analysis reveals different sensitivity profiles, especially from scenario C to F. One could be tempted to choose the modeling assumptions based on the quality of the experimental data. As an example one could argue to choose scenario F or G, since the uncertainty of $k_{\rm eff}$ is much lower and the almost sole dependence of the $k_{\rm eff}$ uncertainty is on the critical water height. Following this argumentation, one could construct modeling assumptions based on the given experimental data to reduce uncertainties and to circumvent possible gaps in the data. But one has to be very careful with these options, and give very good arguments, why one chooses one scenario over another. The resulting covariance matrices directly influence the bias and its uncertainty, and thus the resulting upper sub-critical limit \cite{ Hoefer2015,ICNCBias, Sobes}.

Using the covariance or correlation matrices for the purpose of validation or the determination of the upper sub-critical limit of an application case, the results can vary strongly, depending on the scenario. Following the argumentation of \cite{ICNCBias, Sobes}, a rule of thumb is that the higher the correlation coefficient, the lesser information is available, and thus the upper sub-critical limit decreases. This means, that being not able to distinguish between the different scenarios and identify the correct one based on the available data, one would in this case take the results associated with the highest correlation coefficient to get a more conservative estimate of the bias in code validation or the upper sub-critical limit.

Note that the underlying data for the work presented is partly constructed and fictive as it is a part of a calculation benchmark exercise \cite{UACSABMIV}. From the given data, any modeling assumption from scenario A to H could be justified. For further determination of the scenarios one would need to know e.g. if the fuel content and geometric description for each fuel rod was identical or if it varied. The statements presented above thus are only valid for the combinations of code and experiments discussed here.

To derive more general statements, further investigations have to be carried out. On the other hand, it may be problem dependent if and to what extent the regard for correlations between benchmark experiments could influence the bias determination. A sufficient number of statistical independent data sets, e.g. for experiments conducted in different laboratories using different materials, can always circumvent the problem of the correct determination of integral experimental covariance data. However, the accurate consideration of correlated data seems to be inevitable if the experimental data in a validation procedure is limited.

But even if one can avoid the determination of the accurate integral experimental covariance data due to statistically independent data sets, the selected modeling scenario should always be justified. The modeling assumptions have the potential to decrease the uncertainty of the resulting $k_{\rm eff}$ significantly.

\section{Acknowledgement}
Our work was financed by the German Federal Office for Radiation Protection under grant No.3614R03331. We thank W. J. Marshall and B. T. Rearden for the various fruitful discussions of our results and findings. We further would like to thank D. Mennerdahl for the discussions concerning the different scenarios.

\section{References}
\label{Ref}