Density Estimation on a Network

TL;DR

This paper introduces a new nonparametric density estimation method on networks that adaptively models potential discontinuities at vertices, improving bias correction and scalability over existing methods.

Contribution

It proposes a data-driven, two-step local polynomial regression approach that handles density discontinuities at network vertices, with theoretical bias-variance analysis and practical extensions.

Findings

Reduces bias near vertices with discontinuities

Scales sub-linearly with sample size

Effective in simulation and real dendrite network data

Abstract

This paper develops a novel approach to density estimation on a network. We formulate nonparametric density estimation on a network as a nonparametric regression problem by binning. Nonparametric regression using local polynomial kernel-weighted least squares have been studied rigorously, and its asymptotic properties make it superior to kernel estimators such as the Nadaraya-Watson estimator. When applied to a network, the best estimator near a vertex depends on the amount of smoothness at the vertex. Often, there are no compelling reasons to assume that a density will be continuous or discontinuous at a vertex, hence a data driven approach is proposed. To estimate the density in a neighborhood of a vertex, we propose a two-step procedure. The first step of this pretest estimator fits a separate local polynomial regression on each edge using data only on that edge, and then tests for equality of the estimates at the vertex. If the null hypothesis is not rejected, then the second step re-estimates the regression function in a small neighborhood of the vertex, subject to a joint equality constraint. Since the derivative of the density may be discontinuous at the vertex, we propose a piecewise polynomial local regression estimate to model the change in slope. We study in detail the special case of local piecewise linear regression and derive the leading bias and variance terms using weighted least squares theory. We show that the proposed approach will remove the bias near a vertex that has been noted for existing methods, which typically do not allow for discontinuity at vertices. For a fixed network, the proposed method scales sub-linearly with sample size and it can be extended to regression and varying coefficient models on a network. We demonstrate the workings of the proposed model by simulation studies and apply it to a dendrite network data set.

Tables

Table 1. Table 1: Bias, standard deviation, and mean squared error of the proposed local (piecewise) polynomial regression compared to all existing methods. The result is based on 100 100 100 simulations of Cases I, III and III with 1000 1000 1000 data points on each edge. Bias, standard deviation and mean squared error are reported at the vertex as it is approached from e 2 subscript 𝑒 2 e_{2} in Figure  2 .

	Case I			Case II			Case III
	Bias	SD	MSE	Bias	SD	MSE	Bias	SD	MSE
LPLR	$-$0.0621	0.0593	0.0074	$-$0.0629	0.0423	0.0057	$-$0.0971	0.0123	0.0096
ESDK	$-$0.4481	0.0196	0.2011	$-$0.2711	0.0229	0.0741	$-$0.1062	0.0121	0.0114
ESCK	$-$0.4956	0.0184	0.2459	$-$0.2312	0.0301	0.0542	$-$0.1324	0.0100	0.0176
DE	$-$0.4670	0.0172	0.2180	$-$0.2438	0.0173	0.0597	$-$0.1140	0.0116	0.0131

Table 2. Table 2: Type II error, bias, standard deviation and mean squared error at the vertex approached from the right edge, based on 3000 3000 3000 simulations of 1000 1000 1000 data points per edge in each scenario.

$({\beta }_l,{\beta }_r)$	Type II error rate	Bias	SD	MSE
(3.5, 4.5)	0.0056	$-$0.0796	0.1480	0.0209
(3.55, 4.45)	0.0227	$-$0.0786	0.1482	0.0199
(3.6, 4.4)	0.0540	$-$0.0811	0.1458	0.0215
(3.65, 4.35)	0.1283	$-$0.0833	0.1453	0.0207
(3.7, 4.3)	0.2453	$-$0.0861	0.1487	0.0210
(3.75, 4.25)	0.3906	$-$0.0925	0.1454	0.0203
(3.8, 4.2)	0.6030	$-$0.1017	0.1398	0.0183
(3.85, 4.15)	0.7633	$-$0.1040	0.1245	0.0158
(3.9, 4.1)	0.8863	$-$0.0956	0.1135	0.0136
(3.95, 4.05)	0.9460	$-$0.0781	0.1038	0.0127
(4, 4)	NA	$-$0.0631	0.1013	0.0140

Equations

$$\hat{f}(x)=\frac{1}{N}\sum _{i=1}^N{K}_h(x_i-x),$$

$$E(y_i|x_i)=E\left(\frac{c_i}{N\omega }\right)=\frac{1}{\omega }E({\hat{p}}_i)=\frac{1}{\omega }{p}_i\approx m(x_i),$$

$$Var(y_i|x_i)=Var\left(\frac{c_i}{N\omega }\right)=\frac{1}{{\omega }^2}Var({\hat{p}}_i)=\frac{p_i-p_i^2}{N{\omega }^2}\approx \frac{1}{N\omega }m(x_i)+\frac{1}{N}m(x_i{)}^2.$$

$$m(x)$$

$$=m(x_0)+{\beta }_1(x-x_0)+{\beta }_2(x-x_0{)}^2+\cdots +{\beta }_p(x-x_0{)}^2,$$

$$\displaystyle\sum^{n}_{i=1}\big{\{}y_{i}-\beta_{0}-\beta_{1}(x_{i}-x_{0})-\dots-\beta_{p}(x_{i}-x_{0})^{p}\big{\}}^{2}K_{h}\left({x_{i}-x_{0}}\right),$$

$$\mathbf{Y}=\left[\begin{array}{c}{y}_1\\ ⋮\\ y_n\end{array}\right],\mathrm{and}\mathbf{X}=\left[\begin{array}{cccc}1& x_1-x& \dots & (x_1-x{)}^p\\ ⋮& ⋮& \ddots & ⋮\\ 1& x_n-x& \dots & (x_n-x{)}^p\end{array}\right].$$

$$\mathbf{X}=\left[\begin{array}{cc}1& x_1-x\\ ⋮& ⋮\\ 1& x_n-x\end{array}\right].$$

$$E\left\{\hat{m}(x)-m(x)\right\}=\frac{1}{2}{h}^2\left\{\frac{{\left({\sigma }_K^2\right)}^2-\left({\sigma }_K^3\right){\sigma }_K^1}{{\sigma }_K^2-{\left({\sigma }_K\right)}^2}\right\}{m}^{(2)}(x)+O\left(\omega \right)+o\left(h^2\right)$$

$$\mathrm{Var}\left\{\hat{m}(x)\right\}=C(N,h,\omega ,x)Q_K+o\left(\frac{1}{Nh}\right),$$

$$C(N,h,\omega ,x)=\frac{1}{Nh}m(x)+\frac{\omega }{Nh}m(x{)}^2\mathrm{and}{Q}_K=\frac{R_K^0{\left({\sigma }_K^2\right)}^2-2R_K^1{\sigma }_K^2{\sigma }_K^1+R_K^2{\left({\sigma }_K^1\right)}^2}{{\left\{{\sigma }_K^2-{\left({\sigma }_K^1\right)}^2\right\}}^2}.$$

$$\hat{m}(x)-m(x)-\frac{1}{2}{h}^2\left\{\frac{{\left({\sigma }_K^2\right)}^2-\left({\sigma }_K^3\right){\sigma }_K^1}{{\sigma }_K^2-{\left({\sigma }_K^1\right)}^2}\right\}{m}^{(2)}(x)\stackrel{d}{\to }N\left\{0,C(N,h,\omega ,x)Q_K\right\}.$$

$$E\left\{\hat{m}(x)-m(x)\right\}=\frac{1}{2}{h}^2{\sigma }_K^2{m}^{(2)}(x)+O\left(\omega \right)+o\left(h^2\right)$$

$$\mathrm{Var}\left\{\hat{m}(x)\right\}=C(N,h,\omega ,x)R_K^0+o\left(\frac{1}{Nh}\right).$$

$$H_0:m_{e_1}(x)=m_{e_2}(x)=\cdots =m_{e_J}(x)$$

$$H_1:\mathrm{Not}{H}_0.$$

$$(\mathbf{C}\hat{\mathbf{m}}{)}^T{\mathbf{\Sigma }}^{-1}\mathbf{C}\hat{\mathbf{m}}\stackrel{a}{\sim }{\chi }_{J-1}^2,$$

$$m_{e_2}(x_2)$$

$$=m_{e_1}(v)+m_{e_2}^{(1)}(x)(x_2-v)$$

$$\approx m_{e_1}(x)+m_{e_1}^{(1)}(x)(v-x)+m_{e_2}^{(1)}(v)(x_2-v),$$

$$m_{e_3}(x_3)\approx m_{e_1}(x)+m_{e_1}^{(1)}(x)(v-x)+m_{e_3}^{(1)}(v)(x_3-v).$$

$$\displaystyle\sum_{x_{il}\in e_{l}}\big{\{}y_{il}-\beta_{0}+\beta_{1}(e_{l})(x_{il}-x)\big{\}}^{2}K_{il}+\sum_{j\neq l}\sum_{x_{ij}\in e_{j}}\big{\{}y_{ij}-\beta_{0}-\beta_{1}(e_{l})(v-x)+\beta_{1}(e_{j})(x_{ij}-v)\big{\}}^{2}K_{ij},$$

$${\left(\mathbf{Y}-\mathbf{X}\mathbf{\beta }\right)}^T\mathbf{W}\left(\mathbf{Y}-\mathbf{X}\mathbf{\beta }\right),$$

$$\hat{m}(x)={\mathbf{e}}_1^T{\left({\mathbf{X}}^T\mathbf{W}\mathbf{X}\right)}^{-1}{\mathbf{X}}^T\mathbf{W}\mathbf{Y}.$$

