# Global Tests for Smoothed Functions in Mean Field Variational Additive
Models

Mark J. Meyer and Junyi Wei
Department of Mathematics and Statistics, Georgetown University
Washington, DC USA

###### Abstract

Variational regression methods are an increasingly popular tool for their
efficient estimation of complex. Given the mixed model representation of
penalized effects, additive regression models with smoothed effects and
scalar-on-function regression models can be fit relatively efficiently in a
variational framework. However, inferential procedures for smoothed and
functional effects in such a context is limited. We demonstrate that by using
the Mean Field Variational Bayesian (MFVB) approximation to the additive model
and the subsequent Coordinate Ascent Variational Inference (CAVI) algorithm,
we can obtain a form of the estimated effects required of a Frequentist test
for semiparametric curves. We establish MFVB approximations and CAVI
algorithms for both Gaussian and binary additive models with an arbitrary
number of smoothed and functional effects. We then derive a global testing
framework for smoothed and functional effects. Our empirical study
demonstrates that the test maintains good Frequentist properties in the
variational framework and can be used to directly test results from a
converged, MFVB approximation and CAVI algorithm. We illustrate the
applicability of this approach in a wide range of data illustrations.

##  1 Introduction

Variational approximations cover a useful class of algorithms for calculating,
in a deterministic fashion, estimates from Bayesian posterior distributions.
One such class, the Mean Field Variational Bayesian (MFVB) approximation, has
been well studied with many standard statistical models derived for use in the
class, see for example Ormerod and Wand (2010) or Blei et al. (2017). Recent
research suggests that estimates arising from MFVB share similar properties
with their Frequentist counter-parts including the consistency and the
asymptotic normality of variational estimators (Bickel et al., 2013; You et
al., 2014; Wang and Blei, 2019; Westling and McCormick, 2019; Yang et al.,
2020; Zhang and Zhou, 2020). Such evaluations demonstrate that variational
models have good Frequentist properties despite originating from a Bayesian
estimation paradigm.

One understudied area, however, is in hypothesis testing. While percentile-
based intervals can be constructed for parameters based on the conditional
posterior, underestimation of the variance can result in misleading
inference—particularly for non-normally distributed posteriors (Gelman et al.,
2013; Giordano et al., 2018). Some authors have proposed a bootstrap-based
approach, see Chen et al. (2018), while others convert back to a fully
Bayesian model when inference is of interest, for example Goldsmith and Kitago
(2016). For complicated variational regression models like additive models,
performing both the modeling and inference in the same framework would be
preferable to a more computationally burdensome approach of bootstrapping or
running a separate, fully Bayesian estimation procedure. One such model is the
additive model which typically implies a regression model with a scalar
outcome and an additive smoothed or semiparametric effect or potentially a set
of such effects in addition to scalar covariates

Additive models with smoothed effects have been broadly studied, see texts by
Ruppert et al. (2003) or Wood (2017), for example. Variational Bayesian
approaches that incorporate smoothing include work by Luts and Wand (2015);
Lee and Wand (2016); Wand (2017); Hui et al. (2019); Yang and Yang (2024) and
references therein. Methods that incorporate both smoothed and functional
effects into non-variational additive models include the work of McLean et al.
(2014); Scheipl et al. (2015, 2016); Li et al. (2016); Ma et al. (2021), among
others. However, to our knowledge, there is no work examining variational
additive models that accommodate an arbitrary number of both smoothed and
functional effects. And while there is work on global testing for smoothed
effects, see for example Zhang et al. (2000); Zhang and Lin (2003);
Crainiceanu et al. (2005); Greven et al. (2008); Wood (2017), such tests have
not been considered in the variational setting.

Global tests of the functional components in regression models are also
limited. Much of the work has emphasized point-wise testing, see for example
Malloy et al. (2010); Goldsmith et al. (2013) and Meyer et al. (2015). Some
authors have proposed global tests including Meyer et al. (2015) who examine
an approach in a fully Bayesian function-on-function setting based on what the
authors refer to as simultaneous band scores. Another example is Ivanescu et
al. (2015) who note the connection between penalized functional effects and
smoothed effects and suggest performing inference using an approach similar to
Crainiceanu and Ruppert (2004) or Greven et al. (2008). Neither authors
explore the operating characteristics of their global tests for functional
effects and such testing in the variational context has not been extensively
studied.

In this manuscript, we explore MFVB approximations to variational additive
models for both Gaussian and binary outcomes, developing Coordinate Ascent
Variational Inference (CAVI) algorithms to perform estimation (Section 2).
These models can incorporate an arbitrary number of both smoothed and
functional effects and thus include semiparametric regression models (in the
Ruppert et al. (2003) sense), scalar-on-function regression models (of the
form described by Ramsay and Dalzell (1991) and Ramsay and Silverman (2005)),
and the combination of the two. The target of inference we consider is either
a smoothed effect or a functional effect, both of which can be represented as
smoothed effects. Zhang et al. (2000) and Zhang and Lin (2003) propose a
Frequentist-based test for semiparametric curves which we refer to as the ZLS
test. We demonstrate that the CAVI algorithms for fitting MFVB approximations
to both Gaussian and binary variational additive models admit forms required
for the ZLS test which, we argue, can be used for a global test (Section 3).
We explore the proprieties of the ZLS test, including the type I error rate
and power, for global inference on both smoothed and functional effects
arising from converged CAVI algorithms (Section 4). In Section 5, we
demonstrate the use of the ZLS test in several data examples. Finally, we
provide a discussion of the method in Section 6.

##  2 Statistical Framework

###  2.1 Mean Field Variational Bayesian Approximations

We begin with a brief overview of the MFVB approximation. Let $\theta$
represent a vector of parameters and Y represent a vector of observed data.
The posterior distribution of $\theta$ given the data, Y, is then
$p({\mbox{\boldmath$\theta$}}|\textbf{Y})=p(\textbf{Y},{\mbox{\boldmath$\theta$}})\big{/}p(\textbf{Y})$.
The marginal distribution in the denominator, $p(\textbf{Y})$, is typically
intractable. Consequently, the posterior estimates must be obtained
algorithmically. In a fully Bayesian approach, this would involve a Markov
Chain Monte Carlo simulation or some other iterative random sampling technique
like Hamiltonian Monte Carlo, for example. The MFVB approach utilizes a
density transformation approach that relies on the observation that, for an
arbitrary density $q$, $p(\textbf{Y})$ is bounded below by
$\underline{p}(\textbf{Y};q)$ where
$\underline{p}(\textbf{Y};q)=\exp\left[\int
q({\mbox{\boldmath$\theta$}})\log\left\\{\frac{p(\textbf{Y},{\mbox{\boldmath$\theta$}})}{q({\mbox{\boldmath$\theta$}})}\right\\}d{\mbox{\boldmath$\theta$}}\right].$
The algorithms used to determine the approximations work by maximizing
$\underline{p}(\textbf{Y};q)$ over a class of tractable densities, $q$. As
$\underline{p}(\textbf{Y};q)$ is maximized, the Evidence Lower Bound (ELBO)
will be minimized. The ELBO is a metric measuring the distance between the
approximation, $q({\mbox{\boldmath$\theta$}})$, and the posterior,
$p({\mbox{\boldmath$\theta$}}|\textbf{Y})$, and is an equivalent, up to an
additive constant, to the Kullback-Leibler or K-L Divergence.

For a partition of the parameter space $\theta$ of size $M$,
$\\{{\mbox{\boldmath$\theta$}}_{1},\ldots,{\mbox{\boldmath$\theta$}}_{M}\\}$,
the MFVB approximation constructs $q({\mbox{\boldmath$\theta$}})$ using a
product density transformation of the form
$q({\mbox{\boldmath$\theta$}})=\prod_{m=1}^{M}q_{m}({\mbox{\boldmath$\theta$}}_{m})$.
Thus the resulting approximation is the product of the $q$-densities which are
similar to the conditional posterior densities in a Gibbs sampling framework
(Gelman et al., 2013). Optimization occurs via an iterative process with
convergence defined to occur when changes in the variational lower bound,
$\underline{p}(\textbf{Y};q)$, become minimal. Minimizing the changes in
$\underline{p}(\textbf{Y};q)$ in turn minimizes the ELBO, and therefore the KL
divergence, between $q({\mbox{\boldmath$\theta$}})$ and
$p({\mbox{\boldmath$\theta$}}|\textbf{Y})$. To perform this optimization, the
coordinate ascent variational inference (CAVI) algorithm is commonly used.
Ormerod and Wand (2010) and Blei et al. (2017) provide detailed reviews and
introductions to variational Bayesian techniques. Chapter 13 of Gelman et al.
(2013) also gives a good introduction to the topic. While there are other
variational approaches to approximating the posterior, our focus is on MFVB
approximations obtained via CAVI. The relationship to the Gibbs sampling
framework, where one component of $\theta$ is updated at a time, means the
algorithm will ultimately admit a useful form for the global test.

###  2.2 Variational Additive Models

####  2.2.1 Gaussian Outcome

Let $y_{i}$ be the continuous outcome of interest for subjects $i=1,\ldots,n$.
Each subject $i$ has a scalar covariate response vector, $\textbf{x}_{i}$,
that is $p\times 1$. The intercept, $\alpha$, and vector of non-smoothed, non-
functional coefficients is $\beta$ which is also $p\times 1$. Finally, we
denote the model errors with $\epsilon_{i}$ and assume
$\epsilon_{i}\stackrel{{\scriptstyle iid}}{{\sim}}N(0,\sigma_{e}^{2})$. We
place similar priors on $\alpha$, $\beta$, and $\sigma_{e}^{2}$ regardless of
the other model components. Thus, $\alpha\sim N(0,\sigma_{a}^{2})$,
${\mbox{\boldmath$\beta$}}\sim N(0,\sigma_{b}^{2}I_{p\times p})$, and
$\sigma_{e}^{2}\sim IG(a_{e},b_{e})$ where $\sigma_{a}^{2}$ and
$\sigma_{b}^{2}$ are fixed and set to something large, $I_{p\times p}$ denotes
a $p\times p$ identity matrix, and $IG$ denotes the inverse gamma
distribution. We take $a_{e}$ and $b_{e}$ to be fixed as well and set them to
something small, 0.01 for example.

For $M$ total smoothed effects and $F$ total functional effects, a general
additive model might have the form

$\displaystyle y_{i}=\alpha$
$\displaystyle+{\mbox{\boldmath$\beta$}}^{\prime}\textbf{x}_{i}+\sum_{m=1}^{M}s_{m}(z_{im})$
$\displaystyle+\sum_{f=1}^{F}\int_{t\in\mathcal{T}}w_{if}(t)\gamma_{f}(t)dt+\epsilon_{i},$
(1)

where $s_{m}(z_{im})$ is a smoothed effect and $\gamma_{f}(t)$ is a functional
effect. The functional effects, as written, are over the same time domain
although this does not have to be the case, i.e. $t$ and $\mathcal{T}$ can
vary with $f$. The covariates $z_{im}f$ and $w_{1f}(t)$ are not part of
$\textbf{x}_{i}$ which is a vector of scalar covariates. The vector $\beta$
contains the scalar effects and $\alpha$ is the model intercept.

Using basis expansions to represent each smoothed and functional effect, the
matrix version of Equation (1) is

$\displaystyle\textbf{Y}=\alpha\textbf{1}$
$\displaystyle+X{\mbox{\boldmath$\beta$}}+\sum_{m=1}^{M}\Xi_{m}{\mbox{\boldmath$\zeta$}}_{m}$
$\displaystyle+\sum_{f=1}^{F}W_{f}\Theta_{f}{\mbox{\boldmath$\lambda$}}_{f}+{\mbox{\boldmath$\epsilon$}},$
(2)

