# A Penalized Functional Linear Cox Regression Model for Spatially-defined
Environmental Exposure with an Estimated Buffer Distance

Jooyoung Lee Corresponding author. Email<EMAIL_ADDRESS>Department of
Applied Statistics, Chung-Ang University, Seoul Korea  Zhibing He  Department
of Biostatistics, Center for Methods on Implementation and Prevention Science,
CT USA  Department of Biostatistics, Yale School of Public Health, CT USA
Charlotte Roscoe  Channing Division of Network Medicine, Department of
Medicine, Brigham and Women’s Hospital and Harvard Medical School, MA USA
Department of Environmental Health, Harvard T.H. Chan School of Public Health,
MA USA  Division of Population Sciences, Dana Farber Cancer Institute, MA USA
Peter James  Division of Chronic Disease Research Across the Life Course
(CoRAL), Department of Population Medicine, Harvard Medical School and Harvard
Pilgrim Health Care Institute, MA USA  Department of Environmental Health,
Harvard T.H. Chan School of Public Health, MA USA  Li Xu  Department of
Epidemiology, Harvard T.H. Chan School of Public Health, MA USA  Donna
Spiegelman  Department of Biostatistics, Yale School of Public Health, CT USA
Center for Methods on Implementation and Prevention Science, Yale School of
Public Health, CT USA  David Zucker  Department of Statistics and Data
Science, Hebrew University, Jerusalem Israel  Molin Wang  Channing Division of
Network Medicine, Department of Medicine, Brigham and Women’s Hospital and
Harvard Medical School, MA USA  Department of Epidemiology, Harvard T.H. Chan
School of Public Health, MA USA  Department of Biostatistics, Harvard T.H.
Chan School of Public Health, MA USA

###### Abstract

In environmental health research, it is of interest to understand the effect
of the neighborhood environment on health. Researchers have shown a protective
association between green space around a person’s residential address and
depression outcomes. In measuring exposure to green space, distance buffers
are often used. However, buffer distances differ across studies. Typically,
the buffer distance is determined by researchers _a priori_. It is unclear how
to identify an appropriate buffer distance for exposure assessment. To address
geographic uncertainty problem for exposure assessment, we present a domain
selection algorithm based on the penalized functional linear Cox regression
model. The theoretical properties of our proposed method are studied and
simulation studies are conducted to evaluate finite sample performances of our
method. The proposed method is illustrated in a study of associations of green
space exposure with depression and/or antidepressant use in the Nurses’ Health
Study.

Keywords: Functional linear Cox regression model; B-spline basis functions;
Group bridge Lasso; neighborhood environment; buffer distance.

##  1 Introduction

Spatial uncertainty in environmental exposure assessment is a critically
important issue in environmental epidemiology, as it can bias health effect
estimates (Spiegelman, 2010). In environmental exposure assessment,
researchers often query spatial datasets using a geographic information system
(GIS) to quantify the environment surrounding residential addresses, with the
objective of assessing associations of environmental exposure and health
across many residents. To assess spatial exposures, spatial data are typically
quantified (e.g., spatially averaged) based on circular distance buffers
surrounding addresses, which are created by specifying multiple distances
(radii) and drawing concentric circles at the specified distances around the
address points. Typically, the radii of circular buffers are selected by
researchers _a priori_ based on hypothesized mechanistic pathways. For
example, a larger buffer radius compared to a smaller buffer radius may be
deemed appropriate to capture the environment relevant for physical activity,
whereas the reverse may be hypothesized if the researcher aims to assess
acoustic properties directly around participants’ homes. Selection of buffers
_a priori_ , however, may lead to inconsistency in exposure assessment across
epidemiological studies. If we use a smaller buffer radius than the true
buffer radius, we ignore the contribution of environmental exposures in the
effect estimate outside the buffer radius and bias may occur since
environmental exposures in the near distance are correlated. On the other
hand, if we use a larger buffer radius than the correct buffer radius, then
the variance of the effect estimate may increase (Baek et al., 2016). In
addition, the effect estimate may be null in the intermediate area, and the
traditional approach does not work in this setting. Therefore, spatial
uncertainty in exposure assessment may hinder accurate inference of
associations with health outcomes, as well as compromise comparability across
studies. These concerns, in turn, may affect the translation of research
findings into actionable environmental policy. Kwan (2012) defined this issue
as the Uncertain Geographic Context Problem, or “the spatial uncertainty in
the actual areas that exert the contextual influences under study.” In the
context of environmental epidemiology, this can be interpreted as uncertainty
in selecting the correct buffer surrounding residential addresses to observe
the influence of an environmental exposure on a health outcome. Previous
studies have used multiple discrete buffer radii for spatial environmental
exposures and shown that associations between exposure and health measures are
sensitive to different buffer radii (James et al., 2014). However, finding the
limit beyond which the environment has no influence on health is difficult
since the additional effect of exposure by increasing buffer radius cannot be
estimated.

In this paper, we measure environmental exposures in concentric rings with
multiple radii, investigate the association between exposure on health
outcomes as a function of radius, and identify a radial distance beyond which
no association between exposure and health is observed. The environmental
exposures can be statistically expressed using functional covariates as a
function of radius within the framework of functional regression models. Then,
finding areas in which environmental exposure is associated with the health
outcome can be translated into the problem of identifying the domain of
functional covariates in which they may affect the response in the functional
regression models (James et al., 2009). In other words, by solving the domain
selection problem in the functional regression models with survival data, we
can determine relevant areas for a particular exposure-health association.

Functional linear regression models were originally developed by Ramsay &
Silverman (1997) to address the problem of estimating regression coefficients
when covariates are functions. Functional linear Cox regression models have
been recently proposed to address the relationship between functional
covariates and survival outcomes. Gellar et al. (2015) maximized the penalized
partial likelihood using B-splines smooth functions and Qu et al. (2016)
formulated the functional Cox regression model through a reproducing kernel
Hilbert space. Kong et al. (2018) adopted the functional principle component
analysis framework in the time-to-event settings. To our knowledge, however,
the domain selection problem has not been addressed in functional linear Cox
regression models. In functional linear regression models, a method for
obtaining smooth and locally sparse estimators based on the penalized least
squares method was developed by Lin et al. (2017) to identify the null regions
where the coefficient functions are zero; that is, the regions where no
association of exposure and outcome is observed. Additionally, to find a
truncation point in a functional linear regression model, Hall & Hooker (2016)
proposed a penalized least squares method with a penalty on a truncation point
and Guan et al. (2020) adopted a nested group bridge penalty to estimate a
truncation point, which can be used to identify the optimal/inflection point
across exposures. While the extension of Lin et al. (2017) to the Cox
regression model can be implemented to find regions with null-effect, their
method relies on the value of a pre-determined threshold to select variables
based on the local quadratic approximation. The method of Guan et al. (2020)
is restricted to the special case, in which the coefficient functions can be
zero only in the tail region of the distribution, and uses a nested penalty
structure. However, it is often inadequate to impose increasing penalties as
the buffer distance increases in environmental exposure applications.

In this paper, we develop methods for identifying the areas corresponding to
non-null regions of coefficient functions in survival data analysis settings.
We build on a combination of the functional Cox model literature and the
sparse functional data analysis literature and use a more flexible method for
allowing for regions where the coefficient function is equal to zero. We
consider a penalized scalar-on-function linear Cox regression model to
estimate the coefficient functions, as well as to select important domains. We
use a B-spline basis to represent the coefficient function. We adopt a
smoothness penalty to enforce smoothness of the coefficient function and a
group bridge penalty for the selection of spline coefficients (Huang et al.,
2009, 2014). Our aim is to address geographic uncertainty via data-informed
identification of optimal exposure areas as well as to estimate the
association of environmental exposure and survival outcomes within these non-
null regions. In addition, by investigating how exposure-health association
changes with the distance from the residence, it is possible to find an area
where the strongest association is observed. Our proposed method is, in
general, to identify non-null regions. However, it can be also applied to the
special case where the objective is to identify a certain location after which
there is no association by identifying the largest non-null location.

This article is organized as follows. In Section 2, we describe the motivating
data set used in our case study including the green space exposure and the
Nurses’ Health Study. In Section 3, we present the proposed methods for
functional linear Cox regression models with the smooth and group bridge
penalties. In addition, we derive the theoretical properties of the proposed
estimators. In Section 4, we examine the performance of our method using
simulation studies. In Section 5, we apply our method to an illustrative
example to examine the association between greenspace exposure and depression
incidence in the Nurses’ Health Study (NHS). We provide concluding remarks in
Section 6.

##  2 Motivating Data Set

Depressive disorders have prevailed as one of the leading causes of
disability, globally, for the last 30 years (James et al., 2018). They affect
more than 264 million people worldwide (World Health Organization, 2017).
Biological predisposition and psychosocial factors can contribute to
depression onset (World Health Organization, 2017); however, modifiable
environmental exposures may also play a role in increasing risk for, or
protecting against, depression (van den Bosch & Meyer-Lindenberg, 2019).
Researchers have shown a protective association between urban greenspace
surrounding the residential addresses and depression or depressive symptoms
across multiple contexts (Sarkar et al., 2018; Pun et al., 2018; Liu et al.,
2019).

Our investigation is motivated by a study by Banay et al. (2019), which showed
that women in the Nurses’ Health Study (NHS) who lived in the most-green
versus least green exposure quintile showed a lower risk of clinically-
diagnosed depression and/or antidepressant use. The authors used two distance
buffers (radii: 250m, 1250m) to assess greenness exposure, which were selected
_a priori_ , and which, as the authors acknowledged, may introduce bias due to
geographic uncertainty. Statistically significant associations between
greenspace and depression were observed in the smaller buffer (250m), though
not the larger buffer (1250m); however, the limit at which the exposure no
longer exerts an effect on the outcome is unknown. In the present paper, using
higher resolution greenness data, we aim to identify the radial distance
beyond which greenness exposure has no influence on depression incidence.

The NHS is a prospective cohort study of U.S. female nurses living across 11
states. It was established in 1976, when a total of 121,701 nurses aged 30-55
were recruited. Nurses receive biennial questionnaires on lifestyle risk
factors and health outcomes. In 2000, depression diagnosis statuses and
antidepressant prescription items were included on the questionnaire. The
current analysis included all women who returned questionnaires in 2000 and
had no diagnosis of depression before 2000. The 35,283 eligible participants
contributed 296,436 person-years of follow-up, and 3,250 incident depression
cases occurred between 2000 and 2010.

Figure 1: An example of Normalized Difference Vegetation Index (NDVI) data
surrounding an example residential address with multiple radial distance
buffers.

