# Bayesian Hierarchical Modeling and Analysis for Physical Activity
Trajectories Using Actigraph Data

Pierfrancesco Alaimo Di Loro
Dpt. of Statistical Sciences
University of Rome "La Sapienza"
<EMAIL_ADDRESS>&Marco Mingione
Dpt. of Statistical Sciences
University of Rome "La Sapienza"
<EMAIL_ADDRESS>&Jonah Lispsitt
Fielding School of Public Health
University of California, Los Angeles
<EMAIL_ADDRESS>&Christina M. Batteate
Center of Occupational and Environmental Health
University of California, Los Angeles
<EMAIL_ADDRESS>&Michael Jerrett
Fielding School of Public Health
University of California, Los Angeles
<EMAIL_ADDRESS>&Sudipto Banerjee
Fielding School of Public Health
University of California, Los Angeles
<EMAIL_ADDRESS>
###### Abstract

Rapid technological developments in accelerometers have generated substantial
interest in monitoring human activity. Wearable devices, such as wrist-worn
sensors that monitor gross motor activity (actigraph units) continuously
record the activity levels of a subject, producing massive amounts of high-
resolution measurements. Analyzing actigraph data needs to account for spatial
and temporal information on trajectories or paths traversed by subjects
wearing such devices. Inferential objectives include estimating a subject’s
physical activity levels along a given trajectory; identifying trajectories
that are more likely to produce higher levels of physical activity for a given
subject; and predicting expected levels of physical activity in any proposed
new trajectory for a given set of health attributes. We devise a Bayesian
hierarchical modeling framework for spatial-temporal actigraphy data to
deliver fully model-based inference on trajectories while accounting for
subject-level health attributes and spatial-temporal dependencies. We
undertake a comprehensive analysis of an original dataset from the Physical
Activity through Sustainable Transport Approaches in Los Angeles (PASTA-LA)
study to formally ascertain spatial zones and trajectories exhibiting
significantly higher levels of physical activity.

Keywords: Bayesian hierarchical models; Directed acyclic graphs; Gaussian
processes; Physical activity; Sparsity; Spatial-temporal statistics.

##  1 Introduction

Promoting a healthy lifestyle continues to stoke substantial research
activities in public health. The “Physical Activity Guidelines for Americans”
(2nd edition) suggests that most individuals, depending on age and body
composition, receive 150-300 minutes of moderate to vigorous physical activity
(MVPA) weekly (Piercy et al., 2018). In general, the scientific community
agrees that regular physical activity (PA) can have immediate and long-term
health benefits (Reiner et al., 2013; Bull et al., 2020). Despite these well-
known benefits, most Americans fail to meet recommended requirements (Piercy
et al., 2018). Specifically, only 1 in 5 high-school adolescents and 1 in 4
adults meet recommended levels of physical activity (PA). Given the well-
established relationships between lack of PA and several obesity-related
chronic conditions such as heart disease, type 2 diabetes, and cancer, as well
as many physical and mental health benefits, an urgent need exists to improve
monitoring of PA and to establish public health programs that promote more PA
111More details at
https://www.cdc.gov/chronicdisease/resources/publications/factsheets/physical-
activity.htm.

Technologies for monitoring spatial energetics (James et al., 2016; Drewnowski
et al., 2020) and promoting physical activity continue to emerge. _Actigraphy_
broadly refers to the monitoring of human rest and activity cycles using
wearable devices. Actigraphy data are gathered directly from wearable sensors
or indirectly through smart-phone mobile applications and record repeated
measurements at very high resolution. Accelerometers, in particular, are
motion sensors that measure acceleration along different axes and are able to
collect large amounts of data (Plasqui and Westerterp, 2007; Sikka et al.,
2019). They are increasingly conspicuous because of their affordability,
accuracy, and availability in smart-phones, smart-watches and other wearable
devices. Many devices include Global Positioning System (GPS) sensors that
reference measurements with location tracking along trajectories, or paths,
traversed by the subject. Collected data can be quickly downloaded and
promptly analyzed to obtain insights into their pattern and structure.

We pursue a comprehensive analysis of an original actigraphy data set from the
Physical Activity through Sustainable Transport Approaches in Los Angeles
(PASTA-LA) study. Analyzing such data is sought for several reasons: (i)
estimating a subject’s physical activity levels along a given trajectory; (ii)
identifying trajectories that are more likely to produce higher levels of
physical activity for a given subject; and (iii) predicting expected levels of
physical activity in any proposed new trajectory for a given set of health
attributes. Researchers find actigraphy tracking especially attractive as it
allows for a better understanding of what behavioral and environmental factors
influence population and individual health and, hence, aid in public health
recommendation and policy.

Given that actigraphs generate large amounts of spatial-temporal data, it is
natural choose from the rich classes of such models (Cressie and Wikle, 2011;
Gelfand et al., 2010). However, actigraph data present some notable challenges
(Kestens et al., 2017): they exhibit dependence along trajectories and must be
accounted while predicting PA along arbitrary (unobserved) trajectories. We
argue against a customary spatiotemporal process over $\mathbb{R}^{2}$ and
disentangle spatial effects from dependence along trajectories. The balance of
the paper is organized as follows. Section 2 introduces the PASTA-LA data set
with insights into accelerometry data. The model for the temporal correlation
is introduced in Section 3, while spatial effects are discussed in Section
3.4. An extensive simulation study validating our model is proposed in Section
3.5. Data analysis, including model assessment and comparisons, are presented
in Section 4. Finally, we conclude with a discussion in Section 5.

##  2 Data

###  2.1 Data collection

Our dataset is compiled from the original Physical Activity through
Sustainable Transport Approaches in Los Angeles (PASTA-LA) study conducted on
a cohort of 460 individuals monitored between May 2017 and June 2018. Data
were collected through different sources: online questionnaires, a smartphone
app named MOVES, a GPS device (GlobalSat DG-500), and a wearable Actigraph
unit (Actigraph GT3X+). We focus on the values recorded by Actigraph and GPS
devices worn by $134$ subjects for two one-week periods (one in 2017 and one
in 2018). Study protocol for safeguarding participant information received
necessary institutional review board (IRB) approval. The data were stored on a
secure computer and a redacted version was created for purposes of data
sharing.

##### Questionnaires

The online questionnaires included two baseline and four follow-up surveys:
one baseline and two follow-ups for each collection period of the actigraph
and GPS data. Each survey consisted of responses pertaining to the
participant’s demographics and transportation habits. This survey is the first
baseline questionnaire and contains personal information such as sex, age,
BMI, ethnicity and other socioeconomic factors. A user ID was assigned to each
survey response data and a redacted master key was generated using all ID
types for joining with other study data.

##### Actigraph

The Actigraph unit is an accelerometer roughly the same size and weight of the
average wrist-watch that can be worn on the wrist, hip, and thigh and measure
the directional acceleration at a pre-specified time frequency (generally 10
Hz to 30 Hz). The Actigraph GT3X+ model used for the PASTA-LA study can detect
movement in up to three orthogonal planes (anteroposterior, mediolateral, and
vertical). Data are stored in an internal memory and can be downloaded to
other hardware for analysis through a proprietary software. During download,
the software converts the raw acceleration information to activity counts,
step counts, caloric expenditure, and activity levels, aggregated at the level
of sample epochs that can be specified by the user. The proprietary software
precludes recovering the raw data once they have been processed in the
download phase.

Our study asked the $134$ participants to keep the Actigraph unit on them
(wrist) at all times other than during bathing and sleeping (awake time was
assumed approximately from 7am to 11pm). Observations recorded outside this
daily time-window were excluded. When participants arrived at the research
offices to drop off devices, some described issues of efficacy in the ability
to keep the device on or charged. Indeed, while the actigraphs were supposed
to hold a charge long enough to last the whole week, this was not always the
case (possibly due to external conditions affecting the battery life or
variations in manfacturing). This resulted in a large amount of missing data.
During download the data were aggregated in sample epochs spanning $10$
seconds. Measurements included the activity counts for the three axles and
step counts. Time-stamps for final measurements (hour, minute, and second)
were referenced by the mid-point between the beginning and the end of the
epoch.

(a)

(b)

(c)

(d)

Figure 1: Derived missing data patterns for actigraph data during the daytime
for 4 individuals arising from querying the inclinometer and eliminating low
values: blank spots are missing data, black spots are observed values.

The Actigraph GT3X+ automatically records a measure of inclinometer values on
how many of the $10$-second epochs have been spent by the individual lying
down (inclinometer.lying), sitting (inclinometer.sitting), standing
(inclinometer.standing) or without wearing the accelerometer
(inclinometer.off). These variables were not of primary interest to the PASTA-
LA study and were not directly addressed in their data collection protocol for
quality assurance. Inclinometers222see https://actigraphcorp.com for details,
have reported sizable error rates up to $30\%$ depending upon where they are
worn and are likely less accurate when worn on the wrist (Peterson et al.,
2015). We sought to exploit convergence of accelerometry and inclinometer data
to derive periods of inactivity. We checked that large values of
inclinometer.off corresponded to low ($\sim 0$) values of activity in all the
possible endpoints. Then, we dropped observations with inclinometer.off larger
than $5s$ (i.e., the accelerometer was inactive for more than half of the
epoch). This yields “missing data” from querying the inclinometer and
eliminating low values. These resulted in $\approx 6.4\times 10^{6}$ scattered
observations exhibiting missingness patterns as in Figure 1.

##### GPS

The GlobalSat DG-500 was provided to the $134$ subjects, which recorded the
subject’s location (latitude and longitude) every $15$ seconds and comprised
date and time of localization and speed in kilometers per hour (computed as
distance over time through linear interpolation). In order to avoid a
geographical imbalance that could bias and invalidate the model estimates, for
the current analysis we restricted attention to subjects living and working in
the Westwood neighborhood of Los Angeles.

##### Joining

GPS and accelerometer data were all assigned a participant ID aligned with the
questionnaires’ master-key to facilitate joining across all ID types
(including email) while redacting and encrypting user data. The first baseline
questionnaire, Actigraph and GPS were available for our group of $134$
individuals. Henceforth, we refer to this specific group of units. We then
build two different data sets:

  * •

The first dataset, $D_{1}$, comprising $N\simeq 5\times 10^{6}$ measurements
is obtained by joining the first baseline questionnaire with Actigraph data
and includes the MAG (Section 2.2) at the different timestamps and all the
individual predictors, but no spatial information.

  * •

The second dataset, $D_{2}$, consists of $N\simeq 5\times 10^{5}$ measurements
(Figure 2) and is obtained by joining $D_{1}$ with GPS data. Actigraph and GPS
data are temporally misaligned so we joined them by linearly interpolating the
GPS locations on the temporal grid from the actigraph, keeping only those
interpolated values wth subsequent GPS measurements less than 30 seconds
apart. This is a reasonable assumption, given that the individual’s trajectory
is well-approximated by a piece-wise linear GPS trajectory.

(a)

(b)

Figure 2: (a): Observed locations over the Westwood area. (b): Observed
physical activity level over the Westwood area on a subset of $10$ individuals

###  2.2 Measure of physical activity

Our primary endpoint of analysis is the MAG (Magnitude of Acceleration)
defined as:

$MAG_{kt}=\sqrt{x_{kt}^{2}+y_{kt}^{2}+z_{kt}^{2}},\qquad k=1,\dots,K,\quad
t=t_{k1},\dots,t_{kT}\;,$ (1)

where $t_{kj}$ is the $j$-th time point for the $k$-th individual; and $x$,
$y$ and $z$ are the activity counts of the three axes (Ott et al., 2000).
There are substantial investigations into the statistical relationships
between accelerometer measurements and energy expenditure measures (EE)
(Freedson et al., 2012; Taraldsen et al., 2012) and, in particular, the
Metabolic Equivalent of Task (METs) which is currently the standard measure of
rate of activity intensity (Crouter et al., 2006; Hall et al., 2013; Lyden et
al., 2014). mig (2017) offers an extensive review of proposed accelerometer
measurement cut-points and transformation into physical activity metrics. More
specifically, Sasaki et al. (2011), Santos-Lozano et al. (2013) and Kamada et
al. (2016) investigated axis counts and vector magnitude resulting from the
GT3X+ accelerometer in both controlled and free-living environments, while
Aguilar-Farias et al. (2019) investigated the accuracy of relationships
between MAGs and METs in comparison to those based on the vertical axis counts
only and validated the results with the EE and METs as quantified by a
portable calorimeter.

The MAG-to-MET relationship expounded in Sasaki et al. (2011) is expressed as
a function of the MAG per minute, which we rescale to our $10$ seconds
aggregated counts as:

$\text{MET}_{kt}=(0.000863\times 6)\cdot MAG_{kt}+0.668876,$ (2)

and perform the same to the corresponding cut points for different PA
intensity level (Table 1).

Activity intensity | MET range | MAG
---|---|---
Sedentary or light | $[0,\,3)$ | $[0,\,493)$
Moderate | $[3,\,6)$ | $[493,\,1029)$
Hard | $[6,\,9)$ | $[1029,\,1608)$
Very hard | $[9,\,\infty)$ | $[1608,\,\infty)$
Table 1: MAG activity count cut-points for different PA intensity levels

Based upon the aforementioned literature, inference for the MAG is transformed
into METs through (2) to interpret results from a physical activity
perspective. Nevertheless, equations directly relating accelerometer
measurements with physical activity metrics in free-living studies must be
interpreted with caution. Relationships between MAG and MET have been posited
in controlled studies and validated while patients are performing specific
tasks (i.e. walking on a treadmill, gardening etc.). The relationship between
the recorded movement (acceleration) and the corresponding energy expenditure,
can vary significantly across different tasks affecting the reliability of
acceleration-based energy expenditure metrics (Lyden et al., 2011; Freedson et
al., 2012; Montoye et al., 2018).

##  3 The model

