Pacific Symposium on Biocomputing 2024

Scalar-Function Causal Discovery for Generating Causal Hypotheses with
Observational Wearable Device Data

Valeriya Rogovchenko, Austin Sibu, and Yang Ni†
Department of Statistics, Texas A&M University,
College Station, TX 77843, USA† E-mail: yni@stat.tamu.edu
Digital health technologies such as wearable devices have transformed health data analytics, providing continuous, high-resolution functional data on various health metrics, thereby
opening new avenues for innovative research. In this work, we introduce a new approach forgenerating causal hypotheses for a pair of a continuous functional variable (e.g., physical activities recorded over time) and a binary scalar variable (e.g., mobility condition indicator).
Our method goes beyond traditional association-focused approaches and has the potential to reveal the underlying causal mechanism. We theoretically show that the proposed
scalar-function causal model is identifiable with observational data alone. Our identifiability
theory justifies the use of a simple yet principled algorithm to discern the causal relationshipby comparing the likelihood functions of competing causal hypotheses. The robustness andapplicability of our method are demonstrated through simulation studies and a real-world
application using wearable device data from the National Health and Nutrition Examination
Survey.
Keywords: Causal identifiability, digital health, NHANES, observational data, wearable device.

1. Introduction
The rise of wearable devices has revolutionized the way we collect and analyze health data,offering an unprecedented wealth of information about human health and behavior. These devices such as accelerometers and continuous glucose monitors allow for frequent measurementof various variables over time including physical activities, sleep patterns, electrocardiogramsignals, and blood glucose levels. The availability of these measurements enables researchers to
ask questions that previously could not be answered, e.g., how to quantify the effect of physicalactivities on all-cause mortality? In these types of scenarios, often, one variable (e.g., physicalactivities) is longitudinal/functional and the other (e.g., mortality) is a scalar. Thus, manystatistical methods such as scalar-on-function regression models10,17 have been successfully
deployed to estimate the association of the scalar-function pair.
The focus of this paper is, however, different from the existing literature for modelingwearable device data. Instead of association, we investigate whether it is possible to discern thecausal relationship between a scalar and a function. More specifically, we aim to identify whichof the scalar-function pair is more likely to be the cause or effect given observational data alone.
© 2023 The Authors. Open Access chapter published by World Scientific Publishing Company anddistributed under the terms of the Creative Commons Attribution Non-Commercial (CC BY-NC)
4.0 License.

201

Pacific Symposium on Biocomputing 2024

We introduce a novel scalar-function causal discovery method to generate data-driven causal
hypotheses. Revealing the causality underlying observed data can deepen our understandingof the physical mechanism involved in the data-generating process and potentially pave theway for better health interventions and policy-making.
The field of causal discovery has seen a significant surge in interest and development
over recent years due to wide-ranging applicability across various domains.4,5,12–14,16,22 While
traditional causal discovery methods are typically tailored to handle either continuous ordiscrete variables exclusively, real-world scenarios are often far more complex. For example,
in the fields of social and health sciences, data frequently comprise a mix of different types ofvariables, necessitating more versatile approaches.
In such scenarios, one may either discard discrete data or convert continuous data into a
discrete form;9,20 either way, a lot of information contained in the original data is lost. In lightof these limitations, there have been some recent developments to discover causality for mixed
data.19,21 However, these methods have only been developed for scalar variables, which cannot
be used for functional data. To deal with functional data, some very recent works6,23 have been
proposed, which, however, cannot accomodate scalar and/or discrete data. In summary, to thebest of our knowledge, there are no existing methods that can identify causality between acontinuous functional variable and a binary scalar variable.
This paper, therefore, aims to fill this critical gap in the causal discovery literature so that
digital health researchers will have a powerful tool to identify causality in a wide range of
observational wearable device data. Our approach is based on a probabilistic causal modelthat quantifies the likelihood of each possible causal direction (from function to scalar or fromscalar to function). We theoretically establish the causal identifiability property of our model
under common causal assumptions. Equipped with the identifiability property, we can simply
identify causal directions based on likelihood functions.
We conduct simulation studies to assess the empirical identifiability of the proposed
method. In addition, to validate our method in real-world scenarios, we present an application with two variables that have a clear causal relationship. Specifically, we consider mobility
conditions and physical activities. Since it is clear that mobility issues may lead to reduced
activities, we will test whether our method can correctly identify such causal relationship without prior knowledge using the National Health and Nutrition Examination Survey (NHANES)data.
The rest of the paper is organized in the following way. In Section 2, we describe the proposed scalar-function causal discovery model, theoretically prove that the causal relationshipis identifiable, and develop a likelihood-based estimation procedure. In Section 3, we evaluate
the proposed method through various simulations as well as a real wearable device dataset
from NHANES, demonstrating its capability to correctly identify the true causal relationship.
We conclude our paper with a brief discussion in Section 4.

