Close
About
FAQ
Home
Collections
Login
USC Login
Register
0
Selected
Invert selection
Deselect all
Deselect all
Click here to refresh results
Click here to refresh results
USC
/
Digital Library
/
University of Southern California Dissertations and Theses
/
Bayesian models for a respiratory biomarker with an underlying deterministic model in population research
(USC Thesis Other)
Bayesian models for a respiratory biomarker with an underlying deterministic model in population research
PDF
Download
Share
Open document
Flip pages
Copy asset link
Request this asset
Request accessible transcript
Transcript (if available)
Content
Bayesian models for a respiratory biomarker with an underlying deterministic model in
population research
By
Jingying Weng
A Dissertation Presented to the
FACULTY OF THE USC GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
(BIOSTATISTICS)
August 2022
Copyright 2022 Jingying Weng
ii
Table of Contents
List of Tables .................................................................................................................... vii
List of Figures .................................................................................................................. viii
Abbreviations ...................................................................................................................... x
Abstract .............................................................................................................................. xi
Chapter 1. Overview ........................................................................................................... 1
1.1 Background ........................................................................................................... 1
Fractional concentration of exhaled Nitric Oxide (FeNO) ................................. 1
Multiple flow FeNO ............................................................................................ 2
Two compartment model (2CM) for FeNO ........................................................ 3
Estimating airway and alveolar nitric oxide in the 2CM .................................... 3
Bayesian methods for hierarchical nonlinear models. ........................................ 5
1.2 Dissertation overview ........................................................................................... 5
Project 1 (Chapter 2). Hierarchical Bayesian estimation of covariate effects on
airway and alveolar nitric oxide .......................................................................... 6
Project 2 (Chapter 3). Longitudinal Hierarchical Bayesian models of covariate
effects on airway and alveolar nitric oxide ......................................................... 7
iii
Project 3. (Chapter 4) Statistical modeling of multiple flow exhaled nitric
oxide in population research ............................................................................... 8
Reference .................................................................................................................... 9
Chapter 2: Hierarchical Bayesian estimation of covariate effects on airway and alveolar
nitric oxide ........................................................................................................................ 11
Abstract ..................................................................................................................... 11
Introduction ............................................................................................................... 13
Methods..................................................................................................................... 16
Deterministic, steady-state two-compartment model for FeNO ....................... 16
Estimating NO parameters from the 2CM ........................................................ 17
Estimating associations of covariates with NO parameters .............................. 17
Prior distributions.............................................................................................. 23
Calculation of posterior via MCMC: ................................................................ 23
Simulation study ............................................................................................... 24
CHS data analysis ............................................................................................. 27
Results ....................................................................................................................... 29
Simulation study: .............................................................................................. 29
CHS data analysis: ............................................................................................ 35
Discussion ................................................................................................................. 38
iv
Software .................................................................................................................... 42
References ................................................................................................................. 43
Supplement: .............................................................................................................. 46
Section 1: Convergence failure rates in the simulation study ........................... 47
Section 2: Bias, Coverage, Power and Confidence Interval Length figures for
𝛽 , by scenario .................................................................................................... 49
Section 3: Bias, Coverage, Power and Confidence Interval Length figures for
𝛼 , by scenario .................................................................................................... 64
Section 4. Additional details on CHS results .................................................... 77
For CHS participants, the estimated participant-level NO parameters were
highly correlated (Pearson’s R>0.92) across U-HB, U-NLME and TS-NLME
methods under the Caw parameterization (Supplementary Figure 4.1) and also
under the J’aw parameterization (Supplementary Figure 4.2). As shown in
Supplementary Figure 4.3, calculated logJaw estimates from U-HB (Caw
parameterization) were extremely highly correlated (Pearson’s R>0.999) with
direct logJaw estimates from U-HB (Jaw-parameterization). ............................. 78
Section 5. Simulation study results using all available estimates (not limited to
datasets for which all methods converged). ...................................................... 84
Chapter 3: Longitudinal Hierarchical Bayesian models of covariate effects on airway and
alveolar nitric oxide .......................................................................................................... 86
Abstract ..................................................................................................................... 86
v
Introduction ............................................................................................................... 88
Methods..................................................................................................................... 90
Deterministic two compartment model for FeNO ............................................ 90
Estimating NO parameters in the 2CM ............................................................. 90
Estimating associations of covariates with longitudinally assessed NO
parameters ......................................................................................................... 91
Simulation study ............................................................................................... 98
CHS data analysis ........................................................................................... 100
Results ..................................................................................................................... 102
Simulation study ............................................................................................. 102
CHS data analysis ........................................................................................... 107
Discussion ............................................................................................................... 110
References ............................................................................................................... 114
Supplementary materials ......................................................................................... 118
Chapter 4: Statistical modeling of multiple flow exhaled nitric
oxide in population research ........................................................................................... 148
Abstract: .................................................................................................................. 148
Introduction ............................................................................................................. 149
Methods................................................................................................................... 152
vi
Multiple flow FeNO data ................................................................................ 152
Two compartment model (2CM) .................................................................... 152
Estimating NO parameters for a single participant ......................................... 153
Extensions to longitudinal study designs ........................................................ 160
Strengths and weaknesses of the six estimation methods ............................... 161
R package ........................................................................................................ 162
Results ..................................................................................................................... 165
Discussion ............................................................................................................... 170
Appendix ................................................................................................................. 173
Appendix 1: Example R code for population multiple flow FeNO studies .... 173
Appendix 2: Customized U-NLME and U-HB: ............................................. 178
References ............................................................................................................... 183
Chapter 5. Summary and future directions ..................................................................... 191
Summary ................................................................................................................. 191
Two stage versus unified approaches. ............................................................ 196
Future work and applications .......................................................................... 197
vii
List of Tables
Table 1: Seven scenarios considered in the cross-sectional data simulation study ...26
Table 2: Analysis of CHS data...................................................................................37
Table 3: Parameter values used to generate simulation study data ............................98
Table 4: Four scenarios considered in the longitudinal data simulation study. .........99
Table 5: Highlights of simulation studies. ...............................................................161
viii
List of Figures
Figure 1: Examples of real-time FeNO and flow .........................................................2
Figure 2: Two compartment model of FeNO ..............................................................3
Figure 3: Hierarchical model structure relating FeNO measurements at multiple flow
rates to NO parameters that are a function of a potential determinant X (e.g., air
pollution). ...........................................................................................................21
Figure 4: Relative bias (a), coverage (b), power, (c) and CI length (d) of the selected
estimation methods from Scenario 1..................................................................31
Figure 5: Type I error
of the selected estimation methods .........................................33
Figure 6: Analysis of CHS data: estimated associations of traffic-related air
pollution with CA, logCaw, and logDaw (𝛽 and 95% CI). ...................................35
Figure 7: Diagram of the hierarchical model structure relating Longitudinal FeNO
measurements at multiple flow rates to NO parameters that are a function of a
potential determinant X (e.g., air pollution) .....................................................94
Figure 8: Comparison of method performance. .......................................................102
Figure 9: Power and Type I errors for the 6 methods ..............................................105
Figure 10: Estimated associations between NO parameters and standardized height
in the CHS using 6 methods.............................................................................108
ix
Figure 11: Multiple flow FeNO data used to estimate NO parameters from the two-
compartment model for a single participant. ..................................................152
Figure 12: Steady state two-compartment model and key NO parameters. .........153
Figure 13: Conceptual diagram of two-stage estimation analysis in multiple flow
exhaled nitric oxide population research. ........................................................156
Figure 14: Conceptual diagram of unified estimation analysis in multiple flow
exhaled nitric oxide population research. ........................................................159
Figure 15: Estimated differences in NO parameters comparing asthma/allergy
groups to a reference (healthy) group of children ............................................167
x
Abbreviations
NO parameters: CA, logCaw, logDaw
CA: concentration of NO in the alveolar region, (ppb)
Caw: concentration of NO in the airway tissue, (ppb)
Daw: airway tissue diffusion capacity, pL· s
-1
· ppb
-1
J’aw: maximum flux of NO in the airway
FeNO: fractional concentration of exhaled nitric oxide, ppb
NLME: nonlinear mixed effect model
HMA: Högman and Merilӓinen Algorithm
xi
Abstract
This dissertation is the culmination of three projects, two statistical and one focused on
disseminating the technical results to the breath research community along with
companion software. All projects are motivated by the multiple flow FeNO data collected
in the Southern California Children’s Health Study (CHS), introduce this dataset more
fully in Chapter 2. Chapter 1 introduces the background and motivations. In Chapter 2,
we develop a novel Unified Hierarchical Bayesian (U-HB) method. We then conduct an
extensive simulation study to evaluate its properties and compare it to alternative unified
and two step approaches. Next, we then re-analyze cross-sectional data from the CHS to
estimate the association of traffic-related air pollution exposure with alveolar and airway
inflammation. In Chapter 3, we introduce a longitudinal extension of U-HB and
demonstrate that HB methods outperform alternative methods for estimating associations
of covariates with NO parameters. A manuscript summarizing this work is currently
under review (submitted March 2022). Chapter 4 aims to provide a high-level overview
of the technical work on cross-sectional and longitudinal methods in Chapters 2 and 3,
targeting researchers conducting population studies with multiple flow FeNO. We review
key statistical issues, provide an overview of relevant statistical methods while
highlighting their strengths and weaknesses, and demonstrate how to implement these
methods using R code that we have developed and made freely available on GitHub.
Using cross sectional multiple flow data from the CHS, we apply six methods to estimate
differences in NO parameters, comparing healthy children to those with asthma only,
xii
allergic rhinitis only, or both asthma and allergy. Chapter 5 summarizes the dissertation
and discuss existing problems and future development of HB methods
1
Chapter 1. Overview
1.1 Background
As early warning systems of health, biomarkers are quantitative measures informing what
is currently or has happened in an organism. Biomarkers have been widely applied in
human studies, both in clinical and epidemiological applications. Exhaled breath
biomarkers have many advantageous features, including non-invasive assessment that is
easily repeatable.[1] The fractional concentration of exhaled nitric oxide, or FeNO, is one
of the most well-studied biomarkers in exhaled breath but there remain important
challenges associated with the assessment of FeNO which can be addressed through the
development of novel statistical methods. The theme of this dissertation is the
development of statistical methods to address unmet needs in population-level research
using a cutting-edge method which partitions FeNO into its proximal and distal sources.
Fractional concentration of exhaled Nitric Oxide (FeNO)
In 1991, it was discovered that NO is produced in airway tissue and measurable in
exhaled breath.[2] Since then, FeNO has gained widespread use as a marker of airway
inflammation with applications in clinical settings (e.g., asthma[3]) and epidemiological
research (e.g., studies of inhaled environmental exposures[4, 5]). Early in the study of
FeNO, it was discovered that FeNO is inversely related to expiratory flow rate[6, 7].
2
Conventional assessment of FeNO (FeNO50) is conducted at a fixed 50 ml/s flow rate[3,
8, 9], essentially treating flow dependency as a nuisance.
Multiple flow FeNO
Measuring FeNO at multiple (rather than single)
flow rates provides additional information that
can be leveraged to gain insights about the
proximal versus distal sources of NO. At high
flow rates, FeNO primarily reflects lower
concentrations in the alveolar region while at low
flow rates, it primarily reflects higher
concentrations in the airway tissue[10, 11]. The
growing interest in multiple flow FeNO has
resulted in its inclusion in the most recent update
to the guidelines for FeNO assessment.[12]
Several mathematical models have been
developed to describe the dynamics of NO in the
lower respiratory tract, using systems of
differential equations with parameters describing
NO sources in the airway and alveolar compartments.[13]
Multiple flow FeNO is collected using real-time measurement, as shown in Figure 1. The
participant inhales and exhales several times in a closed mask or spirometer using
pollution free air to clear the NO in the airway and lung, and then takes a deep breath to
exhale slowly and constantly to collect the real-time FeNO concentration and flow rate.
Figure 1: Examples of real-time
FeNO and flow profiles from four
maneuvers at different target flow
rates, collected as part of a multiple
flow FeNO protocol.
3
As shown in Figure 2, time-series of FeNO concentrations from an expiratory maneuver
start at zero and increase rapidly toward a steady state. The maneuver is summarized by
averaging the NO concentration and flow rate during a plateau within the steady state for
several seconds.
Two compartment model (2CM) for FeNO
The statistical models in this dissertation are based on
the simple steady-state two-compartment model
(henceforth referred to as the 2CM) of FeNO, which
assumes a cylindrically shaped airway compartment
and an expansile alveolar compartment (see Figure
2). Under the 2CM, FeNO measured at the mouth can
be related to parameters quantifying airway and
alveolar NO sources as follows (CA: the concentration of NO in the alveolar region,
(ppb); Caw: the concentration of NO in the airway tissue, (ppb); and Daw: the airway tissue
diffusion capacity, pL· s
-1
· ppb
-1
):
𝐹𝑒𝑁𝑂 = 𝐶 𝑎𝑤
+ ( 𝐶 𝐴 − 𝐶 𝑎𝑤
)× 𝑒 −𝐷 𝑎𝑤
flow ⁄
These so-called NO parameters are of increasing interest in clinical research.[14]
Estimating airway and alveolar nitric oxide in the 2CM
Figure 2: Two compartment model of FeNO
4
To estimate NO parameters (CA, Caw, Daw) of the 2CM model for a single participant at a
single visit, researchers have relied on linear regression-based methods—under various
linearizing assumptions. Examples include the Taylor expansion based first order linear
approximation methods[15], linT and linP[13]; second order quadratic approximation
methods, quadT and quadP[16], and a third order approximation method, the Hö gman
and Merilӓinen algorithm (HMA). Other approaches include nonlinear regression[17] and
an extension in which a natural log transformation is applied to both sides,
acknowledging the increased variation in residuals that occurs as flow rate increases.[16]
In biomedical and epidemiological research studies, multiple flow FeNO is typically
measured in multiple participants, potentially at multiple visits. Using data from this
study population, it is often of interest to estimate how NO parameters are associated
with disease, treatment, or environmental exposures. Researchers have typically used ad
hoc two stage methods for these studies. In these two stage (TS) methods, NO parameters
are estimated separately for each participant in Stage 1 and then treated as known
outcomes in Stage 2. But previous studies found that TS methods can suffer from issues
related to small sample sizes in Stage I (e.g., < 10 FeNO maneuvers are used to estimate
3 NO parameters for each participant) which can result in Stage I models failing to
converge or producing unreasonable values for NO parameters. Participants whose data
poses problems in Stage I are often excluded from Stage II, potentially biasing Stage II
results. An appealing alternative is unified models in which we simultaneously estimate
NO parameters and their associations with covariates. This dissertation introduces a new
unified approach in which we estimate NO parameters and their associations with
5
covariates together using a Hierarchical Bayesian approach. The new unified Hierarchical
Bayesian approach is called U-HB. In cross-sectional and longitudinal settings, U-HB is
compared to two stage approaches using extensive simulation studies and multiple flow
FeNO data collected from the southern California Children’s Health Study.
Bayesian methods for hierarchical nonlinear models.
The Bayesian approach allows for full characterization of uncertainty in unknown
parameters in the model. It has been applied widely in many fields, including
pharmacokinetics modeling, which is an application closely related to the multiple flow
FeNO problem, in which deterministic models are linked with hierarchical statistical
models to estimate biologically meaningful parameters (and their associations with
covariates). In the U-HB model, FeNO is modeled with a two-level hierarchy for cross-
sectional data: maneuver and participant levels and with a three-level hierarchy for
longitudinal data: maneuver, visit, and participant levels. In each level, we define a set of
prior distributions for unknown parameters and their hyperparameters and any constraints
of parameter spaces (e.g., positivity). Due to the normal conjugacy of the distributions in
our particular hierarchical structure, we were able to apply Gibbs sampling to accelerate
the MCMC fitting process using the JAGS program in R. In particular, we used the
R2jags package.
1.2 Dissertation overview
6
This dissertation is the culmination of three projects, two statistical and one focused on
disseminating the technical results to the breath research community along with
companion software. All projects are motivated by the multiple flow FeNO data collected
in the Southern California Children’s Health Study (CHS). The CHS is a population-
based prospective cohort study of schoolchildren originally designed to study the long-
term respiratory health effects of exposure to air pollution. The CHS includes one of the
largest and longest studies of FeNO in children to date. We introduce this dataset more
fully in Chapter 2.
Project 1 (Chapter 2). Hierarchical Bayesian estimation of covariate effects on airway
and alveolar nitric oxide
In cross-sectional population studies with multiple flow FeNO data, traditionally NO
parameters have been estimated separately for each participant (Stage I) and then related
to covariates (Stage II). Our goal was to develop an estimation method to simultaneously
estimate NO parameters and their association with covariates. We developed a novel
Unified Hierarchical Bayesian (U-HB) method, that we fitted using Gibbs sampling. An
extensive simulation study was conducted to evaluate the properties of U-HB and
compare it to alternative unified and two step approaches. We then re-analyzed cross-
sectional data from the CHS to estimate the association of traffic-related air pollution
exposure with alveolar and airway inflammation. This dissertation includes a preprint
version of work, which has already been published:
7
Weng, J., Molshatzki, N., Marjoram, P., Gauderman, W.J., Gilliland, F.D. and Eckel,
S.P., 2021. Hierarchical Bayesian estimation of covariate effects on airway and alveolar
nitric oxide. Scientific Reports, 11(1), pp.1-12.
Project 2 (Chapter 3). Longitudinal Hierarchical Bayesian models of covariate effects on
airway and alveolar nitric oxide
In studies of chronic respiratory disease, such as asthma, longitudinal studies of within-
participant changes in biomarkers can be particularly informative. This is also the case
with NO parameters as assessed by multiple flow FeNO. While longitudinal population
studies of multiple flow FeNO are of interest, the statistical properties of ad hoc two stage
methods used to date by researchers have not been evaluated and unified methods have
not been developed. Our goal was to extend the cross-sectional U-HB model to studies
with longitudinal assessments of multiple flow FeNO. We developed a longitudinal
version L-U-HB and its two stage alternative L-TS-HB. Longitudinal unified and two
stage methods were compared in extensive simulation studies and in application to the
CHS evaluating associations of height with NO parameters. In the simulation studies, we
demonstrated that HB methods outperformed alternative methods for estimating
associations of covariates with NO parameters. Using CHS data, L-U-HB estimated
positive associations of height with airway and alveolar NO concentrations. A manuscript
summarizing this work is currently under review (submitted March 2022).
8
Project 3. (Chapter 4) Statistical modeling of multiple flow exhaled nitric
oxide in population research
This chapter aims to provide a high-level overview of the technical work on cross-
sectional and longitudinal methods in Chapters 2 and 3, targeting researchers conducting
population studies with multiple flow FeNO. We review key statistical issues, provide an
overview of relevant statistical methods while highlighting their strengths and
weaknesses, and demonstrate how to implement these methods using R code that we have
developed and made freely available on GitHub. Using cross sectional multiple flow data
from the CHS, we apply six methods to estimate differences in NO parameters,
comparing healthy children to those with asthma only, allergic rhinitis only, or both
asthma and allergy.
In summary, this dissertation proposes novel statistical methods for population studies
using multiple flow FeNO, which were evaluated in simulation studies and applied to
CHS data. We have developed publicly available R code with the goal of facilitating the
adoption of the methods by researchers working with multiple flow FeNO data.
9
Reference
1. Braun PX, Gmachl CF, Dweik RA. Bridging the Collaborative Gap: Realizing the
Clinical Potential of Breath Analysis for Disease Diagnosis and Monitoring–Tutorial.
Sensors Journal, IEEE. 2012;12(11):3258-3270.
2. Gustafsson L, Leone A, Persson M, Wiklund N, Moncada S. Endogenous nitric
oxide is present in the exhaled air of rabbits, guinea pigs and humans. Biochem Biophys
Res Commun. 1991;181(2):852-857.
3. Dweik RA, Boggs PB, Erzurum SC, et al. An official ATS clinical practice
guideline: interpretation of exhaled nitric oxide levels (FENO) for clinical applications.
Am J Respir Crit Care Med. 2012;184(5).
4. La Grutta S, Ferrante G, Malizia V, Cibella F, Viegi G. Environmental effects on
fractional exhaled nitric oxide in allergic children. J Allergy (Cairo). 2012;2012:916926.
5. Scarpa MC, Kulkarni N, Maestrelli P. The role of non-invasive biomarkers in
detecting acute respiratory effects of traffic-related air pollution. Clin Exp Allergy.
2014;44(9):1100-1118.
6. Hogman M, Stromberg S, Schedin U, Frostell C, Hedenstierna G, Gustafsson LE.
Nitric oxide from the human respiratory tract efficiently quantified by standardized single
breath measurements. Acta Physiol Scand. 1997;159(4):345-346.
7. Silkoff PE, McClean PA, Slutsky AS, et al. Marked flow-dependence of exhaled
nitric oxide using a new technique to exclude nasal nitric oxide. Am J Respir Crit Care
Med. 1997;155(1):260-267.
8. ATS. Recommendations for standardized procedures for the on-line and off-line
measurement of exhaled lower respiratory nitric oxide and nasal nitric oxide in adults and
children-1999. This official statement of the American Thoracic Society was adopted by
the ATS Board of Directors, July 1999. Am J Respir Crit Care Med. 1999;160(6):2104-
2117.
9. ATS/ERS. ATS/ERS recommendations for standardized procedures for the online
and offline measurement of exhaled lower respiratory nitric oxide and nasal nitric oxide,
2005. Am J Respir Crit Care Med. 2005;171(8):912-930.
10. George SC. How accurately should we estimate the anatomical source of exhaled
nitric oxide? J Appl Physiol. 2008;104(4):909-911.
11. Eckel SP, Salam MT. Single high flow exhaled nitric oxide is an imperfect proxy for
distal nitric oxide. Occup Environ Med. 2013;70(7):519-520.
12. Horvath I, Barnes PJ, Loukides S, et al. A European Respiratory Society technical
standard: exhaled biomarkers in lung disease. Eur Respir J. 2017;49(4).
13. George SC, Hogman M, Permutt S, Silkoff PE. Modeling pulmonary nitric oxide
exchange. J Appl Physiol. 2004;96(3):831-839.
10
14. Lehtimä ki L, Karvonen T, Hö gman M. Clinical values of nitric oxide parameters
from the respiratory system. Current Medicinal Chemistry. 2020;27(42):7189-7199.
15. SC G, M H, S P, PE S. Modeling pulmonary nitric oxide exchange. Journal of
applied physiology (Bethesda, Md : 1985). 2004;96(3).
16. Silkoff PE, Sylvester JT, Zamel N, Permutt S. Airway nitric oxide diffusion in
asthma: Role in pulmonary function and bronchial responsiveness. Am J Respir Crit Care
Med. 2000;161(4 Pt 1):1218-1228.
17. Eckel SP, Linn WS, Berhane K, et al. Estimation of parameters in the two-
compartment model for exhaled nitric oxide. PLoS ONE. 2014;9(1):e85471.
11
Chapter 2: Hierarchical Bayesian estimation of covariate
effects on airway and alveolar nitric oxide
Abstract
Exhaled breath biomarkers are an important emerging field. The fractional concentration
of exhaled nitric oxide (FeNO) is a marker of airway inflammation with clinical and
epidemiological applications (e.g., air pollution health effects studies). Systems of
differential equations describe FeNO—measured non-invasively at the mouth—as a
function of exhalation flow rate and parameters representing airway and alveolar sources
of NO in the airway. Traditionally, NO parameters have been estimated separately for
each study participant (Stage I) and then related to covariates (Stage II). Statistical
properties of these two-step approaches have not been investigated. In simulation studies,
we evaluated finite sample properties of existing two-step methods as well as a novel
Unified Hierarchical Bayesian (U-HB) model. The U-HB is a one-step estimation method
developed with the goal of properly propagating uncertainty as well as increasing power
and reducing type I error for estimating associations of covariates with NO parameters.
We demonstrated the U-HB method in an analysis of data from the southern California
Children’s Health Study relating traffic-related air pollution exposure to airway and
alveolar airway inflammation.
12
Notation:
NO parameters: CA, logCaw, logDaw
CA: concentration of NO in the alveolar region, (ppb)
Caw: concentration of NO in the airway tissue, (ppb)
Daw: airway tissue diffusion capacity, pL· s
-1
· ppb
-1
J’aw: maximum flux of NO in the airway
FeNO: fractional concentration of exhaled nitric oxide, ppb
NLME: nonlinear mixed effect model
HMA: Högman and Merilӓinen Algorithm
13
Introduction
Exhaled breath biomarkers are an important emerging field, since they can be measured
non-invasively and repeatedly.[1] The fractional concentration of exhaled nitric oxide
(FeNO) is a marker of airway inflammation with applications in clinical settings (e.g.,
asthma[3]) and epidemiological research (e.g., studies of inhaled environmental
exposures[4, 5]). The nature of the non-invasive assessment of FeNO results in
challenges that can be at least partially addressed by innovations in statistical
methodology.
FeNO concentrations are inversely related to expiratory flow rate, suggesting both an
airway and an alveolar source of nitric oxide (NO).[6, 7] High flow FeNO primarily
reflects low NO concentrations in the alveolar region while low flow FeNO primarily
reflects higher NO concentrations in the airway tissue.[10, 11] Several mathematical
models have been developed to describe the dynamics of NO in the lower respiratory
tract, using systems of differential equations with parameters describing NO sources in
the airway and alveolar compartments.[13] Conventional standardized assessment of
FeNO at a fixed 50 ml/s flow rate[3, 8, 9] treats flow dependency as a nuisance.
However, measurement of FeNO at multiple expiratory flow rates (“multiple flow
FeNO”) is a promising technique that takes advantage of the information across flow
rates to non-invasively assess airway and alveolar inflammation using estimated
parameters quantifying airway and alveolar source of NO. The growing interest in
multiple flow FeNO has resulted in its inclusion in the most recent update to the
guidelines for FeNO assessment.[12]
14
Non-invasive assessment of airway and alveolar NO has also sparked interest in relating
estimated airway and/or alveolar NO to covariates, such as respiratory diseases[18, 19]
and environmental exposures, including: ambient air pollution,[20, 21] traffic-related air
pollution [22] and pollen.[23] All of these studies have used a “Two-Stage” (TS)
approach to estimating the associations of covariates with NO parameters. Stage I
consists of estimating NO parameters for each participant (using a variety of methods,
which we will describe later). Stage II consists of relating each estimated NO parameter
(treated as a known outcome) to covariates in a standard linear regression model.
However, the statistical properties of the TS approach have not been thoroughly
investigated.
The TS approach ignores the uncertainty in the Stage I estimated NO parameters when
they are used in Stage II, thereby failing to propagate the statistical uncertainty.
Additionally, TS methods may fail to produce NO parameter estimates in Stage I for
some participants due to small sample size (e.g., in our motivating example there are ~9
observations/participant).[16] Subsequent exclusion of these subjects in Stage II may bias
the estimated association in Stage II. Another issue is that the Stage I model (i.e., the
model estimating NO parameters) is only an approximation of reality and may include
various assumptions aimed at simplification which impact parameter estimation (e.g.,
alveolar NO estimation has been shown to be sensitive to applying a first, second, or third
order approximation to an exponential function in the Stage I model.)[16] To overcome
the weakness of TS approaches, we had previously explored a nonlinear mixed effect
15
(NLME) unified model to simultaneously estimate NO parameters and their associations
with covariates. However, the method suffered from convergence issues.[24] For NLME
models, which are widely-used in the field of pharmacokinetics, convergence issues are
not uncommon and can sometimes be overcome through the careful choice of starting
values.[16, 25]
This paper presents a novel Unified Hierarchical Bayesian (U-HB) model for NO
parameter estimation, based on previous statistical methodology developments in
pharmacokinetics[25] and implemented using Gibbs sampling.[26] Our U-HB model
simultaneously estimates NO parameters and their associations with covariates, thereby
propagating uncertainty throughout the entire analysis. The U-HB model is a novel
application of Hierarchical Bayesian methods to the field of multiple flow FeNO. We
present simulation studies evaluating the statistical properties of existing two-step
methods and the novel U-HB approach. We then show results from an application of
these methods to data from the southern California Children’s Health Study (CHS). The
CHS is a landmark cohort study on the effects of air pollution exposures on children’s
respiratory health.[27-29] The methodological work in this paper is motivated by the
need for better statistical methods to address CHS research questions. In particular, our
group has been interested in evaluating the effects of traffic-related air pollution
exposures on airway and alveolar inflammation. Traffic is a major source of
anthropogenic air pollution which is increasingly recognized as impacting human
health.[30-36]
16
Methods
We begin by introducing the deterministic steady-state two-compartment model for
FeNO and then discuss a variety of TS approaches which have been used to estimate both
the alveolar and airway NO parameters from this model and associations with covariates.
We then introduce our novel U-HB approach.
Deterministic, steady-state two-compartment model for FeNO
Our modeling work is based on the simple steady-state two-compartment model
(henceforth referred to as the 2CM) which assumes a cylindrically-shaped airway
compartment and an expansile alveolar compartment.[13, 37] Under the 2CM, FeNO can
be related to parameters quantifying airway and alveolar NO sources as follows:
FeNO = 𝐶 𝑎𝑤
+ ( 𝐶 𝐴 − 𝐶 𝑎𝑤
)× 𝑒 −𝐷 𝑎𝑤
flow ⁄
(1)
Equation 1 describes how FeNO (ppb) measured at the mouth is related to expiratory
flow rate, “flow” (mL/s), and three “NO parameters”: CA, the concentration of NO in the
alveolar region (ppb); Caw, the concentration of NO in the airway tissue (ppb); and Daw,
the airway tissue diffusion capacity (pL· s
-1
· ppb
-1
). Note that another widely used
parameter J’aw , the maximum flux of NO in the airway, is the product of Caw and
Daw.[13] A number of common methods for estimating NO parameters in the 2CM use an
alternative J’aw parameterization of Equation 1, but here we use the Caw parameterization
since it allows for direct estimation of a more interpretable parameter.
17
Estimating NO parameters from the 2CM
The model in Equation 1 is deterministic and nonlinear. To estimate 2CM parameters
from stochastically observed multiple flow FeNO data (steady state summaries of
repeated FeNO maneuvers at a range of target flow rates), researchers have relied on
linear regression methods under various linearizing assumptions or nonlinear regression
methods. For example, previous studies have used approximations based on Taylor
expansions to the exponential function, including first order linear approximation
methods (linT and linP),[13] second order quadratic approximation methods (quadT and
quadP),[16] and a third order approximation method, the Högman and Merilӓinen
algorithm (HMA).[18, 38] Nonlinear approaches include standard nonlinear least squares
regression[17] which essentially adds a random normal error to Equation 1, as well as a
natural log transform-both-sides extension.[16] The log transform-both-sides model
acknowledges the increased variation in residuals that occurs as flow rate (and hence
FeNO concentration) increases.
Estimating associations of covariates with NO parameters
TS methods. We selected three existing TS approaches to be evaluated in a simulation
study, along with the new U-HB method. In each of the following TS methods, Stage I
estimates of the three NO parameters for participant i (𝐶 𝐴 ̂
𝑖 , 𝑙𝑜𝑔𝐶 𝑎𝑤
̂
𝑖 , and 𝑙𝑜𝑔𝐷 𝑎𝑤
̂
𝑖 ) are
treated as known values and used as the outcomes in three separate Stage II linear
regression models relating each NO parameter to a covariate X:
18
𝐶 𝐴 ̂
𝑖 = 𝛼 𝐶 𝐴 + 𝛽 𝐶 𝐴 𝑋 𝑖 + 𝜀 𝐶 𝐴 𝑖 𝑙𝑜𝑔𝐶 𝑎𝑤
̂
𝑖 = 𝛼 𝑙𝑜𝑔𝐶 𝑎𝑤
+ 𝛽 𝑙𝑜𝑔𝐶 𝑎𝑤
𝑋 𝑖 + 𝜀 𝑙𝑜𝑔𝐶 𝑎𝑤
𝑖 𝑙𝑜𝑔𝐷 𝑎𝑤
̂
𝑖 = 𝛼 𝑙𝑜𝑔𝐷 𝑎𝑤
+ 𝛽 𝑙𝑜𝑔𝐷 𝑎𝑤
𝑋 𝑖 + 𝜀 𝑙𝑜𝑔𝐷 𝑎𝑤
𝑖 (2)
Henceforth, we will refer to the three TS approaches by the names of their Stage I
models, which are:
1. The TS Högman and Merilӓinen Algorithm (TS-HMA). HMA[18, 38] is a
widely-applied method for estimating NO parameters using an iterative algorithm
involving a third-order approximation to the 2CM. In a study with N participants,
Stage I consists of fitting N HMA models, one per participant. Stage I input data
for each participant consists of 3 observations: average FeNO measured at low,
medium, and high flow rates. In the CHS we use 30, 100, and 300 mL/s
respectively.[16]
2. The TS Nonlinear Least Squares model (TS-NLS). The natural log transform-
both-sides nonlinear least squares model proposed previously by our group,[16]
and henceforth referred to as NLS, is implemented using standard nonlinear-least
squares software (“nls” from the nlme package in R). The log transform-both-
sides approach better satisfies the assumption of normally distributed errors while
maintaining the physiological interpretation of the NO parameters.[16] For N
participants, Stage I consists of fitting N log transform-both-sides NLS models of
the following form, based on Equation 1:
log ( FeNO) = log ( exp( 𝑙𝑜𝑔𝐶 𝑎𝑤
)+ ( 𝐶 𝐴 − exp( 𝑙𝑜𝑔𝐶 𝑎𝑤
) )× 𝑒 −exp( 𝑙𝑜𝑔𝐷 𝑎𝑤
) 𝑓𝑙𝑜𝑤 𝑗 ⁄
)+ 𝜀 𝑗
(3)
19
where j=1…n, indexes the observations for each participant. In the CHS, the
protocol asked each participant to perform 9 maneuvers, but in practice
participants had between 4-12 valid maneuvers each. Each maneuver is
summarized by a paired observation of FeNO concentration and flow.[39] Note
that in Equation 3, two NO parameters are expressed (and directly estimated) as
log𝐶 𝑎𝑤
and log𝐷 𝑎𝑤
to enforce positivity of 𝐶 𝑎𝑤
and 𝐷 𝑎𝑤
and to better satisfy
the normality assumptions in the TS-NLME method (below). This is common
practice in pharmacokinetics modelling.[40]
3. The TS Nonlinear Mixed Effect model (TS-NLME): In a nonlinear mixed effects
(NLME) approach, we use FeNO maneuvers from all participants to estimate a
single NLME model of the form:
lo g( FeNO
𝑖𝑗
) = log( exp( 𝑙𝑜𝑔𝐶 𝑎𝑤
𝑖 )+ ( 𝐶 𝐴 𝑖 − exp( 𝑙𝑜𝑔𝐶 𝑎𝑤
𝑖 ) )× 𝑒 −exp( 𝑙𝑜𝑔𝐷 𝑎𝑤
𝑖 ) 𝑓𝑙𝑜𝑤 𝑖𝑗
⁄
)+
𝜀 𝑖𝑗
(4)
with participant-level random effects which follow a multivariate normal
distribution. Participant-level estimates of NO parameters are obtained by
combining fixed effect parameter estimates with empirical Bayes estimates of the
random effects. Because TS-NLME fits a single model (rather than N models) in
Stage I, it does not suffer from the small sample size issues that affect TS-HMA
or TS-NLS (each of which often fails to produce estimates for a subset of the
population at Stage I). In TS-NLME, which has been applied previously,[22] NO
parameters estimated in a Stage I NLME are then related to covariates in a
separate Stage II.
20
In secondary analyses, we evaluated two alternative TS-NLS approaches: a constrained
version of TS-NLS (in which we assume 𝐶 𝐴 >0.001) and a version of TS-NLS with an
inverse-variance weighted Stage II.
Unified methods. In our unified approaches, we aim to model the measured outcome
FeNO conditionally on both latent NO parameters and measured environmental factors
using hierarchical modeling. Unified methods, in contrast to the TS methods,
simultaneously estimate the NO parameters and their associations with the covariate 𝑋 𝑖
in a single model. Below we present two unified methods, one based on the frequentist
NLME approach and our novel method using a Bayesian approach:
1. The Unified Nonlinear Mixed Effect model (U-NLME): In U-NLME, rather than
estimate separate Stage II models relating estimated NO parameters to 𝑋 𝑖 as in
TS-NLME, we incorporate 𝑋 𝑖 into the mean function for each NO parameter in
the NLME model (equations are conceptually similar to those in the U-HB,
presented in Equations 5-7). Recall that NLME is recognized to have convergence
issues and is sensitive to starting values.[16, 25] Early work on the U-NLME
indeed demonstrated problems with convergence.[24] However, for completeness
we include U-NLME for comparison purposes.
2. The Unified Hierarchical Bayesian model (U-HB): The U-HB method can be
described as a two-level hierarchical Bayesian model, with j indexing FeNO
maneuvers nested in participant i, as displayed in Figure 3.
21
Level 1: Maneuver
Similar to Equation 4, we assume that log(FeNO) for participant i at maneuver j is
normally distributed:
log ( 𝐹𝑒𝑁𝑂 𝑖𝑗
)∼ 𝑁 ( 𝑓 ( 𝜃 𝑖 , 𝑋 𝑖 , 𝑓𝑙𝑜𝑤 𝑖𝑗
) , 𝜎 2
) (5)
with a mean function:
𝑓 ( 𝜃 𝑖 , 𝑋 𝑖 , 𝑓𝑙𝑜𝑤 𝑖𝑗
)= log( exp( 𝑙𝑜𝑔𝐶 𝑎𝑤
𝑖 )+ ( 𝐶 𝐴 𝑖 − exp( 𝑙𝑜𝑔𝐶 𝑎𝑤
𝑖 ) )×
𝑒 −exp( 𝑙𝑜𝑔𝐷 𝑎𝑤
𝑖 ) 𝑓𝑙𝑜𝑤 𝑖𝑗
⁄
) (6)
Figure 3: Hierarchical model structure relating FeNO measurements at multiple flow rates to NO
parameters that are a function of a potential determinant X (e.g., air pollution).
22
that depends on the 3-dimensional vector of participant-specific NO parameters 𝜃 𝑖 =
( 𝐶 𝐴 𝑖 , 𝑙𝑜𝑔𝐶 𝑎𝑤
𝑖 , 𝑙𝑜𝑔𝐷 𝑎𝑤
𝑖 ) ′, 𝑋 𝑖 (the role of 𝑋 𝑖 becomes clear in Level 2), and flow rates. The
variance (σ
2
) is assumed constant across participants and flow rates.
Level 2: Subject
We assumed the participant-level NO parameters are each a linear function of 𝑋 𝑖 :
𝜃 𝑖 = [
𝐶 𝐴 𝑖 𝑙𝑜𝑔𝐶 𝑎𝑤
𝑖 𝑙𝑜𝑔𝐷 𝑎𝑤
𝑖 ] = [
𝛼 𝐶 𝐴 + 𝛽 𝐶 𝐴 𝑋 𝑖 + 𝜖 𝐶 𝐴 𝑖 𝛼 𝑙𝑜𝑔𝐶 𝑎𝑤
+ 𝛽 𝑙𝑜𝑔𝐶 𝑎𝑤
𝑋 𝑖 + 𝜖 𝑙𝑜𝑔𝐶 𝑎𝑤
𝑖 𝛼 𝑙𝑜𝑔𝐷 𝑎𝑤
+ 𝛽 𝑙𝑜𝑔𝐷 𝑎𝑤
𝑋 𝑖 + 𝜖 𝑙𝑜𝑔𝐷 𝑎𝑤
𝑖 ] (7)
where the vector 𝜖 𝑖 = ( 𝜖 𝐶 𝐴 𝑖 , 𝜖 𝑙𝑜𝑔𝐶 𝑎𝑤
𝑖 , 𝜖 𝑙𝑜𝑔𝐷 𝑎𝑤
𝑖 ) ′ of participant-level random effects is
assumed to have a multivariate normal distribution with mean zero and variance-
covariance matrix Σ
𝜃 . Based on the previous data analyses,[41] we assumed the joint
distribution of the NO parameters can be reasonably modeled using a multivariate normal
(MVN) distribution, and the equations above are equivalent to a formulation of:
𝜃 𝑖 ~𝑀𝑉𝑁 ( 𝜇 𝜃 𝑖 , Σ
𝜃 ) where 𝜇 𝜃 𝑖 represents a mean vector with, for example, the first
element equal to 𝛼 𝐶 𝐴 + 𝛽 𝐶 𝐴 𝑋 𝑖 , and Σ
𝜃 is the same variance-covariance matrix of NO
parameters. Since the concentration of NO in the alveolar region should be non-negative,
we imposed a non-negative constraint on 𝐶 𝐴 . We implemented our U-HB model using
the Gibbs sampling program JAGS.[42] But JAGS can only specify univariate truncated
normal distributions. To achieve this constraint in our multivariate context, we
constructed the truncated MVN distribution into two steps. First, we sampled from a
bivariate normal distribution for (𝑙𝑜𝑔𝐶 𝑎𝑤
, 𝑙𝑜𝑔𝐷 𝑎𝑤
), then we sampled from a zero-
23
truncated normal distribution for 𝐶 𝐴 with the conditional expectation and variance given
(𝑙𝑜𝑔𝐶 𝑎𝑤
, 𝑙𝑜𝑔𝐷 𝑎𝑤
).
Prior distributions
We assumed the vectors of regression intercepts (𝛼 ) and slopes (𝛽 ) in Equation 7 each
had independent normal prior distributions (with I indicating a square identity matrix):
𝛼 ∼ 𝑀𝑉𝑁 ( 𝜇 𝛼 , 𝐼 𝜎 𝛼 2
)
𝛽 ∼ 𝑀𝑉𝑁 ( 𝜇 𝛽 , 𝐼 𝜎 𝛽 2
)
(8)
where 𝜇 𝛼 = ( 𝜇 𝛼 𝐶 𝐴 , 𝜇 𝛼 𝑙𝑜𝑔𝐶 𝑎𝑤
, 𝜇 𝛼 𝑙𝑜𝑔𝐷 𝑎𝑤
)
′
= (2, 4.2, 2.5) ′ and 𝜇 𝛽 =
( 𝜇 𝛽 𝐶 𝐴 , 𝜇 𝛽 𝑙𝑜𝑔𝐶 𝑎𝑤
, 𝜇 𝛽 𝑙𝑜𝑔𝐷 𝑎𝑤
) ′= (0, 0, 0). Non-informative values were assumed for the
variance hyperparameters 𝜎 𝛼 2
= 𝜎 𝛽 2
= 10
3
. Similarly, the variance of the residuals (σ
2
)
was assigned a non-informative Inv-Gamma(10
−3
, 10
−3
) prior distribution and the
variance-covariance matrix of NO parameters (Σ
𝜃 ) was assigned a non-informative
inverse-Wishart prior distribution.
Calculation of posterior via MCMC:
Taking a Bayesian perspective, the general posterior distribution can be written as:
𝑝 ( 𝜇 𝛼 , 𝜎 2
𝛼 , 𝜇 𝛽 , 𝜎 2
𝛽 , Σ
Θ
|𝑙𝑜𝑔𝐹𝑒𝑁𝑂 𝑖𝑗
, 𝑋 𝑖 , 𝑓𝑙𝑜𝑤 𝑖𝑗
)∝
𝑝 ( log𝐹𝑒𝑁𝑂 𝑖𝑗
|𝑓 ( θ
𝑖 , 𝑓𝑙𝑜𝑤 𝑖𝑗
) , 𝜎 2
) 𝑝 ( θ
𝑖 |𝜇 θ
𝑖 , Σ
θ
) 𝑝 ( 𝜇 θ
𝑖 |𝛼 , 𝛽 , 𝑋 𝑖 ) 𝑝 ( 𝛼 |𝜇 𝛼 , 𝜎 2
𝛼 ) 𝑝 ( 𝛽 |𝜇 𝛽 , 𝜎 2
𝛽 ) 𝑝 ( 𝜎 2
)
(9)
24
This density cannot be directly calculated. Instead, we rely on MCMC methods to
provide samples from the density. Specifically, we implemented Gibbs sampling using
JAGS, taking advantage of normal conjugate distributions for this hierarchical model.
While convergence is sometimes non-trivial to obtain, Gibbs sampling generally had
good performance in the models considered here.
To conduct U-HB analyses, we used three parallel MCMC chains in JAGS. Once the
Gelman-Rubin 𝑅 ̂
convergence criteria[43] reached ≤1.1, we drew 12,000 additional
samples and checked that the 𝑅 ̂
remained ≤1.1 (second checkpoint). If so, we used those
12,000 samples to construct the posterior distributions of the parameters. Otherwise, we
continued sampling as long as necessary to satisfy the convergence criteria.
Simulation study
We conducted an extensive simulation study to compare the performance of the five
methods. To produce realistic synthetic data, we modelled our data generation scenarios
on the CHS. Specifically, we simulated data for 1000 individuals and assumed each
individual performed 8 FeNO maneuvers, two at each of four flow rates: 30, 50, 100, and
300 ml/s. These particular flow rates were selected for three reasons: (1) they match the
flow restrictors provided with the FeNO sampling instrument, (2) they are within a
reasonable range of flows for schoolchildren, and (3) optimal flow rate sampling design
has been studied in detail previously, and studies with these flow rates balanced
theoretical performance with feasibly.[39] Based on previously described distributions in
25
the CHS,[41] we assumed NO parameters had a multivariate normal distribution with an
additional non-negative constraint for 𝐶 𝐴 in data generation step (the randomly sampled
vector of NO parameters was discarded if it had a negative 𝐶 𝐴 value). The variance-
covariance matrix was set to be similar to values observed in preliminary analyses of
CHS data:
Σ
Θ
= (
0.44 0.13 −0.14
0.13 0.62 −0.15
−0.14 −0.15 0.36
) (10)
We also assumed that a standard normal covariate X had a potentially different linear
relationship with each NO parameter. In all data generating scenarios, we set the fixed
effect intercepts (mean NO parameters for participants with mean X) to mirror those used
in previous statistical methods work on FeNO in the CHS (i.e., 1, 3.5, and 2.5 for CA,
logCaw and logDaw respectively). Our primary focus was on the associations of X with
NO parameters (𝛽 𝐶 𝐴 , 𝛽 𝑙𝑜𝑔𝐶 𝑎𝑤
, and 𝛽 𝑙𝑜𝑔𝐷 𝑎𝑤
).
We considered seven different scenarios, in which one or more of the 𝛽 ’s varied (ranging
from 0 to 0.2, with a step size of 0.02) as summarized in Table 1. In Scenario 1, the
magnitudes of the associations between X and each NO parameters were assumed to be
equal (𝛽 𝐶 𝐴 = 𝛽 𝑙𝑜𝑔𝐶 𝑎𝑤
= 𝛽 𝑙𝑜𝑔𝐷 𝑎𝑤
) and all 𝛽 ’s varied together. In Scenarios 2-4, only one of
the NO parameters was associated with X; in Scenarios 5-7, only one NO parameter was
not associated with X. The scenarios where at least one association is truly zero permit
estimation of Type I error rates. For each setting of each scenario, we conducted 1000
replicates. Using these seven scenarios, we compared the five methods in terms of the
26
following criteria. Bias was calculated as (estimate – true value) and relative bias as
(estimate– true value)/true value. Considering 95% confidence intervals or credible
intervals (henceforth called 95% CI for simplicity), we calculated their length and
coverage. Power was defined as the proportion of the 95% CIs that did not include 0
when the true value was not 0. We used Scenario 1 for primary results, where all NO
parameters were assumed to be equally affected by X. We used Scenarios 5-7 for
primary results on type I error rates, which we define as the proportion of 95% CIs that
did not include 0 when the true value was 0. All results were presented only for simulated
datasets on which all methods satisfied convergence criteria.
Table 1: Seven scenarios considered in the simulation study, where the relation of X to
each NO parameter (𝜷 𝑪 𝑨 , 𝜷 𝒍𝒐𝒈𝑪 𝒂𝒘
, 𝜷 𝒍𝒐𝒈𝑫 𝒂𝒘
) varied from 0 to 0.2, with a step size of 0.02.
𝜷 𝑪 𝑨 𝜷 𝒍𝒐𝒈𝑪 𝒂𝒘
𝜷 𝒍𝒐𝒈𝑫 𝒂𝒘
Scenario 1
*
0 to 0.2 0 to 0.2 0 to 0.2
Scenario 2 0.02 to 0.2 0
†
0
Scenario 3 0 0.02 to 0.2 0
Scenario 4 0 0 0.02 to 0.2
Scenario 5
*
0 0.02 to 0.2 0.02 to 0.2
Scenario 6
*
0.02 to 0.2 0 0.02 to 0.2
Scenario 7
*
0.02 to 0.2 0.02 to 0.2 0
*
In Scenarios 1 and 5-7, the non-zero β values are identical (e.g., the first three settings of Scenario 1 have 𝛽 𝐶 𝐴 =
𝛽 𝑙𝑜𝑔𝐶 𝑎𝑤
= 𝛽 𝑙𝑜𝑔𝐷 𝑎𝑤
= 0, 𝛽 𝐶 𝐴 = 𝛽 𝑙𝑜𝑔𝐶 𝑎𝑤
= 𝛽 𝑙𝑜𝑔𝐷 𝑎𝑤
= 0.02, 𝛽 𝐶 𝐴 = 𝛽 𝑙𝑜𝑔𝐶 𝑎𝑤
= 𝛽 𝑙𝑜𝑔𝐷 𝑎𝑤
= 0.04).
†
Cells marked 0 indicate that X had no effect on the corresponding NO parameter.
27
In total, we generated 71 sets of simulated datasets (7 scenarios with 10 values each of
the varying parameters, as well as one dataset with all 𝛽 ’s=0), each replicated 1,000
times. Analyses were conducted using the clusters at the University of Southern
California’s High-Performance Computing Center using R version 3.5 and JAGS 4.0.
The median time needed to produce estimates for 5 models (TS-NLS, TS-HMA, TS-
NLME, U-NLME, U-HB) for each dataset was 6 hours.
CHS data analysis
In a previous publication using data from CHS participants,[41] we investigated the
association between traffic-related air pollution (TRAP) and NO parameters. Briefly, we
conducted a cross-sectional analysis of multiple flow FeNO measured in 1635
schoolchildren ages 12-15, using TS-NLME and TS-HMA. By design, CHS participants
had ~9 FeNO maneuvers each (3 at 50 ml/s, 2 each at 30, 100, and 300 ml/s flow rates).
Exposure to TRAP in the indoor schoolroom microenvironment at the time of the FeNO
test was assessed using approximately concurrent room air measurements of NO (ppb). In
this paper, we re-analyzed these data using both “unadjusted” models (with TRAP
exposure as the only covariate) and minimally “adjusted” models (with TRAP exposure
plus adjustment for: sex, age, asthma). The original publication adjusted for additional
potential confounders or study design variables (race/ethnicity, rhinitis history, use of
asthma medications in the past 12 months, secondhand tobacco smoke exposure, parental
education, time of day of FeNO test, FeNO analyzer, and CHS community), but results
from minimally adjusted models were similar to those of the fully adjusted model.[41]
28
Note that the original publication used the alternative J’aw parameterization of the two-
compartment model:
log ( FeNO) = log ( exp( 𝑙𝑜𝑔𝐽 ‘
𝑎𝑤
− 𝑙𝑜𝑔𝐷 𝑎𝑤
)+ ( 𝐶 𝐴 − exp( 𝑙𝑜𝑔𝐽 ‘
𝑎𝑤
− 𝑙𝑜𝑔𝐷 𝑎𝑤
) )∗
exp( −exp( 𝑙𝑜𝑔𝐷 𝑎𝑤
) flow ⁄ )+ 𝜀 (8)
whereas in this paper we used the “Caw” parameterization (Equation 2) due to the direct
estimation of the more interpretable airway wall concentration parameter. Thus, here we
fit TS-NLME models to CHS data with both the “J’aw” and “Caw” parameterization to
facilitate direct comparison with the previously published work.
29
Results
Simulation study:
Computational Time and Convergence. TS-NLS and TS-HMA models always converged
in Stage II and their run times were typically less than 1 minute, but they failed to
estimate NO parameters in Stage I for some participants due to small sample size. In
average, 31.6% of the Stage I models failed to converge for TS-NLS while it was <1%
for TS-HMA. Hence these participants were left out of the Stage II model. TS-NLME
models had a convergence rate of 90% (average time: ~10 minutes) while U-NLME had
a convergence rate of only 70% (average time: ~20 minutes). For U-HB models, the
median time to achieve 𝑅 ̂
<1.1 was 5.5 hours (~90% of models converged within 16
hours and 1% never converged in the 48 hours allowed). For the rest of this paper, to aid
comparison, we only show results obtained from the subset of datasets for which all
methods satisfied converged criteria. Results for all datasets are shown in Supplementary
Section 5.
Bias. As shown in Figure 4a (Scenario 1), U-HB consistently showed the smallest bias.
For other methods, the directions and relative magnitudes of the bias varied across NO
parameters. Trends in bias were consistent across the range of 𝛽 , so we report average
relative bias. All methods produced negatively biased estimates of 𝛽 𝐶 𝐴 , with the average
relative bias being the smallest for U-HB ( -4.1%) and much larger for the other methods
(approximately -12.3% for TS-HMA and -50.4% for TS-NLS, -42.6% for TS-NLME and
-42.3% for U-NLME). For 𝛽 𝑙𝑜𝑔𝐶 𝑎𝑤
, average relative bias was negative for U-HB (-7.7%)
30
and TS-HMA (-67.6%) but positive for the other methods (11.1% for U-NLME, 11.4 %
for TS-NLS, and 45.1% for TS-NLME). Conversely, for 𝛽 𝑙𝑜𝑔 𝐷 𝑎𝑤
, the average relative
bias was positive for U-HB (8.8%) and TS-HMA (64.1%) but negative for the other
methods (-10.2% for U-NLME, -53.4 for TS-NLS, and -55.2% for TS-NLME). For
Scenarios 2-7 (Supplementary Figures 2.2-2.7), U-HB also had smaller bias than all other
methods. In the secondary analyses, constrained TS-NLS had smaller bias than the
standard TS-NLS for estimating 𝐶 𝐴 , but, it considerably underestimated logCaw. Inverse-
variance weighted TS-NLS had even worse performance (results not shown). Thus, we
only report standard TS-NLS results.
CI coverage and length. As shown in Figure 4b and 4d (Scenario 1), U-HB was the only
method to produce CIs with appropriate coverage (~95%) for associations with all NO
parameters (𝛽 𝐶 𝐴 , 𝛽 𝑙𝑜𝑔𝐶 𝑎𝑤
, 𝛽 𝑙𝑜𝑔 𝐷 𝑎𝑤
). TS-NLME consistently had the worst coverage
while its one-step analog, U-NLME, performed much better for 𝛽 𝑙𝑜𝑔𝐶 𝑎𝑤
and 𝛽 𝑙𝑜𝑔 𝐷 𝑎𝑤
but
only moderately better for 𝛽 𝐶 𝐴 . To further understand this difference in TS-NLME vs U-
NLME coverage, we note that the CIs for TS-NLME were very narrow (the narrowest of
all methods) and TS-NLME estimates were considerably more biased than for U-NLME
(except for 𝛽 𝐶 𝐴 , where the bias was similar for the two methods). TS-NLS had poor
coverage for 𝛽 𝐶 𝐴 (due to large bias, despite wide CI) and 𝛽 𝑙𝑜𝑔 𝐷 𝑎𝑤
(due to large bias).
TS-HMA had the widest CIs for 𝛽 𝑙𝑜𝑔 𝐶𝑎𝑤 and 𝛽 𝑙𝑜𝑔 𝐷𝑎𝑤
but the bias for these parameters
was also large, reducing the CIs’ coverage. For Scenarios 2-7, results were similar
(Supplementary Figures 2.2-2.7).
31
Figure 4: Relative bias (a), coverage (b), power, (c) and CI length (d) of the selected estimation
methods from Scenario 1 of the simulation study (𝜷 𝑪 𝑨 = 𝜷 𝒍𝒐𝒈𝑪 𝒂𝒘
= 𝜷 𝒍𝒐𝒈𝑫 𝒂𝒘
).
4a.
4b.
4c.
4d.
32
Type I error. Figure 5 shows type I error rates from Scenarios 5, 6 and 7, all of which had
only one NO parameter not affected by the covariate X (e.g., type I error for 𝛽 𝐶 𝐴 was
obtained from Scenario 5 where 𝛽 𝐶 𝐴 = 0 while 𝛽 𝑙𝑜𝑔𝐶 𝑎𝑤
= 𝛽 𝑙𝑜𝑔𝐷 𝑎𝑤
varied). Results of
Scenarios 2-4, in which only one NO parameter association was non-zero, are shown in
Supplementary Figures 2.2-2.4. Of all the methods, U-HB generally had the lowest type I
error rates, with values of ~0.05 for 𝛽 𝐶 𝐴 , 𝛽 𝑙𝑜𝑔𝐶 𝑎𝑤
, and 𝛽 𝑙𝑜𝑔 𝐷 𝑎𝑤
but that increased slightly
as the effect of X on the other NO parameters increased. As might be expected from the
earlier analyses of bias and CI coverage/length, the other methods performed less well.
TS-NLME had poor type I errors for 𝛽 𝑙𝑜𝑔𝐶 𝑎𝑤
and 𝛽 𝑙𝑜𝑔 𝐷 𝑎𝑤
, and the type I errors
increased dramatically as the effect of X increased. U-NLME had the highest type I error
for 𝛽 𝐶 𝐴 , and the type I error increased more modestly as the effect of X increased. TS-
NLS and TS-HMA had similar patterns in type I error rates, which were inflated slightly
for 𝛽 𝐶 𝐴 and 𝛽 𝑙𝑜𝑔𝐶 𝑎𝑤
but ~0.05 for 𝛽 𝑙𝑜𝑔 𝐷 𝑎𝑤
.
33
Power. As shown in Figure 4c (Scenario 1), power results were complex. U-HB had
power >0.99 for effect sizes of 0.12 for 𝛽 𝐶 𝐴 , 0.18 for 𝛽 𝑙𝑜𝑔 𝐶𝑎𝑤 , and 0.16 for 𝛽 𝑙𝑜𝑔 𝐷𝑎𝑤
.
Although U-HB never had the best power relative to the other methods, the apparently
good power curves of some methods require careful interpretation. TS-NLME appeared
to have the best power across all NO parameters, but this was due to its narrow
confidence intervals and came at the cost of inflated type-I error rates (e.g., Figure 4c for
Scenario 1 when 𝛽 𝐶 𝐴 = 𝛽 𝑙𝑜𝑔𝐶 𝑎𝑤
= 𝛽 𝑙𝑜𝑔𝐷 𝑎𝑤
=0). TS-HMA had a similar power curve to U-
HB for 𝛽 𝐶 𝐴 (the parameter for which TS-HMA had the least bias), but TS-HMA had low
power curves for 𝛽 𝑙𝑜𝑔 𝐶𝑎𝑤 and 𝛽 𝑙𝑜𝑔 𝐷𝑎𝑤
(parameters for which TS-HMA had displayed
considerable bias and wide CI). Relative to the other methods, U-NLME power was high
for 𝛽 𝑙𝑜𝑔 𝐶𝑎𝑤 and 𝛽 𝑙𝑜𝑔 𝐷𝑎𝑤
, but poor for 𝛽 𝐶 𝐴 (the parameter for which U-NLME had
Figure 5: Type I error
*
of the selected estimation methods from the simulation study.
34
considerable bias). TS-NLS had poor power for 𝛽 𝐶 𝐴 and 𝛽 𝑙𝑜𝑔 𝐷𝑎𝑤
but better power for
𝛽 𝑙𝑜𝑔 𝐶𝑎𝑤 (the parameter for which TS-NLS had the least bias).
Overall summary of simulation study results. Across our simulation study scenarios, U-
HB had good statistical properties and consistently outperformed U-NLME and the three
TS approaches. The other unified approach, U-NLME, performed well for 𝛽 𝑙𝑜𝑔𝐶 𝑎𝑤
and
𝛽 𝑙𝑜𝑔 𝐷 𝑎𝑤
, but not for 𝛽 𝐶 𝐴 (biased, poor coverage, high type I error) and suffered from
serious convergence issues, failing to converge for 30.2% of simulated datasets
(Supplementary Figure 1). The worst bias was observed for TS methods. Interestingly,
TS-HMA had a different direction of bias compared to TS-NLS, TS-NLME and U-
NLME for 𝛽 𝑙𝑜𝑔𝐶 𝑎𝑤
and 𝛽 𝑙𝑜𝑔 𝐷 𝑎𝑤
. The correlation between NO parameters had an
influence on bias across all methods. We assumed a positive correlation between CA and
logCaw, a negative correlation between CA and logDaw, and a negative correlation between
logCaw and logDaw to mimic the distributions previously observed in CHS participants.
Neither the magnitude nor the direction of the estimated relative bias for 𝛽 𝐶 𝐴 from all the
models were affected by the values of 𝛽 𝑙𝑜𝑔𝐶 𝑎𝑤
or 𝛽 𝑙𝑜𝑔 𝐷 𝑎𝑤
(comparing bias figures for
Scenario 1 to Scenarios 2, 6 and 7, recall all Scenarios are included in Supplementary
Figures 2.1-2.7). But when estimating 𝛽 𝑙𝑜𝑔𝐶 𝑎𝑤
and 𝛽 𝑙𝑜𝑔 𝐷 𝑎𝑤
, the relative bias changed
both in magnitude and direction, except for U-HB (comparing bias figures for Scenario 1
to Scenarios 3, 5 and 7; comparing Scenario 1 to Scenarios 4, 5 and 7). We focused the
simulation study results on 𝛽 ’s, the association between X and NO parameters, as our
main interest. Results for 𝛼 ’s, the population mean NO parameters when X=0, are
presented in Supplementary Figures 3.1-3.7.
35
CHS data analysis:
As shown in Figure 6, we confirmed the previously published result (using the J’aw
parameterization TS-NLME)[41] of a statistically significant positive association
between TRAP and alveolar NO (𝛽 𝐶 𝐴 ) using all estimation methods. Results were quite
similar with and without adjustment for age, sex, and asthma status (Table 2). Using U-
HB, we found that a 10 ppb higher TRAP concentration was associated with a 0.15 (95%
CI: 0.07, 0.24) ppb increase in CA, after adjusting for sex, age, and asthma. This U-HB
Figure 6: Analysis of CHS data: estimated associations of traffic-related air pollution with
CA, logCaw, and logDaw (𝜷 ̂
and 95% CI) using the selected methods,
*
with no adjustments for
covariates.
*
TS-NLME_J denotes the previously published model using a J’aw parameterization of TS-NLME.
36
estimated 𝛽 𝐶 𝐴 was nearly twice the size of the analogous estimate from TS-NLME (0.08,
95% CI: 0.03, 0.14 from both parameterizations). As in the previous publication, there
was little evidence for effects of TRAP on airway wall inflammation (𝛽 𝑙𝑜𝑔𝐶 𝑎𝑤
and
𝛽 𝑙𝑜𝑔 𝐷 𝑎𝑤
). Based on the CHS findings, Scenario 2 of the simulation study (𝛽 𝐶 𝐴 varies,
𝛽 𝑙𝑜𝑔𝐶 𝑎𝑤
= 0, 𝛽 𝑙𝑜𝑔𝐷 𝑎𝑤
= 0, Supplemental Figure 2.2) might be a closer comparator to the
CHS data than Scenario 1. The finding in the CHS data that the estimated 𝛽 𝐶 𝐴 was larger
for U-HB than for TS-NLME agrees with results from both simulation Scenarios 1 and 2,
where U-HB had a modest negative bias for 𝛽 𝐶 𝐴 compared to the more considerable
negative bias for TS-NLME. While we were motivated by the problem of estimating
TRAP associations with NO parameters (Stage II), it is interesting to note that all
methods produce participant-level NO parameter estimates (Stage I) and these estimates
were most similar when comparing U-HB to TS-NLME and U-NLME (Supplemental
Figures 4.1-4.3). However, TS-NLS and TS-HMA failed to estimate NO parameters in
Stage I for subsets of participants due to small sample size of flow rate (31.6% for TS-
NLS and <1% for TS-HMA in average).
37
Table 2: Analysis of CHS data: estimated associations of a 10 ppb increase in traffic-
related air pollution with CA, logCaw, and logDaw (𝜷 ̂
and 95% CI) using the selected
methods, without and with adjustments for covariates (age, sex, asthma).
Estimation
method
Model Estimated 𝜷 𝑪 𝑨
(95% CI)
Estimated 𝜷 𝒍𝒐𝒈𝑪 𝒂𝒘
(95% CI)
Estimated 𝜷 𝒍𝒐𝒈𝑫 𝒂𝒘
(95% CI)
TS-NLS
Unadjusted 0.24 (0.06, 0.41)
**
0.04 (-0.05, 0.13) 0.02 (-0.05, 0.10)
Adjusted 0.23 (0.06, 0.40)
**
0.04 (-0.05, 0.13) 0.03 (-0.05, 0.10)
TS-HMA
Unadjusted 0.33 (0.19, 0.47)
***
0.17 (0.02, 0.32)
*
-0.14 (-0.31, 0.02)
Adjusted 0.32 (0.18,0.46)
***
0.16 (0.02, 0.31)
*
-0.14 (-0.31 0.02)
TS-NLME
Unadjusted 0.08 (0.03, 0.14)
**
0.02 (-0.04, 0.07) 0.04 (-0.01, 0.08)
Adjusted 0.08 (0.03, 0.14)
**
0.01 (-0.04, 0.07) 0.03 (-0.01, 0.08)
TS-NLME (J’aw)
†
Unadjusted 0.08 (0.03, 0.14)
**
0.07 (-0.01, 0.15) 0.03 (0.01, 0.08)
Adjusted 0.08 (0.03, 0.14)
**
0.07 (-0.01, 0.14) 0.03 (-0.01, 0.07)
U-NLME
Unadjusted 0.13 (0.04, 0.22)
**
-0.01 (-0.09, 0.06) 0.07 (-0.01, 0.15)
Adjusted 0.13 (0.05, 0.22)
**
-0.03 (-0.09, 0.06) 0.07 (-0.01, 0.15)
U-HB
Unadjusted 0.15 (0.07, 0.24) -0.02 (-0.11, 0.07) 0.07 (-0.03, 0.18)
Adjusted 0.15 (0.07, 0.24) -0.02 (-0.12, 0.07) 0.07 (-0.04, 0.18)
†
Previously published model: J’aw parameterization of TS-NLME.
*
: P<0.05,
**
: P<0.01,
***
: P<0.001
38
Discussion
In this paper, we performed an extensive simulation study to evaluate the finite sample
performance of a variety of statistical methods used to estimate the associations of
covariates with NO parameters from the steady-state two compartment model of lower
respiratory tract NO. One of these methods was a novel Unified Hierarchical Bayesian
(U-HB) model which simultaneously estimates participant-level NO parameters and their
associations with covariates. In the simulation studies, U-HB outperformed four other
methods: three two-stage methods and one unified method, all implemented using
frequentist approaches. When applying U-HB to the motivating data example—
investigating associations of traffic-related pollution exposure with NO parameters in a
cross-sectional sample of southern California schoolchildren—we confirmed the
previously published positive association of traffic with alveolar inflammation (CA).
However, the estimated association was nearly two times larger (0.15 vs 0.08) using U-
HB as compared to the two-stage method (TS-NLME) applied in the previous
publication. This result was consistent with the simulation study finding that estimates of
𝛽 𝐶 𝐴 tended to be higher (less negative bias) in U-HB than in TS-NLME.
Numerous applied publications use the two-stage approach, with NO parameters
estimated in Stage I and then treated as observed outcomes in Stage II analyses relating
NO parameters to covariates. Our paper is, to our knowledge, the first evaluation of
statistical methods used to estimate the associations of covariates with NO parameters. A
number of publications have previously studied statistical methods for estimating
participant-level NO parameters (i.e., only Stage I) [13, 16, 18] In our own previous
39
comparison of Stage I methods, we found that NLS outperformed linear, quadratic, and
third order approximations (i.e., HMA), especially in terms of reduced bias for estimates
of D
𝑎𝑤
𝑁𝑂 and more appropriate CI coverage for all NO parameters.
19
In this paper
where we compare TS methods, we found that TS-NLS also generally outperformed TS-
HMA. The direction of bias for TS-NLS and TS-NLME were always the same, while TS-
HMA usually differed from the other methods. TS-HMA also had the widest CI while
TS-NLME had the shortest CI (resulting in high power but low coverage). The worst bias
was observed for TS methods, which we suspect was a consequence of not propagating
uncertainty in the second stage. U-NLME, which had the second best performance after
U-HB, had serious convergence issues, as we had reported in our early work on U-
NLME.[24]
This paper has several strengths. First, our novel U-HB model takes a unified approach,
which combines the estimation of NO parameters and their associations with covariates
into a single step. This allows for the propagation of uncertainty across stages and avoids
the exclusion of some participants in two-stage models due to failure to estimate their
Stage I model. Second, we used a Bayesian approach for U-HB which allowed us to:
clearly specify a hierarchical model structure, fully characterize the posterior
distributions of the parameters of interest, and to have the flexibility to use diffuse priors
(as we do here) or more informative priors if appropriate. Unlike TS-NLS or TS-NLME,
U-HB does not use inferential approaches that rely on large sample size/normality
assumptions that might be violated with multiple flow FeNO datasets. Third, we used the
Caw parameterization of the 2CM (Equation 1) for all the methods presented in the paper
40
since it directly estimates what we believe to be a more interpretable NO parameter—the
airway wall concentration of NO (Caw, ppb)—rather than J’aw (the airway wall tissue
diffusion capacity, pL· s
-1
· ppb
-1
). Based on CHS data analysis (Supplementary Figure
4.3), had we used the J’aw parameterization, our results would have been similar to what
we observed using the Caw parameterization. An additional advantage of U-HB is that,
since it is implemented using MCMC, we can fully characterize the posterior of J’aw from
a Caw parameterization model by simply calculating exp(logCaw)/exp(logDaw) at each
MCMC iteration. Fourth, we compared U-HB to four alternative methods, which spanned
the range applied in practice to study the associations of covariates with NO parameters.
We did exclude one common two-stage method (“TS-linT”) in which Stage I is a simple
linear regression for each participant, based on the linear approximation proposed by
Tsoukias et al.[44] We excluded TS-linT because its simplified assumption of the
approximation cannot estimate Caw (it only estimates J’aw and CA). Fifth, we performed
an extensive set of simulation studies across a range of effect sizes, which were of
particular concern in situations where the covariate is related to only one (or two) of the
three correlated NO parameters. We examined not only bias and CI length but also power
and type I error rates in those scenarios to understand how each model behaved. Finally,
we implemented U-HB via Gibbs sampling in readily available software (the R interface
to JAGS). This implementation of the U-HB worked well across simulation scenarios
(especially compared to other methods) and converged considerably better than the
NLME methods.
41
Limitations of our work include the following three issues. First, the current
implementation of U-HB in JAGS is computationally intensive and requires hours to fit
instead of seconds (TS-NLS, TS-HMA), minutes (TS-NLME) or tens of minutes (U-
NLME). To overcome this challenge in the exploratory or interactive model-building
stage of an applied data analysis, an analyst may wish to use a two-stage version of HB.
One can run U-HB with no covariates to obtain participant-level NO parameters in Stage
1, and then in Stage 2, run the usual separate linear regression models relating NO
parameters to covariates. Penultimate and final models could then be estimated with U-
HB. Running the final analyses using U-HB, over the course of several hours, is a
minimal commitment of resources compared to the months or years that most studies
have devoted to planning and data collection. Future work may consider alternative
MCMC methods, with the goal of reducing run time. Second, our simulation studies
though extensive, did not cover all possible scenarios. Here we generated data using
magnitudes of association and NO parameter distributions similar to those observed in
the CHS data, under a single covariate scenario. However, U-HB can be readily applied
to scenarios with multiple covariates or even interactions, with some additional
computational cost. For example using CHS data, U-HB converged after 4.8 hours for the
model with only one covariate, but 9.2 hours for a model with two covariates and 9.5
hours for a model with four covariates. Third, we used the conservative approach of only
comparing methods using datasets for which all the methods converged. If we included
all datasets (dropping a dataset for a method only if that method failed to converge on
that dataset), the relative performance the methods was similar (Supplemental Figure 5).
42
In summary, we have performed the first evaluation of statistical methods used to
estimate the associations of covariates with NO parameters. We found the best
performance using a novel unified estimation method (U-HB), which requires longer
computation time (hours rather than minutes) but which can be readily implemented
using standard statistical software.
Software
An example dataset and R code for implementing all 5 methods are available on Github
https://github.com/USCbiostats/FeNO_UHB
43
References
1. Braun PX, Gmachl CF, Dweik RA. Bridging the Collaborative Gap: Realizing the
Clinical Potential of Breath Analysis for Disease Diagnosis and Monitoring–Tutorial.
Sensors Journal, IEEE. 2012;12(11):3258-3270.
2. Dweik RA, Boggs PB, Erzurum SC, et al. An official ATS clinical practice
guideline: interpretation of exhaled nitric oxide levels (FENO) for clinical applications.
Am J Respir Crit Care Med. 2012;184(5).
3. La Grutta S, Ferrante G, Malizia V, Cibella F, Viegi G. Environmental effects on
fractional exhaled nitric oxide in allergic children. J Allergy (Cairo). 2012;2012:916926.
4. Scarpa MC, Kulkarni N, Maestrelli P. The role of non-invasive biomarkers in
detecting acute respiratory effects of traffic-related air pollution. Clin Exp Allergy.
2014;44(9):1100-1118.
5. Hogman M, Stromberg S, Schedin U, Frostell C, Hedenstierna G, Gustafsson LE.
Nitric oxide from the human respiratory tract efficiently quantified by standardized single
breath measurements. Acta Physiol Scand. 1997;159(4):345-346.
6. Silkoff PE, McClean PA, Slutsky AS, et al. Marked flow-dependence of exhaled
nitric oxide using a new technique to exclude nasal nitric oxide. Am J Respir Crit Care
Med. 1997;155(1):260-267.
7. George SC. How accurately should we estimate the anatomical source of exhaled
nitric oxide? J Appl Physiol. 2008;104(4):909-911.
8. Eckel SP, Salam MT. Single high flow exhaled nitric oxide is an imperfect proxy for
distal nitric oxide. Occup Environ Med. 2013;70(7):519-520.
9. George SC, Hogman M, Permutt S, Silkoff PE. Modeling pulmonary nitric oxide
exchange. J Appl Physiol. 2004;96(3):831-839.
10. ATS. Recommendations for standardized procedures for the on-line and off-line
measurement of exhaled lower respiratory nitric oxide and nasal nitric oxide in adults and
children-1999. This official statement of the American Thoracic Society was adopted by
the ATS Board of Directors, July 1999. Am J Respir Crit Care Med. 1999;160(6):2104-
2117.
11. ATS/ERS. ATS/ERS recommendations for standardized procedures for the online
and offline measurement of exhaled lower respiratory nitric oxide and nasal nitric oxide,
2005. Am J Respir Crit Care Med. 2005;171(8):912-930.
12. Horvath I, Barnes PJ, Loukides S, et al. A European Respiratory Society technical
standard: exhaled biomarkers in lung disease. Eur Respir J. 2017;49(4).
13. Hogman M, Merilainen P. Extended NO analysis in asthma. J Breath Res.
2007;1(2):024001.
44
14. Hogman M. Extended NO analysis in health and disease. J Breath Res.
2012;6(4):047103.
15. Modig L, Dahgam S, Wass K, Nyberg F, Olin AC. Effects of short-term exposure to
air pollution on the levels of exhaled nitric oxide among adults – Results from the
ADONIX study in Gothenburg, Sweden. Eur Respir J. 2012;40(s56):609s.
16. Rosa MJ, Perzanowski MS, Divjan A, et al. Association of recent exposure to
ambient metals on fractional exhaled nitric oxide in 9–11year old inner-city children.
Nitric Oxide. 2014;40:60-66.
17. Eckel S, Zhang Z, Habre R, et al. Traffic-related air pollution and alveolar nitric
oxide in southern California children. Eur Respir J. 2016;47(5):1348-1356.
18. Lopez V, Prieto L, Perez-Frances C, Barato D, Marin J. Natural pollen exposure
increases the response plateau to adenosine 5'-monophosphate and bronchial but not
alveolar nitric oxide in sensitized subjects. Respiration. 2011;83(3):225-232.
19. Eckel SP, Linn WS, Berhane K, et al. Estimation of parameters in the two-
compartment model for exhaled nitric oxide. PLoS ONE. 2014;9(1):e85471.
20. Eckel S, Berhane K, Linn W, et al. Statistical Methods For Studying Determinants
Of Airway And Alveolar Nitric Oxide. Am J Respir Crit Care Med. 2013;187:A5092.
21. Davidian M, Giltinan DM. Nonlinear models for repeated measurement data: an
overview and update. J Agric Biol Envir S. 2003;8(4):387-419.
22. Casella G, George EI. Explaining the Gibbs sampler. The American Statistician.
1992;46(3):167-174.
23. Gauderman WJ, Avol E, Gilliland F, et al. The effect of air pollution on lung
development from 10 to 18 years of age. N Engl J Med. 2004;351(11):1057-1067.
24. McConnell R, Berhane K, Yao L, et al. Traffic, susceptibility, and childhood asthma.
Environ Health Perspect. 2006;114(5):766-772.
25. Chen Z, Salam MT, Eckel SP, Breton CV, Gilliland FD. Chronic effects of air
pollution on respiratory health in Southern California children: findings from the
Southern California Children’s Health Study. J Thorac Dis. 2015;7(1):46.
26. Salam MT, Islam T, Gilliland FD. Recent evidence for adverse effects of residential
proximity to traffic sources on asthma. Curr Opin Pulm Med. 2008;14(1):3-8.
27. Health Effects Institute. Panel on the Health Effects of Traffic-Related Air Pollution.
Traffic-related air pollution: a critical review of the literature on emissions, exposure,
and health effects. Health Effects Institute; 2010.
28. Brook RD, Rajagopalan S, Pope CA, et al. Particulate matter air pollution and
cardiovascular disease. Circulation. 2010;121(21):2331-2378.
29. Benbrahim-Tallaa L, Baan RA, Grosse Y, et al. Carcinogenicity of diesel-engine and
gasoline-engine exhausts and some nitroarenes. Lancet Oncol. 2012;13(7):663-664.
45
30. Anderson HR, Favarato G, Atkinson RW. Long-term exposure to air pollution and
the incidence of asthma: meta-analysis of cohort studies. Air Qual Atmos Health.
2013;6(1):47-56.
31. Kelly FJ, Fussell JC. Air pollution and public health: emerging hazards and
improved understanding of risk. Environ Geochem Health. 2015;37(4):631-649.
32. Schultz ES, Litonjua AA, Melé n E. Effects of long-term exposure to traffic-related
air pollution on lung function in children. Curr Allergy Asthma Rep. 2017;17(6):41.
33. Tsoukias NM, George SC. A two-compartment model of pulmonary nitric oxide
exchange dynamics. J Appl Physiol. 1998;85(2):653-666.
34. Hogman M, Holmkvist T, Wegener T, et al. Extended NO analysis applied to
patients with COPD, allergic asthma and allergic rhinitis. Respir Med. 2002;96(1):24-30.
35. Silkoff PE, Sylvester JT, Zamel N, Permutt S. Airway nitric oxide diffusion in
asthma: Role in pulmonary function and bronchial responsiveness. Am J Respir Crit Care
Med. 2000;161(4 Pt 1):1218-1228.
36. Molshatski N, Eckel SP. Optimal flow rate sampling designs for studies with
extended exhaled nitric oxide analysis. J Breath Res. 2017;11(1):016012.
37. Gelman A, Bois F, Jiang J. Physiological pharmacokinetic analysis using population
modeling and informative prior distributions. Journal of the American Statistical
Association. 1996;91(436):1400-1412.
38. Eckel SP, Zhang Z, Habre R, et al. Traffic-related air pollution and alveolar nitric
oxide in southern California children. European Respiratory Journal. 2016.
39. Plummer M. JAGS: A program for analysis of Bayesian graphical models using
Gibbs sampling. Paper presented at: Proceedings of the 3rd international workshop on
distributed statistical computing2003.
40. Brooks SP, Gelman A. General methods for monitoring convergence of iterative
simulations. Journal of computational and graphical statistics. 1998;7(4):434-455.
41. Tsoukias NM, Tannous Z, Wilson AF, George SC. Single-exhalation profiles of NO
and CO2 in humans: effect of dynamically changing flow rate. J Appl Physiol.
1998;85(2):642-652.
46
Supplement:
Table of Contents
Section 1: Convergence failure rates in the simulation study ............................................47
Section 2: Bias, Coverage, Power and Confidence Interval Length figures for 𝛽 , by
scenario ..............................................................................................................................49
Section 3: Bias, Coverage, Power and Confidence Interval Length figures for 𝛼 , by
scenario. .............................................................................................................................64
Section 4. Additional details on CHS results .....................................................................77
Section 5. Simulation study results using all available estimates (not limited to datasets
for which all models converged). .......................................................................................84
47
Section 1: Convergence failure rates in the simulation study
Supplementary Figure 1: Rates of convergence failure in the simulated datasets, by
method (denoted by line type) by scenario (denoted by line color, see paper for
definitions), as a function of effect size, using simulated datasets.
There was no trend in converge rate failures by Scenario or effect size, so we summarized
converge failure rates as the average across Scenarios and effect sizes. There were
notable difference in convergence failure rates across methods: 9% for TS-NLME, 30%
for U-NLME and only 2% for U-HB.
48
49
Section 2: Bias, Coverage, Power and Confidence Interval Length figures for 𝛽 , by
scenario
In this section, we provide full results for our simulation study for all seven Scenarios.
Results for Scenario 1 were included in the main body of the manuscript, but we repeat
them below for completeness.
Supplementary Figure 2.1. Relative bias (a), coverage (b), and power (c) CI length (d) of
the selected estimation methods from Scenario 1 of the simulation study (𝛽 𝐶 𝐴 = 𝛽 𝑙𝑜𝑔𝐶 𝑎𝑤
=
𝛽 𝑙𝑜𝑔𝐷 𝑎𝑤
).
50
We see that 𝛽 𝐶 𝐴 was underestimated by all models except U-HB, with the degree of
negative bias increasing with effect size. U-HB and U-NLME showed least bias for
𝛽 𝑙𝑜𝑔𝐶 𝑎𝑤
and 𝛽 𝑙𝑜𝑔𝐷 𝑎𝑤
but in different directions. TS-NLS, TS-HMA and TS-NLME
showed good power but it is because they have biased estimation, short confidence
interval length, or low converage. TS-NLME performed well for estimating 𝛽 𝑙𝑜𝑔 𝐶 𝑎𝑤
but
51
not for 𝛽 𝐶 𝐴 and 𝛽 𝑙𝑜𝑔𝐷 𝑎𝑤
. U-NLME performed well for 𝛽 𝑙𝑜𝑔𝐶 𝑎𝑤
and 𝛽 𝑙𝑜𝑔𝐷 𝑎𝑤
but not for
𝛽 𝐶 𝐴 .
52
Supplementary Figure 2.2. Relative bias (a), coverage (b), and power (c) CI length (d) of
the selected estimation methods from Scenario 2 of the simulation study (𝛽 𝐶 𝐴 varied,
𝛽 𝑙𝑜𝑔𝐶 𝑎𝑤
= 0, and 𝛽 𝑙𝑜𝑔𝐷 𝑎𝑤
= 0).
53
The directions of bias were changed in this scenario for 𝛽 𝑙𝑜𝑔𝐶 𝑎𝑤
and 𝛽 𝑙𝑜𝑔𝐷 𝑎𝑤
. U-HB had
good performance and was very stable, unlike U-NLME, which had high type I error due
to its bias and narrow CI length.
54
Supplementary Figure 2.3. Relative bias (a), coverage (b), and power (c) CI length (d) of
the selected estimation methods from Scenario 3 of the simulation study (𝛽 𝐶 𝐴 = 0,
𝛽 𝑙𝑜𝑔𝐶 𝑎𝑤
varied, and 𝛽 𝑙𝑜𝑔𝐷 𝑎𝑤
= 0).
55
The direction of bias flipped again for 𝛽 𝑙𝑜𝑔𝐶 𝑎𝑤
and 𝛽 𝑙𝑜𝑔𝐷 𝑎𝑤
. U-HB was still the best and
most stable model.
56
Supplementary Figure 2.4. Relative bias (a), coverage (b), and power (c) CI length (d) of
the selected estimation methods from Scenario 4 of the simulation study (𝛽 𝐶 𝐴 = 0,
𝛽 𝑙𝑜𝑔𝐶 𝑎𝑤
= 0, and 𝛽 𝑙𝑜𝑔𝐷 𝑎𝑤
varied).
57
Estimation of 𝛽 𝑙𝑜𝑔𝐶 𝑎𝑤
was affected greatly by the increasing size of 𝛽 𝑙𝑜𝑔𝐷 𝑎 𝑤 for all
models except U-HB. While estimation of 𝛽 𝐶 𝐴 was hardly affected.
58
Supplementary Figure 2.5. Relative bias (a), coverage (b), and power (c) CI length (d) of
the selected estimation methods from Scenario 5 of the simulation study (𝛽 𝐶 𝐴 = 0,
𝛽 𝑙𝑜𝑔𝐶 𝑎𝑤
= 𝛽 𝑙𝑜𝑔𝐷 𝑎𝑤
varied).
59
60
Overall, the U-HB model again had the best performance.
Supplementary Figure 2.6. Relative bias (a), coverage (b), and power (c) CI length (d) of
the selected estimation methods from Scenario 6 of the simulation study (𝛽 𝑙𝑜𝑔𝐶 𝑎𝑤
= 0
while 𝛽 𝐶 𝐴 = 𝛽 𝑙𝑜𝑔𝐷 𝑎𝑤
varied).
61
62
Overall, the U-HB model again had the best performance.
63
Supplementary Figure 2.7. Relative bias (a), coverage (b), and power (c) CI length (d) of
the selected estimation methods from Scenario 7 of the simulation study (𝛽 𝑙𝑜𝑔𝐷 𝑎𝑤
= 0
while 𝛽 𝐶 𝐴 = 𝛽 𝑙𝑜𝑔𝐶 𝑎𝑤
varied).
64
Overall, the U-HB model again had the best performance.
Section 3: Bias, Coverage, Power and Confidence Interval Length figures for 𝛼 , by
scenario
Supplementary Figure 3.1. Relative bias (a), coverage (b), and power (c) CI length (d) of
the selected estimation methods from Scenario 1 of the simulation study (𝛽 𝐶 𝐴 = 𝛽 𝑙𝑜𝑔𝐶 𝑎𝑤
=
𝛽 𝑙𝑜𝑔𝐷 𝑎𝑤
).
65
66
Supplementary Figure 3.2. Relative bias (a), coverage (b), and power (c) CI length (d) of
the selected estimation methods from Scenario 2 of the simulation study (𝛽 𝐶 𝐴 varied,
𝛽 𝑙𝑜𝑔𝐶 𝑎𝑤
= 0, and 𝛽 𝑙𝑜𝑔𝐷 𝑎𝑤
= 0).
67
68
69
Supplementary Figure 3.3. Relative bias (a), coverage (b), and power (c) CI length (d) of
the selected estimation methods from Scenario 3 of the simulation study (𝛽 𝐶 𝐴 = 0,
𝛽 𝑙𝑜𝑔𝐶 𝑎𝑤
varied, and 𝛽 𝑙𝑜𝑔𝐷 𝑎𝑤
= 0).
70
71
Supplementary Figure 3.4. Relative bias (a), coverage (b), and power (c) CI length (d) of
the selected estimation methods from Scenario 4 of the simulation study (𝛽 𝐶 𝐴 = 0,
𝛽 𝑙𝑜𝑔𝐶 𝑎𝑤
= 0, and 𝛽 𝑙𝑜𝑔𝐷 𝑎𝑤
varied).
72
73
Supplementary Figure 3.5. Relative bias (a), coverage (b), and power (c) CI length (d) of
the selected estimation methods from Scenario 5 of the simulation study (𝛽 𝐶 𝐴 = 0,
𝛽 𝑙𝑜𝑔𝐶 𝑎𝑤
= 𝛽 𝑙𝑜𝑔𝐷 𝑎𝑤
varied).
74
75
Supplementary Figure 3.6. Relative bias (a), coverage (b), and power (c) CI length (d) of
the selected estimation methods from Scenario 6 of the simulation study (𝛽 𝑙𝑜𝑔𝐶 𝑎𝑤
= 0
while 𝛽 𝐶 𝐴 = 𝛽 𝑙𝑜𝑔𝐷 𝑎𝑤
varied).
76
Supplementary Figure 3.7. Relative bias (a), coverage (b), and power (c) CI length (d) of
the selected estimation methods from Scenario 7 of the simulation study (𝛽 𝑙𝑜𝑔𝐷 𝑎𝑤
= 0
while 𝛽 𝐶 𝐴 = 𝛽 𝑙𝑜𝑔𝐶 𝑎𝑤
varied).
77
Section 4. Additional details on CHS results
78
For CHS participants, the estimated participant-level NO parameters were highly
correlated (Pearson’s R>0.92) across U-HB, U-NLME and TS-NLME methods under the
Caw parameterization (Supplementary Figure 4.1) and also under the J’aw
parameterization (Supplementary Figure 4.2). As shown in Supplementary Figure 4.3,
calculated logJaw estimates from U-HB (Caw parameterization) were extremely highly
correlated (Pearson’s R>0.999) with direct logJaw estimates from U-HB (Jaw-
parameterization).
Supplementary Figure 4.1. Correlation of estimated participant-level NO parameters in
the CHS across estimation methods using Caw parameterization: (left) CA, (middle)
logCaw, (right) logDaw.
79
Supplementary Figure 4.2. Correlation of estimated participant-level NO parameters in
the CHS across estimation methods using J’aw parameterization: (left) CA, (middle)
logJ’aw, (right) logDaw.
80
Supplementary Figure 4.3. Comparison of estimated participant-level logJ’aw in the CHS
from J’aw and Caw parameterizations of U-HB.
81
Supplementary Figure 4.4. Analysis of CHS data: estimated population-level mean NO
parameters when traffic-related air pollution was 0 (𝛼 ̂ and 95% CI) using the selected
methods,
*
with no adjustments for covariates.
82
*
TS-NLME(J’ aw) is the previously published model using a J’ aw parameterization of TS-NLME.
83
Supplementary Table 4.1. Analysis of CHS data: estimated population-level parameters
(and 95% CI) using the selected methods,
*
in models that include traffic-related air
pollution, age, sex, and asthma.
Estimated
parameter
(95% CI)
TS-NLS TS-HMA TS-NLME(J’ aw) TS-NLME U-NLME U-HB
𝛼 C
𝐴 1.48 (1.32, 1.64) 1.58 (1.44, 1.71) 1.30 (1.24,1.35) 1.29 (1.24,1.34) 1.25 (1.17,1.33) 1.51 (1.43,1.60)
𝛽 traffic
C
𝐴
0.23 (0.06, 0.40) 0.32 (0.18, 0.46) 0.08 (0.03,0.14) 0.08 (0.03,0.14) 0.13 (0.05,0.22) 0.15 (0.66,0.24)
𝛽 age
C
𝐴
0.52 (0.28, 0.76) 0.65 (0.43, 0.87) -0.01 (-0.05,0.03) -0.01 (-0.04,0.02) -0.03 (-0.09,0.02) -0.02 (-0.07,0.04)
𝛽 male
C
𝐴
0.17 (-0.03, 0.37) 0.23 (0.05, 0.40) 0.10 (0.03,0.17) 0.10 (0.03,0.17) 0.14 (0.03,0.25) 0.17 (0.07,0.28)
𝛽 asthma
C
𝐴
-0.001 (-0.10, 0.10) -0.03 (-0.12, 0.06) 0.24 (0.15,0.33) 0.24 (0.15,0.32) 0.29 (0.14,0.43) 0.41 (0.26,0.55)
𝛼 logC
𝑎𝑤
3.92 (3.84, 4.00) 4.39 (4.25, 4.54) 3.71 (3.65,3.76 3.68 (3.63,3.74) 3.70 (3.62,3.77) 4.13 (4.02,4.23)
𝛽 traffic
logC
𝑎𝑤
0.04 (-0.05, 0.13) 0.16 (0.02, 0.31) 0.02 (-0.04,0.07) 0.01 (-0.04,0.07) -0.03 (-0.09,0.06) -0.02 (-0.12,0.07)
𝛽 age
logC
𝑎𝑤
0.46 (0.34, 0.59) 0.42 (0.20, 0.65) 0.05 (0.01,0.08) 0.05 (0.02,0.08) 0.06 (0.01,0.11) 0.07 (0.01,0.13)
𝛽 male
logC
𝑎𝑤
0.09 (-0.01, 0.20) 0.08 (-0.11, 0.26) 0.12 (0.05,0.19) 0.13 (0.05,0.19) 0.14 (0.04,0.23) 0.13 (0.02,0.25)
𝛽 asthma
logC
𝑎𝑤
0.07 (0.02, 0.13) 0.02 (-0.07, 0.11) 0.36 (0.28,0.44) 0.37 (0.28,0.45) 0.38 (0.27,0.50) 0.33 (0.19,0.48)
𝛼 logD
𝑎𝑤
2.72 (2.65, 2.79) 2.13 (1.96, 2.29) 2.89 (2.84,2.93) 2.93 (2.87,2.96) 2.91 (2.82,2.99) 2.39 (2.27,2.51)
𝛽 traffic
logD
𝑎𝑤
0.03 (-0.05, 0.10) -0.14 (-0.31, 0.02) 0.03 (-0.01,0.07) 0.03 (-0.01,0.08) 0.07 (-0.01,0.15) 0.07 (-0.04,0.18)
𝛽 age
logD
𝑎𝑤
-0.03 (-0.05, 0.10) 0.06 (-0.20,0.32) 0.01 (-0.02,0.04) 0.01 (-0.02,0.03) 0.000 (-0.05,0.05) -0.01 (-0.08,0.05)
𝛽 male
logD
𝑎𝑤
0.09 (0.01, 0.18) 0.07 (-0.14, 0.28) 0.02 (-0.03,0.08) 0.03 (-0.03,0.08) -0.001 (-0.10,0.10) 0.01 (-0.12,0.14)
𝛽 asthma
logD
𝑎𝑤
-0.03 (-0.07, 0.02) 0.05 (-0.05,0.15) 0.12 (0.06,0.19) 0.11 (0.04,0.17) 0.09 (-0.04,0.21) 0.15 (-0.01,0.31)
*
TS-NLME(J’ aw) is the previously published model using a J’ aw parameterization of TS-NLME.
84
Section 5. Simulation study results using all available estimates (not limited to datasets
for which all methods converged).
Supplementary Figure 5: Relative bias (a), coverage (b), and power (c) CI length (d) for
𝛽 ′s using all available estimates (not just the subset of estimates from datasets for which
all methods converged) for the selected estimation methods in simulation study Scenario
1 (𝛽 𝐶 𝐴 = 𝛽 𝑙𝑜𝑔𝐶 𝑎𝑤
= 𝛽 𝑙𝑜𝑔𝐷 𝑎𝑤
).
85
Previously presented results used the subset of simulated datasets on which all methods
converged. Since convergence failure rates varied across methods (e.g., ~30% for U-
NLME), we present Supplementary Figure 5 including results from simulation Scenario 1
on all available estimates (from all datasets rather than the subset of datasets where all
methods converged) to evaluate the impact of this choice on our conclusions.
Conclusions for Scenario 1 from Supplementary Figure 5 (all available estimates) are
similar to those in Figure 2 (subset of estimates). For other scenarios, conclusions are also
same (data not shown).
86
Chapter 3: Longitudinal Hierarchical Bayesian models of
covariate effects on airway and alveolar nitric oxide
Abstract
Biomarkers such as exhaled nitric oxide (FeNO), a marker of airway inflammation, have
applications in the study of chronic respiratory disease where longitudinal studies of
within-participant changes in the biomarker are particularly relevant. A cutting-edge
approach to assessing FeNO, called multiple flow FeNO, repeatedly assesses FeNO
across a range of expiratory flow rates at a single visit and combines these data with a
deterministic model of lower respiratory tract NO to estimate parameters quantifying
airway wall and alveolar NO sources. Previous methodological work for multiple flow
FeNO has focused on methods for data from a single participant or from cross-sectional
studies. Performance of existing ad hoc two-stage methods for longitudinal multiple flow
FeNO in cohort or panel studies has not been evaluated. In this paper, we present a novel
longitudinal extension to a unified hierarchical Bayesian (L-U-HB) model relating
longitudinally assessed multiple flow FeNO to covariates. In several simulation study
scenarios, we compare the L-U-HB method to other unified and two-stage frequentist
methods. In general, L-U-HB produced unbiased estimates, had good power, and its
performance was not sensitive to the magnitude of the association with a covariate and
correlations between NO parameters. In an application relating height to longitudinal
multiple flow FeNO in schoolchildren without asthma, unified analysis methods
estimated positive, statistically significant associations of height with airway and alveolar
87
NO concentrations and negative associations with airway wall diffusivity while estimates
from two-stage methods were smaller in magnitude and sometimes non-significant.
88
Introduction
Sophisticated statistical methods are needed to link longitudinal assessments of a
biomarker to patient-level characteristics and time-varying exposures in the context of a
deterministic mathematical model describing the production and dynamics of the
biomarker within the human body. This paper presents statistical methods developed for
such longitudinal assessments of the fractional concentration of exhaled nitric oxide,
FeNO, a biomarker of airway inflammation used in clinical [45-47] and epidemiological
research[48-50].
FeNO is an exhaled breath biomarker conventionally assessed at the target expiratory
flow rate of 50 ml/s (FeNO50)[51]. A cutting-edge approach, called multiple flow FeNO,
repeatedly assesses FeNO across a range of expiratory flow rates and combines these data
with a deterministic model of NO in the lower respiratory tract to estimate parameters
quantifying the effect of airway wall and alveolar sources. Literature on the modeling of
multiple flow FeNO data has focused on methods for data collected from one person at a
single visit[15, 16, 52, 53]. Multiple flow FeNO data present statistical challenges which
require sophisticated statistical methods. Many studies of multiple flow FeNO conduct
analyses using a two-stage approach: (1) estimate airway and alveolar NO parameters and
(2) treat the estimated NO parameters as observed outcomes in linear regressions relating
NO parameters to factors of interest (asthma medication use, air pollution exposures,
etc.). In a previous paper, we presented a novel unified hierarchical Bayesian (U-HB)
model for estimating cross-sectional associations of covariates with NO parameters using
data from a single multiple flow FeNO test session for each study participant. We also
89
found, in an extensive simulation study, that the U-HB method was less biased and had
better power/type I error compared to conventional two-stage methods[54].
There is a need for longitudinal data analysis methods for FeNO. Biomarkers like FeNO
are particularly promising for tracking within-person changes over time since they tend to
be relatively stable within persons, despite considerable heterogeneity across people[55-
58]. Longitudinal trends in study populations with repeated measures of the conventional
FeNO50 are generally modeled using standard longitudinal data analysis techniques, such
as linear mixed effects models (LMM) or generalized additive linear mixed effects
models (GAMM) [59]. Panel studies or longitudinal cohort studies with repeated
measures of multiple flow FeNO across multiple visits can highlight within-person trends
in proximal and distal inflammation. However, there have been no methods proposed for
longitudinal multiple flow FeNO data, and the performance of existing ad hoc two-stage
methods has not been evaluated.
Here, we present a novel extension of the U-HB model to longitudinal data (L-U-HB)
using a Bayesian implementation of nonlinear mixed effects models. L-U-HB takes as
input longitudinal measurements of multiple-flow FeNO on a group of participants to
estimate associations of NO parameters with time-varying or time-constant covariates, as
well as to quantify within- and between-participant variation. Our work is motivated by
longitudinal multiple flow FeNO data collected as part of the Southern California
Children’s Health Study (CHS) [60]. The CHS, originally designed to study impacts of
long-term air pollution exposures on children’s respiratory health, included repeated
measurements of multiple flow FeNO in the most recent cohort.
90
Methods
First, we introduce the mathematical deterministic model for FeNO in the respiratory
tract and then we describe statistical methods for estimating NO parameters from this
model using multiple flow FeNO data, including two-stage (TS) approaches and unified
(U) approaches for longitudinally assessed multiple flow FeNO.
Deterministic two compartment model for FeNO
Our work is based on the simple steady-state two-compartment model (2CM), which
assumes a cylindrically-shaped airway compartment with related NO parameters: Caw, the
concentration of NO in the airway tissue (ppb); Daw, the airway tissue diffusion capacity
(pL· s
-1
· ppb
-1
) , and an expansile alveolar compartment with related NO parameter CA, the
concentration of NO in the alveolar region (ppb) [15]. Under the 2CM, FeNO (ppb) at the
mouth is deterministically related to expiratory flow rate (ml/s) and the three NO
parameters quantifying airway and alveolar sources of NO, as shown below:
𝐹𝑒𝑁𝑂 = 𝐶 𝑎𝑤
+ ( 𝐶 𝐴 − 𝐶 𝑎𝑤
)× 𝑒 −
𝐷 𝑎𝑤
𝑓𝑙𝑜𝑤 (1)
Estimating NO parameters in the 2CM
The 2CM model for FeNO in Equation 1 is deterministic and nonlinear. In practice,
multiple flow FeNO data is measured with error. Researchers have developed various
methods to estimate 2CM NO parameters using multiple flow data from a given
participant, typically using linear regression approaches with an underlying linearization
assumption, a third order approximation method such as the Högman and Merilӓinen
91
algorithm (HMA) [61, 62], and nonlinear regression which essentially adds an error term
to the right hand side of Equation 1 [17]. Here and in previous work[16, 54], we use the
following fundamental nonlinear statistical model for multiple flow FeNO measured
repeatedly across a range of flow rates for a single participant, with maneuvers indexed
by k:
𝑙𝑜𝑔 ( 𝐹𝑒𝑁𝑂 𝑘 ) = 𝑙𝑜𝑔 ( 𝐶 𝑎𝑤
+ ( 𝐶 𝐴 − 𝐶 𝑎𝑤
)× 𝑒 −
𝐷 𝑎𝑤
𝑓𝑙𝑜𝑤 𝑘 )+ 𝜀 𝑘 (2)
This model formulation includes a “transform-both-sides” [63] approach using the natural
log to acknowledge the increased variation in error that occurs as flow rate (and hence
FeNO concentration) increases w hile maintaining the interpretability of the 2CM NO
parameters. On the logFeNO scale, the error (ε) can be reasonably assumed to be
normally distributed. Henceforth, we will refer to a model estimating NO parameters
using Equation 2 with standard nonlinear-least squares software (e.g., “nls” from the
nlme package in R) as NLS [16]. So far, we have discussed only estimation of NO
parameters from a single multiple flow FeNO test session for one participant. When
multiple flow FeNO data are assessed longitudinally in a study population the data have
three levels of variation: across-participant, within-participant (across visits), and within-
visit (across maneuvers).
Estimating associations of covariates with longitudinally assessed NO parameters
Two Stage (TS) methods
92
In the existing literature, most researchers use ad hoc two stage approaches to relate
estimated NO parameters to covariates. Two stage methods for cross-sectional studies
were discussed in our previous work[54]. A typical longitudinal two stage method
proceeds as follows. In Stage I, NO parameters for each participant at each visit are
estimated via separate models. For example, a separate HMA or NLS model is fit to the
multiple flow FeNO data from each participant at each visit. In Stage II, three linear
mixed effect models (LMMs) are fit, one for each NO parameter, to relate the
longitudinal estimates of the NO parameters to a covariate(s) of interest, denoted
generically as Xij. A participant-level random intercept is included in each LMM to
account for the within-participant correlation in the longitudinal NO parameter data.
Below, we introduce 4 longitudinal two stage (L-TS) methods, differentiated by the name
of the method employed in Stage I:
93
1. L-TS-NLS: Stage I consists of N (participants) x M (visits) separate NLS [16]
models, each NLS model fit to the typically small multiple flow FeNO dataset at
that visit.
2. L-TS-HMA: Similarly, Stage I consists of N x M separate HMA [61, 62] models.
3. L-TS-NLME: Stage I consists of a single longitudinal nonlinear least square
mixed effect (NLME) model, an extension of the approach using N x M separate
NLS models. In the longitudinal NLME, we specified participant-level and visit-
level random intercepts for each NO parameter. At each level, these random
effects follow a multivariate normal distribution, allowing for correlation of NO
parameters. For example, participant-level correlation allows for participants with
high CA to also tend to have high Caw. We implemented NLME using the nlme
package in R (version 3.1-152) [64].
4. L-TS-HB: Stage I consists of a single longitudinal Hierarchical Bayesian (HB)
analog of the longitudinal NLME model, implemented using JAGS (Just another
Gibbs sampler [65]) similar to the U-HB cross-sectional model in our previous
publication[54], but with no covariate X and a partitioning of variance in NO
parameters at the participant- and visit-levels through specification of variance-
covariance matrices, the population mean of the NO parameters, and the
measurement error in Stage I. This model is described in greater detail in the L-U-
HB section below, where the model includes X.
All these TS approaches use the same LMM approach in Stage II.
Unified approaches
Unified methods, in contrast to TS methods, simultaneously estimate NO parameters and
their associations with the covariate 𝑋 𝑖𝑗
in a single model. In this longitudinal version,
we estimated the between/within-participant variation at the same time, which L-TS-NLS
94
and L-TS-HMA were not able to obtain because their estimation in the stage I ignored the
grouping effect.
5. L-U-HB
Figure 7: Diagram of the hierarchical model structure relating Longitudinal FeNO
measurements at multiple flow rates to NO parameters that are a function of a potential
determinant X (e.g., air pollution)
Our novel U-HB model for longitudinal data (U-HBL) has three levels: maneuver, visit,
participants, as described below and displayed in Figure 7:
Level 1: Maneuver
𝑙𝑜𝑔 ( 𝐹𝑒𝑁𝑂 𝑖𝑗𝑘 )∼ 𝑁 ( 𝑓 ( 𝜃 𝑖𝑗
, 𝑓𝑙𝑜𝑤 𝑖𝑗𝑘 ) , 𝜎 𝜖 2
) (3)
𝑓 ( 𝜃 𝑖𝑗
, 𝑓𝑙𝑜𝑤 𝑖𝑗𝑘 ) = 𝑙𝑜𝑔 (𝑒𝑥𝑝 ( 𝑙𝑜𝑔𝐶 𝑎 𝑤 𝑖𝑗
)+ ( 𝐶 𝐴 𝑖𝑗
− 𝑒𝑥𝑝 ( 𝑙𝑜𝑔𝐶 𝑎𝑤
𝑖𝑗
) )× 𝑒 −
𝑒𝑥𝑝 ( 𝑙𝑜𝑔𝐷 𝑎𝑤
𝑖𝑗
)
𝑓𝑙𝑜𝑤 𝑖𝑗𝑘 )
(4)
95
In the first level, log FeNO for participant i at visit j and maneuver k is assumed to be
normally distributed with a mean that is a function of NO parameters: 𝜃 𝑖𝑗
=
( 𝐶 𝐴 𝑖𝑗
, 𝑙𝑜𝑔𝐶 𝑎𝑤
𝑖𝑗
, 𝑙𝑜𝑔𝐷 𝑎𝑤
𝑖𝑗
) ′ and expiratory flow, 𝑓𝑙𝑜𝑤 𝑖𝑗𝑘 . The variance of the
unexplained error in logFeNO, 𝜎 𝜖 2
, was assumed to be the same across flow rates, visits,
and participants.
Level 2: Visit (time)
𝜃 𝑖𝑗
= 𝐴 𝛩 𝑖 + 𝛽 1
𝛩 𝑋 𝑖𝑗
+ 𝛼 𝛩 𝑖𝑗
(5)
In the second level, NO parameters for participant i at visit j ( 𝜃 𝑖𝑗
) are modeled as a
linear function of 𝐴 𝛩 𝑖 , a vector of participant-level mean NO parameter values for
participant i when the covariate 𝑋 𝑖𝑗
= 0. Key parameters of interest include 𝛽 𝛩 =
( 𝛽 𝐶𝑎
, 𝛽 𝑙𝑜𝑔𝐶𝑎𝑤 , 𝛽 𝑙𝑜𝑔𝐷𝑎𝑤 ) ′, the regression coefficients on 𝑋 𝑖𝑗
. Otherwise unexplained
within-participant variation in the NO parameters is represented by the visit-level random
intercepts 𝛼 𝛩 𝑖𝑗
= ( 𝛼 𝐶𝑎
𝑖𝑗
, 𝛼 𝑙𝑜𝑔𝐶𝑎𝑤 𝑖𝑗
, 𝛼 𝑙𝑜𝑔𝐷𝑎𝑤 𝑖𝑗
) assumed to have no correlations and
follow a multivariate normal distribution (MVN) with variance-covariance matrix 𝛴 𝜎 =
𝑑𝑖𝑎𝑔 ( 𝜎 𝐶 𝐴 𝑖𝑗
, 𝜎 𝑙𝑜𝑔 𝐶𝑎𝑤 𝑖𝑗
, 𝜎 𝑙𝑜𝑔 𝐷𝑎𝑤 𝑖𝑗
)
, i.e., 𝛼 𝛩 𝑖𝑗
~MVN (0, 𝛴 𝜏 ) .
Level 3: Participant
𝐴 𝛩 𝑖 = 𝛽 0
𝛩 + 𝛼 𝛩 𝑖 (6)
In the third level, 𝐴 𝛩 𝑖 is decomposed into 𝛽 0
𝛩 , the overall population-mean NO
parameters when 𝑋 𝑖𝑗
= 0 and the otherwise unexplained between-participant variation in
96
the NO parameters, represented by the participant-level random intercepts 𝛼 𝛩 𝑖 =
( 𝛼 𝐶𝑎
𝑖 , 𝛼 𝑙𝑜𝑔𝐶𝑎𝑤 𝑖 , 𝛼 𝑙𝑜𝑔𝐷𝑎𝑤 𝑖 ) ′, assumed to follow a MVN with variance-covariance matrix
𝛴 𝜏 , i.e., 𝛼 𝛩 𝑖 ~MVN (0, 𝛴 𝜏 ) .
Levels 2 and 3 are combined in the following equation with two random intercepts: one
at the visit level and the other at the participant level.
𝜃 𝑖𝑗
= 𝛽 0
𝛩 + 𝛽 1
𝛩 𝑋 𝑖𝑗
+ 𝛼 𝛩 𝑖 + 𝛼 𝛩 𝑖𝑗
(7)
Prior distributions for L-U-HB are specified to be relatively non-informative. We assume
the regression coefficients each have independent multivariate normal priors (with I
indicate a square identity matrix):
𝛽 0
𝛩 ~𝑀𝑉𝛮 ( 𝜇 𝛽 0
, 𝐼𝜎
𝛽 0
𝛩 2
) (8)
𝛽 1
𝛩 ~𝑀𝑉𝛮 ( 𝜇 𝛽 1
, 𝐼𝜎
𝛽 1
𝛩 2
) (9)
with 𝜇 𝛽 = ( 𝜇 𝛽 𝐶𝐴
, 𝜇 𝛽 𝑙𝑜𝑔𝐶𝑎𝑤 , 𝜇 𝛽 𝑙𝑜𝑔𝐷𝑎𝑤 ) , a non-informative prior distribution using large
variances 𝜎 𝛽 0
2
= 𝜎 𝛽 1
2
= 10
3
. The random intercept variance-covariance matrix 𝛴 𝜐 is
assumed to have a non-informative inverse-Wishart prior distribution with non-
informative diagonal matrix D=diag (0.001, 0.001, 0.001). And the variances of within-
participants are sampled from non-informative independent inverse gamma distributions.
Finally, the residual variance 𝜎 2
is assumed to have a non-informative inverse-Gamma
distribution, Inv-Gamma (0.001, 0.001).
97
The L-U-HB and L-TS-HB’s first stage was simulated via JAGS as mentioned above.
The simulation process includes an adaptive mode phase (“burn-in”) and a long enough
updating phase [66]where the adaptive mode is turned off.
6. L-U-NLME
The longitudinal version of the unified NLME model is similar to the cross-sectional one
in which the covariate is linked to the mean function for NO parameters, except that it
also specifies the variance-covariance matrix for the visit level subgroup. We also specify
the diagonal matrix for the visit level variations for our simplified model.
Constraint on CA
CA must be non-negative since it represents the concentration, in ppb, of the NO in the
alveolar compartment [54]. There are similar constraint considerations for Caw and Daw
which we satisfy by modeling logCaw and logDaw since Caw and Caw to have
approximately log-normal population-level distributions. CA tends to have more of a
normal or truncated normal distribution, so a different approach is necessary. In the
simulation study data generation step (described below), we discard samples if their CA
was negative. We also enforce the non-negative constraint in the HB models by using a
truncated distribution function. We didn’t apply the constraint on HMA since our
previous papers [16, 54] proved that it had poor performance, probably due to the large
number of failed to converge in Stage I. The NLME models implemented in the nlme
package in R were fitted without such constraints since there are no readily available
constraint options. Constrained versions of TS-NLS proved to be more biased in our
previous study[54], thus we implement only the unconstrained NLS here.
98
Simulation study
We compare the above methods in an extensive simulation study, roughly based on the
CHS study design. Each simulated dataset consists of 500 participants with 3 visits each
and each visit includes 8 multiple flow FeNO maneuvers (2 each at: 30, 50, 100, and 300
ml/s), which we simulate under a given “scenario” of underlying true associations of the
NO parameters with a standard normal covariate Xij (independent across and within
participants). For a given scenario, 100 replicate datasets (each of N* M =500*3) are
generated. Data-generating values, shown in Table 3, of population-level mean NO
parameters (𝛽 0
𝛩 ) , the between-participant variance-covariance matrix (𝛴 𝛩 𝜐 ), and the
residual variance 𝜎 𝜖 2
are based on values estimated in a preliminary L-TS-NLME
analysis of CHS data and are similar to the values in the previous cross-sectional
study[54]. We set the with-participant covariances to be zero for simplicity.
Table 3: Parameter values used to generate data in the simulation study
NO
parameter
Population-level (across-
participant)
Within-
participant
(across-visit)
Population-level (across-
participant)
Mean (𝜇 ) Standard
Deviation (𝜏 )
Standard
Deviation (𝜎 )
Correlation (𝜌 )
CA 1.5 0.45 0.50 CA, logCaw 0.66
logCaw 3.5 0.65 0.40 CA, logDaw -0.38
logDaw 2.5 0.55 0.33 logCaw, logDaw -0.35
99
The different scenarios for the simulation study are described in Table 4. The
regression coefficients 𝛽 1
𝛩 relating Xij to each NO parameter take a range of values: 0.01,
0.05 to 0.1. In Scenario 1, the reference scenario, all NO parameters have the same
association with the covariate Xij (𝛽 𝐶 𝐴 = 𝛽 𝑙𝑜𝑔𝐶 𝑎 𝑤 = 𝛽 𝑙𝑜𝑔𝐷 𝑎𝑤
), and the value of this
association is either 0.01, 0.05 or 0.1.
Table 4: The 4 simulation study scenarios are each repeated at 3 effect sizes and
replicated 100 times, for 1200 simulated datasets in total.
𝛽 𝐶 𝐴 𝛽 𝑙𝑜𝑔𝐶 𝑎𝑤
𝛽 𝑙𝑜𝑔𝐷 𝑎𝑤
Scenario 1 0.01,0.05,0.1 0.01,0.05,0.1 0.01,0.05,0.1
Scenario 2 0.01,0.05,0.1 0
†
0
†
Scenario 3 0
†
0.01,0.05,0.1 0
†
Scenario 4 0
†
0
†
0.01,0.05,0.1
†Cells marked 0 indicate that X had no effect on the corresponding NO parameter.
In this simulation study we compare performance of the methods based on several
metrics: percent bias (% bias), 95% Confidence/Credible (CI) length and coverage,
power and type I. Percent bias is calculated as: (Estimate-True value)/True value for non-
zero parameters. The 95% CI length and coverage are calculated using the given 95%
confidence interval for frequentist approaches and 95% credible interval (using the
returned posterior distribution, the center portion contains 95% of the values) for
Bayesian approaches. In primary analyses, Scenario 1 is used to calculate bias, power,
95% CI coverage, 95% CI length and the other scenarios (S2-S4) are used to calculate
type I errors. In secondary analyses, Scenarios 2-4 are used to calculate bias, power &
type I error rate, 95% CI coverage, 95% CI length.
100
CHS data analysis
We analyzed data from a CHS cohort originally recruited in kindergarten/1
st
grade [67]
with FeNO50 assessed at 6 study visits over 8 years (spanning ages 8-16) and multiple
flow FeNO assessed at the last 2 study visits, when most children were ages 13-14 and
15-16. The CHS multiple flow FeNO protocol called for 9 maneuvers at each of four
expiratory flow rates (3 at 50 ml/s and 2 each at: 30, 100, and 300 ml/s) collected
using chemiluminescence analyzers (model CLD88-SP with DeNOx accessory to provide
NO-scrubbed air; EcoMedics, Duernten, Switzerland/Ann Arbor, MI, USA) as described
in detail elsewhere [68, 69]. FeNO data processing was based on the ATS/ERS guidelines
for FeNO at 50 mL·s−1 [51] with a search window based on airway turnover[70]. Each
CHS child participant provided informed assent and a parent/guardian provided informed
consent. The CHS data were collected using a protocol approved by the University of
Southern California Institutional Review Board, the analyses in this paper were
conducted under HS-13–00150, and all methods were carried out in accordance with
relevant guidelines and regulations.
In a previous longitudinal analysis, FeNO50 was found to have a strong positive linear
association with height across this age range in children without asthma [71]. To
complement the previous analysis relating longitudinally assessed FeNO50 to height, here
we relate longitudinally assessed NO parameters (from the up to 2 repeated assessments
of multiple flow FeNO) to standardized height (population-mean centered: 162.7
centimeters and population-SD scaled: 8.75 centimeters). The analyses included 1004
children who never reported a doctor diagnosis of asthma and had multiple flow FeNO
101
data available at both visits. The average number of valid multiple flow maneuvers at the
first and second visits were 9.68 and 8.97, respectively.
102
Results
Simulation study
Computation time was longer for unified methods than for two-stage methods, as
expected. The computation time for a given method was similar across scenarios and
slightly shorter for larger β coefficients (Supplementary Table 2). Average computation
times on a high-performance computing platform (3 CPU, 12GB memory) for a single
simulated dataset (500 participants, 3 visits each, 8 maneuvers per visit) were: 30 hours
for L-U-HB, 23 hours for L-TS-HB, 14.7 minutes for L-U-NLME, 11.4 minutes for L-
TS-NLME, 4.6 seconds for L-TS-HMA and 6.2 seconds for L-TS-NLS. Most methods
had reasonable convergence rates (99% for L-TS-HB, 93% for L-U-HB, 94% for L-TS-
NLME) except for L-U-NLME (51%). L-TS-NLS and L-TS-HMA converged for all
datasets, however Stage I of L-TS-NLS had 29% of participant models fail to converge
on average (resulting in the exclusion of these participants’ results in Stage II) while L-
TS-HMA had only 0.016% failures in Stage I (Supplementary Figure 2). The following
simulation study results are a summary of all available converged results, since the
intersection of datasets which converged under all methods is relatively small due to the
convergence issue for L-U-NLME.
Figure 8: Comparison of method performance in terms of: percent bias (a), 95% CI
length (b), and 95% CI coverage (c) for estimating associations with NO parameters
103
(𝛽 𝐶𝑎
: black square, 𝛽 𝑙𝑜𝑔𝐶𝑎𝑤 : red circle, 𝛽 𝑙𝑜𝑔𝐷𝑎𝑤 : blue triangle) in simulation study
Scenario 1, replicated at 3 different effect sizes (𝛽 = 0.01, 0.05 or 0.1).
Figure 8 compares the percent bias and 95% CI interval properties for estimation of
𝛽 𝐶𝑎
, 𝛽 𝑙𝑜𝑔𝐶𝑎𝑤 , 𝛽 𝑙𝑜𝑔𝐷𝑎𝑤 across methods. For many methods, percent bias tended to
decrease as the effect size of the covariate increased. Among all methods, L-U-HB had
the lowest absolute values of percent bias for all three NO parameter associations (all <
4%), however at a small magnitude effect size (true β of 0.01) there was 52% bias for
𝛽 𝑙𝑜𝑔𝐷𝑎𝑤 , equivalent to a 0.052 bias on the original scale. L-U-NLME also had good
performance in the subset of datasets where the method converged. Two stage methods
104
(L-TS-HB, L-TS-NLME, L-TS-NLS) tended to have negative percent bias for
𝛽 𝐶𝑎
, 𝛽 𝑙𝑜𝑔𝐷𝑎𝑤 and positive percent bias for 𝛽 𝑙𝑜𝑔𝐶𝑎𝑤 . For a given method and effect size,
percent bias was smaller for 𝛽 𝐶𝑎
than for associations with other NO parameters, perhaps
because CA is in the linear part of the 2CM for FeNO. In simulation Scenario 3, when
only logCaw had an association with X, the directions of bias for two stage methods were
different from other scenarios (Supplemental Figure 1.C1). But L-U-HB and L-U-NLME
still had the smallest bias. An alternative version of L-TS-HMA with the constraint that
Stage I CA>0 in Scenario 1 (Supplementary Figure 3) resulted in ~30% of Stage I
estimates being dropped, on average, and patterns in percent bias similar to other two
stage methods.
For a given method, the 95% CI lengths (Figure 8b) were similar across NO parameters
and these patterns did not vary based on magnitude of the true β. L-TS-NLME had the
shortest 95% CI lengths followed by L-TS-HB. For 95% CI coverage (Figure 8c),
coverage declined as the magnitude of the true β increased, for all methods except L-U-
HB. It was the only method to produce 95% CIs with approximately appropriate coverage
(91 ~ 95%) for all 𝛽 s. L-U-NLME had slightly larger biases and shorter 95% CI lengths,
resulting in 75 ~ 80% coverage when the true β was 0.1. The higher coverage for these
two unified methods was due to a combination of low bias and longer 95% CI. L-TS-HB
also had reasonable coverage, but it differed across NO parameters, ranging from 63 ~
81%. For L-TS-NLME, L-TS-HMA, and L-TS-NLS, low coverage was due to a
combination of large bias and short CI lengths.
105
Figure 9: Power and Type I errors for the 6 methods (distinguished by color) for true β of
0.1, with power for a given NO parameter’s association calculated from the simulation
scenario where only that association is non-zero (e.g., power for 𝛽 𝐶𝑎
from Scenario 2
where 𝛽 𝐶𝑎
= 0.1, 𝛽 𝑙𝑜𝑔𝐶𝑎𝑤 = 0, 𝛽 𝑙𝑜𝑔𝐷𝑎𝑤 = 0) and Type I error calculated under two
scenarios (e.g., Type I error for 𝛽 𝐶𝑎
under S3: 𝛽 𝐶𝑎
= 0, 𝛽 𝑙𝑜 𝑔𝐶𝑎𝑤 = 0.1, 𝛽 𝑙𝑜𝑔𝐷𝑎𝑤 = 0 and
under S4: 𝛽 𝐶𝑎
= 0, 𝛽 𝑙𝑜𝑔𝐶𝑎𝑤 = 0, 𝛽 𝑙𝑜𝑔𝐷𝑎𝑤 = 0.1), with the non-zero NO parameter
association denoted by shape.
Power curves and Type I error rates for 𝛽 𝐶𝑎
, 𝛽 𝑙𝑜𝑔𝐶𝑎𝑤 , 𝛽 𝑙𝑜𝑔𝐷𝑎𝑤 across all simulation
scenarios are shown in Supplementary Figure 1. For simplicity, here Figure 9 displays a
subset of these data: the power for each NO parameter association versus two versions of
type I error, based on scenarios 2-4 at the largest magnitude effect size considered (true β
of 0.1). Ideally, a method will produce 2 values in the upper left-hand corner of this plot,
106
which indicates high power and low type I error regardless of which another NO
parameter had a non-zero association. Indeed, L-U-HB had relatively high power (1.00
for 𝛽 𝐶𝑎
, 0.97 for 𝛽 𝑙𝑜𝑔𝐶𝑎𝑤 , and 0.89 for 𝛽 𝑙𝑜𝑔𝐷𝑎𝑤 ) and low type I error rates from 0.02 to
0.07 for all three NO parameter associations except 0.12 for 𝛽 𝑙𝑜𝑔𝐶𝑎𝑤 in Scenario 4 where
only 𝛽 𝑙𝑜𝑔𝐷𝑎𝑤 was non-zero. Other methods generally also had good power, but L-TS-
NLS had low power for 𝛽 𝑙𝑜𝑔𝐷𝑎𝑤 and L-TS-HMA had low power for both 𝛽 𝑙𝑜𝑔𝐶𝑎𝑤 and
𝛽 𝑙𝑜𝑔𝐷𝑎𝑤 . Except for L-U-HB, most methods had inflated type I error for 𝛽 𝑙𝑜𝑔𝐶𝑎𝑤 when
𝛽 𝑙𝑜𝑔𝐷𝑎𝑤 was non-zero, or vice versa. For example, L-TS-HB had excellent power for
𝛽 𝐶𝑎
(0.98) and low type I error (0.02) when 𝛽 𝑙 𝑜𝑔𝐶𝑎𝑤 was nonzero but higher type I error
(0.17) when 𝛽 𝑙𝑜𝑔𝐷𝑎𝑤 was nonzero. This issue became more pronounced for 𝛽 𝑙𝑜𝑔𝐶𝑎𝑤 ,
where L-TS-HB again had excellent power for 𝛽 𝑙𝑜𝑔𝐶𝑎𝑤 (1.00) and low type I error (0.01)
when 𝛽 𝐶𝑎
had a nonzero association but very high type I error (0.69) when 𝛽 𝑙𝑜𝑔𝐷𝑎𝑤 was
nonzero. In summary, L-U-HB had high power and low type I error rates for all NO
parameter associations while other methods had good power and type I error for 𝛽 𝐶𝑎
and
good power but inflated type I error rates for 𝛽 𝑙𝑜𝑔𝐶𝑎𝑤 and 𝛽 𝑙𝑜𝑔𝐷𝑎𝑤 . Exceptions to this
pattern were L-TS-NLS and L-TS-HMA, both of which had low power for 𝛽 𝑙𝑜𝑔𝐷𝑎𝑤 due
to their large negative biases.
While the primary focus of the simulation study was on the estimation of 𝛽 𝐶𝑎
, 𝛽 𝑙𝑜𝑔𝐶𝑎𝑤 ,
and 𝛽 𝑙𝑜𝑔𝐷𝑎𝑤 , we also studied estimation of random effect variances at the participant and
visit levels. Participant-level variances, particularly for CA and logDaw, tended to be
underestimated by most methods, though L-U-HB had the overall lowest bias
107
(Supplemental Figure 1) Participant-level correlations also tended to be underestimated
(Supplemental Figure 1).
In summary, this simulation study demonstrated that L-U-HB generally had the smallest
bias, appropriate 95% CI coverage, and good power with low type I error rates for all
three NO parameter associations across scenarios. L-TS-HB, the two-step version of L-
HB, had greatly reduced computation time and maintained good power at the expense of
introducing some bias, poorer 95% CI coverage at larger magnitude effects, and high
type I errors for 𝛽 𝑙𝑜𝑔𝐶𝑎𝑤 and 𝛽 𝑙𝑜𝑔𝐷𝑎𝑤 . Compared to L-TS-HB, L-U-NLME had similar
performance, and even less inflated type I error rates but failed to converge in ~50% of
the simulated datasets. Compared to L-TS-HB, L-TS-NLME had bias in the same
direction but of larger magnitude, resulting in lower coverage and more inflated type I
errors. L-TS-NLS, on average, had ~40% of participants fail to have Stage I estimates
(Supplementary Table 1). Despite this, L-TS-NLS performed well for estimation of 𝛽 𝐶 𝐴 ,
but for 𝛽 𝑙𝑜𝑔𝐶𝑎𝑤 and 𝛽 𝑙𝑜𝑔𝐷𝑎𝑤 had considerable bias, low power and inflated type I error.
CHS data analysis
Applying the 6 methods to a CHS analysis relating NO parameters to height, we observed
associations with height that were: positive or null for CA, positive for logCaw, and
negative or null for logDaw (Figure 4). The two unified methods (L-U-HB and L-U-
NLME) both estimated similar statistically significant associations between height and all
three NO parameters. Specifically, from the L-U-HB model we estimated that, a within-
108
participant increase in height of 8.79 centimeters was associated, on average, with a
0.079 (95% CI: 0.034, 0.125) ppb increase in CA, a 0.158 (95% CI: 0.106, 0.212) increase
in logCaw, and a 0.106 (95% CI: 0.044, 0.171) decrease in logDaw. These latter two
estimates are equivalent to a 17% increase in Caw and a 10% decrease in Daw. From the L-
U-NLME model, analogous estimates were similar: 0.092 (95% CI: 0.046, 0.138) for CA,
0.149 (95% CI: 0.104, 0.194) for logCaw, -0.104 (95% CI: -0.157, -0.052) for logDaw).
Two-stage methods produced lower estimates for 𝛽 𝐶𝑎
and 𝛽 𝑙𝑜𝑔𝐶𝑎𝑤 , and higher estimates
for 𝛽 𝑙𝑜𝑔𝐷𝑎𝑤 . Furthermore, two-stage method estimates were all approximately null and
non-significant for 𝛽 𝑙𝑜𝑔𝐷𝑎𝑤 . Estimates of 𝛽 𝐶𝑎
were not statistically significant for L-TS-
NLS and L-TS-HMA. There was a clear pattern when comparing a unified method to its
two-stage counterpart (i.e., L-U-HB vs L-TS-HB and L-U-NLME vs L-TS-NLME) in
that the unified method had larger magnitude estimates (farther from zero) than their
corresponding two-stage versions. The finding that the unified methods produced more
significant estimates than two-stage ones (especially the L-TS-NLS and L-TS-HMA),
which agreed with results in our simulation. For L-TS-NLS, a Stage I outlier (-1325) for
CA resulted in an extremely wide 95% CI for 𝛽 𝐶𝑎
in Stage II, so we excluded this outlier
for Figure 4. Stage I convergence failures were observed for 528 and 30 out of 1004
models for L-TS-NLS and L-TS-HMA, respectively.
Figure 10: Estimated associations between NO parameters and standardized height in the
CHS using 6 methods.
109
110
Discussion
In this paper, we proposed a novel unified hierarchical Bayesian model, L-U-HB, for
relating longitudinally assessed NO parameters to covariates, extending our previous
cross-sectional unified hierarchical Bayesian model, and performed the first evaluation of
the statistical properties of various two-stage methods for longitudinal analysis of NO
parameters. In a simulation study, L-U-HB performed well for estimating associations of
NO parameters with covariates, with small bias, appropriate 95% CI coverage, good
power, and low type I error rates. The two-step version, L-TS-HB, had greatly reduced
computation time and maintained good power at the expense of introducing some bias,
poorer 95% CI coverage at larger magnitude effects, and high type I errors for
𝛽 𝑙𝑜𝑔𝐶𝑎𝑤 and 𝛽 𝑙𝑜𝑔𝐷𝑎𝑤 . The other unified method L-U-NLME had good performance
when it converged, but it had serious convergence issues. Other two-stage methods had
drawbacks in terms of bias, inflated type I error rates, etc.
In a previous simulation study comparing the performance of methods estimating NO
parameter associations with a covariate in a cross-sectional study [54], U-HB had the best
performance across all simulation scenarios, similar to our findings in this longitudinal
study. L-U-NLME had the second-best performance across all three NO parameters in
the longitudinal study, while its cross-sectional version (U-NLME) had large bias in
estimating 𝛽 𝐶𝑎
. L-TS-NLS also had much better performance in estimating 𝛽 𝐶𝑎
and
𝛽 𝑙𝑜𝑔𝐶𝑎𝑤 in the longitudinal study than TS-NLS for the cross-sectional study but still had
large bias for 𝛽 𝑙𝑜𝑔𝐷𝑎𝑤 . Both L-TS-NLME and the cross-sectional TS-NLME had large
bias and poor coverage.
111
Computation time is an issue. L-U-HB had superior statistical properties to many
competitor methods, albeit at additional computational expense. The cross-sectional U-
HB model had an average computation time of 5.5 hours for N=1000 participants (ref)
while for the longitudinal version (L-U-HB) average computation time was 30 hours for
N=500 participants, 3 visits each. While the longer L-U-HB computation time was
burdensome in a simulation study with 1000s of datasets, it is less of an issue when
analyzing a single dataset. However, note that the computational cost will increase as
more covariates are added, or a more complex model is used in Stage II. The other
competitive method was L-U-NLME, which had the second-best estimation performance
and ran faster, but it had poor convergence. Our results suggest that for further
applications in which there are more variables or more complex models, where the
unified model may become computational intractable, a two-stage version of hierarchical
Bayesian model will have reasonably good performance and therefore be used in iterative
model building, with a single run of L-U-HB for the final results, using L-TS-HB
estimates as starting values to speed convergence.
For biologically plausibility, we constrained CA to be non-negative. In our previous
cross-sectional study, we encountered a problem constraining CA to be non-negative in
the standard JAGs software. We solved it by sampling logCaw and logDaw as they were
bivariate normal distributed, and then sample zero-truncated CA conditioned on them. In
that case, we sampled the variances and correlations one by one and set boundaries for
the last correlation to ensure a 3× 3 positive definite variance covariance matrix. The same
solution was also used for L-U-HB and L-TS-HB when we set NO parameters to be
correlated in both levels. But in our simplified model which assumed no correlations in
112
the visit level, the randomness and constraint on CA were easily specified separately
without considering the validity of the variance-covariance matrix.
When applying L-U-HB to study the association between height and NO parameters
using longitudinal multiple flow FeNO data on healthy schoolchildren in the CHS, we
found positive associations of CA and logCaw with height and a negative association of
logDaw with height. Had we applied only L-TS-HMA or L-TS-NLS methods, we would
have failed to detect associations of height with CA or logDaw. Our findings add to the
limited literature on associations of NO parameters with height/age. A previous analysis
using longitudinal FeNO50 data from the same cohort over a longer follow-up period,
from ages 8-16, found that FeNO50 increased approximately linearly with height and
FeNO50 increased nonlinearly with age [71]. Limited cross-sectional data on trends in NO
parameters by age (for participants less than 20 years old) suggests non-significant
increases of Daw and Caw but a decrease in CA though some influential values for the
oldest participants in the sample may have impacted these results [72]. Additionally,
differences with our findings may be due to a different age range or due to the difference
between a cross-sectional (between-person) vs longitudinal (within-person) design.
There are several directions for future work. To reduce computation time observed when
implementing L-U-HB with JAGS, alternative GIBBS sampling software such as Rstan
could be considered. Several components of the model implementation appeared to affect
convergence rates, such as the number of parameters or the length of adaptation phase.
The length of adaptation was the most important factor, but adaptation is in itself a
computationally intensive operation. Given the computational costs involved, we chose to
base our simulation study on smaller dataset sizes (500 participants, 3 visits each) and
113
simplify our model so that it could converge within two days on average. The simplified
version of our longitudinal model ignored the correlation between the NO parameters
within the participant level (V0 model). L-U-HB and L-TS/U-NLME were able to
specify the diagonal matrix for visit-level variation. While the NLS and the HMA
approaches fitted the FeNO models for each observation, thus ignoring any correlations
between or within participants. In the CHS data analysis, we only selected participants
who have visited both in year 8 and 10 as a more direct comparison to the simulation
study, but unified models are capable of handling unbalanced data.
In conclusion, in this paper we presented a longitudinal extension of the unified
hierarchical Bayesian model for analyzing nonlinear data (e.g., FeNO data). Despite the
long computation time required for achieving convergence, L-U-HB had the best
performance estimating covariate coefficients as well as variance-covariance
components. Its two-stage version L-TS-HB served as a decent alternative to obtain a
faster pilot estimation for much more complicated models.
114
References
1. Jain, N. and J.L. Hill, The National Heart Lung and Blood Institute guidelines finally
say "yes" to fractional exhaled nitric oxide. Ann Allergy Asthma Immunol, 2022.
2. Khatri, S.B., et al., Use of Fractional Exhaled Nitric Oxide to Guide the Treatment of
Asthma: An Official American Thoracic Society Clinical Practice Guideline. Am J
Respir Crit Care Med, 2021. 204(10): p. e97-e109.
3. Dweik, R.A., et al., An official ATS clinical practice guideline: interpretation of
exhaled nitric oxide levels (FENO) for clinical applications. Am J Respir Crit Care Med,
2011. 184(5): p. 602-15.
4. La Grutta, S., et al., Environmental Effects on Fractional Exhaled Nitric Oxide in
Allergic Children. Journal of Allergy, 2011. 2012.
5. Scarpa, M.C., N. Kulkarni, and P. Maestrelli, The role of non-invasive biomarkers in
detecting acute respiratory effects of traffic-related air pollution. Clin Exp Allergy, 2014.
44(9): p. 1100-18.
6. Oliver, A.C., et al., Effects of Reduced Nicotine Content Cigarettes on Fractional
Exhaled Nitric Oxide and Self-Reported Respiratory Health Outcomes Among Smokers
With Psychiatric Conditions or Socioeconomic Disadvantage. Nicotine Tob Res, 2022.
24(1): p. 135-140.
7. Society, A.T., European Respiratory Society. ATS/ERS recommendations for
standardized procedures for the online and offline measurement of exhaled lower
115
respiratory nitric oxide and nasal nitric oxide, 2005. Am J Respir Crit Care Med, 2005.
171: p. 912-930.
8. SC, G., et al., Modeling pulmonary nitric oxide exchange. Journal of applied
physiology (Bethesda, Md. : 1985), 2004. 96(3).
9. Roy, K., et al., Use of different exhaled nitric oxide multiple flow rate models in
COPD. Eur Respir J, 2007. 29(4): p. 651-9.
10. Eckel, S.P., et al., Estimation of parameters in the two-compartment model for
exhaled nitric oxide. PLoS One, 2014. 9(1): p. e85471.
11. Karvonen, T., et al., Comparison of feasibility and estimates of central and peripheral
nitric oxide parameters by different mathematical models. J Breath Res, 2017. 11(4): p.
047102.
12. Weng, J., et al., Hierarchical Bayesian estimation of covariate effects on airway and
alveolar nitric oxide. Sci Rep, 2021. 11(1): p. 17180.
13. Karrasch, S., et al., Accuracy of FENO for diagnosing asthma: a systematic review.
Thorax, 2017. 72(2): p. 109-116.
14. HL, P., et al., Exhaled nitric oxide levels to guide treatment for adults with asthma.
The Cochrane database of systematic reviews, 2016. 9(9).
15. Silkoff, P.E., et al., Longitudinal stability of asthma characteristics and biomarkers
from the Airways Disease Endotyping for Personalized Therapeutics (ADEPT) study.
Respir Res, 2016. 17: p. 43.
116
16. Kharitonov, S.A., et al., Reproducibility of exhaled nitric oxide measurements in
healthy and asthmatic adults and children. Eur Respir J, 2003. 21(3): p. 433-8.
17. Garcia, E., et al., Patterns and determinants of exhaled nitric oxide trajectories in
schoolchildren over a 7-year period. Eur Respir J, 2020. 56(1).
18. Linn, W.S., et al., Exhaled nitric oxide in a population-based study of Southern
California Schoolchildren, in Respir Res. 2009. p. 28.
19. Hö gman, M. and P. Merilä inen, Extended NO analysis in asthma. J Breath Res,
2007. 1(2): p. 024001.
20. Hö gman, M., et al., Extended NO analysis applied to patients with COPD, allergic
asthma and allergic rhinitis. Respir Med, 2002. 96(1): p. 24-30.
21. Silkoff, P.E., et al., Airway nitric oxide diffusion in asthma: Role in pulmonary
function and bronchial responsiveness. Am J Respir Crit Care Med, 2000. 161(4 Pt 1): p.
1218-28.
22. Latif, A.H. and S.G. Gilmour, Transform-both-sides nonlinear models for in vitro
pharmacokinetic experiments. Stat Methods Med Res, 2015. 24(3): p. 306-24.
23. Pinheiro, J., et al., Nlme: Linear and Nonlinear Mixed Effects Models. R package
version 31-110, 2013. 3: p. 1-113.
24. Plummer, M., JAGS: A program for analysis of Bayesian graphical models using
Gibbs sampling. undefined, 2003.
25. Plummer, M. JAGS: A program for analysis of Bayesian graphical models using
Gibbs sampling. 2003.
117
26. McConnell, R., et al., Traffic, susceptibility, and childhood asthma. Environ Health
Perspect, 2006. 114(5): p. 766-72.
27. Linn, W.S., et al., Multiple-flow exhaled nitric oxide, allergy, and asthma in a
population of older children. Pediatr Pulmonol, 2013. 48(9): p. 885-96.
28. Linn, W.S., et al., Extended exhaled nitric oxide analysis in field surveys of
schoolchildren: a pilot test. Pediatr Pulmonol, 2009. 44(10): p. 1033-42.
29. Puckett, J.L., et al., Impact of analysis interval on the multiple exhalation flow
technique to partition exhaled nitric oxide. Pediatr Pulmonol, 2010. 45(2): p. 182-91.
30. Garcia, E., et al., Patterns and determinants of exhaled nitric oxide trajectories in
schoolchildren over a 7-year period. 2020.
31. Hö gman, M., et al., Effects of growth and aging on the reference values of
pulmonary nitric oxide dynamics in healthy subjects - IOPscience. 2017.
118
Supplementary materials
We now show a number of further details of results of our simulation study. We begin by
exploring convergence rates across a range of methods and scenarios.
Table 1a: Number of simulation study datasets (out of 100) with successful convergence
for each method and scenario.
Method Scenario 1 Scenario 2 Scenario 3 Scenario 4
𝛽 𝐶𝑎
= 𝛽 𝑙𝑜𝑔𝐶𝑎𝑤 =
𝛽 𝑙𝑜𝑔𝐷𝑎𝑤
𝛽 𝑙𝑜𝑔𝐶𝑎𝑤 = 𝛽 𝑙𝑜𝑔𝐷𝑎𝑤 = 0 𝛽 𝐶𝑎
= 𝛽 𝑙𝑜𝑔𝐷𝑎𝑤 = 0 𝛽 𝐶𝑎
= 𝛽 𝑙𝑜𝑔𝐶𝑎𝑤 = 0
𝛽 𝐶𝑎
𝛽 𝑙𝑜𝑔𝐶𝑎𝑤 𝛽 𝑙𝑜𝑔𝐷𝑎𝑤
0.01 0.05 0.1 0.01 0.05 0.10 0.01 0.05 0.10 0.01 0.05 0.10
L-U-HB 95 100 89 89 97 98 91 88 95 90 97 92
L-TS-HB 100 99 99 97 99 100 98 100 100 100 100 99
L-U-NLME 54 55 49 49 48 48 46 45 47 56 55 56
L-TS-NLME 96 96 93 96 95 93 95 96 92 93 92 96
L-TS-HMA 100 100 100 100 99 100 100 100 100 100 100 99
L-TS-NLS 100 100 100 100 100 100 100 100 100 100 100 100
All methods
converged
53 54 41 45 45 47 41 39 44 51 54 50
119
Table 1b: Average Stage I convergence failure rates* for select two-stage methods: L-TS-
HMA and L-TS-NLS.
1
Method Scenario 1 Scenario 2 Scenario 3 Scenario 4
𝛽 𝐶𝑎
= 𝛽 𝑙𝑜𝑔𝐶𝑎𝑤 =
𝛽 𝑙𝑜 𝑔 𝐷𝑎𝑤
𝛽 𝑙𝑜𝑔𝐶𝑎𝑤 = 𝛽 𝑙𝑜𝑔𝐷𝑎𝑤 = 0 𝛽 𝐶𝑎
= 𝛽 𝑙𝑜𝑔𝐷𝑎𝑤 = 0 𝛽 𝐶𝑎
= 𝛽 𝑙𝑜𝑔𝐶𝑎𝑤 = 0
𝛽 𝐶𝑎
𝛽 𝑙𝑜𝑔𝐶𝑎𝑤 𝛽 𝑙𝑜𝑔𝐷𝑎𝑤
0.01 0.05 0.1 0.01 0.05 0.10 0.01 0.05 0.10 0.01 0.05 0.10
L-TS-HMA 0.2 0.2 0.3 0.2 0.2 0.3 0.3 0.3 0.2 0.3 0.3 0.2
L-TS-NLS 435.3 435.6 439.3 432.4 434.5 435.1 433.5 433.5 433.6 438.2 434.2 436.3
Table 2: Average Computation time for each method and simulation scenario.
Average computation times on a high-performance computing platform (3 CPU, 8GB
memory) for a single simulated dataset (500 participants, 3 visits each, 8 maneuvers per
visit) were: 29 hours for L-U-HB, 24 hours for L-TS-HBL, 14.7 minutes for L-U-NLME,
11.6 minutes for L-TS-NLME, 4.7 seconds for L-TS-HMA and 6.1 seconds for L-TS-
1
Stage I convergence failure rate is the number of participant’s Stage I models that failed to converge divided by the total number of participants (500)
120
NLS. However, the computation times differed slightly across scenarios. Table 2
summarizes the computation time of each scenario for all methods.
Scenario
Scenario 1 Scenario 2 Scenario 3 Scenario 4 overa
ll
𝛽 𝐶𝑎
= 𝛽 𝑙𝑜𝑔𝐶𝑎𝑤 =
𝛽 𝑙𝑜𝑔𝐷𝑎𝑤
𝛽 𝑙𝑜𝑔𝐶𝑎𝑤 = 𝛽 𝑙𝑜𝑔𝐷𝑎𝑤
= 0
𝛽 𝐶𝑎
= 𝛽 𝑙𝑜𝑔𝐷𝑎𝑤 = 0 𝛽 𝐶𝑎
= 𝛽 𝑙𝑜𝑔𝐶𝑎𝑤 = 0
𝛽 𝐶𝑎
𝛽 𝑙𝑜𝑔𝐶𝑎𝑤 𝛽 𝑙𝑜𝑔𝐷𝑎𝑤
0.01 0.05 0.1 0.01 0.05 0.10 0.01 0.05 0.10 0.01 0.05 0.10
L-U-HB
(hours ）
28 32.3 28.4 33.4 31.1 31.5 31.6 30.2 30.2 30.4 28 30
L-TS-HB
(hours)
25.4 22.7 23.7 22.6 23 21.8 23 21.7 23.3 21.3 22.7 22.4
L-U-NLME
(minutes)
18.5 13.7 11.7 15.2 17.3 14.7 16.5 14.6 14.4 12.3 18.2 15.3
L-TS-NLME
(minutes)
12.4 11.4 11 12.2 12.1 12.6 12.7 11.1 11.9 12.9 13.4 12.5
L-TS-HMA
(seconds)
4.9 4.6 4.7 4.7 4.4 4.5 4.4 4.5 4.5 4.5 4.5 4.6
L-TS-NLS
(seconds)
6.1 6.2 5.8 6.2 5.8 5.9 5.7 5.9 6.1 5.8 6 5.8
121
Figure 1: Extensive comparison of method performance in Simulation Scenarios with
subplots for the different parameters
The performance of methods were compared in terms of bias, 95% CI length and 95% CI
coverage for coefficients, populational means, participant level standard deviations,
participant level correlations and visit level correlations of NO parameters (𝛽 𝐶𝑎
: black
square, 𝛽 𝑙𝑜𝑔𝐶𝑎𝑤 : red circle, 𝛽 𝑙𝑜𝑔𝐷𝑎𝑤 : blue triangle), replicated at 3 different effect sizes
(𝛽 = 0.01, 0.05 or 0.1). Power and type I error rates for coefficients. These results are
shown below.
A: Scenario 1
A1: Estimated Associations β
122
A2: Population level mean ɑ
123
A3: Participant level standard deviations τ
124
A4: Participant level correlations ⍴
125
A5: Visit level standard deviations σ
126
B: Scenario 2
B1: Estimated Associations β (Bias instead of % Bias)
127
128
B2: Population level mean ɑ
B3: Participant level standard deviations τ
129
B4: Participant level correlations ⍴
130
B5: Visit level standard deviations σ
131
C: Scenario 3
C1: Estimated Associations β (Bias instead of % Bias)
132
133
C2: Population level mean ɑ
C3: Participant level standard deviations τ
134
C4: Participant level correlations ⍴
135
C5: Visit level standard deviations σ
136
D: Scenario 4
D1: Estimated Associations β (Bias instead of % Bias)
137
138
D2: Population level mean ɑ
D3: Participant level standard deviations τ
139
D4: Participant level correlations ⍴
140
D5: Visit level standard deviations σ
141
142
Figure 2: Summary on all converged datasets (44% - 64%) for Scenario 1 Estimated
Associations β
2
Figure 2.a: The percentage bias (% bias), 95 % confidence interval length and the
coverage of 95 % confidence intervals
2
Instead of using all available converged results, we summarize those datasets converged
by all methods.
143
Figure 2.b: The power vs type I error rates:
144
Figure 3: Constrain the HMA Stage I to be non-negative for NO parameters and ENO.
We previously used all estimations from Stage I for HMA analysis. However, if we use
only non-negative estimations, the bias patterns were similar to other two-stage methods
but at a cost of 26.4% drop rate on average.
Table 3: CHS summary
145
The effects of height to the NO parameters among CHS study data were summarized in
Table 3 in the form of Mean (95% CI). Significant associations were highlighted in
yellow.
Table 3.a: Mean (95% CI) for coefficient of standardized height
Models 𝛽 𝐶𝑎
𝛽 𝑙𝑜𝑔𝐶𝑎𝑤 𝛽 𝑙𝑜𝑔𝐷𝑎𝑤
L-U-HB 0.079 (0.034, 0.125) 0.158 (0.106, 0.212) -0.106 (-0.171, -0.044)
L-TS-HB 0.046 (0.014, 0.078) 0.05 (0.027, 0.072) 0.006 (-0.015, 0.028)
L-U-NLME 0.092 (0.046, 0.138) 0.149 (0.104, 0.194) -0.104 (-0.157, -0.052)
L-TS-NLME 0.038 (0.008, 0.067) 0.057 (0.033, 0.081) 0 (-0.021, 0.021)
L-TS-HMA -0.019 (-0.098, 0.06) 0.063 (0.009, 0.118) -0.006 (-0.058, 0.046)
L-TS-NLS -0.145 (-1.83, 1.54) 0.08 (0.033, 0.127) -0.016 (-0.057, 0.025)
* On stage I, HMA: 401 failed; NLS: 542 failed
Full results
146
parameter
L-U-HB L-TS-HB L-U-NLME L-TS-NLME L-TS-HMA L-TS-NLS
ɑ
𝐶𝑎
1.943
(1.892,
1.995)
1.939
(1.902,
1.976)
1.732
(1.682,
1.782)
1.726
(1.692,
1.761)
1.914
(1.831,
1.997)
1.006 (-
0.704,
2.716)
ɑ
𝑙𝑜𝑔𝐶𝑎𝑤
4.412
(4.341,
4.485)
4.392
(4.359,
4.426)
4.011
(3.958,
4.064)
4.004
(3.97,
4.039)
4.152
(4.095,
4.209)
4.129
(4.078,
4.18)
ɑ
𝑙𝑜𝑔𝐷𝑎𝑤
2.167
(2.086,
2.246)
2.187
(2.146,
2.227)
2.629
(2.569,
2.69)
2.639
(2.604,
2.674)
2.574
(2.521,
2.628)
2.635
(2.592,
2.677)
𝜏 𝐶𝑎
0.507
(0.464,
0.557)
0.528
(0.443,
0.599)
0.629
(0.588,
0.674)
0.458
(0.427,
0.491)
NA NA
𝜏 𝑙𝑜𝑔𝐶𝑎𝑤
0.709
(0.653,
0.77)
0.696
(0.644,
0.754)
0.52
(0.461,
0.586)
0.522
(0.497,
0.55)
NA NA
𝜏 𝑙𝑜𝑔𝐷𝑎𝑤
0.861
(0.788,
0.942)
0.843
(0.776,
0.927)
0.702
(0.645,
0.764)
0.546
(0.521,
0.573)
NA NA
147
𝜌 𝐶𝑎 ,𝑙𝑜𝑔𝐶𝑎𝑤
0.816
(0.76,
0.883)
0.849
(0.773,
0.905)
0.834
(0.705,
0.91)
NA NA NA
𝜌 𝐶𝑎 ,𝑙𝑜𝑔𝐷𝑎𝑤
-0.172 (-
0.308, -
0.044)
-0.208 (-
0.338, -
0.058)
-0.448 (-
0.561, -
0.318)
NA NA NA
𝜌 𝑙𝑜𝑔𝐶𝑎𝑤 ,𝑙𝑜𝑔𝐷𝑎𝑤
-0.706 (-
0.759, -
0.648)
-0.69 (-
0.746, -
0.631)
-0.571 (-
0.633, -
0.502)
NA NA NA
𝜎 𝐶𝑎
0.659
(0.619,
0.698)
0.656
(0.615,
0.699)
0.369
(0.345,
0.394)
0.468
(0.448,
0.488)
NA NA
𝜎 𝑙𝑜𝑔𝐶𝑎𝑤
0.359
(0.331,
0.388)
0.366
(0.339,
0.392)
0.615
(0.573,
0.66)
0.304
(0.291,
0.317)
NA NA
𝜎 𝑙𝑜𝑔𝐷𝑎𝑤
0.37
(0.333,
0.408)
0.367
(0.331,
0.402)
0.355
(0.323,
0.39)
0.243
(0.232,
0.253)
NA NA
148
Chapter 4: Statistical modeling of multiple flow exhaled nitric
oxide in population research
Abstract:
Multiple flow FeNO is of increasing interest and appears in an increasing number of
research studies. However, understanding of appropriate statistical methods to quantify
associations of patient characteristics (e.g., asthma severity), treatment, and other factors
with estimated NO parameters remains limited. Researchers typically use an ad hoc two-
stage method to: (1) estimate NO parameters and (2) treat the NO parameters as known
quantities and relate them to factors of interest. We summarize various statistical methods
for estimating NO parameters and relating them to factors of interest and highlight the
strengths/weaknesses of these methods. Compared to two-stage methods, unified
methods will likely have statistical advantages such as reduced bias, increased power, and
reduced issues with inflated Type I error due to correlations between NO parameters.
However, computation is more intensive. We have developed a freely available R
package to implement both two-stage and unified methods. Using data from a population-
based study of school children, we illustrate use of the software to estimate differences in
NO parameters comparing healthy children to those with asthma only, allergic rhinitis
only, or both asthma and allergy using a variety of statistical methods and compare the
results.
149
Introduction
FeNO (fractional concentration of exhaled nitric oxide) is a marker of airway
inflammation, measured non-invasively in exhaled breath, with clinical[73, 74] and
epidemiological applications.[4, 75] FeNO has most commonly been assessed at the
standardized expiratory flow rate of 50 ml/s (FeNO50),[9, 76] treating flow dependency as
a nuisance. However, new insights are being gained by using FeNO measured across a
range of flow rates to partition FeNO into its proximal (airway wall) and distal (alveolar)
contributions.[14] This “multiple flow FeNO” technique is noninvasive, relying instead
on a combination of mathematical and statistical models to estimate biologically
meaningful airway wall and alveolar parameters. The methodology of multiple flow
FeNO assessment poses technical challenges, at least some of which can be addressed by
innovations in statistical methodology.
Deterministic mathematical models have been developed to describe the dynamics of NO
in the lower respiratory tract, using systems of differential equations with parameters
describing NO sources in airway and alveolar compartment.[37, 77, 78] The steady-state
two compartment model (2CM) is simple and robust, and hence widely used.[13] A
variety of statistical approaches have been developed to estimate NO parameters in the
2CM from multiple flow FeNO data, ranging from simple linear regression to more
complex nonlinear regression.[13, 16, 79, 80]
While there has been considerable
150
attention to the underlying deterministic physiological mathematical model and to the
statistical models estimating NO parameters from these models using FeNO data
measured with error, there has been little discussion on the issue of relating estimated NO
parameters to factors of interest (e.g., disease status, asthma severity, treatment),
henceforth referred to as “covariates”.
This paper focuses on the statistical issues in multiple flow FeNO research in population
studies, wherein the goal is to relate NO parameters to covariates. Typically, researchers
have used an ad hoc two stage method. In Stage I, NO parameters are estimated, often in
a separate regression model for each participant (for each visit, in the case of longitudinal
data). In Stage II, the estimated NO parameters are treated as known and used as the
outcomes in regression models relating the NO parameters to covariates. There are some
potential issues and missed opportunities with this approach. In particular, the two-stage
approach ignores the (potentially unequal) uncertainty of the estimated NO parameters
and ignores correlation across NO parameters. In previous work we have developed and
evaluated an alternative unified estimation approach, implemented using Hierarchical
Bayesian methodology, for both cross-sectional[81] and longitudinal study designs.
Unified estimation methods overcome some of the issues of two stages approaches, albeit
at the expense of additional computational complexity and time.
151
The purpose of this paper is threefold: (1) to review the statistical issues in multiple flow
FeNO research in population studies, (2) to provide an overview of relevant statistical
analysis approaches while highlighting their strengths and weaknesses, and (3) to
demonstrate implementation of these analyses using freely available statistical software,
facilitating the translation of these statistical methods developments to the breath research
community. Specifically, we introduce a new R package
(https://github.com/USCbiostats/PopFeNO) for estimating NO parameters in population
research studies which enables simple and convenient implementation and comparison of
two-stage and unified estimation approaches. Using multiple flow FeNO data from the
Southern California Children’s Health Study, we illustrate use of this software to estimate
differences in NO parameters by asthma/allergy status using a variety of two stage and
unified approaches.
152
Methods
Multiple flow FeNO data
Multiple flow FeNO studies require data from
each participant on FeNO across a range of
expiratory flow rates.[39] In a single FeNO
maneuver, the concentration of NO at the mouth
(FeNO) is measured in real-time along with the
expiratory flow rate while the participant exhales
steadily and slowly, aiming to achieve a target
flow rate.[76] The maneuver is summarized
using a “plateau” mean FeNO and the
corresponding mean expiratory flow rate. Figure
1 shows an example of the paired FeNO and flow
data from 9 maneuvers for a single participant.
Two compartment model (2CM)
Figure 11: Multiple flow FeNO data used to
estimate NO parameters* from the two-
compartment model for a single participant.
* For this example CHS participant, estimated NO
parameters are: 𝐶 ̂
𝐴 𝑁𝑂 = 6.73 ( 95% 𝐶𝐼 : 4.01, 9.39) ,
𝐶 ̂
𝑎𝑤
𝑁𝑂 = 276.60 ( 95% 𝐶𝐼 : 199.23, 570.62) , and
𝐷 ̂
𝑎𝑤
𝑁𝑂 = 10.51 ( 95% 𝐶𝐼 : 4.48, 16.59) from the NLS
estimation method.
153
At high flow rates, FeNO primarily reflects lower
concentrations in the alveolar region, while at low
flow rates it primarily reflects higher concentrations in
the airway wall tissue. This flow dependency can be
leveraged to gain insights about the proximal versus
distal sources of NO. The popular steady-state two
compartment model (2CM) assumes a cylindrical
airway compartment and an expansile alveolar
compartment (Figure 2).[13] In the 2CM model,
FeNO is related to the following parameters
quantifying airway and alveolar NO sources: 𝐶 𝐴 𝑁𝑂 , the concentration of NO in the
alveolar region, (ppb); 𝐶 𝑎𝑤
𝑁𝑂 , the concentration of NO in the airway wall tissue, (ppb);
and 𝐷 𝑎𝑤
𝑁𝑂 , the airway wall tissue diffusion capacity, pL· s
-1
· ppb
-1
) and expiratory flow
rate (ml/s) in the following equation:
𝐹𝑒𝑁𝑂 = 𝐶 𝑎𝑤
+ ( 𝐶 𝐴 − 𝐶 𝑎𝑤
)× 𝑒 −
𝐷 𝑎𝑤
𝑓𝑙𝑜𝑤 (1)
where we have dropped the “NO” from the “NO parameter” notation for simplicity.
Estimating NO parameters for a single participant
Figure 12: Steady state two-
compartment model (2CM) and key
NO parameters.
154
NO parameters from the 2CM can be estimated for a study participant at a given visit
using the multiple flow FeNO data described earlier (e.g., 3-10 paired flow/FeNO
values[39]). The statistical methods for estimating NO parameters range in complexity
from simple linear regression to nonlinear regression.[13, 16, 79, 80] Simple linear
regression approaches rely on a first order linear approximation to the exponential
function in Equation 1, which is only valid for high enough flow rates such that the ratio
of 𝐷 𝑎𝑤
𝑁𝑂 to flow is small.[13, 44, 82, 83] An implication of these simplifications and
assumptions is that not all NO parameters can be estimated using linear approximation
methods (specifically, 𝐶 𝑎𝑤
𝑁𝑂 and 𝐷 𝑎𝑤
𝑁𝑂 cannot be estimated), hence we do not
consider these methods further here. A widely used estimation method is the Hö gman and
Merilӓinen algorithm (HMA),[18, 38] which consists of a custom iterative algorithm
using a third order approximation to the exponential function in Equation 1, which takes
as input 3 paired values: average FeNO from “low”, “medium”, and “high” target flow
rates, and includes an internal data consistency check. Rather than approximating the
exponential function in Equation 1, there are statistical methods which can readily
estimate parameters from so called “nonlinear” models. While it has been proposed to use
a nonlinear regression model directly estimating 2CM parameters using the most direct
translation of Equation 1 (adding a normally distributed error term to the righthand
side),[17] it has been demonstrated that a modification introducing a natural log
155
transformation of both sides of the 2CM equation has improved statistical properties
because it accounts for the increasing variation of FeNO with flow:[16]
𝑙𝑜𝑔 ( 𝐹𝑒𝑁𝑂 𝑗 ) = 𝑙𝑜𝑔 (𝐶 𝑎𝑤
+ ( 𝐶 𝐴 − 𝐶 𝑎𝑤
)× 𝑒 −
𝐷 𝑎𝑤
𝑓𝑙𝑜𝑤 𝑗 )+ 𝜀 𝑗 (2)
where the errors 𝜀 𝑗 are normally distributed. We henceforth refer to this log transform
both sides nonlinear regression method as “NLS” for nonlinear least squares.
Estimating NO parameters and their association with covariates in a population study
There are two general approaches to estimating associations of NO parameters with
covariates: two stage methods and unified methods, henceforth denoted by “TS” or “U”,
respectively.
Two stage estimation methods. The conceptual framework of two stage methods,
displayed in Figure 3, has been described previously.[81] Briefly, in Stage I, NO
parameters are estimated for each participant (and for each visit in longitudinal study). In
Stage II, the estimates from Stage I of the three NO parameters for participant i (𝐶 ̂
𝑖 𝐴 ,
𝑙𝑜𝑔 𝐷 ̂
𝑖 𝑎𝑤
, and 𝑙𝑜𝑔 𝐷 ̂
𝑖 𝑎𝑤
) are treated as known values and used as the outcomes in three
156
separate linear regression models relating each NO parameter to a covariate, referred to
generically as X.
Stage I models typically consist of separate NLS or HMA models for each participant.
One downside of this approach is that due to small sample size (e.g., 9 maneuvers or 3
“average” FeNO values, respectively), Stage I models can fail to converge for a subset of
Figure 13: Conceptual diagram of two-stage estimation analysis in multiple flow exhaled nitric
oxide population research, where associations of a covariate X with NO parameters are of
interest (indicated by grey box). In Stage I, NO parameters are estimated and then in Stage II
estimated NO parameters are related to covariate X, all in a population of N study participants,
each participant indexed by i.
157
participants. Participants without Stage I results are then dropped from Stage II analyses,
causing potential selection bias in the estimated associations of NO parameters with
covariates. In previous work, this problem was particularly pronounced for NLS in Stage
I.[81] Even if Stage I models do converge, due to small sample sizes the estimated NO
parameters may be extreme (influenced by outliers) and/or estimated with varying
degrees of uncertainty.
In population studies, the Stage I small sample size issues can be partially remedied by
using a mixed effect model in Stage I, where a single Stage I model is fit to the pooled set
of multiple flow datasets for all participants. A mixed model extension to NLS is a
nonlinear mixed effects (NLME) model[25] of the form:
log ( 𝐹𝑒𝑁𝑂 𝑖𝑗
) = 𝑓 ( 𝜃 𝑖 , flow
𝑖𝑗
)+ 𝜀 𝑖𝑗
where:
𝑓 ( 𝜃 𝑖 , flow
𝑖𝑗
)= 𝑙𝑜𝑔 (exp( 𝑙𝑜𝑔 𝐶 𝑎𝑤
𝑖 )+ ( 𝐶 𝐴 𝑖 − exp( 𝑙𝑜𝑔 𝐶 𝑎𝑤
𝑖 ) )× 𝑒 −
exp( 𝑙𝑜𝑔 𝐷 𝑎𝑤
𝑖 )
𝑓𝑙𝑜𝑤 𝑖𝑗
) (3)
and i indexes participant and j indexes maneuver. Note that again modeling natural log
transformations, this time of certain NO parameters (𝐶 𝑎𝑤
𝑁𝑂 and 𝐷 𝑎𝑤
𝑁𝑂 ), improves the
statistical properties of this modeling approach by both better reflecting the
approximately log normal (right skewed) distribution of these parameters typically
observed in study populations while simultaneously imposing non-negativity constraints
on these parameters.[81] Special considerations for the distribution of 𝐶 𝐴 𝑁𝑂 will be
158
discussed later for more sophisticated models, but note that attempts to constrain
𝐶 𝐴 𝑁𝑂 > 0 within NLS consistently underestimate 𝐶 𝑎𝑤
𝑁𝑂 compared to unconstrained
versions.[81] While all participants will have estimated NO parameters if the NLME
model converges, standard methods for estimating NLME models suffer from well-
documented convergence problems. One alternative is to use a Hierarchical Bayesian
(HB) analog (described in greater detail later) to the frequentist NLME in Stage I.
However, both the NLME and HB approaches are most naturally suited to a unified
estimation approach rather than a two-stage approach.
Unified estimation methods. Unified methods simultaneously estimate both NO
parameters and their associations with covariates under a hierarchical model structure,
with maneuvers nested within participants, as shown in the conceptual diagram in Figure
4. For more compact notation, we define the vector of NO parameters for each participant
to be 𝜃 𝑖 = ( 𝐶 𝐴 𝑖 , log 𝐶 𝑎𝑤
𝑖 , log 𝐷 𝑎𝑤
𝑖 ) ′. Conceptually, unified methods are appealing
because they inherently account for unequal uncertainty in estimating NO parameters
across participants (e.g., from varying numbers of valid FeNO maneuvers) and account
for correlation of NO parameters, rather than ignoring it, when estimating associations
with the covariate X. This last feature reduces the potential for Type I error. For example,
159
if X had a true effect on only 𝐶 𝑎𝑤
𝑁𝑂 , but 𝐶 𝑎𝑤
𝑁𝑂 and 𝐶 𝐴 𝑁𝑂 were correlated, ignoring
the correlation could result in finding a spurious association between X and 𝐶 𝐴 𝑁𝑂 .
Previously, we have proposed two unified methods, U-NLME and U-HB.[81] The
frequentist method, U-NLME, is the unified version of TS-NLME, where the NO
parameters are modeled as a function of the covariate X, as implemented in the function
nlme() in the R package nlme[84] (Appendix 2). For the Bayesian method, U-HB, the
hierarchical model structure is similar to that of U-NLME except that we also incorporate
prior distributions for the parameters (e.g., the prior for 𝛽 1
( 𝐶 𝐴 )
is a normal distribution
with hyperparameters, i.e., mean and variance, selected to be noninformative, as
described in greater detail elsewhere[81]). Results of U-HB are approximations of the full
Figure 14: Conceptual diagram of unified estimation analysis in multiple flow exhaled
nitric oxide population research, where associations of a covariate X with NO parameters
are of interest (indicated by grey box) as estimated from a population of N study
participants, with each participant indexed by i and each maneuver by j.
160
posterior distributions of the parameters of interest rather than point MLE estimates and
standard errors. We have implemented the U-HB method using a Gibbs sampling Markov
chain Monte Carlo (MCMC) algorithm via the R2jags package interface to JAGS in R
(Appendix 2). While the hierarchical model in Figure 4 describes a basic model for
multiple flow FeNO population studies that are cross-sectional in design, it can also be
modified to introduce additional complexities (e.g., additional levels, adjustment to
confounders).
Extensions to longitudinal study designs
Similar to FeNO50, some of the greatest insights from multiple flow FeNO data may from
come from the study of within-participants changes in NO parameters over time in
longitudinal study designs. The methods for multiple flow FeNO population studies can
be readily extended to longitudinal data, as we have previously described [cite]. Briefly,
two stage methods now consist of Stage I: estimate NO parameters separately for each
participant at each visit and Stage II: relate these estimated NO parameters to covariates
in standard models for longitudinal data (e.g., linear mixed models). For unified methods,
an additional level is added to the hierarchical model so that now there are three levels:
maneuvers nested in visits nested in participants. Since the estimation methods had
similar performance in cross-sectional and longitudinal studies (with U-HB having clear
advantages but requiring longer computation time), we focus the remainder of this paper
161
on cross-sectional studies for simplicity, but the example code for longitudinal studies is
also provided in https://github.com/USCbiostats/PopFeNO.
Strengths and weaknesses of the six estimation methods
In extensive simulation studies, we have previously compared the performance of six
methods: TS-NLS, TS-HMA, TS-NLME, U-NLME, TS-HB, U-HB. In Table 1 we
provide a summary of the strengths and weaknesses of the various approaches for
estimating associations of NO parameters with a covariate X in a cross-sectional study
design. In brief, unified methods, particularly U-HB, tended to outperform the other
methods, at the expense of longer computation times.[81]
Table 5: Highlights of prior simulation studies, summarizing the strengths and
weaknesses of six methods for multiple flow FeNO population studies.
Method Degree
of Bias
95% CI
coverage
Power Type I
error
Convergence/
missing data
issues
Compute
time*
TS-NLS Severe Very low High Inflated Frequent < 10
seconds
TS-
HMA
Largely Low Low
(CawNO)
Inflated Frequent < 10
seconds
TS-
NLME
Largely Very low High Highly
Inflated
(CawNO,
DawNO)
Rare < 10
minutes
U-
NLME
Moderate Sometimes
low
High Inflated
(DawNO)
Frequent < 20
minutes
162
TS-HB
†
Little High High Inflated
(CawNO,
DawNO)
Rare < 5 hours
U-HB Very
little
High High Rarely
Inflated
Rare < 6 hours
* Computation times are general examples for a cross-sectional multiple flow FeNO population study with 1000
participants, 8 FeNO maneuvers per participant, and one covariate X.
†
TS-HB was not included in the simulation study for cross-sectional FeNO data in the (cite), we used the results from
longitudinal simulation study in the paper
R package
We have developed the R package https://github.com/USCbiostats/PopFeNO which
enables easy implementation of the six methods described here for multiple flow FeNO
population studies. Two-stage method functions take as input the multiple flow FeNO
data and the covariate(s) and return a dataset of estimated NO parameters with covariates
as the Stage I outcome, enabling the user to readily customize Stage II regression models
(e.g., readily try adjustment for different confounders or evaluate nonlinearity of the
association with X). Unified method functions produce all estimates in one function, so
the U-NLME and U-HB functions require additional syntax to specify the form of the
associations between NO parameters and covariates (Appendix 2).
163
Application to CHS data. The Southern California Children’s Health Study (CHS) is a
population-based prospective cohort study which includes one of the largest and longest
studies of FeNO in schoolchildren to date,[27, 28, 85-87] including multiple flow FeNO
assessments at 2 visits for many participants.[69] Here, we demonstrate two-stage and
unified estimation techniques for multiple flow FeNO studies using cross-sectional data
on multiple flow FeNO from 1618 CHS participants collected at the March 2010—June
2010 study visit. We aim to evaluate differences in NO parameters by allergy/asthma
status, similar to a prior analysis.[68] Responses to questionnaires were used to assess
whether a doctor had ever diagnosed the child with asthma (“asthma”) or whether in the
past 12 months, “your child had a problem with sneezing or a runny or blocked nose
when he/she did not have a cold or the flu?” (“allergy”). Participants were categorized
into four groups based on these questionnaire responses: “healthy” reference (N=943: no
asthma and never/not current allergy), allergy only (N=163), asthma only (N=355) and
having both asthma and allergy (N=157). Multiple flow FeNO was collected using
chemiluminescence analyzers (Model 1080i, EcoMedics, Ann Arbor, MI/Duernten,
Switzerland), as described previously,[68, 69] and processed using guideline-based
methods tailored to the CHS, described previously.[22] The protocol[68, 69] called for 9
maneuvers per participant (3 at 50 ml/s, followed by 2 each at: 30, 100, and 300 ml/s),
164
but in practice not all participants had precisely 9 valid FeNO maneuvers. Our primary
analysis related NO parameters to the four allergy/asthma groups without additional
covariate adjustment, as shown below for 𝐶 𝐴 𝑁𝑂 :
𝐶 𝐴 = 𝛽 0
+ 𝛽 1
ⅹ𝐼 𝑎𝑙𝑙𝑒𝑟𝑔𝑦 + 𝛽 2
ⅹ𝐼 𝑎𝑠𝑡 ℎ𝑚𝑎
+ 𝛽 3
ⅹ𝐼 𝑎𝑙𝑙𝑒𝑟𝑔𝑦 & 𝑎𝑠𝑡 ℎ𝑚𝑎
+ 𝜀
where the covariates are indicators for inclusion in the allergy only, asthma only, or
allergy and asthma groups. Hence, the parameters of interest are 𝛽 1
, 𝛽 2
, and 𝛽 3
which
represent the average difference in the NO parameter (here, 𝐶 𝐴 𝑁𝑂 ) between a given
group (e.g., for 𝛽 1
the allergy only group) and the reference “healthy” group. In a
secondary analysis, we additionally adjusted for age (mean centered at 13.4 years), height
(mean centered at 160 cm), and gender. Estimation of NO parameters and their
associations with asthma/allergy groups was conducted via six methods: four two-stage
methods (TS-NLS, TS-HMA, TS-NLME, TS-HB) and two unified methods (U-NLME,
U-HB). Note that to implement the HMA method in the CHS, we use average FeNO
values from 30, 100 and 300 ml/s flow rates and did not drop NAs or negative estimated
CAs to keep consistent with the approaches used by other frequentist methods in earlier
analyses.[16] Uncertainties are reported as 95% confidence intervals for frequentist
methods and 95% credible intervals for Bayesian methods, both denoted 95% CI. For
Bayesian methods, the reported “estimate” is the mean of posterior distribution.
165
Results
The 1618 participants ranged in age from 12 to 15 years old, 48.7% were male, and
93.88% had at least 9 valid maneuvers (range: 4 to 14). Figure 5 shows the estimated
differences in NO parameters in each of the asthma/allergy groups, compared to the
reference “healthy” group, using the six estimation methods. Results from U-HB model
showed an approximately additive effect of allergy and asthma on 𝐶 𝐴 and logCaw,
whereas for logDaw values were similarly elevated in children with allergy, asthma, or
both asthma and allergy. On average, 𝐶 𝐴 𝑁𝑂 was 0.30 ppb higher in children with allergy
(95 % CI: 0.10, 0.53), 0.18 higher in children with asthma (95 % CI: 0.05, 0.34), and 0.58
higher in children with allergy and asthma (95 % CI: 0.37, 0.80). On average, 𝐶 𝑎𝑤
𝑁𝑂
was 16% higher in children with allergy (95 % CI: -5%, 42%), 17% higher in children
with asthma (95 % CI: 1%, 37%), and 67% higher in children with allergy and asthma
(95 % CI: 37%, 105%). 𝐷 𝑎𝑤
𝑁𝑂 was 29% higher in children with allergy (95 % CI: 1%,
61%), 17% higher in children with asthma (95 % CI: 1%, 37%), and 21% higher in
children with allergy and asthma (95 % CI: -6%, 52%). Estimated differences were
slightly lower after adjusting for age, sex, and height.
Generally, estimates of differences in 𝐶 𝐴 𝑁𝑂 by asthma/allergy status were the most
stable across the 6 estimation methods (e.g., 5 of 6 methods estimated increasingly large
differences for allergy, asthma, allergy and asthma) though only 4 methods (TS-HMA,
166
TS-NLME, TS-HB, U-HB) found all 3 allergy/asthma groups to be statistically
significantly different from the reference group. Results for 𝐶 𝑎𝑤
𝑁𝑂 and 𝐷 𝑎𝑤
𝑁𝑂 were
less consistent across estimation methods. For 𝐶 𝑎𝑤
𝑁𝑂 , all methods estimated the largest
difference to be between the asthma and allergy group and the reference group (all
statistically significant), though this estimated difference was quite large for TS-NLME.
Additionally, for 𝐶 𝑎𝑤
𝑁𝑂 4 methods (TS-HMA, U-NLME, TS-HB, U-HB) estimated a
similar magnitude difference between the allergy vs reference group as compared to the
difference for the asthma vs reference group. For 𝐷 𝑎𝑤
𝑁𝑂 , 3 methods (TS-NLS, TS-
HMA, U-HB) estimated similar differences with largely overlapping CI for all 3
asthma/allergy groups while 3 methods (TS-NLME, U-NLME, TS-HB) estimated a
lower difference for the allergy group.
When comparing results across the estimation methods, we found similar patterns of the
differences across methods as in previous simulation studies. [cite]. Results for the two
unified methods, U-NLME and U-HB were most similar for 𝐶 𝑎𝑤
𝑁𝑂 , while the modest
variability across unified methods in results for 𝐶 𝐴 𝑁𝑂 and 𝐷 𝑎𝑤
𝑁𝑂 might be partly
explained by the greater uncertainty (wider 95% CI) for these parameters. Results for U-
HB and its two-stage analog TS-HB were most similar for 𝐶 𝐴 𝑁𝑂 and 𝐶 𝑎𝑤
𝑁𝑂 , while the
95% CI for differences in 𝐷 𝑎𝑤
𝑁𝑂 were considerably wider for U-HB resulting in fewer
statistically significant differences for 𝐷 𝑎𝑤
𝑁𝑂 . Comparing U-NLME and its two-stage
167
analog TS-NLME, the estimates from TS-NLME tended to be larger. In previous
simulation studies, we had found that two-stage analogs of unified methods often
produced overly narrow 95% CI. For the simplest (and most commonly used) two-stage
methods, TS-NLS and TS-HMA, results for 𝐶 𝐴 𝑁𝑂 were similar though estimated
differences from TS-HMA were slightly larger. TS-HMA had particularly long 95% CI
for differences in 𝐷 𝑎𝑤
𝑁𝑂 resulting in none of the estimated differences being
statistically significant, which is as reported in previous simulation studies [cite].
Figure 15: Estimated differences in NO parameters comparing asthma/allergy groups to a
reference (healthy) group of children, at the 2010 study visit in the Southern California
Children’s Health Study. Six estimation methods were used: TS-NLS, TS-HMA, TS-
NLME, U-NLME, TS-HB, U-HB.
168
Computation times for the six estimation methods were: 2.13 seconds for TS-NLS, 1.53
seconds for TS-HMA, 2.6 minutes for TS-NLME, 6.9 minutes for U-NLME, 1.90 hours
for TS-HB, and 3 hours for U-HB. Due to convergence failures/inadequate data, Stage I
estimates were unavailable for 363 participants using TS-NLS and 24 using TS-HMA
resulting in the exclusion of these participants from Stage II of these methods.
When considering the suitability of two stage versus unified methods for these data, it is
interesting to note that there was variation in the uncertainty of estimated NO parameters
across study participants, as indicated by the following Stage I results from TS-NLS: the
169
SE of estimated 𝐶 𝐴 𝑁𝑂 ranging from 0.022 to 33.77 (median: 0.54, mean: 0.94), the SE
of estimated log 𝐶 𝑎𝑤
𝑁𝑂 ranging from 0.046 to 10.42 (median: 0.46, mean: 0.79), and the
SE of estimated log 𝐷 𝑎𝑤
𝑁𝑂 ranging from 0.063 to 10.85 (median: 0.64, mean: 0.98).
There was also correlation in NO parameters. U-HB estimated the correlations were 0.59
between 𝐶 𝐴 𝑁𝑂 and log 𝐶 𝑎𝑤
𝑁𝑂 , - 0.16 between 𝐶 𝐴 𝑁𝑂 and log 𝐷 𝑎𝑤
𝑁𝑂 and - 0.55
between log 𝐶 𝑎𝑤
𝑁𝑂 and log 𝐷 𝑎𝑤
𝑁𝑂 .
170
Discussion
In this paper, we describe a variety of statistical methods to quantify associations of
covariates (e.g., patient characteristics, treatment, exposures) with estimated NO
parameters, highlighting their relative strengths and weaknesses, with the aim of helping
data analysts in the breath research community to select appropriate statistical methods
that provide the best power to detect associations of interest. We demonstrated the
application of six methods to estimate differences in NO parameters by asthma/allergy
status in a cross-sectional, population-based study of schoolchildren. While ad hoc two
stage methods have been almost exclusively applied in practice, unified methods,
particularly a unified Hierarchical Bayesian method (U-HB), have appealing statistical
properties and can be readily implemented using standard statistical software, as
demonstrated in the R code and accompanying tutorial we provide at:
https://github.com/USCbiostats/PopFeNO. Given that the U-HB method is more
computationally intensive than alternative methods, the two-stage analog (TS-HB)—
which has slightly worse statistical performance but much faster computation time—
could be used in interactive model building and in more complex modeling scenarios,
ideally followed by a final run with the U-HB model.
171
Literature on statistical methods for population multiple flow FeNO studies is limited[81]
and already summarized in this paper. However, there is a large related literature in
population pharmacokinetic modeling,[88, 89] where it is of interest to relate patient
characteristics to pharmacokinetic parameters. In population pharmacokinetics, two stage
methods have been generally supplanted by unified modeling approaches, as we suggest
here, due to the advantages of unified methods in handling sparse participant-level
data[90] (analogous to a small number of FeNO maneuvers per participant) and in
appropriately estimating between-participant variances.[91] Hierarchical Bayesian
models are well established population pharmacokinetics methodology, having been in
use for decades.[92] By adapting the methodological developments in this related field,
we can drive forward progress in multiple flow FeNO methodology.
For investigators conducting population multiple flow FeNO studies, we provide R code
and an accompanying tutorial at: https://github.com/USCbiostats/PopFeNO to aid in the
translation of these statistical methods to practical data analyses. We provide code for
both two-step and unified approaches, both frequentist and Bayesian. In general, the U-
HB model is recommended due to its advantageous statistical properties (including
superior power and reduced type I errors), however it can have a long runtime (i.e., time
to convergence) for large datasets and/or complex relationships with covariates. When
172
such runtimes present problems, one might first use the comparatively quick-fitting TS-
HB model to obtain rough estimates of NO parameters and variance-covariance matrices.
Intermediate model building that is more iterative or interactive (e.g., determining
whether covariates associations are linear or nonlinear, identifying potential confounders)
could use the participant-level NO parameter estimates from the HB Stage I of TS-HB as
outcomes in standard Stage II models. One could then use the resultant TS-HB estimates
to provide good ‘ball-park’ starting values for the full U-HB model specification to fit the
final model, thereby improving likely runtime for that final analysis.
In conclusion, multiple flow FeNO is an area of breath research with great potential[14]
but it is still under development. Growing interest is reflected in the inclusion of multiple
flow FeNO in the most recent update to the guidelines for FeNO assessment.[76]
Multiple flow FeNO has been included in a growing number of clinical studies, for
example both cross-sectional studies[93-97] and studies with repeated visits for each
participant,[98-106] as well as several epidemiological studies.[22, 107-109] All these
studies have used two stage analyses, which do not take full advantage of the rich data
structure and can suffer from biases and inflated type I errors. Researchers should
consider alternative methods, such as U-HB, in their analyses of population data on
multiple flow FeNO to maximize the potential of their research.
173
Appendix
Appendix 1: Example R code for population multiple flow FeNO studies
Here we provide example code for analyzing population multiple flow FeNO data with
both two-stage and unified approaches. Details of data manipulation were simplified, the
complete code and results is in the demo file for our package, it can be found in Github
(https://github.com/USCbiostats/PopFeNO).
Set up workspace and load required packages
set.seed(2022)
rm(list=ls())
setwd(“##”)
library(MASS)
library(lme4)
library(nlme)
library(reshape2)
library(plotrix)
library(latex2exp)
Prepare dataset:
load("FeNO.RData") # dat is the dataset name
Print the head of FeNO dataset:
174
Gather the variables according for each participant, separate by the visit entries (only
use Year 8)
AdjustX<-aggregate(
Enodata[,c("id","yr","heightC","maleC","ageC",
"aNO","bmipctl","race2n","Allergy/Asthma",
"pflowlev")],
by=list(Enodata$yr_id),FUN = unique)
AdjustX<-AdjustX[order(AdjustX$id,AdjustX$yr),]
AdjustX8<-AdjustX[AdjustX$yr==8,]
Run simulations:
Two-stage methods:
Fit the two-stage NLS model for example:
1. Estimate NO parameters
NLS_S1_Y8<-TS_NLS_S1(Enodata8,AdjustX8)
2. Fit three linear regressions for NO them
175
Ca_Yr8_NLS<-lme(Ca~AA2, random=~1|id,
data = NLS_S1_Y8,na.action = na.omit)
logCaw_Yr8_NLS<-lme(logCaw~AA2, random=~1|id,
data = NLS_S1_Y8,na.action = na.omit)
logDaw_Yr8_NLS<-lme(logDaw~AA2, random=~1|id,
data = NLS_S1_Y8,na.action = na.omit)
3. Summary the results in NLSY8 dataset
R code files of function for each model:
TS_NLS_s1.R: NLS model for Stage I
TS_HMA_s1.R: HMA model for Stage I
TS_NLME_s1.R: NLME model for Stage I
TS_HB_s1.R: NLME model for Stage I
Unified methods: Estimate NO parameters and coefficients simultaneously.
U-NLME:
1. Fit the TS model as initial model. (Optional, to get initial values and it was easier to
converge)
2. Update the TS model with covariates:
fitU.nlme_Y8 <- update(UNLME_S1_fit_Y8,
fixed = list(CaNO ~ AA2,
176
logCawNO ~ AA2,
logDawNO ~ AA2),
start = c(startFix_Y8[1],0,0,0,
startFix_Y8[2],0,0,0,
startFix_Y8[3],0,0,0),
control=list(tolerance=0.01)
)
U-HB:
For different models, it is required to customize initial values, input dataset, JAGS model
functions.
For example:
Ca_mean[i] <-beta0_Ca + beta1_Ca * X1[i]
+ beta2_Ca * X2[i]
+ beta3_Ca * X3[i]
logCaw_mean[i]<-beta0_logCaw + beta1_logCaw * X1[i]
+ beta2_logCaw * X2[i]
+ beta3_logCaw * X3[i]
logDaw_mean[i]<-beta0_logDaw + beta1_logDaw * X1[i]
+ beta2_logDaw * X2[i]
+ beta3_logDaw * X3[i]
inits_U <-
list(
list(
tau_c =runif(1,10,200),
beta0_Ca =runif(1,0,5),
beta0_logCaw =runif(1,0,5),
beta0_logDaw =runif(1,0,5),
177
beta1_Ca =runif(1,-0.5,0.5),
beta1_logCaw =runif(1,-0.5,0.5),
beta1_logDaw =runif(1,-0.5,0.5),
beta2_Ca =runif(1,-0.5,0.5),
beta2_logCaw =runif(1,-0.5,0.5),
beta2_logDaw =runif(1,-0.5,0.5),
beta3_Ca =runif(1,-0.5,0.5),
beta3_logCaw =runif(1,-0.5,0.5),
beta3_logDaw =runif(1,-0.5,0.5),
tau_sub =solve(SubM())
),
list(
same as the first list
),
list(
same as the first list
)
)
Summarize and compare the estimates from different models
I used the plotCI function here:
plotCI(x=xvals,err="y",y=c(0,NLSY8$Ca_est[2:4],
0,NLSY8$logCaw_est[2:4],
0,NLSY8$logDaw_est[2:4])+10.5,
li=c(0,NLSY8$Ca_lb[2:4],
0,NLSY8$logCaw_lb[2:4],
0,NLSY8$logDaw_lb[2:4])+10.5,
ui=c(0,NLSY8$Ca_ub[2:4],
0,NLSY8$logCaw_ub[2:4],
178
0,NLSY8$logDaw_ub[2:4])+10.5,
yaxt="n",xaxt="n",xlim=c(0,38.5),ylim=c(0,12),
xlab="",ylab="",pch=19,lty=1,cex=1.5)
...
Same function to add results for other methods
…
Appendix 2: Customized U-NLME and U-HB:
R code for customized U-NLME:
TS_NLME_S1<-function(dat,tol1,tol2){
mydat <- subset(dat,select=c(yr_id,logeno,flow,id))
dat2 <- groupedData(logeno~ flow | id, data=mydat)
fitnlsList <- nlsList(logeno~SSnonLinLogPlogJwoGrad(
flow,logCawNO,CaNO,logDawNO),
data=dat2, control=list(tolerance = 0.1))
cf <- coef(fitnlsList)
startfix <- unlist(lapply(cf, median, na.rm = TRUE))
idREstart <- matrix(0,nrow=length(unique(mydat$id)),ncol=3)
colnames(idREstart) <- names(startfix)
rownames(idREstart) <- unique(mydat$id)
idvisitmat <- matrix(as.integer(unlist(
strsplit(rownames(coef(fitnlsList)),"/"))),
byrow=T,ncol=1)
for(i in 1:3){
fillcol <- tapply(cf[,i],
as.factor(idvisitmat[,1]),mean,na.rm=T)
idREstart[which(!is.nan(fillcol)),i] <-
fillcol[which(!is.nan(fillcol))]-startfix[i]
}
179
fit1.nlme <- nlme(logeno~SSnonLinLogPlogJwoGrad(
flow,logCawNO,CaNO,logDawNO),
data=mydat,
fixed=logCawNO + CaNO +logDawNO ~ 1,
random= list(pdLogChol(
logCawNO + CaNO +logDawNO ~ 1)),
groups=~id,
start=list( fixed=c(1.82,4.12,2.63),
random=list(id=idREstart)
),
verbose=T,
control=list(tolerance =
tol1,opt="nlminb",natural=TRUE,niterEM=50))
fit2.nlme <- update(fit1.nlme,
control=list(tolerance =tol2,maxIter=50),
verbose=T)
return(fit2.nlme)
}
UNLME_S1_fit_Y8<-TS_NLME_S1(Enodata8,0.1,0.01)
fitU.nlme_Y8 <- update(UNLME_S1_fit_Y8,
fixed = list(
CaNO ~ AA2, logCawNO ~ AA2, logDawNO ~ AA2),
start = c(startFix_Y8[1],0,0,0,
startFix_Y8[2],0,0,0,
startFix_Y8[3],0,0,0),
control=list(tolerance=0.01)
)
The argument of fixed effects and starting values should be modified according to the
linear regression models for NO parameters and covariates.
U-HB:
180
The customized JAGS model function:
Offset: used to localize participant for unbalanced data
Truncated distribution: ;T(0,) added after the distribution function
Code in Blue: Conditional distribution of CA|(logCaw, logDaw). And the sampling
strategy to assure the positive definitive of the variance – covariance matrix
Code in Yellow: Customized for the linear relationship between NO parameters and
covariates.
Jags_U_CHS <- function() {
for( i in 1 : Nid) {
for( j in offset[i]:(offset[i+1]-1)) {
logeno[j] ~ dnorm(logmu[j],tau_c)
mu[j] <- exp(logCaw[i])+
(Ca[i]-exp(logCaw[i]))*
exp(-1*exp(logDaw[i])/flow[j])
logmu[j]<-log(max(mu[j] ,0.001))#logFeNO>0.001
}
Ca[i] <- param_subject[i, 1]
logCaw[i] <- param_subject[i, 2]
logDaw[i] <- param_subject[i, 3]
param_subject[i,1] ~ dnorm(mu1_22[i],tau1_22);T(0,)
mu1_22[i]<-Ca_mean[i]+part1[1,1]*(param_subject[i,2]
-logCaw_mean[i])
+part1[1,2]*(param_subject[i,3]
-logDaw_mean[i])
param_subject[i,2:3] ~ dmnorm(
c(logCaw_mean[i],logDaw_mean[i]),
tau22[1:2,1:2])
181
#customize your Xs
Ca_mean[i] <-beta0_Ca + beta1_Ca * X1[i]
+ beta2_Ca * X2[i]
+ beta3_Ca * X3[i]
logCaw_mean[i]<-beta0_logCaw + beta1_logCaw * X1[i]
+ beta2_logCaw * X2[i]
+ beta3_logCaw * X3[i]
logDaw_mean[i]<-beta0_logDaw + beta1_logDaw * X1[i]
+ beta2_logDaw * X2[i]
+ beta3_logDaw * X3[i]
##**************##
}
#---sub---#
tau1_22<-1/cov1_22
############# variance of A Conditioned on B
cov1_22[1,1]<-varCa-(covCalogCaw^2*tau22[1,1]
+2*covCalogDaw*tau22[2,1]*covCalogCaw
+covCalogDaw^2*tau22[2,2])
############# SigmaAB * tauB
part1[1,1]<-covCalogCaw*tau22[1,1]+covCalogDaw*tau22[2,1]
part1[1,2]<-covCalogCaw*tau22[1,2]+covCalogDaw*tau22[2,2]
############# Sigma B
tau22[1:2,1:2]<-inverse(cov22[1:2,1:2])
cov22[1,1] <- varlogCaw
cov22[1,2] <- covlogCawlogDaw
cov22[2,1] <- covlogCawlogDaw
cov22[2,2] <- varlogDaw
############# var-cov
covCalogCaw <-corlogCawCa*sqrt(varCa*varlogCaw)
covCalogDaw <-corlogDawCa*sqrt(varCa*varlogDaw)
covlogCawlogDaw <-corlogCawlogDaw*sqrt(varlogCaw*varlogDaw)
varCa<-1/tauCa
varlogCaw<-1/taulogCaw
varlogDaw<-1/taulogDaw
sdCa <- sqrt(varCa)
182
sdlogCaw <- sqrt(varlogCaw)
sdlogDaw <- sqrt(varlogDaw)
sigma_c<-sqrt(1/tau_c)
############# sample correlation
corlogDawCa ~ dunif(L_corlogDawCa,U_corlogDawCa)
L_corlogDawCa <- corlogCawCa*corlogCawlogDaw
-sqrt((corlogCawlogDaw^2-1)*(corlogCawCa^2-1))
U_corlogDawCa <- corlogCawCa*corlogCawlogDaw
+sqrt((corlogCawlogDaw^2-1)*(corlogCawCa^2-1))
corlogCawCa ~ dunif(-1,1)
corlogCawlogDaw ~ dunif(-1,1)
############# sample tau
tauCa ~ dgamma(1.0E-3, 1.0E-3)
taulogCaw ~ dgamma(1.0E-3, 1.0E-3)
taulogDaw ~ dgamma(1.0E-3, 1.0E-3)
#---error---#
tau_c ~ dgamma(1.0E-3, 1.0E-3)#sigma_C
############# sample mean
beta0_Ca ~ dnorm(beta0_prior[1],1e-6)
beta0_logCaw ~ dnorm(beta0_prior[2],1e-6)
beta0_logDaw ~ dnorm(beta0_prior[3],1e-6)
#customize your Xs
beta1_Ca ~ dnorm(betaC_prior[1], 1e-6)
beta1_logCaw ~ dnorm(betaC_prior[2], 1e-6)
beta1_logDaw ~ dnorm(betaC_prior[3], 1e-6)
beta2_Ca ~ dnorm(betaC_prior[1], 1e-6)
beta2_logCaw ~ dnorm(betaC_prior[2], 1e-6)
beta2_logDaw ~ dnorm(betaC_prior[3], 1e-6)
beta3_Ca ~ dnorm(betaC_prior[1], 1e-6)
beta3_logCaw ~ dnorm(betaC_prior[2], 1e-6)
beta3_logDaw ~ dnorm(betaC_prior[3], 1e-6)
##**************##
}
183
References
1. Braun, P.X., C.F. Gmachl, and R.A. Dweik, Bridging the Collaborative Gap:
Realizing the Clinical Potential of Breath Analysis for Disease Diagnosis and
Monitoring–Tutorial. Sensors Journal, IEEE, 2012. 12(11): p. 3258-3270.
2. Gustafsson, L., et al., Endogenous nitric oxide is present in the exhaled air of
rabbits, guinea pigs and humans. Biochem Biophys Res Commun, 1991. 181(2):
p. 852-857.
3. Dweik, R.A., et al., An official ATS clinical practice guideline: interpretation of
exhaled nitric oxide levels (FENO) for clinical applications. Am J Respir Crit
Care Med, 2012. 184(5).
4. La Grutta, S., et al., Environmental effects on fractional exhaled nitric oxide in
allergic children. J Allergy (Cairo), 2012. 2012: p. 916926.
5. Scarpa, M.C., N. Kulkarni, and P. Maestrelli, The role of non-invasive biomarkers
in detecting acute respiratory effects of traffic-related air pollution. Clin Exp
Allergy, 2014. 44(9): p. 1100-1118.
6. Hogman, M., et al., Nitric oxide from the human respiratory tract efficiently
quantified by standardized single breath measurements. Acta Physiol Scand,
1997. 159(4): p. 345-6.
7. Silkoff, P.E., et al., Marked flow-dependence of exhaled nitric oxide using a new
technique to exclude nasal nitric oxide. Am J Respir Crit Care Med, 1997. 155(1):
p. 260-7.
8. ATS, Recommendations for standardized procedures for the on-line and off-line
measurement of exhaled lower respiratory nitric oxide and nasal nitric oxide in
adults and children-1999. This official statement of the American Thoracic
Society was adopted by the ATS Board of Directors, July 1999. Am J Respir Crit
Care Med, 1999. 160(6): p. 2104-17.
9. ATS/ERS, ATS/ERS recommendations for standardized procedures for the online
and offline measurement of exhaled lower respiratory nitric oxide and nasal nitric
oxide, 2005. Am J Respir Crit Care Med, 2005. 171(8): p. 912-30.
10. George, S.C., How accurately should we estimate the anatomical source of
exhaled nitric oxide? J Appl Physiol, 2008. 104(4): p. 909-11.
11. Eckel, S.P. and M.T. Salam, Single high flow exhaled nitric oxide is an imperfect
proxy for distal nitric oxide. Occup Environ Med, 2013. 70(7): p. 519-20.
12. Horvath, I., et al., A European Respiratory Society technical standard: exhaled
biomarkers in lung disease. Eur Respir J, 2017. 49(4).
184
13. George, S.C., et al., Modeling pulmonary nitric oxide exchange. J Appl Physiol,
2004. 96(3): p. 831-9.
14. Lehtimä ki, L., T. Karvonen, and M. Hö gman, Clinical values of nitric oxide
parameters from the respiratory system. Current Medicinal Chemistry, 2020.
27(42): p. 7189-7199.
15. SC, G., et al., Modeling pulmonary nitric oxide exchange. Journal of applied
physiology (Bethesda, Md. : 1985), 2004. 96(3).
16. Eckel, S.P., et al., Estimation of parameters in the two-compartment model for
exhaled nitric oxide. PLoS ONE, 2014. 9(1): p. e85471.
17. Silkoff, P.E., et al., Airway nitric oxide diffusion in asthma: Role in pulmonary
function and bronchial responsiveness. Am J Respir Crit Care Med, 2000. 161(4
Pt 1): p. 1218-28.
18. Hogman, M. and P. Merilainen, Extended NO analysis in asthma. J Breath Res,
2007. 1(2): p. 024001.
19. Hogman, M., Extended NO analysis in health and disease. J Breath Res, 2012.
6(4): p. 047103.
20. Modig, L., et al., Effects of short-term exposure to air pollution on the levels of
exhaled nitric oxide among adults – Results from the ADONIX study in
Gothenburg, Sweden. Eur Respir J, 2012. 40(s56): p. 609s.
21. Rosa, M.J., et al., Association of recent exposure to ambient metals on fractional
exhaled nitric oxide in 9–11year old inner-city children. Nitric Oxide, 2014. 40:
p. 60-66.
22. Eckel, S., et al., Traffic-related air pollution and alveolar nitric oxide in southern
California children. Eur Respir J, 2016. 47(5): p. 1348-56.
23. Lopez, V., et al., Natural pollen exposure increases the response plateau to
adenosine 5'-monophosphate and bronchial but not alveolar nitric oxide in
sensitized subjects. Respiration, 2011. 83(3): p. 225-232.
24. Eckel, S., et al., Statistical Methods For Studying Determinants Of Airway And
Alveolar Nitric Oxide. Am J Respir Crit Care Med, 2013. 187: p. A5092.
25. Davidian, M. and D.M. Giltinan, Nonlinear models for repeated measurement
data: an overview and update. J Agric Biol Envir S, 2003. 8(4): p. 387-419.
26. Casella, G. and E.I. George, Explaining the Gibbs sampler. The American
Statistician, 1992. 46(3): p. 167-174.
27. Gauderman, W.J., et al., The effect of air pollution on lung development from 10
to 18 years of age. N Engl J Med, 2004. 351(11): p. 1057-67.
28. McConnell, R., et al., Traffic, susceptibility, and childhood asthma. Environ
Health Perspect, 2006. 114(5): p. 766-72.
185
29. Chen, Z., et al., Chronic effects of air pollution on respiratory health in Southern
California children: findings from the Southern California Children’s Health
Study. J Thorac Dis, 2015. 7(1): p. 46.
30. Salam, M.T., T. Islam, and F.D. Gilliland, Recent evidence for adverse effects of
residential proximity to traffic sources on asthma. Curr Opin Pulm Med, 2008.
14(1): p. 3-8.
31. Health Effects Institute. Panel on the Health Effects of Traffic-Related Air
Pollution, Traffic-related air pollution: a critical review of the literature on
emissions, exposure, and health effects. 2010: Health Effects Institute.
32. Brook, R.D., et al., Particulate matter air pollution and cardiovascular disease.
Circulation, 2010. 121(21): p. 2331-2378.
33. Benbrahim-Tallaa, L., et al., Carcinogenicity of diesel-engine and gasoline-
engine exhausts and some nitroarenes. Lancet Oncol, 2012. 13(7): p. 663-664.
34. Anderson, H.R., G. Favarato, and R.W. Atkinson, Long-term exposure to air
pollution and the incidence of asthma: meta-analysis of cohort studies. Air Qual
Atmos Health, 2013. 6(1): p. 47-56.
35. Kelly, F.J. and J.C. Fussell, Air pollution and public health: emerging hazards
and improved understanding of risk. Environ Geochem Health, 2015. 37(4): p.
631-649.
36. Schultz, E.S., A.A. Litonjua, and E. Melé n, Effects of long-term exposure to
traffic-related air pollution on lung function in children. Curr Allergy Asthma
Rep, 2017. 17(6): p. 41.
37. Tsoukias, N.M. and S.C. George, A two-compartment model of pulmonary nitric
oxide exchange dynamics. J Appl Physiol, 1998. 85(2): p. 653-66.
38. Hogman, M., et al., Extended NO analysis applied to patients with COPD,
allergic asthma and allergic rhinitis. Respir Med, 2002. 96(1): p. 24-30.
39. Molshatski, N. and S.P. Eckel, Optimal flow rate sampling designs for studies
with extended exhaled nitric oxide analysis. J Breath Res, 2017. 11(1): p. 016012.
40. Gelman, A., F. Bois, and J. Jiang, Physiological pharmacokinetic analysis using
population modeling and informative prior distributions. Journal of the American
Statistical Association, 1996. 91(436): p. 1400-1412.
41. Eckel, S.P., et al., Traffic-related air pollution and alveolar nitric oxide in
southern California children. European Respiratory Journal, 2016.
42. Plummer, M. JAGS: A program for analysis of Bayesian graphical models using
Gibbs sampling. in Proceedings of the 3rd international workshop on distributed
statistical computing. 2003. Vienna, Austria.
43. Brooks, S.P. and A. Gelman, General methods for monitoring convergence of
iterative simulations. Journal of computational and graphical statistics, 1998. 7(4):
p. 434-455.
186
44. Tsoukias, N.M., et al., Single-exhalation profiles of NO and CO2 in humans:
effect of dynamically changing flow rate. J Appl Physiol, 1998. 85(2): p. 642-52.
45. Jain, N. and J.L. Hill, The National Heart Lung and Blood Institute guidelines
finally say "yes" to fractional exhaled nitric oxide. Ann Allergy Asthma Immunol,
2022.
46. Khatri, S.B., et al., Use of Fractional Exhaled Nitric Oxide to Guide the
Treatment of Asthma: An Official American Thoracic Society Clinical Practice
Guideline. Am J Respir Crit Care Med, 2021. 204(10): p. e97-e109.
47. Dweik, R.A., et al., An official ATS clinical practice guideline: interpretation of
exhaled nitric oxide levels (FENO) for clinical applications. Am J Respir Crit
Care Med, 2011. 184(5): p. 602-15.
48. La Grutta, S., et al., Environmental Effects on Fractional Exhaled Nitric Oxide in
Allergic Children. Journal of Allergy, 2011. 2012.
49. Scarpa, M.C., N. Kulkarni, and P. Maestrelli, The role of non-invasive biomarkers
in detecting acute respiratory effects of traffic-related air pollution. Clin Exp
Allergy, 2014. 44(9): p. 1100-18.
50. Oliver, A.C., et al., Effects of Reduced Nicotine Content Cigarettes on Fractional
Exhaled Nitric Oxide and Self-Reported Respiratory Health Outcomes Among
Smokers With Psychiatric Conditions or Socioeconomic Disadvantage. Nicotine
Tob Res, 2022. 24(1): p. 135-140.
51. Society, A.T., European Respiratory Society. ATS/ERS recommendations for
standardized procedures for the online and offline measurement of exhaled lower
respiratory nitric oxide and nasal nitric oxide, 2005. Am J Respir Crit Care Med,
2005. 171: p. 912-930.
52. Roy, K., et al., Use of different exhaled nitric oxide multiple flow rate models in
COPD. Eur Respir J, 2007. 29(4): p. 651-9.
53. Karvonen, T., et al., Comparison of feasibility and estimates of central and
peripheral nitric oxide parameters by different mathematical models. J Breath
Res, 2017. 11(4): p. 047102.
54. Weng, J., et al., Hierarchical Bayesian estimation of covariate effects on airway
and alveolar nitric oxide. Sci Rep, 2021. 11(1): p. 17180.
55. Karrasch, S., et al., Accuracy of FENO for diagnosing asthma: a systematic
review. Thorax, 2017. 72(2): p. 109-116.
56. HL, P., et al., Exhaled nitric oxide levels to guide treatment for adults with
asthma. The Cochrane database of systematic reviews, 2016. 9(9).
57. Silkoff, P.E., et al., Longitudinal stability of asthma characteristics and
biomarkers from the Airways Disease Endotyping for Personalized Therapeutics
(ADEPT) study. Respir Res, 2016. 17: p. 43.
187
58. Kharitonov, S.A., et al., Reproducibility of exhaled nitric oxide measurements in
healthy and asthmatic adults and children. Eur Respir J, 2003. 21(3): p. 433-8.
59. Garcia, E., et al., Patterns and determinants of exhaled nitric oxide trajectories in
schoolchildren over a 7-year period. Eur Respir J, 2020. 56(1).
60. Linn, W.S., et al., Exhaled nitric oxide in a population-based study of Southern
California Schoolchildren, in Respir Res. 2009. p. 28.
61. Hö gman, M. and P. Merilä inen, Extended NO analysis in asthma. J Breath Res,
2007. 1(2): p. 024001.
62. Hö gman, M., et al., Extended NO analysis applied to patients with COPD,
allergic asthma and allergic rhinitis. Respir Med, 2002. 96(1): p. 24-30.
63. Latif, A.H. and S.G. Gilmour, Transform-both-sides nonlinear models for in vitro
pharmacokinetic experiments. Stat Methods Med Res, 2015. 24(3): p. 306-24.
64. Pinheiro, J., et al., Nlme: Linear and Nonlinear Mixed Effects Models. R package
version 31-110, 2013. 3: p. 1-113.
65. Plummer, M., JAGS: A program for analysis of Bayesian graphical models using
Gibbs sampling. undefined, 2003.
66. Plummer, M. JAGS: A program for analysis of Bayesian graphical models using
Gibbs sampling. 2003.
67. McConnell, R., et al., Traffic, susceptibility, and childhood asthma. Environ
Health Perspect, 2006. 114(5): p. 766-72.
68. Linn, W.S., et al., Multiple-flow exhaled nitric oxide, allergy, and asthma in a
population of older children. Pediatr Pulmonol, 2013. 48(9): p. 885-96.
69. Linn, W.S., et al., Extended exhaled nitric oxide analysis in field surveys of
schoolchildren: a pilot test. Pediatr Pulmonol, 2009. 44(10): p. 1033-42.
70. Puckett, J.L., et al., Impact of analysis interval on the multiple exhalation flow
technique to partition exhaled nitric oxide. Pediatr Pulmonol, 2010. 45(2): p. 182-
91.
71. Garcia, E., et al., Patterns and determinants of exhaled nitric oxide trajectories in
schoolchildren over a 7-year period. 2020.
72. Hö gman, M., et al., Effects of growth and aging on the reference values of
pulmonary nitric oxide dynamics in healthy subjects - IOPscience. 2017.
73. Louis, R., et al., European Respiratory Society Guidelines for the Diagnosis of
Asthma in Adults. European Respiratory Journal, 2022.
74. Khatri, S.B., et al., Use of Fractional Exhaled Nitric Oxide to Guide the
Treatment of Asthma: An Official American Thoracic Society Clinical Practice
Guideline. American journal of respiratory and critical care medicine, 2021.
204(10): p. e97-e109.
188
75. Chen, X., et al., The association between short-term exposure to ambient air
pollution and fractional exhaled nitric oxide level: A systematic review and meta-
analysis of panel studies. Environmental Pollution, 2020. 265: p. 114833.
76. Horvá th, I., et al., A European Respiratory Society technical standard: exhaled
biomarkers in lung disease. Eur Respir J, 2017. 49(4): p. 1600965.
77. Condorelli, P., et al., A simple technique to characterize proximal and peripheral
nitric oxide exchange using constant flow exhalations and an axial diffusion
model. J Appl Physiol, 2007. 102(1): p. 417-25.
78. Suresh, V., et al., Effect of heterogeneous ventilation and nitric oxide production
on exhaled nitric oxide profiles. J Appl Physiol, 2008. 104(6): p. 1743-52.
79. Roy, K., et al., Use of different exhaled nitric oxide multiple flow rate models in
COPD. Eur Respir J, 2007. 29(4): p. 651.
80. Karvonen, T., et al., Comparison of feasibility and estimates of central and
peripheral nitric oxide parameters by different mathematical models. Journal of
Breath Research, 2017. 11(4): p. 047102.
81. Weng, J., et al., Hierarchical Bayesian estimation of covariate effects on airway
and alveolar nitric oxide. Scientific reports, 2021. 11(1): p. 1-12.
82. Pietropaoli, A.P., et al., Simultaneous measurement of nitric oxide production by
conducting and alveolar airways of humans. J Appl Physiol, 1999. 87(4): p. 1532-
42.
83. Karvonen, T. and L. Lehtimä ki, Effect of exhalation flow rates and level of nitric
oxide output on accuracy of linear approximation of pulmonary nitric oxide
dynamics. Journal of Breath Research, 2021. 15(3): p. 036003.
84. Pinheiro, J., et al., nlme: Linear and Nonlinear Mixed Effects Models. R package
version 3.1-118. R Foundation for Statistical Computing, Vienna, 2014.
85. Zhang, Y., et al., Long-term exposures to air pollutants affect FeNO in children: a
longitudinal study. European Respiratory Journal, 2021. 58(5).
86. Garcia, E., et al., Patterns and determinants of exhaled nitric oxide trajectories in
schoolchildren over a 7-year period. European Respiratory Journal, 2020. 56(1).
87. Chen, Z., et al., Chronic effects of air pollution on respiratory health in Southern
California children: findings from the Southern California Children’s Health
Study. J Thorac Dis, 2015. 7(1): p. 46-58.
88. Kiang, T.K., et al., Fundamentals of population pharmacokinetic modelling.
Clinical pharmacokinetics, 2012. 51(8): p. 515-525.
89. Racine-Poon, A. and J. Wakefield, Statistical methods for population
pharmacokinetic modelling. Statistical Methods in Medical Research, 1998. 7(1):
p. 63-84.
189
90. Roe, D.J., et al., Comparison of population pharmacokinetic modeling methods
using simulated data: results from the population modeling workgroup. Statistics
in medicine, 1997. 16(11): p. 1241-1262.
91. Mould, D. and R.N. Upton, Basic concepts in population modeling, simulation,
and model‐based drug development—part 2: introduction to pharmacokinetic
modeling methods. CPT: pharmacometrics & systems pharmacology, 2013. 2(4):
p. 1-14.
92. Wakefield, J., The Bayesian analysis of population pharmacokinetic models.
Journal of the American Statistical Association, 1996. 91(433): p. 62-75.
93. Thornadtsson, A., et al., Extended nitric oxide analysis may improve personalized
anti-inflammatory treatment in asthmatic children with intermediate FENO50.
Journal of Breath Research, 2015. 9(4): p. 047114.
94. Asano, T., et al., Small airway inflammation is associated with residual airway
hyperresponsiveness in Th2-high asthma. Journal of Asthma, 2020. 57(9): p. 933-
941.
95. Sardó n, O., et al., Alveolar nitric oxide and its role in pediatric asthma control
assessment. BMC Pulmonary Medicine, 2014. 14(1): p. 1-7.
96. Linkosalo, L., et al., Relation of bronchial and alveolar nitric oxide to exercise‐
induced bronchoconstriction in atopic children and adolescents. Pediatric allergy
and immunology, 2012. 23(4): p. 360-366.
97. Lehtimaki, L., et al., Extended exhaled NO measurement differentiates between
alveolar and bronchial inflammation. American journal of respiratory and critical
care medicine, 2001. 163(7): p. 1557-1561.
98. Karvonen, T., et al., Onset of action of inhaled glucocorticoids on bronchial and
alveolar nitric oxide output. Journal of Breath Research, 2020. 15(1): p. 016008.
99. Nicolini, G., et al. Both bronchial and alveolar exhaled nitric oxide are reduced
with extrafine beclomethasone dipropionate in asthma. in Allergy and Asthma
Proceedings. 2010. OceanSide Publications.
100. Karampitsakos, T., et al., The effect of bronchodilation and spirometry on
fractional exhaled nitric oxide (FeNO50), bronchial NO flux (JawNO) and
alveolar NO concentration (CANO) in children and young adults. Journal of
Asthma, 2018. 55(8): p. 882-889.
101. Cattoni, I., et al., Mechanisms of decrease in fractional exhaled nitric oxide
during acute bronchoconstriction. Chest, 2013. 143(5): p. 1269-1276.
102. Mason, P., et al., Exhaled nitric oxide dynamics in asthmatic reactions induced by
diisocyanates. Clinical & Experimental Allergy, 2016. 46(12): p. 1531-1539.
103. Lehtimä ki, L., et al., Inhaled fluticasone decreases bronchial but not alveolar
nitric oxide output in asthma. European Respiratory Journal, 2001. 18(4): p. 635-
639.
190
104. Akamatsu, T., et al., Switching from salmeterol/fluticasone to
formoterol/budesonide combinations improves peripheral airway/alveolar
inflammation in asthma. Pulmonary Pharmacology & Therapeutics, 2014. 27(1):
p. 52-56.
105. Fritscher, L.G., et al., The effect of montelukast on exhaled nitric oxide of alveolar
and bronchial origin in inhaled corticosteroid-treated asthma. Respiratory
medicine, 2009. 103(2): p. 296-300.
106. Gelb, A.F., et al., Alveolar and airway sites of nitric oxide inflammation in treated
asthma. American journal of respiratory and critical care medicine, 2004. 170(7):
p. 737-741.
107. Dahgam, S., et al., Effects of Short-Term Exposure to Ozone on Alveolar Nitric
Oxide. Epidemiology, 2009. 20(6): p. S145-S146.
108. Rosa, M.J., et al., Association of recent exposure to ambient metals on fractional
exhaled nitric oxide in 9–11 year old inner-city children. Nitric Oxide, 2014. 40:
p. 60-66.
109. Rosa, M.J., et al., Fractional exhaled nitric oxide exchange parameters among 9‐
year‐old inner‐city children. Pediatric pulmonology, 2011. 46(1): p. 83-91.
1.
191
Chapter 5. Summary and future directions
Summary
Multiple flow FeNO is a promising technique for assessing airway and alveolar
inflammation but there is room for improvement in the methodology through the
development of appropriate statistical methods, particularly for analyses of multiple flow
FeNO data collected in population studies. My dissertation introduced novel unified
hierarchical Bayesian (HB) models tailored to this problem. The HB model introduced in
Chapter 2 for cross-sectional studies is the base for several extensions, including the
longitudinal and two-stage versions in Chapters 3 and 4. Unified HB approaches
outperformed the two-stage methods currently in wide use as well as unified frequentist
methods (U-NLME and L-U-NLME) which, as implemented in the standard R package,
often suffered from convergence issues.
Unified UB methods, as other Bayesian simulation approaches, consumed large
computation resource and required much longer computation time. In our previous cross-
sectional simulation study of N=1000 (participants) with correlated NO parameters, U-
HB took 5.5 hours to run on median (Chapter 3). We tested several other simulated data
scenarios and found that U-HB converged slightly more quickly when the NO parameters
were not correlated or when datasets were smaller. However, with considerably smaller
192
datasets (e.g., N=100), sometimes the U-HB model took longer to converge. In
longitudinal datasets of N=500, M=3 (visit times) with correlated NO parameters took 30
hours to run L-U-HB, on median (Chapter 4). We also found that the L-U-HB converged
much more quickly when the NO parameters were only correlated in either participant or
visit level, or when the datasets were smaller in either N or M.
The key contributions of this dissertation were the following:
Constrained multivariate normal distribution in JAGS. NO parameters have biological
interpretations and hence a range of reasonable values. To better satisfy distributional
assumptions and to impose a positivity constraint, 𝐶 𝑎𝑤
𝑁𝑂 and 𝐷 𝑎𝑤
𝑁𝑂 were directly
modeled as their natural log transformation and exponentiated in the nonlinear 2CM
function. In preliminary data, 𝐶 𝐴 𝑁𝑂 did not have a lognormal distribution but more of a
truncated normal distribution with a biological requirement of a positivity constraint,
since concentrations are non-negative. Hence in HB models, the participant-level vector
of three NO parameters was modeled using a truncated multivariate normal distribution
with a non-negative constraint on 𝐶 𝐴 𝑁𝑂 . The JAGS software did not permit truncated
multivariate normal distributions, so to solve this problem I separately sampled 𝐶 𝐴 𝑁𝑂
from a truncated univariate normal distribution conditional on log 𝐶 𝑎𝑤
𝑁𝑂 and
193
log 𝐷 𝑎𝑤
𝑁𝑂 sampled from a bivariate normal distribution. The prior distributions for the
means, standard deviations, and correlations for 𝛽 𝐶 𝐴 , 𝛽 𝑙𝑜𝑔 𝐶𝑎𝑤 , 𝛽 𝑙𝑜𝑔 𝐷𝑎𝑤
were sampled
independently. However, a variance-covariance matrix with more than two variables
constructed by randomly sampled standard deviations and correlations can still
potentially fail the positive definite requirement. The positive definitive variance-
covariance structure was only affected by the correlations, so I first sampled the
correlations 𝜌 𝐶 𝐴 ,𝑙𝑜𝑔 𝐶𝑎𝑤 , 𝜌 𝐶 𝐴 ,𝑙𝑜𝑔 𝐷𝑎𝑤
and then sampled the third correlation 𝜌 𝑙𝑜𝑔 𝐶𝑎𝑤 ,𝑙𝑜𝑔 𝐷𝑎𝑤
within the range required to fulfill the positive definitive requirement conditioned on the
sampled other two. (Appendix ##)
Correlated NO parameters. Biologically, there are likely correlations between NO
parameters at the participant level and at the visit level within participant. For example,
participants with high 𝐶 𝑎𝑤
𝑁𝑂 may also tend to have higher 𝐶 𝐴 𝑁𝑂 as the concentrations
of NO in airway and alveolar region are likely not completely independent. These
correlations also relate to the issue of positive definiteness because if there are no
correlations the variance-covariance matrices for NO parameters in each level are
ensured to be positive definite. However, in the longitudinal data, the positive definite
problem can be avoided for the simplified model by assuming independence between NO
parameters in the visit level. Note that in population pharmacokinetics, a similar
194
assumption of within-participant independence in parameters is often made.[91] The non-
negative constraint on 𝐶 𝐴 𝑁𝑂 was only specified in the visit level, which was assumed to
have independent variance matrix that is therefore already positive definite. Also, the
variance-covariance matrix for the subject level was only used to generate the participant
randomness and included no non-negative constraints.
Another problem related to correlation is convergence issues. For example, a longitudinal
HB model fitted to simulated data for N=500 participants, with 3 visits per participant,
converged in ~5 hours assuming NO parameters were independent at both levels and
converged in more than 24 hours if there were only visit-level correlations or only
participant-level correlations. However, after 4-5 days it failed to converge when
correlations were permitted at both levels. To accelerate the convergence, we found that
using a long adaptation phase and carefully specified initial values for both covariates
associations and variance-covariance matrices were the most helpful strategies. In future
studies, we wish to use a complete model with no simplification of the NO parameter
correlation structure for longitudinal study designs, which might be facilitated by
considering simulations with more repeated visits per participant.
195
Although NO parameter correlations brought technical challenges to our modelling
strategies, it also revealed interesting information. In Chapters 2 - 4, we focused on
Simulation Scenario 1, which assumed all NO parameters were affected by the covariate
and only discussed the other simulation scenarios in the context of type I errors. In the
supplementary materials for Chapters 2 and 3, results were presented for all scenarios.
Unified HB methods proved to be the most stable across all scenarios, but other methods
were more sensitive to scenarios in which only select NO parameters were affected by the
covariate, likely due to correlation of NO parameters. Other explanations for differences
in performance across the methods include whether uncertainty in Stage I estimates of
NO parameters is appropriately propagated to Stage II.
Extensive simulation studies evaluating the performance of several methods. The
simulation studies of Chapter 2 explored many different scenarios with different
combinations of covariate coefficients (i.e., one scenario had the covariate affect all NO
parameters, other scenarios had the covariate affect only one or two of the three NO
parameters). We did this to mimic the way that different regions in the respiratory system
might be affected by certain diseases or drugs. In Chapter 3, due to the increased
computation time for longitudinal data, we shank the sample size from 1000 participants
to 500 and only considered 3 visits. However, only having 3 visits will typically result in
196
more uncertainty in estimates of the visit-level correlation of NO parameters.
Furthermore, we found that some U-HB models required more than 5 days to converge.
Due to the 48-hour time limit on USC’s high performance computing cluster, we
simplified the Longitudinal U-HB model to assume no visit-level correlation between NO
parameters beyond the participant level correlations. We also reduced the number of
scenarios, only considering Scenario 1 where all NO parameters were affected and
scenarios where only one of the NO parameters was affected. In future studies, with more
efficient algorithms, we could test longitudinal U-HB models on simulated data sets with
larger sample sizes or more visits and allowing for correlation of NO parameters at both
the participant and visit levels.
Two stage versus unified approaches.
Two stage methods suffer in studies with sparse data for Stage I (e.g., a small number of
FeNO maneuvers per participant), ignore potential differences in the uncertainty of Stage
I NO parameter estimates across participants, and fail to appropriately account for the
correlation of NO parameter estimates. The unified methods considered in this
dissertation overcome these issues and are novel to the field of population studies of
multiple flow FeNO.
197
Bayesian versus Frequentist Approaches. The Bayesian approach has been increasingly
popular in many areas of application. It provides an alternative view that assumes all
parameters are random variable whose behavior is best captured by a probability
distribution rather than being fixed, single values that are to be estimated. This allows
them to describe uncertainty more fully in parameter estimation, as captured by the
resulting parameter posterior distributions. The Bayesian approach typically imposes a
high computational burden, and as such its popularity has increased in recent decades as
the growth of computational ‘horsepower’ has made it more feasible to fit such models.
A unified Bayesian approach such as our U-HB method further increases the
computational burden, since everything is being estimated simultaneously, but in doing
so allows us to fully characterize and capture parameter uncertainty, and correlations
between parameters. As such we consider this extra computational time to be well spent.
While run-times were indeed long, given enough iterations, convergence was generally
obtained.
Future work and applications
This dissertation has demonstrated U-HB to be the least biased and most powerful
approach, at the cost of computational time. A major area of future work would be to
investigate alternative algorithms to implement the U-HB method. In early pilot studies, I
198
had implemented U-HB models using both Metropolis-Hasting algorithms with the R
package “mcmc” and Gibbs sampling using JAGS, and found that JAGS converged much
faster and had better flexibility to specify additional features like the non-negative
constraint for CANO and to modify the correlation structure for NO parameters. The
implementation of HB models in this dissertation uses JAGS, however there are other
alternatives using off-the-shelf software like Rstan or NIMBLE or custom-written
algorithms. Efficiency might also be improved using parallel MCMC. In this work, I used
conjugate prior distributions for both NO parameters and covariates. The normal
conjugate prior distributions of NO parameters can be replaced by other conjugate priors
like beta or Poisson. Alternative software could also be adapted from the population
pharmacokinetics field. For example, NONMEM, the standard software for
pharmacokinetics and pharmacodynamic data, can deal with nonlinear mixed effect
models with varieties of frequentist and Bayesian approaches but is expensive, not open
source and has a limited set of potential prior distributions. Although this dissertation is
focused on the 2CM for multiple flow FeNO, the methods and lessons learned are also
applicable to other pharmacokinetic phenomenon described with other closed form,
nonlinear deterministic models. In summary, my work has demonstrated the potential of
applying unified estimation methods, especially using a Bayesian perspective, to data
analysis fields that traditionally rely on two-stage frequentist approaches.
199
Asset Metadata
Creator
Weng, Jingying (author)
Core Title
Bayesian models for a respiratory biomarker with an underlying deterministic model in population research
Contributor
Electronically uploaded by the author
(provenance)
School
Keck School of Medicine
Degree
Doctor of Philosophy
Degree Program
Biostatistics
Degree Conferral Date
2022-08
Publication Date
07/26/2022
Defense Date
05/16/2022
Publisher
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
Air pollution,airway inflammation,asthma,Bayesian,FeNO,Gibbs sampling,hierarchical model,JAGS,MCMC,multilevel model,OAI-PMH Harvest
Format
application/pdf
(imt)
Language
English
Advisor
Eckel, Sandrah (
committee chair
), D'Argenio, David (
committee member
), Gauderman, William (
committee member
), Gilliland, Frank (
committee member
), Marjoram, Paul (
committee member
)
Creator Email
wengjing@usc.edu,wengjingying@gmail.com
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-oUC111375229
Unique identifier
UC111375229
Legacy Identifier
etd-WengJingyi-11002
Document Type
Dissertation
Format
application/pdf (imt)
Rights
Weng, Jingying
Type
texts
Source
20220728-usctheses-batch-962
(batch),
University of Southern California
(contributing entity),
University of Southern California Dissertations and Theses
(collection)
Access Conditions
The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the author, as the original true and official version of the work, but does not grant the reader permission to use the work if the desired use is covered by copyright. It is the author, as rights holder, who must provide use permission if such use is covered by copyright. The original signature page accompanying the original submission of the work to the USC Libraries is retained by the USC Libraries and a copy of it may be obtained by authorized requesters contacting the repository e-mail address given.
Repository Name
University of Southern California Digital Library
Repository Location
USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA
Repository Email
uscdl@usc.edu
Abstract (if available)
Abstract
This dissertation is the culmination of three projects, two statistical and one focused on disseminating the technical results to the breath research community along with companion software. All projects are motivated by the multiple flow FeNO data collected in the Southern California Children’s Health Study (CHS), introduce this dataset more fully in Chapter 2. Chapter 1 introduces the background and motivations. In Chapter 2, we develop a novel Unified Hierarchical Bayesian (U-HB) method. We then conduct an extensive simulation study to evaluate its properties and compare it to alternative unified and two step approaches. Next, we then re-analyze cross-sectional data from the CHS to estimate the association of traffic-related air pollution exposure with alveolar and airway inflammation. In Chapter 3, we introduce a longitudinal extension of U-HB and demonstrate that HB methods outperform alternative methods for estimating associations of covariates with NO parameters. A manuscript summarizing this work is currently under review (submitted March 2022). Chapter 4 aims to provide a high-level overview of the technical work on cross-sectional and longitudinal methods in Chapters 2 and 3, targeting researchers conducting population studies with multiple flow FeNO. We review key statistical issues, provide an overview of relevant statistical methods while highlighting their strengths and weaknesses, and demonstrate how to implement these methods using R code that we have developed and made freely available on GitHub. Using cross sectional multiple flow data from the CHS, we apply six methods to estimate differences in NO parameters, comparing healthy children to those with asthma only, allergic rhinitis only, or both asthma and allergy. Chapter 5 summarizes the dissertation and discuss existing problems and future development of HB methods
Tags
airway inflammation
Bayesian
FeNO
Gibbs sampling
hierarchical model
JAGS
MCMC
multilevel model
Linked assets
University of Southern California Dissertations and Theses