The outcomes corresponding to the $K$ subjects are referenced with respect to
the time at which they are recorded and the position in the trajectory. While
it is tempting to work with a spatiotemporal process, dependence introduced by
such processes may not be appropriate. An individual can visit the same
location numerous times in his/her trajectory. These revisits need not occur
at regular intervals and can be at distant time points. This suggests that
proximity of two spatial locations in a trajectory need not result in strongly
dependent MAGs recorded there. It appears more reasonable to model dependence
among MAG measurements through a temporal process. In fact, such temporal
processes can be motivated by the position vectors defining the trajectories
as we describe below.

Let $Z_{k}(\cdot):\mathbb{R}^{2}\rightarrow\mathbb{R}$ be a spatial process
corresponding to individual $k$. The domain of $Z_{k}(\cdot)$ is restricted to
the trajectories
$\gamma_{k}(t)=\left(\gamma^{x}_{k}(t),\gamma^{y}_{k}(t)\right)$, where
$k=1,\dots,K$ and $t\in\mathbb{R}^{+}$, which defines the movements of the
$k$-th individual along time. As shown in Figure 3, the process actually
belongs to a one-dimensional space, for which we define a proper distance
measure $d(t_{ki},t_{kj})=\|\gamma_{k}(t_{kj})-\gamma_{k}(t_{ki})\|$, where
$t_{ki}$ is the $i$-th recorded time point from individual $k$. We approximate
such distances as the elapsed time between the two points
$d(t_{ki},t_{kj})=|t_{kj}-t_{ki}|$, which would result in a good approximation
of the spatial distance (especially if the subject is moving at constant
speed). More generally, the elapsed separation across time will reflect
dependence better than the spatial distance. The faster an individual is
moving from one point to the other, the shorter the time elapsed, and higher
the correlation between the two measurements. Hence, we model our measurements
as $Y_{k}(\cdot)\equiv
Z_{k}\circ\gamma_{k}(\cdot):\mathbb{R}^{+}\rightarrow\mathbb{R}$, which, by
construction, is a valid stochastic process.

(a)

(b)

Figure 3: Example of observed points (a) and trajectory (b): black dots are
realizations, grey line is domain of the process

This will form the edifice of the model in Section 3.1, where we are modeling
the dependence by solely considering stochastic evolution through time. How
should spatial information be introduced in the model? Two individuals at the
same spatial coordinate experience the same spatial effect but different
temporal effects because their physical activities are a function of their
trajectory’s temporal evolution. An added complication is that trajectories
intersect and overlap and, in practice, can have multiple observations at the
same location. Even more flexible spatiotemporal covariance kernels (e.g.,
nonseparable or nonstationarity) will struggle to recognize the above
features. Hence, we introduce the spatial effect in the mean using spline
regression (see Section 3.4).

###  3.1 Temporal model

Let $\mathcal{T}=\cup_{k=1}^{K}\mathcal{T}_{k}$ where
$\mathcal{T}_{k}=\left\\{t_{ki}\right\\}_{i=1}^{T_{k}}$ and
$t_{ki}\in\mathbb{R}^{+}$ be the set of the $n=\sum_{k=1}^{K}T_{K}$ observed
time points. We model $\boldsymbol{Y}(\mathcal{T})$ as the finite realization
of a $K$-variate process $\boldsymbol{Y}(\cdot)$ over $\mathbb{R}^{+}$:

$\boldsymbol{Y}(t)=\boldsymbol{X}(t,\boldsymbol{\gamma}(t))^{\top}\bbeta+\boldsymbol{w}(t)+\boldsymbol{\varepsilon}(t),\quad
t\in\mathbb{R}^{+},$ (3)

where $\boldsymbol{Y}(t)=(Y_{1}(t),Y_{2}(t),\ldots,Y_{K}(t))^{\top}$ is a
$K\times 1$ vector of measurements at time $t$ on the $K$ individuals,
$\boldsymbol{X}(t,\boldsymbol{\gamma}(t))$ is a $p\times K$ matrix, each row
being the values of a covariate for the $K$ individuals,
$\boldsymbol{w}(t)=(w_{1}(t),w_{2}(t),\ldots,w_{K}(t))^{\top}$ is a $K\times
1$ vector comprising a temporal process for each individual, and
$\boldsymbol{\varepsilon}(t)\sim\mathcal{N}_{K}(0,\tau^{2}\boldsymbol{I}_{K}),\;\tau^{2}\in\mathbb{R}^{+},$
is a white noise process for measurement error. Each element of
$\boldsymbol{w}(t)$ is specified as $w_{k}(t)\stackrel{{\scriptstyle
ind}}{{\sim}}\mathcal{GP}\left(0,c_{\boldsymbol{\theta}}(\cdot,\cdot)\right)$,
where $c_{\boldsymbol{\theta}}(\cdot,\cdot)$ is a covariance function with
parameters $\boldsymbol{\theta}\in\Theta$.

Let $y_{ki}$ and $\boldsymbol{x}_{ki}$ be the outcome and covariates for
individual $k$ at time point $t_{ki}$, respectively, so
$\left\\{\left(y_{ki},\boldsymbol{x}_{ki}\right):\;k=1,\dots,K,\;i=1,\dots
T_{k}\right\\}$ is the observed data. Let $\boldsymbol{y}_{k}$ and
$\boldsymbol{w}_{k}$ be $T_{k}\times 1$ vectors comprising all measurements
and random effects on patient $k$, respectively. Forming the $n\times 1$
vectors
$\boldsymbol{y}=\begin{bmatrix}\boldsymbol{y}_{1\cdot}^{\top}&\boldsymbol{y}_{2\cdot}^{\top}&\cdots&\boldsymbol{y}_{K\cdot}^{\top}\end{bmatrix}^{\top}$
and
$\boldsymbol{w}=\begin{bmatrix}\boldsymbol{w}_{1\cdot}^{\top}&\boldsymbol{w}_{2\cdot}^{\top}&\cdots&\boldsymbol{w}_{K\cdot}^{\top}\end{bmatrix}^{\top}$,
and the $n\times p$ matrix
$\boldsymbol{X}=\begin{bmatrix}\boldsymbol{X}_{1\cdot}^{\top}&\boldsymbol{X}_{2\cdot}^{\top}&\ldots&\boldsymbol{X}_{K\cdot}^{\top}\end{bmatrix}^{\top}$,
where $\boldsymbol{X}_{k}$ is the $T_{k}\times p$ matrix of predictors
corresponding to $\boldsymbol{y}_{k}$, we extend (3) to a hierarchical model
with posterior distribution

$p(\bbeta,\boldsymbol{w},\boldsymbol{\theta},\tau^{2}\,|\,\boldsymbol{y})\propto
p(\boldsymbol{\theta},\tau^{2})\times
N(\bbeta\,|\,\boldsymbol{\mu}_{\beta},\boldsymbol{V}_{\beta})\times
N(\boldsymbol{w}\,|\,\mathbf{0},\boldsymbol{C}_{\boldsymbol{\theta}})\times
N(\boldsymbol{y}\,|\,\boldsymbol{X}\bbeta+\boldsymbol{w},\tau^{2}\boldsymbol{I}_{n})\;.$
(4)

The covariance matrix
$\boldsymbol{C}_{\boldsymbol{\theta}}=\mbox{diag}\left(\boldsymbol{C}_{\boldsymbol{\theta}_{1},1},\boldsymbol{C}_{\boldsymbol{\theta}_{2},2},\ldots,\boldsymbol{C}_{\boldsymbol{\theta}_{K},K}\right)$
is $n\times n$ block-diagonal with
$\boldsymbol{C}_{\boldsymbol{\theta}_{k},k}=[c_{\boldsymbol{\theta}}(t_{ki},t_{kj})]$
as the $T_{k}\times T_{k}$ temporal covariance matrix corresponding to
individual $k$. Each individual is allowed its own covariance parameters,
$\boldsymbol{\theta}_{k}$, and
$\boldsymbol{\theta}=\\{\boldsymbol{\theta}_{1},\boldsymbol{\theta}_{2},\ldots,\boldsymbol{\theta}_{K}\\}$
in (4) is the collection of all the covariance kernel parameters. Applying (4)
involves the determinant and inverse of
$\boldsymbol{C}_{\boldsymbol{\theta}}$, which require $\mathcal{O}(n^{2})$
storage space and $\mathcal{O}(n^{3})$ floating point operations (flops). The
block-diagonal structure of $\boldsymbol{C}_{\boldsymbol{\theta}}$
considerably alleviates this burden since
$\mbox{det}(\boldsymbol{C})=\prod_{k=1}^{K}\mbox{det}(\boldsymbol{C}_{\boldsymbol{\theta}_{k},k})$
and
$\boldsymbol{C}^{-1}=\mbox{diag}\left(\boldsymbol{C}_{\boldsymbol{\theta}_{1},1}^{-1},\boldsymbol{C}_{\boldsymbol{\theta}_{2},2}^{-1},\ldots,\boldsymbol{C}_{\boldsymbol{\theta}_{K},K}^{-1}\right)$.
This reduces the flop count from
$\mathcal{O}(n^{3})=\mathcal{O}((\sum_{k=1}^{K}T_{k})^{3})$ to
$\mathcal{O}(K\sum_{k=1}^{K}(T_{k})^{3})$, with a significant saving of
calculations especially when the $T_{k}$’s are reasonably small ($<10^{4}$).
Furthermore, each $\boldsymbol{C}_{\boldsymbol{\theta}_{k},k}$ can be computed
in parallel rendering further further scalability to the algorithm.

However, analyzing the Actigraph data in Section 2 will involve $T_{k}>10^{5}$
measurements from some individuals. Full inference will be impractical without
any exploitable structure for each
$\boldsymbol{C}_{\boldsymbol{\theta}_{k},k}$. Analyzing massive spatiotemporal
data has witnessed burgeoning interest and a comprehensive review is beyond
the scope of this work (see, e.g., Heaton et al., 2019, and references
therein). We will pursue an approximation due to Vecchia (Vecchia, 1988) that
has generated substantial recent interest (Datta et al., 2016b, a; Katzfuss et
al., 2020; Katzfuss and Guinness, 2021; Peruzzi et al., 2020) in scalable
Bayesian modeling.

###  3.2 Independent DAG models over individuals

We adapt Vecchia’s likelihood approximation (Vecchia, 1988) to the random
effects $\boldsymbol{w}_{k}$ for each $k=1,2,\ldots,K$. Beginning with the
observed time points $\\{t_{k1}<t_{k2}<\cdots<t_{kT_{k}}\\}$ for individual
$k$ and the directed acyclic graphical (DAG) representation
$p(\boldsymbol{w}_{k})=p(w_{k1})\prod_{i=2}^{T_{k}}p(w_{ki}\,|\,w_{k1},\ldots,w_{k(i-1)})$,
we define

$p(\boldsymbol{w}_{k})\approx\widetilde{p}(\boldsymbol{w}_{k})=p(w_{k1})\prod_{i=2}^{T_{k}}p(w_{ki}|\boldsymbol{w}_{k,N(i)})\;,$
(5)

where $\widetilde{p}(\cdot)$ is the joint density derived from
$p(\boldsymbol{w}_{k})$ by restricting the parents (conditional sets) of each
$w_{ki}$ in the DAG to a set $w_{kN(i)}=\\{w_{kj}:j\in N(i)\\}$, where $N(i)$
is a set of prefixed size $m$ comprising the $m$ nearest neighbors of $t_{ki}$
from the past. Thus, $N(i)=\\{t_{k(i-m)},\ldots,<t_{k(i-1)}\\}$ for $i>m$ and
$N(i)=\\{t_{k1},\ldots,t_{k(i-1)}\\}$ for $i\leq m$. Such approximations yield
valid probability likelihoods (Lauritzen, 1996; Stein et al., 2004; Murphy,
2012) and can be extended to stochastic processes (Datta et al., 2016b) for
inference on arbitrary time points.