Normalized Difference Vegetation Index (NDVI) values were used to quantify
greenness surrounding residential addresses of NHS participants. NDVI is a
satellite-derived indicator of quantity of photosynthesizing vegetation; it is
widely used as a marker for exposure to greenspaces in epidemiological studies
(Rhew et al., 2011). NDVI data derived from Landsat-7 satellite imagery were
ascertained and processed via Google Earth Engine (GEE). We used the
LandsatSimpleComposite algorithm to compute a Landsat Top of Atmosphere
composite from 30m resolution raw scenes and selected the least-cloudy scenes.
We computed the normalized difference between the near-infrared and red bands
to calculate NDVI. From the 30m NDVI composite, we used the GEE
reduceNeighborhood algorithm to create 30m resolution surfaces for the
contiguous US in which each pixel represented the average of the NDVI values
within a radial buffer surrounding that pixel in the input image (radii: 90m,
150m, 270m, 510m, 750m, 990m, 1230m, 1500m or 2100m) in 2000 (Figure 1).
Pixels were included in the calculation of average NDVI if 50% of their
surface was within the radial buffer. NDVI values range from 0 to 1, with
higher values indicating greater cover of photosynthesizing vegetation. NDVI
values below 0 indicate water and were set to 0. We spatially joined this NDVI
dataset to the location of the residential addresses of NHS participants using
the raster package (https://rspatial.org/raster) in R. NDVI reaches its
maximum and highest level of geographic variation during the summer. We
therefore used the mean NDVI value of summer months (July 1 – September 30) in
2000 (study baseline) for exposure assessment.

##  3 Methods

###  3.1 Functional linear Cox regression model

We let $\tilde{T}_{i}$ be time to event of interest, $C_{i}$ the censoring
time for individual $i=1,\ldots,n$, and denote $\tau$ as the end of follow-up
time. We define $X_{i}(s)$ as the average of environmental exposure on a
concentric ring at radius $s$, centered on the residence for individual $i$.
We let $\mathcal{S}$ be a closed and bounded interval of $\mathbb{R}$ and
$\\{X_{i}(s):s\in\mathcal{S}\\}$ be a random function expressing the exposure.
We take $\mathcal{S}$ to be $[0,1]$ without loss of generality. We observe the
environmental exposures at the pre-specified radii of concentric rings
$\\{s_{r}\in\mathcal{S},1\leq r\leq\mathcal{R}\\}$, where $\mathcal{R}$ is the
number of radii of observed concentric rings. Let
$Z_{i}=(Z_{i1},\ldots,Z_{ip})^{T}$ be a $p$-vector of scalar covariates that
we adjust for, and define $W_{i}=(X_{i},Z_{i})^{T}$. Under the non-informative
right-censoring setting, we observe $n$ i.i.d. samples
$(T_{i},\Delta_{i},Z_{i},\\{X_{i}(s_{r}),1\leq r\leq\mathcal{R}\\})$ for
$i=1,\ldots,n$, where $T_{i}=\min(\tilde{T}_{i},C_{i})$, and
$\Delta_{i}=I(T_{i}\leq C_{i})$ is an indicator for occurrence of the event.
The functional linear Cox regression model specifies the hazard function as

$\lambda(t)=\lambda_{0}(t)\exp\left(\int_{0}^{1}X_{i}(s)\beta(s)ds+\theta^{T}Z_{i}\right),$
(1)

where $\lambda(t)=d\Lambda(t)$, $\lambda_{0}(t)=d\Lambda_{0}(t)$ is the
baseline hazard function, $\beta(s)$ represents the effect of average
environmental exposures $X(s)$ measured on the concentric rings with radius
$s$ on the survival outcome, and we denote
$\theta=(\theta_{1},\ldots,\theta_{p})^{T}$ for a vector of covariate effects.

Our approach utilizes the B-spline basis expansion of $\beta(s)$ (De Boor,
2001; James et al., 2009; Lin et al., 2017). We let $\mathcal{S}_{dM_{n}}$
denote the linear space of spline functions on $\mathcal{S}$ spanned by the
B-spline basis functions $\\{B_{j}(s):j=1,\ldots,M_{n}+d+1\\}$ with degree $d$
and $M_{n}$ equally spaced inner knots
$0=\kappa_{0}<\kappa_{1}<\cdots<\kappa_{M_{n}}<\kappa_{M_{n}+1}=1$, and let
$B(s)=(B_{1}(s),\ldots,B_{M_{n}+d+1}(s))^{T}$. Over each of subintervals
$[s_{m},s_{m+1}]$, each B-spline basis function $B_{m}(s)$ is a piecewise
polynomial of degree $d$ and nonzero over no more than $d+1$ adjacent
subintervals. Under suitable smoothness assumptions, we can approximate
$\beta(s)$ as a linear combination of B-spline basis functions
$\beta(s)\approx\beta^{*}(s)=\sum_{k=1}^{L_{n}}b_{k}B_{k}(s)=B^{T}(s)b$, where
$L_{n}=M_{n}+d+1$, $\beta^{*}(s)\in\mathcal{S}_{dM_{n}}$, and
$b=(b_{1},\ldots,b_{L_{n}})^{T}$ is a vector of coefficients.

###  3.2 Smoothness and sparsity of $\beta(s)$

The estimators of coefficient functions $\beta(s)$ can be wiggly when $M_{n}$
becomes large. To enforce smoothness on $\beta(s)$, we impose the smoothness
penalty $b^{T}Jb$ where
$J_{i,j}=\int_{0}^{1}\frac{\partial^{2}{B}_{i}(s)}{\partial
s^{2}}\frac{\partial^{2}{B}_{j}(s)}{\partial s^{2}}ds$ (Cardot et al., 2003;
Yuan & Cai, 2010).

In order to identify the buffer distance, we need a sparse representation of
$\beta(s)$. By obtaining the sparsity of $\beta(s)$, it is possible to
identify non-null regions $\mathcal{S}_{1}=\\{s:\beta(s)\neq
0,s\in\mathcal{S}\\}$ and null regions
$\mathcal{S}_{0}=\\{s:\beta(s)=0,s\in\mathcal{S}\\}$. By the B-spline
approximation, the sparsity of $\beta(s)$ can be characterized by the sparsity
of the coefficient vector $b$. Specifically, if $\beta(s)=0$ for
$s\in[\kappa_{j-1},\kappa_{j}]$, then $b_{j}=\cdots=b_{j+d}=0$ since the
B-spline basis functions have the local support property (Wang & Kai, 2015).
This overlapping group structure of coefficients of $b$ requires the selection
of groups of variables where the groups can be overlapped. In this regard, we
impose the group bridge penalty (Huang et al., 2009, 2014) to select at both
group-level and individual-level. Therefore, by estimating $b$, we can
estimate $\beta(s)$ and identify the buffering distance.

We estimate $\alpha=(b^{T},\theta^{T})^{T}$ by maximizing the penalized log
partial likelihood

$\mathcal{L}_{n}(\alpha)=n^{-1}l_{n}(\alpha)-\lambda_{1}\sum_{j=1}^{M_{n}+1}\|b_{A_{j}}\|_{1}^{\gamma}-\lambda_{2}b^{T}Jb,$
(2)

where $l_{n}(\alpha)$, the log partial likelihood, is given by

$\displaystyle l_{n}(\alpha)=$
$\displaystyle\sum_{i=1}^{n}\Delta_{i}\bigg{\\{}\sum_{k=1}^{L_{n}}b_{k}\int_{0}^{1}X_{i}(s)B_{k}(s)ds+\theta^{T}Z_{i}$
$\displaystyle-\log\sum_{T_{j}\geq
T_{i}}\exp\left(\sum_{k=1}^{L_{n}}b_{k}\int_{0}^{1}X_{j}(s)B_{k}(s)ds+\theta^{T}Z_{j}\right)\bigg{\\}},$
(3)

$A_{j}=\\{j,j+1,\ldots,j+d\\}$, $b_{A_{j}}=(b_{k},k\in A_{j})^{T}$, and
$\|b_{A_{j}}\|_{1}=|b_{j}|+\ldots+|b_{j+d}|$ for $0<\gamma<1$ (Wang & Kai,
2015). We note that $\lambda_{1}$ is a regularization parameter that controls
the local sparsity of $\beta(s)$, and $\lambda_{2}$ enforces smoothness on
$\beta(s)$. We denote the maximizers as $\hat{\theta}$ and $\hat{b}$, then the
spline estimator of $\beta$ is $\hat{\beta}(s)=B(s)^{T}\hat{b}$.

###  3.3 Computational algorithm

Since direct maximization of (2) is difficult due to the non-convexity of the
group bridge penalty for $0<\gamma<1$, we follow the approach of Huang et al.
(2009) to maximize (2). For $0<\gamma<1$, we let

$\mathcal{L}_{n}(\alpha)=n^{-1}l_{n}(\alpha)-\sum_{j=1}^{M_{n}+1}\mu_{j}^{1-1/\gamma}\|b_{A_{j}}\|_{1}-\zeta\sum_{j=1}^{M_{n}+1}\mu_{j}-\lambda_{2}b^{T}Jb,$

where
$\mu_{j}=\\{(1-\gamma)/(\gamma\zeta)\\}^{\gamma}\|b_{A_{j}}\|_{1}^{\gamma}$
with $\zeta=\lambda_{1}^{1/(1-\gamma)}\gamma^{\gamma/(1-\gamma)}(1-\gamma).$
Define $\tilde{l}_{n}(\alpha)=-l_{n}(\alpha)$ and denote the gradient vector
$\bigtriangledown\tilde{l}_{n}(\alpha)=\partial\tilde{l}_{n}(\alpha)/\partial\alpha$
and the Hessian matrix
$\bigtriangledown^{2}\tilde{l}_{n}(\alpha)=\partial^{2}\tilde{l}_{n}(\alpha)/\partial\alpha\partial\alpha^{T}$.
Consider the Cholesky decomposition of
$\bigtriangledown^{2}\tilde{l}_{n}(\alpha^{\circ}),\bigtriangledown^{2}\tilde{l}_{n}(\alpha^{\circ})=V^{T}V$
($V$ is an upper triangle matrix), and set the pseudo response vector
$Y=(V^{T})^{-1}\\{\bigtriangledown^{2}\tilde{l}_{n}(\alpha^{\circ})\alpha^{\circ}-\bigtriangledown\tilde{l}_{n}(\alpha^{\circ})\\}$.
Following the arguments in Zhang & Lu (2007) and Huang et al. (2014), the
partial log-likelihood can be approximated by a quadratic form using a Taylor
expansion around $\alpha^{\circ}$, i.e., $\tilde{l}_{n}(\alpha)\approx
1/2(Y-V\alpha)^{T}(Y-V\alpha)$ (see Section 1 of the Supplemental Materials
for the derivation). Therefore, we propose to minimize

$\frac{1}{2n}(Y-V\alpha)^{T}(Y-V\alpha)+\sum_{j=1}^{M_{n}+1}\mu_{j}^{1-1/\gamma}\|b_{A_{j}}\|_{1}+\zeta\sum_{j=1}^{M_{n}+1}\mu_{j}+\lambda_{2}b^{T}Jb.$
(4)

Let $V=(v_{1}\cdots v_{L_{n}},v_{L_{n}+1}\cdots v_{L_{n}+p})^{T}$, and define
an augmented data set $(\overline{Y},\overline{V})$ (Zou & Hastie, 2005) as

$\overline{Y}=\begin{pmatrix}Y\\\
\mathbf{0}\end{pmatrix},\quad\quad\overline{V}=\begin{pmatrix}v_{1}\cdots
v_{L_{n}}&v_{L_{n}+1}\cdots v_{L_{n}+p}\\\
\sqrt{2n\lambda_{2}}D&\mathbf{0}^{T}\end{pmatrix},$

where $J=D^{T}D.$ It follows that, (4) can be written as

$\mathcal{Q}_{n}(\alpha)=\frac{1}{2n}(\overline{Y}-\overline{V}\alpha)^{T}(\overline{Y}-\overline{V}\alpha)+\sum_{j=1}^{M_{n}+1}\mu_{j}^{1-1/\gamma}\|b_{A_{j}}\|_{1}+\zeta\sum_{j=1}^{M_{n}+1}\mu_{j}.$
(5)

We propose the following algorithm:

  * Step 1.

For a given choice of $\lambda_{2}$, obtain an initial value of $\alpha$,
denoted by $\alpha^{(0)}$, by solving
$\operatorname*{arg\,min}\limits_{\alpha}\frac{1}{2n}(Y-V\alpha)^{T}(Y-V\alpha)+\lambda_{2}b^{T}Jb$
iteratively until convergence where an initial value is a $\mathbf{0}$ vector
(Gellar et al., 2015). The solution at each iteration is given by

$\hat{\alpha}=(V^{T}V+2n\lambda_{2}J^{*})^{-1}V^{T}Y$

where $J^{*}$ is a block diagonal matrix with $J$ and the zero matrix of size
$p\times p$ on the diagonal.

  * Step 2.

At the $m$th iteration, $m=1,2,\ldots$, compute

$\mu_{j}^{(m)}=\left(\frac{1-\gamma}{\gamma\zeta}\right)^{\gamma}\|b_{A_{j}}^{(m-1)}\|_{1}^{\gamma},\quad
j=1,\ldots,M_{n}+1,$

and

$\xi_{k}^{(m)}=\sum_{j:k\in A_{j}}(\mu_{j}^{(m)})^{1-1/\gamma},\quad
k=1,\ldots,L_{n}.$

  * Step 3.

Compute $\overline{Y}$ and $\overline{V}$ based on $\alpha^{(m-1)}$.

  * Step 4.

Obtain the new estimate by coordinate descent algorithm as

$\alpha^{(m)}=\operatorname*{arg\,min}_{\alpha}\left\\{\frac{1}{2n}(\overline{Y}-\overline{V}\alpha)^{T}(\overline{Y}-\overline{V}\alpha)+\sum_{k=1}^{L_{n}}\xi_{k}^{(m)}|b_{k}|\right\\}.$

  * Step 5.

Stop with convergence.

The algorithm always converges to a local minimizer of $\mathcal{Q}(\alpha)$
given fixed tuning parameters depending on the initial value of
$\alpha^{(0)}$, since the group penalty is not convex. To avoid the effect of
this issue, we suggest to use the multiple initial values and select the one
which can achieve the smallest value of $\mathcal{Q}(\alpha)$.

###  3.4 Assumptions and asymptotic results

In this section, we discuss the asymptotic properties of the estimators
defined in the previous section. We first introduce notations. Let
$(\beta_{0},\theta_{0})^{T}$ denote the true parameter value. We let
$\|\cdot\|_{2}$ denote the Euclidean norm and $\|\cdot\|_{\infty}$ denote the
supremum norm. In the Hilbert space $H=L^{2}(\mathcal{S})$, containing all
square integrable functions on $\mathcal{S}$, we define inner product $\langle
x,y\rangle=\int_{\mathcal{S}}x(s)y(s)ds$, $\forall x,y\in H$ and norm
$\|x\|=\langle x,x\rangle^{1/2}$. For two sequences of positive numbers
$a_{n}$ and $b_{n}$, $a_{n}\lesssim b_{n}$ indicates $a_{n}/b_{n}$ is
uniformly bounded and $a_{n}\asymp b_{n}$ if $a_{n}\lesssim b_{n}$ and
$b_{n}\lesssim a_{n}$. We define the covariance operator $\Gamma$ of the
random process $X$ as $\Gamma x(t)=\int_{0}^{1}EX(t)X(s)x(s)ds$, and
$\Gamma_{\Delta}x(t)=\int_{0}^{1}E\Delta X(t)\Delta
X(s)x(s)ds=\int_{0}^{1}E\Delta X(t)X(s)x(s)ds$. We denote $\Gamma_{n}$ as the
empirical version of $\Gamma$ of X; that is,
$\Gamma_{n}x(t)=n^{-1}\sum_{i=1}^{n}\int_{0}^{1}X_{i}(t)X_{i}(s)x(s)ds.$ Let
$H$ be the $n\times M_{n}$ matrix with element
$h_{i,j}=\langle\Gamma_{n}B_{i},B_{j}\rangle$. Let $\rho_{min}(A)$ and
$\rho_{max}(A)$ denote the minimum and maximum eigenvalues of a symmetric
matrix $A$, respectively. For a $s\times s$ matrix $A$,
$\|A\|_{\infty}=\max\\{\sum_{k=1}^{s}|(A)_{jk}|:1\leq j\leq s\\}$ and
$\|A\|_{2}=\\{\rho_{\max}(A^{T}A)\\}^{1/2}$. We define the seminorm
$\|f\|_{\Gamma}=\left\\{\int_{0}^{1}(\Gamma f)(s)f(s)ds\right\\}^{1/2}$. Let
$\mathcal{H}(r)$ be the collection of functions defined on $[0,1]$ whose $v$th
derivative exists and satisfies the Hölder condition of order $\delta$ with
$r\equiv v+\delta>2$, $|f^{(v)}(s_{1})-f^{(v)}(s_{2})|\leq
D|s_{1}-s_{2}|^{\delta}$, for $0\leq s_{1},s_{2}\leq 1$ and some constant
$D>0$.

We make the following assumptions:

  1. (A1)

$\beta\in\mathcal{H}(r)$.

  2. (A2)

Assume that the exposure function and the additional covariates have been
centered so that $E(\Delta Z)=0$ and $E(\Delta X(s))=0,s\in[0,1]$.

  3. (A3)

Assume $\tilde{T}\perp C|W$.

  4. (A4)

The observed event time $T_{i}$, $1\leq i\leq n$ is in a finite interval, say
$[0,\tau]$, and the baseline cumulative hazard function
$\int_{0}^{\tau}\lambda_{0}(s)ds<\infty$, and there exists a small positive
constant $\epsilon$ such that: (i) $P(\Delta=1|W)>\epsilon$, and (ii)
$P(C>\tau|W)>\epsilon$ almost surely with respect to the probability measure
of $W$.

  5. (A5)

The covariate $Z$ takes values in a bounded subset of $R^{p}$, and the
$L_{2}$-norm $\|X\|_{2}$ is bounded almost surely.

To have a sparse estimator of $\beta(s)$, we partition the set of number of
knots $\\{1,\ldots,M_{n}+1\\}$ into three sets (Wang & Kai, 2015)

$\displaystyle\mathcal{A}_{1}$
$\displaystyle=\\{j:\beta(s)=0,s\in[\kappa_{j-1},\kappa_{j}]\\},$
$\displaystyle\mathcal{A}_{2}$
$\displaystyle=\\{j:0<\max_{s\in[\kappa_{j-1},\kappa_{j}]}|\beta(s)|\leq
D^{*}M_{n}^{-r}\text{ for some }D^{*}>D\\},$ $\displaystyle\mathcal{A}_{3}$
$\displaystyle=\\{1,\ldots,M_{n}+1\\}-\mathcal{A}_{1}\cup\mathcal{A}_{2}.$

If (A1) holds, there is a vector of $\beta(s)^{*}=B^{T}(s)b^{*}$ where
$\left\|\beta_{0}(s)-\beta^{*}(s)\right\|_{\infty}\leq D_{0}M_{n}^{-r}$ for
some constant $D_{0}$. As suggested in Wang & Kai (2015), we introduce a
sparse modification of $\beta^{*}(s)\in\mathcal{S}_{dM_{n}}$ denoted by
$\beta_{0}^{*}(s)=B^{T}(s)b_{0n}$, where
$b_{0n}=(b_{0,1},\ldots,b_{0,L_{n}})^{T}$ is a sparse modification of $b^{*}$
with $b_{0,k}=b_{k}^{*}I(k\not\in\mathcal{B}_{1})$, and
$\mathcal{B}_{1}=\cup_{j\in\mathcal{A}_{1}\cup\mathcal{A}_{2}}A_{j}$. Note
that under (A1), $\|\beta_{0}-B^{T}b_{0n}\|_{\infty}=O(M_{n}^{-r})$ by Lemma 1
in Wang & Kai (2015).

  1. (A6)

$M_{n}=o(n^{1/4})$, $\lambda_{1}\nu_{n}=O(n^{-1/2}M_{n}^{1/2})$ with
$\nu_{n}=\left(\sum_{j\in\mathcal{A}_{3}}\|b_{0,A_{j}}\|_{1}^{2\gamma-2}\right)^{1/2}$
and $\lambda_{2}=o(n^{-1/2}M_{n}^{1/2})$.

The assumptions (A7) – (A10) are in the Supplementary Materials. The
assumptions (A1) – (A5), and (A7) – (A8) are the same as those in Qu et al.
(2016) and the assumptions (A6), (A9) and (A10) are similar to Huang et al.
(2014).

######  Theorem 1 (Convergence rate).

Define $D_{n}=M_{n}^{1/2}n^{-1/2}+M_{n}^{-r}$. Under conditions (A1) – (A10),
there exists a local maximizer $(\hat{b},\hat{\theta})$ such that

$\|\hat{\theta}-\theta_{0}\|_{2}=O_{p}(D_{n}),$
$\|\hat{\beta}-\beta_{0}\|=O_{p}(D_{n}),$

and

$\|\hat{\beta}-\beta_{0}\|_{\Gamma}=O_{p}(D_{n}).$

By Theorem 1, the optimal choice of $M_{n}\asymp O(n^{1/(2r+1)})$. Then,
$\|\hat{\theta}-\theta_{0}\|_{2}=O_{p}(n^{-r/(2r+1)})$ and
$\|\hat{\beta}-\beta_{0}\|_{\Gamma}=O_{p}(n^{-r/(2r+1)})$. The details of the
proof are provided in the Supplementary Materials, which show that there
exists a local maximum in the ball
$\\{\alpha_{0}+\delta_{n}u:\|u\|_{2}=c_{0}\\}$ with probability at least
$1-\epsilon$, where $\alpha_{0}=(b_{0n},\theta_{0})$.

######  Theorem 2 (Sparsistency).

Under conditions (A1) – (A10), we have
$(\hat{b}_{A_{j}}:j\in\mathcal{A}_{1}\cup\mathcal{A}_{2})=\bm{0}$ with
probability tending to 1.

This theorem guarantees that the non-null regions are identified by a set of
$\mathcal{A}_{1}\cup\mathcal{A}_{2}$, in which $\hat{\beta}(s)=0$ with large
probabilities when
$s\in\cup_{j:\mathcal{A}_{1}\cup\mathcal{A}_{2}}[\kappa_{j-1},\kappa_{j}]$.
The details of the proof are provided in the Supplementary Materials.

###  3.5 Choice of tuning parameters

The choice of tuning parameters determines the performance of the proposed
methods. We choose the number of interior knots $M_{n}$ to be relatively
large, following Marx & Eilers (1999); Lin et al. (2017). Alternatively, we
could use K-fold cross-validation to choose $M_{n}$ (Huang et al., 2004; Wang
& Kai, 2015). As suggested in Huang et al. (2010), we set $\gamma=0.5.$ Using
the techniques similar to those in Guan et al. (2020), we compute the
effective number of parameters by the Newton-Raphson method, given by

$df(\lambda_{1},\lambda_{2})=\text{tr}[\\{H_{0}+n\lambda_{2}J_{0}\\}^{-1}H_{0}],$

where $H_{0}$ denote the submatrix of $\bigtriangledown^{2}l(\hat{\alpha})$
with rows and columns corresponding to the nonzero $\hat{b}$, and $J_{0}$
denote the submatrix of $J$ with rows and columns corresponding to the nonzero
$\hat{b}$. We select tuning parameters $\lambda_{1}$ and $\lambda_{2}$ using a
BIC-type criterion,

$\text{BIC}(\lambda_{1},\lambda_{2})=-2\log
l_{n}(\hat{\alpha})+\log(n)df(\lambda_{1},\lambda_{2}).$ (6)

The BIC-type criterion is better suited to identify the true model structure
because it penalizes a model more for its complexity compared to AIC (Bishop &
Nasrabadi, 2006).

###  3.6 Variance estimation using two-stage approach

In this section, we introduce a two-stage approach to obtain the interval
estimates for $\beta(t)$ given the non-null regions. In the first stage, we
define the non-null region indicators as
$I_{0}=\bigcup_{j=1}^{M_{n}+1}\\{j:\hat{\beta}(s)\neq 0\ \text{for}\ \forall
s\in[\kappa_{j-1},\kappa_{j}]\\}$, and then non-null regions are
$\mathcal{S}_{s}=\bigcup_{j\in I_{0}}[\kappa_{j-1},\kappa_{j}]$, where
$\hat{\beta}(s)$ on $\mathcal{S}$ are estimated based on the algorithm in
Section 3.3. We define $A_{s}=\bigcup_{j\in J_{0}}A_{j}$. In the second stage,
given the identified non-null regions $\mathcal{S}_{s}$, we re-estimate
$\beta(s)$ in $s\in\mathcal{S}_{s}$, denoted by $\hat{\beta}^{*}(s)$, based on
the B-spline basis functions $\\{B_{j}(s):j\in A_{s}\\},$ by maximizing the
penalized partial likelihood only with the smoothness penalty, written as

$\mathcal{L}^{\ast}_{n}(\tilde{\alpha})=n^{-1}l_{n}^{*}(\tilde{\alpha})-\tilde{\lambda}_{2}\tilde{b}^{T}J_{0}\tilde{b},$
(7)

where the log partial likelihood is given by

$\displaystyle l_{n}^{\ast}(\tilde{\alpha})=$
$\displaystyle\sum_{i=1}^{n}\Delta_{i}\bigg{\\{}\sum\limits_{k\in
A_{s}}b_{k}\int_{0}^{1}X_{i}(s)B_{k}(s)ds+\theta^{T}Z_{i}$
$\displaystyle-\log\sum_{T_{j}\geq T_{i}}\exp\left(\sum\limits_{k\in
A_{s}}b_{k}\int_{0}^{1}X_{j}(s)B_{k}(s)ds+\theta^{T}Z_{j}\right)\bigg{\\}},$

$\tilde{b}=(b_{k},k\in A_{s})^{T}$, and
$\tilde{\alpha}=(\tilde{b}^{T},\theta^{T})^{T}$. Note that coefficients
corresponding to the B-spline basis functions are included only in the non-
null regions. Here, the smoothness penalty tuning parameter
$\tilde{\lambda}_{2}$ is chosen by a BIC-type criterion of the form (6).

Based on the Condition (C2) in Hao et al. (2021), we simultaneously
diagonalize
$H=-\partial^{2}l^{*}_{n}(\tilde{\alpha})/\partial\tilde{\alpha}\partial\tilde{\alpha}^{T}$
and

$\mathcal{P}=\begin{bmatrix}J_{0}&0_{q\times p}\\\ 0_{p\times q}&0_{p\times
p}&\end{bmatrix}$

by producing a matrix $R$ and a diagonal matrix $\mathcal{D}$ such that
$R^{T}HR=I$ and $R^{T}\mathcal{P}R=\mathcal{D}$. Now, we write $R$ and
$\mathcal{D}$ as

$R=\begin{bmatrix}R_{11}&R_{12}\\\
R_{21}&R_{22}\end{bmatrix},\quad\mathcal{D}=\begin{bmatrix}\Pi&0_{q\times
p}\\\ 0_{p\times q}&0_{p\times p}\end{bmatrix},$

where $\Pi$ is a diagonal matrix with diagonal elements
$\pi_{1},\ldots,\pi_{q}$ with an increasing order. Then, we define

$\psi_{l}(s)=\sum_{u=1}^{q}R_{12,ul}B_{u}(t),$
$\phi_{\nu}(s)=\sum_{u=1}^{q}R_{11,u\nu}B_{u}(t).$

By applying the Theorem 3.3 in Hao et al. (2021),

$\widehat{\operatorname{Var}}(\hat{\beta}^{*}(s))=\frac{1}{n}\left[\sum_{l=1}^{p}\psi_{l}^{2}(s)+\sum_{\nu=1}^{q}\frac{\varphi_{\nu}^{2}(s)}{\left(1+\tilde{\lambda}_{2}\pi_{\nu}\right)^{2}}\right],$
(8)

where $q=|A_{s}|$.

It is often of interest to estimate cumulative effects
$\int_{s\in\mathcal{S}_{s}}\hat{\beta}^{*}(s)ds$, which summarizes the overall
effect of exposure within a circular distance buffer. We can calculate the
variance of the cumulative effect as

$\widehat{\operatorname{Var}}\left(\int_{s\in\mathcal{S}_{s}}\hat{\beta}^{*}(s)ds\right)=\frac{1}{n}\left[\sum_{l=1}^{p}|\psi_{l}^{0}|^{2}+\sum_{\nu=1}^{q}\frac{|\varphi_{\nu}^{0}|^{2}}{\left(1+\tilde{\lambda}_{2}\pi_{\nu}\right)^{2}}\right],$

where $\psi_{l}^{0}=\int_{s\in\mathcal{S}_{s}}\psi_{l}(s)ds$, and
$\phi_{\nu}^{0}=\int_{s\in\mathcal{S}_{s}}\phi_{\nu}(s)ds$.

Note that this variance formula is only approximate because we do not take
into account the variance associated with the selection of the nonzero
$\beta(s)$ in the first stage. In our simulation studies, we observed that
this approximation performs well.

##  4 Simulation Studies

We conducted simulation studies to evaluate the finite sample performance of
the proposed method (denoted by Spline-Gbridge), compared with other methods
including the Spline method (Gellar et al., 2015), the Lasso method, the
Gbridge method (Huang et al., 2009) and the Spline-Lasso method (Guo et al.,
2016). The Spline method has the smoothness penalty and no sparseness penalty
in (2) by setting $\lambda_{1}=0$. The Lasso method has the lasso penalty
where the group bridge penalty is replaced with
$\lambda_{1}\sum_{j=1}^{L_{n}}|b_{j}|$ and no smoothness penalty in (2) by
setting $\lambda_{2}=0$. The Gbridge method has the group bridge penalty and
no smoothness penalty in (2) by setting $\lambda_{2}=0$. The Spline-lasso
method has the smoothness penalty and the lasso penalty in (2).

We generated failure time $\tilde{T}_{i}$ from an exponential distribution
with $\lambda_{0}(s)=1$ and a hazard function

$\lambda(t)=\lambda_{0}(t)\exp\left(\int_{0}^{1}X_{i}(s)\beta(s)ds+\theta_{1}Z_{i1}+\theta_{2}Z_{i2}\right),$

where the covariate function $X_{i}(s)$ was generated by
$X_{i}(s)=\sum_{j=1}^{52}c_{ij}B_{j}(s)$, where $B_{j}(s)$ are cubic B-spline
basis functions with 48 equally spaced inner knots, and $c_{ij}$ was generated
from the standard normal distribution. Two scalar covariates were generated by
$Z_{1}\sim N(1,0.5^{2})$ and $Z_{2}\sim N(0,1^{2})$, respectively. We set
$\theta_{1}=\log(0.8)$, and $\theta_{2}=\log(1.2).$ We consider three
scenarios for $\beta(s)$ with $s\in[0,1]$:

  * •

Scenario I: $\beta_{1}(s)=0$,

  * •

Scenario II: $\beta_{2}(s)=2\sin(2\pi s)I_{0\leq s<0.5}$,

  * •

Scenario III: $\beta_{3}(s)=-2\sin(\pi(s-0.5))I_{0\leq s<0.5}$,

where $I_{(\cdot)}$ denotes the indicator function. For the scenarios II and
III, $\beta(s)=0$ on $[0.5,1]$. The censoring time was generated from an
exponential distribution with a censoring rate of $10\%$.

In the fitted models, we set $M_{n}=26$ and $d=3$. We considered sample sizes
of $n=500$ and $1000$, and we independently generated $1000$ datasets for each
case. The performance of estimators $\hat{\beta}(s)$ was evaluated by the
integrated mean squared errors according to IMSE =
$\int_{0}^{1}\left(\hat{\beta}(s)-\beta(s)\right)^{2}ds/\int_{0}^{1}\beta(s)^{2}ds$
for Scenario II and III, and we consider the IMSE
$\int_{0}^{1}\left(\hat{\beta}(s)-\beta(s)\right)^{2}ds$ for Scenario I. We
also examined the performance of estimators $\hat{\theta}_{1}$ and
$\hat{\theta}_{2}$ with percent bias and empirical standard errors (in
parentheses). To assess the performance of $\hat{\beta}(s)$ on identifying the
non-null regions, equivalently the supremum of non-null regions in our
scenarios, we obtained the average of the supremum of non-null regions. The
tuning parameters $\lambda_{1}$ and $\lambda_{2}$ were selected based on BIC.

Table 1: Simulation results for three scenarios. The mean and standard
deviation (in parentheses) of IMSE for $\hat{\beta}(s)$ and a percent bias and
empirical standard error for $\hat{\theta}_{1}$ and $\hat{\theta}_{2}$. Tuning
parameters are selected based on BIC. We set a sample size of $n=1000$, and
the results are based on 1000 simulated data sets. The results for $n=500$ are
shown in the Supplementary Materials.

Scenario | Parameter | Spline | Lasso | Gbridge | Spline-Lasso | Spline-Gbridge
---|---|---|---|---|---|---
I | $\beta_{1}(s)$ | 0.154 (0.155) | 0.010 (0.056) | 0.084 (1.199) | 0.010 (0.062) | 0.046 (0.131)
$\theta_{1}$ | -0.145 (0.076) | -0.287 (0.076) | -0.325 (0.076) | -0.301 (0.076) | -0.293 (0.076)
$\theta_{2}$ | -0.539 (0.038) | -0.704 (0.038) | -0.675 (0.038) | -0.718 (0.038) | -0.642 (0.038)
II | $\beta_{2}(s)$ | 1.160 (0.483) | 0.907 (0.162) | 0.947 (1.265) | 0.527 (0.348) | 0.483 (0.344)
$\theta_{1}$ | -0.188 (0.076) | -0.840 (0.077) | -0.662 (0.078) | -0.385 (0.078) | -0.221 (0.077)
$\theta_{2}$ | -0.886 (0.038) | -1.585 (0.038) | -1.343 (0.038) | -1.074 (0.038) | -0.909 (0.038)
III | $\beta_{3}(s)$ | 0.237 (0.155) | 0.931 (0.134) | 0.971 (0.512) | 0.458 (0.378) | 0.416 (0.387)
$\theta_{1}$ | -0.027 (0.077) | -0.840 (0.078) | -0.714 (0.078) | -0.366 (0.077) | -0.265 (0.077)
$\theta_{2}$ | -0.757 (0.038) | -1.585 (0.038) | -1.426 (0.038) | -1.091 (0.038) | -1.011 (0.038)

Table 1 reports the mean and the standard deviation (in parentheses) of IMSE
when the sample size is $n=1000$. The estimators based on the Spline-Gbridge
model outperform over those based on the Gbridge model. For Scenario I where
no effect exists, the Lasso and Spline-Lasso models give smaller IMSE than the
other models. However, for Scenario II and III, the Spline-Gbridge estimators
have the smallest IMSE than the other methods. The simulation results for
$n=500$ are in the Supplementary Material. Overall, our proposed model Spline-
Gbridge gives better performance in terms of IMSE.

Table 2 reports the mean and standard deviation (in parentheses) of the
supremum of the non-null region, i.e., the buffer distance. In the true model,
the buffer distance is 0 for Scenario I and 0.5 for Scenario II and III. The
Spline model performs well in terms of IMSE but it cannot identify the buffer
distance. The performance of identifying the buffer distance based on the
Lasso and Spline-Lasso for Scenario II and III is poor. The simulation results
for $n=500$ are in the Supplementary Materials. The Spline-Gbridge estimators
outperform the other methods in terms of estimation accuracy as well as
identifying the buffer distance by exploiting the local support property of
B-spline basis functions (Wang & Kai, 2015).

Table 2: Simulation results for three scenarios. The average of the supremum of the non-null regions. Tuning parameters are selected based on BIC. We set a sample size of $n=1000$, and the results are based on 1000 simulated data sets. The results for $n=500$ are shown in the Section 3 in the Supplementary Materials. Scenario | Lasso | Gbridge | Spline-Lasso | Spline-Gbridge
---|---|---|---|---
I | 0.058 (0.199) | 0.040 (0.166) | 0.041 (0.177) | 0.109 (0.260)
II | 0.136 (0.208) | 0.208 (0.223) | 0.683 (0.429) | 0.563 (0.298)
III | 0.091 (0.172) | 0.124 (0.192) | 0.683 (0.445) | 0.597 (0.342)

We conducted additional simulation studies to assess the accuracy of the
estimate introduced in Section 3.6 of the standard error of
$\hat{\beta}^{*}(s)$. Therefore, we conducted the simulation studies with
Scenarios II and III given the true non-zero areas $[0,0.5]$ as follows. Given
the true non-zero areas $[0,0.5]$, we estimate $\beta(s)$ in $[0,0.5]$,
denoted by $\hat{\beta}^{*}(s)$, by maximizing (7) and estimate the standard
error of $\hat{\beta}^{*}(s)$ using the formula (8). The smoothness penalty
tuning parameter was chosen by the BIC-type criterion.

Figure 2: Simulation results for two-stage variance estimation for Scenarios
II and III on the true non-zero areas $[0,0.5]$. The black solid lines
represent the average of the estimated standard error (ASE), and the blue
dotted lines represent the empirical standard error (ESE) for Scenario II (the
left penal) and Scenario III (the right penal). Figure 3: Simulation results
for two-stage variance estimation for scenarios II and III on the true non-
zero areas $[0,0.5]$. The left penal shows the coverage probability for
$\beta(s)$ for Scenario II and the right penal shows the the coverage
probability for $\beta(s)$ for Scenario III.

Figure 2 displays the average of the estimated standard error (ASE) of
$\beta(s)$ and the empirical standard error (ESE) of $\beta(s)$. The average
of the estimated standard errors and the empirical standard erros are in good
agreement in both Scenarios. Figure 3 illustrates the coverage probability of
$\beta(s)$. Scenario II shows a slight undercoverage at the end and the center
of the interval, however, overall, it shows good performance. For Scenario
III, we see a generally good performance over the interval.

##  5 Case Study

The aim of the case study is to investigate the radial buffer distance of NDVI
after which there is no effect on depression incidence and to estimate the
effects of green space on depression incidence in NHS within the selected
buffer distance. The objective measurements of green space consist of the NDVI
data assessed within circular distance buffers at nine radii ($s_{1}$=90m,
$s_{2}$=150m, $s_{3}$=270m, $s_{4}$=510m, $s_{5}$=750m, $s_{6}$=990m,
$s_{7}$=1230m, $s_{8}$=1500m, and $s_{9}$=2100m). We obtained the average NDVI
values in the distance buffers between radii $s_{j-1}$ and $s_{j}$, denoted by
$X_{i}(s_{j})$, for $j=1,\ldots,9$ and we let $s_{0}=0$. Figure 4 shows the
pattern of the average NDVI values in concentric rings based on circular
distance buffers at radii of 0m, 90m, 150m, 270m, 510m, 750m, 990m, 1230m,
1500m, and 2100m for a total of 35,283 eligible participants in NHS. The range
of the average NDVI was $(0.48,0.49)$ and the maximum NDVI of 0.49 was
observed at radius 270m and the minimum NDVI of 0.48 was observed at radius
90m. Since the NDVI values were obtained at few radii, we use linear
interpolation to estimate individual NDVI values as continuous NDVI functions
$X_{i}(s)$.

Figure 4: The average NDVI values in concentric rings based on circular
distance buffers at radii of 0m, 90m, 150m, 270m, 510m, 750m, 990m, 1230m,
1500m, and 2100m for a total of 35,283 eligible participants in NHS.

In a naive analysis, a Cox regression model was fitted for each circular
distance buffer separately. To apply our proposed method, we used the two-
stage approach described in Section 3.6, where in the first stage, we fit the
functional linear Cox model with the sparseness and smoothness penalties to
identify non-null regions and find the buffer distance by the supremum of a
set of non-null regions, and in the second stage, we re-fit the functional
linear Cox model with the smoothness penalty within the buffer distance
obtained from the first stage. In both the naive analysis and the analysis
using our method, we stratified by age in months and additionally adjusted for
potential risk factors including race, smoking status, pack-years of smoking,
alcohol consumption, marital status, educational attainment, husband’s
educational attainment, baseline Mental Health Inventory-5 score, and area-
level characteristics such as Census tract population density, Census tract
median income, and Census tract median home value. We set the B-spline basis
$\\{B_{j}(s):j=1,\ldots,9\\}$ with degree 3, $M_{n}=5$, and inner knots at
150m, 270m, 510m, 990m, 1500m, and $\mathcal{S}$ as $[90,2100]$m. The tuning
parameters were selected based on BIC.

Radius(m) | Hazard ratio | 95% CI | p-value
---|---|---|---
$90$ | 0.952 | (0.928, 0.976) | 0.0001
$150$ | 0.950 | (0.926, 0.975) | $<0.0001$
$270$ | 0.951 | (0.926, 0.976) | 0.0001
$510$ | 0.955 | (0.930, 0.980) | 0.0005
$750$ | 0.955 | (0.930, 0.980) | 0.0005
$990$ | 0.956 | (0.932, 0.982) | 0.0008
$1230$ | 0.958 | (0.933, 0.983) | 0.0011
$1500$ | 0.958 | (0.934, 0.983) | 0.0011
$2100$ | 0.959 | (0.935, 0.983) | 0.0011
Table 3: The estimated hazard ratios, confidence intervals, and p-value for
0.1 unit increase in the average NDVI within each of nine circular distance
buffers (90m, 150m, 270m, 510m, 750m, 990m, 1230m, 1500m, and 2100m) based on
separate stratified Cox regression models.

Table 3 shows the results of naive analysis including the hazard ratios (HRs),
confidence intervals, and p-value of 0.1 unit increase in the average NDVI
within nine circular distance buffers, respectively, using the stratified Cox
regression model. The NDVI decreased the risk of depression in all nine
circular distance buffers, and the risk of depression was smallest at 150m
circular distance buffer (HR: 0.950, 95% CI: 0.926–0.975, p-value$<0.001$) and
the risk of depression increased as the buffer distance increased. In the
naive analysis, we were not able to estimate the effect of NDVI as a function
of radius; therefore, it was not possible to identify the buffer distance.

In our proposed model, the non-null regions of $[0,510]$m were identified;
thus, the buffer distance of 510m was chosen. Given the buffer distance of
510m, $\hat{\beta}^{*}(s)$ and its variance were estimated using the
functional Cox model with the smoothness penalty. Figure 5 shows the hazard
ratio of 0.1 unit increase in the average NDVI at a radius in meters of $s$
given the 510m buffer distance, i.e., $\exp\left(0.1\hat{\beta}^{*}(s)\right)$
while controlling for the NDVI values at different radii and confounders. The
estimated cumulative effect for 0.1 unit increase in the average NDVI within
the 510m buffer distance
$\left(\exp\left(\int_{0}^{510}0.1\hat{\beta}^{*}(s)\right)\right)$ was 0.946
(95% CI: 0.920–0.973, p-value=0.0001), which was comparable to the results
based on the naive analysis using the 510m circular distance buffer (HR:
0.955, 95% CI: 0.930–0.980, p-value=0.0005).

Figure 5: The solid black line is the estimated hazard ratio of 0.1 unit
increase in the average NDVI and the shaded area (gray) is its 95% confidence
interval as a function of radius within the buffer distance of 510m using the
proposed two-stage method.

##  6 Discussion

In our greenness and depression study, we investigated the risk of depression
associated with the greenness of the neighborhood environment and identified
the data-informed optimal buffer distance used to measure green space
exposure. To overcome the limitations of traditional analysis, which uses
radial buffer distances selected a priori, we propose the two-stage functional
Cox regression model, in which the optimal buffer distance is identified in
the first stage, and given this optimal buffer distance, we estimate the
functional association between incident depression and residential greenness
exposure. Our analysis result using this proposed method suggests that the
greenness within a 510m buffer distance from the residence had a protective
effect on the risk of developing depression, and the protective effect was
strongest at 90m and decreased as the distance from the residence increased.

In our study based on the NHS, NDVI exposure was sparsely measured at only
nine different circular distance buffers and we used the linear interpolated
NDVI exposures to obtain regular and dense functional data. The identification
of buffer distance may be sensitive to the choice of number and locations of
knots in the B-spine basis functions; therefore, dense NDVI exposure may be
required for future research. In addition, the statistical inference in the
second stage does not take into account the variability induced by the
identification of non-null regions. The development of appropriate statistical
models to account for post-selection inference in functional linear models is
a topic for our future research. Furthermore, our proposed method is based on
finite basis functions such as B-spline functions, and we applied the
theoretical results in the functional linear Cox models under the reproducing
kernel Hilbert space (RKHS) (Hao et al., 2021) to the second stage in our
proposed two-stage approach. It may be of interest to apply the RKHS framework
to the selection problem in functional coefficients for future research.

## Supplementary Materials

The Supplementary Materials contain detailed proofs of all theorems,
additional simulation results, and R code for the simulation studies and the
data example.

## Acknowledgements

This project was supported by the National Research Foundation of Korea grant
funded by the Korean government (MSIT) (No. 2022R1C1C1011806), the NIH/NIEHS
grant R01ES026246, P30 ES00000, UM1 CA186107, the NIH/NHLBI grant R01
HL150119, and the DoD grant W81XWH2210030.

## References

  * (1)
  * Baek et al. (2016) Baek, J., Sánchez, B. N., Berrocal, V. J. & Sanchez-Vaznaugh, E. V. (2016), ‘Distributed lag models: examining associations between the built environment and health’, Epidemiology (Cambridge, Mass.) 27(1), 116.

  * Banay et al. (2019) Banay, R. F., James, P., Hart, J. E., Kubzansky, L. D., Spiegelman, D., Okereke, O. I., Spengler, J. D. & Laden, F. (2019), ‘Greenness and depression incidence among older women’, Environmental Health Perspectives 127(2), 027001.
  * Bishop & Nasrabadi (2006) Bishop, C. M. & Nasrabadi, N. M. (2006), Pattern recognition and machine learning, Vol. 4, Springer.
  * Cardot et al. (2003) Cardot, H., Ferraty, F. & Sarda, P. (2003), ‘Spline estimators for the functional linear model’, Statistica Sinica pp. 571–591.
  * De Boor (2001) De Boor, C. (2001), A practical guide to splines, Springer-Verlag New York.
  * Gellar et al. (2015) Gellar, J. E., Colantuoni, E., Needham, D. M. & Crainiceanu, C. M. (2015), ‘Cox regression models with functional covariates for survival data’, Statistical Modelling 15(3), 256–278.
  * Guan et al. (2020) Guan, T., Lin, Z. & Cao, J. (2020), ‘Estimating truncated functional linear models with a nested group bridge approach’, Journal of Computational and Graphical Statistics 29(3), 620–628.
  * Guo et al. (2016) Guo, J., Hu, J., Jing, B.-Y. & Zhang, Z. (2016), ‘Spline-lasso in high-dimensional linear regression’, Journal of the American Statistical Association 111(513), 288–297.
  * Hall & Hooker (2016) Hall, P. & Hooker, G. (2016), ‘Truncated linear models for functional data’, Journal of the Royal Statistical Society: Series B (Statistical Methodology) 78(3), 637–653.
  * Hao et al. (2021) Hao, M., Liu, K.-y., Xu, W. & Zhao, X. (2021), ‘Semiparametric inference for the functional cox model’, Journal of the American Statistical Association 116(535), 1319–1329.
  * Huang et al. (2010) Huang, J., Horowitz, J. L. & Wei, F. (2010), ‘Variable selection in nonparametric additive models’, Annals of Statistics 38(4), 2282.
  * Huang et al. (2014) Huang, J., Liu, L., Liu, Y. & Zhao, X. (2014), ‘Group selection in the cox model with a diverging number of covariates’, Statistica Sinica pp. 1787–1810.
  * Huang et al. (2009) Huang, J., Ma, S., Xie, H. & Zhang, C.-H. (2009), ‘A group bridge approach for variable selection’, Biometrika 96(2), 339–355.
  * Huang et al. (2004) Huang, J. Z., Wu, C. O. & Zhou, L. (2004), ‘Polynomial spline estimation and inference for varying coefficient models with longitudinal data’, Statistica Sinica pp. 763–788.
  * James et al. (2009) James, G. M., Wang, J., Zhu, J. et al. (2009), ‘Functional linear regression that’s interpretable’, The Annals of Statistics 37(5A), 2083–2108.
  * James et al. (2014) James, P., Berrigan, D., Hart, J. E., Hipp, J. A., Hoehner, C. M., Kerr, J., Major, J. M., Oka, M. & Laden, F. (2014), ‘Effects of buffer size and shape on associations between the built environment and energy balance’, Health & Place 27, 162–170.
  * James et al. (2018) James, S. L., Abate, D., Abate, K. H., Abay, S. M., Abbafati, C., Abbasi, N., Abbastabar, H., Abd-Allah, F., Abdela, J., Abdelalim, A. et al. (2018), ‘Global, regional, and national incidence, prevalence, and years lived with disability for 354 diseases and injuries for 195 countries and territories, 1990–2017: a systematic analysis for the global burden of disease study 2017’, The Lancet 392(10159), 1789–1858.
  * Kong et al. (2018) Kong, D., Ibrahim, J. G., Lee, E. & Zhu, H. (2018), ‘Flcrm: Functional linear cox regression model’, Biometrics 74(1), 109–117.
  * Kwan (2012) Kwan, M.-P. (2012), ‘The uncertain geographic context problem’, Annals of the Association of American Geographers 102(5), 958–968.
  * Lin et al. (2017) Lin, Z., Cao, J., Wang, L. & Wang, H. (2017), ‘Locally sparse estimator for functional linear regression models’, Journal of Computational and Graphical Statistics 26(2), 306–318.
  * Liu et al. (2019) Liu, Y., Wang, R., Grekousis, G., Liu, Y., Yuan, Y. & Li, Z. (2019), ‘Neighbourhood greenness and mental wellbeing in guangzhou, china: What are the pathways?’, Landscape and Urban Planning 190, 103602.
  * Marx & Eilers (1999) Marx, B. D. & Eilers, P. H. (1999), ‘Generalized linear regression on sampled signals and curves: a p-spline approach’, Technometrics 41(1), 1–13.
  * Pun et al. (2018) Pun, V. C., Manjourides, J. & Suh, H. H. (2018), ‘Association of neighborhood greenness with self-perceived stress, depression and anxiety symptoms in older us adults’, Environmental Health 17, 1–11.
  * Qu et al. (2016) Qu, S., Wang, J.-L., Wang, X. et al. (2016), ‘Optimal estimation for the functional cox model’, The Annals of Statistics 44(4), 1708–1738.
  * Ramsay & Silverman (1997) Ramsay, J. & Silverman, B. (1997), Functional data analysis, Springer.
  * Rhew et al. (2011) Rhew, I. C., Vander Stoep, A., Kearney, A., Smith, N. L. & Dunbar, M. D. (2011), ‘Validation of the normalized difference vegetation index as a measure of neighborhood greenness’, Annals of Epidemiology 21(12), 946–952.
  * Sarkar et al. (2018) Sarkar, C., Webster, C. & Gallacher, J. (2018), ‘Residential greenness and prevalence of major depressive disorders: a cross-sectional, observational, associational study of 94879 adult uk biobank participants’, The Lancet Planetary Health 2(4), e162–e173.
  * Spiegelman (2010) Spiegelman, D. (2010), ‘Approaches to uncertainty in exposure assessment in environmental epidemiology’, Annual Review of Public Health 31, 149–163.
  * van den Bosch & Meyer-Lindenberg (2019) van den Bosch, M. & Meyer-Lindenberg, A. (2019), ‘Environmental exposures and depression: biological mechanisms and epidemiological evidence’, Annual Review of Public Health 40, 239–259.
  * Wang & Kai (2015) Wang, H. & Kai, B. (2015), ‘Functional sparsity: Global versus local’, Statistica Sinica pp. 1337–1354.
  * World Health Organization (2017) World Health Organization (2017), Depression and other common mental disorders: global health estimates, Technical report, World Health Organization.
  * Yuan & Cai (2010) Yuan, M. & Cai, T. T. (2010), ‘A reproducing kernel hilbert space approach to functional linear regression’, The Annals of Statistics 38(6), 3412–3444.
  * Zhang & Lu (2007) Zhang, H. H. & Lu, W. (2007), ‘Adaptive lasso for cox’s proportional hazards model’, Biometrika 94(3), 691–703.
  * Zou & Hastie (2005) Zou, H. & Hastie, T. (2005), ‘Regularization and variable selection via the elastic net’, Journal of the Royal Statistical Society: Series B (statistical methodology) 67(2), 301–320.


# Protecting IoT Servers Against Flood Attacks with the Quasi Deterministic
Transmission Policy

Erol Gelenbe  Institute of Theoretical & Applied Informatics
Polish Academy of Science, 44-100 Gliwice, PL
& CNRS I3S Lab., Université Côte d’Azur, FR
& Yaşar University, Bornova, Izmir, TR
ORCID:0000-0001-9688-2201  Mohammed Nasereddin  Institute of Theoretical &
Applied Informatics
Polish Academy of Science, 44-100 Gliwice, PL
ORCID: 0000-0002-3740-9518

###### Abstract

IoT Servers that receive and process packets from IoT devices should meet the
QoS needs of incoming packets, and support Attack Detection software that
analyzes the incoming traffic to identify and discard packets that may be part
of a Cyberattack. Since UDP Flood Attacks can overwhelm IoT Servers by
creating congestion that paralyzes their operation and limits their ability to
conduct timely Attack Detection, this paper proposes and evaluates a simple
architecture to protect a Server that is connected to a Local Area Network,
using a Quasi-Deterministic Transmission Policy Forwarder (SQF) at its input
port. This Forwarder shapes the incoming traffic, sends it to the Server in a
manner which does not modify the overall delay of the packets, and avoids
congestion inside the Server. The relevant theoretical background is briefly
reviewed, and measurements during a UDP Flood Attack are provided to compare
the Server performance, with and without the Forwarder. It is seen that during
a UDP Flood Attack, the Forwarder protects the Server from congestion allowing
it to effectively identify Attack Packets. On the other hand, the resulting
Forwarder congestion can also be eliminated at the Forwarder with “drop”
commands generated by the Forwarder itself, or sent by the Server to the
Forwarder.

###### Index Terms:

Internet of Things (IoT), Cyberattack Detection, Congestion, Traffic Shaping,
Quasi-Deterministic Transmission Policy (QDTP), Quality of Service

##  I Introduction

With some $30$ Billion devices on the Internet [1], many types of anomalies
have been observed as a result of cyberattacks [2, 3, 4, 5], including Denial
of service (DoS) attacks that disable target systems by flooding them with
huge streams of requests [6]. While many such attacks go unreported when they
occur, just one Distributed DoS attack in $2017$ targeting Google, compromised
$180,000$ web servers which flooded Google servers at overall bitrates of
$2.54$ Tera-bits/sec [7]. Other attacks aim mainly at the IoT [8, 9, 10, 11],
while Botnet attacks [12] are particularly vicious since they spread by
inducing their victims to become attackers [13, 14, 15].

UDP Flood attacks [16] are also exploited by Botnets to create massive
congestion that overwhelms network nodes and ports. Using spoofed-source-
address UDP packets, they cause their victims to crash due to high traffic
volumes, creating denial of service, causing lost data and resulting in
missing and incomplete readings of the data carried by legitimate IoT traffic
[17, 18].

###  I-A Prior Work

Because of the concern about cybersecurity, there is a large literature on
cyberattacks and Attack Detection (AD) methods [19, 20, 6]. These methods are
typically evaluated for accuracy using statistical methods [11, 21], and
various Machine Learning based AD algorithms are often tested under ideal
conditions on general purpose computers [22, 23, 24, 25, 26], where attack
traffic is treated as data, and the attack’s overload on the processing
capacity and performance of the victim node is not taken into account.

Various AD test-beds [27] for cyber-physical and IoT networks are presented in
[28, 29, 30]. Experiments on windfarms under SYN attacks are discussed in [31]
and other experimental IoT security studies can be found in [32, 33, 34]. Data
collection and display for flood attacks are discussed in [35] while in [36]
real-time data collection for IoT DNS attacks is presented. Denial of Service
(DoS) attacks against Software Defined Networks (SDN) that support the IoT
have also been studied [37]. However this prior work relates to attack
emulation environments that do not include the overload caused by attacks, as
recently discussed for autonomous vehicles[38] and IoT servers [39].

###  I-B Motivation and Research Plan

Thus the present paper is motivated by the need to:

  * •

Experimentally evaluate the effect of IoT Server overload during an ongoing
UDP Flood attack,

  * •

Understand the attack’s impact on the Server’s capacity to carry out Attack
Detection (AD) and other useful processing functions,

  * •

Demonstrate a system architecture, and a traffic shaping policy [40] that was
initially proposed to mitigate the IoT’s Massive Access Problem (MAP) [41, 42,
43], to guarantee that in the presence of attacks that create large packet
flows, the Server can operate seamlessly and accomplish its role for AD and
other useful IoT processing functions,

  * •

Experimentally demonstrate that mitigation actions can be triggered to rapidly
eliminate the long-term effects of such UDP Flood attacks from the system as a
whole.

Thus in Section II, we provide new measurements on the experimental test-bed
shown in Figure 1, to illustrate the effect of a UDP Flood attack emanating
from an IoT traffic source against the IoT Server that receives packets from
different IoT devices. These measurements show that the Server is
significantly impacted during an attack and is impeded from conducting its AD
functions in a timely fashion.

Based on this observation, Section III proposes and evaluates a novel system
architecture shown in Figure 6 where the Server is preceded by a Smart “Quasi-
Deterministic Forwarding Policy (QDP)” Forwarder (SQF) that shapes the traffic
that is forwarded to the Server. Our results show that if we select the SQF
parameters based on mathematical principles [40], then the SQF effectively
limits the undesirable effects of an attack against the Server. However,
attack packets accumulate at the SQF which protects the Server, and mitigation
actions may discard the accumulated attack traffic.

##  II Initial Measurement Results

As a first step, we have conducted experiments on the Local Area Network (LAN)
test-bed for a system shown in Figure 1, in which IoT devices represented by
several Raspberry Pi machines, send UDP traffic to the Server. One of the Pi
machines is also programmed to generate attack traffic either at predetermined
instants or at random.

These Raspberry Pi $4$ Model B Rev $1.2$ machines (RPi$1$ and RPi$2$)
machines, each have a $1.5$GHz ARM Cortex-A$72$ quad-core processor and $2$GB
LPDDR$4-3200$ SDRAM, running Raspbian GNU/Linux $11$ (bullseye), a Debian-
based operating system optimized for Raspberry Pi hardware. The normally
operating (uncompromised) Raspberry PIs periodically send UDP Protocol packets
that simply contain the measurements of the temperature of the Raspberry Pi,
to the Server that is shown in Figure 1. The choice of this particular data is
simply in order to provide an example of real data that the Raspberry Pis can
send.

The Server is an Intel 8-Core i$7-8705$G, with a Linux $5.15.0-60-$ generic
$66-$Ubuntu SMP based operating system. It is equipped with $16$GB of RAM, it
runs at $3.1$Ghz and has a $500$GB hard drive. It communicates with the
Raspberry Pis via the Ethernet Local Area Network (LAN) shown in Figure 1, and
receives IoT traffic from them via the UDP protocol.

As shown in Figure 2, the Server supports the UDP protocol with SNMP for
incoming packets. It operates the accurate AD algorithm reported in [26], and
supports the other normal processing needs of incoming UDP packets. The UDP
protocol’s simplicity fits the needs of the simple IoT devices that we use,
since UDP does not establish a connection before transmitting and does not use
ACKs or error recovery for communications [44].

While many datasets can be used to generate attack traffic, including the
KDD$99$ dataset or its improved version NSL-KDD, UNSW-NB15, CICDS$2017$, and
the Bot-IoT dataset [36], in this work we use MHDDoS [45] containing $56$
recent real-world DoS attacks with $56$ different techniques, for attack
traffic emulation in this Ethernet based test-bed.

Figure 1: The test-bed composed of an Ethernet network with Raspberry Pi
machines that generate normal traffic, as well as possible traffic. An IoT
Server receives the IoT traffic via the same network. Figure 2: Internal
software architecture of the Server, containing a SNMP network protocol
manager, an AD system [26] that identifies attack packets, and software for
processing the incoming packet contents.

In Figure 3 we show measurements of the effect of a $60$ sec Flood Attack,
which overwhelms the 8-Core Server with some $400,000$ packets that accumulate
at its input buffer. The Server’s activities, including AD, are paralyzed by
the attack and the packet backlog takes nearly $300$ minutes to clear. Thus we
see that the attack significantly impairs its capability for AD, its ability
to discard attack packets and to subsequently process benign packets.

Figure 3: Experimentally measured queue length (the $y$-axis is in number of
packets) over time (the $x$-axis is in seconds) at the Server input prior to
processing at the AD module, during a $60$ second UDP Flood Attack launched
from one of the Rasberry Pis of Figure 1 against the Server. The backlog of
packets at the Server initially rises rapidly to $400,000$ packets, and
without human intervention the congestion at the Server lasts far longer than
the attack itself, up to several hours, due to the fact that the Server is
paralyzed and stops its AD processing packets for long time intervals. These
long interruptions in AD processing time are observed as the large outliers in
AD processing times in Figure 5.

The detailed measurements of the Server’s AD processing times per packet, when
there is no attack, and when a UDP Flood Attack occurs, are reported in Figure
4 and Figure 5.

We observe that the Server’s AD processing time per packet, when no attacks
occur, has an average value of $2.98$ milliseconds (ms). On the other hand,
when the Server is targeted by a UDP Flood Attack, we observe a substantial
increase in the AD algorithm’s average processing time to $4.82$ ms. Moreover,
the AD processing time per packet when the Server is under attack, exhibits
some very large outliers, as shown in Figure 5. We observe that these
“outlier” processing times are close to $10^{3}$ times larger than the typical
values, showing that during a UDP Flood Attack the Server’s AD procesing of
packets is repeatedly paralyzed for a substantial amount of time, as also
shown in Figure 3.

Figure 4: In the Upper figure, we show the histogram of measurements of the
Server’s AD processing time per packet, when there is no attack, exhibiting an
average processing time of $2.98$ ms and variance $0.0055~{}ms^{2}$. In the
Lower figure an attack is occurring: the Server’s measured average AD
processing time of packets rises substantially to $4.82$ ms with a variance
$0.51~{}ms^{2}$,

Figure 5: In the Upper figure, we show successive measurements of the Server’s
AD processing time per packet during a UDP Flood Attack (in the absence of the
QDTP Forwarder SQF), showing large outliers that initially become more severe,
and gradually become less frequent over time. In the Lower figure, the AD
processing time of packets that is measured after the UDP Attack begins,
reveals very large outliers in AD processing times, indicating that the AD is
intermittently paralyzed or unable to operate.  Figure 6: The figure shows the
modified system architecture where a Smart QDTP Forwarder (SQF) constantly
acts as a traffic shaping interface between the Ethernet LAN and the Server.
The effect of the SQF is to eliminate the paralyzing effect of the packet
flood at the Server, buffering packets within the SQF and forwarding in a
manner which allows the Server to conduct its AD processing and other work in
a timely fashion.

Figure 7: The queue length at the Server in the presence of a $60$ second UDP
Flood Attack. The figure Above shows the case without the SQF, and we see that
the queue length peaks to $400,000$ packets and descends slowly over some
$15,000$ seconds. The figure Below compares the queue length in logarithmic
scale, with SQF in Blue using the parameter $D=3~{}ms$, against the case
without SQF in Red, with the same UDP Flood Attack which lasts $60$ seconds.
Note that because the value of $D$ we use is very close to the average value
of $T_{n}$ measured to be $2.98$ ms in the absence of an attack, as shown in
Figure , the fluctuations in the values of $T_{n}$ will cause a small queue
buildup (in the order of a few packets), as seen in te Blue plot in the figure
Below.

Figure 8: We measure the Server queue length (repesented logarithmically) when
the Server is targeted by a UDP Flood Attack that lasts $30-sec$ (Above) and
$10-sec$ (Below). The Red curves show the case without the SQF traffic shaping
Forwarder, while the Blue Curves show the effect of the use of the SQF which
uses QDTP. We observe the huge difference in queue length. For both the $30$
and $10$ second attacks, we have set $D=3~{}ms$.

Figure 9: AD Processing Time at the Server when the SQF with the QDTP Policy
is installed and the parameter $D=2.7~{}ms$ is used. We observe that the AD
processing time $T_{n}$ has an average value of $2.97~{}ms$ and variance of
$0.0041~{}sec^{2}$ in the absence of an attack (Above). In the presence of a
UDP Flood Attack (Below) the average processing time of the AD per packet is
higher by roughly 10% on average, at $3.28~{}ms$ with a variance of
$0.0023~{}sec^{2}$ so that the SQF is effective in protecting the Server from
paralysis and excessive slowdown.

Figure 10: AD Processing Time at the Server when the SQF with the QDTP Policy
is installed and the parameter $D=3.20~{}ms$ is used. We observe that the AD
processing time $T_{n}$ has an average value of $3.00~{}ms$ and variance of
$0.0036~{}sec^{2}$ in the absence of an attack (Above). In the presence of a
UDP Flood Attack (Below) the average processing time of the AD per packet is
quasi-identical on average, at $2.99~{}ms$ with a variance of
$0.0067~{}sec^{2}$ so that in this case too, the SQF is effective in
protecting the Server from paralysis and excessive slowdowns. Figure 11: SQF
queue length ($y$-axis in number of packets) against time ($x$-axis in
seconds) when a UDP Flood attack lasts for $60$ seconds. We have used
$D=3~{}ms$, and no mitigation action takes place.

###  II-A Lindley’s Equation when the SQF is not Used

If the SQF module is not being used as shown in Figure 1:

  * •

Let $0=a_{0}\leq a_{1}\leq a_{2},~{}...~{}$, be the successive packet arrival
instants at the Server through the Ethernet LAN from any of the IoT devices
connected to the LAN. We also define the interarrival time
$A_{n+1}=a_{n+1}-a_{n}$.

  * •

Let $T_{n}$ denote the Server’s AD processing time for the $n-th$ packet, and
assume that the Server processes packets in First Come First Served (FCFS)
order.

Then the total waiting time $L_{n+1}$ of the $n+1$-th incoming packet, between
the instant $a_{n}$ and the start of the AD processing time of the Server, is
given by the well known Lindley’s equation [46, 47]:

$L_{n+1}=[L_{n}+T_{n}-A_{n+1}]^{+},~{}n\geq 0,~{}L_{0}=0~{},$ (1)

where for a real number $X$, we use the notation:

$[X]^{+}=X~{}if~{}X>0,~{}and~{}[X]^{+}=0~{}if~{}X\leq 0.$ (2)

Note that $L_{0}=0$ because the first incoming packet encounters an empty
queue infront of the AD. Note also that whenever we have $T_{n}>A_{n+1}$ then
$L_{n+1}>L_{n}$, i.e. the waiting time increases.

During a Flood Attack, the values of $A_{n}$ and $T_{n}$ will be modified, as
we see from Figure 3, indicating that packet arrival rates have considerably
increased so that the values of $A_{n}$ are much smaller, while Figure 4 shows
that the values of $T_{n}$ are also larger. However the form of (1) does not
change.

##  III Effect of the Smart QDTP Forwarder (SQF)

In Figure 6, we present our proposed modified architecture where the Server,
whose role is to process incoming IoT packets – including operating the AD
module in order to detect attacks – is “protected” by a Smart QDTP Forwarder
(SQF) which is placed between the Ethernet based sources of IoT traffic, and
the Server’s input port. The role of the SDF is to shape the incoming traffic
directed at the Server using the Quasi-Deterministic Transmission Policy
(QDTP) [48, 40].

QDTP is a simple policy that delays some of the packets it receives, by
forwarding them to the Server at time $t_{n}\geq a_{n}$, where $a_{n}$ is the
$n$-th packet’s arrival instant to the SQF, and $t_{n}$ is the instant at
which SQF forwards the packet to the Server, and is defined by:

$\displaystyle t_{n+1}=\max\\{t_{n}+D,a_{n+1}\\},~{}t_{0}=a_{0},~{}n\geq 0,$
(3) $\displaystyle so~{}that~{}t_{n+1}-t_{n}\geq D~{},$ (4)

where $D>0$ is a constant parameter of the QDTP algorithm that needs to be
fixed.

When the $n-th$ packet is transmitted by the SQF, we assume that it arrives
instantaneously at the Server’s input queue for AD processing. Here we are in
fact assuming that the physical transmission time from the SQF to the Server,
and the network protocol service time inside the Server, are tiny compared to
the AD processing duration $T_{n}$ at the Server. Thus the total delay $Q_{n}$
experienced by the $n$-th packet due to the action of the SQF, that elapses
from the arrival of the $n$-th packet to the SQF at $a_{n}$, until its arrival
to the AD at the Server at $t_{n}$, is:

$\displaystyle Q_{0}=t_{0}-a_{0}=0,$ (5) $\displaystyle
Q_{n+1}=t_{n+1}-a_{n+1},$ $\displaystyle=\max\\{t_{n}+D,a_{n+1}\\}-a_{n+1},$
$\displaystyle=0,~{}if~{}t_{n}+D\leq a_{n+1},~{}and$
$\displaystyle=t_{n}+D-a_{n+1},~{}otherwise.$ (6)

Since $t_{n}=Q_{n}+a_{n}$, we obtain the recursive expression:

$\displaystyle Q_{n+1}=[t_{n}+D-a_{n+1}]^{+},$
$\displaystyle=[Q_{n}+D-A_{n+1}]^{+},~{}n\geq 0,$ (7)

which is also an instance of Lindley’s equation (1).

On the other hand, the Server’s AD module also acts as a FCFS queue and we can
exploit Lindley’s equation again to compute $W_{n},~{}n\geq 0$ the waiting
time of the $n$-th packet that arrives to the Server to be processed for
attack detection, which is:

$\displaystyle W_{n+1}$ $\displaystyle=$
$\displaystyle[W_{n}+T_{n}-(t_{n+1}-t_{n})]^{+},~{}W_{0}=0,$ (8)
$\displaystyle\leq$ $\displaystyle W_{n}+T_{n}-(t_{n}-t_{n+1}),$ (9)

since the $n$-th packet’s AD service time is $T_{n}$ and the $n+1$-th
interarrival interval to the Server’s AD queue is $t_{n+1}-t_{n}$.

Therefore using equations (9) and (4) we obtain:

$W_{n+1}\leq W_{n}+T_{n}-D,$ (10)

and we have the following key insight into how to choose $D$:

Result 1. If we fix the parameter $D$ in the QDTP policy for the SQF to a
value so that $D>T_{n}$ for all $n\geq 0$, then the waiting time $W_{n}$ at
the Server will remain at the value $W_{n}=0$ for all $n\geq 0$.

We now present experiments showing the usefulness of Result 1. Noting from
Figure 4 that the measured average value of $T_{n}$ is $2.98~{}ms$ when there
is no attack, we first select $D=3~{}ms$ which is just above that value.

Figure 7 compares the case without SQF (Above) and with SQF (Below) during a
$60~{}sec$ UDP Flood Attack. Note that the figure Above represents the Server
queue length varying over time, without the SQF. The figure Below is in
logarithmic scale for the Server queue length, and compares the cases without
SQF (in Red) and with SQF (in Blue) for the Server queue length varying over
time. Since $D=3~{}ms$ is very close to the average of $T_{n}$, the
fluctuations in the values of $T_{n}$ cause a small queue buildup of a few
packets, as seen in the Blue plot in the lower part of the shown Below.

Figure 8 shows the results of four experiments where we measure the queue
length at the Server when a UDP Flood Attack lasts $30$ (Above) and $10$
(Below) seconds, without (Red) and with (Blue) the SQF Forwarder. Without SQF,
the Server’s AD processing time increases significantly. In the $30~{}sec$
attack, approximately $470,000$ packets are received at the Server and without
SQF it takes $44.45$ minutes for the Server to return to normal process them,
while in the $10~{}sec$ attack $153,667$ packets are received and it takes the
Server roughly $15$ minutes to process them. Note that in these curves it
takes some $99$ seconds for the compromised RPi to launch the attack.

Figure 9 shows that when we use the SQF based system with $D=2.7~{}ms$, which
is smaller than the value recommended by Result 1, when there is no attack
this choice of $D$ has very little effect. However, when a UDP Flood Attack
occurs, the Server’s AD processing is somewhat slowed down and the average
value of $T_{n}$ increases by roughly 10%.

On the other hand, Figure 10 confirms Result 1 since it shows that, if we take
$D=3.2~{}ms$ which guarantees that $D>T_{n}$ most of the time, then the
measured average value of $T_{n}$ remains at around $3~{}ms$ showing that it
has not been slowed down by the attack’s overload effect. Of course the same
is seen when no attack occurs.

##  IV SQF Queue Buildup and Attack Mitigation

When a Flood Attack occurs, the SQF accumulates packets in its input queue,
and forwards them to the Server using the QDTP algorithm with $D=3~{}ms$, so
that the Server does not experience any AD slowdown, ensuring that the Server
continues to operate as usual. Figure 11 shows the sudden increase and then
slow decrease of the SQF input queue when a UDP Flood attack lasts for $60$
seconds, and the SQF uses $D=3~{}ms$. Since both the SQF (and the Server) do
not drop packets, the attack packets will accumulate at the input queue of the
SQF.

Thus in this section we test a possible mitigating action that the SQF can
take. Since Flood Attacks are characterized by an unusually high packet
arrival rate, and this is also one of the attack detection parameters used by
the AD used in this work [26], we now test an additional feature, as follows:

  1. 1.

If the SQF receives more than $N$ packets in a time interval that is less than
or equal to $D$, it drops all incoming packets for the next $K.D$ time units.

  2. 2.

Here $N$ and $K$ are parameters of the mitigating action.

  3. 3.

The action is repeated as long as the condition 1) (above) on $N$ persists.

