Regression analysis of mixed sparse synchronous and asynchronous longitudinal covariates with varying-coefficient models

To estimate $\beta(t)$ and $\gamma(t)$ simultaneously, we consider a local linear estimator (Fan and Gijbels, 1996; Cao et al., 2018):

$$\begin{split}&\{\hat{\beta}_{w}(t),\hat{\dot{\beta}}_{w}(t),\hat{\gamma}_{w}(t),\hat{\dot{\gamma}}_{w}(t)\}\\ =&\mathop{\rm arg\min}_{\{\beta_{0},\beta_{1},\gamma_{0},\gamma_{1}\}}\sum_{i=1}^{n}\iint K_{h_{1},h_{2}}(t_{1}-t,t_{2}-t)\left\{Y_{i}(t_{1})-X_{i}(t_{1})^{T}\beta_{0}-X_{i}(t_{1})^{T}\beta_{1}(t_{1}-t)\right.\\ &\left.\ \ \quad\quad\quad\quad\quad\quad\quad-Z_{i}(t_{2})^{T}\gamma_{0}-Z_{i}(t_{2})^{T}\gamma_{1}(t_{2}-t)\right\}^{2}dN^{*}_{i}(t_{1},t_{2}).\end{split}$$

(2.2)

For any fixed time point $t,$ denote $R(t_{1},t_{2},t)=\{X(t_{1})^{T},X(t_{1})^{T}(t_{1}-t),Z(t_{2})^{T},Z(t_{2})^{T}(t_{2}-t)\}^{T}$, and $\rho_{0}(t)\triangleq\{\beta_{0}(t)^{T},\dot{\beta}_{0}(t)^{T},\gamma_{0}(t)^{T},\dot{\gamma}_{0}(t)^{T}\}^{T},$ where $\beta_{0}(t)$ and $\gamma_{0}(t)$ are the true non-parametric regression functions and $\dot{\beta}_{0}(t)$ and $\dot{\gamma}_{0}(t)$ are the corresponding derivative functions. Assume that the conditional variance function $\mbox{var}\left\{Y(t)\mid X(t),Z(t)\right\}=\sigma\left\{t,X(t),Z(t)\right\}^{2}<\infty.$ We need the following assumptions.

1. (A1)

   $N^{*}_{i}(t,s)$ is independent of $\{X_{i}(t),Z_{i}(t),Y_{i}(t)\}$ and $E\left\{dN^{*}_{i}(t,s)\right\}=\eta(t,s)dtds,$ $i=1,\ldots,n,$ where $\eta(t,s)$ is a twice continuously differentiable function for any $(t,s)\in[0,1]^{\otimes 2}.$ For $t_{1}\neq s_{1},t_{2}\neq s_{2},{P}(dN^{*}(t_{1},t_{2})=1\mid N^{*}(s_{1},s_{2})-N^{*}(s_{1}-,s_{2}-)=1)=f(t_{1},t_{2},s_{1},s_{2})$ $dt_{1}dt_{2}$ where $f(t_{1},t_{2},s_{1},s_{2})$ is continuous for $t_{1}\neq s_{1},t_{2}\neq s_{2}$ and
   $f(t_{1}\pm,t_{2}\pm,s_{1}\pm,s_{2}\pm)$ exist.