202

Pacific Symposium on Biocomputing 2024

2. Method
2.1. Notations
We use capital letters to denote random variables and small letters to denote their realized
values. We use boldface to denote vectors or matrices and non-boldface to denote scalars.
With a slight abuse of notation, we use P (·) to denote both probability mass and density
functions, which can be understood from the context as it is determined by the type of therandom variable under consideration. Let Mn×n be the cone of n × n positive definite matrices.
2.2. Causal Probability Model
We are interested in identifying the causal relationship between two statistically dependent
random variables: a random binary variable Y ∈ {0, 1} and a random function measured on ntime points X = (X(t1 ), ..., X(tn ))⊤ ∈ Rn . One can view these functional measurements as afinite realization of an infinite stochastic process X(·) such as the Gaussian process.18
We consider two competing causal hypothesesa ,H0 : X → Y or X causes Y

vs
H1 : Y → X or Y causes X
Under each hypothesis, we will set up a probability model. Specifically, let PX→Y (X = x, Y =y) denote the probability model of H0 and PY →X (X = x, Y = y) denote the probability modelof H1 . Using the probability chain rule, we have
PX→Y (X = x, Y = y) = PX→Y (Y = y | X = x) · PX→Y (X = x),
PY →X (X = x, Y = y) = PY →X (X = x | Y = y) · PY →X (Y = y),(1)

where PX→Y (Y = y | X = x) and PX→Y (X = x) are respectively the conditional and marginal
probability distributions under H0 : X → Y and similarly PY →X (X = x | Y = y) andPY →X (Y = y) are those under H1 : Y → X . Next, we will discuss the choice of these four
probability distributions.
For the marginal distribution of Y , we assume it to be a Bernoulli distribution with successprobability ρ ∈ (0, 1),
PY →X (Y = y) = ρy (1 − ρ)1−y .

(2)For the marginal distribution of X , we assume it to be a multivariate Gaussian distribution
with mean µ ∈ Rn and covariance matrix Σ ∈ Mn×n ,PX→Y (X = x) = N (x | µ, Σ),(3)where N (x | ·, ·) is the Gaussian probability density function evaluated at x.
To model the conditional distribution of Y given X , we adopt a linear logistic regression,
log

PX→Y (Y = 1|X = x)
= α0 + x ⊤ α1 ,PX→Y (Y = 0|X = x)a Note that we are not performing null hypothesis testing. Our method is exploratory.

203

Pacific Symposium on Biocomputing 2024

where α0 is the intercept and α1 ̸= 0 ∈ Rn are the slopes. That is, Y is conditionally Bernoulli,
PX→Y (Y = y|X = x) = ϕyx (1 − ϕx )1−y(4)with the success probability depending on X through a sigmoid transformation,
ϕx =

1
.
1 + e−α0 −x⊤ α

To specifiy PY →X (X | Y ), we employ a multivariate linear regression model,
X = β0 + Y β1 + ϵ

where β 0 ∈ Rn are the intercepts, β 1 ̸= 0 ∈ Rn are the slopes, and ϵ ∈ Rn are Gaussian errors
with mean zero and covariance Ω. The multivariate linear regression model above implies,
PY →X (X = x|Y = y) = N (x | θ y , Ω)(5)with θy = β 0 + yβ 1 .
Putting (1)-(5) together, we have
PX→Y (X = x, Y = y) = ϕyx (1 − ϕx )1−y N (x | µ, Σ)
PY →X (X = x, Y = y) = N (x | θ y , Ω)ρy (1 − ρ)1−y(6)