$$E\left\{\hat{m}(x)-m(x)\right\}=\frac{1}{2}{\mathbf{e}}_1^T{\mathbf{U}}_0^{-1}\left\{h^2{\mathbf{R}}_0+h(v-x)\left({\mathbf{R}}_1-{\mathbf{U}}_1{\mathbf{U}}_0^{-1}{\mathbf{R}}_0\right)\right\}+o\left(h(v-x)\right)+o\left(h^2\right)+O(\omega ),$$

$$and$$

$$\mathrm{Var}(\hat{m}(x))=\frac{1}{Nh}C(\omega ,x){\mathbf{e}}_1^T{\mathbf{U}}_0^{-1}\left[{\mathbf{M}}_0+\left(\frac{v-x}{h}\right)\left({\mathbf{M}}_0{\mathbf{U}}_0^{-1}{\mathbf{U}}_1+{\mathbf{M}}_1-{\mathbf{U}}_1{\mathbf{U}}_0^{-1}{\mathbf{M}}_0\right)\right]{\mathbf{U}}_0^{-1}{\mathbf{e}}_1$$

$$+o\left(\frac{1}{Nh}\right)+o\left(\frac{\omega }{Nh}\right)+o\left(\frac{v-x}{N{h}^2}\right)+o\left(\frac{\omega (v-x)}{N{h}^2}\right),$$

$$m_{e_2}(x_2)$$

$$=m_{e_1}(v)+m_{e_2}^{(1)}(v)(x_2-v)$$

$$\approx m_{e_1}(x)+m_{e_1}^{(1)}(x)(v-x)+m_{e_2}^{(1)}(v)(x_2-v)$$

$$\approx m_{e_1}(x)+m_{e_1}^{(1)}(x)(v-x)+m_{e_1}^{(1)}(v)(x_2-v).$$

$$m_{e_1}(x_1)\approx m_{e_1}(x)+(x_1-x)m_{e_1}^{(1)}(x)$$

$$m_{e_2}(x_2)\approx m_{e_2}(v_1)+(x_2-v_1)m_{e_2}^{(1)}(v_1)$$

$$m_{e_3}(x_3)\approx m_{e_3}(v_2)+(x_3-v_2)m_{e_3}^{(1)}(v_2).$$

$$m_{e_2}(x_2)\approx m_{e_1}(x)+(v_1-x)m_{e_1}^{(1)}(x)+(x_2-v_1)m_{e_2}^{(1)}(v_1)$$

$$m_{e_3}(x_3)\approx m_{e_2}(v_1)+(v_2-v_1)m_{e_2}^{(1)}(v_1)+(x_3-v_2)m_{e_3}^{(1)}(v_2)$$

$$\approx m_{e_1}(x)+(v_1-x)m_{e_1}^{(1)}(x)+(v_2-v_1)m_{e_2}^{(1)}(v_1)+(x_3-v_2)m_{e_3}^{(1)}(v_2).$$

Full text

Keywords: Asymptotic bias and variance, discontinuous density, kernel density estimation, local piecewise linear estimation, pretest estimation

1 Introduction

In this paper, we analyze spatial patterns of points that lie on a network. There are many types of events that can occur along networks. The lines or curves that form the network may be roads, rivers, subway lines, spider webs, blood vessels, or nerve fibers. The events may be vehicles, street crimes, car accidents, retail stores, or neuroanatomical features such as dendritic spines. For a recent summary of applications that involve the modelling of point patterns on a network of lines, see [2]. In order to carry out analyses of those events, researchers need a range of specific techniques.

One of the frequently demanded tasks is density estimation on a network, but few statistical methods had been developed to address this need until recently. A natural first attempt at analyzing such data is to take the kernel density estimate (KDE) on the one-dimensional real line

[TABLE]

and apply it directly to network data by defining $|x_{i}-x|$ as the network distance, where $x$ is any location on the network, and $x_{1},\dots,x_{N}$ are the observed data locations. $K_{h}(\cdot)=K(\cdot/h)/h$, where $K$ is a kernel function, and $h$ is the bandwidth. The network distance is defined as the shortest path distance between two points on a network. However, under this approach, estimate (1) does not conserve mass. This happens because that the induced kernel $K_{h}(\cdot)$ is not a probability density on the network at points within the $h$-neighborhood of a vertex, and so (1) is not a probability density. As a result, it will overestimate the true density. [8] summarized widely used methods for density estimation on a network. Many published papers, such as [13], mentioned by Okabe et al. computed a kernel density estimate on a network, but most have used (1).

Although there is currently no consensus view on how to perform density estimation on a network, the two most popular heuristic techniques, namely the equal-split discontinuous kernel (ESDK) estimator and the equal-split continuous kernel (ESCK) estimator ([9], [8], [12]) seem to produce reasonable results in applications. However their methods suffer from a lack of theoretical grounding, and their computational cost is high (see Section 7.2). Furthermore, ESCK does not allow for discontinuity at vertices, which is found in many applications, such as traffic network, and the example in Figure 1 and Section 7.4. Although ESDK produces discontinuous density at vertices, the estimation in the small neighborhood of vertices always uses data from all edges, which can lead to large bias. This is demonstrated in Case I of the simulation study in Section 7.1.

More recently [7] proposed a diffusion density estimator (DE) on networks using the connection between kernel smoothing and diffusion [3]. DE is based on numerical solution of the heat equation, and it is not only faster than ESDK and ESCK, but also provides a sound statistical rationale and helps establish theoretical properties. The estimate is asymptotically unbiased with rate $h^{2}$ away from the boundary region, where $h$ is the bandwidth. At a terminal vertex, it has the standard boundary rate $h$ of KDE. Moreover, similar to ESCK, DE does not allow for discontinuity at vertices, which consequently, can lead to large bias.

In this paper, we take a different approach. We propose local piecewise linear regression (LPLR) on a network, and establish theoretical properties for the proposed method. The piecewise linear structure is needed to model the discontinuity in the regression function or the first derivative of the regression function at a vertex. By binning, we form a histogram of observations (bin centers and counts) on the network, which is then smoothed by local polynomial regression (LPR). The advantage of local polynomial estimation is reduction of the order of the bias, especially when the evaluation point is near or on the boundary vertex, and, as a practical matter, the amount of computation is reduced. To accommodate densities with discontinuous derivatives at vertices, we introduce piecewise local regression. As far as we know, local polynomial regression on networks has not be proposed and investigated and piecewise local linear regression is also new. We note that the proposed method assumes that the data comes from a Poisson process on a network, but McSwiggan et al is applicable to more general point processes. However we want to point out that our approach has better order of the bias near vertices, and it only imposes continuity at vertices when there is evidence suggested by the data that the density is continuous there. Our approach is not limited to density estimation as it can be extended to handle other types of analysis of network data, such as regression analysis. The implementation of the proposed method is fast, the computation time scales only with data sub-linearly.

To motivate our research and illustrate its contribution, we introduce the dendrite data collected by the Kosik Lab, UC Santa Barbara, and first analyzed by [1] and [6]. In this example, the events are dendritic splines, which are of clinical importance. Cognitive disorders such as ADHD and autism may result from abnormalities in dendritic spines, especially the number of spines and their maturity [10]. The events on the network are the locations of 566 spines observed on one branch of the dendritic tree of a rat neuron, as shown in Figure 1. The density of the spines shows clear discontinuity at vertices A, C and D; see Figures 7, 9, 11 and 13. We will need a density estimation method that allows for multiple levels of smoothness at vertices – discontinuous, continuous with discontinuous derivative, and continuous with continuous derivative. We show that our proposed method will remove the bias that has been noted for existing methods such as [8] and [7] due to their inflexibility at vertices. KDE, ESDK, ESCK and DE’s inflexibility arises because the smoothness of the density function at a vertex needs to be decided before choosing an estimator: if discontinuity is desired at a vertex, KDE is applied to each edge, otherwise ESDK, ESCK or DE should be applied. In contrast, we propose a data-driven approach and determine the smoothness of the density function at a vertex by statistical testing. In fact, estimation at vertices is a major contribution of this paper.

The paper is organized as follows. Section 2 gives basic definitions. In Section 3, a binned local polynomial regression estimator is constructed individually on each edge, and asymptotic properties of the estimator are established. As the evaluation point approaches the vertex from each edge, we obtain limits of the regression functions. In Section 4, we construct an asymptotic test for the equality of those limits. If the null hypothesis is not rejected, then the regression function is re-estimated for evaluation points that are within the $h$-neighborhood of a vertex by binned local piecewise polynomial regression using data from all neighboring edges. These estimators are constructed and their asymptotic properties are studied in Section 5. In Section 6, we discuss networks with loops and related computational issues when the bandwidth is large. In Section 7, we demonstrate the workings of the proposed method with simulated data and a real example, and discuss and analyze practical issues. We end with a discussion in Section 8.

2 Preliminaries

In this section, we introduce basic definitions and formulate density estimation as a local linear regression problem.

2.1 Network and density function on a network