2. (A2)

   $\rho_{0}(t)$ is twice continuously differentiable for $\forall t\in[0,1].$ For any fixed time point $t,$ $E\left\{R(t_{1},t_{2},t)R(s_{1},s_{2},t)^{T}\right\}\in{\mathbb{R}}^{2(p+q)\times 2(p+q)}$ and ${E}\left\{R(t_{1},t_{2},t)\right\}\in{\mathbb{R}}^{2(p+q)}$ are twice continuously differentiable for $(t_{1},t_{2},s_{1},s_{2})\in[0,1]^{\otimes 4}.$ $E\left\{C(t,t)C(t,t)^{T}\right\}\eta(t,t)\in{\mathbb{R}}^{(p+q)\times(p+q)}$ is positive definite for $\forall t\in[0,1],$ where $C(t,t)=\left\{X(t)^{T},Z(t)^{T}\right\}^{T}$. Moreover,

   $$\|E\left[C(t,t)C(t,t)^{T}\sigma\left\{t,X(t),Z(t)\right\}^{2}\right]\|_{\infty}\eta(t,t)<\infty,\quad\forall t\in[0,1],$$

   where for a square matrix $A,$ $\|A\|_{\infty}=\max_{1\leq i\leq n}\sum_{j=1}^{n}|a_{ij}|.$

3. (A3)

   $K(\cdot,\cdot)$ is a bivariate density function satisfying $\iint z_{1}K(z_{1},z_{2})dz_{1}dz_{2}=\iint z_{2}K(z_{1},z_{2})$
   $dz_{1}dz_{2}=0,\iint K(z_{1},z_{2})^{2}dz_{1}dz_{2}<\infty,\iint z_{1}^{2}K(z_{1},z_{2})dz_{1}dz_{2}<\infty,\iint z_{2}^{2}K(z_{1},z_{2})dz_{1}dz_{2}<\infty$ and $\iint z_{1}^{2}z_{2}^{2}K(z_{1},z_{2})dz_{1}dz_{2}<\infty.$

4. (A4)

   $(nh_{1}h_{2})^{1/2}(h_{1}^{2}+h_{2}^{2})\to 0$ and $nh_{1}h_{2}\to\infty.$

Condition (A1) requires that the observation process be independent of the covariate processes. Analogous assumptions have been made in Lin and Ying (2001) and Cao et al. (2015). Condition (A2) ensures the identifiability of $\rho_{0}(t)$ at any fixed time point $t,$ and some smoothness assumptions on the expectation of certain functional of $R(t_{1},t_{2},t)$ and $\rho_{0}(t).$ Conditions (A3) and (A4) specify valid kernels and bandwidths. For the bivariate kernel function, we recommend kernel functions with bounded support, such as the product of the univariate Epanechnikov kernel.

The following theorem establishes the asymptotic results of $\hat{\beta}_{w}(t)$ and $\hat{\gamma}_{w}(t).$ The proofs are relegated in the Appendix.

Theorem 1 derives the joint asymptotic distribution of $\hat{\beta}_{w}(t)$ and $\hat{\gamma}_{w}(t).$ If we let $h_{\alpha}=h_{1}=h_{2}$, the bias of $\hat{\beta}_{w}(t)$ and $\hat{\gamma}_{w}(t)$ in Theorem 1 is $O(h_{\alpha}^{2})$, which is the same as that in Theorem 2 of Cao et al. (2015). The reason why the $O(h_{\alpha})$ term vanishes is that we use symmetric density functions as the kernel functions, which are specified in (A3). When the sample size $n$ goes to infinity, under (A4), the bias terms vanish. In general, the off-diagonal part of the limiting variance covariance matrix is non-zero. In the special case that $E\{X(t)Z(t)^{T}\}=0_{p\times q},$ $\hat{\beta}_{w}(t)$ and $\hat{\gamma}_{w}(t)$ are independent. Moreover, if we assume that $\sigma\{t,X(t),Z(t)\}^{2}=\sigma(t)^{2},$ then $\mbox{Var}\{\hat{\beta}_{w}(t)\}=[E\{X(t)X(t)^{T}\}\eta(t,t)]^{-1}\big{\{}\iint K(z_{1},z_{2})^{2}dz_{1}dz_{2}\big{\}}\sigma(t)^{2}$ and $\mbox{Var}\{\hat{\gamma}_{w}(t)\}=[E\{Z(t)Z(t)^{T}\}\eta(t,t)]^{-1}\big{\{}\iint K(z_{1},z_{2})^{2}dz_{1}dz_{2}\big{\}}\sigma(t)^{2},$ respectively.

