# A fully Bayesian semi-parametric scalar-on-function regression (SoFR) with
measurement error using instrumental variables.

Roger S. Zoh, Yuanyuan Luan, Carmen Tekwe

###### Abstract

Wearable devices such as the ActiGraph are now commonly used in health studies
to monitor or track physical activity. This trend aligns well with the growing
need to accurately assess the effects of physical activity on health outcomes
such as obesity. When accessing the association between these device-based
physical activity measures with health outcomes such as body mass index, the
device-based data is considered functions, while the outcome is a scalar-
valued. The regression model applied in these settings is the scalar-on-
function regression (SoFR). Most estimation approaches in SoFR assume that the
functional covariates are precisely observed, or the measurement errors are
considered random errors. Violation of this assumption can lead to both under-
estimation of the model parameters and sub-optimal analysis. The literature on
a measurement corrected approach in SoFR is sparse in the non-Bayesian
literature and virtually non-existent in the Bayesian literature. This paper
considers a fully nonparametric Bayesian measurement error corrected SoFR
model that relaxes all the constraining assumptions often made in these
models. Our estimation relies on an instrumental variable (IV) to identify the
measurement error model. Finally, we introduce an IV quality scalar parameter
that is jointly estimated along with all model parameters. Our method is easy
to implement, and we demonstrate its finite sample properties through an
extensive simulation. Finally, the developed methods are applied to the
National Health and Examination Survey to assess the relationship between
wearable-device-based measures of physical activity and body mass index among
adults living in the United States.

##  1 Introduction

Regression settings where the number of covariates $p$ exceeds the sample
size, $n$, coined as “large-p–small-n” scenarios, have gained considerable
attention recently. Two main estimation approaches are employed in such
settings. The first approach assumes that no a priori order exists for the
covariates and attempts to estimate the regression parameters for each
important covariate while driving regression parameters for unimportant
covariates to zero, effectively performing variable selection. See for example
Hastie et al. (2015) for a frequentist overview, and Van Erp et al. (2019) for
the Bayesian counterpart. The second approach to large-p–small-n regression
analysis falls under the functional data analysis functional data analysis
umbrella. Functional data analysis fundamentally assumes that there exists a
functional regression parameter, often with a smooth functional form, that
describes the effects of potentially large functional covariates on the
response (Ramsay, 2006). Scalar-on-function (SoFR), an example of functional
data analysis, is increasingly employed in biomedical studies to evaluate the
relationship between function-valued covariates and scalar-valued outcomes. In
these settings, the functional covariates are observed on a dense or sparse
grid. There exists an extensive literature devoted to SoFR, for example,
Ramsay (2006); Wang et al. (2016), provide a general introduction to
functional data analysis data analysis and SoFR. Additionally, Morris (2015)
provides a great review on functional regression. Most of these prior works
focus on cases where non-functional or scalar-valued covariates are ignored.
Here we assume a broader definition of SoFR to allow for both function and
scalar-valued covariates. Namely, let $\\{X(t):t=1,\cdots,T<\infty\\}$, $Z$,
and $Y$ denote the functional covariate, scalar-valued covariate, the scalar-
valued response, respectively. Then the SoFR model is represented as:

$\displaystyle Y_{i}$ $\displaystyle=$
$\displaystyle\beta_{0}+\mbox{\boldmath$\beta$}_{z}^{\rm
T}{\mathbf{Z}}_{i}+\int_{S}\beta(s)X_{i}(s)\,dx+\epsilon_{i}\quad\forall
i=1,\cdots,n.$ (1)

In (1), we assume ${\mathbf{Z}}\in\mathbb{R}^{n\times p}$, where $p<<n$.
Estimation of model (1) has extensively been discussed in the non-Bayesian
framework using techniques developed in the FDA literature. Under these
approaches, estimation is mostly done using a pre-selected set of basis
functions such as B-splines, thin plate, polynomial splines, and wavelet or
based on data driven set of basis functions obtained by functional principal
components combined with some penalty term to impose various desired levels of
smoothness on the functional solution. See for example Reiss et al. (2017) for
a nice review of more recent advances in SoFR approaches.

This manuscript is motivated by the need to accurately evaluate the effects of
physical activity patterns on the body mass index (BMI) for a sample of adults
living in the United States. In this motivating example, Actigraph, a wearable
device used to monitor physical activity is used to approximate the individual
patterns of physical activity of the study participants, while BMI is used as
an indicator for adiposity or obesity. Devices such as Actigraph collect data
frequently such as at $60$ second epochs while the device is worn, resulting
in high dimensional longitudinal data or function data with complex
heterogeneous covariance structures. With the rise in the use of wearable
devices to track physical activity, several authors have recently questioned
their accuracy in monitoring physical activity (Feito et al., 2012; O’Driscoll
et al., 2020; An et al., 2017). Functional data analysis, specifically, SoFR
can be applied to assess how the function-valued data obtained from wearable
devices influence scalar-valued outcomes such as BMI. However, current SoFR
approaches are limited in their ability to correct for potential biases due to
measurement errors when estimating the association between imprecisely
observed functional covariates prone to complex error structures and scalar-
valued outcomes.

While, there is a growing interest in addressing measurement error in SoFR Yao
et al. (2005); Cardot et al. (2007); Goldsmith et al. (2013); Morris (2015);
Reiss et al. (2017), most of these approaches have assumed that the function-
valued covariates are observed at discrete time points with independence error
structures. In contrast, Chakraborty and Panaretos (2017) used a multivariate
regression calibration approach while allowing for a correlated measurement
error. More recently, Tekwe et al. (2019) proposed a measurement error-
adjusted SoFR model that relies on an instrumental variable to identify the
model under the assumption of heteroscedastic error structures for the
measurement error. The model parameters were estimated using generalized
method of moments (GMM), and bootstrap confidence intervals were obtained for
the estimated functional parameter. However, some limitations of the GMM
approach include the requirement of two-stage estimation and their proposed
GMM method does not allow the inclusion of error free covariates. Finally,
another limitation in current SoFR models is that the error term is often
assumed to have a normal or symmetric distributional form, because interest is
primarily on the mean portion of the model.

Here, our contribution differs from previously proposed methods for addressing
measurement error in SoFR in several ways. First, we adjust for measurement
error using a fully Bayesian semi-parametric method, which allows for
arbitrary measurement error covariance structures. Our Bayesian approach
eliminates the need for two-stage estimation in the presence of error-free
covariates, while allowing for an accurate uncertainty quantification in the
model parameters estimates. Because Bayesian inference relies on the posterior
distribution, quantification of the uncertainty is automatic. Second, we
propose the use of a function-valued instrumental variable to obtain a
reasonable estimate of the covariance matrix for the measurement error. The
chosen instrumental variable can be potentially biased for the true unknown
functional in a functional fashion. We provide a way to estimate this
additional scaling functional parameter. introduce an additional parameter in
the model to quantify the fit of the chosen instrumental variable as an
instrument for the true or latent covariate. We term this parameter a “global
instrumental variable quality parameter”.

The remainder of the paper is structured as follows. In section 2, we specify
the model along with all the priors, model assumptions, and the posterior
distributions. We evaluate the performance of our approach in section 3
through an extensive simulation study. Our proposed method is applied to the
National Health and Examination Survey (NHANES) data set to unbiasedly assess
the association between wearable-device-based measures of physical activity
intensity and body mass index (BMI) among adults living in the U.S. in section
4. We conclude with some concluding remarks and a discussion of some potential
extensions in section 5.

##  2 Model

###  2.1 Model Specification

Suppose $\left\\{\mbox{\boldmath$W$}_{i},Y_{i}\right\\}$ is the data pair
collected for each individual $i=1,\cdots,n$, where
$\mbox{\boldmath$W$}_{i}=\\{W_{i}(s),s\in\mbox{\boldmath$S$}\\}$ and $S$ is
the collection of indexes. The SoFR model with functional covariate subject to
measurement error is formulated as follows

$\displaystyle Y_{i}$ $\displaystyle=$
$\displaystyle\alpha_{0}+{\mbox{\boldmath$\beta$}}_{z}^{\rm
T}\mbox{\boldmath$Z$}_{i}+\int_{\mbox{\boldmath$S$}}\mbox{\boldmath$\beta$}(s)X_{i}(s)\,ds+\epsilon_{i}$
(2) $\displaystyle W_{i}(s)$ $\displaystyle=$ $\displaystyle
X_{i}(s)+U_{i}(s),$ (3)

similar to the model considered by (Tekwe et al., 2019). In model 2,
$\mbox{\boldmath$X$}=\\{X(s),s\in\mbox{\boldmath$S$}\\}$ is the true but
unobserved (latent) functional covariate for which we would to objective
estimate its effect $\mbox{\boldmath$\beta$}(s)$ on the response $Y$. However,
we only observe $W$, a contaminated (error prone) but unbiased proxy of $X$.
$\epsilon_{i}$ is the error for the response model and $U_{i}(s)$ the error
process for the measurement error part of the model. Without further
assumptions on the error structure $U$ (see Carroll et al., 2006; Yao et al.,
2005; Cardot et al., 2007; Goldsmith et al., 2013; Morris, 2015; Reiss et al.,
2017)) or the existence of replicate proxies of the functional
$\mbox{\boldmath$X$}(s)$, estimation of model 2 parameters is virtually
impossible. Recently, (Tekwe et al., 2019) provided an approach to estimating
model 2 assuming the existence of an additional replicate of $W$ in the form
of an instrumental variable (IV). However, this assumption may reveal too
constraining since it imposes that the IV be an unbiased proxy of $X$ in
additions to the usual conditions required of an IV. Here we propose the
following model

$\displaystyle Y_{i}$ $\displaystyle=$
$\displaystyle\alpha_{0}+\mbox{\boldmath$\beta$}_{z}^{\rm
T}\mbox{\boldmath$Z$}_{i}+\int_{\mbox{\boldmath$S$}}\beta(s)X_{i}(s)ds+\epsilon_{i}$
(4) $\displaystyle W_{i}(s)$ $\displaystyle=$ $\displaystyle
X_{i}(s)+U_{i}(s),$ (5) $\displaystyle M_{i}(s)$ $\displaystyle=$
$\displaystyle\mbox{\boldmath$\delta$}(s)X_{i}(s)+\omega_{i}(s);$ (6)

where we have an additional model for
$\mbox{\boldmath$M$}_{i}=\\{M_{i}(s),s\in S\\}$ of the observed value for the
IV. We assume $\epsilon_{i}\stackrel{{\scriptstyle
iid}}{{\sim}}\mbox{\boldmath$G$}$, for some distribution $G$ with mean zero
and finite second moment; $E\left\\{U_{i}(s)\right\\}=0$,
$E\left\\{U_{i}(s)U_{i}(t)\right\\}=\sigma_{u,s,t}<\infty$;
$E\\{\omega_{i}(s)\\}=0$, and
$E\\{\omega_{i}(t)\omega_{i}(s)\\}=\sigma_{\omega,t,s}$. We assume all the
model error terms are independent. Estimating $\mbox{\boldmath$\delta$}(s)$
reveals challenging during our Bayesian estimation. In what follows, we assume
$\mbox{\boldmath$\delta$}(s)\neq 0$ is known for $s\in\mbox{\boldmath$S$}$ and
provide an approach to estimating $\mbox{\boldmath$\delta$}(s)$ later . Given
the value of $\mbox{\boldmath$\delta$}(t)\neq 0$, estimation is based on the
scaled version of $M^{*}_{i}(s)=M_{i}(s)/\mbox{\boldmath$\delta$}(s)$ and we
have the following model

$\displaystyle M^{*}_{i}(s)$ $\displaystyle=$ $\displaystyle
X_{i}(s)+\omega^{*}_{i}(s);$ (7)

The number of parameters in models 5 and 6 can be extremely large making model
estimation nearly impossible. We proceed to perform dimension reduction using
a basis expansion. Let $K$ be a given set of orthogonal basis functions
stacked as $\mbox{\boldmath$B$}(s)=\left(b_{s,1},\cdots,b_{s,n}\right)^{\rm
T}\in\mathbb{R}^{n\times K}$. Applying the basis expansion to the functional
covariates, the models in 4, 5, and 6 become