2.3. Causal Identifiability
Since we only have access to observational data, the two competing causal hypotheses may
not be identifiable, i.e., PX→Y (X = x, Y = y) = PY →X (X = x, Y = y) for all x and y . For
example, if both X and Y are Gaussian, they are not identifiable. Consequently, even with an
infinite amount of data, one cannot tell these two causal models apart – clearly an undesirable
feature. Fortunately, we will show, both theoretically and empirically, that the proposed modelis identifiable.
Definition 1 (Causal Identifiability). We say H0 and H1 are identifiable if one cannot findany values of {α0 , α, µ, Σ}b and {β 0 , β 1 , Ω, ρ}c such that PX→Y (X = x, Y = y) = PY →X (X =x, Y = y) for all x and y .
Under the causal sufficiency assumption (i.e., there is no unmeasured confounder) commonly adopted in the literature,2,4,11,13,15,22,23 we have the following identifiability theorem.
Theorem 1 (Causal Identifiability). Assuming causal sufficiency, the causal hypotheses
H0 and H1 are identifiable under model (6).
Proof. We will show by contradiction. Suppose,
PX→Y (X = x, Y = y | α0 , α, µ, Σ) = PY →X (X = x, Y = y | β 0 , β 1 , Ω, ρ)b The parameters of P
c The parameters of P

X→Y (X = x, Y = y)

Y →X (X = x, Y = y)

204

(7)

Pacific Symposium on Biocomputing 2024

for all x ∈ Rn and y ∈ {0, 1}. Summing up both sides of (7) over y from 0 to 1, we have
1
X

PX→Y (Y = y | X = x, α0 , α)PX→Y (X = x | µ, Σ)

y=0

=

1
X

(8)
PY →X (X = x | Y = y, β 0 , β 1 , Ω)PY →X (Y = y | ρ)

y=0

The left-hand side of (8) is given by
1
X

PX→Y (Y = y | X = x, α0 , α)PX→Y (X = x | µ, Σ)

y=0= PX→Y (X = x | µ, Σ)

1
X

PX→Y (Y = y | X = x, α0 , α)y=0

= PX→Y (X = x | µ, Σ)(9)= N (x | µ, Σ),

where the second equality is due to the law of total probability.
The right-hand side of (8) is given by
1
X

PY →X (X = x | Y = y, β 0 , β 1 , Ω)PY →X (Y = y | ρ)y=0

= ρ · PY →X (X = x | Y = 1, β 0 , β 1 , Ω) + (1 − ρ) · PY →X (X = x | Y = 0, β 0 , β 1 , Ω)
= ρN (x | β0 + β1 , Ω) + (1 − ρ)N (x | β0 , Ω).

(10)Note that (9) is a Gaussian distribution whereas (10) is a mixture of Gaussian distribution.
Therefore, for them to be equivalent, we must have ρ = 0, ρ = 1, or β 1 = 0, which aredegenerated cases (i.e., either Y is deterministically 0 or 1, or X and Y are independent).
Although our theorem relies on the causal sufficiency assumption, the experiments in Section 3.1.3 empirically show that the proposed method is relatively robust to the presence of
unmeasured confounders.
2.4. Estimation
Theorem 1 establishes a property of the probability model and therefore is a population-level
result. It implies that for a large enough sample size, one can correctly identify the correctcausal hypothesis even with observational data alone. For a finite sample, our identifiability
result paves the way for a simple, yet useful, causal discovery algorithm based on the maximum
likelihood estimation (MLE). We aim to determine whether X causes Y or vice versa by
quantifying the respective likelihoods. Therefore, when provided with a dataset of N subjects,(x1 , y1 ), . . . , (xN , yN ), we conclude H0 : X → Y if

205

Pacific Symposium on Biocomputing 2024