Our method depends on the proper choice of bandwidths. If we let $h_{\alpha}=h_{1}=h_{2},$ according to condition (A4), a valid bandwidth is larger than $O(n^{-1/2}),$ and the bias of $\hat{\beta}_{w}(t)$ and $\hat{\gamma}_{w}(t)$ in Theorem 1 is of order $h_{\alpha}^{2},$ so we should choose bandwidth $h_{\alpha}=o(n^{-1/6}).$ With this choice of $h_{\alpha}$, we obtain a rate of convergence $O(n^{1/3}),$ which is the same as that achieved in Cao et al. (2015). This is slower than the non-parametric rate of convergence $n^{2/5}$ due to the asynchronous longitudinal data structure.

In Theorem 1, $A^{*}(t)$ and $\Sigma^{*}(t)$ are not directly computable. In practice, we use the estimating equation derived from (2.2) and estimate the variance covariance matrix of $\hat{\rho}_{w}(t)$ with the sandwich formula:

$$\begin{split}&\widehat{\mathrm{var}}\{\hat{\rho}_{w}(t)\}\\ =&\left\{\sum_{i=1}^{n}\iint K_{h_{1},h_{2}}(t_{1}-t,t_{2}-t)R_{i}(t_{1},t_{2},t)R_{i}(t_{1},t_{2},t)^{T}dN^{*}_{i}(t_{1},t_{2})\right\}^{-1}\\ &\times\sum_{i=1}^{n}\left[\iint K_{h_{1},h_{2}}(t_{1}-t,t_{2}-t)R_{i}(t_{1},t_{2},t)\left\{Y_{i}(t_{1})-R_{i}(t_{1},t_{2},t)^{T}\hat{\rho}_{w}(t)\right\}dN^{*}_{i}(t_{1},t_{2})\right]^{\otimes 2}\\ &\times\left\{\sum_{i=1}^{n}\iint K_{h_{1},h_{2}}(t_{1}-t,t_{2}-t)R_{i}(t_{1},t_{2},t)R_{i}(t_{1},t_{2},t)^{T}dN^{*}_{i}(t_{1},t_{2})\right\}^{-1}.\end{split}$$

The variance covariance matrices of $\hat{\beta}_{w}(t)$ and $\hat{\gamma}_{w}(t)$ are the upper $p\times p$ submatrix of the upper left $(2p)\times(2p)$ submatrix and the upper $q\times q$ submatrix of the lower right $(2q)\times(2q)$ submatrix of $\widehat{\mathrm{var}}\{\hat{\rho}_{w}(t)\}$, respectively.

To get uncertainty quantification and evaluate the overall magnitude of variation of the non-parametric regression functions, we construct simultaneous confidence bands (SCB). Here, we use $\gamma(t)$ to illustrate. Specifically, for a pre-specified confidence level $1-\alpha$ and a given smooth function of $q$ dimensions $a(t),$ we aim to find smooth random functions $L(t)$ and $U(t)$, such that

$$\mathbb{P}\{L(t)\leq a(t)^{T}\gamma(t)\leq U(t),\forall t\in[a,b]\}\rightarrow 1-\alpha$$

as number of subjects $n\rightarrow\infty$ for a pre-specified closed interval $[a,b].$ Theoretical results on SCB construction rely on extreme value theory with slow rate of convergence (Fan and Zhang, 2000; Cao et al., 2018). We shall adopt a wild bootstrap approach to resample the residuals (Wu, 1986; Liu, 1988). The wild bootstrap is computationally efficient as it multiplies a random disturbance on the residuals without repeated estimation of the non-parametric regression function. Specific steps are as follows.

Step 1. For each subject $i,i=1,2,\ldots,n,$ calculate the kernel-weighted residual