To illustrate the effect of this simple policy, in a first experiment we set
$N=10$ and $K=3$ and implement the suggested drop-based mitigation policy.

In this experiment, an RPi launches a $10$ second Flood Attack, and the
resulting queue length at the input of the SQF is shown in Figure 12, where we
see that the SQF input buffer reaches a small value of $12$ packets. The
attack starts at the $34$-th second and lasts $10$ seconds, but thanks to the
mitigation policy there is no accumulation of packets. After the attack ends
the SQF can continue to operate normally.

Figure 13 displays the queue length of the SQF input buffer in a second
experiment, when the attack lasts $60$ seconds, showing similar results to the
first experiment. Both measurements show the importance of having a simple
mitigating action to deal with high volume Flood Attacks.

However, although this policy appears attractive, it comes at the cost of
dropping legitimate (non-attack) packets that come from non-compromised IoT
devices.

Figure 12: SQF Smart QDTP Forwarder queue length when the attack lasts for
$10$ seconds, with the mitigation action, and $D=~{}3ms$. Figure 13: SQF queue
length with $D=3~{}ms$ when the attack lasts for $60$ seconds, and the
mitigation action is also applied.

##  V Conclusions

This paper has considered the effect of UDP Flood attacks on an IoT Server
that processes incoming traffic from IoT devices via a local area network. The
Server also incorporates an AD module.