$\displaystyle Y_{i}$ $\displaystyle=$
$\displaystyle\alpha_{0}+\mbox{\boldmath$\beta$}_{z}\mbox{\boldmath$Z$}_{i}+\sum^{K}_{k=1}\gamma_{k}\widetilde{X}_{i,k}+\widetilde{\epsilon}_{i}$
(8) $\displaystyle\widetilde{W}_{i,k}$ $\displaystyle=$
$\displaystyle\widetilde{X}_{i,k}+\widetilde{U}_{i,k},\mbox{where}\;k=1,\cdots,K$
(9) $\displaystyle\widetilde{M}_{i,k}$ $\displaystyle=$
$\displaystyle\widetilde{X}_{i,k}+\widetilde{\omega}^{*}_{i,k},\mbox{where}\;k=1,\cdots,K$
(10)

where $\widetilde{X}_{i,k}=\int_{S}\mbox{\boldmath$B$}_{.k}(s)X_{i}(s)ds$,
with $\mbox{\boldmath$B$}_{.k}$ denoting the $k^{th}$ (column $k$) basis
function; and $\gamma_{k}=\int_{S}\mbox{\boldmath$B$}_{.k}(s)\beta(s)ds$;
$\widetilde{W}_{i,k}=\int_{S}\mbox{\boldmath$B$}_{.k}(s)W_{i}(s)ds$;
$\widetilde{M}_{i,k}=\int_{S}\mbox{\boldmath$B$}_{.k}(s)M^{*}_{i}(s)ds$;
$\widetilde{U}_{i,k}=\int_{S}\mbox{\boldmath$B$}_{.k}(s)U_{i}(s)ds$;
$\widetilde{\omega}^{*}_{i,k}=\int_{S}\mbox{\boldmath$B$}_{.k}(s)\omega^{*}_{i}(s)ds$.

###  2.2 Model Assumptions and Priors

To estimate the proposed model parameters, we adopt a fully Bayesian approach
which conveys many advantages over non-Bayesian approaches, including the
estimation for the model parameters and their related uncertainties completely
based on the posterior distribution.

We now discuss the choices of the model assumptions and prior distributions.
Let’s first provide a stack version of the data and model parameter. We define
for each $i=1,\cdots,n$,
$\widetilde{\mbox{\boldmath$W$}}_{i}=(\widetilde{W}_{i,1},\cdots,\widetilde{W}_{i,K})^{\rm
T}$;
$\widetilde{\mbox{\boldmath$X$}}_{i}=(\widetilde{X}_{i,1},\cdots,\widetilde{X}_{i,K})^{\rm
T}$;
$\widetilde{\mbox{\boldmath$\omega$}}^{*}_{i}=(\widetilde{\omega}^{*}_{i,1},\cdots,\widetilde{\omega}^{*}_{i,K})^{\rm
T}$ ;
$\widetilde{\mbox{\boldmath$U$}}_{i}=(\widetilde{U}_{i,1},\cdots,\widetilde{U}_{i,K})^{\rm
T}$. To allow for flexibility in the error terms of the model, we will assume
a truncated Dirichlet Process mixture (tDPMM) prior similar to the one
considered in Sarkar et al. (2018). This provides considerable flexibility and
help to capture departure from normal (symmetric) distributions of the error
terms. The parameters will be estimated conditional on $K$, the number of
basis functions. We also consider the following conjugate priors which permit
straightforward implementation of the Gibbs steps for obtaining samples from
our complex joint posterior distributions.
$\epsilon_{i}\stackrel{{\scriptstyle
iid}}{{\sim}}\sum^{K_{\epsilon}}_{k=1}\pi_{k}\hbox{Normal}(.|\mu_{k},\sigma^{2}_{k})$,
where $\sum^{K_{\epsilon}}_{k=1}\pi_{k}\mu_{k}=0$, and
$(\mu_{k},\sigma^{2}_{k})\sim\mbox{NIG}(\mu_{0},\nu_{0})$, for
$k=1,\cdots,K_{\epsilon}$;
$\mbox{\boldmath$\pi$}_{\epsilon}=(\pi_{1},\cdots,\pi_{K_{\epsilon}})^{\rm
T}\sim\hbox{Dirch}(\alpha_{\epsilon}/K_{\epsilon},\cdots,\alpha_{\epsilon}/K_{\epsilon})$,
where $K_{\epsilon}$ denotes the class membership probabilities for the error
$\epsilon_{i}$.

For the measurement error vector,
$\widetilde{\mbox{\boldmath$U$}}_{i}|\widetilde{\mbox{\boldmath$X$}}_{i}\stackrel{{\scriptstyle
iid}}{{\sim}}\sum^{K_{u}}_{k=1}\pi_{u,k}\hbox{MVN}_{K}(\mbox{\boldmath$\mu$}_{u,k},\Sigma_{k,u}),$
with
$\mbox{\boldmath$\mu$}_{u,k}\sim\hbox{MVN}_{K}(\mbox{\boldmath$\mu$}_{0},\Sigma_{0})$,
$\Sigma_{k,u}\sim\mbox{Inv-W}(\nu_{0},\Psi_{0,u})$ constraint to
$\sum^{K_{\epsilon}}_{k=1}\pi_{u,k}\mbox{\boldmath$\mu$}_{u,k}={\mathbf{0}}$;
with $\mbox{\boldmath$\pi$}_{u}=(\pi_{1},\cdots,\pi_{K_{u}})^{\rm
T}\sim\hbox{Dirch}(\alpha_{u}/K_{u},\cdots,\alpha_{u}/K_{u})$. For the error
model in the instrumental variable (IV) equation,
$\widetilde{\mbox{\boldmath$\omega$}}^{*}_{i}|\mbox{\boldmath$X$}_{i}\stackrel{{\scriptstyle
iid}}{{\sim}}\sum^{K_{\omega}}_{k=1}\pi_{\omega,k}\hbox{Normal}_{k_{\omega}}(\mbox{\boldmath$\mu$}_{\omega,k},\Sigma_{\omega,u})$,
$\mbox{\boldmath$\mu$}_{\omega,k}\sim\hbox{Normal}(\mbox{\boldmath$\mu$}_{0},\Sigma_{0})$,
$\Sigma_{\omega,k}\sim\mbox{Inv-W}(\nu_{0},\Psi_{0,\omega})$ constrained to
$\sum^{K_{\epsilon}}_{k=1}\pi_{\omega,k}\mbox{\boldmath$\mu$}_{\omega,k}={\mathbf{0}}$.
For the cluster probability, we further assume
$\mbox{\boldmath$\pi$}_{\omega}=(\pi_{1},\cdots,\pi_{K_{\omega}})^{\rm
T}\sim\hbox{Dirch}(\alpha_{u}/K_{\omega},\cdots,\alpha_{u}/K_{\omega})$.
Additionally, we assume the following distribution for
$\widetilde{\mbox{\boldmath$X$}}_{i}\sim\sum^{K_{x}}_{k=1}\pi_{x,k}\hbox{MVN}_{K}(\mbox{\boldmath$\mu$}_{x,k},\Sigma_{x,k})$,
with
$\mbox{\boldmath$\mu$}_{k}\sim\hbox{MVN}(\mbox{\boldmath$\mu$}_{0},\Sigma_{0}),\quad\Sigma_{x,k}\sim\mbox{Inv-W}(\nu_{0},\Psi_{0})$;
where $\mbox{\boldmath$\pi$}_{x}=(\pi_{1},\cdots,\pi_{K_{x}})^{\rm
T}\sim\hbox{Dirch}(\alpha_{x}/K_{x},\cdots,\alpha_{x}/K_{x})$. We employ a
flat prior $p(\mbox{\boldmath$\beta$}_{z})\propto 1$. Finally, we use the
following smoothing prior for
$\mbox{\boldmath$\gamma$}\propto\exp\left\\{-\frac{1}{2\tau}\mbox{\boldmath$\gamma$}^{\rm
T}\mbox{\boldmath$P$}_{\gamma}\mbox{\boldmath$\gamma$}\right\\}$ and
$\tau\sim\hbox{Inv-Ga}(\alpha_{0},\beta_{0})$, where $P$ is the second order
difference penalty matrix similar to what was proposed in Lang and Brezger
(2004). Note that $\tau$ represents a smoothing parameter where smaller values
of $\tau$ will lead to stronger smoothing. The choice of a penalized prior on
$\gamma$ eliminates the need to precisely choose the number basis $K$ in the
model since unnecessary basis functions will be heavily penalized and their
coefficients driven to zero (Ruppert et al., 2003; Lang and Brezger, 2004).

###  2.3 Posterior

With our choice of prior distributions described above, sampling from the
joint posterior distribution can be done using a sequence of straightforward
Gibbs steps with the full conditional posterior distribution completely
defined and easy to simulate from. We provide the form of the full conditional
posterior distributions in section 3 of the Online supplementary material
along with the choice of prior hyper-parameters and starting values for the
MCMC. A software in the form of R code, together with a sample input data set
and complete documentation will be available on https://github.com/rszoh/SoFR-
ME. Our Bayesian model makes many significant advantages over the non-Bayesian
approach. All inferences about $\gamma$ (hence $\mbox{\boldmath$\beta$}(s)$)
can be done simply based on the draws from the posterior distributions and no
additional bootstrap approach is needed. Our Bayesian model also provides a
way to extract potential subgroups that differ in term of their true
unobserved (latent) $X$ function value since the truncated Dirichlet process
prior amount to performing clustering Li et al., 2019.

###  2.4 Estimation of $\mbox{\boldmath$\delta$}(s)$

Although in our discussion above, we assumed that
$\mbox{\boldmath$\delta$}(t)$ is known. In application $\delta$ in unknown and
should be estimated from the data. Here we proposed a simple approach to
estimating $\delta$. Note base on Equations 4 and 5, we have
$E(W_{i}(s))=E(X_{i}(s))$ and
$E(M_{i}(s))=\mbox{\boldmath$\delta$}(s)E(X_{i}(s))$. Therefore ,
$\mbox{\boldmath$\delta$}(s)=E(M_{i}(s))/E(W_{i}(s))$ provided
$E(X_{i}(s))\neq 0$. So an estimate of $\mbox{\boldmath$\delta$}(s)$, denoted
$\widehat{\mbox{\boldmath$\delta$}}(s)$ is obtained as

$\displaystyle\widehat{\mbox{\boldmath$\delta$}}(s)$ $\displaystyle=$
$\displaystyle\frac{\sum^{n}_{i=1}M_{i}(s)}{\sum^{n}_{i=1}W_{i}(s)}$ (11)

If we further assume that $\mbox{\boldmath$\delta$}(s)$ is a smooth function,
we can further smooth the point-wise estimates obtained in Equation 11 using
any of the many local polynomial smoothing approaches available in R.

##  3 Simulation & Results

###  3.1 Simulation Set-up

We use the following set up for our simulation. We simulate, for each
$i=1,\cdots,n$ independently, $X_{i}(t)$ from a Gaussian process with
$E\\{X_{i}(t)\\}=\sin(2\pi t)$, $\hbox{var}\\{X_{i}(t)\\}=\sigma^{2}_{x}$, and
$\hbox{cov}\\{X_{i}(t_{l}),X_{i}(t_{j})\\}=\rho_{x}\sigma^{2}_{x}$; simulate
$W_{i}(t)=X_{i}(t)+U_{i}(t)$, where $U(t)$ is a Gaussian process with
$E\\{U(t)\\}=0$, $\hbox{var}\\{U_{i}(t)\\}=\sigma^{2}_{u}$, and
$\hbox{cov}\\{U_{i}(t_{l}),U(t_{j})\\}=\rho_{u}\sigma^{2}_{u}$;
$M_{i}(t)=\mbox{\boldmath$\delta$}(t)X_{i}(t)+\omega_{i}(t)$ where
$\omega_{i}(t)$ is a Gaussian process with $E\\{\omega_{i}(t)\\}=0$,
$\hbox{var}\\{\omega_{i}(t)\\}=\sigma^{2}_{\omega}$, and
$\hbox{cov}\\{\omega(t_{l}),\omega(t_{j})\\}=\rho_{\omega}\sigma^{2}_{\omega}$.
Finally, we simulate the response for each unit as
$Y_{i}=\mbox{\boldmath$Z$}_{i}^{\rm
T}\mbox{\boldmath$\beta$}_{z}+\int^{1}_{0}\mbox{\boldmath$\beta$}(t)X_{i}(t)dt+\epsilon_{i}$,
where $\epsilon_{i}\sim\mbox{\boldmath$G$}$, for some distribution $G$. The
data is simulated assuming $\mbox{\boldmath$\beta$}_{z}={\mathbf{0}}$. In each
simulation exercise, we will assess the performance of each estimator of
$\mbox{\boldmath$\beta$}(t)$ using the mean squared integrated error (MSIE)
which equals the sum of the average bias squared ($ABias^{2}$) and the average
sample variance (AVAR) and are expressed respectively as:

$\displaystyle ABias^{2}(\widehat{\beta})$ $\displaystyle=$
$\displaystyle\frac{1}{n_{grid}}\sum^{n_{grid}}_{l=1}\left\\{\bar{\beta}(t_{l})-\beta(t_{l})\right\\}^{2}\quad\mbox{and}$
(12) $\displaystyle Avar(\widehat{\beta})$ $\displaystyle=$
$\displaystyle\frac{1}{n_{s}n_{grid}}\sum^{n_{r}}_{i=1}\sum^{n_{grid}}_{j=1}\left\\{\widehat{\beta}(t_{j})-\bar{\beta}(t_{j})\right\\}^{2}$
(13) $\displaystyle MSIE$ $\displaystyle=$ $\displaystyle ABias^{2}+Avar,$
(14)

where $n_{grid}$ represents the number of equally selected grid points between
0 and 1;
$\bar{\beta}(t_{j})=\frac{1}{n_{r}}\sum^{n_{r}}_{j=1}\widehat{\beta}^{r}(t_{j})$
represents the point-wise average, averaged over the $n_{r}$ replicates of
$\beta(t_{j})$ at a specific time point $t_{j}$.

In what follows, we denote our approach as
$\widehat{\mbox{\boldmath$\beta$}}_{BIV}$ (Bayesian approach using IV),
$\widehat{\mbox{\boldmath$\beta$}}_{IV,\delta}$ (scaling $M$ by
$\mbox{\boldmath$\delta$}(t)$) and $\widehat{\mbox{\boldmath$\beta$}}_{IV}$
(not scaling $M$) are both variants of the approach by Tekwe et al., 2019;
Finally, $\widehat{\mbox{\boldmath$\beta$}}_{ReF.W}$ denotes the estimate of
$\mbox{\boldmath$\beta$}(t)$ using the function pfr from the R package refund
(Goldsmith et al., 2019) using $W$ as the functional covariate using the
penalized B-splines options. We described in details the options we used for
the pfr function in our simulation in the Supplementary Material. Results are
based on $N=600$ replications. Finally, we choose
$\mbox{\boldmath$\delta$}(t)=\frac{\delta}{2}(1+\sin(2\pi
t))+\min(\delta,0.02)$.

###  3.2 Simulation Results

Our simulation studies are divided into four parts. In the part 1, we compare
the performance of all competing approaches. Data are simulated assuming
$\sigma_{x}=4.0$, $\sigma_{u}=4.0$, $\sigma_{e}=1.0$, $\sigma_{\omega}=1.0$,
$\rho_{x}=\rho_{u}=\rho_{\omega}=0.25$, with $n=500$ and $1000$. The
functional covariate $X(t)$ is evaluated at $t=50$ different but equidistant
time points on the $(0,1)$ interval. We simulate the error
$\epsilon_{i}\stackrel{{\scriptstyle
iid}}{{\sim}}\hbox{Normal}(0,\sigma^{2}_{e})$ and $\sigma_{e}=1$. The results
are provided in Table 1. We find that $\mbox{\boldmath$\beta$}_{BIV}$,
$\mbox{\boldmath$\beta$}_{IV}$, and $\mbox{\boldmath$\beta$}_{IV,\delta}$
perform very similarly and tended to have very low bias. However,
$\mbox{\boldmath$\beta$}_{ReF.X}$ has a very high biased and overall large
$MISE$ compared to the others. It is important to note that
$\mbox{\boldmath$\beta$}_{ReF.X}$ only uses $W$ and this might clearly explain
the low performance. For the remaining of the paper, we will drop
$\mbox{\boldmath$\beta$}_{ReF.X}$ from of our discussion since we believe.

In the second part of the simulation, we investigated the impact of various
levels of the variances $\sigma_{x},\sigma_{u},\sigma_{\omega}$ on the
estimates of $\mbox{\boldmath$\beta$}(t)$. Varying values of
$\sigma_{x},\sigma_{u},\sigma=0.5,1.0,16$ were considered. In each setting
only one parameter is chancged and the others are held constant at the value
of 1.0. The data were simulated assuming $\delta=2$, $\sigma_{e}=1.0$,
$\rho_{x}=\rho_{u}=\rho_{\omega}=0.25$, with $n=500$. We report the results in
Table 2. For all approaches, we note that the estimated MSIE increases with
$\sigma_{u}$ and $\sigma_{\omega}$ but decreases with increasing $\sigma_{x}$.
However $\mbox{\boldmath$\beta$}_{BIV}$ tended to have the lowest MSIE in all
cases. We also note that all approaches tended to have higher bias with
increasing $\sigma_{u}$ and $\sigma_{\omega}$ although the increase in the
bias seems larger in the $\mbox{\boldmath$\beta$}_{BIV}$ compared to the
others.

In the third part of our simulations, we considered varying level of
dependency in the error processes for $X(s)$, $U(s)$, and $\omega(s)$ and
assessed their impact the estimate of $\mbox{\boldmath$\beta$}(s)$. The
correlation coefficient values considered for $\rho_{x};\rho_{u};\rho_{w}$
were $\\{0,0.25,0.75\\}$. The data were simulated assuming $\sigma_{x}=4.0$,
$\sigma_{u}=4.0$, $\sigma_{e}=1.0$, $\sigma_{\omega}=1.0$ and $n=500$. The
results are reported in Table 3. Overall, we note that MSIE tended to be
similar for varying values of the correlations values for $\rho_{u};\rho_{w}$.
However, lower values of $\rho_{x}$ ($\rho_{x}=0$) tended to have higher MSIE.
But our approach ($\mbox{\boldmath$\beta$}_{BIV}$) tended to have smaller MISE
in all cases too. So far in our discussion, we have observed that the Bayesian
estimate ($\mbox{\boldmath$\beta$}_{BIV}$) seems to outperform its competitors
($\mbox{\boldmath$\beta$}_{IV},\mbox{\boldmath$\beta$}_{IV,\delta}$). We also
note that the original estimate of Tekwe et al. (2019) and its
$\mbox{\boldmath$\delta$}(t)$ adjusted version do not differ in term of MISE,
at least for the choice of $\mbox{\boldmath$\delta$}(t)$ we considered so far.
In the fourth part of simulation, we investigate the performance of our
approach with respect to varying values of $\mbox{\boldmath$\delta$}(t)$(i.e.,
$\delta$). Note that in the function $\mbox{\boldmath$\delta$}(t)$, $\delta$
is a scaling factor and large value of $\delta$ are associated with wider
range of the function. Result of the simulation are reported Table 4. Again,
data are simulated assuming $\sigma_{x}=4.0$, $\sigma_{u}=4.0$,
$\sigma_{e}=1.0$, $\sigma_{\omega}=1.0$ and $n=500$,
$\rho_{x};\rho_{u};\rho_{w}=0.25$. We note again that
$\mbox{\boldmath$\beta$}_{BIV}$ has overall smaller MSIE and bias across all
versions of the function $\mbox{\boldmath$\delta$}(t)$ (varying $\delta$)
considered. However, both
$\mbox{\boldmath$\beta$}_{IV},\mbox{\boldmath$\beta$}_{IV,\delta}$ perform
virtually the same albeit the $\mbox{\boldmath$\delta$}(t)$ corrected
($\mbox{\boldmath$\beta$}_{IV,\delta}$) seems to have smaller MSIE compared to
the original approach $\mbox{\boldmath$\beta$}_{IV}$.

##  4 Application

We now discuss the application of our proposed methods to the motivating data.
The Center for Disease Control established the NHANES program in the 1960s to
describe the health of Americans. This cross-sectional survey employs a
complex multi-stage cluster probability sampling scheme and assigns
participants sampling weights so that the resulting data represent the U.S.
population as a whole. The data are based on in-person interviews and physical
examinations. Approximately $5,000$ participants from 15 geographic regions
are examined annually. Data collected include demographics, socioeconomic
parameters, dietary data, health behaviors, medical conditions, physiological
examination results, and laboratory assessments of blood samples. We analyzed
data on $1893$ NHANES participants, 48.5% women, median age = 43 years, with
an age range of 20–69 years who have accelerator data for four or more days
from the 2005-2006 cycles of the survey. Participants included in the
application wore the wearable devices for physical activity around their hips
(ActiGraph AM-7164) to track their daily physical activity for at least four
days during waking hours (Matthews et al., 2008; Troiano et al., 2008).
However, the device was taken off during water activities and sleep. The
device recorded data on physical activity intensity and step counts every 60
seconds while worn. True physical activity intensity is considered the true
function-valued covariate, $X(t)$. The device-based measure of true physical
activity intensity or the observed physical activity intensity is considered
$W(t)$ in this application. While the device-based measures of step counts is
considered the instrumental variable, $M(t)$. We assumed that conditional on
true physical activity intensity, step counts is independent of the
measurement error. We also assumed that the effects of step count on the
device-based measure of physical activity intensity is only through true
physical activity intensity. The error free covariates included in the
analyses were age, sex, and ethnicity. Sample weights were applied to the data
to account for the over sampling of racial groups using analytic guidelines
provided for the NHANES (Johnson et al., 2013). The device-based data were
grouped in $30$-minutes (half-hour) intervals and we averaged over $30$
minutes across the multiple days of wear time. We standardized the data by
considering our baseline $t=0$ to be the beginning time at which the device is
placed on, while $t=24$ represents the end of the day. We show the plots of
the data: half-hour average step counts, half-hour average intensity, and the
value of BMI in Supplementary Figure $S1$. We apply a $\hbox{log}2$
transformation on $M$ and $W$. Then we obtain an estimate
$\hat{\mbox{\boldmath$\delta$}}(t)$ as discussed above. Using the default
priors described in Subsection 2.2, we ran our MCMC for $10,000$ iterations,
discarded $2000$ as burn-in, and selected values of $3$ for thinning. The
remaining MCMC draws were then used for the final inference. Additional MCMC
convergence diagnosis tools were used but omitted here for brevity. We show
the plot of $\hbox{log}2$ transformed data for physical activity intensity
($\mbox{\boldmath$W$}(t)$) and steps counts $\mbox{\boldmath$W$}(t)$ in Figure
1a and 1b along with the mean curved of the two groups of physical activity
obtained from the posterior.