Let $v_{1},v_{2}\in\mathbb{R}^{d}$ be vertices and define an edge connecting $v_{1}$ and $v_{2}$ as a curve $e:[0,1]\to\mathbb{R}^{d}$ where $e(0)=v_{1}$ and $e(1)=v_{2}$. We call $L=\cup^{J}_{i=j}e_{j}$, $0<J<\infty$ a network in $\mathbb{R}^{d}$. A path between two points $x,y\in L$ is a sequence $\pi=(x,v_{1},\dots,v_{n},y)$, where $v_{1},\dots,v_{n}$ are vertices, such that $e_{j}=[v_{j},v_{j+1}]$ are edges, and $x\in e_{0}$, $y\in e_{n}$. The path length is Length$(\pi)=\int_{x}^{v_{1}}|e_{0}^{\prime}(t)|dt+\sum_{j=1}^{n-1}\int^{v_{j+1}}_{v_{j}}|e_{j}^{\prime}(t)|dt+\int_{v_{n}}^{y}|e_{n}^{\prime}(t)|dt$, where integration is with respect to arc length. In the case where there are multiple paths between two points, we will use the shortest path to define length. If there no paths from $x$ to $y$, we define the network distance between $x$ and $y$ to be infinity. Often, the network is embedded in a two-dimensional space, such as a street network (if there are no overpasses). However all of our results apply easily to a network of curves embedded in a higher dimensional space, such as the dendrites which are embedded in $\mathbb{R}^{3}$. In the theory section (Section 3, 4 and 5), we will study in deal a simple network, where there is one vertex connecting $J$ edges. We discuss application to complex networks in Section 6 and 7.

Our objective is to construct an estimate of an unobservable underlying probability density function $f(x)$ of event locations based on observed data $x_{1},\cdots,x_{n}\in L$. The density satisfies $f(x)\geq 0$ for all $x\in L$, and $\int_{L}f(x)dx=\sum_{j=1}^{J}\int_{e_{j}}f(x)dx=1$. Stated differently, let $\boldsymbol{x}=\{x_{1},\dots,x_{n}\}$ be a realization of a Poisson process $\boldsymbol{X}$ on the network. We want to estimate the intensity function $\lambda(x)$ of $\boldsymbol{X}$. The intensity function $\lambda(x)$ is the expected number of random points per unit length of network in the small neighborhood of $x$. The estimation of $f$ and of $\lambda$ are equivalent because, for fixed $N$, $\lambda(x)=Nf(x)$ for all $x\in L$.

2.2 Problem formulation

In this paper, we consider density estimation via local polynomial regression. This is achieved by way of binning. Here we use the “simple binning” discussed in [5]. For $i=1,\dots,n$, $x_{i}$ are the bin centers, $c_{i}$ are bin counts, $y_{i}$ are bin heights, and $\omega$ is bin width. Letting $N$ be the total number of observations on the linear network, we have $y_{i}={c_{i}}/{N\omega}$. These rectangles form a histogram with total area $1$, since the area of the $i$th bin is $c_{i}/N$ and $\sum_{i=1}^{n}c_{i}/N=1$. For a chosen $\omega$, we consider the regression model $y_{i}=m(x_{i})+\epsilon(x_{i})$. We have

[TABLE]

where $p_{i}$ is the expected proportion of points in bin $i$, and $m(\cdot)$ is the regression function. We note that this assumption would still hold for more general point processes. For the variance, $\hat{p}_{i}=c_{i}/N$ is the sample proportion so

[TABLE]

The variance expression depends on the Poisson process assumption. The approximate equality symbol is due to binning induces a bias, but the bias can be made negligible by a small bin width. Hall and Wand studied the accuracy of binned kernel density estimators [5], and their results show that, with the commonly-used Epanechnikov kernel, we require $\omega$ to go to zero faster than $h$ the bandwidth, if binning is not to have a significant effect on the bias of the estimator. Therefore, the error term $\epsilon(x_{i})$ has approximately mean zero and variance $m(x_{i})/N\omega+m(x_{i})^{2}/N$.

Assuming that $\omega$ is sufficiently small, we can smooth the histogram by local polynomial regression using $(x_{i},y_{i})$, $i=1,\dots,n$, as data. We propose a “pretest” estimation procedure that consists of the following steps:

Local polynomial regression on each edge individually. 2. 2.

Test joint equality of intercepts at the vertex. 3. 3.

If joint equality is not rejected, then the regression function is re-estimated by local piecewise polynomial regression using data from all neighboring edges for evaluation points that are within the $h$-neighborhood of a vertex.

An asymptotic test for joint equality of slopes at the vertex can be constructed similarly.

3 Local polynomial regression

First, we consider each edge individually. This is equivalent to a fixed equally-spaced design model, where the $x$-variables are the bin centers $x_{1},\dots,x_{n}$, and the $y$-variables are the bin heights $y_{1},\dots,y_{n}$. We want to estimate the regression function $m(x)=E(Y|X=x)$. Using a Taylor expansion, we can approximate $m(x)$, where $x$ is close to a point $x_{0}$, by a $p$ degree polynomial:

[TABLE]

provided that all the required derivatives exist. The local polynomial regression estimator minimizes with respect to $\beta_{0},\beta_{1},\dots,\beta_{p}$ the function

[TABLE]

where $K_{h}\left({x_{i}-x_{0}}\right)=K(\frac{x_{i}-x_{0}}{h})/h$. Let $\hat{\boldsymbol{\beta}}=\{\hat{\beta_{0}},\hat{\beta_{1}},\dots,\hat{\beta_{p}}\}$ denote the minimizer of (2). Then $\hat{\beta_{0}}$ estimates $m(x)$ and $s!\hat{\beta_{s}}$ estimates $m^{(s)}(x)$, the $s$th derivative of $m(x)$, for $s=1,\dots,p$. Conveniently, (2) is a standard weighted least-squares regression problem. Let \boldsymbol{W}={\rm diag}\big{\{}K_{h}(x_{1}-x),\dots,K_{h}(x_{n}-x)\big{\}}. Define

[TABLE]

Assuming invertibility of $\boldsymbol{X}^{\prime}\boldsymbol{W}\boldsymbol{X}$, then $\hat{\boldsymbol{\beta}}=\left(\boldsymbol{X}^{T}\boldsymbol{W}\boldsymbol{X}\right)^{-1}\boldsymbol{X}^{T}\boldsymbol{W}\boldsymbol{Y}$. The estimate of the regression function at $x$ is $\hat{m}(x)=\hat{\beta}_{0}=\boldsymbol{e}^{T}_{1}\boldsymbol{\hat{\beta}}$, where $\boldsymbol{e}$ is the $(p+1)\times 1$ vector with $1$ being the first entry and zero elsewhere. [11] studied the leading conditional bias and variance of the above estimator, but here we consider in detail only the local linear $(p=1)$ binned estimator. We make the following assumptions:

1. A1.

The function $m^{(2)}(\cdot)$ is continuous. 2. A2.

The kernel $K$ is symmetric and supported on $(-1,1)$. Also, $K$ has a bounded first derivative. 3. A3.

As $N\to\infty$ and $\omega\to 0$, $Nh\to\infty$ and $\omega=o(h^{2})$, where $N$ is sample size, $h$ is bandwidth and $\omega$ is bin width.

It follows from the definition of the estimator that $E\{\hat{m}(x)\}\approx\boldsymbol{e}^{T}_{1}\left(\boldsymbol{X}^{T}\boldsymbol{W}\boldsymbol{X}\right)^{-1}\boldsymbol{X}^{T}\boldsymbol{W}\boldsymbol{M}$, where the vector $\boldsymbol{M}=\left\{m(x_{1}),\dots,m(x_{n})\right\}^{T}$ contains the true regression function values at the each of the $x_{i}$’s. For local linear regression we have that

[TABLE]

Also, ${\rm Var}\{\hat{m}(x)\}=\boldsymbol{e}^{T}_{1}\left(\boldsymbol{X}^{T}\boldsymbol{W}\boldsymbol{X}\right)^{-1}\left(\boldsymbol{X}^{T}\boldsymbol{W}\boldsymbol{V}\boldsymbol{W}\boldsymbol{X}\right)\left(\boldsymbol{X}^{T}\boldsymbol{W}\boldsymbol{X}\right)^{-1}\boldsymbol{e}_{1}$, where $\boldsymbol{V}={\rm Var}(\boldsymbol{Y})$ is a diagonal matrix with diagonal entries $m(x_{i})/\omega N-m(x_{i})^{2}/N$, $i=1,\dots,n$. Note that if $m$ is a linear function then the local linear estimator is exactly unbiased. To find the leading bias term for general function $m$, we review the result of local linear regression [4] with the binning procedure discussed in [5]. This estimator uses only data on a single edge.

Theorem 1**.**

Suppose that $x$ is a point on the line segment of interest, and that A1, A2 and A3 hold. Let $\hat{m}(x)=\boldsymbol{e}^{T}_{1}\left(\boldsymbol{X}^{T}\boldsymbol{W}\boldsymbol{X}\right)^{-1}\boldsymbol{X}^{T}\boldsymbol{W}\boldsymbol{Y}$, then

[TABLE]

where $\sigma^{i}_{K}=\int u^{i}K(u)du$, $R^{i}_{K}=\int u^{i}K(u)^{2}du$,

[TABLE]

The estimator has the following asymptotic distribution:

[TABLE]

For an interior point $x$, we can simplify (4) and (5). Note that in this case $\sigma^{i}_{K}=0$ for odd $i$, we have

[TABLE]

The local linear estimator has better properties at the boundary than the kernel density estimator (KDE). Asymptotically, the local linear estimator’s bias is of the same order, $h^{2}$, at a boundary point as at an interior point, where the bias would be of order $h$ for KDE. Even if $x$ is at the boundary of the density’s support, since the local linear estimator fits a weighted least squares line through data near the boundary, if the true relationship is linear, this estimator will be exactly unbiased.