max

α0 ,α,µ,Σ

N
Y

PX→Y (X = xi , Y = yi | α0 , α, µ, Σ) >

i=1

max

β 0 ,β 1 ,Ω,ρN
Y

PY →X (X = xi , Y = yi | β 0 , β 1 , Ω, ρ),

i=1

and H1 : Y → X otherwise. Note that the two competing hypotheses have the same model
complexity (i.e., the same number of parameters) and hence a model complexity penalty,
which is typically needed for model selection, is not necessary here. The factorized form of
the proposed model (6) allows us to separately find the MLE of each of its four components
using existing standard techniques.
However, we note that in our motivating application, X is high-dimensional (n = 1,440).
For better statistical and computational efficiency, we choose to reduce its dimensionality
before finding the MLE. Specifically, the functional principal component analysis (FPCA) isused, which can reduce the functional data into a few uncorrelated functional principal components (FPCs) that explain the most variation among all the functional bases. We decompose
the covariance function of a stochastic process X(·) as,
Cov(X(s), X(t)) =∞X

λk ψk (s)ψk (t),k=1where λk ’s are the nonnegative eigenvalues in descending order and ψk (·)’s are the corresponding orthogonal eigenfunctions. By the Karhunen-Loève theorem,
X(t) = µ(t) +∞X

Zk ψk (t),k=1∞where µ(t) = E[X(t)], {ψk (t)}∞k=1 is referred to as the FPCs, and {Zik (t)}k=1 denotes thecorresponding FPC scores. In practice, we would choose the first K ≪ n FPC scores Z =(Z1 , . . . , ZK )⊤ that explain 99% variance and replace X by Z in the proposed model whenfinding the MLEs.
Finally, to assess the uncertainty of our approach, we use the bootstrap3 technique inour real data application. We first generate B bootstrap samples by resampling subjects
with replacement. Each bootstrap sample has the same size as the original dataset. Then weapply our method to each bootstrap sample and record our choice between H0 and H1 . The
proportion of times that we choose H0 or H1 reflects our confidence toward each hypothesis.

3. Experiments
We first tested our model through various simulation scenarios on synthetic data where there is
known ground truth. After confirming its effectiveness, we then applied our method to a realworld mobility-activity dataset, demonstrating its practical capability in generating plausible
causal hypotheses.
3.1. Simulations
To assess the efficacy of the proposed model, we performed simulations on three different
synthetic datasets including one with unmeasured confounders. Each simulation was repeated
500 times, measuring the accuracy by the frequency at which we correctly identified the true206

Pacific Symposium on Biocomputing 2024

hypothesis. By considering varying sample sizes N between 50 and 200, we investigated theasymptotic behavior of our method. Furthermore, we examined the performance of the model
under various signal strengths. For ease of exposition, δ i always denote the standard Gaussian
white noises hereafter, i.e., δ i ∼ N (0, I) where I is an identity matrix.
3.1.1. Case 1: True Direction X → Y
For each subject i = 1, . . . , N , the functional data xi were created by first sampling their meanm from a centered Gaussian process at n = 30 evenly spaced time points,m ∼ GP(0, K)

with the powered exponential covariance function,
K(t, s) = exp{−|t − s|κ },of which the power κ = 1.9, and then setting
xi = m + δ i .