where Y and $\epsilon$ are both $n\times 1$ vectors and $X$ is an
$n\times(P-1)$ matrix of scalar effects. Considering the smoothed model
components first: for $K$ knots, $\Xi_{m}$ is the $n\times K$ matrix
containing the basis-expanded representation of the $z_{im}$, $i=1,\ldots,n$,
and ${\mbox{\boldmath$\zeta$}}_{m}$ is the $K\times 1$ vector of basis
coefficients. For the functional components, let $\Theta_{f}$ be a $T\times L$
matrix of basis functions and let $\lambda$ be a $L\times 1$ vector of basis
coefficients. Then $W_{f}$ is the $n\times T$ matrix containing the functional
covariate. Both smoothed and functional effects use a B-spline basis
expansion. Under empirical testing, we found that $K=8$ knots works well for
smoothed effects and $L=12$ works best for functional effects. In Equation
(2), we use the basis expansion ${\mbox{\boldmath$\gamma$}}=\Theta\lambda$
where $\gamma$ is the $T\times 1$ vector of functional coefficients,
$\gamma(t)$. Goldsmith et al. (2011) implement a similar approach for
variational scalar-on-function regression.

All basis expanded effects, both smoothed and functional, require penalization
to avoid overfitting. It is well known that penalized regression models can be
represented as mixed models (Ruppert et al., 2003, Chapter 4). In the
Frequentist context, this means the components one wishes to penalize are
treated as random effects. However, there is no distinction between random and
fixed effects in the Bayesian context as any unknown parameter is treated as
if it is random and has a prior placed on it. Thus, the equivalent is a prior
specification that introduces the penalty. Because of the relationship between
mixed effects models and penalized smoothing, we can leverage existing
algorithms for mean field variational mixed effect models to estimate Equation
(2). Ormerod and Wand (2010) provide one such algorithm which forms the basis
of our estimation procedure for the variational additive model.

We place a mean-zero shrinkage prior on the elements of
${\mbox{\boldmath$\zeta$}}_{m}$ of the form ${\mbox{\boldmath$\zeta$}}_{m}\sim
N(\textbf{0},\omega_{m}\mathcal{P}_{m})$ where $\omega_{m}$ is a tuning
parameter and 0 is a $K\times 1$ vector of zeros. The matrix $\mathcal{P}_{m}$
is a $K\times K$ penalty matrix which is the difference operator in matrix
form, see Eilers and Marx (1996) for more details. The prior on the basis
expanded functional effects is similar: ${\mbox{\boldmath$\lambda$}}_{f}\sim
N(\textbf{0},\eta_{f}\Delta_{f})$ where $\eta_{f}$ is a tuning parameter and 0
is an $L\times 1$ vector of zeros. We distinguish the penalty matrix for the
functional effect, $\Delta_{f}$, because it has a different form than in the
smoothed case: $\Delta_{f}=\xi\Delta_{0}+(1-\xi)\Delta_{2}$ where $\Delta_{0}$
is the zeroth derivative matrix and $\Delta_{2}$ is the second derivative
matrix of the B-spline basis function. This penalty matrix is consistent with
other penalized functional regression settings including Goldsmith and Kitago
(2016) and Meyer et al. (2022).

Regardless of effect type, the tuning parameters need to be estimated to
induce shrinkage. Thus, we place conditionally conjugate inverse-gamma priors
on all $\omega_{m}$ and $\eta_{f}$: $\omega_{m}\sim
IG(a_{\omega,m},b_{\omega,m})$ and $\eta_{m}\sim IG(a_{\eta,m},b_{\eta,m})$.
We set the hyper-parameters $a_{\omega,m},b_{\omega,m},a_{\eta,m},$ and
$b_{\eta,m}$, to something small, 0.01 or less. Since these parameters are
tuning parameters, it is reasonable to place a weakly informative prior on
each. However, our formulation does allow for differing amounts of information
to passed to different components when the model contains multiple smoothed
and functional effects. The Directed Acyclic Graph (DAG) describing the Markov
Blanket for this model is in Figure 1.

Y$\alpha$$\beta$$\zeta$$\lambda$$\omega$$\eta$ Figure 1: DAG for Gaussian
variational additive model

Combining the above model specifications, the resulting $q$ densities for MFVB
approximation are

$\displaystyle q({\mbox{\boldmath$\theta$}})$ $\displaystyle\sim
N\left[{\mbox{\boldmath$\mu$}}_{q({\mbox{\boldmath$\theta$}})},{\mbox{\boldmath$\Sigma$}}_{q({\mbox{\boldmath$\theta$}})}\right]$
$\displaystyle q(\sigma_{e}^{2})$ $\displaystyle\sim
IG\left[a_{e}+\frac{N}{2},B_{q(\sigma_{e}^{2})},\right]$ $\displaystyle
q(\omega_{1})$ $\displaystyle\sim
IG\left[a_{\omega,1}+\frac{K}{2},B_{q(\omega_{1})}\right],$
$\displaystyle\cdots$ $\displaystyle q(\omega_{M})$ $\displaystyle\sim
IG\left[a_{\omega,M}+\frac{K}{2},B_{q(\omega_{M})}\right],$ $\displaystyle
q(\eta_{1})$ $\displaystyle\sim
IG\left[a_{\eta,1}+\frac{L}{2},B_{q(\eta_{1})}\right],$ $\displaystyle\cdots$
$\displaystyle q(\eta_{F})$ $\displaystyle\sim
IG\left[a_{\eta,F}+\frac{L}{2},B_{q(\eta_{F})}\right],$

depending on the number of smoothed ($M$) and functional ($F$) components in
the model. The subscript-${q(\cdot)}$ notation indicates the parameter to
which the quantity belongs under the mean field approximation. The vector
$\theta$ contains all mean model parameters, i.e.
${\mbox{\boldmath$\theta$}}=\left[\begin{array}[]{cccccccc}\alpha&{\mbox{\boldmath$\beta$}}&{\mbox{\boldmath$\zeta$}}_{1}&\cdots&{\mbox{\boldmath$\zeta$}}_{M}&{\mbox{\boldmath$\lambda$}}_{1}&\cdots&{\mbox{\boldmath$\lambda$}}_{F}\end{array}\right]^{\prime}$.
Let C denote the design matrix for Equation (2), that is
$\textbf{C}=\left[\begin{array}[]{cccccccc}\textbf{1}&X&\Xi_{1}&\cdots&\Xi_{M}&W_{1}&\cdots&W_{F}\end{array}\right]$.
Algorithm 1 describes the estimation procedure. Convergence is determined when
changes in $\log[\underline{p}(\textbf{Y};q)]$ become negligible, where

$\displaystyle\log[$
$\displaystyle\underline{p}(\textbf{Y};q)]=\frac{1}{2}(P+K+L)-\frac{N}{2}\log(2\pi)$
$\displaystyle-\frac{P}{2}\log(\sigma_{b}^{2})+\frac{1}{2}\log\left(|{\mbox{\boldmath$\Sigma$}}_{q({\mbox{\boldmath$\theta$}})}|\right)$
$\displaystyle-\frac{1}{2\sigma_{b}^{2}}\left[{{\mbox{\boldmath$\mu$}}_{q(\alpha,{\mbox{\boldmath$\beta$}})}}^{\prime}{\mbox{\boldmath$\mu$}}_{q(\alpha,{\mbox{\boldmath$\beta$}})}+\text{tr}\left\\{{\mbox{\boldmath$\Sigma$}}_{q(\alpha,{\mbox{\boldmath$\beta$}})}\right\\}\right]$
$\displaystyle-
a_{e}\log(b_{e})-\left(a_{e}+\frac{N}{2}\right)\log(B_{q(\sigma^{2})})$
$\displaystyle+\log\left\\{\Gamma\left(a_{e}+\frac{N}{2}\right)\right\\}-\log\left\\{\Gamma(a_{e})\right\\}$
$\displaystyle+\sum_{m=1}^{M}\bigg{[}a_{\omega,m}\log(b_{\omega,m})$
$\displaystyle\hskip
30.0pt-\left(a_{\omega,m}+\frac{K}{2}\right)\log(B_{q(\omega_{m})})$
$\displaystyle\hskip
30.0pt+\log\left\\{\Gamma\left(a_{\omega,m}+\frac{K}{2}\right)\right\\}$
$\displaystyle\hskip 30.0pt-\log\\{\Gamma(a_{\omega,m})\\}\bigg{]}$
$\displaystyle+\sum_{f=1}^{F}\bigg{[}a_{\eta,f}\log(b_{\eta,f})$
$\displaystyle\hskip
30.0pt-\left(a_{\eta,f}+\frac{L}{2}\right)\log(B_{q(\eta_{f})})$
$\displaystyle\hskip
30.0pt+\log\left\\{\Gamma\left(a_{\eta,f}+\frac{L}{2}\right)\right\\}$
$\displaystyle\hskip 30.0pt-\log\\{\Gamma(a_{\eta,f})\\}\bigg{]}.$

The algorithm is written generally and can be used to estimate parameters from
Equation (1) for an arbitrary number of smoothed and functional effects.

Algorithm 1 Estimate Equation (1) components. The abbreviations ‘bd’ is for
block diagonal.

1:$B_{q(\omega_{1})},\ldots,B_{q(\omega_{M})},B_{q(\eta_{1})},\ldots,B_{q(\eta_{F})}>0$,
and $\epsilon>0$, small

2:while $\Delta\log[\underline{p}(\textbf{Y};q)]>\epsilon$ do

3: $\mathcal{D}\leftarrow\text{bd}\left\\{(\sigma_{b}^{2})^{-1}I_{p\times
p},\frac{a_{\omega,1}+\frac{K}{2}}{B_{q(\omega_{1})}}\mathcal{P}_{1},\cdots,\frac{a_{\omega,M}+\frac{K}{2}}{B_{q(\omega_{M})}}\mathcal{P}_{M},\frac{a_{\eta,1}+\frac{L}{2}}{B_{q(\eta_{1})}}\Delta_{1},\cdots,\frac{a_{\eta,F}+\frac{L}{2}}{B_{q(\eta_{F})}}\Delta_{M}\right\\}$

4:
${\mbox{\boldmath$\Sigma$}}_{q({\mbox{\boldmath$\theta$}})}\leftarrow\left[\frac{a_{e}+\frac{N}{2}}{B_{q(\sigma_{e}^{2})}}\textbf{C}^{\prime}\textbf{C}+\mathcal{D}\right]^{-1}$

5:
${\mbox{\boldmath$\mu$}}_{q({\mbox{\boldmath$\theta$}})}\leftarrow\left(\frac{a_{e}+\frac{N}{2}}{B_{q(\sigma_{e}^{2})}}\right){\mbox{\boldmath$\Sigma$}}_{q({\mbox{\boldmath$\theta$}})}\textbf{C}^{\prime}\textbf{Y}$

6: $B_{q(\sigma^{2})}\leftarrow
b_{e}+\frac{1}{2}\left[\left\\{\textbf{Y}-\textbf{C}{\mbox{\boldmath$\mu$}}_{q({\mbox{\boldmath$\theta$}})}\right\\}^{\prime}\left\\{\textbf{Y}-\textbf{C}{\mbox{\boldmath$\mu$}}_{q({\mbox{\boldmath$\theta$}})}\right\\}+\text{tr}\left\\{\textbf{C}^{\prime}\textbf{C}\right\\}{\mbox{\boldmath$\Sigma$}}_{q({\mbox{\boldmath$\theta$}})}\right]$