Extension to higher order $p$ is straightforward, see [11]. For $p$ odd, the bias is of order $h^{p+1}$ everywhere. When $p$ is even, the bias is of order $h^{p+2}$ away from the boundary but $h^{p+1}$ in the boundary region. By increasing the polynomial order from even to the next odd number, the order of the bias in the interior remains unchanged, but the bias simplifies. Kernel estimation corresponds to $p=0$. On the other hand, by increasing the polynomial order from odd to the next even number, the bias order increases in the interior. This effect is analogous to the bias reduction achieved by higher-order kernels.

4 Test for joint equality at vertex

When multiple edges meet at a vertex, there are as many estimates for $m(v)$. Consider a simple network consisting of three edges $e_{1},e_{2}$ and $e_{3}$ and a vertex $v$ (see Figure 2). One estimate is $\hat{m}_{e_{1}}(v)$, the local linear regression estimate of $m_{e_{1}}(v):=\lim_{x\to v,\ x\in e_{1}}m(x)$ using only data on $e_{1}$. The estimates $\hat{m}_{e_{2}}(v)$ and $\hat{m}_{e_{3}}(v)$ are constructed similarly. In this way, for a vertex connecting $J$ edges, we obtain $J$ estimates at that vertex. In this section, we construct a test to study whether the $m_{e_{j}}(v)$ for $j=1,\dots,J$ (or any subset) are equal.

Let $v$ be a vertex and let $e_{l}$, $l=1,\dots,J$, be the edges connected by $v$. We want to test

[TABLE]

Theorem 2**.**

Let $\hat{\boldsymbol{m}}=\{\hat{m}_{e_{1}}(v),\dots,\hat{m}_{e_{J}}(v)\}^{T}$, and $\boldsymbol{C}$ be a contrast matrix such that $\boldsymbol{C1}=\boldsymbol{0}$. Under the null hypothesis, we have

[TABLE]

where $\boldsymbol{\Sigma}=\boldsymbol{C}\,diag(V_{e_{1}},\dots,V_{e_{J}})\,\boldsymbol{C}^{T}$, and $V_{e_{l}}$ is the asymptotic variance of the $i$th estimator. Here $\stackrel{{\scriptstyle a}}{{\sim}}$ denotes “asymptotically distributed as.”

We note that the test statistic (12) is invariant under the choice of contrast matrices. This is easy to see by noticing that the rows of $\boldsymbol{C}$ are linearly independent.

5 Estimation with equal intercepts at the vertex

If $H_{0}$ is not rejected, then evaluation points within the $h$-neighborhood of the vertex can be estimated using data from all neighboring edges, subject to $m_{e_{1}}(v)=\dots=m_{e_{J}}(v)$. The resulting estimator has a lower variance compared to estimating separately on each edge.

For now, we make the following assumptions:

1. B1.

No loop shorter than $2h$ in the network. 2. B2.

All edges are longer than $h$.

If the $h$ neighborhood of the evaluation point $x$ is completely covered by an edge, then the estimator is the same as if we only use data from that edge, and its asymptotic bias and variance are given in Theorem 1. Next, we derive the estimator when the $h$ neighborhood $x$ contains a vertex. Note that, by assumption B1, there is a unique path between the evaluation point and any data point that is within its $h$-neighborhood, and by assumption B2, we need only consider the case where the neighborhood of $x$ contains exactly one vertex. For a fixed network, the above assumptions eventually will hold as $N\to\infty$ and $h\to 0$. We will consider the effect of these assumptions on the implementation of the proposed estimator in Section 6.

We note that our estimator will have an additional bias when the test in Section 4 makes a type II error, that is, when the density is not continuous at the vertex but the test accepts that it is continuous. We will investigate this problem in Section 7.3.

5.1 Deriving the estimator

Let us first consider the simple network in Figure 2, and consider a point $x\in e_{1}$ in the $h$-neighborhood of the vertex $v$, and data points $(x_{1},y_{1})$ on $e_{1}$, $(x_{2},y_{2})$ on $e_{2}$ and $(x_{3},y_{3})$ on $e_{3}$. Directional distance is used here. The evaluation point $x$ is the origin when measuring distance. Note that we will use directional distance in our proposed model. When we consider a path from a point on the network to another, we use the evaluation point $x$ as the origin. For example the distance between $x_{1}$ and $x$ is $x_{1}-x$ and the distance between $x_{1}$ and $v$ is $x_{1}-v$, which are both negative.

Using a Taylor expansion up to first order, we can approximate the regression function on $e_{1}$ at $x_{1}$ as $m_{e_{1}}(x_{1})\approx m_{e_{1}}(x)+m^{(1)}_{e_{1}}(x)(x_{1}-x)$ provided that all the required derivatives exist. For the regression function on $e_{2}$ at $x_{2}$, we have

[TABLE]

where (13) holds because that the regression functions from all edges are assumed to be equal at the vertex, (14) is a Taylor expansion of $m_{e_{1}}(v)$ at $x$ on $e_{1}$, and $m^{(1)}_{e_{j}}(v)$ is defined as $\lim_{x\to v,x\in e_{j}}m_{e_{j}}^{(1)}(x)$. Similarly we have

[TABLE]

For a simple network with one vertex $v$ of degree $J$ and an evaluation point $x$ on $e_{l}$, we have for $x_{i}\in e_{j}$, $j=1,\dots,J$ and $j\neq l$, $m_{e_{j}}(x_{i})\approx m_{e_{l}}(x)+m^{(1)}_{e_{l}}(x)(v-x)+m^{(1)}_{e_{j}}(v)(x_{i}-v)$, and for $x_{i}\in e_{l}$, $m_{e_{l}}(x_{i})\approx m_{e_{l}}(x)+m^{(1)}_{e_{l}}(x)(x_{i}-x)$.

Under the above construction, the regression functions agree at the vertex but can have different slopes there. Letting $\beta_{0}=m_{e_{l}}(x)$, $\beta_{1}(e_{l})=m^{(1)}_{e_{l}}(x)$, and $\beta_{1}(e_{j})=m^{(1)}_{e_{j}}(v)$ for $j\neq l$, these parameters are estimated by minimizing:

[TABLE]

where $K_{ij}=K\left(d_{ij}/h\right)$, and $d_{ij}$ is the network distance, i.e., the arc length, between $x_{ij}$ and $x$. In matrix notation,  we have

[TABLE]

where $\boldsymbol{\beta}=[\beta_{0},\beta_{1}(e_{1}),\dots,\beta_{1}(e_{J})]^{T}$ and $\boldsymbol{X}$ is a $n\times(J+1)$ matrix. The first column of $\boldsymbol{X}$ is identically $1$. In the $j$th column, for $j>1$ and $j\neq l$, the $1+\sum_{i=1}^{j-1}n_{i}$ to $\sum_{i=1}^{j}n_{i}$ entries are $\{x_{ij}-v\}_{i=1,\dots,n_{j}}$, where $n_{i}$ is the number of bin centers on edge $i$, and the remaining entries are zero. In the $l$th column of $\boldsymbol{X}$, the $1+\sum_{i=1}^{l-1}n_{i}$ to $\sum_{i=1}^{l}n_{i}$ entries are $\{x_{il}-v\}_{i=1,\dots,n_{l}}$, and the remaining entries are $v-x$. We also note that $\boldsymbol{Y}=\{y_{ij}\}_{i=1,\dots,n_{j};j=1\dots,J}$, and $\boldsymbol{W}={\rm diag}\{K_{ij}\}_{i=1,\dots,n_{j};j=1\dots,J}$.

The solution to the minimization problem (15) is $\hat{\boldsymbol{\beta}}=\left(\boldsymbol{X}^{T}\boldsymbol{W}\boldsymbol{X}\right)^{-1}\boldsymbol{X}^{T}\boldsymbol{W}\boldsymbol{Y}$, and the estimate for the regression function at $x$ under the constraint $m_{e_{1}}(v)=\dots=m_{e_{J}}(v)$ is given by

[TABLE]

Note that we can extend the above derivation to quadratic and higher order approximation, allowing different $m^{(p)}_{e_{j}}(v)$ for all $j$ under the constraint that $m^{p^{\prime}}_{e_{j}}$ are equal for $p^{\prime}<p$.

5.2 Asymptotic properties – the local linear case ($p=1$)

In this subsection, we make the following assumptions:

1. C1.

The functions $m_{e_{j}}(v)$ are equal for all $j$. 2. C2.

The function $m^{(2)}_{e_{j}}(\cdot)$ for all $j$ are continuous on $e_{j}$, excluding the vertex. 3. C3.

The kernel $K$ is symmetric and supported on $(-1,1)$. Also, $K$ has a bounded first derivative. 4. C4.

The vertex $v$ is always in the $h$-neighborhood of $x$ as $h$ approaches 0, so $x$ must approach $v$ at least as fast as $h$ approaches 0. 5. C5.

As $N\to\infty$ and $\omega\to 0$, $Nh\to\infty$ and $\omega=o(h^{2})$, where $N$ is sample size, $h$ is bandwidth, and $\omega$ is bin width.

Note that if C4 does not hold, then we are in the case already studied by [11].

Theorem 3**.**

Suppose that $x\in e_{l}$ is within the $h$-neighborhood of a vertex that connects $J$ segments, and that C1, C2, C3, C4 and C5 hold. Let $\hat{m}(x)=\boldsymbol{e}^{T}_{1}\left(\boldsymbol{X}^{T}\boldsymbol{W}\boldsymbol{X}\right)^{-1}\boldsymbol{X}^{T}\boldsymbol{W}\boldsymbol{Y}$, where the matrices are defined as above, then