We first show that such attacks, even when they last just a few seconds,
create overload for the IoT Server, so that its normal operations, including
AD, are substantially slowed down. We see, in particular, that a $60$ second
attack may create a backlog of packets at the Server, that may require several
hours to clear out.

Thus we propose that the Server’s input be “protected” by a special SQF front-
end that operates the QDTP policy, in order to allow the timely operation of
the Server even when an attack occurs. This approach requires that an
inexpensive lightweight hardware addition, such as an RPi, be installed
between the local area network that supports the IoT devices and the Server.
Several experiments are used to illustrate the effectiveness of the proposed
approach.

However, the SQF with its QDTP policy requires that a key timing parameter $D$
be chosen. Therefore, we provide a theoretical analysis of how $D$ should be
selected: we show that it must be just larger than the AD processing time of
the Server under normal, i.e. non-attack, conditions. We then validate this
observation with several experiments and show that the SQF can preserve the
Server from congestion and overload, and allow it to operate normally.
However, we note that the congestion that has been eliminated at the Server
may now accumulate at the SQF input, although this in itself does not stop the
RPi based SQF from continuing its normal operation.

Furthermore, when the incoming traffic rate is such that it clearly indicates
an attack, or when the Server informs the SQF that an attack is occuring, we
can implement a mitigating action at the SQF to drop incoming packets in
relatively short successive time intervals. This approach is tested
experimentally and shown to be effective. However, the fact that such a policy
may also drop incoming legitimate packets implies that there will be
circumstances when it cannot be used and a close coupling between AD at the
Server and packet drop actions at the SQF will be needed.