$$\hat{r}_{i}(t_{1},t_{2},t)=K_{h_{1},h_{2}}(t_{1}-t,t_{2}-t)\left\{Y_{i}(t_{1})-R_{i}(t_{1},t_{2},t)^{T}\hat{\rho}_{w}(t)\right\}.$$

Step 2. Obtain an approximation of $\hat{\gamma}_{w}(t)-\gamma_{0}(t)$ through

$$\hat{Q}_{n}(t)=\frac{1}{n}\sum_{i=1}^{n}\left\{f_{n}(t)\right\}^{-1}\iint Z_{i}(t_{2})\hat{r}_{i}(t_{1},t_{2},t)dN^{*}_{i}(t_{1},t_{2}),$$

where $f_{n}(t)=\frac{1}{n}\sum_{i=1}^{n}\iint K_{h_{1},h_{2}}(t_{1}-t,t_{2}-t)Z_{i}(t_{2})Z_{i}(t_{2})^{T}dN^{*}_{i}(t_{1},t_{2}).$

Step 3. Independently generate $\left\{u_{i},i=1,\dots,n\right\}$ from a distribution with mean $0$ and standard deviation $1.$. For example, we can use $N(0,1)$ or the Rademacher distribution which takes values $+1$ and $-1$ with equal probability. Construct the stochastic process

$$\hat{Q}_{n}^{u}(t)=\frac{1}{n}\sum_{i=1}^{n}u_{i}\cdot\left\{f_{n}(t)\right\}^{-1}\iint Z_{i}(t_{2})\hat{r}_{i}(t_{1},t_{2},t)dN^{*}_{i}(t_{1},t_{2}).$$

Step 4. Repeat Step 3 $B$ times to obtain $\{\sup\limits_{t}\mid a(t)^{T}\hat{Q}_{n}^{u(1)}(t)\mid,\dots,\sup\limits_{t}\mid a(t)^{T}\hat{Q}_{n}^{u(B)}(t)\mid\},$ and denote its $1-\alpha$ empirical percentile as $c_{\alpha}.$ The SCB for $a(t)^{T}\gamma(t)$ can be written as

$$\left(a(t)^{T}\hat{\gamma}_{w}(t)-c_{\alpha},a(t)^{T}\hat{\gamma}_{w}(t)+c_{\alpha}\right).$$

We recommend the Rademacher distribution as the multiplier in Step 3 through extensive simulation studies with different random variables. The details are omitted.

We remark that for a given $p$-dimensional smooth function $d(t),$ the construction of SCB for $d(t)^{T}\beta(t)$ is exactly the same except in Step 2, let $\hat{J}_{n}(t)=\frac{1}{n}\sum_{i=1}^{n}\left\{g_{n}(t)\right\}^{-1}\iint X_{i}(t_{1})\hat{r}_{i}(t_{1},t_{2},t)\newline dN^{*}_{i}(t_{1},t_{2}),$ where $g_{n}(t)=\frac{1}{n}\sum_{i=1}^{n}\iint K_{h_{1},h_{2}}(t_{1}-t,t_{2}-t)X_{i}(t_{1})X_{i}(t_{1})^{T}dN^{*}_{i}(t_{1},t_{2}).$ The remaining steps are based on $\hat{J}_{n}(t)$ and we replace $a(t)$ by $d(t).$

In (2.1), $\beta(t)$ is the regression function of the synchronous longitudinal covariate. In classic non-parametric regression analysis of longitudinal data, we get the $O(n^{2/5})$ rate of convergence. Can we achieve the same rate of convergence as the estimation of $\beta(t)$ with mixed synchronous and asynchronous longitudinal covariates? In this subsection, we provide an affirmative answer to this question.