7: if $M\geq 1$ then

8: loop$\ m=1,\ldots,M$

9: $B_{q(\omega_{m})}\leftarrow
b_{\omega,m}+\frac{1}{2}\left[{{\mbox{\boldmath$\mu$}}_{q({\mbox{\boldmath$\zeta$}}_{m})}}^{\prime}{\mbox{\boldmath$\mu$}}_{q({\mbox{\boldmath$\zeta$}}_{m})}+\text{tr}\left\\{\frac{a_{\omega,m}+\frac{K}{2}}{B_{q(\omega_{m})}}\mathcal{P}_{m}\right\\}\right]$

10: end loop

11: end if

12: if $F\geq 1$ then

13: loop$\ f=1,\ldots,F$

14: $B_{q(\eta_{f})}\leftarrow
b_{\eta,f}+\frac{1}{2}\left[{{\mbox{\boldmath$\mu$}}_{q({\mbox{\boldmath$\lambda$}}_{f})}}^{\prime}{\mbox{\boldmath$\mu$}}_{q({\mbox{\boldmath$\lambda$}}_{f})}+\text{tr}\left\\{\frac{a_{\eta,f}+\frac{L}{2}}{B_{q(\eta_{f})}}\Delta_{f}\right\\}\right]$

15: end loop

16: end if

17:end while

####  2.2.2 Binary Outcome

Suppose instead that $y_{i}$ is a binary outcome for subjects $i=1,\ldots,n$.
As in the Gaussian case, $\textbf{x}_{i}$ is a vector of scalar covariates
that is $p\times 1$, $\alpha$ is the intercept, and $\beta$ is $p\times 1$
vector where both $\alpha$ and $\beta$ are non-smoothed, non-functional
effects. Prior specification for these components are the same as in Section
2.2.1. To model the binary outcome, we use the latent variable representation
which has the form