While this paper has focused on an architecture with multiple sources of IoT
traffic represented by several RPis, future work will consider Edge Systems
having multiple IoT servers, as well as multiple IoT devices and packet
sources, and will study the usage of dynamic policies for AD and traffic
routing at the edge for complex IoT Server and SQF gateway architectures.

Another important issue that should be addressed in future work is the energy
consumption of such edge systems [49], so that dynamic management policies may
be used to minimize energy consumption, as well as to optimize Quality of
Service and Cybersecurity.

## Acknowledgements

The support of the EU H2020 IoTAC Research and Innovation Action under GA No.
952684 is gratefully acknowledged.

## References

  * [1] Cisco, _Cisco Annual Internet Report (2018–2023)_ , Mar. 2020. [Online]. Available: https://www.cisco.com/c/en/us/solutions/collateral/executive-perspectives/annual-internet-report/white-paper-c11-741490.html
  * [2] A. Zanella, M. Zorzi, A. F. dos Santos, P. Popovski, N. Pratas, C. Stefanovic, A. Dekorsy, C. Bockelmann, B. Busropan, and T. A. Norp, “M2M massive wireless access: Challenges, research issues, and ways forward,” in _2013 IEEE Globecom Workshops_. IEEE, 2013, pp. 151–156.
  * [3] H. Jin, W. T. Toor, B. C. Jung, and J.-B. Seo, “Recursive pseudo-Bayesian access class barring for M2M communications in LTE systems,” _IEEE Transactions on Vehicular Technology_ , vol. 66, no. 9, pp. 8595–8599, 2017.
  * [4] J. Liu, L. Song _et al._ , “A novel congestion reduction scheme for massive machine-to-machine communication,” _IEEE Access_ , vol. 5, pp. 18 765–18 777, 2017.
  * [5] L. Tello-Oquendo, I. Leyva-Mayorga, V. Pla, J. Martinez-Bauset, J.-R. Vidal, V. Casares-Giner, and L. Guijarro, “Performance analysis and optimal access class barring parameter configuration in LTE-A networks with massive M2M traffic,” _IEEE Transactions on Vehicular Technology_ , vol. 67, no. 4, pp. 3505–3520, 2018.
  * [6] E. Johns and M. Ell, “Cyber security breaches survey 2023,” April 2023. [Online]. Available: https://www.gov.uk/government/statistics/cyber-security-breaches-survey-2023/cyber-security-breaches-survey-2023
  * [7] Cloudflare. [Online]. Available: https://www.cloudflare.com/learning/ddos/famous-ddos-attacks/
  * [8] R. Mahmoud, T. Yousuf, F. Aloul, and I. Zualkernan, “Internet of Things (IoT) security: Current status, challenges and prospective measures,” in _2015 10th International Conference for Internet Technology and Secured Transactions (ICITST)_. IEEE, 2015, pp. 336–341.
  * [9] M. U. Farooq, M. Waseem, A. Khairi, and S. Mazhar, “A critical analysis on the security concerns of Internet of Things (IoT),” _International Journal of Computer Applications_ , vol. 111, no. 7, pp. 1–6, 2015.
  * [10] “Hp study reveals 70 percent of Internet of Things devices vulnerable to attack,” Accessed on 25.01.2022. [Online]. Available: https://www.hp.com/us-en/hp-news/press-release.html%3Fid=1744676
  * [11] S. Evmorfos, G. Vlachodimitropoulos, N. Bakalos, and E. Gelenbe, “Neural network architectures for the detection of SYN flood attacks in IoT systems,” in _Proceedings of the 13th ACM International Conference on PErvasive Technologies Related to Assistive Environments_ , 2020, pp. 1–4.
  * [12] N. Statt, “How an army of vulnerable gadgets took down the web today,” October 2016. [Online]. Available: https://www.theverge.com/2016/10/21/13362354/dyn-dns-ddos-attack-cause-outage-status-explained
  * [13] B. Tushir, H. Sehgal, R. Nair, B. Dezfouli, and Y. Liu, “The impact of dos attacks onresource-constrained iot devices: A study on the mirai attack,” _arXiv preprint arXiv:2104.09041_ , 2021.
  * [14] H. Sinanović and S. Mrdovic, “Analysis of mirai malicious software,” in _2017 25th International Conference on Software, Telecommunications and Computer Networks (SoftCOM)_. IEEE, 2017, pp. 1–5.
  * [15] M. Antonakakis, et al., “Understanding the Mirai Botnet,” in _Proceedings of the 26th USENIX Security Symposium_ , 2017. [Online]. Available: https://www.usenix.org/conference/usenixsecurity17/technical-sessions/presentation/antonakakis
  * [16] Cloudflare. [Online]. Available: https://www.cloudflare.com/learning/ddos/udp-flood-ddos-attack/
  * [17] A. Iqbal, S. Aftab, I. Ullah, M. A. Saeed, and A. Husen, “A classification framework to detect DoS attacks.” _International Journal of Computer Network & Information Security_, vol. 11, no. 9, 2019.
  * [18] T. Chin, X. Mountrouidou, X. Li, and K. Xiong, “Selective packet inspection to detect DoS flooding using software defined networking (sdn),” in _2015 IEEE 35th International Conference on Distributed Computing Systems Workshops_ , 2015, pp. 95–99.
  * [19] C. Douligeris and A. Mitrokotsa, “DDoS attacks and defense mechanisms: classification and state-of-the-art,” _Computer networks_ , vol. 44, no. 5, pp. 643–666, 2004.
  * [20] M. Ell and R. Gallucci, “Cyber security breaches survey 2022,” July 2022. [Online]. Available: https://www.gov.uk/government/statistics/cyber-security-breaches-survey-2022/cyber-security-breaches-survey-2022
  * [21] H. Ahmetoglu and R. Das, “A comprehensive review on detection of cyber-attacks: Data sets, methods, challenges, and future research directions,” _Internet of Things_ , vol. 20, p. 100615, 2022. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S254266052200097X
  * [22] M. Banerjee and S. Samantaray, “Network traffic analysis based iot botnet detection using honeynet data applying classification techniques,” _International Journal of Computer Science and Information Security (IJCSIS)_ , vol. 17, no. 8, 2019.
  * [23] T. A. Tuan, H. V. Long, R. Kumar, I. Priyadarshini, N. T. K. Son _et al._ , “Performance evaluation of Botnet DDoS attack detection using machine learning,” _Evolutionary Intelligence_ , pp. 1–12, 2019.
  * [24] A. I. Al-issa, M. Al-Akhras, M. S. ALsahli, and M. Alawairdhi, “Using machine learning to detect DoS attacks in wireless sensor networks,” in _2019 IEEE Jordan International Joint Conference on Electrical Engineering and Information Technology (JEEIT)_ , 2019, pp. 107–112.
  * [25] E. Y. Güven and Z. Gürkaş-Aydın, “Mirai botnet attack detection in low-scale network traffic,” _Intelligent Automation & Soft Computing_, vol. 37, no. 1, pp. 419–437, 2023.
  * [26] E. Gelenbe and M. Nakıp, “Traffic based sequential learning during botnet attacks to identify compromised iot devices,” _IEEE Access_ , vol. 10, pp. 126 536–126 549, 2022.
  * [27] J. Mirkovic, S. Fahmy, P. Reiher, and R. K. Thomas, “How to test dos defenses,” in _2009 cybersecurity applications & technology conference for homeland security_. IEEE, 2009, pp. 103–117.
  * [28] M. Kaouk, F.-X. Morgand, and J.-M. Flaus, “A testbed for cybersecurity assessment of industrial and IoT-based control systems,” in _Lambda Mu 2018 - 21è Congrè de Maîtrise des Risques et Sûreté de Fonctionnement, Oct 2018, Reims, France_. [Online]. Available: https://hal.science/hal-02074654v1/document
  * [29] M. Annor-Asante and B. Pranggono, “Development of smart grid testbed with low-cost hardware and software for cybersecurity research and education,” _Wireless Personal Communications_ , vol. 101, pp. 1357–1377, 2018.
  * [30] O. A. Waraga, M. Bettayeb, Q. Nasir, and M. A. Talib, “Design and implementation of automated iot security testbed,” _Computers & security_, vol. 88, p. 101648, 2020.
  * [31] V. K. Singh, R. Sharma, and M. Govindarasu, “Testbed-based performance evaluation of attack resilient control for wind farm scada system,” in _2020 IEEE Power & Energy Society General Meeting (PESGM)_. IEEE, 2020, pp. 1–5.
  * [32] A. Ghaleb, S. Zhioua, and A. Almulhem, “Scada-sst: a scada security testbed,” in _2016 World Congress on Industrial Control Systems Security (WCICSS)_. IEEE, 2016, pp. 1–6.
  * [33] A. Tesfahun and D. L. Bhaskari, “A scada testbed for investigating cyber security vulnerabilities in critical infrastructures,” _Automatic Control and Computer Sciences_ , vol. 50, pp. 54–62, 2016.
  * [34] B. Reutimann and I. Ray, “Simulating measurement attacks in a scada system testbed,” in _Critical Infrastructure Protection XV: 15th IFIP WG 11.10 International Conference, ICCIP 2021, Virtual Event, March 15–16, 2021, Revised Selected Papers 15_. Springer, 2022, pp. 135–153.
  * [35] S.-U. Park and S.-M. Hwang, “Test bed construction for apt attack detection,” _International Journal of Control and Automation_ , vol. 11, no. 4, pp. 175–186, 2018.
  * [36] R. Arthi and S. Krishnaveni, “Design and development of iot testbed with ddos attack for cyber security research,” in _2021 3rd International Conference on Signal Processing and Communication (ICPSC)_. IEEE, 2021, pp. 586–590.
  * [37] A. P. Wright and N. Ghani, “A testbed for the evaluation of denial of service attacks in software-defined networks,” in _2019 SoutheastCon_. IEEE, 2019, pp. 1–6.
  * [38] P. V. Sontakke and N. B. Chopade, “Impact and analysis of denial-of-service attack on an autonomous vehicle test bed setup,” in _Proceedings of Third International Conference on Intelligent Computing, Information and Control Systems: ICICCS 2021_. Springer, 2022, pp. 221–236.
  * [39] M. Nasereddin, M. Nakip, and E. Gelenbe, “Measurement based evaluation and mitigation of flood attacks on a LAN test-bed,” in _Submitted for publication_ , 2023.
  * [40] E. Gelenbe and K. Sigman, “Iot traffic shaping and the massive access problem,” in _ICC 2022, IEEE International Conference on Communications, 16–20 May 2022, Seoul, South Korea_ , no. https://zenodo.org/record/5918301. https://zenodo.org/record/5918301, 2022, pp. 1–6.
  * [41] M. Nakip, V. Rodoplu, C. Güzeliş, and D. T. Eliiyi, “Joint forecasting-scheduling for the Internet of Things,” in _2019 IEEE Global Conference on Internet of Things (GCIoT)_ , 2019, pp. 1–7.
  * [42] V. Rodoplu, M. Nakıp, D. T. Eliiyi, and C. Güzeliş, “A multiscale algorithm for joint forecasting–scheduling to solve the massive access problem of iot,” _IEEE Internet of Things Journal_ , vol. 7, no. 9, pp. 8572–8589, 2020\.
  * [43] M. Nakip and E. Gelenbe, “Randomization of data generation times improves performance of predictive IoT networks,” in _2021 IEEE 7th World Forum on Internet of Things (WF-IoT)_ , 2021, pp. 350–355.
  * [44] S. Kumar and S. Rai, “Survey on transport layer protocols: Tcp & udp,” _International Journal of Computer Applications_ , vol. 46, no. 7, pp. 20–25, 2012.
  * [45] “MHDDoS - DDoS Attack Script With 56 Methods,” Online, May 2022, accessed: 2023-02-22. [Online]. Available: https://github.com/MatrixTM/MHDDoS
  * [46] K. Sigman, _Stationary Marked Point Processes: An Intuitive Approach_. Chapman and Hall, CRC, May 1995\.
  * [47] E. Gelenbe and I. Mitrani, _Analysis and Synthesis of Computer Systems, 2nd Ed._ World Scientific, Imperial College Press, 2010. [Online]. Available: https://doi.org/10.1142/p643
  * [48] E. Gelenbe, M. Nakip, D. Marek, and T. Czachorski, “Diffusion analysis improves scalability of IoT networks to mitigate the massive access problem,” pp. 182–189, November 2021.
  * [49] P. Fröhlich, E. Gelenbe, J. Fiołka, J. Checinski, M. Nowak, and Z. Filus, “Smart sdn management of fog services to optimize QoS and energy,” _Sensors_ , vol. 21, no. https://doi.org/10.3390/s21093105, p. 3105, 2021.