The connection between sparsity and conditional independence follows by
writing (5) as a linear model
$\boldsymbol{w}_{k}=\boldsymbol{A}_{k}\boldsymbol{w}_{k}+\boldsymbol{\eta}_{k}$,
where $\boldsymbol{A}_{k}$ is a $T_{k}\times T_{k}$ strictly lower triangular
matrix,
$\boldsymbol{\eta}_{k}\sim\mathcal{N}_{T_{k}}(\mathbf{0},\boldsymbol{D}_{k})$
and $\boldsymbol{D}_{k}$ is the $T_{k}\times T_{k}$ diagonal matrix such that
$\left[\boldsymbol{D}_{k}\right]_{ii}=d_{ii}=\text{Var}\left(w_{ki}|\\{w_{kj},j<i\\}\right)$
for $i=1,\dots,T_{k}$. The DAG imposes the lower-triangular structure on
$\boldsymbol{A}_{k}$ and its $(i,j)$-th entry is allowed to be nonzero only
for $j\in N(i)$. Therefore, each row of $\boldsymbol{A}_{k}$ has at most $m$
nonzero entries so that
$\widetilde{\boldsymbol{C}}_{k}^{-1}=(\boldsymbol{I}_{T_{k}}-\boldsymbol{A}_{k})^{\top}\boldsymbol{D}_{k}^{-1}(\boldsymbol{I}_{T_{k}}-\boldsymbol{A}_{k})$
is sparse, where $\widetilde{\boldsymbol{C}}_{k}^{-1}$ is the precision matrix
corresponding to $\widetilde{p}(\boldsymbol{w}_{k})$. Replacing
$\boldsymbol{C}$ with $\widetilde{\boldsymbol{C}}$ in (4) yields a
computationally efficient hierarchical model with
$N(\prod_{k=1}^{K}N(\boldsymbol{w}_{k}\,|\,\mathbf{0},\widetilde{\boldsymbol{C}}_{k})$
as the prior on $\boldsymbol{w}$.

The key observation is that the nonzero elements of the $i$-th row of
$\boldsymbol{A}_{k}$ is the solution $\boldsymbol{a}_{k}$ of the $m\times m$
linear system
$\boldsymbol{C}_{\boldsymbol{\theta},k}[N(i),N(i)]\boldsymbol{a}_{k}=\boldsymbol{C}_{\boldsymbol{\theta},k}[N(i),i]$,
where $[\cdot,\cdot]$ indicates submatrices defined by the given row and
column index sets. Obtaining the nonzero elements of $\boldsymbol{A}_{k}$ and
$\boldsymbol{D}_{k}$ costs $\mathcal{O}(T_{k}m^{3})$ (scales linearly with
$T_{k}$) instead of $\mathcal{O}(T_{k}^{3})$ as would have been without
sparsity. This cheaply delivers the quadratic form
$\boldsymbol{w}_{k}^{\top}\widetilde{\boldsymbol{C}}_{k}^{-1}\boldsymbol{w}_{k}$
in terms of $\boldsymbol{A}_{k}$ and $\boldsymbol{D}_{k}$ and the determinant
$\text{det}(\widetilde{\boldsymbol{C}}_{k})=\prod_{i=1}^{T_{k}}d_{ii}$ at
almost no additional cost. The lower triangular matrix $\boldsymbol{A}_{k}$ is
not just sparse but also banded, with a lower bandwidth equal to $m$.
Consequently, $\widetilde{\boldsymbol{C}}_{k}^{-1}$ is also banded with lower
and upper bandwidth equal to $m$. This leads to further accrual of
computational benefits. The overall cost is
$\mathcal{O}(\sum_{k=1}^{K}T_{k}m^{3})=\mathcal{O}(nm^{3})$ (linear in $n$)
for computing the posterior for any given values of the parameters.

###  3.3 Implementation using collapsed models

The Bayesian hierarchical model in (4), either with
$\boldsymbol{C}_{\boldsymbol{\theta}}$ or with
$\widetilde{\boldsymbol{C}}_{\boldsymbol{\theta}}$ in the prior for
$\boldsymbol{w}$, allows full posterior inference for
$\\{\bbeta,\boldsymbol{w},\boldsymbol{\theta},\tau^{2}\\}$ using Markov chain
Monte Carlo (MCMC). Gibbs sampling with random walk Metropolis steps provide
full conditional distributions in closed form for
$\\{\bbeta,\boldsymbol{w}\\}$ and also for $\tau^{2}$ with an
$\mathcal{IG}(a_{\tau},b_{\tau})$ prior. However, this convenience is
nullified in practice by strong autocorrelation and poor mixing of the chains
(Liu et al., 1994). Samplers based on spatial DAG-based models have been
devised, explored and compared in Finley et al. (2019). Instead of (4), we
sample from

$p(\bbeta,\boldsymbol{\theta},\tau^{2}\,|\,\boldsymbol{y})\propto
p(\boldsymbol{\theta},\tau^{2})\times
N(\bbeta\,|\,\boldsymbol{\mu}_{\beta},\boldsymbol{V}_{\beta})\times
N(\boldsymbol{y}\,|\,\boldsymbol{X}\bbeta,\widetilde{\boldsymbol{C}}_{\boldsymbol{\theta}}+\tau^{2}\boldsymbol{I}_{n})\;,$
(6)

which is derived from (4) by integrating out $\boldsymbol{w}$, thereby
“collapsing” the parameter space to a much smaller domain without
$\boldsymbol{w}$. This considerably improves mixing and convergence.

We will need to compute the inverse and determinant of
$\widetilde{\boldsymbol{\Lambda}}=\widetilde{\boldsymbol{C}}_{\boldsymbol{\theta}}+\tau^{2}\boldsymbol{I}_{n}$,
which is $n\times n$. While $\widetilde{\boldsymbol{\Lambda}}^{-1}$ does not
share the same convenient factorization of $\widetilde{\boldsymbol{C}}^{-1}$
and is also not guaranteed to be sparse, the Sherman-Woodbury-Morrison
formulas reveal

$\widetilde{\boldsymbol{\Lambda}}^{-1}=\tau^{-2}\boldsymbol{I}-\tau^{-4}\boldsymbol{\Omega}^{-1},\quad\text{with}\quad\boldsymbol{\Omega}=\widetilde{\boldsymbol{C}}^{-1}+\tau^{-2}\boldsymbol{I}\;,$
(7)

where $\boldsymbol{\Omega}$ enjoys the same sparsity as $\boldsymbol{C}^{-1}$.
Moreover,
$\text{det}(\widetilde{\boldsymbol{\Lambda}})=\tau^{2n}\text{det}(\widetilde{\boldsymbol{C}})\text{det}(\boldsymbol{\Omega})$.
The core of the algorithm is therefore to compute
$\widetilde{\boldsymbol{\Lambda}}^{-1}$ through $\boldsymbol{\Omega}$. In our
application, the random effect is assumed to be the realization of $K$
independent temporal processes. As discussed in Section 3.2, this implies a
block-diagonal structure for $\widetilde{\boldsymbol{C}}$ that can be shown to
be shared also by $\boldsymbol{\Omega}$ (see Eq. (7)). Each block
$\boldsymbol{\Omega}_{k}$ of $\boldsymbol{\Omega}$ can be computed
independently for each individual and the same holds for its inverse and its
determinant. This means that the body of the algorithm will consist of a loop
over all the individuals, which allows for straightforward parallelization;
see Algorithm 1. Unlike in spatial DAGs (Datta et al., 2016b; Finley et al.,
2019), we do not need fill-reducing permutation methods since neighbors sets
for temporal processes consist of contiguous observations and
$\left\\{\boldsymbol{\Omega}_{k}\right\\}_{k=1}^{K}$ are banded matrices with
no gaps.

0: Initialization

begin

for  _$k=1,\dots,K$_ do

a: Compute $d^{k}_{ij}=|t_{j}-t_{i}|,\quad\forall
t_{j},t_{i}\in\mathcal{T}_{k}$

b: Find the neighbor sets $\\{N_{k}(i)\\}_{i=1}^{T_{k}}$

end for

end

1: Metropolis-Hastings update for $\\{\boldsymbol{\theta},\tau^{2}\\}$

$p(\boldsymbol{\theta},\tau^{2}|\cdot)\propto
p(\boldsymbol{\theta},\tau^{2})\times\frac{1}{\sqrt{\det{\widetilde{\boldsymbol{\Lambda}}}}}\exp\left(-\frac{1}{2}\left(\boldsymbol{y}-\boldsymbol{X}\bbeta\right)^{\top}\widetilde{\boldsymbol{\Lambda}}^{-1}\left(\boldsymbol{y}-\boldsymbol{X}\bbeta\right)\right)$

begin

for  _$k=1,\dots,K$_ do

a: Compute
$\boldsymbol{L}_{k}=\left(\boldsymbol{I}_{T_{k}}-\boldsymbol{A}_{k}\right)^{\top}$,
$\boldsymbol{d}_{k}=\text{diag}(\boldsymbol{D}_{k})$ and
$\boldsymbol{R}_{k}=\boldsymbol{D}_{k}^{-1}\left(\boldsymbol{I}_{T_{k}}-\boldsymbol{A}_{k}\right)$
using $\boldsymbol{C}_{k}$ and $\\{N_{k}(i)\\}_{i=1}^{T_{k}}$

b: Compute
$\boldsymbol{\Omega}_{k}=\boldsymbol{L}_{k}\cdot\boldsymbol{R}_{k}+\tau^{-2}\boldsymbol{I}_{T_{k}}$
exploiting sparsity

c: Compute
$\boldsymbol{r}_{k}=\boldsymbol{y}_{k\cdot}-\boldsymbol{X}_{k\cdot}\bbeta$ and
$\delta_{\boldsymbol{D}_{k}}=\prod_{i=1}^{T_{k}}d_{k,i}$

d: Compute
$\boldsymbol{v}_{k}=\boldsymbol{\Omega}_{k}^{-1}\boldsymbol{r}_{k}$,
$\boldsymbol{u}_{k}=\boldsymbol{\Omega}_{k}^{-1}\boldsymbol{X}_{k\cdot}$ and
$\delta_{\boldsymbol{\Omega}_{k}}=\text{det}(\boldsymbol{\Omega}_{k})$
exploiting the sparse Cholesky decomposition of $\boldsymbol{\Omega}_{k}$

e: Collect $\boldsymbol{r}_{k}$, $\boldsymbol{v}_{k}$ and $\boldsymbol{u}_{k}$
into $\boldsymbol{r}$, $\boldsymbol{v}$ and $\boldsymbol{u}$, respectively.

end for

f: Compute
$q_{1}=\tau^{2n}\cdot\prod_{k=1}^{K}\delta_{\boldsymbol{D}_{k}}\cdot\prod_{k=1}^{K}\delta_{\boldsymbol{\Omega}_{k}}$
and
$q_{2}=\boldsymbol{r}^{\top}\boldsymbol{r}/\tau^{2}-\boldsymbol{r}^{\top}\boldsymbol{v}/\tau^{4}$

g: Get
$p(\boldsymbol{\theta},\tau^{2}|\cdot)\propto\frac{\exp(-q_{2}/2)}{\sqrt{q_{1}}}\cdotp(\boldsymbol{\theta},\tau^{2})$

end

2: Gibbs’ sampler update for $\bbeta$

$\bbeta|\cdot\sim\mathcal{N}_{p}(\boldsymbol{B}^{-1}\boldsymbol{b},\boldsymbol{B}^{-1})$,
where
$\boldsymbol{B}=\boldsymbol{X}^{\top}\widetilde{\boldsymbol{\Lambda}}^{-1}\boldsymbol{X}+\boldsymbol{V}_{\bbeta}^{-1}$
and
$\boldsymbol{b}=\boldsymbol{X}^{\top}\widetilde{\boldsymbol{\Lambda}}^{-1}\boldsymbol{y}+\boldsymbol{V}_{\bbeta}^{-1}\boldsymbol{\mu}_{\bbeta}$

begin

a: Compute $\boldsymbol{F}=\boldsymbol{V}^{-1}_{\bbeta}$ and
$\bf=\boldsymbol{V}^{-1}_{\bbeta}\boldsymbol{\mu}_{\bbeta}$

b: Compute
$\boldsymbol{b}=\boldsymbol{y}^{\top}\boldsymbol{X}/\tau^{2}-\boldsymbol{y}^{\top}\boldsymbol{v}/\tau^{4}+\bf$
and
$\boldsymbol{B}=\boldsymbol{X}^{\top}\boldsymbol{X}/\tau^{2}-\boldsymbol{X}^{\top}\boldsymbol{v}/\tau^{4}+\boldsymbol{F}$

c: Generate
$\bbeta\sim\mathcal{N}_{p}(\boldsymbol{B}^{-1}\boldsymbol{b},\boldsymbol{B}^{-1})$

end

Repeat steps 1 and 2 to obtain $M$ MCMC samples for
$\\{\bbeta,\boldsymbol{\theta},\tau^{2}\\}$

Algorithm 1 Sampling from the posterior of the collapsed temporal model

We devised a Gibbs sampler with Metropolis random walk updates for (6), where
$\bbeta$ is updated from its full conditional distribution, while
$\\{\boldsymbol{\theta},\tau^{2}\\}$ are updated using an adaptive Metropolis
step based on Haario et al. (2001). Here, after the first few iterations, a
new proposal covariance matrix is regularly computed on the run according to
the empirical covariance of the current chain. Subsequently, a mixture of the
original and adaptive proposal is used as the new proposal. Convergence toward
the desired acceptance rate is assured for an appropriate choice of the
variance terms and of the adaptation rule (Roberts and Rosenthal, 2009). The
algorithm has been coded using the R 4.0.1 statistical environment. All
expensive computations are managed by the Eigen library (version 3.3.7), which
provides efficient routines for numerical linear algebra with an emphasis on
sparse matrices. Our implementation of (6) outperforms the algorithms that
update $\boldsymbol{w}$ in terms of computational speed as it is implemented
in the spNNGP package (Finley et al., 2017). We present these comparisons in
the online supplement (Appendix).

###  3.4 Including spatial effects

Accounting for spatial information in our Actigraph dataset presents some new
considerations. As mentioned in Section 1, spatial information is available to
us in terms of the physical location along the trajectory as well as through
covariates that are functions of space. Considering the discussion in Section
3, the analytical goals of this dataset suggest accounting for spatial
heterogeneity. Here, as argued earlier, modeling $\boldsymbol{w}(\cdot)$ in
(3) as a spatiotemporal process, including scalable versions, has challenges
given that: (i) the trajectory’s domain does not have a positive area; and
(ii) associations among the measurements are more amenable to the temporal
scale. Therefore, we introduce spatial effects into the mean employing a
smooth function of space, $f_{S}(\cdot):\mathbb{R}^{2}\rightarrow\mathbb{R}$,
approximated by a spline basis representation (see, e.g., Goodman and Hardin,
2006; Ramsay and Silverman, 2007). For instance, if $J_{x}$ and $J_{y}$ are
the dimensions of independently defined B-spline basis expansions on the $x$
and $y$ coordinates, respectively, then
$f_{S}\left((x,y)\right)\approx\widetilde{f}_{S}\left((x,y)\right)=\sum_{j_{X}=1}^{J_{X}}\sum_{j_{Y}=1}^{J_{Y}}\beta_{S,(j_{X},j_{Y})}B_{x,j_{x}}(x)B_{y,j_{y}}(y)$,
where $B_{x,j_{X}}=\left[\boldsymbol{B}_{x}\right]_{j_{x}}$ and
$B_{y,j_{Y}}=\left[\boldsymbol{B}_{y}\right]_{j_{y}}$ are the $j_{x}$-th and
$j_{y}$-th element of the B-spline basis along the two axis. For any location
$(x,y)\in\mathbb{R}^{2}$ the elements of the previous sum can be more
compactly expressed through the tensor product basis
$\boldsymbol{B}_{S}(x,y)=\left(\boldsymbol{B}_{x}\otimes\boldsymbol{B}_{y}\right)(x,y)$.
The size of this basis is $J_{S}=J_{x}\cdot J_{y}$ and depends on the size of
the two original spline basis, which in turn depends on the chosen number of
knots ${knots}_{x},{knots}_{y}$ and degree ${deg}_{x},{deg}_{y}$ (namely
$J_{c}={knots}_{c}+{deg}_{c}$ for $c=x,y$). We now modify (3) to include the
spline,

$\boldsymbol{Y}(t)=\boldsymbol{X}(t)\bbeta+\boldsymbol{B}_{S}\left(\boldsymbol{\gamma}(t)\right)\bbeta_{S}+\boldsymbol{w}(t)+\boldsymbol{\varepsilon}(t),\qquad
t\in\mathcal{\mathbb{R}^{+}}\;,$ (8)

where
$\boldsymbol{\gamma}(t)=\left\\{\boldsymbol{\gamma}_{1}(t),\boldsymbol{\gamma}_{2}(t),\ldots,\boldsymbol{\gamma}_{K}(t)\right\\}$,
$\boldsymbol{\gamma}_{k}(\cdot)=(\gamma_{k,x}(t),\gamma_{k,y}(t)):\mathbb{R}^{+}\rightarrow\mathbb{R}^{2}$
is the trajectory function mapping time $t$ for individual $k$ to its position
and $\boldsymbol{B}_{S}\left(\boldsymbol{\gamma}(t)\right)$ is the $K\times
J_{S}$ matrix with row $k$ corresponding to the $J_{S}$ basis elements for the
coordinates at time point $t$ for individual $k$. A proper choice of $J_{S}$
(i.e. knots and degree) is required to fit a spline surface flexible enough to
describe the spatial variations at the scale of interest without incurring
over-fitting. Let us denote with
$\boldsymbol{B}=\boldsymbol{B}_{S}\left(\boldsymbol{\gamma}({\cal T})\right)$
the $n\times J_{S}$ matrix containing the B-spline basis elements evaluated at
the observed location of each individual $\boldsymbol{\gamma}({\cal
T})=\left\\{\gamma_{1}(t_{11}),\gamma_{1}(t_{12}),\dots,\gamma_{K}(t_{KT_{k}})\right\\}$.
Following Equation (6), we sample from the posterior,

$p(\bbeta,\bbeta_{S},\boldsymbol{\theta},\tau^{2}\,|\,\boldsymbol{y})\propto
p(\boldsymbol{\theta},\tau^{2})\times p_{S}(\bbeta_{S})\times
N(\bbeta\,|\,\boldsymbol{\mu}_{\beta},\boldsymbol{V}_{\beta})\times
N(\boldsymbol{y}\,|\,\boldsymbol{X}\bbeta+\boldsymbol{B}\bbeta_{S},\widetilde{\boldsymbol{C}}_{\boldsymbol{\theta}}+\tau^{2}\boldsymbol{I}_{n})\;,$
(9)

where the prior $p_{S}(\cdot)$ shall be accurately chosen. We must consider
that our Actigraph data includes millions of observations in a limited study
area, of which some assume different values in the same location (or in its
immediate vicinity) so over-fitting will not be an issue. However, some areas
present sparsely observed points (trajectories are not uniformly distributed,
as shown in Figure 2). This may cause coefficients corresponding to those
regions to be weakly identified. To control for the balance of all these
components, we may assign ad-hoc priors to the spatial spline regression
coefficients (Eilers and Marx, 1996) for penalizing deviation from a certain
degree of smoothness and favoring identifiability. This yields the Bayesian
P-Spline (Hastie et al., 2000; Lang and Brezger, 2004). While keeping the
Gaussian priors, we effectuate shrinkage by choosing a suitable precision
matrix $\boldsymbol{P}$ and introducing a shrinkage parameter $\lambda$ at a
deeper level of the hierarchy. To be precise,
$\displaystyle\bbeta_{S}\,|\,\lambda\propto\exp\left\\{-\frac{\lambda}{2}\cdot\bbeta_{S}\boldsymbol{P}\bbeta_{S}^{\top}\right\\}$
and $\lambda\sim\mathcal{G}(\alpha_{\lambda},\beta_{\lambda})$. We consider
two possible forms for $\boldsymbol{P}$, which imply different penalization
for the coefficients:

  * •

Ridge-like prior, which is to say
$\boldsymbol{P}=\boldsymbol{P}_{RL}=\boldsymbol{I}_{J_{S}}$;

  * •

First-order random walk prior, which is to say:

$\boldsymbol{P}=\boldsymbol{P}_{RW}:\;[\boldsymbol{P}_{RW}]_{ij}=\begin{cases}n_{i}\quad&i=j\\\
-1\quad&i\sim j\\\ 0\quad&\text{otherwise}\end{cases}$

where $n_{i}$ is the number of neighbors of knot $i$ and $i\sim j$ denotes a
neighboring relationship between the knots.

1: Gibbs’ sampler update for $\boldsymbol{\psi}$

$\boldsymbol{\psi}|\cdot\sim\mathcal{N}_{J}(\boldsymbol{G}^{-1}\boldsymbol{g},\boldsymbol{G}^{-1})$,
where
$\boldsymbol{G}={\boldsymbol{X}^{*}}^{\top}\widetilde{\boldsymbol{\Lambda}}^{-1}{\boldsymbol{X}^{*}}+\boldsymbol{V}_{\boldsymbol{\psi}}^{-1}$
and
$\boldsymbol{g}={\boldsymbol{X}^{*}}^{\top}\widetilde{\boldsymbol{\Lambda}}^{-1}\boldsymbol{y}+\boldsymbol{V}_{\boldsymbol{\psi}}^{-1}\boldsymbol{\mu}_{\boldsymbol{\psi}}$

begin

a: Compute $\boldsymbol{F}=\boldsymbol{V}^{-1}_{\boldsymbol{\psi}}$ and
$\bf=\boldsymbol{V}^{-1}_{\boldsymbol{\psi}}\boldsymbol{\mu}_{\boldsymbol{\psi}}$

b: Compute
$\boldsymbol{g}=\boldsymbol{y}^{\top}{\boldsymbol{X}^{*}}/\tau^{2}-\boldsymbol{y}^{\top}\boldsymbol{v}/\tau^{4}+\bf$
and
$\boldsymbol{G}={\boldsymbol{X}^{*}}^{\top}{\boldsymbol{X}^{*}}/\tau^{2}-{\boldsymbol{X}^{*}}^{\top}\boldsymbol{v}/\tau^{4}+\boldsymbol{F}$

c: Generate
$\boldsymbol{\psi}\sim\mathcal{N}_{p^{*}}(\boldsymbol{G}^{-1}\boldsymbol{g},\boldsymbol{G}^{-1})$

end

2: Gibbs’ sampler update for $\lambda$

$\lambda|\cdot\sim Ga\left(\alpha_{\lambda}^{*},\,\beta_{\lambda}^{*}\right)$,
where $\alpha_{\lambda}^{*}=\alpha_{\lambda}+1/2$ and
$\beta_{\lambda}^{*}=\beta_{\lambda}+\bbeta_{S}^{\top}\boldsymbol{P}\bbeta_{S}$

begin

a: Compute $h=\bbeta_{S}^{\top}\boldsymbol{P}\bbeta_{S}$ and get:
$\alpha_{\lambda}^{*}=\alpha_{\lambda}+1/2$ and
$\beta_{\lambda}^{*}=\beta_{\lambda}+h$

b: Generate
$\lambda\sim\mathcal{G}\left(\alpha_{\lambda}^{*},\,\beta_{\lambda}^{*}\right)$

end

Algorithm 2 $\boldsymbol{\psi}$ and $\lambda$ Gibbs’ update in the collapsed
algorithm with shrinkage

Both precision matrices provide a multivariate Gaussian prior distribution on
the coefficients. However, the latter is improper since
$\text{rank}\left(\boldsymbol{P}_{RW}\right)<J_{S}$. Nevertheless, if we
collect the B-Spline basis elements with the other covariates as
$\boldsymbol{X}^{*}=\left[\boldsymbol{X},\boldsymbol{B}\right]$ and stack the
corresponding coefficients into the joint vector
$\boldsymbol{\psi}=\left[\bbeta,\bbeta_{S}\right]$, then the posterior
distribution of the latter is a proper multivariate Gaussian with full
conditional distribution
$\boldsymbol{\psi}\,|\,\cdot\propto\mathcal{N}_{J}\left(\boldsymbol{\psi}\left|\boldsymbol{G}^{-1}\boldsymbol{g},\,\boldsymbol{G}^{-1}\right.\right)$,
where
$\boldsymbol{G}={\boldsymbol{X}^{*}}^{\top}\widetilde{\boldsymbol{\Lambda}}^{-1}{\boldsymbol{X}^{*}}+\boldsymbol{V}_{\boldsymbol{\psi}}^{-1}$
and
$\boldsymbol{g}={\boldsymbol{X}^{*}}^{\top}\widetilde{\boldsymbol{\Lambda}}^{-1}\boldsymbol{y}+\boldsymbol{V}_{\boldsymbol{\psi}}^{-1}\boldsymbol{\mu}_{\boldsymbol{\psi}}$
with
$\boldsymbol{V}_{\boldsymbol{\psi}}^{-1}=\text{diag}\left(\boldsymbol{V}_{\bbeta}^{-1},\,\lambda\cdot\boldsymbol{P}\right)$
and
$\boldsymbol{g}=\left[\boldsymbol{\mu}_{\bbeta},\,\boldsymbol{\mu}_{\bbeta_{S}}\right]^{\top}=\mathbf{0}^{\top}$.
Moreover, the Gamma prior on $\lambda$ implies a Gamma full-conditional
distribution
$\displaystyle\lambda\,|\,\cdot\propto\mathcal{G}\left(\lambda\left|\alpha_{\lambda}+1/2,\,\beta_{\lambda}+\bbeta_{S}^{\top}\boldsymbol{P}\bbeta_{S}\right.\right)$.

Estimating the model in (8) is achieved through a straightforward extension of
Algorithm 1. We jointly update $\boldsymbol{\psi}$ and $\lambda$ from their
full conditional distributions. Algorithm 2 shows how the Gibbs’ sampling step
of Algorithm 1 can be modified to get full inference also on the spline
coefficients $\bbeta_{S}$ and the shrinkage parameter $\lambda$. In practical
terms, this requires $J_{S}$ additional linear coefficients to be estimated,
whose size $p^{*}=p+J_{S}$ may undermine the efficiency of the algorithm. For
example, calculations in Step 1b are quadratic w.r.t.
$p^{*}\rightarrow\mathcal{O}(n{p^{*}}^{2})$. Steps 1a and 1b (i.e. the most
expensive in $p^{*}$) are executed in the first iteration and subsequently,
only in those iterations where new values of $\boldsymbol{\theta}$ are
accepted. When $\boldsymbol{\theta}$ is rejected, we retain in memory the
previously computed value (which would stay unchanged). Thus, if we attain an
optimal acceptance rate of $\approx 20\%-30\%$ in the Metropolis Hastings step
on $\boldsymbol{\theta}$, the computation is avoided in the majority of cases
with a sensible improvement in computation time and speed.

###  3.5 Simulations

We conducted simulation experiments to evaluate the model described in Section
3.4 and compared the performance of our algorithm in terms of fitting,
prediction error and computational speed with other routines available from
the spNNGP package. Additional comparative experiments are provided in the
supplementary material. We executed our MCMC algorithms on a computing
environment equipped with 12 modern computational nodes with 16 cores each
(bringing the overall number of cores to 192), roughly equivalent to 3
TeraFlop/sec, and $64$Gb of RAM. Each of the following jobs, and the ones from
Section 4, have been executed on a single node exploiting all 16 cores. The
results presented here and in Section 4 are based upon posterior samples that
were retained after diagnosing convergence using visual tools (e.g.,
traceplots, autocorrelation), effective sample sizes, Monte Carlo standard
errors (MCSE) and other diagnostics offered by the coda, mcse and bayesplot
packages in the R computing environment; the online supplement (Appendix)
includes specific details.

We first generated $T_{k}=2\times 10^{5}$ time points for $K=5$ individuals,
where each time point $t_{ki}$ followed exponential waiting times between
observations, i.e. $t_{ki}=\sum_{h=1}^{i-1}\delta_{h}$, and
$\delta_{h}\stackrel{{\scriptstyle iid}}{{\sim}}\mbox{Exp}(5)$. Given the time
points, we constructed spatial trajectories $\boldsymbol{\gamma}_{k}(\cdot)$,
$k=1,\ldots,K$, by simulating
$\boldsymbol{s}_{k}=\left[\gamma_{k}(t_{k1}),\dots,\gamma_{k}(t_{kT_{k}})\right]^{\top}$,
where subsequent components were independent Gaussian random walks over the
square $\mathcal{S}=(1,10)\times(1,10)$, with the variance of each step along
the horizontal and vertical axis proportional to the elapsed time between two
subsequent observations. If the trajectory left the square, it was projected
onto the border and the next step would resume from there. The simulated
trajectories are shown in Figure 4(a).

(a) Trajectories

(b) Points

Figure 4: Observed trajectories (a) and observed points (b) for the simulated dataset. Param. (True) | S-Spline | P-Spline
---|---|---
Point | Interval | Point | Interval
$\beta_{01}\;(-3.76)$ | -3.799 | (-3.846,-3.752) | -3.797 | (-3.844,-3.75)
$\beta_{02}\;(0.65)$ | 0.572 | (0.523,0.62) | 0.575 | (0.526,0.623)
$\beta_{03}\;(-0.60)$ | -0.649 | (-0.697,-0.6) | -0.646 | (-0.693,-0.598)
$\beta_{04}\;(2.36)$ | 2.326 | (2.277,2.374) | 2.328 | (2.28,2.376)
$\beta_{05}\;(-0.33)$ | -0.359 | (-0.408,-0.31) | -0.356 | (-0.404,-0.308)
$\beta_{1}\;(2.59)$ | 2.599 | (2.59,2.608) | 2.599 | (2.59,2.608)
$\beta_{2}\;(2.70)$ | 2.691 | (2.683,2.7) | 2.691 | (2.683,2.7)
$\beta_{3}\;(-0.58)$ | -0.586 | (-0.595,-0.577) | -0.586 | (-0.595,-0.577)
$\sigma^{2}\;(1)$ | 1.001 | (0.973,1.032) | 0.993 | (0.965,1.023)
$\phi\;(1)$ | 0.994 | (0.948,1.04) | 1.01 | (0.964,1.063)
$\tau^{2}\;(1)$ | 1.001 | (0.984,1.018) | 1.001 | (0.984,1.018)
Metric | Out-of-sample | In-sample | Out-of-sample | In-sample
Coverage | 0.95 | 0.99 | 0.95 | 0.99
RMSPE (r) | 0.07 (1.18) | 0.03 (0.84) | 0.07 (1.19) | 0.03 (0.84)
PIW | 4.66 | 4.44 | 4.66 | 4.44
DIC | 115’543 | 115’556
Fitting time (h) | 2.18 | 2.2
Table 2: Parameter estimates, predictive validation and fitting times (hours)
on the simulated dataset for all the considered models.

Given the time points and positions, we generated the latent temporal Gaussian
processes $w_{k}(\cdot)\stackrel{{\scriptstyle
ind}}{{\sim}}\mathcal{GP}(0,c_{\boldsymbol{\theta}}(\cdot,\cdot))$ with an
exponential covariance $\displaystyle
c_{\boldsymbol{\theta}}(t,\,t^{\prime})=\sigma^{2}\exp\\{-\phi\cdot|t-t^{\prime}|\\}$,
where $\sigma^{2}>0$ represents the variance of the process, $\phi>0$ is the
decay in temporal correlation (range) and $\tau^{2}>0$ the residual variance
(nugget). The spatial effects are then introduced through
$f_{S}(\cdot):\mathcal{S}\rightarrow\mathbb{R}$ by considering a tensor
product spline basis of degree $2$ and with $9$ knots over the square domain
(including boundary knots), where the spline coefficients $\bbeta_{S}$ have
been fixed to randomly generated values from
$\mathcal{N}_{81}\left(\mathbf{0},\lambda\boldsymbol{I}_{81}\right)$ with
$\lambda=0.5$. The model also included individual-specific intercepts
$\left\\{\beta_{0k}\right\\}_{k=1}^{5}$ and the effect of 3 covariates with
random values drawn independently at each location from a $\mathcal{N}(0,1)$
distribution, leading to covariate vectors
$\left\\{\boldsymbol{x}_{ki}\right\\}_{i=1}^{T_{k}},\;k=1,\dots,K$. The effect
of the covariates is assumed common across individuals, and set to be
determined by slopes $\bbeta=[\beta_{1},\beta_{2},\beta_{3}]^{\top}$.

We generated values of the outcome for individual $k$ at time $t_{ki}$ and
location $\boldsymbol{s}_{ki}=\gamma_{k}(t_{ki})$ according to the generative
process defined by (8) with parameters fixed as above. This yielded a
simulated dataset
$D_{sim}=\left\\{(\text{Ind}_{j},t_{j},\boldsymbol{s}_{j},y_{j},\boldsymbol{x}_{j}^{\top})\right\\}_{j=1}^{n}$
with $n=10^{5}$ observations, where $\text{Ind}_{j}$ denotes the individual
corresponding to row $j$. Then, we fit the model in (9) on $70\%$ of the total
observations in $D_{sim}$ using Algorithm 1 with the Gibbs’ sampling modified
as in Algorithm 2. The remaining $30\%$ were held out to assess out-of-sample
predictive performances in terms of Relative and Root Mean Squared Prediction
Error (RMSPE), and Coverage, Predictive Interval Width (PIW). Intercept and
slope regression parameters were assigned $\mathcal{N}(0,\,10^{6})$ priors;
the variance components, $\sigma^{2}$ and $\tau^{2}$, were both assigned
inverse Gamma $\mathcal{IG}(2,2)$ priors; and the decay parameter $\phi$
received a Gamma prior $\mathcal{G}(1,1)$. For the spline coefficients, we
considered both the penalized versions in Section 3.4. The first is referred
to as an S-Spline (shrinking splines), and the second as P-Spline (penalized
splines).

Table 2 presents the posterior estimates. We also included the Deviance
Information Criterion (DIC) for both models. Performances in the two settings
are almost identical, but the DIC favors the S-Spline model. This is not
surprising as the data were generated using an analogous shrinkage prior for
the $\bbeta_{S}$’s. Further details, including the estimates of the spline
coefficients are provided in the online supplement (Appendix). Figure 5
presents the posterior estimate of the spatial surface. We compare the true
latent surface with the two (practically identical) estimates.

(a) True splines surface

(b) Estimated S-Splines

(c) Estimated P-Splines

Figure 5: True (top left) and estimated spline surfaces (bottom left and
right), including the point-wise difference between the true one and the
S-Spline estimated (top right).

##  4 Application

We analyze activity levels throughout the “active time”—when the Actigraph
device records the individual as being physically active—excluding epochs when
the device was not worn or when there was no physical activity (e.g., the
individual was sitting or lying down). The data processing and merging of
actigraph data with GPS locations resulted in two final datasets (Section 2).
These are treated separately. In both applications, $70\%$ of the total
observations were used for training the model, while the remaining were
excluded to assess the out-of-sample predictive performances.

###  4.1 Temporal model

We first analyze $D_{1}$. Our predictors include a binary variable indicating
if the measures refer to the period before or after a Bruin Bike Share (BBS)
program was launched in Westwood, Los Angeles to account for the effect of a
new specific policy which aims at improving the physical activity level of the
participants. We account for the daily periodic behavior that characterizes
most human activities by modeling the impact of the hour of the day on the
physical activity level as a non-linear function
$f_{H}(\cdot):[7,23)\rightarrow\mathbb{R}$, which is approximated by a linear
combination of $J_{H}$ spline basis functions $\phi_{j}(\cdot)$ with unknown
coefficients $\beta_{H,j}$’s, $\displaystyle
f_{H}(h)\approx\widetilde{f}_{H}(h)=\sum_{j=1}^{J_{H}}\beta_{H,j}B_{H,j}(h)=\boldsymbol{B}_{H}(h)\bbeta_{H}$.
The full process specification yields

$\boldsymbol{Y}(t)=\boldsymbol{X}(t)\bbeta+\boldsymbol{B}_{H}\left(h(t)\right)\bbeta_{H}+\boldsymbol{w}(t)+\boldsymbol{\varepsilon}(t),\qquad
t\in\mathcal{\mathbb{R}^{+}}$ (10)

where $h(\cdot):\mathbb{R}^{+}\rightarrow[7,23)$ links each time point to the
corresponding hour of the day and
$\boldsymbol{B}_{H}\left(\cdot\right):[7,23)\rightarrow\mathbb{R}^{J_{H}}$
links each hour of the day to the values of the splines at that point.

We employ a second order approximation with $4$ internal knots spread
uniformly over the domain. Collecting the basis elements in the design matrix
and stacking the coefficients as for the spatial model in Section 3.4
introduces $6$ additional columns (hence, only $6$ additional slope
parameters) in the design matrix (i.e. the spline basis functions evaluated at
the observed time-points). The large number of observations in each epoch and
the reduced number of knots subdue any concerns surrounding over-fitting and
unrobust inference.

We use a logarthmic transformation, $lMAG_{k}(t)=\log(MAG_{k}(t))$ for
$k=1,2,\ldots,K$ and $t=t_{k1},\ldots,t_{kT_{k}}$. We denote the parameter
associated with variable “varname” as $\beta_{\text{varname}}$ and the levels
of each categorical covariate as $\text{varname}_{(j)}$ for
$j=1,\dots,J_{\text{varname}}$. Hence,

$lMAG_{k}(t)=\beta_{0}+\beta_{\text{BMI}}\cdot\text{BMI}_{k}+\sum_{j=2}^{J_{\text{Sex}}}\beta_{\text{Sex},j}\cdot\mathbb{I}\left(\text{Sex}_{k}=\text{Sex}_{(j)}\right)\\\
+\beta_{\text{BBS}}\cdot\text{BruinBikeShare}(t)+\sum_{j=1}^{J_{H}}\beta_{H,j}B_{H,j}\left(h(t)\right)+w_{k}(t)+\epsilon_{k}(t)\;,$
(11)

where $\mathbb{I}(\cdot)$ denotes the indicator function, $w_{k}(\cdot)$ is
the DAG-based approximation (Section 3.2) for
$\mathcal{GP}(0,c_{\boldsymbol{\theta}}(\cdot,\cdot))$ with
$c_{\boldsymbol{\theta}}(t,t^{\prime})$ the exponential function and
$\epsilon_{k}(t)\overset{\text{iid}}{\sim}\mathcal{N}\left(0,\tau^{2}\right)$.
For the categorical variables “Age” and “Sex”, the levels start from $2$ as
$J=1$ is taken as the baseline, which corresponds to an Asian female with age
between 20 and 25 years. Other socioeconomic factors (e.g. education and
income level) have been excluded from the analysis as they are strongly
associated with ethnicity and age, while Lux (detecting light exposition) was
excluded after some preliminary analysis revealed its low predictive
power333The PASTA-LA study did not contemplate a rigorous protocol for the
light exposition sensor, and hence this variable is likely to not have been
recorded accurately.. We have assigned priors such as
$\bbeta\sim\mathcal{N}_{J}\left(\mathbf{0},10^{6}\cdot\boldsymbol{I}_{J}\right)$,
$\sigma^{2}\sim\mathcal{IG}(2,2)$ and $\tau^{2}\sim\mathcal{IG}(2,2)$ with $J$
being the total number of $\bbeta$ coefficients. The presence of temporal
dependence was investigated through an individual-specific exploratory
analysis on the residuals from an ordinary least squares linear regression.
Subsequently, it was decided that an exponential covariance function for
temporal dependence (corresponding to an Ornstein-Uhlenbeck process) will be a
parsimonious and effective choice to model the behavior of the underlying
residual process.

We implemented Algorithm 1 for (11). Inference was based on $10,000$ posterior
samples retained after convergence was evinced (discarding an initial $5,000$
iterations as burn-in). Point estimates from the standard linear model were
used as starting values for the regression coefficients. The run time of the
collapsed sampler (Section 3.3) on $D_{1}$ was $\approx 15$ hours, achieving a
desirable acceptance rate of $\approx 28\%$ at convergence.

###  4.2 Results from temporal analysis

Table 3 presents the posterior estimates. The regression coefficients were
slightly different from a Bayesian linear regression model, but presented very
similar inference. African Americans, Latinos and Whites revealed higher
values of $lMAG$ than Asian-Americans as did males over females. As expected,
higher age-groups revealed lower $lMAG$. Unsuprisingly, the introduction of
the temporal process effectuates slightly wider credible intervals for the
regression parameters. This results in BMI being marginally less credible from
the temporal process model than from linear regression. What is more suprising
in both models is the slightly negative effect of BBS on physical activity.
This, however, is likely a consequence of the fact that at least $2\cdot
10^{6}$ data points in $D_{1}$ corresponded to individuals outside of Westwood
without access to the program. The average of the $lMAG$ for the reference
individual is represented by the common intercept, which is estimated $\approx
5.514$ by our model. This implies a $MAG$ per minute count of $1,488$, which
corresponds to hard physical activity and an average MET of $\approx 7$
according to the Table 1 and (2). This value, while large, is not surprising
as we are modeling the epochs corresponding to active time.

Parameter | Collapsed Model | Linear regression
---|---|---
Point | Interval | Point | Interval
Intercept | 5.514 | (5.507, 5.520) | 5.872 | (5.854, 5.888)
Eth. Latin-American | 0.166 | (0.149, 0.183) | 0.136 | (0.131, 0.142)
Eth. White | 0.073 | (0.005, 0.095) | 0.081 | (0.076, 0.086)
Eth. Black or other | 0.203 | (0.184, 0.221) | 0.164 | (0.158, 0.170)
Sex Male | 0.017 | (0.003, 0.033) | 0.023 | (0.019, 0.027)
BMI | 0.004 | (-0.002, 0.01) | 0.003 | (0.002, 0.004)
Age [25-35) | -0.106 | (-0.121, -0.091) | -0.124 | (-0.129, -0.119)
Age [35-50) | -0.110 | (-0.131, -0.09) | -0.123 | (-0.129, -0.117)
Age [50-70] | -0.092 | (-0.121, -0.065) | -0.144 | (-0.152, -0.137)
BBS | -0.051 | (-0.066, -0.037) | -0.067 | (-0.071, -0.064)
$\sigma^{2}$ | 1.537 | (1.528, 1.546) |  |
$\phi$ | 0.315 | (0.312, 0.319) |  |
$\tau^{2}$ | 1.138 | (1.135, 1.141) |  |
Metric | Out-of-sample | In-sample | Out-of-sample | In-sample
Coverage | 0.94 | 0.97 | 0.94 | 0.94
RMSPE (r) | 0.60 (1.24) | 0.34 (0.93) | 1 (1.59) | 1 (1.59)
PIW | 4.80 | 4.62 | 6.24 | 6.24
Table 3: Parameter credible intervals, $95\%(2.5\%,97.5\%)$ and predictive
validation for $15\times 10^{3}$ MCMC iterations on $D_{1}$.

The estimate of the temporal decay parameter $\phi$, suggests a fairly sharp
decline in temporal association, which drops below $0.05$ after $\approx 1$
minute (computed, with the exponential covariance function as
$\frac{1}{3\phi}$). While the estimated variance of the temporal process,
$\sigma^{2}$, is slightly larger than $\tau^{2}$, the latter’s estimate
indicates substantial residual variation beyond the temporal
process—motivating our analysis in Section 4.3. The model in (11) delivers
statistical estimates (with uncertainty quantification) of personalized daily
PA profiles for any individual for any day. The splines and the temporal
process combine to capture subject-specific diurnal variation. Figure 6
presents the posterior estimates of daily MAGs (log) of two such individuals
(number $204$ (a) and $188$ (b)) throughout the day to evince the inter-
subject variation. This figure illustrates the need to accommodate variations
among subjects when predicting their daily physical activities.

(a)

(b)

Figure 6: Personalized PA profiles for two individuals estimated with $95\%$
credible intervals (dashed lines) using the spline daily effect and the
temporal process in (11).

An example of the out-of-sample predictions of the temporal process for a set
of $100$ subsequent points from one specific subject (number $77$) is shown in
Figure 7, which demonstrates the proposed model’s ability to interpolate the
lMAG values at unobserved time-points or intervals. The interpolated values
(red dots) provide a slightly over-smoothed but accurate reconstruction of the
held out lMAG (grey dots), which is always included in the corresponding 95%
predictive bounds. This smoother behavior characterizes both in-sample and
out-of-sample predictions when compared to the true values and is not
necessarily a limitation of our model. Indeed, accelerations recorded by
accelerometers are generally noisy, and the predicted values may be
interpreted as a denoised version of the raw signal.

(a)

(b)

Figure 7: Out-of-sample predictions for two random individuals: black dots
(observed values), grey dots (test set), pink dots (oos predictions), dashed
line ($95\%$ confidence intervals)

###  4.3 Spatial-temporal analysis

We consider $D_{2}$ and fit the model in (11), adding a spatial term (Section
3.4) to exploit GPS information. $D_{2}$ is restricted to those observations
recorded in Westwood, Los Angeles. We introduce spatial splines obtained
through the tensor product of two analogous univariate B-spline basis on
longitude and latitude. We choose two bases of degree $2$ with $9$ equally
spaced knots over a square encompassing Westwood. This sums up to
$J_{S}=(7+2)\times(7+2)=81$ terms for our complete spline basis, including the
boundary knots.

In practice, since locations are functions of time through the trajectory
function $\boldsymbol{\gamma}_{k}(\cdot),k=1,\dots,K$ of each individual,
rewrite the time dependent component of the process mean as
$\mu_{k}(t)=\beta_{\text{BBS}}\cdot\text{BruinBikeShare}(t)+\sum_{j=1}^{J_{H}}\beta_{H,j}B_{H,j}\left(h(t)\right)+\sum_{j=1}^{J_{S}}\beta_{S,j}B_{S,j}\left(\gamma_{k}(t)\right)$,
where $B_{S}=B_{X}\otimes B_{Y}$ is the tensor product bivariate spline. Given
the reduced number of knots and the high spatial density of observations in
several areas of the map, over-fitting is not a concern. However, there are
also areas in Westwood that present sparsely observed data-points and the
model can struggle to identify the spline coefficients referred to those
areas, jeopardizing convergence of the MCMC algorithm. Therefore, we consider
the S-Spline (Ridge-like prior) for this application, where the shrinkage
parameter $\lambda$ has been assigned a $\mathcal{G}(1,1)$ prior. Other
parameters have been assigned the same priors of the temporal application in
Section 4.1. Our posterior inference was based on $5,000$ samples retained
after convergence out of $10,000$ MCMC iterations. Fitting the model to
$D_{2}$ required $\approx 30$ hours. The acceptance rate obtained is $\approx
28\%$, supporting the consistency of our adaptive strategy.

###  4.4 Results from spatial-temporal analysis

Table 4 presents parameter estimates and predictive performances of the model
and compares with a standard linear regression model which includes the
spatial spline terms, but neglects the temporal dependence structure.
Conclusions on the regression coefficients are very similar to those from
Section 4.1. However, accounting for the spatial effects allow for easier
interpretation of the age-group regression coefficients: the older the person,
the lower is the expected physical activity level. Also, both models estimate
the effect of BBS as trending slightly negative. This somewhat surprising
finding can be attributed to a few factors. First, the observations after the
BBS launch are mostly from the winter season (February to April, the coldest
months in L.A. together with December), while the others include summer and
autumn (June to November, the warmest months). Given that physical activity
levels tend to be lower in the colder months, there is indication of some
possible confounding between the BBS effect and seasonality. Second, not all
subjects were exposed to the BBS after its launch and, hence, could not take
advantage of it.

In this application, the intercept is estimated to be $\approx 5.31$ by our
model and, hence, a $MAG$ per minute count of $1214$ (slightly lower than in
the temporal application). This would again correspond to vigorous physical
activity and a MET of $\approx 6.5$.

Param. | Linear regression | S-Spline |
---|---|---|---
Point | Interval | Point | Interval |
Intercept | 5.613 | (5.536, 5.689) | 5.31 | (5.29, 5.33) |
Eth. Latin-American | 0.093 | (0.079, 0.108) | 0.114 | (0.069, 0.159) |
Eth. White | 0.053 | (0.040, 0.067) | 0.053 | (0.004, 0.102) |
Eth. Black or other | 0.066 | (0.052, 0.080) | 0.095 | (0.054, 0.135) |
Sex Male | 0.019 | (0.008, 0.029) | 0.021 | (-0.014, 0.055) |
BMI | 0.005 | (0.003, 0.006) | 0.006 | (-0.007, 0.02) |
Age [25-35) | -0.170 | (-0.183, -0.156) | -0.191 | (-0.227, -0.155) |
Age [35-50) | -0.217 | (-0.233, -0.201) | -0.249 | (-0.298, -0.199) |
Age [50-70] | -0.381 | (-0.404, -0.359) | -0.456 | (-0.528, -0.384) |
BBS | -0.008 | (-0.091, -0.071) | -0.107 | (-0.140, -0.073) |
$\sigma^{2}$ |  |  | 1.489 | (1.461, 1.517) |
$\phi$ |  |  | 0.364 | (0.351, 0.376) |
$\tau^{2}$ |  |  | 0.777 | (0.768, 0.786) |
Metric | Out-of-sample | In-sample | Out-of-sample | In-sample |
Coverage | 0.94 | 0.94 | 0.95 | 0.99 |
RMSPE (r) | 0.95 (1.43) | 0.95 (1.43) | 0.53 (1.07) | 0.26 (0.74) |
PIW | 5.62 | 5.62 | 4.78 | 4.74 |
Table 4: Parameter estimates and predictive validation on $D_{2}$.

The estimate of $\phi$ implies that the dependence drops to $0$ in less than a
minute. Unsurprisingly, including the spatial effect and the temporal process
improves predictive performances (RMSPE or PIW in Table 4) over a model
including only spatial effects (using linear regression with splines). The
spatial-temporal model delivers satisfactory coverage and outperforms its
competitor in all of the other indices for the training and testing dataset.

Figure 8(a) shows the estimated spatial surface, while Figure 8(b) presents
the width of the posterior predictive intervals. The map clearly evinces zones
(darker shades of red highlighted with white contours) that tend to depict
high levels of physical activity. For example, the largest dark red blob in
the north center-left almost perfectly tracks the UCLA campus boundary
reflecting a campus environment with active mobility (walking, running,
biking). Other zones of high activity identify with locations where more
participants in the study live, including those residing in student dorms
(northwest corner) and residential areas immediately around and in the
predefined Westwood/UCLA study area (such as the south central zone) or
Century City shopping center (to the east). Lighter shades (orange) correspond
to areas that are less developed (open space), such as the areas in the north
east; or they are areas with a high degree of transportation infrastructure
and traffic (e.g., toward the western boundary). These correspond to highways
(such as the Interstate-405 highway or other vehicular transportation
corridors) that often have lower levels of activity because they inhibit
outdoor physical activities due to noise, pollution, safety, etc. Our analysis
reveals three additional high activity areas that are not gleaned from non-
spatial models: the Los Angeles National Veteran Park; the Century City
shopping center and the Stone Canyon Park. The color gradient closely follows
the spatial characteristics of the Westwood neighborhood and reveal how
spatial patterns can impact physical activity behavior after accounting for
variation attributable to known explanatory variables.

(a)

(b)

Figure 8: (a) Spatially smoothed estimates from a shrinkage spline over
Westwood, Los Angeles; (b) width of $95\%$ posterior predictive intervals for
the shrinkage spline.

Figure 9 shows two examples of observed (left) and reconstructed (right) MAGs
along trajectories carved out by two subjects. We find a good degree of
agreement between the two plots, and the ability of our model to recover the
$lMAG$ in locations where it has not been observed. The reliability of the
predictions can be proved through different metrics and, unsurprisingly,
including the spatial effect and the temporal process improves predictive
performances either in terms of MSPE or PIW. We deliver these personalized
trajectory plots for every subject in the study and also predict personalized
MAGs for each subject along any new trajectory. This enables personalized
recommendations based upon an individual’s health attributes including
suggestions for more effective paths to follow for optimal physical
activities, while also informing community level interventions in the built
environment.

(a)

(b)

Figure 9: Two randomly chosen $lMAG$ trajectories over Westwood from
individual 204 (a) and individual 566 (b). Observed trajectories with gaps are
seen on the left panels; spatially reconstructed (predicted) trajectories are
seen on the right panels.

##  5 Discussion

We have devised a Bayesian modeling framework to conduct fully model-based
inference for high-resolution accelerometer data over trajectories compiled
from the PASTA-LA study. Our key data analytic developments included (i)
modeling dependence over trajectories; (ii) accounting for subject-specific
spatial-temporal variation for daily mobility; and (iii) predicting or
interpolating PA levels across trajectories; and (iv) identify zones of high
physical activity in Westwood, Los Angeles. Our spatiotemporal analysis offers
richer inference and evinces relationships between physical activity levels
and a variety of factors, both at the subject level (e.g., personal
attributes) and as a function of space and time. The temporal process was able
to effectively glean the features of the data at finer resolutions, while the
spatial splines accounted for residual spatial heterogeneity. Accommodating
both temporal dependence and spatial heterogeneity demonstrably improved
predictive ability and enabled us to effectively delineate zones of high
physical activity. Furthermore, the ability of the model to pool information
across individuals at all time points allows us to infer about those who
present sparsely observed space-time points (due to technical issues or
protocol violation). In particular, given our improved predictive power, we
can fill gaps and infer about PA levels with good accuracy and ensure the
desired coverage by our prediction intervals.

Our analysis also resolves practical difficulties in using actigraph data. It
is not cost-effective to deploy research-grade GlobalSat GPS and Actigraph
units as they are very expensive and continued usage requires heavy staff
involvement. Our methods can be applied to analogous, but less complete, data
derived from smart phones and smart watches, then such devices could be
deployed in much larger studies with much larger sample sizes at a fraction of
the cost. Given the spatiotemporal nature of outdoor PA research, our ability
to predict in areas of data missingness drastically improve inference related
to the impacts of the built and natural environments on physical activity and
active mobility.

While our approach offers trajectory-based inference for actigraph data, we
recognize that there are several avenues for further research. Our DAG-based
approach for scalable temporal processes can be further enriched with recent
developments (Katzfuss and Guinness, 2021; Peruzzi et al., 2020), although any
of the methods reviewed and evaluated by Heaton et al. (2019) can be
incorporated into our framework. Finally, there is possible merit in modeling
the activity counts in each axis jointly and relaxing the assumptions of
Gaussianity using recent developments in multivariate spatiotemporal count
models and for non-Gaussian outcomes (see, e.g. Bradley et al., 2018, 2020).

Recent public health reviews call for interdisciplinary technological advances
to more effectively measure spatiotemporal energetics of activity spaces in
obesity and chronic disease research (James et al., 2016; Kestens et al.,
2017; Drewnowski et al., 2020). Individual-level data, at aggregate, can be
used to identify anchor points for physical activity and reveal causal
pathways between built environment exposures and health. Our work is a novel
contribution demonstrating methodologies for how these pressing research
questions may be answered.

## Supplementary Materials

The online supplement (Appendix) provides additional simulation experiments
expanding upon the analysis presented here and further details on MCMC
specifications. Computer programs executing the models and reproducing the
analysis are also provided.

## Funding

The work of the authors have been supported in part by National Science
Foundation (NSF) under grants NSF/DMS 1916349 and NSF/IIS 1562303, and by the
National Institute of Environmental Health Sciences (NIEHS) under grants
R01ES030210 and 5R01ES027027.

## References

  * mig [2017] Accelerometer data collection and processing criteria to assess physical activity and other outcomes: A systematic review and practical considerations.  _Sports medicine_ , 47(9):1821–1845, 2017.
  * Aguilar-Farias et al. [2019] Nicolas Aguilar-Farias, GMEE Peeters, Robert J Brychta, Kong Y Chen, and Wendy J Brown.  Comparing actigraph equations for estimating energy expenditure in older adults.  _Journal of sports sciences_ , 37(2):188–195, 2019.
  * Bradley et al. [2018] Jonathan R. Bradley, Scott H. Holan, and Christopher K. Wikle.  Computationally efficient multivariate spatio-temporal models for high-dimensional count-valued data (with discussion).  _Bayesian Analysis_ , 13(1):253–310, 2018.
  * Bradley et al. [2020] Jonathan R. Bradley, Scott H. Holan, and Christopher K. Wikle.  Bayesian hierarchical models with conjugate full-conditional distributions for dependent data from the natural exponential family.  _Journal of the American Statistical Association_ , 115(532):2037–2052, 2020.  doi:10.1080/01621459.2019.1677471.  URL https://doi.org/10.1080/01621459.2019.1677471.
  * Bull et al. [2020] Fiona C Bull, Salih S Al-Ansari, Stuart Biddle, Katja Borodulin, Matthew P Buman, Greet Cardon, Catherine Carty, Jean-Philippe Chaput, Sebastien Chastin, Roger Chou, et al.  World health organization 2020 guidelines on physical activity and sedentary behaviour.  _British journal of sports medicine_ , 54(24):1451–1462, 2020.
  * Cressie and Wikle [2011] Noel Cressie and Christopher K Wikle.  _Statistics For Spatio-Temporal Data_.  John Wiley & Sons, 2011.
  * Cressie and Wikle [2015] Noel Cressie and Christopher K Wikle.  _Statistics For Spatio-Temporal Data_.  John Wiley & Sons, 2015.
  * Crouter et al. [2006] Scott E Crouter, Kurt G Clowers, and David R Bassett Jr.  A novel method for using accelerometer data to predict energy expenditure.  _Journal of applied physiology_ , 100(4):1324–1331, 2006.
  * Datta et al. [2016a] A. Datta, S. Banerjee, A. O. Finley, N. A. S. Hamm, and M. Schaap.  Non-separable dynamic nearest-neighbor gaussian process models for large spatio-temporal data with an application to particulate matter analysis.  _Annals of Applied Statistics_ , 10:1286–1316, 2016a.  URL http://dx.doi.org/10.1214/16-AOAS931.
  * Datta et al. [2016b] Abhirup Datta, Sudipto Banerjee, Andrew O Finley, and Alan E Gelfand.  Hierarchical nearest-neighbor gaussian process models for large geostatistical datasets.  _Journal of the American Statistical Association_ , 111(514):800–812, 2016b.
  * Drewnowski et al. [2020] Adam Drewnowski, James Buszkiewicz, Anju Aggarwal, Chelsea Rose, Shilpi Gupta, and Annie Bradshaw.  Obesity and the built environment: A reappraisal.  _Obesity_ , 28(1):22–30, 2020.  doi:https://doi.org/10.1002/oby.22672.  URL https://onlinelibrary.wiley.com/doi/abs/10.1002/oby.22672.
  * Eilers and Marx [1996] Paul HC Eilers and Brian D Marx.  Flexible smoothing with b-splines and penalties.  _Statistical science_ , pages 89–102, 1996.
  * Finley et al. [2017] A Finley, A Datta, and S Banerjee.  Spnngp: Spatial regression models for large datasets using nearest neighbor gaussian processes.  _R package version 0.1_ , 1, 2017.
  * Finley et al. [2013] Andrew O Finley, Sudipto Banerjee, and Alan E Gelfand.  spbayes for large univariate and multivariate point-referenced spatio-temporal data models.  _arXiv preprint arXiv:1310.8192_ , 2013.
  * Finley et al. [2019] Andrew O Finley, Abhirup Datta, Bruce D Cook, Douglas C Morton, Hans E Andersen, and Sudipto Banerjee.  Efficient algorithms for bayesian nearest neighbor gaussian processes.  _Journal of Computational and Graphical Statistics_ , pages 1–14, 2019.
  * Freedson et al. [2012] Patty Freedson, Heather R Bowles, Richard Troiano, and William Haskell.  Assessment of physical activity using wearable monitors: Recommendations for monitor calibration and use in the field.  _Medicine and science in sports and exercise_ , 44(1 Suppl 1):S1, 2012.
  * Gelfand et al. [2010] Alan E Gelfand, Peter Diggle, Peter Guttorp, and Montserrat Fuentes.  _Handbook of Spatial Statistics_.  CRC press, 2010.
  * Goodman and Hardin [2006] Tim Goodman and Doug Hardin.  Refinable multivariate spline functions.  In _Studies in Computational Mathematics_ , volume 12, pages 55–83. Elsevier, 2006.
  * Haario et al. [2001] Heikki Haario, Eero Saksman, Johanna Tamminen, et al.  An adaptive metropolis algorithm.  _Bernoulli_ , 7(2):223–242, 2001.
  * Hall et al. [2013] Katherine S Hall, Cheryl A Howe, Sharon R Rana, Clara L Martin, and Miriam C Morey.  Mets and accelerometry of walking in older adults: Standard versus measured energy cost.  _Medicine and science in sports and exercise_ , 45(3):574, 2013.
  * Hastie et al. [2000] Trevor Hastie, Robert Tibshirani, et al.  Bayesian backfitting (with comments and a rejoinder by the authors.  _Statistical Science_ , 15(3):196–223, 2000.
  * Heaton et al. [2019] Matthew J Heaton, Abhirup Datta, Andrew O Finley, Reinhard Furrer, Joseph Guinness, Rajarshi Guhaniyogi, Florian Gerber, Robert B Gramacy, Dorit Hammerling, Matthias Katzfuss, et al.  A case study competition among methods for analyzing large spatial data.  _Journal of Agricultural, Biological and Environmental Statistics_ , 24(3):398–425, 2019.
  * James et al. [2016] Peter James, Marta Jankowska, Christine Marx, Jaime E. Hart, David Berrigan, Jacqueline Kerr, Philip M. Hurvitz, J. Aaron Hipp, and Francine Laden.  “spatial energetics”: Integrating data from gps, accelerometry, and gis to address obesity and inactivity.  _American Journal of Preventive Medicine_ , 51(5):792–800, 2016.  ISSN 0749-3797.  doi:https://doi.org/10.1016/j.amepre.2016.06.006.  URL https://www.sciencedirect.com/science/article/pii/S0749379716302276.
  * Kamada et al. [2016] Masamitsu Kamada, Eric J Shiroma, Tamara B Harris, and I-Min Lee.  Comparison of physical activity assessed using hip-and wrist-worn accelerometers.  _Gait & posture_, 44:23–28, 2016.
  * Katzfuss and Guinness [2021] Matthias Katzfuss and Joseph Guinness.  A general framework for vecchia approximations of gaussian processes.  _Statist. Sci._ , 36(1):124–141, 02 2021.  doi:10.1214/19-STS755.  URL https://doi.org/10.1214/19-STS755.
  * Katzfuss et al. [2020] Matthias Katzfuss, Joseph Guinness, Wenlong Gong, and Daniel Zilber.  Vecchia approximations of gaussian-process predictions.  _Journal of Agricultural, Biological and Environmental Statistics_ , 25:383–414, 2020.
  * Kestens et al. [2017] Yan Kestens, Rania Wasfi, Alexandre Naud, and Basile Chaix.  “contextualizing context”: Reconciling environmental exposures, social networks, and location preferences in health research.  _Current Environmental Health Reports_ , 4:51–60, 2017\.  doi:https://doi.org/10.1007/s40572-017-0121-8.
  * Lang and Brezger [2004] Stefan Lang and Andreas Brezger.  Bayesian p-splines.  _Journal of computational and graphical statistics_ , 13(1):183–212, 2004.
  * Lauritzen [1996] Steffen L Lauritzen.  _Graphical Models_ , volume 17.  Clarendon Press, 1996.
  * Liu et al. [1994] Jun S Liu, Wing Hung Wong, and Augustine Kong.  Covariance structure of the gibbs sampler with applications to the comparisons of estimators and augmentation schemes.  _Biometrika_ , 81(1):27–40, 1994.
  * Lyden et al. [2011] Kate Lyden, Sarah L Kozey, John W Staudenmeyer, and Patty S Freedson.  A comprehensive evaluation of commonly used accelerometer energy expenditure and met prediction equations.  _European journal of applied physiology_ , 111(2):187–201, 2011.
  * Lyden et al. [2014] Kate Lyden, Sarah Kozey Keadle, John Staudenmayer, and Patty S Freedson.  A method to estimate free-living active and sedentary behavior from an accelerometer.  _Medicine and science in sports and exercise_ , 46(2):386, 2014.
  * Montoye et al. [2018] Alexander HK Montoye, Rebecca W Moore, Heather R Bowles, Robert Korycinski, and Karin A Pfeiffer.  Reporting accelerometer methods in physical activity intervention studies: A systematic review and recommendations for authors.  _British journal of sports medicine_ , 52(23):1507–1516, 2018.
  * Murphy [2012] Kevin P Murphy.  _Machine Learning: A Probabilistic Perspective_.  MIT press, 2012.
  * Ott et al. [2000] Alise E Ott, Russell R Pate, Stewart G Trost, Dianne S Ward, and Ruth Saunders.  The use of uniaxial and triaxial accelerometers to measure children’s “free-play” physical activity.  _Pediatric Exercise Science_ , 12(4):360–370, 2000.
  * Peruzzi et al. [2020] Michele Peruzzi, Sudipto Banerjee, and Andrew O. Finley.  Highly scalable bayesian geostatistical modeling via meshed gaussian processes on partitioned domains.  _Journal of the American Statistical Association_ , 0(0):1–14, 2020.  doi:10.1080/01621459.2020.1833889.  URL https://doi.org/10.1080/01621459.2020.1833889.
  * Peterson et al. [2015] Neil E Peterson, John R Sirard, Pamela A Kulbok, Mark D DeBoer, and Jeanne M Erickson.  Inclinometer validation and sedentary threshold evaluation in university students.  _Research in nursing & health_, 38(6):492, 2015\.
  * Piercy et al. [2018] Katrina L Piercy, Richard P Troiano, Rachel M Ballard, Susan A Carlson, Janet E Fulton, Deborah A Galuska, Stephanie M George, and Richard D Olson.  The physical activity guidelines for americans.  _Jama_ , 320(19):2020–2028, 2018.
  * Plasqui and Westerterp [2007] Guy Plasqui and Klaas R Westerterp.  Physical activity assessment with accelerometers: An evaluation against doubly labeled water.  _Obesity_ , 15(10):2371–2379, 2007.
  * Ramsay and Silverman [2007] James O Ramsay and Bernard W Silverman.  _Applied Functional Data Analysis: Methods and Case Studies_.  Springer, 2007.
  * Reiner et al. [2013] Miriam Reiner, Christina Niermann, Darko Jekauc, and Alexander Woll.  Long-term health benefits of physical activity–a systematic review of longitudinal studies.  _BMC public health_ , 13(1):1–9, 2013.
  * Roberts and Rosenthal [2009] Gareth O Roberts and Jeffrey S Rosenthal.  Examples of adaptive mcmc.  _Journal of Computational and Graphical Statistics_ , 18(2):349–367, 2009.
  * Santos-Lozano et al. [2013] Alejandro Santos-Lozano, F Santin-Medeiros, Greet Cardon, Gema Torres-Luque, R Bailon, C Bergmeir, Jonatan R Ruiz, Alejandro Lucía Mulas, Nuria Garatachea, et al.  Actigraph gt3x: validation and determination of physical activity intensity cut points.  2013\.
  * Sasaki et al. [2011] Jeffer E Sasaki, Dinesh John, and Patty S Freedson.  Validation and comparison of actigraph activity monitors.  _Journal of science and medicine in sport_ , 14(5):411–416, 2011.
  * Sikka et al. [2019] Robby S Sikka, Michael Baer, Avais Raja, Michael Stuart, and Marc Tompkins.  Analytics in sports medicine: Implications and responsibilities that accompany the era of big data.  _JBJS_ , 101(3):276–283, 2019.
  * Stein et al. [2004] Michael L Stein, Zhiyi Chi, and Leah J Welty.  Approximating likelihoods for large spatial data sets.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 66(2):275–296, 2004.
  * Taraldsen et al. [2012] Kristin Taraldsen, Sebastien FM Chastin, Ingrid I Riphagen, Beatrix Vereijken, and Jorunn L Helbostad.  Physical activity monitoring by use of accelerometer-based body-worn sensors in older adults: A systematic literature review of current knowledge and applications.  _Maturitas_ , 71(1):13–19, 2012.
  * Vecchia [1988] Aldo V Vecchia.  Estimation and model identification for continuous spatial processes.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 50(2):297–312, 1988.

## Appendix

##  6 Additional simulation experiments

We carried out two additional experiments to test the reliability of our
algorithm and verify comparative performances with the Sequential NNGP as it
is implemented in the spNNGP package [Finley et al., 2017]. We did not
consider the Response NNGP because it does not recover the latent component.
The first one is described in Section 6.1 and includes simulated observations
for one single individual; the second one includes simulated observations for
multiple individuals and is described in Section 6.2.

###  6.1 Experiment 1

We generated observations $\left\\{y(t_{j})\right\\}_{j=1}^{\top}$ for $K=1$
individual, using $T=10^{5}$ time-points, where each
$t_{i}=\sum_{h=1}^{i-1}\delta_{h}$, and $\delta_{h}\sim Exp(5),\,\forall\,h$.
The model included an intercept $\beta_{0}$ and 3 covariates, $x_{1}$, $x_{2}$
and $x_{3}$ all drawn from $\mathcal{N}(0,1)$, with associated slopes
$\beta_{1}$, $\beta_{2}$ and $\beta_{3}$. We modeled the covariance structure
between any two simulations at time-points $t$ and $t^{\prime}$ using the
exponential covariance function:

$\mathbf{Cov}_{\boldsymbol{\theta}}\left[Y(t),Y(t^{\prime})\right]=c_{\boldsymbol{\theta}}(t,t^{\prime})=\sigma^{2}e^{(\phi|t-t^{\prime}|)},\quad\sigma^{2},\phi\in\mathbb{R}^{+},$
(12)

where $\sigma^{2}$ represents the variance of the process (sill), $\phi$ is
the decay in temporal correlation (range) and $\tau^{2}$ the residual variance
(nugget). In this data generation step the parameters have been set to the
following values: $\beta_{0}=-1.878$, $\beta_{1}=0.326$, $\beta_{2}=-0.302$,
$\beta_{3}=1.182$, $\sigma^{2}=\phi=\tau^{2}=1$. A chunk of the simulated
trajectory and its density can be observed as an example in Figures 10(a) and
10(b), respectively.

(a)

(b)

Figure 10: Simulated uni-dimensional Gaussian process (a) and its density (b).

We fitted the model on the simulated data using our Collapsed NNGP
implementation, specifically optimized for the temporal setting, while fitting
the Sequential NNGP using the spNNGP package. The latter, while generally used
for fitting spatial (i.e. two-dimensionals) models, can be adapted to the
temporal (uni-dimensional) case by providing a set of locations where $t$ is
one of the coordinates and the other is fixed to a constant value (e.g.
$\left\\{\widetilde{\boldsymbol{s}}_{j}\right\\}_{j=1}^{\top}=\left\\{(t_{j},0)\right\\}_{j=1}^{\top}$).
In our implementation, the intercept and slope regression parameters were
given a vague normal prior distribution $\mathcal{N}(0,10^{6})$. The variance
components, $\sigma^{2}$ and $\tau^{2}$, were both assigned an inverse Gamma
prior $\mathcal{IG}(2,2)$, and the decay parameter $\phi$ was ascribed a
$\mathcal{G}(1,1)$. On the other hand, the spNNGP assumes a flat prior on the
intercept and slope coefficients and a uniforma $\mathcal{U}(a,b)$ prior on
the decay parameter $\phi$. In this experiment we fixed $a=0.5$ and $b=30$.
All the models were trained on the same random sample composed of the $70\%$
of the total observations, while the the remaining $30\%$ have been excluded
to assess the out-of-sample predictive performances in terms Relative Mean
Squared Prediction Error (RMSPE), Root Mean Squared Prediction Error (rMSPE),
Coverage, Predictive Interval Width (PIW).

Param. (True) | Collapsed NNGP | Sequential NNGP
---|---|---
Point | Interval | ESS | Point | Interval | ESS
$\beta_{0}\;(-1.88)$ | -1.87 | (-1.89, -1.85) | 4999 | -1.87 | (-1.89, -1.85) | 57
$\beta_{1}\;(0.33)$ | 0.33 | (0.32, 0.34) | 4999 | 0.33 | (0.32, 0.34) | 1285
$\beta_{2}\;(-0.30)$ | -0.30 | (-0.31, -0.29) | 4999 | -0.30 | (-0.31, -0.3) | 1365
$\beta_{3}\;(1.18)$ | 1.18 | (1.17, 1.19) | 4999 | 1.18 | (1.17, 1.19) | 1342
$\sigma^{2}\;(1)$ | 1.00 | (0.97, 1.03) | 472 | 1.00 | (0.97, 1.03) | 294
$\phi\;(1)$ | 0.99 | (0.95, 1.04) | 496 | 0.99 | (0.95, 1.04) | 65
$\tau^{2}\;(1)$ | 1.01 | (0.99, 1.03) | 457 | 1.01 | (0.99, 1.03) | 165
Metric | Out-of-sample | In-sample |  | Out-of-sample | In-sample |
Coverage | 0.95 | 0.99 |  | 0.96 | 0.99 |
RMSPE (r) | 0.39 (1.19) | 0.20 (0.85) |  | 0.39 (1.19) | 0.20 (0.85) |
PIW | 4.68 | 4.46 |  | 4.78 | 4.47 |
Run time (h) | 1.77 | 1.86
Table 5: Parameter estimates, predictive validation and fitting times (hours)
on the simulated dataset for all the considered models.

We ran the $10000$ MCMC iterations, fixing the number of neighbours $m=10$.
The first $5000$ simulations have been dropped as burn-in, while the last
$5000$ have been retained for estimation and prediction purposes. No thinning
has been considered. Results are summarized in Table 5. The two approaches
provide identical outputs, both in terms of estimation and prediction.
However, our implementation is faster than its competitor (at least in the
context of the temporal setting) and provides way better performances in terms
of Effective Sample Size (ESS).

####  6.1.1 Computation time evaluation

In order to delve more into the computational aspect, we quantified the
linearity of all the algorithms: by construction, the fitting time should
increase linearly with the sample size. We split observations in $l=1,\dots,5$
different fitting windows $\left\\{t_{1},\dots,t_{T_{l}}\right\\}$ with
increasing sizes $T_{l}=\\{\\{2^{l}\\}_{l=0}^{6}\cup\\{100\\}\\}\times 10^{3}$
and the computation time of one sampler iteration, fixing $m=30$ for all the
considered algorithms has been recorded for $\widetilde{M}=100$ times. Figure
11 shows that all algorithms scale linearly with the sample size. However, our
implementation of the collapsed NNGP, whilst pointed out as generally less
efficient than its competitors in Finley et al. [2019], scales with a rate of
$\approx 0.376\cdot 10^{-4}$ per data-point, while the Sequential NNGP scale
with a rate equal to $\approx 4.5736\cdot 10^{-4}$, which is sensibly higher.
For an exact numerical analysis, results are reported in Table 6.

Figure 11: Time elapsed (in seconds) for 1 MCMC iteration for the two considered algorithms with increasing sample size T and fixed $m=30$. $T\times 10^{3}$ | Algorithm | Min | $q_{025}$ | Median | Mean | $q_{975}$ | Max
---|---|---|---|---|---|---|---
$1$ | Collapsed | 0.01 | 0.01 | 0.02 | 0.02 | 0.03 | 0.03
Sequential | 0.12 | 0.12 | 0.13 | 0.14 | 0.18 | 0.18
$2$ | Collapsed | 0.03 | 0.03 | 0.03 | 0.04 | 0.06 | 0.07
Sequential | 0.25 | 0.25 | 0.26 | 0.27 | 0.31 | 0.34
$4$ | Collapsed | 0.06 | 0.06 | 0.07 | 0.09 | 0.16 | 0.16
Sequential | 0.50 | 0.50 | 0.52 | 0.64 | 1.21 | 1.21
$8$ | Collapsed | 0.13 | 0.14 | 0.30 | 0.26 | 0.32 | 0.41
Sequential | 1.01 | 1.01 | 2.34 | 1.99 | 2.41 | 2.56
$16$ | Collapsed | 0.27 | 0.28 | 0.60 | 0.46 | 0.63 | 0.64
Sequential | 2.02 | 2.02 | 4.65 | 3.46 | 4.75 | 4.77
$32$ | Collapsed | 0.55 | 0.56 | 1.23 | 1.17 | 1.28 | 1.37
Sequential | 4.08 | 4.09 | 9.40 | 8.87 | 9.59 | 10.16
$64$ | Collapsed | 1.01 | 1.03 | 2.46 | 1.90 | 2.79 | 2.85
Sequential | 2.51 | 7.49 | 18.71 | 14.43 | 20.13 | 20.69
$100$ | Collapsed | 1.60 | 1.61 | 1.67 | 1.68 | 1.87 | 1.99
Sequential | 11.68 | 11.74 | 11.93 | 12.01 | 12.80 | 13.87
Table 6: Time (in seconds) of one MCMC iteration for the two considered
algorithms with increasing sample size (T) and fixed $m=30$. Figure 12: Time
elapsed (in seconds) for 1 MCMC iteration for the Collapsed NNGP and the Full
GP with increasing sample size T.

Additionally, we wanted to quantify the computational advantage of the
proposed NNGP-based collapsed algorithm over the standard MCMC update of the
full GP model [Cressie and Wikle, 2015], where the latter is already
implemented in R through spLM() function of the spBayes package [Finley et
al., 2013]. Since Datta et al. [2016b] proved that the NNGP approximation with
$30$ neighbours provides almost exactly the same inference of the full GP, we
fixed $m=30$ and built again different sets of data with increasing number of
data-points, but this times with $T_{l}=100,1000,5000,10000$ (sizes have been
reduced to comply with the slow update of the Full GP). Results, which are
summarized in Figure 12, showed that for $T_{l}=100$ the computational time
difference between the full GP and the collapsed NNGP is negligible. However,
as the size increases, the saving of time increases exponentially: 15 seconds
(per iteration) when $n=5000$, $122$ seconds per iteration when $n=10000$. For
the last scenario, which is the more realistic in a MCMC inference context,
this means that the collapsed NNGP will provide us with the same results 14
days in advance with respect to the Full GP model.

###  6.2 Experiment 2

The aim of this experiment is to verify the ability of our algorithm in
recovering the true parameters and to determine if pooling information from
multiple individuals can help in improving the accuracy of the estimates.
Comparison with the Sequential NNGP is not feasible, since it does not allow
the contemporary fitting of multiple Gaussian processes with common
parameters. Thus, we compare performances of the Pooled NNGP (that’s how we
will refer to the proposed collapsed in what follows) with the single models
estimated separately for each individual.

|  |  | Individuals
---|---|---|---
Param. | True | Pooled NNGP | 1 | 2 | 3 | 4 | 5
$\beta_{01}$ | -9.39 | -9.41 (-9.46, -9.37) | -9.41 (-9.46, -9.37) |  |  |  |
$\beta_{02}$ | 1.63 | 1.59 (1.54, 1.64) |  | 1.59 (1.54, 1.64) |  |  |
$\beta_{03}$ | -1.51 | -1.53 (-1.57, -1.48) |  |  | -1.54 (-1.58, -1.49) |  |
$\beta_{04}$ | 5.91 | 5.91 (5.86, 5.96) |  |  |  | 5.91 (5.86, 5.96) |
$\beta_{05}$ | -0.82 | -0.80 (-0.85, -0.76) |  |  |  |  | -0.81 (-0.86, -0.76)
$\beta_{1}$ | 6.48 | 6.48 (6.47, 6.49) | 6.48(6.46,6.5) | 6.48(6.46,6.50) | 6.48(6.46,6.50) | 6.46(6.44,6.48) | 6.49(6.47,6.51)
$\beta_{2}$ | 6.76 | 6.75 (6.74, 6.76) | 6.74(6.72,6.76) | 6.75(6.74,6.77) | 6.75(6.73,6.77) | 6.75(6.73,6.77) | 6.76(6.74,6.78)
$\beta_{3}$ | -1.46 | -1.46 (-1.47, -1.45) | -1.47(-1.49,-1.45) | -1.48(-1.50,-1.46) | -1.47(-1.48,-1.45) | -1.45(-1.47,-1.43) | -1.47(-1.49,-1.45)
$\sigma^{2}$ | 1 | 0.98 (0.96, 1.01) | 0.94(0.88,1.00) | 1.02(0.953,1.085) | 0.93(0.88,0.99) | 1.03(0.97,1.1) | 1.01(0.94,1.08)
$\phi$ | 1 | 1.01 (0.97, 1.06) | 1.00(0.91,1.11) | 1.03(0.93,1.14) | 1.06(0.95,1.17) | 0.94(0.85,1.04) | 0.99(0.9,1.11)
$\tau^{2}$ | 1 | 1 (0.99, 1.02) | 1.04(1.00,1.08) | 0.98(0.94,1.02) | 1.00(0.967,1.04) | 1.00(0.96,1.04) | 0.99(0.95,1.03)
Coverage |  | 0.95 (0.99) | 0.95 (0.99) | 0.95 (0.99) | 0.95 (0.99) | 0.96 (0.99) | 0.95 (0.99)
RMSPE |  | 0.012 (0.006) | 0.02 (0.01) | 0.02 (0.01) | 0.02 (0.01) | 0.02 (0.01) | 0.02 (0.01)
rMSPE |  | 1.22 (0.84) | 1.24 (0.83) | 1.24 (0.83) | 1.22 (0.85) | 1.22 (0.85) | 1.23 (0.84)
PIW |  | 4.67 (4.44) | 4.99 (4.44) | 4.94 (4.43) | 4.93 (4.43) | 4.95 (4.44) | 4.94 (4.44)
Fitting time |  | $0.59$ | $0.33$ | $0.30$ | $0.35$ | $0.34$ | $0.37$
Table 7: Parameter estimates ($95\%$ credible intervals), out-of-sample
prediction error and fitting times (hours) on the simulated dataset for the
pooled and the single models.

We generated $2\cdot 10^{4}$ observations for $K=5$ individuals, using the
same scheme of Experiment 1 (total of $10^{5}$ data-points). Results are
presented in Table 7. The model also included 3 covariates and an intercept
for each individual drawn from independent $N(0,1)$. Observations were then
generated as described in Section 6.1. The simulated data was split into two
sets: $70\%$ composed the train set for estimation purposes, while the
remaining $30\%$ was used to asses model predictive performances. RMSPE,
coverage of the predictive $95\%$ credible intervals and their mean width were
used as measures of the goodness of fit. For all the models, the intercept and
slope regression parameters were given a flat normal prior distribution
$\mathcal{N}(0,10^{6})$. The variance components, $\sigma^{2}$ and $\tau^{2}$,
were both assigned an inverse Gamma $\mathcal{IG}(2,2)$ priors, and the decay
parameter $\phi$ received a Gamma prior $\mathcal{G}(1,1)$. The advantage of
pooling information from multiple individuals for the estimation of common
parameters, while the independence assumption among them still holds, is
evident according to all criteria. First of all, there is a sensible gain in
the estimation accuracy of the common parameters. Indeed, while the true value
of the parameters are included in the intervals also considering one single
individual at a time, the widths of $95\%$ credible intervals are sensibly
smaller when we pool information together. Furthermore, some slight advantage
is also visible for prediction purposes, where the Pooled NNGP provides larger
coverage and smaller RMSPE. Additionally, thanks to parallelization of the
code, there is almost no loss in terms of the computational time required for
the fitting: $\approx 40$ minutes to fit one individual VS $\approx 55$
minutes to fit the pooled model.