Specifically, we propose a two-step approach to estimate $\beta(t)$ and $\gamma(t),$ respectively. In the first step, we aim to obtain unbiased estimation of $\beta(t)$ from (2.1) using only $X(t)$ and $Y(t)$. There are two strategies to handle the part that involves $\gamma(t).$ We can either remove it or take it into account by a non-parametric function. In the second step, we regress residuals from the first step on the asynchronous longitudinal covariates $Z(t)$ to get regression function estimation of $\gamma(t)$ similar to that in Cao et al. (2015).

We need a critical assumption

$$E\{Z(t)\mid X(t)\}=E\{Z(t)\}.$$

(2.3)

Condition (2.3) states that $X(t)$ and $Z(t)$ are uncorrelated when $X(t)\neq 0$ for any $t.$ This condition is plausible when $X(t)$ indicates a time-invariant treatment, which is randomized and independent of other covariates. This condition means that there is no unmeasured confounder between treatment and outcome. This assumption can be strong in observational studies, where unmeasured confounders abound. In classic linear models, when the covariates are orthogonal with each other, the estimation of regression coefficients remains unbiased with a larger variance if some covariates are omitted. Unlike linear models, bias occurs when orthogonal covariates are omitted for longitudinal data if a non-parametric intercept is not added. If $X(t)$ and $Z(t)$ are correlated, there will be bias using the proposed two-step method.

A centering approach to estimate $\beta(t)$.  From (2.1), taking a conditional expectation on $X(t),$ by (2.3), we have

$$E\{{Y}(t)\mid X(t)\}={X}(t)^{T}\beta(t)+\{EZ(t)\}^{T}\gamma(t).$$

(2.4)

Taking another expectation, we get unconditional expectation of $Y(t)$ as

$$E\{Y(t)\}=\{EX(t)\}^{T}\beta(t)+\{EZ(t)\}^{T}\gamma(t).$$

(2.5)

By taking the difference between (2.4) and (2.5), we obtain

$$E\{\tilde{Y}(t)\mid\tilde{X}(t)\}=\tilde{X}(t)^{T}\beta(t),$$

(2.6)

where $\tilde{Y}(t)=Y(t)-E\{Y(t)\}$ and $\tilde{X}(t)=X(t)-E\{X(t)\}.$ Unbiased estimation of $\beta(t)$ can be obtained through (2.6). In (2.6), $m_{Y}(t)=E\{Y(t)\}$ and $m_{X}(t)=E\{X(t)\}$ are not directly observed. We propose Nadaraya-Watson type of estimators as follows (Nadaraya, 1964; Watson, 1964). Let $\hat{m}_{Y}(t_{0})=\{\sum\limits_{i=1}^{n}\int K_{h}(t-t_{0})Y_{i}(t)dN_{i}(t)\}/\{\sum\limits_{i=1}^{n}\int K_{h}(t-t_{0})dN_{i}(t)\},$ and $\hat{m}_{X}(t_{0})=\{\sum\limits_{i=1}^{n}\int K_{h}(t-t_{0})X_{i}(t)dN_{i}(t)\}/\{\sum\limits_{i=1}^{n}\int K_{h}(t-t_{0})dN_{i}(t)\},$ where $K(\cdot)$ is a kernel function with bounded support, $K(\cdot)\geq 0,\int_{-\infty}^{\infty}K(t)dt=1$ and $K_{h}(t)=h^{-1}K(t/h).$ The bandwidth $h\rightarrow 0$ and $nh\rightarrow\infty.$ Denote $\hat{Y}_{i}(t)=Y_{i}(t)-\hat{m}_{Y}(t)$ and $\hat{X}_{i}(t)=X_{i}(t)-\hat{m}_{X}(t).$ We aim to find

$$\begin{split}&\{\hat{\beta}_{c}(t),\hat{\dot{\beta}}_{c}(t)\}\\ =&\mathop{\rm arg\min}_{\beta_{0},\beta_{1}\in R^{p}}\Big{[}\frac{1}{n}\sum_{i=1}^{n}\int K_{h}(t_{1}-t)\left\{\hat{Y}_{i}(t_{1})-\hat{X}_{i}(t_{1})^{T}\beta_{0}-\hat{X}_{i}(t_{1})^{T}\beta_{1}(t_{1}-t)\right\}^{2}dN_{i}(t_{1})\Big{]}.\end{split}$$