[TABLE]

where the matrices are defined as follows:

$C(\omega,x)=m(x)-\omega m(x)^{2}$. 2. 2.

$\boldsymbol{U}_{0}$* is symmetric, and its first row is $\sum_{j}\mu^{(j)}_{0},\mu^{(1)}_{1},\dots,\mu^{(J)}_{1}$ and its diagonal is $\sum_{j}\mu^{(j)}_{0},\mu^{(1)}_{2},$ $\dots,\mu^{(J)}_{2}$, where $\mu^{(j)}_{i}=\int_{e_{j}}u^{i}K(u)du$. All other entries are zero.* 3. 3.

$\boldsymbol{U}_{1}$* is symmetric, and its first row is $0,-\mu^{(1)}_{0},\cdots,-\mu^{(l-1)}_{0},$ $\sum_{j\neq l}\mu^{(j)}_{0},-\mu^{(l+1)}_{0},\cdots,-\mu^{(J)}_{0}$. Its *$(l+1)$th row is $\sum_{j\neq l}\mu^{(j)}_{0},\mu^{(1)}_{1},\cdots,\mu^{(l-1)}_{1},0,\mu^{(l+1)}_{1},\cdots,\mu^{(J)}_{1}$ and its diagonal is $0,-2\mu^{(1)}_{1},$ $\cdots,-2\mu^{(l-1)}_{1},0,-2\mu^{(l+1)}_{1},\cdots,-2\mu^{(J)}_{1}$. All other entries of $\boldsymbol{U}_{1}$ are zero. 4. 4.

$\boldsymbol{R}_{0}$* is a column vector, its first entry is $\sum_{j\neq l}m^{(2)}_{e_{j}}(v)\mu^{(j)}_{2}+m^{(2)}_{e_{l}}(x)\mu^{(l)}_{2}$, and the rest are $m^{(2)}_{e_{j}}(v)\mu^{(j)}_{3}$, for $j=1,\cdots,J$.* 5. 5.

$\boldsymbol{R}_{1}$* is a column vector, its first entry is $-2\sum_{j\neq l}m^{(2)}_{e_{j}}(v)\mu^{(j)}_{1}$, and the rest, except the *$(l+1)$*th entry, are $-3m^{(2)}_{e_{j}}(v)\mu^{(j)}_{2}$, for $j\neq l+1$, and the *$(l+1)$th entry is $\sum_{j\neq l}m^{(2)}_{e_{j}}(x)\mu^{(j)}_{2}$, for $j=l+1$. 6. 6.

$\boldsymbol{M}_{0}$* is symmetric, its first row is $\sum_{j}R_{0}^{(j)},R_{1}^{(1)},\cdots,R_{1}^{(J)}$, and its diagonal entries are $\sum_{j}R_{0}^{(j)},R_{2}^{(1)},\cdots,R_{2}^{(J)}$. All the other entries of $\boldsymbol{M}_{0}$ are zero.* 7. 7.

$\boldsymbol{M}_{1}$* is symmetric, its first row is $0,-R^{(1)}_{0},\cdots,-R^{(l-1)}_{0},\sum_{j}R^{(j)}_{0},-R^{(l+1)}_{0},\cdots,-R^{(J)}_{0}$, its *$(l+1)$th row is

$\sum_{j}R^{(j)}_{0},R^{(1)}_{1},\cdots,R^{(l-1)}_{1},0,$* $R^{(l+1)}_{1},\cdots,R^{(J)}_{1}$, and its diagonal is $0,-R^{(1)}_{1},\cdots,$ $-R^{(l-1)}_{1},0,-R^{(l+1)}_{1},\cdots,-R^{(J)}_{1}$. All other entries of $\boldsymbol{M}_{1}$ are zero.*

Note that for $v$ to be in the $h$-neighborhood of the evaluation point $x$, we need $|v-x|<h$, and for $v$ to stay in the $h$-neighborhood of the evaluation point asymptotically, we need $x$ to approach $v$ at least as fast as $h$ approaches [math]. So $(v-x)/h$ is small and approaches zero as $N\to\infty$.

5.3 Estimation with equal first derivatives at the vertex

We have studied the asymptotic properties of local piecewise linear regression estimator at evaluation points that are within the $h$-neighborhood of a vertex, under the assumption that ${m}_{e_{j}}(v)$, $j=1,\dots,J$, agree at the vertex. Similarly, we can construct an asymptotic test for the joint equality of the limits of ${m}^{(1)}_{e_{j}}(x)$, $j=1,\dots,J$, as $x\to v$ along multiple edges. If we accept the null hypothesis that the first derivatives are equal at the vertex, this constraint should be added to the estimation procedure. For example, in Figure 2, consider an evaluation point $x\in e_{1}$ and a data point $x_{2}\in e_{2}$. Using local linear approximations, we have

[TABLE]

The regression function at $x_{3}\in e_{3}$ has a similar expansion. In the last line we used $m^{(1)}_{e_{2}}(v)=m^{(1)}_{e_{1}}(v)$, and expand $m^{(1)}_{e_{1}}(v)$ around the evaluation point $x$ on $e_{1}$. However to estimate first derivatives, one should use at least local quadratic polynomials.

6 Practical Issues

6.1 Estimation with large bandwidth

In this section, we will describe a more complex situation. When there are loops in the network, we use the shortest path between the evaluation point and a data point, and the distance between the points is defined by the shortest path distance. When there are multiple vertices in the $h$-neighborhood of an evaluation point, and vertices have different results from the joint equality tests, we will show that only data points that have direct access to the evaluation point contribute to the estimation. If $x$ is the evaluation point, $x_{i}$ is a data point that is within the $h$-neighborhood of $x$, and $v_{1},\dots,v_{n}$ are the vertices on the shortest path between $x$ and $x_{i}$, we say that $x_{i}$ has direct access to $x$ if the H0 (Theorem 2) is accepted at all $v_{i}$, $i=1,\cdots,n$. We will demonstrate this via the following example.

Here $x$ is the evaluation point, and $x_{1}$, $x_{2}$ and $x_{3}$ are data points on edges $e_{1}$, $e_{2}$ and $e_{3}$ respectively. We consider two scenarios:

$H_{0}$ is accepted at $v_{1}$ and $v_{2}$, 2. 2.

$H_{0}$ is accepted at $v_{1}$ and rejected at $v_{2}$.

In the first scenario, by Taylor expansions, we have

[TABLE]

Now, since $H_{0}$ is accepted at $v_{1}$ and $v_{2}$, we have $m_{e_{1}}(v_{1})=m_{e_{2}}(v_{1})$ and $m_{e_{2}}(v_{2})=m_{e_{3}}(v_{2})$, and (19) and (20) become

[TABLE]

The problem of minimizing $\sum_{i,j}\{y_{ij}-m_{e_{j}}(x_{i})\}^{2}K_{ij}$ can be set up in the same way as before. However, in the second scenario, since we do not have $m_{e_{2}}(v_{2})=m_{e_{3}}(v_{2})$, the problem becomes minimization of $\sum_{i,j=1,2}\{y_{ij}-m_{e_{j}}(x_{i})\}^{2}K_{ij}+\sum_{i}\{y_{i3}-m_{e_{3}}(x_{i})\}^{2}K_{i3}$. We see that the estimation of $m_{e_{1}}(x)$ only depends on the first summation. In other words, only data points that have direct access to the evaluation point contribute to the estimation of $m(x)$.

6.2 Bin width

For fixed sample size, we will let the bin width $\omega\to 0$. For simplicity, we show how this affects the calculation of $\left(\omega\boldsymbol{X}^{T}\boldsymbol{W}\boldsymbol{X}\right)^{-1}\omega\boldsymbol{X}^{T}\boldsymbol{W}\boldsymbol{Y}$. In the case where the regression function is estimated on each edge individually, when $\omega\to 0$, to calculate $\omega\boldsymbol{X}^{T}\boldsymbol{W}\boldsymbol{X}$ for $x$ near the vertex, we use that the Riemann sum converges to the integral as $\omega\to 0$:

[TABLE]

where $x=ch$, $0\leq c\leq 1$, and $\mu^{p}(c)=\int^{1}_{-c}u^{p}K(u)du$ is the $p$th truncated moment. Note that when the evaluation point is the vertex itself, $c=0$. The calculation of the term $\omega\boldsymbol{X}^{T}\boldsymbol{W}\boldsymbol{Y}$ requires

[TABLE]

Here $y$ is bar height and $\omega$ is bar width, so $y\,\omega$ is bar area, which equals to $c/N$, where $c$ is bar count and $N$ is sample size. Let $z_{i}$ be the location of $i$th observation. As $\omega\to 0$, eventually each bin is either empty or has exactly one observation in it. Therefore bin count $c_{i}$ is [math] for most bins, and $1$ for bins that contain one observation, and the location of that observation becomes the bin center in the limit.

When the joint equality constraint is added, to approximate $\omega\boldsymbol{X}^{T}\boldsymbol{W}\boldsymbol{X}$ for $x\in e_{l}$ near the vertex, we also need to calculate, for $j\neq l$, $x_{ij}\in e_{j}$,

[TABLE]

Binning with $\omega$ fixed might still be used with very large sample size.

7 Implementation

7.1 Simulation studies

Let us consider the simple network shown in Figure 2: three edges meet at a vertex. We propose three interesting cases on this network using the beta distribution. The three cases that we report here are