The estimated functional parameter $\mbox{\boldmath$\beta$}(t)$ under each of
the approaches ($\mbox{\boldmath$\beta$}_{BIV}$,
$\mbox{\boldmath$\beta$}_{IV}$ described in the simulations are illustrated in
Figure 1a. Note that for $\mbox{\boldmath$\beta$}_{IV}$, we obtain estimate of
$\mbox{\boldmath$\beta$}(t)$ based on the residuals obtained after we regress
$Y$ (BMI) on the error free covariates $Z$.

We report the $90\%$ point-wise credible intervals along with the posterior
mean estimate of $\mbox{\boldmath$\beta$}(t)$ based on the posterior draws.
The posterior mean of $\mbox{\boldmath$\beta$}(t)$ tends to increase with
decreasing amount of the physical activity and decrease with with increase
physical activity (increased step counts and physical activity intensity). We
also observe that $\mbox{\boldmath$\beta$}(t)$ is significant for peak hour
physical intensity (between 12 and 22 hrs). This is followed by an increase of
the $\mbox{\boldmath$\beta$}(t)$ function at later hours. We find that these
results are consistent with the belief that increasing physical activities are
associated with lower BMI (refs). The estimate of
$\mbox{\boldmath$\beta$}_{IV}(t)$ seems to closely resemble our estimate
$\mbox{\boldmath$\beta$}_{BIV}(t)$ but seems slightly attenuated towards the
zero vertical line.

In addition to estimating the function parameter $\mbox{\boldmath$\beta$}(t)$,
it might be of interest to extract potential individual groups that differ in
term of their true (latent) physical activity pattern. Our Bayesian semi-
parametric model allows one to infer the potential grouping of the true latent
physical activity pattern along with the group mean functions. Base on our
posterior distribution of the group membership indicator, we estimate two main
physical activity pattern along with their respective mean functions (see Fig
(b)). The size of these groups are: 1653 (cluster 1), 0 (cluster 2), and 240
(cluster 3). We estimate the posterior mean difference BMI between two similar
individuals (race and age) in cluster 1 and cluster 3 based on their mean
cluster physical activity as
$\int^{1}_{0}\beta(t)\\{\hat{\mbox{\boldmath$X$}}_{\text{cluster
1}}-\hat{\mbox{\boldmath$X$}}_{\text{\text{cluster 3} }}\\}dt$ which we
estimated to be $-0.958\;(-1.560,-0.354)$ along with a $90\%$ posterior
interval for the difference - which does not contain zero. This suggests that
these two clusters are significantly different in term of their effect on BMI.
Also, we also report the estimates for the error free covariates (see Table )

##  5 Conclusion

We proposed a fully Bayesian approach to estimate the functional parameter in
the SoFR in the presence of error free covariates and a sparsely observed
functional covariate contaminated by measurement error. We make use of a
functional instrumental variable to estimate the measurement error covariance.
Through extensive simulations, we demonstrated that overall, our proposed
estimator outperformed its competitors. Our model simultaneously allows for
the estimation of a global instrumental variable quality parameter, $\delta$,
which may indicate the quality of the estimator of
$\mbox{\boldmath$\beta$}(t)$. Note also that $\delta$ allows for the use of an
IV that has potentially a different scale than the functional covariate
$X(t)$. This confers more flexibility in the choice of an IV. We believe that
our work, although the first to provide a fully Bayesian approach to a
measurement error SoFR with an instrumental variable, can be substantially
extended. The assumption of a single global quality parameter for the
functional instrumental variable, $\mbox{\boldmath$M$}(t)$, could be relaxed
to allow for a functional quality parameter. Additionally, the MCMC proposed
here, although fast can be also be improved using a tailored variational
approach. These points will be the subject of future work.

Figure 1: Plot the analysis summary. (a) Show the estimated
$\mbox{\boldmath$\beta$}(t)$ for all the participants (n=1893)
$\mbox{\boldmath$\beta$}_{BIV}$ (continuous line),
$\mbox{\boldmath$\beta$}_{ReF.W}$ (longdash line), and
$\mbox{\boldmath$\beta$}_{IV}$(in dot dash line). (b) Plot the log physical
intensity per cluster. More active (continuous line) and less active (dashed
line). (c) Estimated $\mbox{\boldmath$\beta$}(t)$ using the less-active group
data (n=1400). $\mbox{\boldmath$\beta$}_{BIV}$ (continuous line),
$\mbox{\boldmath$\beta$}_{ReF.W}$ (longdash line), and
$\mbox{\boldmath$\beta$}_{IV}$(in dot dash line). (d) Estimated
$\mbox{\boldmath$\beta$}(t)$ using the more active group data (n=493).
$\mbox{\boldmath$\beta$}_{BIV}$ (continuous line),
$\mbox{\boldmath$\beta$}_{ReF.W}$ (longdash line), and
$\mbox{\boldmath$\beta$}_{IV}$(in dot dash line).

(a) $\hbox{log}2$ Physical activity Intensity

(b) $\hbox{log}2$ Steps Counts

(c) Estimated $\mbox{\boldmath$\beta$}(t)$

(d) Posterior Mean of the latent Physical activity for each of the 3 latent
groups identified.

Table 1: Estimated $MISE$ based on $500$ replicates. Data are simulated assuming no effect for the error-free covariates. Method | n | Bias2 | Var | MISE
---|---|---|---|---
$\mbox{\boldmath$\beta$}_{BIV}$ | 500 | 0.0006 | 0.0007 | 0.0013
1000 | 0.0004 | 0.0003 | 0.0007
$\mbox{\boldmath$\beta$}_{IV,\delta}$ | 500 | 0.0051 | 0.0007 | 0.0058
1000 | 0.0001 | 0.0003 | 0.0004
$\mbox{\boldmath$\beta$}_{IV}$ | 500 | 0.0061 | 0.0007 | 0.0068
1000 | 0.0001 | 0.0004 | 0.0004
$\mbox{\boldmath$\beta$}_{IReF.W}$ | 500 | 0.3076 | 0.0004 | 0.3080
1000 | 0.2369 | 0.0002 | 0.2371
Table 2:  Estimated MISE for estimating $\mbox{\boldmath$\beta$}(t)$ when data
are simulated for various values of
$\sigma_{\omega},\sigma_{u},\sigma_{x}\in\\{.5,1,16\\}$. Note that the
parameters are changing one at a time, while the other parameters are held
constant at $1$.

Method | $\sigma_{u}$ | Bias2 | Var | MSIE | $\sigma_{w}$ | Bias2 | Var | MSIE | $\sigma_{x}$ | Bias2 | Var | MSIE
---|---|---|---|---|---|---|---|---|---|---|---|---
$\mbox{\boldmath$\beta$}_{BIV}$ | 0.5 | 0.0002 | 0.0006 | 0.0008 | 0.5 | 0.0004 | 0.0006 | 0.0011 | 0.5 | 0.0208 | 0.0396 | 0.0604
1.0 | 0.0002 | 0.0006 | 0.0008 | 1.0 | 0.0002 | 0.0006 | 0.0008 | 1.0 | 0.0026 | 0.0182 | 0.0208
16 | 0.0036 | 0.0031 | 0.0067 | 16 | 0.0041 | 0.0043 | 0.0083 | 16 | 0.0001 | 0.0001 | 0.0003
$\mbox{\boldmath$\beta$}_{IV,\delta}$ | 0.5 | 0.0051 | 0.0005 | 0.0056 | 0.5 | 0.0053 | 0.0005 | 0.0058 | 0.5 | 0.0536 | 0.0111 | 0.0647
1.0 | 0.0051 | 0.0005 | 0.0056 | 1.0 | 0.0051 | 0.0005 | 0.0056 | 1.0 | 0.0096 | 0.0103 | 0.0199
16 | 0.0066 | 0.0019 | 0.0085 | 16 | 0.0058 | 0.0064 | 0.0122 | 16 | 0.0051 | 0.0000 | 0.0051
$\mbox{\boldmath$\beta$}_{IV}$ | 0.5 | 0.0057 | 0.0005 | 0.0062 | 0.5 | 0.0058 | 0.0005 | 0.0064 | 0.5 | 0.0524 | 0.0117 | 0.0641
1.0 | 0.0064 | 0.0005 | 0.0069 | 1 | 0.0062 | 0.0005 | 0.0067 | 1.0 | 0.0124 | 0.0102 | 0.0226
16 | 0.0074 | 0.0019 | 0.0093 | 16 | 0.0066 | 0.0069 | 0.0135 | 16 | 0.0062 | 0.0000 | 0.0062

Table 3: Estimated MISE when data simulated assuming dependent error terms
$\rho_{u},\rho_{\omega},\rho_{x}\in 0.0,0.25,0.5,0.75$. Note that for varying
correlation, the other correlations are held at $0.25$.

Method | $\rho_{u}$ | Bias2 | Var | MSIE | $\rho_{\omega}$ | Bias2 | Var | MSIE | $\rho_{x}$ | Bias2 | Var | MSIE
---|---|---|---|---|---|---|---|---|---|---|---|---
$\mbox{\boldmath$\beta$}_{BIV}$ | 0.0 | 0.0013 | 0.0006 | 0.0019 | 0.0000 | 0.0007 | 0.0007 | 0.0014 | 0.0 | 0.0006 | 0.0135 | 0.0141
0.25 | 0.0008 | 0.0007 | 0.0016 | 0.25 | 0.0007 | 0.0006 | 0.0014 | 0.25 | 0.0006 | 0.0006 | 0.0013
0.75 | 0.0003 | 0.0008 | 0.0010 | 0.75 | 0.0009 | 0.0008 | 0.0017 | 0.75 | 0.0010 | 0.0003 | 0.0013
$\mbox{\boldmath$\beta$}_{BIV,\delta}$ | 0.0 | 0.0052 | 0.0006 | 0.0058 | 0.0 | 0.0053 | 0.0007 | 0.0060 | 0.0 | 0.0068 | 0.0093 | 0.0160
0.25 | 0.0052 | 0.0008 | 0.0059 | 0.25 | 0.0052 | 0.0007 | 0.0059 | 0.25 | 0.0052 | 0.0006 | 0.0058
0.75 | 0.0052 | 0.0009 | 0.0061 | 0.75 | 0.0051 | 0.0008 | 0.0059 | 0.75 | 0.0052 | 0.0002 | 0.0055
$\mbox{\boldmath$\beta$}_{BIV}$ | 0.0 | 0.0058 | 0.0006 | 0.0064 | 0.0 | 0.0062 | 0.0008 | 0.0069 | 0.0 | 0.0085 | 0.0107 | 0.0192
0.25 | 0.0061 | 0.0008 | 0.0068 | 0.25 | 0.0066 | 0.0007 | 0.0073 | 0.25 | 0.0061 | 0.0006 | 0.0067
0.75 | 0.0061 | 0.0009 | 0.0070 | 0.75 | 0.0063 | 0.0008 | 0.0071 | 0.75 | 0.0064 | 0.0002 | 0.0067

Table 4: Estimated MSIE for various values for various values of the function $\mbox{\boldmath$\delta$}(t)$. Note that increasing values of $\delta$ are associated with increasing range for the function $\mbox{\boldmath$\delta$}(t)$. Method | $\delta$ | ABias2 | AVar | MSIE
---|---|---|---|---
$\mbox{\boldmath$\beta$}_{BIV}$ | 0.005 | 0.0004 | 0.0008 | 0.0012
0.05 | 0.0005 | 0.0009 | 0.0014
5 | 0.0006 | 0.0008 | 0.0014
| 20 | 0.0005 | 0.0008 | 0.0014
$\mbox{\boldmath$\beta$}_{IV,\delta}$ | 0.005 | 0.0052 | 0.0008 | 0.0060
0.05 | 0.0053 | 0.0008 | 0.0061
5 | 0.0052 | 0.0008 | 0.0060
| 20 | 0.0052 | 0.0009 | 0.0060
$\mbox{\boldmath$\beta$}_{IV}$ | 0.0050 | 0.0053 | 0.0008 | 0.0062
0.05 | 0.0068 | 0.0008 | 0.0076
5 | 0.0060 | 0.0009 | 0.0069
| 20 | 0.0064 | 0.0009 | 0.0073
Table 5: Posterior summary estimates for the error free covariates using the
Bayesian Instrumental variable approach. Age was median centered by
subtracting the median age (43 years).

| Non-His_WHite | Black | Hispanic | Asian | Others | Gender(Male) | Age (Median-43)
---|---|---|---|---|---|---|---
Mean | 31.65 | 31.77 | 32.23 | 33.04 | 30.03 | 0.06 | 0.08
5% | 29.80 | 29.75 | 29.93 | 31.08 | 27.84 | -0.45 | 0.05
95% | 33.33 | 33.67 | 34.41 | 34.82 | 32.10 | 0.58 | 0.10

## acknowledgement

Zoh’s research was supported by National Cancer Institute Supplemental under
award U01-CA057030-29S1. Tekwe’s research was supported by National Cancer
Institute Supplemental under award U01-CA057030-29S2. This research was also
supported in part by Lilly Endowment, Inc., through its support for the
Indiana University Pervasive Technology Institute.

Conflict of Interest

The authors have declared no conflict of interest.

## References

  * An et al. (2017) Hyun-Sung An, Gregory C Jones, Seoung-Ki Kang, Gregory J Welk, and Jung-Min Lee.  How valid are wearable physical activity trackers for measuring steps?  _European Journal of Sport Science_ , 17(3):360–368, 2017.
  * Cardot et al. (2007) Herve Cardot, Christophe Crambes, Alois Kneip, and Pascal Sarda.  Smoothing splines estimators in functional linear regression with errors-in-variables.  _Computational Statistics & Data Analysis_, 51(10):4832–4848, 2007.
  * Carroll et al. (2006) R. J. Carroll, D. Ruppert, L. A. Stefanski, and C. M. Crainiceanu.  _Measurement Error in Nonlinear Models: A Modern Perspective, Second Edition_.  Chapman and Hall, 2006.
  * Chakraborty and Panaretos (2017) Anirvan Chakraborty and Victor M Panaretos.  Regression with genuinely functional errors-in-covariates.  _arXiv preprint arXiv:1712.04290_ , 2017.
  * Feito et al. (2012) Yuri Feito, David R Bassett, Dixie L Thompson, and Brian M Tyo.  Effects of body mass index on step count accuracy of physical activity monitors.  _Journal of Physical Activity and Health_ , 9(4):594–600, 2012.
  * Goldsmith et al. (2013) Jeff Goldsmith, Sonja Greven, and CIPRIAN Crainiceanu.  Corrected confidence bands for functional data using principal components.  _Biometrics_ , 69(1):41–51, 2013.
  * Goldsmith et al. (2019) Jeff Goldsmith, Fabian Scheipl, Lei Huang, Julia Wrobel, Chongzhi Di, Jonathan Gellar, Jaroslaw Harezlak, Mathew W. McLean, Bruce Swihart, Luo Xiao, Ciprian Crainiceanu, and Philip T. Reiss.  _refund: Regression with Functional Data_ , 2019.  URL https://CRAN.R-project.org/package=refund.  R package version 0.1-21.
  * Hastie et al. (2015) Trevor Hastie, Robert Tibshirani, and Martin Wainwright.  _Statistical learning with sparsity: the lasso and generalizations_.  CRC press, 2015.
  * Johnson et al. (2013) CL Johnson, R Paulose-Ram, CL Ogden, MD Carroll, D Kruszon-Moran, SM Dohrmann, and LR Curtin.  National health and nutrition examination survey: analytic guidelines, 1999-2010.  _Vital and health statistics. Series 2, Data evaluation and methods research_ , (161):1, 2013.
  * Lang and Brezger (2004) Stefan Lang and Andreas Brezger.  Bayesian p-splines.  _Journal of computational and graphical statistics_ , 13(1):183–212, 2004.
  * Li et al. (2019) Yuelin Li, Elizabeth Schofield, and Mithat Gönen.  A tutorial on dirichlet process mixture modeling.  _Journal of mathematical psychology_ , 91:128–144, 2019\.
  * Matthews et al. (2008) Charles E Matthews, Kong Y Chen, Patty S Freedson, Maciej S Buchowski, Bettina M Beech, Russell R Pate, and Richard P Troiano.  Amount of time spent in sedentary behaviors in the ’united states’, 2003–2004.  _American Journal of Epidemiology_ , 167(7):875–881, 2008.
  * Morris (2015) Jeffrey S Morris.  Functional regression.  _Annual Review of Statistics and Its Application_ , 2:321–359, 2015.
  * O’Driscoll et al. (2020) Ruairi O’Driscoll, Jake Turicchi, Kristine Beaulieu, Sarah Scott, Jamie Matu, Kevin Deighton, Graham Finlayson, and James Stubbs.  How well do activity monitors estimate energy expenditure? a systematic review and meta-analysis of the validity of current technologies.  _British Journal of Sports Medicine_ , 54(6):332–340, 2020.
  * Ramsay (2006) James O Ramsay.  _Functional data analysis_.  Wiley Online Library, 2006.
  * Reiss et al. (2017) Philip T Reiss, Jeff Goldsmith, Han Lin Shang, and R Todd Ogden.  Methods for scalar-on-function regression.  _International Statistical Review_ , 85(2):228–249, 2017.
  * Ruppert et al. (2003) D. Ruppert, M. A. Wand, and R. J. Carroll.  _Semiparametric Regression_.  Cambridge University Press, 2003.
  * Sarkar et al. (2018) Abhra Sarkar, Debdeep Pati, Antik Chakraborty, Bani K Mallick, and Raymond J Carroll.  Bayesian semiparametric multivariate density deconvolution.  _Journal of the American Statistical Association_ , 113(521):401–416, 2018.
  * Tekwe et al. (2019) Carmen D Tekwe, Roger S Zoh, Miao Yang, Raymond J Carroll, Gilson Honvoh, David B Allison, Mark Benden, and Lan Xue.  Instrumental variable approach to estimating the scalar-on-function regression model with measurement error with application to energy expenditure assessment in childhood obesity.  _Statistics in medicine_ , 38(20):3764–3781, 2019\.
  * Troiano et al. (2008) Richard P Troiano, David Berrigan, Kevin W Dodd, Louise C Mâsse, Timothy Tilert, Margaret McDowell, et al.  Physical activity in the United States measured by accelerometer.  _Medicine and Science in Sports and Exercise_ , 40(1):181, 2008.
  * Van Erp et al. (2019) Sara Van Erp, Daniel L Oberski, and Joris Mulder.  Shrinkage priors for bayesian penalized regression.  _Journal of Mathematical Psychology_ , 89:31–50, 2019.
  * Wang et al. (2016) Jane-Ling Wang, Jeng-Min Chiou, and Hans-Georg Müller.  Functional data analysis.  _Annual Review of Statistics and Its Application_ , 3:257–295, 2016.
  * Yao et al. (2005) Fang Yao, Hans-Georg Müller, and Jane-Ling Wang.  Functional data analysis for sparse longitudinal data.  _Journal of the American statistical association_ , 100(470):577–590, 2005.


on $p$. Then, for each $i$,

$S(p_{i},r_{i})-S(q_{i},r_{i})=\left[\tilde{S}(p_{i},r_{i})+d(r_{i})\right]-\left[\tilde{S}(q_{i},r_{i})-d(r_{i})\right]=\tilde{S}(p_{i},r_{i})-\tilde{S}(q_{i},r_{i}).$
(39)

Combining (38) and (39) gives us the result.

###  B.2 Proof of Theorem 2

The proof closely resembles the proof of Theorem 4 in Howard et al., (2021).
Fix $c=1$ without loss of generality, so that $\hat{\delta}_{i}$,
$\delta_{i}$, and $\gamma_{i}$ are each bounded by
$\left[-\frac{1}{2},\frac{1}{2}\right]$ and consequently the differences
$\hat{\delta}_{i}-\delta_{i}$ and $\hat{\delta}_{i}-\gamma_{i}$ are bounded by
$\left[-1,1\right]$.

Our goal is to show that

$L_{t}(\lambda)=\exp\left\\{\lambda
S_{t}-\psi_{E}(\lambda)\sum_{i=1}^{t}\left(\hat{\delta}_{i}-\gamma_{i}\right)^{2}\right\\}$
(40)

is a supermartingale for each $\lambda\in[0,1)$, where
$\psi_{E}(\lambda)=\psi_{E,1}(\lambda)=-\log(1-\lambda)-\lambda$.

Let $Y_{i}=\hat{\delta}_{i}-\delta_{i}$ and $Z_{i}=\gamma_{i}-\delta_{i}$,
such that $Y_{i}-Z_{i}=\hat{\delta}_{i}-\gamma_{i}\in[-1,1]$. Using the fact
that $\exp\left\\{\lambda\xi-\psi_{E}(\lambda)\xi^{2}\right\\}\leq
1+\lambda\xi$ for all $\lambda\in[0,1)$ and $\xi\geq-1$ (Fan et al.,, 2015),
we have

$\exp\left\\{\lambda\left(Y_{i}-Z_{i}\right)-\psi_{E}(\lambda)\left(Y_{i}-Z_{i}\right)^{2}\right\\}\leq
1+\lambda\left(Y_{i}-Z_{i}\right).$

Multiplying each side by $e^{\lambda Z_{i}}$, we get

$\displaystyle\exp\left\\{\lambda
Y_{i}-\psi_{E}(\lambda)\left(Y_{i}-Z_{i}\right)^{2}\right\\}$
$\displaystyle\leq e^{\lambda Z_{i}}\left(1-\lambda Z_{i}\right)+\lambda
e^{\lambda Z_{i}}Y_{i}.$

Now, conditional on $\mathcal{G}_{i-1}$, we know that
$\operatorname*{\mathbb{E}}\left[{Y_{i}}\;\middle|\;{\mathcal{G}_{i-1}}\right]=\operatorname*{\mathbb{E}}\left[{\hat{\delta}_{i}-\delta_{i}}\;\middle|\;{\mathcal{G}_{i-1}}\right]=0$
by the first condition. We also know that $(Z_{i})_{i=1}^{\infty}$ is
predictable w.r.t. $(\mathcal{G}_{i})_{i=1}^{\infty}$, as
$(\gamma_{i})_{i=1}^{\infty}$, $(p_{i})_{i=1}^{\infty}$,
$(q_{i})_{i=1}^{\infty}$, $(r_{i})_{i=1}^{\infty}$ are all predictable w.r.t.
$(\mathcal{G}_{i})_{i=1}^{\infty}$. Thus, by taking conditional expectations
on both sides, and substituting back $Y_{i}=\hat{\delta}_{i}-\delta_{i}$ and
$Z_{i}=\gamma_{i}-\delta_{i}$ on the left-hand side, we get

$\displaystyle\operatorname*{\mathbb{E}}\left[{\exp\left\\{\lambda\left(\hat{\delta}_{i}-\delta_{i}\right)-\psi_{E}(\lambda)\left(\hat{\delta}_{i}-\gamma_{i}\right)^{2}\right\\}}\;\middle|\;{\mathcal{G}_{i-1}}\right]$
$\displaystyle\leq e^{\lambda Z_{i}}\left(1-\lambda Z_{i}\right)+\lambda
e^{\lambda
Z_{i}}\operatorname*{\mathbb{E}}\left[{Y_{i}}\;\middle|\;{\mathcal{G}_{i-1}}\right]$
$\displaystyle=e^{\lambda Z_{i}}\left(1-\lambda Z_{i}\right)\leq 1,$

using $1-x\leq e^{-x}$ for the final inequality.

Finally, observe that the left-hand side is precisely the expected increment
of $(L_{t}(\lambda))_{t=1}^{\infty}$ at time $t$. It follows that
$\operatorname*{\mathbb{E}}\left[{L_{t}(\lambda)}\;\middle|\;{\mathcal{G}_{t-1}}\right]\leq
L_{t-1}(\lambda)$, i.e., $(L_{t}(\lambda))_{t=1}^{\infty}$ is a
supermartingale. Thus, by (36) and (37), we have

$\operatorname*{\mathbb{P}}\left({\exists t\geq 1:S_{t}\geq
u(\hat{V}_{t})}\right)\leq\alpha$ (41)

for a sub-$\psi_{E}$ (i.e., sub-exponential) uniform boundary $u$ with
crossing probability $\alpha/2$ and scale $c=1$. Using the fact that
$\frac{1}{t}S_{t}=\frac{1}{t}\sum_{i=1}^{t}\hat{\delta}_{i}-\frac{1}{t}\sum_{i=1}^{t}\delta_{i}=\hat{\Delta}_{t}-\Delta_{t}$,
we can divide each side of the inequality by $t$ to obtain the upper
confidence bound. Finally, by applying the same argument to
$(-S_{t},\hat{V}_{t})$ and combining the upper and lower bounds with a union
bound, we obtain the CS:

$\operatorname*{\mathbb{P}}\left({\forall t\geq
1:\left\lvert\hat{\Delta}_{t}-\Delta_{t}\right\rvert<\frac{u(\hat{V}_{t})}{t}}\right)\geq
1-\alpha.$ (42)

###  B.3 Proof of Theorem 3

Let
$S_{t}=\sum_{i=1}^{t}(\hat{\delta}_{i}-\delta_{i})=t(\hat{\Delta}_{t}-\Delta_{t})$
and $\hat{V}_{t}=\sum_{i=1}^{t}(\hat{\delta}_{i}-\gamma_{i})^{2}$ where
$(\gamma_{i})_{i=1}^{\infty}$ is predictable w.r.t.
$(\mathcal{G}_{i})_{i=1}^{\infty}$. As shown in the proof of Theorem 2
(Appendix B.2), we have that, for any $\lambda\in[0,\lambda_{\mathrm{max}})$,

$L_{t}(\lambda):=\exp\left\\{\lambda
S_{t}-\psi_{E}(\lambda)\hat{V}_{t}\right\\}$ (43)

is a supermartingale w.r.t. $(\mathcal{G}_{t})_{t=1}^{\infty}$, that is,
$\operatorname*{\mathbb{E}}\left[{L_{t}(\lambda)}\;\middle|\;{\mathcal{G}_{t-1}}\right]\leq
L_{t-1}(\lambda)$ for each $t$.

Now, under $H_{0}^{\mathsf{w}}(p,q)$, we have that
$\exp\left\\{-\lambda\sum_{i=1}^{t}\delta_{i}\right\\}\geq 1$, so

$\displaystyle L_{t}(\lambda)$
$\displaystyle=\exp\left\\{\lambda\sum_{i=1}^{t}\hat{\delta}_{i}-\psi_{E}(\lambda)\hat{V}_{t}\right\\}\exp\left\\{-\lambda\sum_{i=1}^{t}\delta_{i}\right\\}$
$\displaystyle\geq\exp\left\\{\lambda\sum_{i=1}^{t}\hat{\delta}_{i}-\psi_{E}(\lambda)\hat{V}_{t}\right\\}=E_{t}(\lambda).$
(44)

Note that $(E_{t}(\lambda))_{t=1}^{\infty}$ is adapted to
$(\mathcal{G}_{t})_{t=1}^{\infty}$. Because $L_{t}(\lambda)$ is a
supermartingale w.r.t. $(\mathcal{G}_{t})_{t=1}^{\infty}$, it then follows
that, under $H_{0}^{\mathsf{w}}$, at any stopping time $\tau$,

$\operatorname*{\mathbb{E}}\left[{E_{\tau}(\lambda)}\right]\leq\operatorname*{\mathbb{E}}\left[{L_{\tau}(\lambda)}\right]\leq
1$ (45)

where the last inequality follows from Doob’s optional stopping theorem
(Doob,, 1940) and that $L_{0}(\lambda)=1$ vacuously.

###  B.4 Proof of Lemma 3

Let $\tilde{S}$ be a linear equivalent of $S$. Define $\tilde{S}_{\mathbf{w}}$
as follows:

$\tilde{S}_{\mathbf{w}}(\mathbf{p}_{t},r_{t}):=\sum_{k=1}^{K}w_{k}\tilde{S}(p_{t}^{(k)},r_{t}).$
(46)

Then, $\tilde{S}_{\mathbf{w}}$ is a linear equivalent of $S_{\mathbf{w}}$, as
for any $\mathbf{p}$ and $r$, $\tilde{S}(\mathbf{p},r)$ is a linear function
of $r$ and
$S(\mathbf{p},r)-\tilde{S}(\mathbf{p},r)=\sum_{k=1}^{K}w_{k}\left[S(p^{(k)},r)-\tilde{S}(q^{(k)},r)\right]$
does not depend on $\mathbf{p}$.

We now show the analogs of (38) and (39) using the linearity of sums in
$S_{\mathbf{w}}$ and $\tilde{S}_{\mathbf{w}}$. First note that, because
$p_{t}^{(k)}$ is a forecast made at the beginning of round $t-k+1$, it only
uses information available to forecasters up to time $t-k$, i.e.,
$p_{t}^{(k)}\in\mathcal{F}_{t-k}$. (For example,
$p_{t}^{(1)}\in\mathcal{F}_{t-1}$, i.e., 1-step ahead forecasts are
predictable w.r.t. $(\mathcal{F}_{t})_{t=1}^{\infty}$, consistent with our
previous formulation in Section 3.2.) It then follows that
$(p_{t}^{(k)})_{t=1}^{\infty}$ is predictable w.r.t.
$(\mathcal{F}_{t})_{t=1}^{\infty}$ and $(\mathcal{G}_{t})_{t=1}^{\infty}$ for
each $k$, as
$p_{t}^{(k)}\in\mathcal{F}_{t-k}\subseteq\mathcal{F}_{t-1}\subsetneq\mathcal{G}_{t-1}$.
This gives us the analog of (38):

$\operatorname*{\mathbb{E}}\left[{\tilde{S}_{\mathbf{w}}(\mathbf{p}_{t},y_{t})}\;\middle|\;{\mathcal{G}_{t-1}}\right]=\sum_{k=1}^{K}w_{k}\operatorname*{\mathbb{E}}\left[{\tilde{S}(p_{t}^{(k)},y_{t})}\;\middle|\;{\mathcal{G}_{t-1}}\right]=\sum_{k=1}^{K}w_{k}\tilde{S}(p_{t}^{(k)},r_{t})=\tilde{S}_{\mathbf{w}}(\mathbf{p}_{t},r_{t}).$
(47)

Also, analogous to (39), we have that

$\displaystyle
S_{\mathbf{w}}(\mathbf{p}_{t},r_{t})-S_{\mathbf{w}}(\mathbf{q}_{t},r_{t})$
$\displaystyle=\sum_{k=1}^{K}w_{k}\left[S(p_{t}^{(k)},r_{t})-\tilde{S}(q_{t}^{(k)},r_{t})\right]$
$\displaystyle=\sum_{k=1}^{K}w_{k}\left[\left(\tilde{S}(p_{t}^{(k)},r_{t})+d(r_{t})\right)-\left(\tilde{S}(q_{t}^{(k)},r_{t})-d(r_{t})\right)\right]$
$\displaystyle=\sum_{k=1}^{K}w_{k}\left[\tilde{S}(p_{t}^{(k)},r_{t})-\tilde{S}(q_{t}^{(k)},r_{t})\right]$
$\displaystyle=\tilde{S}_{\mathbf{w}}(\mathbf{p}_{t},r_{t})-\tilde{S}_{\mathbf{w}}(\mathbf{q}_{t},r_{t})$
(48)

where $d(r):=S(p,r)-\tilde{S}(p,r)$ is the function that does not depend on
$p$.

Given (47) and (B.4), the proof of Theorem 2 can readily be extended to derive
the analogous confidence sequences for average forecast score differentials
under $S_{\mathbf{w}}$. We state this extension as a corollary.

###  B.5 Proof of Lemma 4

First, observe that $(T(p_{t},q_{t}))_{t=1}^{\infty}$ is predictable w.r.t.
$(\mathcal{G}_{t})_{t=1}^{\infty}$, since
$T(p_{t},q_{t})\in\mathcal{F}_{t-1}\subset\mathcal{G}_{t-1}$. This means that
the equivalent of (38) holds for (32):

$\displaystyle\operatorname*{\mathbb{E}}\left[{w(p_{t},q_{t},y_{t})}\;\middle|\;{\mathcal{G}_{t-1}}\right]$
$\displaystyle=\frac{\operatorname*{\mathbb{E}}\left[{S(p_{t},y_{t})-S(q_{t},y_{t})}\;\middle|\;{\mathcal{G}_{t-1}}\right]}{T(p_{t},q_{t})}$
$\displaystyle=\frac{\operatorname*{\mathbb{E}}\left[{\tilde{S}(p_{t},y_{t})-\tilde{S}(q_{t},y_{t})}\;\middle|\;{\mathcal{G}_{t-1}}\right]}{T(p_{t},q_{t})}$
$\displaystyle=\frac{\tilde{S}(p_{t},r_{t})-\tilde{S}(q_{t},r_{t})}{T(p_{t},q_{t})}$
$\displaystyle=\frac{S(p_{t},r_{t})-S(c_{t},r_{t})}{T(p_{t},q_{t})}=w(p_{t},q_{t},r_{t}),$

where we use (39) in the second and fourth equalities and (38) in the third
equality. Note that the score differences are bounded by a value that depends
on both $q_{0}$ and $S$. For the Brier score, this bound is given by
$\left[-\frac{2}{q_{0}},\frac{2}{q_{0}}\right]$.

##  Appendix C Computing the Gamma-Exponential Mixture

Here, we derive the gamma-exponential mixture corresponding to the e-process
for the weak null, as discussed in Section 3.5. By Proposition 9 in Howard et
al., (2021), the expression for the gamma-exponential mixture is given by

$m(s,v):=\frac{\left(\frac{\rho}{c^{2}}\right)^{\frac{\rho}{c^{2}}}}{\Gamma\left(\frac{\rho}{c^{2}}\right)P\left(\frac{\rho}{c^{2}},\frac{\rho}{c^{2}}\right)}\frac{\Gamma\left(\frac{v+\rho}{c^{2}}\right)P\left(\frac{v+\rho}{c^{2}},\frac{cs+v+\rho}{c^{2}}\right)}{\left(\frac{cs+v+\rho}{c^{2}}\right)^{\frac{v+\rho}{c^{2}}}}\exp\left\\{\frac{cs+v}{c^{2}}\right\\},$
(49)

where $\Gamma(a,z):=\int_{z}^{\infty}u^{a-1}e^{-u}du$ is the upper incomplete
gamma function, $\Gamma(a):=\Gamma(a,0)$ is the gamma function, and $P$ is the
regularized lower incomplete gamma function defined as

$P(a,z):=\frac{1}{\Gamma(a)}\int_{0}^{z}u^{a-1}e^{-u}du$ (50)

for $a,z>0$. Note that both $\Gamma$ and $P$ can computed efficiently in
standard scientific computing software. (E.g., $P$ can be computed using
boost::math::gamma_p in C++ and scipy.special.gammainc in Python.)

In the following, we derive the expression (49) and show its upper bound in
the case where $s+v+\rho<0$. For notational simplicity, we assume $c=1$, but
an analogous proof can be made for any $c>0$.

Recall that $\psi_{E}(\lambda)=-\log(1-\lambda)-\lambda$ for
$\lambda\in[0,1)$. Using the gamma density
$f_{\rho}(\lambda)=C(\rho)(1-\lambda)^{\rho-1}e^{-\rho(1-\lambda)}$, where
$C(\rho)=\frac{\rho^{\rho}}{P(\rho,\rho)\Gamma(\rho)}$, we have that, for
$\rho>0$,

$\displaystyle m(s,v)$ $\displaystyle=C(\rho)\int_{0}^{1}\exp\left\\{\lambda
s-\psi_{E}(\lambda)v\right\\}\cdot(1-\lambda)^{\rho-1}e^{-\rho(1-\lambda)}d\lambda$
$\displaystyle=C(\rho)\int_{0}^{1}e^{\lambda(s+v)}(1-\lambda)^{v}\cdot(1-\lambda)^{\rho-1}e^{-\rho(1-\lambda)}d\lambda$
$\displaystyle=C(\rho)\int_{0}^{1}(1-\lambda)^{v+\rho-1}e^{\lambda(s+v)-\rho(1-\lambda)}d\lambda$
$\displaystyle=C(\rho)\left(\int_{0}^{1}(1-\lambda)^{v+\rho-1}e^{-(s+v+\rho)(1-\lambda)}d\lambda\right)e^{s+v},$
(51)

where in the last equality we used

$\lambda(s+v)-\rho(1-\lambda)=(s+v)-(1-\lambda)(s+v)-(1-\lambda)\rho=-(s+v+\rho)(1-\lambda)+(s+v).$

Now, let $a=v+\rho$ and $z=s+v+\rho$, and note that $a>0$.

If $z>0$, then using the change-of-variable formula
$u=(s+v+\rho)(1-\lambda)=z(1-\lambda)$, we have that

$\displaystyle m(s,v)$
$\displaystyle=C(\rho)\left(\int_{z}^{0}\left(\frac{u}{z}\right)^{a-1}e^{-u}\frac{du}{-z}\right)e^{s+v}$
$\displaystyle=C(\rho)\cdot\frac{1}{z^{a}}\left(\int_{0}^{z}u^{a-1}e^{-u}du\right)e^{s+v}$
(52) $\displaystyle=C(\rho)\frac{\Gamma(a)P(a,z)}{z^{a}}e^{s+v},$ (53)

where we use the fact that the integral in (52) corresponds to the numerator
of the lower incomplete gamma function $P(a,z)$ in (50). The expression (53)
can be computed in closed-form.

On the other hand, if $z<0$, then using the change-of-variable formula
$u=-(s+v+\rho)(1-\lambda)=-z(1-\lambda)$, we obtain

$\displaystyle m(s,v)$
$\displaystyle=C(\rho)\left(\int_{-z}^{0}\left(\frac{u}{-z}\right)^{a-1}e^{u}\frac{du}{z}\right)e^{s+v}$
$\displaystyle=C(\rho)\cdot\frac{1}{(-z)^{a}}\left(\int_{0}^{-z}u^{a-1}e^{u}du\right)e^{s+v}$
$\displaystyle=C(\rho)\cdot\frac{1}{|z|^{a}}\left(\int_{0}^{|z|}u^{a-1}e^{u}du\right)e^{s+v}.$
(54)

Although the integral in (54) is no longer a regularized lower gamma function,
we can still show that $m(s,v)$ as a whole is bounded by $e$. Since $e^{u}\leq
e^{|z|}=e^{-z}$ for $u\leq|z|$, we have that

$\displaystyle m(s,v)$ $\displaystyle\leq
C(\rho)\cdot\frac{1}{|z|^{a}}\left(\int_{0}^{|z|}u^{a-1}du\right)e^{-z}\cdot
e^{s+v}$
$\displaystyle=C(\rho)\cdot\frac{1}{|z|^{a}}\left(\int_{0}^{|z|}u^{a-1}du\right)e^{-\rho}$
(55)
$\displaystyle=C(\rho)\cdot\frac{1}{|z|^{a}}\left.\left(\frac{u^{a}}{a}\right)\right|_{0}^{|z|}e^{-\rho}$
$\displaystyle=\frac{C(\rho)e^{-\rho}}{v+\rho},$ (56)

where in (55) we used $-z+(s+v)=-(s+v+\rho)+(s+v)=-\rho$, and in (56) we
substituted in $a=v+\rho$. We can further bound this value, using the fact
that $v>0$ and substituting back in $C(\rho)$:

$\displaystyle m(s,v)\leq\frac{C(\rho)e^{-\rho}}{v+\rho}$
$\displaystyle\leq\frac{C(\rho)e^{-\rho}}{\rho}$
$\displaystyle=\rho^{\rho-1}e^{-\rho}\cdot\left(\int_{0}^{\rho}u^{\rho-1}e^{-u}du\right)^{-1}$
$\displaystyle\leq\rho^{\rho-1}e^{-\rho}\cdot\left(e^{-\rho}\int_{0}^{\rho}u^{\rho-1}du\right)^{-1}$
(57)
$\displaystyle=\rho^{\rho-1}\cdot\left[\left.\left(\frac{u^{\rho}}{\rho}\right)\right|_{0}^{\rho}\right]^{-1}$
$\displaystyle=1,$ (58)

where in (57) we used the fact that $e^{-\rho}\leq e^{-u}$ for $u\in[0,\rho]$.
It makes sense that $m(s,v)$ is upper-bounded by 1 when $z=s+v+\rho<0$,
because $s<-(v+\rho)<0$ would imply that the sum of score differentials is
negative, which is an evidence that supports the weak null. In our
implementation, we use the upper bound (56), which can be computed easily and
gets substantially smaller than 1 when $v+\rho\gg 0$.

##  Appendix D Comparing CS Widths on IID Means

Figure 10: (Left) Histogram of
$\delta_{i}\overset{\text{IID}}{\sim}\mathrm{Beta}(30,10)-\mathrm{Beta}(10,30)$
for $i=1,\dotsc,10,000$. (Right) Plot of the intrinsic time
$\hat{V}_{t}=\sum_{i=1}^{t}(\hat{\delta}_{i}-\hat{\Delta}_{i-1})^{2}$, where
$\hat{\Delta}_{i-1}=\sum_{j=1}^{i-1}\hat{\delta}_{j}$. Note that the
horizontal axis $t$ is drawn in log-scale. Also note that the hyperparameter
$v_{\text{opt}}$ determines the intrinsic time $\hat{V}_{t}$ at which the
uniform boundary is the tightest.

To validate and finetune our CSs, we compare the widths of various time-
uniform CS applicable to the mean of independent and identically distributed
(IID) random variables. We compare both the Hoeffding-style (H) CS (Theorem 1
and the empirical-Bernstein (EB) CS (Theorem 2) using both the conjugate-
mixture (CM) and the polynomial stitching uniform boundaries, as described in
Section 3.4. We also include predictable-mixed confidence sequences, namely
the predictably-mixed Hoeffding (PM-H) CS and empirical Bernstein (PM-EB) CS
from Waudby-Smith and Ramdas, (2020), which are time-uniform CS that are also
applicable for means of IID random variables. We also include the asymptotic
CS described in Waudby-Smith et al., (2021). As for the data, we use the
difference between two i.i.d. Beta random variables, as a proxy for score
differentials between two sets of forecasts: for $i=1,\dotsc,10,000$,

$\delta_{i}\overset{\text{IID}}{\sim}\mathrm{Beta}(30,10)-\mathrm{Beta}(10,30).$
(59)

Note that $-1\leq\delta_{i}\leq 1$ a.s. and that
$\operatorname*{\mathbb{E}}\left[{\delta_{i}}\right]=\frac{30}{30+10}-\frac{10}{10+30}=\frac{1}{2}$.
Figure 10 illustrates the data sampled according to (59) (left) as well as the
intrinsic time
$\hat{V}_{t}=\sum_{i=1}^{t}(\hat{\delta}_{i}-\hat{\Delta}_{i-1})^{2}$, where
$\hat{\Delta}_{i-1}=\sum_{j=1}^{i-1}\hat{\delta}_{j}$, over time (right). Note
that, when choosing the uniform boundary for the CSs, the hyperparameter
$v_{\text{opt}}$ determines the intrinsic time $\hat{V}_{t}$ at which the
width of the uniform boundary is the tightest (Proposition 3, Howard et al.,
(2021)).

Figure 11: (Left) Hyperparameter tuning for the width of CM-EB CS by adjusting
the optimal intrinsic time parameter $v_{\text{opt}}$. The choice
$v_{\text{opt}}=10$ gives the smallest width overall, although
$v_{\text{opt}}=100$ is comparable for $t>10^{3}$. The width of stitching EB
CS (with $v_{\text{opt}}=1$) is also drawn as a point of comparison. (Right)
Comparing the widths of different types of time-uniform CS. Overall, the
asymptotic CS is the tightest, followed by the EB variants and then by the
Hoeffding variants.

In Figure 11 (left), we present the hyperparameter tuning results of the CM-EB
CS with respect to its optimal intrinsic time parameter $v_{\text{opt}}$,
which corresponds to the point in the variance process where the boundary is
the tightest. As a point of comparison, we also plot the polynomial stitching
EB CS with hyperparameters $v_{\text{opt}}=10$, $s=1.4$, and $\eta=2$.
Comparing the values of $v_{\text{opt}}\in\\{0.1,1,10,100,1000\\}$, we find
that the CM-EB CS is the tightest across time (up to $T=10^{4}$) with
$v_{\text{opt}}=10$. Based on these results, we use the CM-EB CS with
$v_{\mathrm{opt}}=10$ for our other experiments throughout the paper (except
for the baseball experiments, which we conduct as a post-hoc analysis using
the sample variance from the data).

In Figure 11 (right), we plot the widths of the different kinds of time-
uniform CS, optimized for the intrinsic time $v_{\mathrm{opt}}=10$ when
applicable. Generally speaking, we observe that the CSs are the tightest for
the asymptotic CS, followed by the EB CS variants and the Hoeffding CS
variants. This is consistent with our intuition, as the EB CS additionally
makes use of the estimated variance to achieve smaller widths than the
Hoeffding CS, and the asymptotic CS is the “limit” of EB CS in terms of the
width while sacrificing the non-asymptotic guarantee (Waudby-Smith et al.,,
2021). For the EB CS variants, the CM-EB CS is tighter towards the beginning
($t<10^{3}$) while the PM-EB CS becomes slightly tighter afterwards, and the
stitching CS is not as tight as the other two. This is also as expected, as
the CM-EB and PM-EB CS have similar widths (up to differences determined by
the choice of hyperparameters) (Waudby-Smith and Ramdas,, 2020) and the
stitching CS tends to be looser in practice (Howard et al.,, 2021). This trend
is also analogous for the Hoeffding CS variants, although the stitching
variant does become tighter for larger $t$ in this case.


the $k$ paths in $\Pi$ arbitrarily. For each candidate $k$-shredder $S$ in
$\mathcal{C}$, define the $k$-tuple

$(\delta_{\pi_{1}}(s_{1}),\delta_{\pi_{2}}(s_{2}),\dots,\delta_{\pi_{k}}(s_{k}))$

where $s_{i}=\pi_{i}(S)$ is the vertex in $S$ that lies along the $i$-th path
$\pi_{i}$ in $\Pi$. This $k$-tuple can be represented as a $k$-digit number in
base $\ell$, inducing the lexicographical order between candidate
$k$-shredders in $\mathcal{C}$. We sort the candidates by their $k$-tuple and
organize the sorted list in a $|\mathcal{C}|\times k$ table (see Figure 7) by
non-decreasing order from top to bottom. This can be done via any stable
sorting algorithm to sort each $k$-tuple for each of the $k$ digits, starting
from the least significant digit to the most.

Figure 7: Candidate $k$-shredders sorted by their $k$-tuple numbers. Each row
represents the $k$-tuple number $T_{S}$ for some candidate $k$-shredder
$S\in\mathcal{C}$. Candidates are sorted in non-decreasing order from top to
bottom. Here we have fixed two rows for two candidate $k$-shredders $S,T$ that
straddle each other.

Using this sorted table, we can detect candidate $k$-shredders that are not
straddled by others. The main idea is as follows. Suppose that a candidate
$k$-shredder $S$ straddles another candidate $k$-shredder $T$. Let
$s_{i}=\pi_{i}(S)$ and $t_{i}=\pi_{i}(T)$. By Definition 4.7, there exist two
paths $\pi_{i},\pi_{j}\in\Pi$ (with $i<j$) such that
$\delta_{\pi_{i}}(s_{i})<\delta_{\pi_{i}}(t_{i})$ and
$\delta_{\pi_{j}}(s_{j})>\delta_{\pi_{j}}(t_{j})$. Without loss of generality,
let $i$ be the smallest path index such that the above holds. Observe that
because $i<j$ and $\delta_{\pi_{i}}(s_{i})<\delta_{\pi_{i}}(t_{i})$, the row
for $S$ must be sorted above the row for $T$. On the other hand, in the $j$-th
column of this table, we will see that the value in the row for $S$
($\delta_{\pi_{j}}(s_{j})$) is greater than the value in the row for $T$
($\delta_{\pi_{j}}(t_{j})$). This inversion is the essential marker that we
use to detect straddling. To formalize, a candidate $k$-shredder on row $s$ is
straddled by another candidate if and only if there exists a column of the
table with top-to-bottom-ordered entries $(e_{1},\ldots,e_{|\mathcal{C}|})$
that contains an entry $e_{t}$ with $t<s$ and $e_{t}>e_{s}$, or there exists
an entry $e_{t}$ with $t>s$ and $e_{t}<e_{s}$.

Detecting which rows are involved in an inversion can be done as follows. For
each column of the table, we scan through the column twice. In the first run,
we scan from top to bottom, keeping track of the maximum value seen so far. If
at any point, the entry on row $s$ is less than the maximum value seen so far,
then we can infer that the $s$-th candidate $k$-shredder in the table was
straddled. In the second run, we scan from bottom to top, keeping track of the
minimum value seen so far. If the next entry on row $s$ is greater than the
minimum, we can infer that the $s$-th candidate $k$-shredder in the table was
straddled.

We will prove why this scanning process suffices. Consider a candidate
$k$-shredder sorted on row $s$. Suppose there exists a column of the table
with top-to-bottom-ordered entries $(e_{1},\ldots,e_{|\mathcal{C}|})$ such
that there exists an entry $e_{t}$ with $t<s$ and $e_{t}>e_{s}$, or there
exists an entry $e_{t}$ with $t>s$ and $e_{t}<e_{s}$. If $t<s$ and
$e_{t}>e_{s}$, this implies that $e_{s}<\max\\{e_{1},\ldots,e_{s-1}\\}$. We
will detect this when scanning the column from top to bottom. Similarly, if
$t>s$ and $e_{t}<e_{s}$, then
$e_{s}>\min\\{e_{s+1},\ldots,e_{|\mathcal{C}|}\\}$, which is detected when
scanning the column from bottom to top. Conversely, suppose that when scanning
a column from top to bottom, we detect that row $s$ is less than the maximum
entry seen so far. Then, by definition, we have detected an inversion
involving row $s$. When scanning from bottom to top, if we detect that row $s$
is greater than the minimum entry seen so far, we again have detected an
inversion involving row $s$. Hence, this scanning process precisely captures
any inversions if they occur.

Now, we calculate the time complexity. Firstly, we preprocess $\Pi$ so that
for a vertex $v\in\Pi$, we can query the distance of $v$ from $x$ in constant
time. This can be done in $\mathcal{O}(\ell)$ time by traversing each path in
$\Pi$ starting from $x$ and marking the distance from $x$ for each traversed
vertex. Next, we build the $k$-tuple numbers for each candidate in
$\mathcal{C}$. This takes $\mathcal{O}(k|\mathcal{C}|)$ time. Since there are
$|\mathcal{C}|$ $k$-tuple numbers, each of which are $k$ digits in base
$\ell$, lexicographically sorting them by their digits takes
$\mathcal{O}(k|\mathcal{C}|\log|\mathcal{C}|)$ time. The column scanning
process can be completed in $\mathcal{O}(k|\mathcal{C}|)$ operations by
iterating through the table twice. Hence, the entire algorithm runs in
$\mathcal{O}(k|\mathcal{C}|\log|\mathcal{C}|+\ell)$ time. ∎

###### Proof of Lemma 4.9.

Firstly, we apply Lemma B.1 to eliminate candidate $k$-shredders that are
straddled by other candidates in
$\mathcal{O}(k|\mathcal{C}|\log|\mathcal{C}|+\ell)$ time. We are left with a
set of bridges $\mathcal{B}$ and a set $\mathcal{C}$ of candidate
$k$-shredders sorted by their $k$-tuple numbers in nondecreasing order.
Because we have eliminated pairwise candidate straddling, the sorted list of
remaining candidate $k$-shredders is well-ordered under $\preceq,\succeq$.
That is, for any pair of candidate $k$-shredders $S_{i},S_{j}$ ($i<j$), we
have $S_{i}\preceq S_{j}$. Hence, if $S_{j}\preceq\Gamma$ for a bridge
$\Gamma$ of $\Pi$, then $S_{i}\preceq\Gamma$ for all $i<j$. The goal is to
develop a procedure that takes a bridge $\Gamma$ and determines all candidate
$k$-shredders straddled by $\Gamma$ in $\mathcal{O}(\mathrm{vol}(\Gamma))$
time.

Let $\pi\in\Pi$. Given a vertex $z\in\pi$, we can query for the greatest index
$lo$ such that $S_{lo}$ satisfies
$\delta_{\pi}(\pi(S_{lo}))\leq\delta_{\pi}(z)$. If no such index exists, then
we set $lo=0$. Similarly, we can query the smallest index $hi$ such that
$S_{hi}$ satisfies $\delta_{\pi}(\pi(S_{hi}))\geq\delta_{\pi}(z)$. If no such
index exists, we set $hi=|\mathcal{C}|+1$. Using the sorted list of candidate
$k$-shredders, we can preprocess the vertices of $\Pi$ in
$\mathcal{O}(k|\mathcal{C}|+\ell)$ time to return $lo$ and $hi$ in constant
time. To extend this a bit further, for each bridge $\Gamma\in\mathcal{B}$ we
can compute the greatest index $lo$ and the least index $hi$ (if they exist)
such that $S_{lo}\preceq\Gamma$ and $S_{hi}\succeq\Gamma$ in
$\mathcal{O}(\mathrm{vol}(\Gamma))$ time. This can be done by simply scanning
through the attachment set of $\Gamma$ and performing the queries described
above. This is especially useful because every candidate $k$-shredder $S_{i}$
with $i\in(lo,hi)$ satisfies $S_{i}\not\preceq\Gamma$ and
$S_{i}\not\succeq\Gamma$, which means $\Gamma$ straddles $S_{i}$. With these
steps, we can compute an open interval $(lo,hi)$ containing the indices of all
candidate $k$-shredders which are straddled by $\Gamma$. To find all candidate
$k$-shredders which are straddled by any bridge in $\mathcal{B}$, we can
simply merge all the $(lo,hi)$ intervals for each bridge
$\Gamma\in\mathcal{B}$. Since $lo$ and $hi$ are nonnegative values that are on
the order of $|\mathcal{C}|$, we can merge all intervals in
$\mathcal{O}(|\mathcal{C}|)$ time. Tallying up all the steps, we have that the
entire algorithm runs in
$\mathcal{O}(k|\mathcal{C}|\log|\mathcal{C}|+\ell)+\mathcal{O}(\mathrm{vol}(\mathcal{B}))+\mathcal{O}(|\mathcal{C}|)=\mathcal{O}(k|\mathcal{C}|\log|\mathcal{C}|+\ell+\mathrm{vol}(\mathcal{B}))$
time. ∎

See 4.1

###### Proof of 4.1.

Given a pair of vertices $(x,y)$, obtaining a set $\Pi$ of $k$ openly-disjoint
simple paths from $x$ to $y$ can be done in $\mathcal{O}(km)$ time. Finding
all the bridges of $\Pi$ and obtaining the list of candidate $k$-shredders can
be done in $\mathcal{O}(m+n)$ time using DFS/BFS. Note that the number of
candidate $k$-shredders is $\mathcal{O}(n)$ because there are at most
$\mathcal{O}(n)$ components of $G\setminus\Pi$.

In [CT99], Cheriyan and Thurimella actually used radix sort to obtain the
sorted table described in figure Figure 7. Specifically, they noted that each
digit in a $k$-tuple is bounded by $n$ because all simple paths from $x$ to
$y$ are of length at most $n$. Hence, sorting at most $n$ $k$-digit numbers in
base $n$ can be done in $\mathcal{O}(k(n+n))=\mathcal{O}(kn)$ time using radix
sort. After obtaining the sorted table, they followed the remaining steps in
Lemma 4.9. This entire pruning process therefore takes
$\mathcal{O}(k(n+n)+m)=\mathcal{O}(kn+m)$ time. Hence, by Lemma 4.8, we have
obtained all $k$-shredders separating $x$ and $y$. In total, the algorithm
takes $\mathcal{O}(km)+\mathcal{O}(m+n)+\mathcal{O}(kn+m)=\mathcal{O}(km)$
time. ∎

##  Appendix C Omitted Proofs from Section 5

See 5.11

###### Proof.

We can assume that $|U|=k$ as otherwise $U=\\{x\\}$, which is caught by line
13. Firstly, we have $\delta_{\pi}(\pi(U))<\delta_{\pi}(z)$ because $U$ does
not contain the endpoints of any path by line 17. The proof now becomes
identical to that of Lemma 5.9.

Assume for the sake of contradiction that $x$ and $z$ are connected in
$G\setminus U$. Then, there exists a path $P$ from $x$ to $z$ in $G\setminus
U$. Let $u$ denote the last vertex along $P$ such that $u\in\pi_{1}$ for some
path $\pi_{1}\in\Pi$ and $\delta_{\pi_{1}}(u)<\delta_{\pi_{1}}(\pi_{1}(U))$.
Such a vertex must exist because $x\in P$ and satisfies
$\delta_{\pi}(x)<\delta_{\pi}(\pi(U))$. Similarly, let $v$ denote the first
vertex along $P$ after $u$ such that $v\in\pi_{2}$ for some path
$\pi_{2}\in\Pi$ and $\delta_{\pi_{2}}(v)>\delta_{\pi_{2}}(\pi_{2}(U))$. Such a
vertex must exist because $z\in P$ and satisfies
$\delta_{\pi}(\pi(U))<\delta_{\pi}(z)$.

Consider the open interval $P(u{\rightsquigarrow}v)$ of $P$. We claim that
every inner vertex $w$ in the open interval $P(u{\rightsquigarrow}v)$ must not
be in $\Pi$. To see this, suppose that there exists a vertex $w\in
P(u{\rightsquigarrow}v)$ such that $w\in\pi_{3}$ for a path $\pi_{3}\in\Pi$.
Then, $\delta_{\pi_{3}}(w)=\delta_{\pi_{3}}(\pi_{3}(U))$ as otherwise, we
contradict the definition of $u$ or $v$. This implies $w=\pi_{3}(U)$ because
the vertex distance $\delta_{\pi_{3}}(\pi_{3}(U))$ from $x$ along $\pi_{3}$ is
unique. But this contradicts the definition of $P$ as $P$ is a path in
$G\setminus U$.

Hence, the vertices of $P(u{\rightsquigarrow}v)$ exist as a bridge $\Gamma$ of
$\Pi$ with attachment set $\\{u,v,\ldots\\}$. If there are no vertices in
$P(u{\rightsquigarrow}v)$, then $\Gamma$ is the edge $(u,v)$. Furthermore, we
have that $\delta_{\pi_{1}}(u)<\delta_{\pi_{1}}(\pi_{1}(U))$ and
$\delta_{\pi_{2}}(v)>\delta_{\pi_{2}}(\pi_{2}(U))$. Therefore, $\Gamma$
straddles $U$. Lemma 5.10 implies $U$ would have been pruned. This is checked
by line 17, which implies the algorithm would have returned $\varnothing$ for
$U$, giving our desired contradiction. ∎

See 5.12

###### Proof.

This was omitted from the pseudocode for simplicity. However, we can easily
compute $\mathrm{vol}(Q_{x})$ at the end of the algorithm by iterating through
each path and exploring the bridges of $\Pi$ connected to the vertices
preceding the unverified vertex. This takes $\mathcal{O}(\nu)$ time for each
path because the unverified vertex $u$ is defined such the bridges attached to
vertices preceding $u$ have total volume at most $\nu$. Hence, over all paths,
this can be done in $\mathcal{O}(k\nu)$ time. ∎

##  Appendix D Proof of Lemma 6.7

See 6.7

###### Proof.

If we sample an edge $(u,v)$ uniformly at random, we have

$\displaystyle\Pr[u\in Q]$
$\displaystyle\geq\frac{1}{2}\cdot\frac{\mathrm{vol}(Q)}{m}$
$\displaystyle\rightarrow$ $\displaystyle\Pr[u\notin Q]$ $\displaystyle\leq
1-\frac{\mathrm{vol}(Q)}{2m}.$

If we independently sample $T$ edges $(u,v)$, we have

$\bigcap_{(u,v)}\Pr[u\notin
Q]\leq\left(1-\frac{\mathrm{vol}(Q)}{2m}\right)^{T}.$

If we set $T\geq\frac{2m}{\mathrm{vol}(Q)}$, then

$\left(1-\frac{\mathrm{vol}(Q)}{2m}\right)^{T}\leq\left(1-\frac{\mathrm{vol}(Q)}{2m}\right)^{\frac{2m}{\mathrm{vol}(Q)}}\leq\frac{1}{e}.$

Notice that the failure probability is now a constant which we can boost by
scaling $T$ by a polylogarithmic factor. Because
$\frac{1}{2}\nu<\mathrm{vol}(Q)$, we have

$\frac{2m}{\mathrm{vol}(Q)}\leq\frac{2m}{\frac{1}{2}\nu}\leq\frac{2m}{\frac{1}{2}\nu}\cdot
100\log n=400\frac{m}{\nu}\log n.$

Setting $T=400\frac{m}{\nu}\log n$, we get

$\displaystyle\bigcap_{(u,v)}\Pr[u\notin
Q]\leq\left(1-\frac{\mathrm{vol}(Q)}{2m}\right)^{T}=\left(1-\frac{\mathrm{vol}(Q)}{2m}\right)^{\frac{2m}{\frac{1}{2}\nu}\cdot
100\log n}$ $\displaystyle\leq\left(\frac{1}{e}\right)^{100\log n}$
$\displaystyle\leq\left(\frac{1}{2}\right)^{100\log n}$ $\displaystyle\leq
n^{-100}.$

Therefore, with probability $1-n^{-100}$, we will sample an edge $(u,v)$ such
that $u\in Q$. ∎