(2.7)

Define

$$\begin{split}S_{n,l}(t)&=\frac{1}{nh}\sum\limits_{i=1}^{n}\int K_{h}(t_{1}-t)\hat{X}_{i}(t_{1})\hat{X}_{i}(t_{1})^{T}\{(t_{1}-t)/h\}^{l}dN_{i}(t_{1}),l=0,1,2\\ \mbox{and}\quad q_{n,l}(t)&=\frac{1}{nh}\sum\limits_{i=1}^{n}\int K_{h}(t_{1}-t)\hat{X}_{i}(t_{1})\hat{Y}_{i}(t_{1})\{(t_{1}-t)/h\}^{l}dN_{i}(t_{1}),l=0,1.\end{split}$$

Then

$$\begin{pmatrix}\hat{\beta}_{c}(t)\\ h\hat{\dot{\beta}}_{c}(t)\end{pmatrix}=\begin{pmatrix}S_{n,0}&S_{n,1}\\ S_{n,1}&S_{n,2}\end{pmatrix}^{-1}\begin{pmatrix}q_{n,0}\\ q_{n,1}\end{pmatrix}.$$

Next, we show the asymptotic properties of $\hat{\beta}_{c}(t)$. Denote $\tilde{W}_{i}(t_{1},t)=\{\tilde{X}_{i}(t_{1})^{T},\tilde{X}_{i}(t_{1})^{T}(t_{1}-t)\}^{T},$ $\theta_{0}(t)=\{\beta_{0}(t)^{T},\dot{\beta}_{0}(t)^{T}\}^{T},$ $\hat{\theta}_{c}(t)=\{\hat{\beta}_{c}(t)^{T},\hat{\dot{\beta}}_{c}(t)^{T}\}^{T},$ and $\mbox{var}\{\tilde{Y}(t)\mid\tilde{X}(t)\}=\sigma\{t,\tilde{X}(t)\}^{2}.$ We need the following assumptions.

1. (A5)

   $N_{i}(t)$ is independent of $\left\{X_{i}(t),Y_{i}(t)\right\},$ and $E\{dN_{i}(t)\}=\lambda(t)dt,i=1,\dots,n,$ where $\lambda(t)$ is twice continuously differentiable for $\forall t\in[0,1].$

2. (A6)

   $\theta_{0}(t)$ is twice continuously differentiable for $\forall t\in[0,1]$. For any fixed time point $t,$ $E\{\tilde{W}(t_{1},t)\}$ and $E\{\tilde{W}(t_{1},t)\tilde{W}(s_{1},t)^{T}\}$ are twice continuously differentiable for $(t_{1},s_{1})\in[0,1]^{\otimes 2}$. Furthermore, $E\{\tilde{X}(t)\tilde{X}(t)^{T}\}\lambda(t)$ is positive definite for $\forall t\in[0,1]$. Furthermore,

   $$\|E[\tilde{X}(t)\tilde{X}(t)^{T}\sigma\{t,\tilde{X}(t)\}^{2}]\|_{\infty}\lambda(t)<\infty,$$

   where for a square matrix $A,$ $\|A\|_{\infty}=\max_{1\leq i\leq n}\sum_{j=1}^{n}|a_{ij}|.$

3. (A7)

   $K(\cdot)$ is a symmetric density function satisfying $\int zK(z)dz=0,$ $\int z^{2}K(z)dz<\infty$ and $\int K(z)^{2}dz<\infty.$

4. (A8)

   $nh\to\infty$ and $nh^{5}\to 0.$