We simulate $500$ points on each edge, with the vertex being the origin, from $Beta(1,2)$, $Beta(1,3)$, and $Beta(1,4)$ respectively. 2. 2.

We simulate $500$ points on each edge, with the vertex being the origin, from $Beta(1,4)$. 3. 3.

For each edge, with the vertex being the origin, we simulate $500$ points from the truncated (from $0.5$ to $1$) $Beta(4,4)$. Then the points are shifted to the vertex by $0.5$ and finally multiplied by $2$, so the support of the true density on each edge is from the vertex to a point that is unit distance away from the vertex.

The true density function over the network is normalized so it integrates to $1$. Our simulation study focuses on demonstrating that the proposed estimator can accommodate various behaviors of the true regression function (the true density function), especially at the vertex, namely, discontinuous (Case I), continuous with discontinuous first derivative (Case II), and continuous with continuous first  derivative (Case III).

Simulation results for one dataset are shown in Figure 4. The result for local piecewise linear density is presented in the last row. The network is in red, and the true density functions over the network are the black dashed lines. The blue lines are our estimates. The first column is case I, where the joint equality test for the regression functions at the vertex is rejected, and our estimate is equivalent to local linear regression on each edge. The second column is case II. We first fitted a local linear regression model on each edge, and tested joint equality at the vertex. Since we fail to reject the null hypothesis, we re-estimated the regression function in the $h$-neighborhood of $v$ using data from all neighboring edges subject to the equality constraint. Turning to case III in the third column, the estimate is also subject to the joint equality constraint.

In case III, in addition to joint equality of the regression functions at the vertex, joint equality of their first derivatives is also established. For better viewing and comparison between different amounts of smoothness at the vertex, Figure 5 zooms in at the vertex. Estimation in left panel of Figure 5 is subject to only one constraint (joint equality of the regression functions at the vertex). However, after testing, by adding the constraint that the first derivatives are jointly equal, we see in the right panel of Figure 5 that the estimated curve is smoother over the vertex. This is desired because the true density function is continuous and has continuous first derivative over the vertex.

We should note that any sufficiently small bin width will do, but selection of the bandwidth takes more care. We note that cross-validation could be used for bandwidth selection, and more computationally efficient bandwidth selector for polynomial regression on networks is under investigation.

7.2 Comparison with existing methods

Let us consider again the three cases described in Section 7.1 and apply the equal-split discontinuous kernel estimator (ESDK) [9], the equal-split continuous kernel estimator (ESCK) [8], and the diffusion estimator (DE) [7] to each case.

ESDK [9] modifies (1), so it conserves mass. The algorithm makes a copy of the kernel function for each evaluation point, confined to the line edge containing that point. At each fork, the remaining tail of the kernel is split equally between the outgoing edges. Suppose there are $J-1$ outgoing edges, each outgoing edge receives a copy of the kernel tail weighted by $1/(J-1)$. ESCK [8] uses another modified version of (1). At each fork, with $J-1$ outgoing edges, each outgoing edge receives a copy of the kernel weighted by $2/J$, while the incoming segment receives a copy with the negative weight $2/J-1$. DE [7] uses the estimator $\hat{f}(u)=\sum^{N}_{i=1}K_{t}(u|x_{i})/N$, where $K_{t}(u|z)$ is the heat kernel on a network [3]. Intuitively, $K_{t}(u|z)du$ is the probability that a Brownian motion on the network, started at location $z\in L_{T}$ at time [math], will fall in the infinitesimal interval of length $du$ around the point $u$ at time $t$. The bandwidth parameters for ESDK, ESCK, and DE are selected by cross validation.

Figure 4 illustrates the performances of ESDK, ESCK and DE, compared to LPLR. Columns one, two and three correspond to cases I, II and III, respectively. For case I (discontinuous density), both ESCK and DE produce a continuous estimate at the vertex. ESDK, however, produces a discontinuous estimate, but the bias is still significant. Although ESDK is not a standard kernel density estimator, one should note that the asymptotic bias of kernel density estimation has order of $h$ on the boundary, whereas local linear estimation has order of $h^{2}$. The significant bias at the vertex can be also due to the fact that ESDK uses data from all neighboring edges as the evaluation point approaches the vertex, so it overestimates low density edges and underestimate high density edges.

For case II (continuous density with discontinuous first derivative), ESDK produces a discontinuous estimate. Although ESCK and DE produce continuous estimates, their bias are considerably higher compared to local piecewise polynomial estimation. This is due to lower order of the asymptotic bias of the kernel estimator near a vertex.

For case III (the density and its first derivative are both continuous) the vertex behaves like an interior point. While ESDK still produces a discontinuous estimate, ESCK and DE are comparable to local piecewise linear estimation, because the asymptotic bias of kernel density estimation and local linear estimation are both of order $h^{2}$ at an interior point. However, ESCK and DE do not estimate the derivatives of the estimated curves as the evaluation point approaching a vertex from multiple directions. Hence, these methods are inadequate when the derivatives are of interest.

Finally we report the bias, standard deviation, and mean squared error of the proposed local (piecewise) polynomial regression compared to all existing methods. In each case (I, II and III), the result is based on simulation of $100$ data sets, and each data set consists of $1000$ data points on each edge. Bias, standard deviation and mean squared error are reported at the vertex as it is approached from $e_{2}$ (see Figure 2). The simulation result is summarised in Table 1. We see that ESDK, ESCK and DE only produce comparable results when the vertex behaves like an interior point, whereas LPLR is superior in all other cases. Simulations are run on a Macbook Pro Mid 2015 with 2.5 GHz Intel Core i7. Average time (Case I, II and III) for LPLR, ESDK, ESCK and DE are 6.0082 secs, 4.3217 hrs, 6.4872 hrs and 6.6266 secs for all $1000$ datasets. ESDK, ESCK and DE are implemented in the R package spatstat. We note that the proposed local linear estimator requires at least moderately large data sets in order to produce good models. But for a fixed network and a fixed number of evaluation points, since the optimal bandwidth is of order $n^{-1/5}$, the computation time only scales sub-linearly with sample size. More generally, the computation time is $O(n^{4/5})$ per evaluation point. If the number of evaluation points is $O(n^{1/5})$, then the computation time grows like $n$.

7.3 Additional bias

The proposed method is subject to additional error because of the limitations of the test described in Section 4. For a large network, we would have to consider a large set of statistical tests simultaneously, and multiple testing problem would likely occur. It means that among the vertices at which the true density is continuous, the test rejected continuity at least one of them. However the main focus of this paper is the estimation of density function at individual locations, and the bias introduced by mistaking a continuous vertex as a discontinuous vertex is just the boundary bias of a local polynomial regression estimator. If one wants to count the number of vertices at which the density function is discontinuous, numerous correction techniques have been proposed in the literature, such as the Holm-Bonferroni method [holm1979simple]. Hence we will focus on type II error.

Our estimator will also have an additional bias when the test that the density is continuous at a vertex makes a type II error, that is, when the density is not continuous but the test accepts that it is continuous. The size of additional bias introduced by a type II error is less clear. We investigate this problem through simulations. For simplicity, let the origin be a vertex, and let the two unit-length edges meet at the vertex. On the left edge, we first sample 1000 points from Beta$(1,\beta_{l})$, then multiply the sample by $-1$, so the data points on the left edge range from $-1$ to [math]. On the right edge, we draw 1000 point from Beta$(1,\beta_{r})$. Pairs of $(\beta_{l},\beta_{r})$ are given in the first column of Table 2. Type II error rate is the probability that the test (Theorem 2) accepts the null hypothesis that the true density is continuous at a vertex when it is not. We perform $3000$ simulations of $1000$ data points per edge for each pair of $(\beta_{l},\beta_{r})$. We report the bias, standard deviation and mean squared error at the vertex approached from the right edge. The result is summarized in Table 2. We see that when type II error rates are large, LPLR under the continuity constraint produces little additional bias due to small gaps between $\beta_{l}$ and $\beta_{r}$. Similarly, when $\beta_{l}$ and $\beta_{r}$ are far apart, there is also little additional bias due to small type II error rates. On the other hand, when type II error rate is in the mid range, for example from $0.6$ to $0.7$, there is considerable amount of additional bias at the vertex. Hence in the bias column, we should see bigger (in magnitude) biases in the middle rows, and as we approach the top and bottom rows, biases become smaller (in magnitude). Finally we note that the top row has bigger bias than the bottom row, despite that the type II error rate is practically zero. This is due to that local linear regression generally has bigger bias for functions with greater curvature (Beta$(1,4.5)$ compared to Beta$(1,4)$).

7.4 Application to real data

In this section we apply the proposed local piecewise polynomial regression to dendrite data. The data was collected by the Kosik Lab, UC Santa Barbara, and first analyzed by [1] and [6]. Dendrites are branching filaments which extend from the main body of a neuron (nerve cell) to propagate electrochemical signals. Spines are small protrusions on the dendrites. The network shown in Figure 1 is one of the ten dendritic trees of this neuron. A dendritic tree consists of all dendrites issuing from a single root branching off the cell body; each neuron typically has 4 to 10 dendritic trees. This example was chosen because it is large enough to demonstrate our techniques clearly, without being too large for graphical purposes. The events on the network are the locations of 566 spines observed on one branch of the dendritic tree of a rat neuron. We will show the result of applying the existing methods and the proposed local (piecewise) polynomial estimator to the dendrite data of Figure 1. We show density estimates at vertices A, B, C and D. We used a fixed bandwidth of $9$ microns (the network has a total length of $1933.653$ microns), and we will leave bandwidth selection and variable bandwidths to future projects. Unlike ESDK, ESCK and DE, the continuity of the LPLR estimates at a vertex is only imposed when there is evidence that the density is continuous there. Since the proposed method is particularly advantageous near vertices, we will mostly focus on that region.