We performed the FPCA24 on x1 , . . . , xN using the R package fdapace, and retained first K
FPCs that explained 99% variance. We denote the standardized FPC scores by z 1 , . . . , z N .
To create the causal dependency of yi on xi through z i , we generated yi from a probit
regression,(yi =1 if yi∗ > 0,
0 if yi∗ ≤ 0

whereyi∗ = 0.5 + z ⊤i γ + ϵiwith ϵi ∼ N (0, 1). Here γ = (γ1 , . . . , γK )⊤ is the direct causal effect (signal), which will bevaried at three levels: weak (γk = ±1), moderate (γk = ±1.5), and strong (γk = ±3).
The simulation results are reported in Table 1, showing an expected trend: the stronger thesignal is, the more accurately the true causal direction can be discerned. Also, the accuracy approaches 100% as the sample size increases for the moderate and strong signal cases. Even withthe weak signal, the accuracy was still good, around 90%. Note that for a non-identifiability
model, the expected accuracy is 50%.
3.1.2. Case 2: True Direction Y → X
Exploring the reverse causal direction, we first generated the binary cause variable yi froma Bernoulli distribution with a success probability of 0.5. Then we generated the functional
effect variable,
xi = myi + δ i ,where my ∼ GP(0, Ky ) for y = 0, 1 with the powered exponential covariance function Ky ofwhich the power κ depends on y . Specifically, κ = 1.9 if y = 1, and κ = 0.3 (strong signal), 1.1(moderate signal), or 1.7 (weak signal) if y = 0.

207

Pacific Symposium on Biocomputing 2024

Table 1: Simulations. Accuracy of the proposed model in determining true causal directions
in synthetic datasets over 500 simulations.
Case

Confounder

Sample size

Signal
50

100

150

200

92.8%
92.8%
93.8%

88.8%
97.2%
98.6%

90%
97.6%
99%

87.8%
99.6%
99.8%

Functional

94.6%

98.2%

99.2%

99.2%

Binary

92.4%

99.2%

100%

100%

82.2%
83.4%
97.8%

89%
90.2%
99%

89.4%
96%
99.4%

93.6%
97.2%
99.8%

Functional

39.8%

61.6%

75%

83.8%

Binary

63.6%

77%

85.2%

85.2%

weak
moderatestrong

None
X→Y

weak
moderatestrongNone
Y →X

As anticipated, our simulation results (Table 1) show that as the signal or the sample size
increases, the accuracy approaches 100%.
3.1.3. Case 3: Hidden Confounders
The Simulation Cases 1&2 above demonstrate the validity of Theorem 1, i.e., causal directions
can be identified even with observational data alone. We now empirically assess the robustnessof the proposed method with respect to the violation of the causal sufficiency assumption, i.e.,
we test whether our method can still identify the correct causal direction in the presence of
unmeasured confounders.
Our methodology hinges on determining the causality between two distinct types of variables, binary scalar and continuous functional. Thus, accordingly, we considered that theunobserved confounders, generically denoted by C , are also either binary scalar or continuousfunctional. Consequently, we investigated four separate scenarios depicted in Fig. 1. We generated data from these four causal graphs and hid C from our method (i.e., only took X andY as the inputs of our algorithm). As before, we recorded the frequency at which we correctly
identified the causal direction between X and Y .
In Fig. 1 (a)&(b) where the confounder is binary, we generated the confounder ci froma Bernoulli distribution with success probability 0.5. In Fig. 1 (a), the mean mci of xi was
generated from a conditional Gaussian process mc ∼ GP(0, Kc ) with the powered exponential
covariance function Kc of which the power κ depends on c. Specifically, κ = 1.9 if c = 1 andκ = 1.5 if c = 0. Then as before, we set xi = mci + δ i . Finally, we generated yi from a probit
regression model, yi = 1 if yi∗ > 0 and yi = 0 otherwise, whereyi∗ = 0.5 + 3 · z ⊤i 1K + 3 · ci + ϵi

208

Pacific Symposium on Biocomputing 2024

C

X

C

Y

(a)

X

C

Y

X

(b)

Y

(c)
Binary

C

X

Y

(d)

Functional

Fig. 1: Four confounding scenarios under consideration in Simulation Case 3.
with ϵi ∼ N (0, 1), z i ’s are the FPC scores of xi , and 1K = (1, . . . , 1)⊤ is a vector of ones withlength K .
In Fig. 1 (b), we generated yi from a Bernoulli distribution with the success probability
ρi dependent on ci . More precisely, ρi = 0.9 if ci = 1 and ρi = 0.1 if ci = 0. Subsequently, themean mci ,yi of xi was generated from a conditional Gaussian process mc,y ∼ GP(0, Kc,y ) withthe powered exponential covariance function Kc,y of which the power κ depends on both c andy . To be specific, κ = 1.9 if c = 1 and y = 1, κ = 0.5 if c = 1 and y = 0, κ = 1.0 if c = 0 and
y = 1, and κ = 1.7 if c = 0 and y = 0. Finally, we set xi = mci ,yi + δ i .
In Fig. 1 (c)&(d), the functional confounder ci was generated in the same way as xi inCase 1 with κ = 1.5. Next, we performed the FPCA on c1 , . . . , cN and retained the first J
FPCs that explained 99% variance. We denote the standardized FPC scores by d1 , . . . , dN . In
Fig. 1 (c), m was first generated from a centered Gaussian process m ∼ GP(0, K) with κ = 1.9.
Then the dependence on the confounder was introduced by setting
xi = 0.5 + 5 · m + 5 · ci + δ i .

Finally, we generated yi from a probit regression model, yi = 1 if yi∗ > 0 and yi = 0 otherwise,where⊤yi∗ = 0.5 + 5 · z ⊤i 1K + 5 · di 1J + ϵiwith ϵi ∼ N (0, 1) and z i ’s being the first K FPC scores of xi .
In Fig. 1 (d), we first generated yi from a probit regression model, yi = 1 if yi∗ > 0 andyi = 0 otherwise, whereyi∗ = 0.5 + 3 · d⊤i 1 J + ϵiwith ϵi ∼ N (0, 1). Then we generated mean processes myi from a conditional Gaussian process
my ∼ GP(0, Ky ). In this setting, the power κ of the powered exponential covariance function
Ky depended on y : κ = 1.9 if y = 1, and κ = 0.5 if y = 0. To introduce the influence of theconfounder, we defined
x i = m yi + 3 · c i + δ i .

The results from these four confounding scenarios (Table 1) demonstrate the robustnessof the proposed method. Particularly, as the sample size increased, our method achieved

209

Pacific Symposium on Biocomputing 2024

increasingly better accuracy and was significantly better than a random guess for large sample
sizes.
3.2. Real Data
Next, we applied the proposed methodology to the data collected by the NHANES. This extensive study, conducted by the Centers for Disease Control, gathered a wide range of healthand nutritional information about the U.S. population, including sociodemographic characteristics and various health conditions. To demonstrate the utility of the proposed method,
we are particularly interested in two variables, physical activities X captured via hip-attachedaccelerometers and an indicator variable of mobility issues Y derived from self-reported household interview data. Given the logical assumption of Y → X in this scenario, we aim to verify
if our method can correctly identify this causal direction, primarily seeking to validate theeffectiveness of our method in accurately determining causation from a known truth.
3.2.1. Data Preprocessing
Utilizing the NHANES dataset, we accessed activity data from hip-worn accelerometers during the 2003–2004 and 2005–2006 study waves. The magnitude of acceleration (movement“intensity”) was captured using the ActiGraph AM-7164, delivering an objective measure ofphysical activity and bypassing the inconsistencies of self-reported data. Participants were
instructed to wear the device for seven consecutive days, excluding swimming and bathing
periods. The raw data were segmented into one-minute intervals or “epochs” with intensity
readings accumulated per epoch and saved in long format (each row is a subject-minute).
The well-formatted data are contained in the R package rnhanesdata.7,8 Following the
preprocessing procedure in their paper,8 we included individuals aged 50 to 85 and omitted
non-compliant individuals who have excessive missing accelerometer data, leaving us with
N = 3, 198 subjects.
The activity data for each individual were aggregated over the 7-day period and transformed via log(1 + x). This dataset is organized in a 7N × 1440 matrix, with one row designatedfor each subject-day across all NHANES waves, where 7 denotes the days each subject wore
the accelerometer, and 1440 corresponds to the total number of minutes in a day.
The presence of any mobility issues was represented as a binary variable, categorized aseither “No difficulty” or “Any difficulty,” based on responses from the Physical Functioning
questionnaire. Individuals were classified under “Any difficulty” if they reported challenges
in climbing 10 stairs, walking a quarter mile, abstained from these activities, or required
special walking equipment. Overall, there are 32.4% subjects in the sample who experienceany mobility of movement problem.
3.2.2. Results
We generated B = 100 bootstrap samples and successfully identified the correct causal direction across all samples from comparing the maximized likelihoods of Y → X and X → Y : themobility issue Y unambiguously impacts an individual’s level of physical activity X with high210

Pacific Symposium on Biocomputing 2024

0.020

density

0.015

0.010

0.005

0.000
80

120

160

200

Log−Likelihood Difference in Competing Hypotheses

Fig. 2: Real data. Histogram depicting the maximized log-likelihood differences between thecompeting hypotheses Y → X and X → Y .
confidence. As depicted in Fig. 2, the histogram illustrates the difference in the maximized
log-likelihoods between these two competing hypotheses (the former minus the latter), which
is noticeably bounded away from zero, meaning that Y → X is far more likely than X → Y ,which matches the presumed truth.
4. Discussion
In this paper, we have presented a new causal model for generating bivariate causal hypotheseswith a continuous functional variable (e.g., physical activities) and a binary scalar variable(e.g., mobility issue indicator) in an exploratory fashion, which can provide insights as towhich variable is more likely the cause. We theoretically proved that the underlying causeeffect relationship is identifiable with purely observational data under the causal sufficiency
assumption. Empirically, we used a likelihood-based inference procedure and demonstrated theutility of the proposed method both under and beyond the causal sufficiency setting throughsimulation studies and a real-world wearable device application.
There are several areas where this paper could be strengthened and extended. First, our
NHANES application has focused on physical activities and mobility issue because of theirclear causal relationship. Having demonstrated it is possible to identify their causal relationship, we plan to analyze other variables in the data to generate causal hypotheses in an
exploratory manner, which is an intended use of the proposed method.
Second, our identifiability theory operates under the assumption that there are no unmeasured confounders. Even though our empirical investigations have indicated a degree ofrobustness to the presence of confounders, a theoretical exploration of identifiability within
this context would be interesting and particularly relevant in observational studies where thepresence of unmeasured confounders is common.
Third, we have focused on the bivariate case and hence an extension to multivariate

211

Pacific Symposium on Biocomputing 2024

cases and leveraging additional publicly available datasets can considerably broaden themethod’s applicability. For example, the brain electroencephalogram dataset1 comprises electroencephalogram signals collected over various trials with distinct stimuli for two groups alcoholics and controls. By viewing the electroencephalogram signals as multivariate functionaldata, a recent paper23 attempts to discern the causal relationships among these functions. The
multivariate extension of our method could potentially enrich this research by providing additional insights into the causal relationships modified by the experimental groups by treating
the group as a binary variable. Moreover, it should be relatively straightforward to extend our
method to incorporate multiple categorical scalar variables.
Finally, a Bayesian inference approach could be adopted especially for multivariate cases
where efficient searching strategies in the causal graph space are required. A Bayesian approach
would make it easier to make finite-sample inferences with natural uncertainty quantification
for complex causal graphs.
Acknowledgment
Ni’s research was partially supported by NIH 1R01GM148974-01 and NSF DMS-2112943.
References
1. H. Begleiter. EEG Database. UCI Machine Learning Repository, 1999. DOI:
https://doi.org/10.24432/C5TS3D.
2. J. Choi, R. Chapkin, and Y. Ni. Bayesian causal structural learning with zero-inflated poisson
Bayesian networks. Advances in neural information processing systems, 33:5887–5897, 2020.
3. B. Efron. Bootstrap methods: another look at the jackknife. In Breakthroughs in statistics:
Methodology and distribution, pages 569–593. Springer, 1992.
4. P. Hoyer, D. Janzing, J. M. Mooij, J. Peters, and B. Schölkopf. Nonlinear causal discovery withadditive noise models. Advances in neural information processing systems, 21, 2008.
5. D. Janzing, J. Mooij, K. Zhang, J. Lemeire, J. Zscheischler, P. Daniušis, B. Steudel, and
B. Schölkopf. Information-geometric approach to inferring causal directions. Artificial Intelligence, 182:1–31, 2012.
6. K.-Y. Lee and L. Li. Functional structural equation model. Journal of the Royal StatisticalSociety Series B: Statistical Methodology, 84(2):600–629, 2022.
7. A. Leroux. rnhanesdata: NHANES Accelerometry Data Pipeline, 2023. R package version 1.02.
8. A. Leroux, J. Di, E. Smirnova, E. J. Mcguffey, Q. Cao, E. Bayatmokhtari, L. Tabacu, V. Zipunnikov, J. K. Urbanek, and C. Crainiceanu. Organizing and analyzing the activity data in nhanes.
Statistics in biosciences, 11:262–287, 2019.
9. S. Monti and G. F. Cooper. A multivariate discretization method for learning bayesian networks
from mixed data. arXiv preprint arXiv:1301.7403, 2013.
10. J. S. Morris. Functional regression. Annual Review of Statistics and Its Application, 2:321–359,
2015.
11. Y. Ni. Bivariate causal discovery for categorical data via classification with optimal label permutation. Advances in Neural Information Processing Systems, 35:10837–10848, 2022.
12. J. Pearl et al. Models, reasoning and inference. Cambridge, UK: CambridgeUniversityPress,
19(2):3, 2000.
13. S. Shimizu, P. O. Hoyer, A. Hyvärinen, A. Kerminen, and M. Jordan. A linear non-gaussian
acyclic model for causal discovery. Journal of Machine Learning Research, 7(10), 2006.

212

Pacific Symposium on Biocomputing 2024

14. P. Spirtes, C. Glymour, and R. Scheines. Causation, Prediction, and Search, 2nd Edition, volume 1 of MIT Press Books. The MIT Press, December 2001.
15. P. Spirtes, C. N. Glymour, and R. Scheines. Causation, prediction, and search. MIT press, 2000.
16. O. Stegle, D. Janzing, K. Zhang, J. M. Mooij, and B. Schölkopf. Probabilistic latent variable
models for distinguishing between cause and effect. Advances in neural information processing
systems, 23, 2010.
17. C. D. Tekwe, R. S. Zoh, M. Yang, R. J. Carroll, G. Honvoh, D. B. Allison, M. Benden, andL. Xue. Instrumental variable approach to estimating the scalar-on-function regression modelwith measurement error with application to energy expenditure assessment in childhood obesity.
Statistics in medicine, 38(20):3764–3781, 2019.
18. J.-L. Wang, J.-M. Chiou, and H.-G. Müller. Functional data analysis. Annual Review of Statistics
and Its Application, 3(1):257–295, 2016.
19. W. Wenjuan, F. Lu, and L. Chunchen. Mixed causal structure discovery with application toprescriptive pricing. In Proceedings of the 27th International Joint Conference on Artificial Intelligence, pages 5126–5134, 2018.
20. M. Yamayoshi, J. Tsuchida, and H. Yadohisa. An estimation of causal structure based on latent
lingam for mixed data. Behaviormetrika, 47:105–121, 2020.
21. Y. Zeng, S. Shimizu, H. Matsui, and F. Sun. Causal discovery for linear mixed data. In
B. Schölkopf, C. Uhler, and K. Zhang, editors, Proceedings of the First Conference on Causal
Learning and Reasoning, volume 177 of Proceedings of Machine Learning Research, pages 994–
1009. PMLR, 11–13 Apr 2022.
22. K. Zhang and A. Hyvärinen. On the identifiability of the post-nonlinear causal model. 647–655,
2009.
23. F. Zhou, K. He, K. Wang, Y. Xu, and Y. Ni. Functional bayesian networks for discovering
causality from multivariate functional data. arXiv preprint arXiv:2210.12832, 2022.
24. Y. Zhou, S. Bhattacharjee, C. Carroll, Y. Chen, X. Dai, J. Fan, A. Gajardo, P. Z. Hadjipantelis,
K. Han, H. Ji, C. Zhu, H.-G. Müller, and J.-L. Wang. fdapace: Functional Data Analysis andEmpirical Dynamics, 2022. R package version 0.5.9.

213