Conditions (A5)-(A8) are similar in spirit to conditions (A1)-(A4) in Section 2.1. (A5) is the univariate version of (A1). (A6) is a modified identifiability condition for $\theta_{0}(t).$ (A7) and (A8) are requirements on the kernel function and the bandwidth.

We state the asymptotic distribution of $\hat{\beta}_{c}(t)$ in the following theorem and relegate the proofs in the Appendix.

We need one bandwidth $h$ and the bias of $\beta_{c}(t)$ in Theorem 2 is $O(h^{2})$, which is the same as in classic nonparametric regression such as Theorem 1 of Fan and Zhang (2008). Our method depends on the selection of the bandwidth, which strikes a balance between bias and variability. From (A8), we shall choose the bandwidth $h=o(n^{-1/5}).$ With this choice of the bandwidth, we achieve a convergence rate $O(n^{2/5}),$ the same as non-parametric regression of classic longitudinal data (Fan and Zhang, 2008) and faster than that obtained in Theorem 1. It is counterintuitive that we get an improved rate of convergence with less information – we do not use any information contained in $Z_{i}(S_{ik}),i=1,\ldots,n;k=1,\ldots,M_{i}.$ The key is the orthogonality assumption imposed in (2.3), which allows us to obtain an unbiased estimation of the non-parametric regression function $\beta_{0}(t).$ When (2.3) is violated, we can regress the residuals of $Y(t)$ on $Z(t)$ on the residuals of $X(t)$ on $Z(t)$ similar to the FWL theorem in linear models (Frisch and Waugh, 1933; Lovell, 1963; Ding, 2021). We conjecture that the results would be similar to the simultaneous estimation of $\beta(t)$ and $\gamma(t),$ with a slower convergence rate for the estimation of $\beta(t).$

For statistical inference, we need to estimate the variance of $\hat{\beta}_{c}(t).$ We use the sandwich formula to estimate the variance covariance matrix of $\{\hat{\beta}_{c}(t),\hat{\dot{\beta}}_{c}(t)\}$ as follows:

$$\begin{split}\widehat{\mathrm{var}}\{\hat{\theta}_{c}(t)\}=&\left\{\sum_{i=1}^{n}\iint K_{h}(t_{1}-t)\hat{W}_{i}(t_{1},t)\hat{W}_{i}(t_{1},t)^{T}dN_{i}(t_{1})\right\}^{-1}\\ &\times\sum_{i=1}^{n}\left[\iint K_{h}(t_{1}-t)\hat{W}_{i}(t_{1},t)\left\{\hat{Y}_{i}(t_{1})-\hat{W}_{i}(t_{1},t)^{T}\hat{\theta}_{c}(t)\right\}dN_{i}(t_{1})\right]^{\otimes 2}\\ &\times\left\{\sum_{i=1}^{n}\iint K_{h}(t_{1}-t)\hat{W}_{i}(t_{1},t)\hat{W}_{i}(t_{1},t)^{T}dN_{i}(t_{1})\right\}^{-1},\end{split}$$

where $\hat{W}_{i}(t_{1},t)=\{\hat{X}_{i}(t_{1})^{T},\hat{X}_{i}(t_{1})^{T}(t_{1}-t)\}^{T},\hat{X}_{i}(t)=X_{i}(t)-\hat{m}_{X}(t),$ and $\hat{Y}_{i}(t)=Y_{i}(t)-\hat{m}_{Y}(t).$ The variance covariance matrix of $\hat{\beta}_{c}(t)$ is the upper left $p\times p$ submatrix of $\widehat{\mathrm{var}}\{\hat{\theta}_{c}(t)\}$ while we treat $\dot{\beta}_{c}(t)$ as a nuisance function. In the literature on density estimation, $\hat{\dot{\beta}}_{c}(t)$ estimated from (2.7) has more desirable properties than conventional kernel density estimation. More details on this can be found in Cattaneo et al. (2020).