Figure 6 shows density estimates near vertex A, using ESCK, ESDK and DE. We used a fixed bandwidth of $9$ microns, and an Epanechnikov kernel in ESCK and SEDK. Both ESCK and DE produced continuous estimates at the vertex. Although ESCK produced discontinuity at the vertex, the estimation clearly used data from all edges. Consequently the estimated density on the background edge near the vertex is overly inflated. In Figure 1, the data suggests that the density function has different behaviours when approaching vertex A from different directions. However none of the exisiting methods can characterise this. Next we apply the proposed piecewise local linear regression approach, and the result is in Figure 7. Since the data does not support continuity at vertex A, no further constraint estimation is required.

Next we look at vertex B in Figure 1. Unlike vertex A, spines around vertex B seem to become more sparse as approaching the vertex from all directions, especially from the two edges on the left. Apply ESCK, ESDK and DE, we get the estimates in Figure 8. ESCK and DE produced continuous estimates at vertex B. Under these two approaches, the estimated density is continuous at the vertex over any pair of edges, i.e. $\hat{f}(x)$ is continuous as $x$ travels from $e_{i}$ to $e_{j}$, where $e_{i}$ and $e_{j}$ are two edges connected by the vertex. However continuity may only exist over a subset of edges, and ESCK, ESDK and DE are unable to provide such flexibility. By the proposed piecewise local linear estimation, we first estimate the density on each edge individually using only the data on that edge. This is panel (a) in Figure 9. By the testing procedure discussed in Section 4, there’s strong evidence in the data that the densities from the left two edges are equal when approaching the vertex. Following the constrained estimation in Section 5, we re-estimate the densities on the two edges by local piecewise estimator using data from all three edges subject to the constraint that the densities on the left two edges are equal at the vertex. The result is in panel (b). Suppose $e_{j}$ for $j=1,\cdots,J$ are connected by vertex $v$, let $f_{e_{j}}(x)$ be the density over $e_{j}$. Using the proposed piecewise local linear estimation, we can test for the equality of $f_{e_{j}}(x)$ for any subset of $j\in\{1,\cdots,J\}$,as $x\to v$, and achieve the desired continuity as suggested by the data. On the other hand, under ESCK, ESDK and DE, the continuity at a vertex is decided by the user prior to model fitting, and the decision making is not data-driven.

Similarly, for vertex C, under piecewise local linear estimation, there’s evidence in the data that densities from the two edges in the foreground are equal when approaching the vertex. Figure 11 shows the estimated densities near vertex C, by ESCK, ESDK and DE, while piecewise local linear estimation is in Figure 11.

For vertex D, under piecewise local linear estimation, there’s evidence in the data that the densities from the two left edges are equal when approaching the vertex. Figure 13 shows the estimated densities near vertex D, by ESCK, ESDK and DE, while piecewise local linear estimation is in Figure 13. From this example we see that the network is embedded in $\mathbb{R}^{3}$. The vertex is on the right of the plot where labeled “vertex D”. There is an overpass on the left of the plot, where the two edges do not intersect.

8 Discussion

As we mentioned in the introduction, there is great potential demand in many fields for estimating the density of events on a network. An easy method for such estimations is to use the ordinary kernel density estimation method that assumes an unbounded plane, or kernel density estimation on the real line with the Euclidean distance replaced by network distance. Many papers in the literature employ this method. However, this method yields a bias in density estimation and so the method is likely to lead to misleading conclusions. We also discussed the equal-split discontinuous kernel estimator, the equal-split continuous kernel estimator, and diffusion estimator. None of those methods allows for discontinuity in the estimates. The first two methods lack theoretical justification and are computationally expensive. The diffusion estimator is mathematically equivalent to an infinite-sum generalization of the equal-split continuous rule applied to the Gaussian density, and it inherits the asymptotic properties of a kernel density estimator [7]. Also the diffusion estimator has a slower rate for the boundary bias.

In this paper, we have formulated a density estimation procedure on a linear network via local piecewise polynomial regression by way of binning. We first apply local polynomial regression on each edge individually, then we test joint equality of the regression functions at the vertex. If the null hypothesis is not rejected, locations within the $h$-neighborhood of the vertex are re-estimated by local piecewise polynomial regression using data from all neighboring edges, subject to the equality constraint. The proposed piecewise linear procedure only imposes continuity at vertices when there is evidence in the data, and its asymptotic bias has the same rate at a vertex as an interior point. We studied the local linear case in detail, while there is a straightforward extension to higher-order polynomial approximation. When applying the proposed method to real data, if there are loops in the network, for simplicity, we only considered the shortest path between points, and we showed that only data points that have direct access to the evaluation point contribute to the estimation of the regression function.

Due to space limitation, we have only considered only fixed-bandwidth smoothing, and we have not considered data-based bandwidth selection nor adaptive smoothing on a network. We proposed a test of equal intercepts at a vertex to decide whether to assume equal intercepts when estimating the density in the neighborhood of the vertex. In the future, we will consider an estimator that, rather than using a pre-test, shrinks the unequal-intercepts estimator towards the equal-intercepts estimator. We have assumed that all locations are measured without error and lie exactly on the linear network. However this is not true for some applications, such as ambulance or taxi, where there are GPS error in their locations. Further study into such measurement error problems is required. In addition to estimating probability density or point process intensity, the proposed procedure is also used for regression problems, such as varying coefficient models on a linear network.

9 Proofs

Proof of Theorem 1.

By Taylor expansion,

[TABLE]

Note that if $m$ is a linear function then $m^{(r)}(x)=0$ for $r\geq 2$ so that the local linear estimator is exactly unbiased when $m$ is a linear function. To find the leading bias term for general function $m$, note that

[TABLE]

where $\sigma^{i}_{K}=\int u^{i}K(u)du$. Some straightforward matrix algebra then leads to the following expression for the leading bias term

[TABLE]

To derive the asymptotic variance of $\hat{m}(x)$ we have

[TABLE]

where $\boldsymbol{V}=diag(f(x_{1}),\dots,f(x_{n}))$, with $f(x_{i})=(1/N\omega)m(x_{i})-(1/N)m(x_{i})^{2}$. We stress that $f(x_{i})$ depends on $N$ and $\omega$, but we drop them from its notation for simplicity. Now use approximation analogous to those used above we have

[TABLE]

where $R^{i}_{K}=\int u^{i}K(u)^{2}du$. These expressions can be combined to obtain

[TABLE]

where $C(N,h,\omega,x)=f(x)\omega/h$, and

[TABLE]

∎

Proof of Theorem 2.

Recall the asymptotic result of local polynomial regression:

[TABLE]

where $B_{e_{l}}$ and $V_{e_{l}}$ are the asymptotic bias and variance of $\hat{m}_{e_{l}}(x)$. Let $\hat{\boldsymbol{m}}=(\hat{m}_{e_{1}}(v),\dots,\hat{m}_{e_{J}}(v))^{T}$, we have asymptotically, $h\to 0$ and $nh\to\infty$,

[TABLE]

Now consider a $(J-1)\times J$ contrast matrix $\boldsymbol{C}$ such that $\boldsymbol{C1}=\boldsymbol{0}$ (i.e. each row sum to zero). One choice of $\boldsymbol{C}$ is

[TABLE]

which simply contrasts edge 1 with edge 2, edge 1 with edge 3. Then

[TABLE]

where

[TABLE]

Under the null hypothesis,

[TABLE]

We want to calculate the probability of generating a point at least as unlikely as the observed data point. To do that, we note that under the null hypothesis, the Mahalanobis distance follows a chi-square distribution

[TABLE]

The test statistic is invariant under the choice of contrast matrices. This is easy to see by noticing that the rows of $\boldsymbol{C}$ are linearly independent. So we have a basis of some vector space $\mathbf{V}$ (and it doesn’t matter if $\mathbf{V}$ is all of $\mathbb{R}^{p+1}$, or some subspace thereof), and two different ordered bases for $\mathbf{V}$, $\boldsymbol{b}_{1}$ and $\boldsymbol{b}_{2}$ (necessarily of the same size, since two bases of the same vector space always have the same size, here they are the transpose of two contrast matrices $\boldsymbol{C}_{1}$ and $\boldsymbol{C}_{2}$):

[TABLE]

A change-of-basis matrix is a matrix that translates from $\boldsymbol{b}_{1}$ coordinates to $\boldsymbol{b}_{2}$ coordinates. That is, $A$ is a change-of-basis matrix (from $\boldsymbol{b}_{1}$ to $\boldsymbol{b}_{2}$) if, given the coordinate vector $[\mathbf{x}]_{\boldsymbol{b}_{1}}$ of a vector $\mathbf{x}$ relative to $\boldsymbol{b}_{1}$, then $A[\mathbf{x}]_{\boldsymbol{b}_{1}}=[\mathbf{x}]_{\boldsymbol{b}_{2}}$ gives the coordinate vector of $\mathbf{x}$ relative to $\boldsymbol{b}_{2}$, for all $\mathbf{x}$ in $\mathbf{V}$.