$\displaystyle y_{i}=\left\\{\begin{array}[]{cc}0&y_{i}^{*}<0\\\
1&y_{i}^{*}\geq 0,\end{array}\right.$ (5)

for the latent variable $y_{i}^{*}\in\mathbb{R}$. That is, the binary outcome
$y_{i}$ only represents the observable part of a continuous, underlying
process. We use the model from Equation (1) in the latent space:

$\displaystyle y_{i}^{*}=\alpha$
$\displaystyle+{\mbox{\boldmath$\beta$}}^{\prime}\textbf{x}_{i}+\sum_{m=1}^{M}s_{m}(z_{im})$
$\displaystyle+\sum_{f=1}^{F}\int_{t\in\mathcal{T}}w_{if}(t)\gamma_{f}(t)dt+\epsilon_{i}^{*}.$
(6)

Assuming that the latent errors, $\epsilon_{i}^{*}$, are Gaussian induces the
probit link, i.e. $\epsilon_{i}^{*}\stackrel{{\scriptstyle
iid}}{{\sim}}N(0,1)$. Since the representation in Equation (5) is scale-
invariant, the variance of the latent variable is taken to be 1, see Albert
and Chib (1993) for additional details on the Bayesian probit model.

Because we model in the latent space, the formulation follows analogously to
that in Section 2.2.1. Thus, the vectorized version of the latent model is

$\displaystyle\textbf{Y}^{*}=\alpha\textbf{1}$
$\displaystyle+X{\mbox{\boldmath$\beta$}}+\sum_{m=1}^{M}\Xi_{m}{\mbox{\boldmath$\zeta$}}_{m}$
$\displaystyle+\sum_{f=1}^{F}W_{f}\Theta_{f}{\mbox{\boldmath$\lambda$}}_{f}+{\mbox{\boldmath$\epsilon$}}^{*},$

where $\textbf{Y}^{*}$ and ${\mbox{\boldmath$\epsilon$}}^{*}$ are both
$n\times 1$ vectors and the remaining model components are as previously
defined. We place the same mean-zero shrinkage priors and hyper-priors on the
probit model as we do for the Gaussian model:
${\mbox{\boldmath$\zeta$}}_{m}\sim N(\textbf{0},\omega_{m}\mathcal{P}_{m})$,
${\mbox{\boldmath$\lambda$}}_{f}\sim N(\textbf{0},\eta_{f}\Delta_{f})$,
$\omega_{m}\sim IG(a_{\omega,m},b_{\omega,m})$, and $\eta_{m}\sim
IG(a_{\eta,m},b_{\eta,m})$. We also set the hyper-parameters
$a_{\omega,m},b_{\omega,m},a_{\eta,m},$ and $b_{\eta,m}$, to something small,
0.01 or less for this model as well.

The addition of the latent variable, $y_{i}^{*}$, to model makes this model
tractable and the resulting full conditionals are fully identifiable. Thus,
the CAVI algorithm is a reasonable approach to performing estimation in the
variational context. The $q$ densities that result from the MFVB approximation
are

$\displaystyle q({\mbox{\boldmath$\theta$}})$ $\displaystyle\sim
N\left[{\mbox{\boldmath$\mu$}}_{q({\mbox{\boldmath$\theta$}})},{\mbox{\boldmath$\Sigma$}}_{q({\mbox{\boldmath$\theta$}})}\right]$
$\displaystyle q(\textbf{Y}^{*}|\textbf{Y}=0)$ $\displaystyle\sim
N_{(-\infty,0)}\left[{\mbox{\boldmath$\mu$}}_{q(\textbf{Y}^{*})},I_{n_{0}\times
n_{0}}\right]$ $\displaystyle q(\textbf{Y}^{*}|\textbf{Y}=1)$
$\displaystyle\sim
N_{[0,\infty)}\left[{\mbox{\boldmath$\mu$}}_{q(\textbf{Y}^{*})},I_{n_{1}\times
n_{1}}\right]$ $\displaystyle q(\omega_{1})$ $\displaystyle\sim
IG\left[a_{\omega,1}+\frac{K}{2},B_{q(\omega_{1})}\right],$
$\displaystyle\cdots$ $\displaystyle q(\omega_{M})$ $\displaystyle\sim
IG\left[a_{\omega,M}+\frac{K}{2},B_{q(\omega_{M})}\right],$ $\displaystyle
q(\eta_{1})$ $\displaystyle\sim
IG\left[a_{\eta,1}+\frac{L}{2},B_{q(\eta_{1})}\right],$ $\displaystyle\cdots$
$\displaystyle q(\eta_{F})$ $\displaystyle\sim
IG\left[a_{\eta,F}+\frac{L}{2},B_{q(\eta_{F})}\right],$

where $n_{0}$ is the number of failures and $n_{1}$ is the number of successes
observed in $y_{i}$. The distributions $q(\textbf{Y}^{*}|\textbf{Y}=0)$ and
$q(\textbf{Y}^{*}|\textbf{Y}=1)$ are conditional on the components of Y that
equal zero and one, respectively, and the subscripts $(-\infty,0)$ and
$[0,\infty)$ denote truncation of the distribution to those ranges. Thus,
$q(\textbf{Y}^{*}|\textbf{Y}=0)$ and $q(\textbf{Y}^{*}|\textbf{Y}=1)$ are
truncated normal distributions. Because the Markov blanket for the penalized
components of the binary model are the same as in the Gaussian model, the $q$
densities for the penalty terms are the same in both, see Figure 2 which
contains the DAG for the binary model.

$\textbf{Y}^{*}$Y$\alpha$$\beta$$\zeta$$\lambda$$\omega$$\eta$ Figure 2: DAG
for probit variational additive model

Algorithm 2 describes the CAVI-based estimation approach. Convergence is
determined when changes in $\log[\underline{p}(\textbf{Y};q)]$ become
negligible:

$\displaystyle\log[$
$\displaystyle\underline{p}(\textbf{Y};q)]=\textbf{Y}^{\prime}\log[\Phi\\{\textbf{C}{\mbox{\boldmath$\mu$}}_{q({\mbox{\boldmath$\theta$}})}\\}]$
$\displaystyle+(\textbf{1}-\textbf{Y})^{\prime}\log[\textbf{1}-\Phi\\{\textbf{C}{\mbox{\boldmath$\mu$}}_{q({\mbox{\boldmath$\theta$}})}\\}]-\frac{1}{2}\log\left|{\mbox{\boldmath$\Sigma$}}_{q({\mbox{\boldmath$\theta$}})}\right|$
$\displaystyle-\frac{1}{2\sigma_{b}^{2}}\left[{{\mbox{\boldmath$\mu$}}_{q(\alpha,{\mbox{\boldmath$\beta$}})}}^{\prime}{\mbox{\boldmath$\mu$}}_{q(\alpha,{\mbox{\boldmath$\beta$}})}+\text{tr}\left\\{{\mbox{\boldmath$\Sigma$}}_{q(\alpha,{\mbox{\boldmath$\beta$}})}\right\\}\right]$
$\displaystyle+\sum_{m=1}^{M}\bigg{[}a_{\omega,m}\log(b_{\omega,m})$
$\displaystyle\hskip
30.0pt-\left(a_{\omega,m}+\frac{K}{2}\right)\log(B_{q(\omega_{m})})$
$\displaystyle\hskip
30.0pt+\log\left\\{\Gamma\left(a_{\omega,m}+\frac{K}{2}\right)\right\\}$
$\displaystyle\hskip 30.0pt-\log\\{\Gamma(a_{\omega,m})\\}\bigg{]}$
$\displaystyle+\sum_{f=1}^{F}\bigg{[}a_{\eta,f}\log(b_{\eta,f})$
$\displaystyle\hskip
30.0pt-\left(a_{\eta,f}+\frac{L}{2}\right)\log(B_{q(\eta_{f})})$
$\displaystyle\hskip
30.0pt+\log\left\\{\Gamma\left(a_{\eta,f}+\frac{L}{2}\right)\right\\}$
$\displaystyle\hskip 30.0pt-\log\\{\Gamma(a_{\eta,f})\\}\bigg{]},$

where $\Phi()$ denotes the cdf of a standard Gaussian distribution. The
algorithm can accommodate an arbitrary number of smoothed and functional
effects when fitting the binary model described by Equations (5) and (6).

Algorithm 2 Estimate components from Equations (5) and (6). The abbreviations
‘bd’ is for block diagonal. $\Phi()$ denotes the cdf of a standard Gaussian
distribution and $\phi()$ denotes the pdf.

1:$B_{q(\omega_{1})},\ldots,B_{q(\omega_{M})},B_{q(\eta_{1})},\ldots,B_{q(\eta_{F})}>0$,
and $\epsilon>0$, small. Set ${\mbox{\boldmath$\mu$}}_{q(\textbf{Y}^{*})}$ to
something reasonable, e.g. 0.

2:while $\Delta\log[\underline{p}(\textbf{Y};q)]>\epsilon$ do

3: $\mathcal{D}\leftarrow\text{bd}\left\\{(\sigma_{b}^{2})^{-1}I_{p\times
p},\frac{a_{\omega,1}+\frac{K}{2}}{B_{q(\omega_{1})}}\mathcal{P}_{1},\cdots,\frac{a_{\omega,M}+\frac{K}{2}}{B_{q(\omega_{M})}}\mathcal{P}_{M},\frac{a_{\eta,1}+\frac{L}{2}}{B_{q(\eta_{1})}}\Delta_{1},\cdots,\frac{a_{\eta,F}+\frac{L}{2}}{B_{q(\eta_{F})}}\Delta_{M}\right\\}$

4:
${\mbox{\boldmath$\Sigma$}}_{q({\mbox{\boldmath$\theta$}})}\leftarrow\left[\textbf{C}^{\prime}\textbf{C}+\mathcal{D}\right]^{-1}$

5:
${\mbox{\boldmath$\mu$}}_{q({\mbox{\boldmath$\theta$}})}\leftarrow{\mbox{\boldmath$\Sigma$}}_{q({\mbox{\boldmath$\theta$}})}\textbf{C}^{\prime}\textbf{Y}$

6:
${\mbox{\boldmath$\mu$}}_{q(\textbf{Y}^{*})}\leftarrow\textbf{C}^{\prime}{\mbox{\boldmath$\mu$}}_{q({\mbox{\boldmath$\theta$}})}+\phi\left[\textbf{C}^{\prime}{\mbox{\boldmath$\mu$}}_{q({\mbox{\boldmath$\theta$}})}\right]\bigg{/}\Phi\left[\textbf{C}^{\prime}{\mbox{\boldmath$\mu$}}_{q({\mbox{\boldmath$\theta$}})}\right]^{\textbf{Y}}\left[\Phi\left\\{\textbf{C}^{\prime}{\mbox{\boldmath$\mu$}}_{q({\mbox{\boldmath$\theta$}})}\right\\}-\textbf{1}\right]^{(\textbf{1}-\textbf{Y})}$

7: if $M\geq 1$ then

8: loop$\ m=1,\ldots,M$

9: $B_{q(\omega_{m})}\leftarrow
b_{\omega,m}+\frac{1}{2}\left[{{\mbox{\boldmath$\mu$}}_{q({\mbox{\boldmath$\zeta$}}_{m})}}^{\prime}{\mbox{\boldmath$\mu$}}_{q({\mbox{\boldmath$\zeta$}}_{m})}+\text{tr}\left\\{\frac{a_{\omega,m}+\frac{K}{2}}{B_{q(\omega_{m})}}\mathcal{P}_{m}\right\\}\right]$

10: end loop

11: end if

12: if $F\geq 1$ then

13: loop$\ f=1,\ldots,F$

14: $B_{q(\eta_{f})}\leftarrow
b_{\eta,f}+\frac{1}{2}\left[{{\mbox{\boldmath$\mu$}}_{q({\mbox{\boldmath$\lambda$}}_{f})}}^{\prime}{\mbox{\boldmath$\mu$}}_{q({\mbox{\boldmath$\lambda$}}_{f})}+\text{tr}\left\\{\frac{a_{\eta,f}+\frac{L}{2}}{B_{q(\eta_{f})}}\Delta_{f}\right\\}\right]$

15: end loop

16: end if

17:end while

##  3 Variational ZLS Tests

Zhang et al. (2000) and Zhang and Lin (2003) discuss hypothesis testing for
comparing two or more semiparametric curves. The test statistic, which we
refer to as the ZLS test, requires the existence of a vector, $\textbf{c}(z)$,
such that the estimate of a sermiparametric function, $s(z)$, can be written
as $\hat{s}(z)=\textbf{c}(z)^{\prime}\textbf{Y}$ at some value
$z\in[\min_{i}(z_{i}),\max_{i}(z_{i})]$. We develop the testing framework
first for the estimated variational smoothed functions and then for
variational scalar-on-function regression curves. Without loss of generality,
each test is formulated for a generic effect and we will perform the tests one
at a time.

###  3.1 Smoothed Effects

Under the null and regardless of outcome type, we assume the smoothed
component has no effect, i.e. $H_{0}:s(z)=0\ \forall\ z$. The alternative is
then that there is a non-zero effect for some values of $z$. To test this
null, we derive a ZLS-like test which requires a vector $\textbf{c}(z)$:

######  Proposition 1.

Upon convergence of Algorithm 1, there exists a vector $\textbf{c}(z)$ such
that $\hat{s}(z)=\textbf{c}(z)^{\prime}Y$.

###### Proof.

Let ${\mbox{\boldmath$\Sigma$}}_{q({\mbox{\boldmath$\zeta$}})}$ denote the
converged covariance matrix for the smoothed effect $\zeta$ resulting from
Algorithm 1. That is,
${\mbox{\boldmath$\Sigma$}}_{q({\mbox{\boldmath$\zeta$}})}$ is the principal
submatrix of ${\mbox{\boldmath$\Sigma$}}_{q({\mbox{\boldmath$\theta$}})}$ that
corresponds to $\zeta$. Upon convergence, the quantity
$\textbf{c}(z)=\left(\frac{a_{e}+\frac{N}{2}}{B_{q(\sigma_{e}^{2})}}\right){\mbox{\boldmath$\Sigma$}}_{q({\mbox{\boldmath$\zeta$}})}\Xi$
is a vector such that $\hat{s}(z)=\textbf{c}(z)^{\prime}\textbf{Y}$. ∎

Zhang and Lin (2003) describe the implementation of the ZLS test for non-
Gaussian outcomes that formulates the vector $\textbf{c}(z)$ in terms of the
working model. In the binary case, the working model is represented by the
latent variable $\textbf{Y}^{*}$.

######  Proposition 2.

Upon convergence of Algorithm 2, there exists a vector $\textbf{c}(z)$ such
that $\hat{s}(z)=\textbf{c}(z)^{\prime}\textbf{Y}^{*}$.

###### Proof.

Let ${\mbox{\boldmath$\Sigma$}}_{q({\mbox{\boldmath$\zeta$}})}$ denote the
converged covariance matrix for the smoothed effect $\zeta$ resulting from
Algorithm 2. Thus, ${\mbox{\boldmath$\Sigma$}}_{q({\mbox{\boldmath$\zeta$}})}$
is the principal submatrix of
${\mbox{\boldmath$\Sigma$}}_{q({\mbox{\boldmath$\theta$}})}$ that corresponds
to $\zeta$. Upon convergence, the quantity
$\textbf{c}(z)={\mbox{\boldmath$\Sigma$}}_{q({\mbox{\boldmath$\zeta$}})}\Xi$
is a vector such that $\hat{s}(z)=\textbf{c}(z)^{\prime}\textbf{Y}^{*}$. ∎

The statistic for the ZLS test of a variational smoothed effect is then

$\displaystyle
G_{s}(\textbf{Y})=\int_{z_{(1)}}^{z_{(N)}}\textbf{Y}^{\prime}\textbf{c}(z)\textbf{c}(z)^{\prime}\textbf{Y}dz=\textbf{Y}^{\prime}U\textbf{Y},$
(7)

where $U=\int_{z_{(1)}}^{z_{(N)}}\textbf{c}(z)\textbf{c}(z)^{\prime}dz$,
$z_{(1)}=\min_{i}(z_{i})$, and $z_{(N)}=\max_{i}(z_{i})$. When $\textbf{c}(z)$
arises from the binary model in Algorithm 2, we replace Y with
$\textbf{Y}^{*}$ in Equation (7). Consistent with Zhang et al. (2000) and
Zhang and Lin (2003), we approximate the distribution of $G(\textbf{Y})$ and
$G_{s}(\textbf{Y}^{*})$ as a scaled chi-squared, using the Satterthwaite
approximation (Satterthwaite, 1943).

Under $H_{0}$, the mean and variance of both $G_{s}(\textbf{Y})$ and
$G_{s}(\textbf{Y}^{*})$ are $e_{s}=\textbf{E}^{\prime}U\textbf{E}+tr(UV)$ and
$\psi_{s}=2tr\left\\{(UV)^{2}\right\\}+4\textbf{E}^{\prime}UVU\textbf{E}$
where E is the mean vector of Y (or $\textbf{Y}^{*}$) and $V$ is the
covariance matrix—in the binary case, $V$ is a working covariance matrix.
Under the null, $U\textbf{E}$ is negligible, so we approximate these
quantities with $e_{s}\approx tr(UV)$ and $\psi_{s}\approx
2tr\left\\{(UV)^{2}\right\\}$. Setting the approximate versions of $e_{s}$ and
$\psi_{s}$ equal to the mean and variance of a scaled chi-squared,
$\kappa_{s}\chi^{2}_{\nu_{s}}$, results in a scaling factor of
$\kappa_{s}=\psi_{s}/(2e_{s})$ and a degrees of freedom of
$\nu_{s}=2e_{s}^{2}/\psi_{s}$. The $p$-value for the test can be approximated
by finding $Pr\left[\chi^{2}_{\nu_{s}}>G_{s}(\textbf{Y})/\kappa_{s}\right]$.

###  3.2 Functional Effects

We now extend the ZLS test to a functional effect. The null for this component
is $H_{0}:\gamma(t)=0\ \forall\ t\in\mathcal{T}=\\{t:t=t_{1},\ldots,t_{T}\\}$.
We assume that the functional effect is measured on a grid that can be equally
spaced, although it is not required to be. The test requires $\textbf{c}(t)$
such that $\hat{\gamma}(t)=\textbf{c}(t)^{\prime}\textbf{Y}$ for some time
point $t\in\mathcal{T}$.

######  Proposition 3.

Upon convergence of Algorithm 1, there exists a vector $\textbf{c}(t)$ such
that $\hat{\gamma}(t)=\textbf{c}(t)^{\prime}\textbf{Y}$.

###### Proof.

Let ${\mbox{\boldmath$\Sigma$}}_{q({\mbox{\boldmath$\gamma$}})}$ denote the
converged covariance matrix for the smoothed functional effect $\gamma$
resulting from Algorithm 1. That is,
${\mbox{\boldmath$\Sigma$}}_{q({\mbox{\boldmath$\gamma$}})}$ is the principal
submatrix of ${\mbox{\boldmath$\Sigma$}}_{q({\mbox{\boldmath$\theta$}})}$ that
corresponds to $\gamma$. Upon convergence, the quantity
$\textbf{c}(t)=\left(\frac{a_{e}+\frac{N}{2}}{B_{q(\sigma_{e}^{2})}}\right){\mbox{\boldmath$\Sigma$}}_{q({\mbox{\boldmath$\gamma$}})}W\Theta$
is a vector such that $\hat{\gamma}(t)=\textbf{c}(t)^{\prime}\textbf{Y}$. ∎

Similar to the Gaussian case, the test in the binary setting requires a vector
$\textbf{c}(t)$ such that
$\hat{\gamma}(t)=\textbf{c}(t)^{\prime}\textbf{Y}^{*}$ where $\textbf{Y}^{*}$
is from the working model.

######  Proposition 4.

Upon convergence of Algorithm 2, there exists a vector $\textbf{c}(t)$ such
that $\hat{\gamma}(t)=\textbf{c}(t)^{\prime}\textbf{Y}^{*}$.

###### Proof.

Let ${\mbox{\boldmath$\Sigma$}}_{q({\mbox{\boldmath$\gamma$}})}$ denote the
converged covariance matrix for the smoothed functional effect $\gamma$
resulting from Algorithm 2. Thus,
${\mbox{\boldmath$\Sigma$}}_{q({\mbox{\boldmath$\gamma$}})}$ is the principal
submatrix of ${\mbox{\boldmath$\Sigma$}}_{q({\mbox{\boldmath$\theta$}})}$ that
corresponds to $\gamma$. Upon convergence, the quantity
$\textbf{c}(t)={\mbox{\boldmath$\Sigma$}}_{q({\mbox{\boldmath$\gamma$}})}W\Theta$
is a vector such that $\hat{\gamma}(t)=\textbf{c}(t)^{\prime}\textbf{Y}^{*}$.
∎

The statistic for the functional effect requires only slight alteration from
Equation 7:

$\displaystyle
G(\textbf{Y})=\int_{t_{1}}^{t_{T}}\textbf{Y}^{\prime}\textbf{c}(t)\textbf{c}(t)^{\prime}\textbf{Y}dt=\textbf{Y}^{\prime}Q\textbf{Y},$

where $Q=\int_{t_{1}}^{t_{T}}\textbf{c}(t)\textbf{c}(t)^{\prime}dt$. A similar
approximation to that used for the smoothed effect results in the $p$-value
having the form
$Pr\left[\chi^{2}_{\nu_{\gamma}}>G_{\gamma}(\textbf{Y})/\kappa_{\gamma}\right]$
where $\kappa_{\gamma}=\psi_{\gamma}/(2e_{\gamma})$ and a degrees of freedom
of $\nu_{\gamma}=2e_{\gamma}^{2}/\psi_{\gamma}$. Once again, when the outcome
is binary, we replace Y with $\textbf{Y}^{*}$ in the formulation of the
statistic. Under $H_{0}$, the approximations to the mean and variance are
$e_{\gamma}\approx tr(QV)$ and $\psi_{\gamma}\approx
2tr\left\\{(QV)^{2}\right\\}$ where $V$ is the covariance of Y or working
covariance of $\textbf{Y}^{*}$, depending on the outcome type.

##  4 Empirical Study

To evaluate the tests proposed in Section 3, we consider testing a single
smoothed effect or functional effect at a time. For smoothed effects, we
generate models of the form

$\displaystyle y_{i}=\alpha+s(z_{i})+\epsilon_{i}\text{ and
}y_{i}^{*}=\alpha+s(z_{i})+\epsilon_{i}.$

In the functional case, we generate models using

$\displaystyle y_{i}$
$\displaystyle=\alpha+\int_{t\in\mathcal{T}}w_{i}(t)\gamma(t)dt+\epsilon_{i}\text{
and }$ $\displaystyle y_{i}^{*}$
$\displaystyle=\alpha+\int_{t\in\mathcal{T}}w_{i}(t)\gamma(t)dt+\epsilon_{i}.$

Under both, $\epsilon_{i}\stackrel{{\scriptstyle iid}}{{\sim}}N(0,1)$ and
$y_{i}^{*}$ is the latent variable from the binary model in Equation (5). We
examine type I error and power for sample sizes of $N=50,100,$ and 200.

Under both outcome types, the functional form of $s()$ is based on either the
negated cdf of a standard normal to give a decreasing sigmoidal effect,
$\Phi$, or the pdf of a gamma distribution with both parameters set to 2,
$\Gamma$, which produces a non-central, skewed peak. Thus,

$\displaystyle s(z_{i})=-\Phi\left(\frac{z_{i}-0.5}{0.5}\right)\text{ or
}s(z_{i})=\frac{2^{2}}{\Gamma(2)}z_{i}e^{-2z_{i}}.$

Both “true” effect types were chosen to mimic the observed smoothed effects in
the lidar and ragweed data, respectively, which we analyze in Section 5. When
generating data using a true $\Phi$ curve, we take $z_{i}$ from a $N(0,1)$ and
when generating under the $\Gamma$ curve, we take $x_{i}$ from a
$\chi_{1}^{2}$.

The functional effects are either a two peak effect constructed by adding two
Gaussian pdfs together or a seasonal effect based on a $\sin$ curve. The
curves are then either

$\displaystyle\gamma(t)$
$\displaystyle=\frac{1}{4}\sqrt{\frac{100}{2\pi}}\exp\left\\{-\frac{100}{2}\left(t-\frac{1}{4}\right)^{2}\right\\}$
$\displaystyle+\frac{1}{8}\sqrt{\frac{100}{2\pi}}\exp\left\\{-\frac{100}{2}\left(t-\frac{3}{4}\right)^{2}\right\\}\text{
or}$ $\displaystyle\gamma(t)$
$\displaystyle=\sin\left\\{\pi\left(4t-1\right)\right\\}\left(t+\frac{1222}{10000}\right).$

The grid is taken to be on the unit interval, $\mathcal{T}=(0,1)$, with
equally spaced points of size $T=50$ or 100. To generate $w_{i}(t)$, we use a
Gaussian process centered at $\gamma(t)$ with an auto-regressive 1 covariance
structure and correlation set to $\rho=0.5$.

###  4.1 Type I Error

Table 1: Type I error for smoothed effect using the nominal $\alpha$ of 0.05. Abbreviations: Out. is for Outcome, Gauss. is for Gaussian, and Bin. is for Binary. Out. | Curve | Sample Size
---|---|---
50 | 100 | 200
Gauss. | $\Phi$ | 0.019 | 0.035 | 0.050
| $\Gamma$ | 0.027 | 0.031 | 0.045
Bin. | $\Phi$ | 0.004 | 0.022 | 0.058
| $\Gamma$ | 0.007 | 0.020 | 0.031

Type I error rates for the various settings and outcome types can be found in
Tables 1 and 2 and are based on 1000 simulated datasets. In general, the
choice of knots varies by effect type (smoothed or functional) and outcome
type (Gaussian or binary). For smoothed effects on Gaussian outcomes,
preliminary testing found that $K=8$ knots worked well across all settings and
produced the well controlled type I error rates from Table 1. Under the binary
outcome, we found that $K=6$ knots was preferable for controlling type I error
across most settings. Error rates tend to increase as $N$ increases and for
the binary case can be quite small for the smallest $N$. Decreasing $K$ when
$N$ is small does give rates closer to nominal, which we study in Section 4.3.

Table 2: Type I error for functional effects using the nominal $\alpha$ of 0.05. Abbreviations: Out. is for Outcome, Gauss. is for Gaussian, and Bin. is for Binary. Out. | Curve | $T$ | Sample Size
---|---|---|---
50 | 100 | 200
Gauss. | Two Peak | 50 | 0.029 | 0.042 | 0.052
|  | 100 | 0.049 | 0.054 | 0.051
| Seas. | 50 | 0.029 | 0.042 | 0.052
|  | 100 | 0.049 | 0.054 | 0.051
Bin. | Two Peak | 50 | 0.009 | 0.016 | 0.048
|  | 100 | 0.005 | 0.012 | 0.055
| Seas. | 50 | 0.009 | 0.016 | 0.048
|  | 100 | 0.005 | 0.012 | 0.055

For the functional effects in Table 2, we observe a similar pattern to the
smoothed case. Type I error rates are smallest when $N$ is small and move
toward nominal as $N$ increases. Once again, when the outcome is binary, we
see the test is quite conservative for the two smaller sample sizes of $N=50$
and 100. Knot choice in the functional case also depends on outcome type with
$L=12$ knots performing well in preliminary testing for Gaussian outcomes and
$L=9$ performing well for binary outcomes across all settings. Lowering $L$
can improve type I error when the sample size is smaller. We explore this
further in Section 4.3.

###  4.2 Power

Table 3: Power for smoothed effects using the nominal $\alpha$ of 0.05. Abbreviations: Out. is for Outcome, Gauss. is for Gaussian, and Bin. is for Binary. Out. | Curve | $N$ | $\xi$
---|---|---|---
1 | 3 | 5
Gauss. | $\Phi$ | 50 | 0.092 | 0.962 | 1.000
|  | 100 | 0.456 | 1.000 | 1.000
|  | 200 | 0.624 | 1.000 | 1.000
| $\Gamma$ | 50 | 0.060 | 0.692 | 0.998
|  | 100 | 0.170 | 0.986 | 1.000
|  | 200 | 0.336 | 1.000 | 1.000
Bin. | $\Phi$ | 50 | 0.018 | 0.152 | 0.256
|  | 100 | 0.103 | 0.516 | 0.672
|  | 200 | 0.309 | 0.862 | 0.938
| $\Gamma$ | 50 | 0.032 | 0.262 | 0.386
|  | 100 | 0.120 | 0.568 | 0.716
|  | 200 | 0.278 | 0.848 | 0.954

To evaluate power, we scale the curves $s(z_{i})$ and $\gamma(t)$ by a factor,
$\xi=1,3,$ or 5. For the smoothed effects, consistent with the type I error
evaluation, we use $K=8$ and 6 knots for Gaussian and binary models,
respectively. For the functional effects, we use $L=12$ and 9 for Gaussian and
binary models, respectively. Tables 3 and 4 contain power for all settings
except for the binary model with a functional effect for $N=50$ since this
setting requires fewer knots.

The tests for smoothed effects gains power as $N$ increases, regardless of
outcome type (Table 3). For the same number of knots used to examine type I
error, power is lower when $N$ is small but increases as $\xi$ gets larger.
The tests tend to perform better in the Gaussian case when testing smoothed
effects, although for large $N$ the binary model achieves reasonable levels of
power. Power is similar between the two curve types.

Table 4: Power for functional effects using the nominal $\alpha$ of 0.05 when $T=50$. Abbreviations: Out. is for Outcome, Gauss. is for Gaussian, and Bin. is for Binary. Out. | Curve | $N$ | $\xi$
---|---|---|---
1 | 3 | 5
Gauss. | Two Peak | 50 | 1.000 | 1.000 | 1.000
|  | 100 | 1.000 | 1.000 | 1.000
|  | 200 | 1.000 | 1.000 | 1.000
| Seas. | 50 | 0.998 | 1.000 | 1.000
|  | 100 | 1.000 | 1.000 | 1.000
|  | 200 | 1.000 | 1.000 | 1.000
Bin. | Two Peak | 100 | 0.160 | 0.996 | 1.000
|  | 200 | 0.874 | 1.000 | 1.000
| Seas. | 100 | 0.100 | 0.912 | 0.996
|  | 200 | 0.608 | 1.000 | 1.000

For functional effects, we see that power increases quickly as $\xi$
increases, regardless of outcome type (Table 4). As with the smoothed case,
power improves as $N$ increases as well. The curve type does not have much of
an impact on power for larger values of $\xi$. The Gaussian case in particular
achieves a high level of power quite quickly. Missing from Table 4 is the
power under the $N=50$ case for the binary outcome. The number of knots that
generally gave good type I error, $L=9$, is too large for the smaller sample
case when evaluating power. Thus, we defer an exploration of this setting to
Section 4.3. Table 4 displays power when $T=50$. We also evaluated the $T=100$
case and found good power for Gaussian outcomes under all $N$ and the binary
outcomes when $N=100$ and 200. Results from that simulation are in Section 1
of the Supplementary Material.

###  4.3 Small Sample Binary Model

When the outcome is binary and $N=50$, Tables 1 and 2 suggest a very
conservative test. Our power evaluation reflects this conservatism in the
smoothed case. For functional effects, we omit the results under $L=9$ knots
due to model instability. As we now demonstrate, the binary case is sensitive
to not only the sample size but the number of knots. When $N=50$, lowering $L$
for functional effects or $K$ for smoothed effects will improve the type I
error rate and power.

Table 5: Power for binary models using the nominal $\alpha=0.05$ when $N=50$. Functional effects were generated under the $T=50$ setting. Abbreviations: Smo. is for smoothed, Func. is for functional, and Kn. is for knots—$K$ for smoothed effects, $L$ for functional effects.  Effect | Curve | Kn. | $\xi$
---|---|---|---
0 | 1 | 3 | 5
Smo. | $\Phi$ | 4 | 0.037 | 0.210 | 0.686 | 0.832
|  | 5 | 0.004 | 0.022 | 0.152 | 0.256
| $\Gamma$ | 4 | 0.018 | 0.152 | 0.582 | 0.774
|  | 5 | 0.007 | 0.032 | 0.262 | 0.386
Func. | DN | 6 | 0.034 | 0.344 | 0.994 | 1.000
|  | 7 | 0.017 | 0.166 | 0.914 | 0.920
| Seas. | 6 | 0.034 | 0.144 | 0.854 | 0.982
|  | 7 | 0.017 | 0.048 | 0.278 | 0.406

Table 5 contains type I error and power for both smoothed and functional
effects when the outcome is binary and the sample size is smaller. We consider
$K=4$ and 5 knots in the smoothed case and $L=6$ and 7 knots in the functional
cases. Type I error is the table value when $\xi=0$, otherwise the power
evaluation is the same as in Tables 3 and 4. All testing is performed at the
$\alpha=0.05$ level.

When the sample size is small, the test obtains closer to nominal type I error
with smaller numbers of knots: $K=4$ for smoothed effects and $K=6$ knots for
functional effects. The test is still conservative, but less so when compared
to the results from Tables 1 and 2. The power gain is noticeable, particularly
for the functional effects where $L=6$ knots achieves power when $N=50$
similar to what we observe for larger values of $N$. This is under the $T=50$
case. As $T$ increases, however, more knots are needed. When $T=100$, for
example, $L=7$ knots produces a closer to nominal type I error while $L=6$
knots can result in inflated type I error—see Section 1 of the Supplementary
Material for additional details.

The number of knots we use in Sections 4.1 and 4.2 work in general for most
problems. But as can be seen from Table 5 as well as the additional results in
the Supplementary Material, care must be taken when selecting the number of
knots for binary outcome models with smaller sample sizes.

##  5 Data Illustrations

Table 6: ZLS test results for the data illustrations. All tests are performed at the $\alpha=0.05$ level. Abbreviations: Rag. denotes ragweed, Stront. denotes strontium, DiS denotes day-in-season. Data | $N$ | Variable |  $\chi^{2}$ $(\nu)$ |  $p$-value
---|---|---|---|---
Lidar | 221 | $s(\text{range}_{i})$ | 6.77 (1.32) | 0.015
Rag. | 335 | $s_{1}(\text{DiS}_{i})$ | 35.7 (1.73) | $<0.001$
|  | $s_{2}(\text{temp}_{i})$ | 1.49 (1.19) | 0.271
Stront. | 106 | $s(\text{age}_{i})$ | 3.19 (1.29) | 0.106
Flight | 26 |  $\Delta$SaO2${}_{i}(t)$ | 0.56 (1.41) | 0.601
DTI | 99 | CCA${}_{i}(t)$ | 9.02 (2.56) | 0.020

We demonstrate the use of the variational ZLS test in five data settings that
cover a wide spectrum of substantive fields. We consider three smoothed
estimation problems and two scalar-on-function regression problems. Two of the
five datasets have binary outcomes while the rest have Gaussian outcomes. Test
results are in Table 6 and some graphical results are provided below.
Additional graphical results are in the Supplementary Material, including
graphs of the changes in the ELBO for each model. Code to implement our models
is available for download as an R package from
https://github.com/markjmeyer/fitVAM.

###  5.1 Smoothing Problems

####  5.1.1 Lidar Data

Figure 3: Estimated semiparametric curve, $\hat{s}(\text{range}_{i})$, for the
lidar data model.

The lidar data is a classic illustration of semiparametric regression used by
Ruppert et al. (2003) to illustrate the concept. It features 221 observations
from a light detection and ranging experiment measuring the log of the ratio
of light received from two laser sources. The predictor of interest in this
example is the distance travelled before the light is reflected back to its
source. The data comes from a text on spectroscopic techniques (Sigrist, 1994)
and is available in the R package SemiPar (Wand, 2018). We fit the model

$\displaystyle\log(\text{ratio}_{i})=\alpha+s(\text{range}_{i})+\epsilon_{i},$

where $\epsilon_{i}\sim N(0,\sigma^{2})$, using Algorithm 1 and use the
variational ZLS test from Section 3 to perform a test of the null hypothesis
$H_{0}:s(\text{range}_{i})=0$ for all values of $\text{range}_{i}$.

The results of the test are in Table 6. The estimated curve is in Figure 3.
Given the smoothed effect and Gaussian outcome, we use $K=8$ knots in
modeling. The resulting curve captures the nature of the relationship in the
data quite well, Figure 3. Based on the ZLS test results in Table 6, there
appears to be a significant association between range and log-ratio
($\chi^{2}_{1.32}=6.77$, $p=0.015$).

####  5.1.2 Ragweed Data

Figure 4: Estimated semiparametric curve, $\hat{s_{1}}(\text{DiS}_{i})$, for
the ragweed data model.

Stark et al. (1997) describe determinants of ragweed levels for 335 days in
the city of Kalamazoo, Michigan, USA. Included in the data are the day in the
current ragweed pollen season the measurement was taken on and temperature,
among others. The former variable gives a sense of how far along in the
pollen, i.e. allergy, season the measurement was taken. The data is available
in the R package SemiPar (Wand, 2018). We jointly estimate smoothed effects
for both variables predicting the square root of ragweed levels:

$\displaystyle\sqrt{\text{ragweed}_{i}}=\alpha+s_{1}(\text{DiS}_{i})+s_{2}(\text{temp.}_{i})+\epsilon_{i},$

where $\epsilon_{i}\sim N(0,\sigma^{2})$, DiS denotes day in season, and temp.
denotes temperature. Estimation is performed using Algorithm 1 on $K=8$ knots
and we test both the null of $H_{0}:s(\text{DiS}_{i})=0$ and
$H_{0}:s(\text{temp.}_{i})=0$.

Table 6 displays the results of the variational ZLS test while Figure 4
contains the estimated curve for day in season, i.e.
$\hat{s}_{1}(\text{DiS}_{i})$. The graph of the estimated effect of
temperature is in the Supplementary Material. Even in the presence of
temperature, day in season is a highly significant predictor of the ragweed
level at the nominal $\alpha$ level ($\chi^{2}_{1.73}=35.7$, $p<0.001$).
Temperature, on the other hand, is not a significant predictor. If we fit the
model with temperature alone, we find moderate significance at the
$\alpha=0.05$ level ($\chi^{2}_{1.11}=4.69$, $p=0.036$).

####  5.1.3 Strontium Data

Figure 5: Estimated semiparametric curve, $\hat{s_{1}}(\text{age}_{i})$, for
the ragweed data model.

We now examine a binary outcome model from a study of fossils looking at the
ratios of strontium isotopes from mid-cretaceous deep-sea sections (Bralower
et al., 1997). The study examined the ratio of two strontium isotopes in
relation to the age (in millions of years) of the sample. We discretize the
outcome to below median strontium ratio (coded as 0) versus above median
strontium ratio (coded as 1) and analyze the effect of age using the latent
model

$\displaystyle y_{i}^{*}=\alpha+s(\text{age}_{i})+\epsilon_{i}^{*},$

where $\epsilon_{i}^{*}\sim N(0,1)$.

We obtain estimates via Algorithm 2 on $K=6$ knots and perform the ZLS to test
the null that $H_{0}:s(\text{age}_{i})=0$ at the nominal level. The results of
the test are in Table 6 while the estimated curve is in Figure 5. While we do
not find evidence of a statistically significant relationship between the
ratios of strontium isotopes and the age of the sample
($\chi^{2}_{1.29}=3.19$, $p=0.106$), the estimated curve fits the discretized
data quite well.

###  5.2 Functional Problems

####  5.2.1 Flight Health Data

Meyer et al. (2019) describes a study of the impact of exposure to altitude in
commercial flight on in-flight markers of heart health finding that exposure
to altitude negatively impacts some markers, including heart rate. Using this
data, we examine the impact of exposure to altitude for 310 minutes on post-
flight heart rate. Exposure to altitudes in commercial flights is known to
decrease blood oxygen saturation SaO2. We use five-minute SaO2 measurements to
examine post-flight heart rate. The study was a cross-over design with
patients having one dat at altitude and one day at sea level in a hypobaric
pressure chamber (the order of treatment arms was block randomized). To
account for the pairing, we pre-difference both the SaO2 curves and post-
experimental condition heart rate. The model is then

$\displaystyle\Delta hr_{i}$
$\displaystyle=\alpha+\int_{t\in\mathcal{T}}\Delta\text{SaO2}_{i}(t)\gamma(t)dt+\epsilon_{i},$

where $\epsilon_{i}\sim N(0,\sigma^{2})$, and $\Delta$ denotes the change in
the variable. For estimation via Algorithm 1, we employ 12 knots. We performed
the differencing by taking active treatment (at altitude) minus placebo (at
sea level). From table 6, we lack evidence to suggest there is an association
between changes in blood oxygen level in flight on changes in heart rate post-
flight, from altitude to sea level ($\chi^{2}_{1.41}=0.56$, $p=0.601$).

####  5.2.2 DTI Data

Our last illustration looks at the impact of a diffusion tensor image
fractional anisotropy tract profiles from the corpus callosum on paced
auditory serial addition test (PASAT) scores in patients with multiple
sclerosis. This data has been considered as a functional data example by
Goldsmith et al. (2011) and Goldsmith et al. (2012), among others. We first
discretize the PASAT score to below median (coded as 0) and above median
(coded as 1). Lower PASAT scores are associated with cognitive impairment.
This binary variable serves as our outcome with the latent model

$\displaystyle y_{i}^{*}$
$\displaystyle=\alpha+\int_{t\in\mathcal{T}}\text{CCA}_{i}(t)\gamma(t)dt+\epsilon_{i}^{*},$

where $\epsilon_{i}^{*}\sim N(0,1)$ and CCA denotes the diffusion tensor image
fractional anisotropy tract profiles from the corpus callosum. Although the
ZLS test on variational binary functional models can be conservative for
$N=100$ (see Table 2), we use $L=9$ knots for this analysis. Table 6 displays
the results which suggests an association between the tract profiles and
whether or not the PASAT score was high or low ($\chi^{2}_{2.56}=9.02$,
$p=0.020$).

##  6 Discusion

As variational techniques gain popularity, additional methods will be required
to ensure that statistical inference can be confidently performed within the
framework. We show that for MFVB approximations to Gaussian and binary
additive models, the CAVI algorithms admit forms that can be used to implement
a global test of a smoothed effect, upon convergence of the algorithm. We
further demonstrate that the penalized spline representation of the functional
effect in these additive models also admits a form via the CAVI algorithm
whereby one can construct a global test of a functional effect.

Our empirical evaluation demonstrates that the testing framework has good
Frequentist properties in terms of type I error rate and power. The data
illustrations show the test is applicable to a wide range of substantive
questions with sample sizes varying from 26 to 335. While the choice of knots
can impact type I error and power, particularly for the binary case, we
identify a set of reasonable choices for the number knots depending on the
outcome, effect type, sample size, and grid (in the functional case).

Our work is applied to MFVB approximations for Gaussian and binary additive
models using the product density transformation. Other transformations may
admit similar forms, for example the tangent transformation. The work of Luts
and Wand (2015), which examines semiparametric models for Poisson and negative
binomially distributed outcomes, may also be amenable. Recent work to improve
variance estimation by Giordano et al. (2018) could also present useable forms
for ZLS-like tests. These other transformations and approaches are of interest
for future work to demonstrate the applicability of the test within the
broader variational framework.

## Supplementary information

Further results from our empirical study and data illustrations are in the
Supplementary Material accompanying this manuscript.

## References

  * Albert and Chib [1993] J. Albert and S. Chib.  Bayesian analysis of binary and polychotomous response data.  _Journal of the American Statistical Association_ , 88:669–679, 1993.  doi: 10.1080/01621459.1993.10476321.
  * Bickel et al. [2013] P. Bickel, D. Choi, X. Chang, and H. Zhang.  Asymptotic normality of maximum likelihood and its variational approximation for stochastic blockmodels.  _The Annals of Statistics_ , 41:1922–1943, 2013.  doi: 10.1214/13-AOS1124.
  * Blei et al. [2017] D. M. Blei, A. Kucukelbir, and J. D. McAuliffe.  Variational inference: a review for statisticians.  _Journal of the American Statistical Association_ , 112:859–877, 2017.
  * Bralower et al. [1997] T. Bralower, P. Fullagar, C. K. Paull, G. S. Dwyer, and R. M. Leckie.  Mid-cretaceous strontium-isotope stratigraphy of deep-sea sections.  _Geological Society of America Bulletin_ , 109:1421–1442, 1997.
  * Chen et al. [2018] Y. H. Chen, B. Mukherjee, S. D. Adar, V. J. Berrocal, and B. A. Coull.  Robust distributed lag models using data adaptive shrinkage.  _Biostatistics_ , 19:461–478, 2018.
  * Crainiceanu et al. [2005] C. Crainiceanu, C. Ruppert, G. Claeskens, and M. P. Wand.  Exact likelihood ratio tests for penalised splines.  _Biometrika_ , 92:91–103, 2005.  doi: 10.1093/biomet/92.1.91.
  * Crainiceanu and Ruppert [2004] C. M. Crainiceanu and D. Ruppert.  Likelihood ratio tests in linear mixed models with one variance component.  _Journal of the Royal Statistical Society: Series B_ , 66:165–185, 2004.  doi: 10.1111/j.1467-9868.2004.00438.x.
  * Eilers and Marx [1996] P. H. C. Eilers and B. D. Marx.  Flexible smoothing with b-splines and penalties.  _Statistical Science_ , 11:89–121, 1996.  doi: 10.1214/ss/1038425655.
  * Gelman et al. [2013] A. Gelman, J. B. Carlin, H. S. Stern, D. B. Dunson, A. Vehtari, and D. B. Rubin.  _Bayesian Data Analysis_.  Chapman and Hall-CRC, Boca Raton, FL, 3${}^{\text{rd}}$ edition, 2013.
  * Giordano et al. [2018] R. Giordano, T. Broderick, and M. I. Jordan.  Covariances, robustness, and variational bayes.  _Journal of Machine Learning Research_ , 19:1–49, 2018\.
  * Goldsmith and Kitago [2016] J. Goldsmith and T. Kitago.  Assessing systematic effects of stroke on motor control using hierarchical function-on-scalar regression.  _Journal of the Royal Statistical Society, Series C_ , 65:215–236, 2016.  doi: 10.1111/rssc.12115.
  * Goldsmith et al. [2011] J. Goldsmith, J. Bobb, C. M. Crainiceanu, B. Caffo, and D. Reich.  Penalized functional regression.  _Journal of Computational and Graphical Statistics_ , 20:830–851, 2011.  doi: 10.1198/jcgs.2010.10007.
  * Goldsmith et al. [2012] J. Goldsmith, C. Crainiceanu, B. Caffo, and D. Reich.  Longitudinal penalized functional regression for cognitive outcomes on neuronal tract measurements.  _Journal of the Royal Statistical Society, Series C_ , 61:453–469, 2012.  doi: 10.1111/j.1467-9876.2011.01031.x.
  * Goldsmith et al. [2013] J. Goldsmith, S. Greven, and C. Crainiceanu.  Corrected confidence bands for functional data using principal components.  _Biometrics_ , 69:41–51, 2013.  doi: 10.1111/j.1541-0420.2012.01808.x.
  * Greven et al. [2008] S. Greven, C. M. Crainiceanu, H. Küchenhoff, and A. Peters.  Restricted likelihood ratio testing for zero variance components in linear mixed models.  _Journal of Computational and Graphical Statistics_ , 17:870–891, 2008.  doi: 10.1198/106186008X386599.
  * Hui et al. [2019] F. K. C. Hui, C. You, H. L. Shang, and S. Müller.  Semiparametric regression using variational approximations.  _Journal of the American Statistical Association_ , 114:1765–1777, 2019.  doi: 10.1080/01621459.2018.1518235.
  * Ivanescu et al. [2015] A. E. Ivanescu, A.-M. Staicu, F. Scheipl, and S. Greven.  Penalized function-on-function regression.  _Computational Statistics_ , 30:539–568, 2015.  doi: 10.1007/s00180-014-0548-4.
  * Lee and Wand [2016] C. Y. Y. Lee and M. P. Wand.  Variational methods for fitting complex bayesian mixed effects models to health data.  _Statistics in Medicine_ , 35:165–188, 2016.  doi: 10.1002/sim.6737.
  * Li et al. [2016] H. Li, S. Kozey-Keadle, V. Kipnis, and R. J. Carroll.  Longitudinal functional additive model with continuous proportional outcomes for physical activity data.  _Stat_ , 5:242–250, 2016.  doi: 10.1002/sta4.121.
  * Luts and Wand [2015] J. Luts and M. P. Wand.  Variational inference for count response semiparametric regression.  _Bayesian Analysis_ , 10:991–1023, 2015.
  * Ma et al. [2021] W. Ma, L. Xiao, and B. Liu.  A functional mixed model for scalar on function regression with application to a functional MRI study.  _Biostatistics_ , 22:439–454, 2021.  doi: 10.1093/biostatistics/kxz046.
  * Malloy et al. [2010] E. J. Malloy, J. S. Morris, S. D. Adar, H. Suh, D. R. Gould, and B. A. Coull.  Wavelet-based functional linear mixed models: an application to measurement error-corrected distributed lag models.  _Biostatistics_ , 11:432–452, 2010.  doi: 10.1093/biostatistics/kxq003.
  * McLean et al. [2014] M. W. McLean, G. Hooker, A.-M. Staicu, F. Scheipl, and D. Ruppert.  Functional generalized additive models.  _Journal of Computational and Graphical Statistics_ , 23:249–269, 2014.  doi: 10.1080/10618600.2012.729985.
  * Meyer et al. [2015] M. J. Meyer, B. A. Coull, F. Versace, P. Cinciripini, and J. S. Morris.  Bayesian function-on-function regression for multi-level functional data.  _Biometrics_ , 71:563–574, 2015.  doi: 10.1111/biom.12299.
  * Meyer et al. [2019] M. J. Meyer, I. Mordukhovich, G. A. Wellenius, M. A. Mittleman, J. P. McCracken, B. A. Coull, and E. McNeely.  Changes in heart rate and rhythm during a crossover study of simulated commercial flight in older and vulnerable participants.  _Frontiers Physiology_ , 10:1339, 2019.  doi: 10.3389/fphys.2019.01339.
  * Meyer et al. [2022] M. J. Meyer, J. S. Morris, R. P. Gazes, and B. A. Coull.  Ordinal probit functional outcome regression with application to computer-use behavior in rhesus monkeys.  _Annals of Applied Statistics_ , 16:537–550, 2022.  doi: 10.1214/21-AOAS1513.
  * Ormerod and Wand [2010] J. T. Ormerod and M. P. Wand.  Explaining variational approximations.  _The American Statistician_ , 64:140–153, 2010.
  * Ramsay and Dalzell [1991] J. O. Ramsay and C. J. Dalzell.  Some tools for functional data analysis.  _Journal of the Royal Statistical Society, Series B_ , 53:539–561, 1991.  doi: 10.1111/j.2517-6161.1991.tb01844.x.
  * Ramsay and Silverman [2005] J. O. Ramsay and B. W. Silverman.  _Functional Data Analysis_.  Springer, New York, 2${}^{\text{nd}}$ edition, 2005.
  * Ruppert et al. [2003] D. Ruppert, M. P. Wand, and R. J. Carroll.  _Semiparametric regression_.  Cambridge University Press, Cambridge, United Kingdom, 2003.
  * Satterthwaite [1943] F. E. Satterthwaite.  An approximate distribution of variance components.  _Biometrika_ , 2:110–114, 1943.
  * Scheipl et al. [2015] F. Scheipl, A.-M. Staicu, and S. Greven.  Functional additive mixed models.  _Journal of Computational and Graphical Statistics_ , 24:477–501, 2015.  doi: 10.1080/10618600.2014.901914.
  * Scheipl et al. [2016] F. Scheipl, J. Gertheiss, and S. Greven.  Generalized functional additive mixed models.  _Electronic Journal of Statistics_ , 10:1455–1492, 2016\.  doi: 10.1214/16-EJS1145.
  * Sigrist [1994] M. Sigrist, editor.  _Air Monitoring by Spectroscopic Techniques (Chemical Analysis Series)_ , volume 197.  Wiley, New York, 1994.
  * Stark et al. [1997] P. C. Stark, L. M. Ryan, J. L. McDonald, and H. A. Burge.  Using meteorologic data to model and predict daily ragweed pollen levels.  _Aerobiologia_ , 13:177–184, 1997.
  * Wand [2018] M. Wand.  _SemiPar: Semiparametic Regression_ , 2018.  URL https://CRAN.R-project.org/package=SemiPar.  R package version 1.0-4.2.
  * Wand [2017] M. P. Wand.  Fast approximate inference for arbitrarily large semiparametric regression models via message passing.  _Journal of the American Statistical Association_ , 112:137–168, 2017.  doi: 10.1080/01621459.2016.1197833.
  * Wang and Blei [2019] Y. Wang and D. M. Blei.  Frequentist consistency of variational bayes.  _Journal of the American Statistical Association_ , 114:1147–1161, 2019.
  * Westling and McCormick [2019] T. Westling and T. H. McCormick.  Beyond prediction: A framework for inference with variational approximations in mixture models.  _Journal of Computational and Graphical Statistics_ , 28:778–789, 2019.  doi: 10.1080/10618600.2019.1609977.
  * Wood [2017] S. N. Wood.  _Generalized Additive Models: An Introduction with R_.  Chapman and Hall/CRCs, New York, 2${}^{\text{nd}}$ edition, 2017.  doi: 10.1201/9781315370279.
  * Yang and Yang [2024] F. Yang and Y. Yang.  A sparse estimate based on variational approximations for semiparametric generalized additive models.  _Computational Statistics_ , 2024.  doi: 10.1007/s00180-024-01485-2.
  * Yang et al. [2020] Y. Yang, D. Pati, and A. Bhattacharya.  $\alpha$-variational inference with statistical guarantees.  _The Annals of Statistics_ , 48:886–905, 2020.  doi: 10.1214/19-AOS1827.
  * You et al. [2014] C. You, J. T. Ormerod, and S. Müller.  On variational bayes estimation and variational information criteria for linear regression models.  _Australian & New Zealand Journal of Statistics_, 56:73–87, 2014.
  * Zhang and Zhou [2020] A. Y. Zhang and H. H. Zhou.  Theoretical and computational guarantees of mean field variational inference for community detection.  _Annals of Statistics_ , 48:2575–2598, 2020.
  * Zhang and Lin [2003] D. Zhang and X. Lin.  Hypothesis testing in semiparametric additive mixed models.  _Biostatistics_ , 4:57–74, 2003.
  * Zhang et al. [2000] D. Zhang, X. Lin, and M. Sowers.  Semiparametric regression for periodic longitudinal hormone data from multiple menstrual cycles.  _Biometrics_ , 56:31–39, 2000.


119–126

# Deriving fundamental parameters of millisecond pulsars via AIC in white
dwarfs

A. Taani1† Y.H. Zhao1 A. Moraghan2 1National Astronomical Observatories,
Chinese Academy of Sciences, Beijing 100012, China †email<EMAIL_ADDRESS>
2Center for Galaxy Evolution Research and Department of Astronomy, Yonsei
University, Seoul 120-749, Republic of Korea

(2012)

###### Abstract

We present a study of the observational properties of Millisecond Pulsars
(MSPs) by way of their magnetic fields, spin periods and masses. These
measurements are derived through the scenario of Accretion Induced Collapse
(AIC) of white dwarfs (WDs) in stellar binary systems, in order to provide a
greater understanding of the characteristics of MSP populations. In addition,
we demonstrate a strong evolutionary connection between neutron stars and WDs
with binary companions from a stellar binary evolution perspective via the AIC
process.

###### keywords:

Neutron stars, white dwarfs, cataclysmic variables, fundamental parameters.

††volume: 290††journal: Feeding compact objects: Accretion on all
scales††editors: C.M. Zhang, T. Belloni, M. Mendez, and S.N Zhang, eds.

##  1 Introduction

Observable parameters of binary Millisecond Pulsars (MSPs), e.g. mass of the
pulsar, mass of the companion, spin period, orbital period, eccentricity,
etc., are used to probe the past accretion history of the MSPs. The purpose of
this proceeding is to demonstrate how to infer some of the observable
quantities (spin period, magnetic field and mass) during the Accretion Induced
Collapse (AIC) of a white dwarf (WD) on its way to become a member of the MSP
family.

Figure 1: Left: The distribution of MSPs in terms of spin-period (observed
data from [Manchester et al.(2005)]. Right: Mass distribution of massive CVs
[Ritter & Kolb (2011), (Ritter & Kolb 2011)].

##  2 Evolution of spin period

We are able to determine the spin period of MSPs originating from WDs. The
process begins with simple Keplerian frequency. We assume that the angular
velocity of the NS is equal to the Keplerian angular velocity, $v_{K}$, of the
magnetosphere, at roughly the Alfv$\acute{\textrm{e}}$n surface,

$v_{K}\propto R_{NS}^{-3/2}\rightarrow P_{MSP}\sim R_{NS}^{3/2}$ (1)

from which we obtain the spin of the MSP, $P_{MSP}$, as a function of the WD
spin, $P_{WD}$.

$P_{MSP}\sim P_{WD,min}\left(\frac{R_{NS}}{R_{WD}}\right)^{3/2}$ (2)

where $P_{WD,min}$ is the minimum spin period of the standard WD [Warner
(1995), (Warner 1995)]. Assuming $P_{WD,min}$ $\sim 30$ s, $R_{NS}=10$ km and
$R_{WD}=1000$ km, we thus obtain, $P_{MSP}\sim 1ms.$ Fig. 1 (left) shows the
observed and derived spin period distributions of MSPs. As for the observed
MSPs, they show a relatively Gaussian distribution. According to these
distributions, the ratio of MSPs originating from CVs is about $\sim 10\%$.
This result agrees with some theoretical predictions such as those by [Warner
(1995)] and [Warner & Woudt (2002)].

##  3 Magnetic field

To investigate the correlation between the magnetic fields of MSPs with bottom
fields of CVs (where the field is partially restructured due to accretion), we
follow the model of [Zhang et al. (2009)]. The magnetic flux is assumed to be
conserved. This corresponds to the magnetic fields produced by AIC.

$B_{NS}=B_{f,WD}\times\left(\frac{R_{WD}}{R_{NS}}\right)^{2}$ (3)

If we adopt $R_{NS}$ = $15\times 10^{5}$ cm and $B_{f,WD}\sim 10^{3}$G in CVs,
the minimum value, $B_{NS}\sim 10^{8.5-9}G.$

##  4 Mass

The sample of CVs whose masses we have considered is the set of binary systems
collected by [Ritter & Kolb (2011)]. Among them we have 26 massive CVs in the
range ${\rm~{}1.0-1.3M_{\odot}}$. Fig. 1 (right) shows the relatively Gaussian
distribution, with mean at ${\rm~{}M_{CV}\sim 1.1M_{\odot}}$. A summary of the
known properties of these systems is given in Table 1 of [Taani et al.
(2012)]. Note that the AIC process leads to a MSP with mass less than
Chandrasekhar limit [Zhang et al. (2011)]. This provides evidence for the AIC
in massive CVs and evolutionary hypotheses of MSP birthrate.

##  5 Summary and Conclusions

  1. 1.

CVs would be invoked via the capability of producing a significant portion of
the MSPs via the AIC process, a regime which may be unattainable by normal
channels.

  2. 2.

We find that the quantitative implications of our calculations are that we
estimate the expected $P$ in the observed of MSPs which could have originated
from CVs to be $\sim 10\%$.

  3. 3.

We further find that the predictions of some parameters after AIC process for
the average levels are consistent with the observed MSP population. Future
work will consider other quantities e.g. orbital period, eccentricity, and
mass ratio (q) using more data sets from SDSS.

We are grateful for the discussion with Cole Miller. This research has been
supported by NBRPC (2009CB824800, 2012CB821800) and the NSFC (10773017,
11173034).

## References

  * [Manchester et al.(2005)] Manchester, R.N., Hobbs, G.B., & Teoh, A., et al. 2005 AJ 129, 1993
  * [Ritter & Kolb (2011)] Ritter, H. & Kolb, U., 2011, VizieR Online Data Catalog 1, 2018
  * [Taani et al. (2012)] Taani, A., Zhang, C.M., & Al-Wardat, M., et al. 2012  Ap&SS 340, 147
  * [Warner (1995)] Warner, B., 1995, Cataclysmic Variable Stars. Cambridge Astrophysics Series vol. 28
  * [Warner & Woudt (2002)] Warner, B. & Woudt, P.A., 2002 The Physics of Cataclysmic Variables and Related Objects. Astronomical Society of the Pacific Conference Series vol. 261
  * [Zhang et al. (2009)] Zhang, C.M., Wickramasinghe, D.T., & Ferrario, L., et al. 2009 MNRAS 397, 2208
  * [Zhang et al. (2011)] Zhang, C.M., Wang, J., & Zhao, Y. H., et al. 2011 A&A 527, 83


* [64] E. Mossel.  Reconstruction on trees: beating the second eigenvalue.  Ann. Appl. Probab., 11(1):285–300, 2001.
  * [65] E. Mossel.  Gaussian bounds for noise correlation of functions and tight analysis of long codes.  In Foundations of Computer Science, 2008 (FOCS 08), pages 156–165. IEEE, 2008.
  * [66] E. Mossel.  Arrow‘s impossibility theorem without unanimity.  Posted on Arxiv 0901.4727, 2009.
  * [67] E. Mossel.  Gaussian bounds for noise correlation of functions.  GAFA, 19:1713–1756, 2010.
  * [68] E. Mossel.  A quantitative arrow theorem.  Probability Theory and Related Fields, 154(1):49–88, 2012.
  * [69] E. Mossel.  Gaussian Bounds for Noise Correlation of Resilient Functions.  Israel Journal of Mathematics, 235:111–137, 2020.
  * [70] E. Mossel and J. Neeman.  Robust optimality of Gaussian noise stability.  J. Eur. Math. Soc. (JEMS), 17(2):433–482, 2015.
  * [71] E. Mossel, J. Neeman, and O. Tamuz.  Majority dynamics and aggregation of information in social networks.  Journal of Autonomous Agents and Multi-Agent Systems, (3):408–429, 2014.
  * [72] E. Mossel and R. O‘Donnell.  On the noise sensitivity of monotone functions.  Random Structures Algorithms, 23(3):333–350, 2003.
  * [73] E. Mossel, R. O‘Donnell, and K. Oleszkiewicz.  Noise stability of functions with low influences: invariance and optimality (extended abstract).  In 46th Annual IEEE Symposium on Foundations of Computer Science (FOCS 2005), 23-25 October 2005, Pittsburgh, PA, USA, Proceedings, pages 21–30. IEEE Computer Society, 2005.
  * [74] E. Mossel, R. O‘Donnell, and K. Oleszkiewicz.  Noise stability of functions with low influences: invariance and optimality.  Annals of Mathematics, 171(1):295–341, 2010.
  * [75] E. Mossel, R. O‘Donnell, O. Regev, J. E. Steif, and B. Sudakov.  Non-interactive correlation distillation, inhomogeneous Markov chains, and the reverse Bonami-Beckner inequality.  Israel J. Math., 154:299–336, 2006.
  * [76] E. Mossel, K. Oleszkiewicz, and A. Sen.  On reverse hypercontractivity.  Geometric and Functional Analysis, 23(3):1062–1097, 2013.
  * [77] E. Mossel and M. Z. Racz.  A quantitative gibbard-satterthwaite theorem without neutrality.  In H. J. Karloff and T. Pitassi, editors, STOC. Proceedings of the 44th Symposium on Theory of Computing Conference, STOC 2012, New York, NY, USA, May 19 - 22, 2012, pages 1041–1060. ACM, 2012.
  * [78] E. Mossel and M. Z. Racz.  A quantitative gibbard-satterthwaite theorem without neutrality.  Combinatorica, 35(3):317–387, 2015.
  * [79] E. Mossel and O. Schramm.  Representations of general functions using smooth functions, 2008.  Unpublished manuscript.
  * [80] J. Neeman et al.  A multidimensional version of noise stability.  Electronic Communications in Probability, 19:1–10, 2014.
  * [81] I. Nehama.  Approximately classic judgement aggregation.  Annals of Mathematics and Artificial Intelligence, 68(1-3):91–134, 2013.
  * [82] E. Nelson.  The free Markoff field.  J. Functional Analysis, 12:211–227, 1973.
  * [83] J. Neyman and E. S. Pearson.  Ix. on the problem of the most efficient tests of statistical hypotheses.  Philosophical Transactions of the Royal Society of London. Series A, Containing Papers of a Mathematical or Physical Character, 231(694-706):289–337, 1933.
  * [84] R. G. Niemi and H. F. Weisberg.  A mathematical solution for the probability of paradox of voting.  Behavioral Science, 13:317–323, 1968.
  * [85] N. Nisan.  Algorithmic game theory / economics. postname: From arrow to fourier, 2009\.
  * [86] R. O‘Donnell.  Hardness amplification within np.  In Proceedings of the 34th ACM Symposium on Theory of Computing, 2002.
  * [87] R. O‘Donnell.  Analysis of boolean functions.  Cambridge University Press, 2014.
  * [88] R. O‘Donnell, R. Servedio, L.-Y. Tan, and A. Wan.  A regularity lemma for low noisy-influences, 2010.  Unpublished manuscript.
  * [89] V. Roychowdhury, K.-Y. Siu, A. Orlitsky, and T. Kailath.  A geometric approach to threshold circuit complexity.  In Proceedings of the fourth annual workshop on Computational learning theory, pages 97–111. Morgan Kaufmann Publishers Inc., 1991.
  * [90] L. Russo.  An approximate zero-one law.  Zeitschrift für Wahrscheinlichkeitstheorie und verwandte Gebiete, 61(1):129–139, 1982.
  * [91] M. A. Satterthwaite.  Strategy-proofness and Arrow‘s Conditions: Existence and Correspondence Theorems for Voting Procedures and Social Welfare Functions.  J. of Economic Theory, 10:187–217, 1975.
  * [92] W. Sheppard.  On the application of the theory of error to cases of normal distribution and normal correlation.  Phil. Trans. Royal Soc. London, 192:101–168, 1899.
  * [93] M. Talagrand.  On Russo‘s approximate 0-1 law.  Annals of Probability, 22:1576–1587, 1994.
  * [94] Wikipedia contributors.  Marquis de condorcet — Wikipedia, the free encyclopedia, 2019.  [Online; accessed 2019].
  * [95] R. Wilson.  Social choice theory without the pareto principle.  Journal of Economic Theory, 5(3):478–486, 1972.