To get a change-of-basis matrix, we write each vector of $\boldsymbol{b}_{1}$ in terms of $\boldsymbol{b}_{2}$, and these are the columns of $A$, for $i=1,\dots,n$, $\mathbf{v}_{i}=a_{1i}\mathbf{w}_{1}+a_{2i}\mathbf{w}_{2}+\cdots+a_{ni}\mathbf{w}_{n}$. We know we can do this because $\boldsymbol{b}_{2}$ is a basis, so we can express any vector (in particular, the vectors in $\boldsymbol{b}_{1}$) as linear combinations of the vectors in $\boldsymbol{b}_{2}$. Then the change-of-basis matrix translating from $\boldsymbol{b}_{1}$ to $\boldsymbol{b}_{2}$ is

[TABLE]

Matrix $A$ is always invertible. This is because just like there is a change-of-basis from $\boldsymbol{b}_{1}$ to $\boldsymbol{b}_{2}$, there is also a change-of-basis from $\boldsymbol{b}_{2}$ to $\boldsymbol{b}_{1}$. Since $\boldsymbol{b}_{1}$ is a basis, we can express every vector in $\boldsymbol{b}_{2}$ using the vectors in $\boldsymbol{b}_{1}$, for $i=1,\dots,n$, $\mathbf{w}_{i}=b_{1i}\mathbf{v}_{1}+b_{2i}\mathbf{v}_{2}+\cdots+b_{ni}\mathbf{v}_{n}$. So the matrix $B$, with

[TABLE]

has the property that given any vector $\mathbf{x}$, if $[\mathbf{x}]_{\boldsymbol{b}_{2}}$ is the coordinate vector of $\mathbf{x}$ relative to $\boldsymbol{b}_{2}$, then $B[\mathbf{x}]_{\boldsymbol{b}_{2}}=[\mathbf{x}]_{\boldsymbol{b}_{1}}$ is the coordinate vector of $\mathbf{x}$ relative to $\boldsymbol{b}_{1}$. Then applying first $A$ and then $B$ translates $\boldsymbol{b}_{1}$ coordinates into $\boldsymbol{b}_{2}$ coordinates and back to $\boldsymbol{b}_{1}$ coordinates, and thus $AB$ must be the identity matrix (likewise for $BA$). So $A$ and $B$ are both invertible, and every change-of-basis matrix is necessarily invertible. ∎

Proof of Theorem 3.

Standard calculation shows that the leading term of the bias is given by

[TABLE]

where $\boldsymbol{S}$ is a column vector with entries being $(v-x)^{2}m^{(2)}_{e_{l}}(x)+(x_{ij}-v)^{2}m^{(2)}_{e_{i}}(v)$, for $i=1,\dots,n_{j}$, and $j=1,\dots,J$. Note that for data points $x_{ij}\in e_{l}$, where $e_{l}$ is the edge that the evaluation point $x$ is located, the entries are simplified to $(x_{il}-x)^{2}m^{(2)}_{e_{l}}(x)$, for $i=1,\dots,n_{l}$. To approximate $\omega^{-1}\left(\boldsymbol{X}^{T}\boldsymbol{W}\boldsymbol{X}\right)^{-1}$, we note that

[TABLE]

where $\mu^{(j)}_{i}=\int_{e_{j}}u^{i}K(u)du$. If follows that

[TABLE]

where $\boldsymbol{A}=diag(1,h,\dots,h)$. The matrix $\boldsymbol{U}_{0}$ is symmetric, and its first row is $\sum_{j}\mu^{(j)}_{0},\mu^{(1)}_{1},\dots,\mu^{(J)}_{1}$ and its diagonal is $\sum_{j}\mu^{(j)}_{0},\mu^{(1)}_{2},\dots,\mu^{(J)}_{2}$. All other entries are zero. We also have that matrix $\boldsymbol{U}_{1}$ is symmetric, and its first row is $0,-\mu^{(1)}_{0},\cdots,-\mu^{(l-1)}_{0},\sum_{j\neq l}\mu^{(j)}_{0},-\mu^{(l+1)}_{0},\cdots,-\mu^{(J)}_{0}$. Its $(l+1)$th row is $\sum_{j\neq l}\mu^{(j)}_{0},\mu^{(1)}_{1},\cdots,\mu^{(l-1)}_{1},0,\mu^{(l+1)}_{1},\cdots,\mu^{(J)}_{1}$ and its diagonal is $0,-2\mu^{(1)}_{1},\cdots,-2\mu^{(l-1)}_{1},$ $0,-2\mu^{(l+1)}_{1},\cdots,-2\mu^{(J)}_{1}$. All other entries of $\boldsymbol{U}_{1}$ are zero. Finally $\boldsymbol{U}_{2}=diag(0,\mu^{(1)}_{0},\dots,\mu^{(J)}_{0})$. Take inverse we have

[TABLE]

To approximate $\omega\boldsymbol{X}^{T}\boldsymbol{W}\boldsymbol{S}$, we note that

[TABLE]

Combining this result with the approximations before, we get

[TABLE]

where $\boldsymbol{R}_{0}$, $\boldsymbol{R}_{1}$, $\boldsymbol{R}_{2}$ and $\boldsymbol{R}_{3}$ are column vectors. The first entry of $\boldsymbol{R}_{0}$ is $\sum_{j\neq l}m^{(2)}_{e_{j}}(v)\mu^{(j)}_{2}+m^{(2)}_{e_{l}}(x)\mu^{(l)}_{2}$, and the rest are $m^{(2)}_{e_{j}}(v)\mu^{(j)}_{3}$, for $j=1,\cdots,J$. The first entry of $\boldsymbol{R}_{1}$ is $-2\sum_{j\neq l}m^{(2)}_{e_{j}}(v)\mu^{(j)}_{1}$, and the rest, except the $(l+1)$th entry, are $-3m^{(2)}_{e_{j}}(v)\mu^{(j)}_{2}$, for $j\neq l+1$, and the $(l+1)$th entry is $\sum_{j\neq l}m^{(2)}_{e_{j}}(x)\mu^{(j)}_{2}$. The first entry of $\boldsymbol{R}_{2}$ is $m^{(2)}_{e_{l}}(x)\sum_{j\neq l}\mu^{(j)}_{0}+\sum_{j\neq l}m^{(2)}_{e_{j}}(x)\mu^{(j)}_{0}$, and the rest, except the $(l+1)$th entry, are $m^{(2)}_{e_{l}}(x)\mu^{(j)}_{1}+3m^{(2)}_{e_{j}}(v)\mu^{(j)}_{1}$, for $j\neq l+1$, and the $(l+1)$th entry is $-2\sum_{j\neq l}m^{(2)}_{e_{j}}(v)\mu^{(j)}_{1}$. Finally, the first entry of $\boldsymbol{R}_{3}$ is [math], and the rest, except the $(l+1)$th entry, are $-m^{(2)}_{e_{l}}(x)\mu^{(j)}_{0}-m^{(2)}_{e_{j}}(v)\mu^{(j)}_{0}$, for $j\neq l+1$, and the $(l+1)$th entry is $m^{(2)}_{e_{l}}(x)\sum_{j\neq l}\mu^{(j)}_{0}+\sum_{j\neq l}m^{(2)}_{e_{j}}(v)\mu^{(j)}_{0}$. Consequently,

[TABLE]

To derive the asymptotic variance of $\hat{m}(x)$ we have

[TABLE]

where $\boldsymbol{V}=diag(f(x_{1}),\dots,f(x_{n}))$, with $f(x_{i})=(1/N\omega)m(x_{i})-(1/N)m(x_{i})^{2}$. We stress that $f(x_{i})$ depends on $N$ and $\omega$, but we drop them from its notation for simplicity. Now use approximation analogous to those used above we have

[TABLE]

It follows that

[TABLE]

where $\boldsymbol{M}_{0}$ is symmetric, its first row is $\sum_{j}R_{0}^{(j)},R_{1}^{(1)},\cdots,R_{1}^{(J)}$, and its diagonal entries are $\sum_{j}R_{0}^{(j)},R_{2}^{(1)},\cdots,R_{2}^{(J)}$. All the other entries of $\boldsymbol{M}_{0}$ are zero. $\boldsymbol{M}_{1}$ is symmetric, its first row is $0,-R^{(1)}_{0},\cdots,-R^{(l-1)}_{0},\sum_{j}R^{(j)}_{0},-R^{(l+1)}_{0},\cdots,-R^{(J)}_{0}$, its $(l+1)$th row is $\sum_{j}R^{(j)}_{0},R^{(1)}_{1},\cdots,R^{(l-1)}_{1},0,$ $R^{(l+1)}_{1},\cdots,R^{(J)}_{1}$, and its diagonal is $0,-R^{(1)}_{1},\cdots,-R^{(l-1)}_{1},0,-R^{(l+1)}_{1},\cdots,-R^{(J)}_{1}$. All other entries of $\boldsymbol{M}_{1}$ are zero. Finally $\boldsymbol{M}_{2}$ is also symmetric, its first row is the zero vector, its $(l+1)$th row is $0,-R^{(0)}_{1},\cdots,-R^{(l-1)}_{1},\sum_{j}R^{(j)}_{0},-R^{(l+1)}_{1},\cdots,-R^{(J)}_{1}$, and its diagonal is $0,R^{(1)}_{0},\cdots,R^{(l-1)}_{0},\sum_{j}R^{(j)}_{0},$ $R^{(l+1)}_{0},\cdots,R^{(J)}_{0}$. All other entries of $\boldsymbol{M}_{2}$ are zero.

Combining the above expressions, we have

[TABLE]

∎

Supplementary Materials

The supplementary materials include an R program containing code to perform the local linear regression method for density estimation on a network as described in this article. The program also contains all codes for simulating datasets used in the article.