A Score Based Test for Functional Linear Concurrent Regression

TL;DR

This paper introduces a new score-based test for assessing the effect of covariates in functional linear concurrent regression, offering improved power and robustness over existing methods, with theoretical and empirical validation.

Contribution

It develops a novel score test for zero variance components in functional regression, with theoretical justification and superior performance demonstrated through simulations and real data applications.

Findings

Test maintains correct type I error rate.

Achieves higher power than bootstrapped F test.

Performs well with sparse and noisy data.

Abstract

We propose a novel method for testing the null hypothesis of no effect of a covariate on the response in the context of functional linear concurrent regression. We establish an equivalent random effects formulation of our functional regression model under which our testing problem reduces to testing for zero variance component for random effects. For this purpose, we use a one-sided score test approach, which is an extension of the classical score test. We provide theoretical justification as to why our testing procedure has the right levels (asymptotically) under null using standard assumptions. Using numerical simulations, we show that our testing method has the desired type I error rate and gives higher power compared to a bootstrapped F test currently existing in the literature. Our model and testing procedure are shown to give good performances even when the data is sparsely observed, and the covariate is contaminated with noise. Applications of the proposed testing method are demonstrated on gait study and a dietary calcium absorption data.

Tables

Table 1. Table 1: Estimated type 1 error along with their standard error at α = 5 % 𝛼 percent 5 \alpha=5\% level.

Simulation Scenario	Scenario A		Scenario B
Sampling Design	$n=100$	$n=300$	$n=100$	$n=300$
Dense	.038 (.006)	.041 (.006)	.049 (.007)	.046 (.007)
Sparse	.064 (.008)	.047 (.007)	.059 (.007)	.042 (.006)

Equations

$$Y_i(t)={\beta }_0(t)+\sum _{k=1}^p{X}_{ik}(t){\beta }_k(t)+{ϵ}_i(t),$$

$$Y_i(t)={\beta }_0(t)+X_i(t){\beta }_1(t)+{ϵ}_i(t),$$

$$H_0:{\beta }_1(t)=0\; for\; alltversus{H}_1:{\beta }_1(t)\ne 0\; for\; somet.$$

$$\sum _{i=1}^n||Y_i(\cdot )-{\beta }_0(\cdot )-X_i(\cdot ){\beta }_1(\cdot )|{|}_{F_2}^2+{\lambda }_0\int \{{\beta }_0^r(t){\}}^{2}dt+{\lambda }_1\int \{{\beta }_1^r(t){\}}^{2}dt,$$

$$\sum _{i=1}^n\sum _{j=1}^m\{Y_i(t_j)-{\mathbf{B}}_0^T(t_j){\mathbf{b}}_0-{{\mathbf{X}}_{\mathbf{i}}^{\ast }}^T(t_j){\mathbf{b}}_1{\}}^2+{\lambda }_0{\mathbf{b}}_0^T\text{amsmathbb}{P}_0{\mathbf{b}}_0+{\lambda }_1{\mathbf{b}}_1^T\text{amsmathbb}{P}_1{\mathbf{b}}_1.$$

$$\sum _{i=1}^n||{\mathbf{Y}}_i-\text{amsmathbb}{B}_0{\mathbf{b}}_0-\text{amsmathbb}{X}_i{\mathbf{b}}_1|{|}_2^2+{\lambda }_0{\mathbf{b}}_0^T\text{amsmathbb}{P}_0{\mathbf{b}}_0+{\lambda }_1{\mathbf{b}}_1^T\text{amsmathbb}{P}_1{\mathbf{b}}_1.$$

$${\mathbf{Y}}_i=\text{amsmathbb}{C}_0{\mathbf{\gamma }}_0+\text{amsmathbb}{Z}_i{\mathbf{\gamma }}_1+{\mathbf{ϵ}}_i,$$

$$\mathbf{Y}=\text{amsmathbb}B{\mathbf{\gamma }}_0+\text{amsmathbb}Z{\mathbf{\gamma }}_1+\mathcal{E},$$

$$L_{ML}({\tau }_0,{\tau }_1)=-1/2(\ln |\text{amsmathbb}V|+{\mathbf{Y}}^T\text{amsmathbb}{V}^{-1}\mathbf{Y}).$$

$$S_{{\tau }_1}({\tau }_0,{\tau }_1)$$

$$=-1/2\{tr(\text{amsmathbb}{Z}^T\text{amsmathbb}{V}^{-1}\text{amsmathbb}Z)-(\text{amsmathbb}{V}^{-1/2}\mathbf{Y}{)}^T\text{amsmathbb}{V}^{-1/2}\text{amsmathbb}Z\text{amsmathbb}{Z}^T\text{amsmathbb}{V}^{-1/2}(\text{amsmathbb}{V}^{-1/2}\mathbf{Y})\},$$

$$\text{amsmathbb}I(\mathbf{\theta })=\left\{\begin{array}{cc}{I}_{11}(\mathbf{\theta })& I_{12}(\mathbf{\theta })\\ I_{21}(\mathbf{\theta }{)}^T& I_{22}(\mathbf{\theta })\end{array}\right\},$$

$$T_S=\frac{S_{{\tau }_1}^2(\tilde{\mathbf{\theta }})}{I_{22}(\tilde{\mathbf{\theta }})-I_{21}(\tilde{\mathbf{\theta }}{)}^T{I}_{11}^{-1}(\tilde{\mathbf{\theta }})I_{12}(\tilde{\mathbf{\theta }})}.$$

$$T_S=\{\begin{array}{l}\frac{S_{{\tau }_1}^2(\tilde{\mathbf{\theta }})}{\Lambda (\tilde{\mathbf{\theta }})}\text{if Sτ1(θ~)≥0}\\ 0\text{if Sτ1(θ~)<0},\end{array}$$

$$(1/2{)}^2\frac{{\left({\sum }_{\ell =1}^{k_1}\widetilde{{\lambda }_{\ell }}{x}_{\ell }^2-{\sum }_{\ell =1}^{k_1}\widetilde{{\lambda }_{\ell }}\right)}^2}{{\Lambda }_n(\tilde{\mathbf{\theta }})}I(\sum _{\ell =1}^{k_1}\widetilde{{\lambda }_{\ell }}{x}_{\ell }^2\ge \sum _{\ell =1}^{k_1}\widetilde{{\lambda }_{\ell }}),$$

$$G(s,t)=\sum _{k=1}^{\infty }{\lambda }_k{\varphi }_k(s){\varphi }_k(t),$$

$$\hat{\Sigma }(s,t)=\sum _{k=1}^K\widehat{{\lambda }_k}\widehat{{\varphi }_k}(s)\widehat{{\varphi }_k}(t)+\widehat{{\sigma }^2}I(s=t),$$

$$L_{ML}({\tau }_0,{\tau }_1,\dots ,{\tau }_p)=-1/2(\ln |\text{amsmathbb}V|+{\mathbf{Y}}^T\text{amsmathbb}{V}^{-1}\mathbf{Y}),$$

$$\text{amsmathbb}I(\mathbf{\theta })=\left\{\begin{array}{cc}\text{amsmathbb}{I}_{11}(\mathbf{\theta }{)}_{p\times p}& {\mathbf{I}}_{12}(\mathbf{\theta }{)}_{p\times 1}\\ {\mathbf{I}}_{21}(\mathbf{\theta }{)}_{1\times p}^T& I_{22}(\mathbf{\theta }{)}_{1\times 1}\end{array}\right\}.$$

$$Y_i(t)={\beta }_0(t)+X_i(t){\beta }_1(t)+{ϵ}_i(t),$$

$${ϵ}_i(t)={\xi }_{i1}\sqrt[]{2}cos(\pi t)+{\xi }_{i2}\sqrt[]{2}sin(\pi t)+N(0,0.9^2{I}_{m_i}),$$

$$Y_i(t)={\beta }_0(t)+X_{i1}(t){\beta }_1(t)+X_{i2}(t){\beta }_2(t)+{ϵ}_i(t),$$

$$\frac{S_{{\tau }_1}^2({\mathbf{\theta }}_0)}{\Lambda ({\mathbf{\theta }}_0)}$$

$$=\frac{{\left[1/2\{tr(\frac{\text{amsmathbb}{Z}^T\text{amsmathbb}{V}^{-1}\text{amsmathbb}Z}{n})-(\text{amsmathbb}{V}^{-1/2}\mathbf{Y}{)}^T\frac{\text{amsmathbb}{V}^{-1/2}\text{amsmathbb}Z\text{amsmathbb}{Z}^T\text{amsmathbb}{V}^{-1/2}}{n}(\text{amsmathbb}{V}^{-1/2}\mathbf{Y})\}\right]}^2}{\frac{\Lambda ({\mathbf{\theta }}_0)}{n^2}}$$

$$=\frac{{\left[1/2\{{\sum }_{\ell =1}^{k_1}{\lambda }_{\ell }-(\text{amsmathbb}W{)}^T(\frac{\text{amsmathbb}{V}^{-1/2}\text{amsmathbb}Z\text{amsmathbb}{Z}^T\text{amsmathbb}{V}^{-1/2}}{n})(\text{amsmathbb}W)\}\right]}^2}{{\Lambda }_n({\mathbf{\theta }}_0)}$$

$$=\frac{{\left[1/2\{{\sum }_{\ell =1}^{k_1}{\lambda }_{\ell }-(\text{amsmathbb}W{)}^T\text{amsmathbb}U\text{amsmathbb}D\text{amsmathbb}{U}^T(\text{amsmathbb}W)\}\right]}^2}{{\Lambda }_n({\mathbf{\theta }}_0)}.$$

$$\frac{S_{{\tau }_1}^2({\mathbf{\theta }}_0)}{\Lambda ({\mathbf{\theta }}_0)}$$

$$=\frac{\{1/2({\sum }_{\ell =1}^{k_1}{\lambda }_{\ell }-\text{amsmathbb}{X}^T\text{amsmathbb}D\text{amsmathbb}X){\}}^2}{{\Lambda }_n({\mathbf{\theta }}_0)}$$

$$=\frac{\{1/2({\sum }_{\ell =1}^{k_1}{\lambda }_{\ell }-{\sum }_{\ell =1}^{k_1}{\lambda }_{\ell }{x}_{\ell }^2){\}}^2}{{\Lambda }_n({\mathbf{\theta }}_0)}.$$

$$T_S({\mathbf{\theta }}_0,\Sigma )=\frac{S_{{\tau }_1}^2({\mathbf{\theta }}_0)}{\Lambda ({\mathbf{\theta }}_0)}I(S_{{\tau }_1}({\mathbf{\theta }}_0)\ge 0)\stackrel{d}{=}(1/2{)}^2\frac{{\left({\sum }_{\ell =1}^{k_1}{\lambda }_{\ell }{x}_{\ell }^2-{\sum }_{\ell =1}^{k_1}{\lambda }_{\ell }\right)}^2}{{\Lambda }_n({\mathbf{\theta }}_0)}I(\sum _{\ell =1}^{k_1}{\lambda }_{\ell }{x}_{\ell }^2\ge \sum _{\ell =1}^{k_1}{\lambda }_{\ell }),$$

Full text

††articletype: Original Paper

1 Introduction

Functional linear concurrent regression model arises when the response and covariates are both functions of time (or any continuous index), and the value of the response at a particular time point is modeled as a linear combination of the covariates at that specific time point, where the coefficients of the functional covariates are functions of time (Ramsay and Silverman 2005). One can view the functional linear concurrent regression model as a series of linear regression for each time point, with the assumption that the coefficient functions are smooth over time. Multiple methods exist in literature for estimation of these regression coefficient functions in functional linear concurrent regression and the closely related varying coefficient model (Hastie and Tibshirani 1993) using basis functions with roughness penalty (Ramsay and Silverman 2005), polynomial spline (Huang et al. 2002, 2004), local polynomial smoothing (Wu et al. 1998; Cai et al. 2000; Fan and Zhang 2000), Bayesian modeling (Gelfand et al. 2003), covariance representation techniques (Şentürk and Nguyen 2011), among many others. While estimation of the regression functions is an important problem, in many cases the primary interest might be finding out whether a specific covariate is truly significant or not, i.e., to test for association between a predictor of interest and the response. For example, in the gait study data (Ramsay and Silverman 2005), where there are longitudinal measurements of hip and knee angles taken on 39 children, the main purpose of the study is to understand how the joints in hip and knee interact during a gait cycle (Theologis 2009). One natural question to ask here would be, whether the knee angles (response) are at all associated with the hip angles (covariate). Simply building a point-wise confidence interval of the estimated regression function does not answer the question of the overall significance of the covariate. Thus there is a need for developing testing methods to find out significant predictors in this setting.

Formally, our primary goal is to test the null hypothesis that the coefficient function corresponding to a predictor of interest is identically zero, versus the alternative hypothesis that the coefficient function is non-zero for some time point. Literature relating to such global testing in functional concurrent linear regression or closely related varying coefficient model can be traced back to Huang et al. (2002), Guo (2002) among many others. Huang et al. (2002) employed a resampling subject bootstrap method on an F-type statistic, whereas Guo (2002) used the connection between linear mixed effects models and smoothing splines, and subsequently used a generalized maximum likelihood ratio test to test for significance of predictors. Kim et al. (2018) extended the bootstrap based test (Huang et al. 2002) to general nonlinear functional concurrent model. Both of these tests rely on a subject-level bootstrap method to obtain the p-values making the test computationally intensive. Recently, Wang et al. (2018) developed a method for pointwise as well as global testing using empirical likelihood ratio tests. Their method, which is extremely general, uses a wild bootstrap procedure to perform the test for both dense and sparse functional data. Besides these global testing procedures one can also use confidence bands based methods e.g., Fan and Zhang (2000) for building simultaneous confidence bands for the underlying coefficient functions.

In this article, our goal is to build a classical likelihood based testing method for testing of the global effect of covariates, which is also computationally cheap. We model the unknown regression functions using B-spline basis functions and derive an equivalent random effects model. Under such a framework, we show that our testing problem reduces to testing for zero variance components for a set of random effects. There are multiple existing methods in the literature for testing for zero variance components. Crainiceanu and Ruppert (2004), Greven et al. (2008), and Staicu et al. (2014) considered testing for variance components using likelihood ratio test (LRT) and restricted likelihood ratio test (RLRT). The main challenge of such tests is that the null distribution is different (Crainiceanu and Ruppert 2004) from the commonly used $0.5{\chi}^{2}_{0}$ : $0.5{\chi}^{2}_{1}$ approximation or such mixtures of two chi-square distributions, which is used in Guo (2002). Crainiceanu and Ruppert (2004) showed the large sample chi-square mixture approximations using the usual asymptotic theory (Self and Liang 1987; Liang and Self 1996) for null hypothesis on the boundary of the parameter space do not hold in general as the inherent assumptions are not valid for linear mixed models. In this article, we propose a score based testing method that is computationally efficient. Our procedure is inspired from the work of Molenberghs and Verbeke (2007) which describes an approach of using a one-sided score test in constrained parameter space. The major advantage of working with the score test is, it does not require computations under the alternative. Zhang and Lin (2008) and Lin (1997) also used such one-sided score tests for variance component testing in generalized linear mixed models and longitudinal data. However, the methods mentioned above assumed that the responses are independent given the random effects and that the variances have some parametric form. In contrast, in our functional regression framework, we assume unknown non-trivial covariance structure and estimate the covariance function nonparametrically. The assumption of non-trivial dependence is crucial in functional data because of complex correlation structures that might be present in real data. We derive the asymptotic distribution of the test statistic under the null hypothesis; we show that the commonly used chi-squared approximation of the score test statistic is not appropriate in our situation. However, the null distribution of our test statistic is easy to simulate from. Thus the calculation of p-value for our testing procedure is computationally efficient. We show that asymptotically our testing procedure has the correct type I error rate. Using numerical simulations, we illustrate that our testing method has the desired type I error rates for finite sample sizes and that our proposed testing procedure has higher power than the bootstrapped F-test of Kim et al. (2018).

The rest of the article is organized as follows. In Section 2, we discuss our model specification, present our testing method and derive theoretical properties related to our test statistic. In Section 3, we present a simulation study under various sampling design scenarios and give the simulation results. In Section 4, we demonstrate our proposed test by applying it to the two real data examples: gait data and calcium absorption study and summarize our findings. We conclude by a discussion about some limitations and some possible extensions of our work in Section 5.

2 Methodology

2.1 Modeling framework

Suppose that the observed data for the $i$th subject, $i=1,\ldots,n$, is $\{Y_{i}(t),X_{i1}(t),\ldots,X_{ip}(t)\}$, where $Y(\cdot)$ is a functional response and $X_{1}(\cdot)$,$\ldots,X_{p}(\cdot)$ are the corresponding functional covariates. In practice, the functions for the $i$th subject are observed only on a finite set of points $t_{ij}$, $j=1,\ldots,m_{i}$. We assume that $t_{ij}\in{\mathcal{T}}$, a bounded and closed set. For the rest of the article, we assume $\mathcal{T}=[0,1]$ without loss of generality. To start with, we will assume that $t_{ij}=t_{j}$, and that the covariates $X_{ik}(\cdot)$ are measured without error. We discuss the cases when the functions are observed on irregularly spaced grid, and with additional measurement errors in Section 2.5. We consider a linear functional concurrent regression model,

[TABLE]

where $\beta_{k}(t)$, $k=0,1,2,\ldots,p$, are smooth functions representing functional intercept and functional slope parameters, respectively. We assume $X_{ik}(\cdot)$ are independent and identically distributed (i.i.d.) copies of $\mathcal{X}_{k}(\cdot)$, $k=1,2,\ldots,p$, where $\mathcal{X}_{k}(\cdot)$ is a stochastic processes with finite second moment. For simplicity we illustrate our testing method for the single covariate model

[TABLE]

which can be easily extended to the multiple covariate situation above and this is discussed in Section 2.6. We further assume $\epsilon_{i}(\cdot)$ are i.i.d. copies of $\epsilon(\cdot)$, which is a mean zero Gaussian process plus some Gaussian white noise, that is, $\epsilon(t)=V(t)+w_{t},$ where $V(\cdot)\sim\mathcal{N}(0,G(\cdot,\cdot))$ and $w_{t}$ are i.i.d. $\mathcal{N}(0,\sigma^{2})$ random errors. Thus the covariance function of the error process is given by $\Sigma(s,t)=cov(\epsilon(s),\epsilon(t))=G(s,t)+\sigma^{2}I(s=t)$. Our primary interest lies in testing,

[TABLE]

In general testing $H_{0}$ is difficult since $\beta_{1}(t)$ is an infinite dimensional parameter. In this article, we show that by modeling the coefficient function with splines and using a random effects model, the testing problem can be reduced to a variance component test. This reduction in parameter dimension not only helps in getting satisfactory performance of our testing method but also is computationally cheaper than doing the bootstrapped F test existing in literature. Subsequently we develop a one sided score test for testing our null hypothesis.

2.2 Equivalent random effects model

An usual method (Ramsay and Silverman 2005) to estimate $\beta_{0}(t)$ and $\beta_{1}(t)$ in model (1) is by minimizing the penalized residual sum of squares,

[TABLE]

where $\lambda_{0}$, $\lambda_{1}$ are unknown penalty parameters penalizing $r$-th derivative of the coefficient functions and $||\cdot||_{F_{2}}$ denotes the functional $L_{2}$ norm. Suppose for $\ell=0,1$ $\{B_{k\ell}(t),k=1,2,\ldots,k_{\ell}\}$ is a set of known basis functions. We approximate the unknown coefficient functions using finite basis function expansion as $\beta_{\ell}(t)=\sum_{k=1}^{K_{\ell}}b_{k\ell}B_{k\ell}(t)=\mathbf{B}_{\ell}^{T}(t)\mathbf{b}_{\ell},\hskip 5.69054pt\ell=0,1$, where $B_{\ell}(t)=[B_{1\ell}(t),B_{2\ell}(t),\ldots,B_{K_{\ell}\ell}(t)]^{T}$ and $b_{\ell}=(b_{1\ell},b_{2\ell},\ldots,b_{k_{\ell}\ell})^{T}$ is a vector of unknown coefficients. In this article, we use cubic B-spline basis functions, however, other basis functions can be used as well. Thus we can write $X_{i}(t)\beta_{1}(t)=\sum_{k=1}^{k_{1}}b_{k1}X_{i}(t)B_{k1}(t)=\mathbf{X_{i}^{*}}^{T}(t)\mathbf{b}_{1}$, where $\mathbf{X}_{i}^{*}(t)=[X_{i}(t)B_{11}(t),X_{i}(t)B_{21}(t),\ldots,X_{i}(t)B_{k_{1}1}(t)]^{T}$. We can then rewrite our model (1) as $Y_{i}(t)=\mathbf{B}_{0}^{T}(t)\mathbf{b}_{0}+\mathbf{X_{i}^{*}}^{T}(t)\mathbf{b}_{1}+\epsilon_{i}(t)$. The unknown basis coefficients can then be estimated by minimizing the penalized error sum of squares $\sum_{i=1}^{n}||Y_{i}(\cdot)-\mathbf{B}_{0}^{T}(\cdot)\mathbf{b}_{0}-{\mathbf{X}_{i}^{*}}^{T}(\cdot)\mathbf{b}_{1}||_{2}^{2}+\lambda_{0}\mathbf{b}_{0}^{T}\amsmathbb{P}_{0}\mathbf{b}_{0}+\lambda_{1}\mathbf{b}_{1}^{T}\amsmathbb{P}_{1}\mathbf{b}_{1},$ where $\amsmathbb{P}_{0}$ and $\amsmathbb{P}_{1}$ are the penalty matrices coming from penalizing the $r$-th derivative of the functions $\beta_{0}(t)$ and $\beta_{1}(t)$. In particular $\amsmathbb{P}_{\ell}=\int\mathbf{B}^{r}_{\ell}(t)\mathbf{B}^{r}_{\ell}(t)^{T}dt$ and thus $\int(\beta_{\ell}^{r}(t))^{2}dt=\mathbf{b}_{\ell}^{T}\amsmathbb{P}_{\ell}\mathbf{b}_{\ell}$. Since we only observe data on a fine regular grid $S=\{t_{1},t_{2},\ldots,t_{m}\}$ in practice, the minimization is carried out by minimizing

[TABLE]

Define $\mathbf{Y}_{i}=[Y_{i}(t_{1}),Y_{i}(t_{2}),\ldots,Y_{i}(t_{m})]^{T}$, $\amsmathbb{B}_{0}=[\mathbf{B}_{0}(t_{1})|\mathbf{B}_{0}(t_{2})|\ldots|\mathbf{B}_{0}(t_{m})]^{T}$,

$\amsmathbb{X}_{i}=[\mathbf{X}_{i}^{*}(t_{1})|\mathbf{X}_{i}^{*}(t_{2})|\ldots|\mathbf{X}_{i}^{*}(t_{m})]^{T}$. Thus the least square criterion becomes

[TABLE]

Now since the matrices $\amsmathbb{P}_{0}$, $\amsmathbb{P}_{1}$ are singular (for $r\geq 1$), the equivalent random effects model corresponding to this minimization problem would be rank deficient. As our primary interest lies in testing, we propose to penalize the coefficient functions directly, namely we use $r=0$ and ridge penalty, consequently $\amsmathbb{P}_{\ell}=\int\mathbf{B}_{\ell}(t)\mathbf{B}_{\ell}(t)^{T}dt$. Our simulations show we are able to maintain correct type I error rates of the proposed testing method using this strategy. It then follows that the normal equations are identical to those from the equivalent random effects model $\mathbf{Y}_{i}=\amsmathbb{B}_{0}\mathbf{b}_{0}+\amsmathbb{X}_{i}\mathbf{b}_{1}+\bm{\zeta}_{i}$, where $\bm{\zeta}_{i}\sim\mathcal{N}_{m}(0,\sigma_{\epsilon}^{2}\amsmathbb{I}_{m})$, $\mathbf{b}_{0}\sim\mathcal{N}_{k_{0}}(0,\sigma_{0}^{2}\mathbb{\Sigma}_{0})$, $\mathbf{b}_{1}\sim\mathcal{N}_{k_{1}}(0,\sigma_{1}^{2}\mathbb{\Sigma_{1}})$ and $\mathbf{b}_{0}$, $\mathbf{b}_{1}$, $\zeta_{i}$ are independent. We denote $\mathbb{\Sigma}_{0}=\amsmathbb{P}_{0}^{-1}$, $\mathbb{\Sigma}_{1}=\amsmathbb{P}_{1}^{-1}$, $\sigma_{0}^{2}=\frac{\sigma_{\epsilon}^{2}}{\lambda_{0}}$, $\sigma_{1}^{2}=\frac{\sigma_{\epsilon}^{2}}{\lambda_{1}}$. Using Cholesky decomposition of $\mathbb{\Sigma}_{0}$, $\mathbb{\Sigma}_{1}$ and appropriately reparameterizing (${\amsmathbb{C}_{0}}={\amsmathbb{B}_{0}}\mathbb{\Sigma}_{0}^{1/2}$, $\amsmathbb{Z}_{i}=\amsmathbb{X}_{i}\mathbb{\Sigma}_{1}^{1/2}$, $\bm{\gamma}_{k}=\mathbb{\Sigma}_{k}^{-1/2}\mathbf{b}_{k}$), the model can be rewritten as $\mathbf{Y}_{i}=\amsmathbb{C}_{0}\bm{\gamma}_{0}+\amsmathbb{Z}_{i}\bm{\gamma}_{1}+{\bm{\zeta}_{i}},$ where $\bm{\zeta}_{i}\sim\mathcal{N}_{m}(0,\sigma_{\epsilon}^{2}\amsmathbb{I}_{m})$, $\bm{\gamma}_{0}\sim\mathcal{N}_{k_{0}}(0,\sigma_{0}^{2}\amsmathbb{I}_{k_{0}})$, $\bm{\gamma}_{1}\sim\mathcal{N}_{k_{1}}(0,\sigma_{1}^{2}\amsmathbb{I}_{k_{1}})$, and all the random effects are independent. Thus our test $H_{0}$ can be carried out via testing of a single variance component, namely testing $H_{0}:\sigma_{1}^{2}=0$ against the alternative $H_{1}:\sigma_{1}^{2}>0$.

Now for our testing problem the errors are not independent and more likely to be temporally correlated. Therefore to get the correct likelihood we need to use the true covariance kernel $\Sigma(s,t)$ which takes into account this temporal dependence of the residual vector $\bm{\epsilon}_{i}$s. This motivates us to use the following random effects model

[TABLE]

where $\bm{\epsilon}_{i}\sim\mathcal{N}_{m}(0,\mathbb{\Sigma}_{m\times m})$, $\bm{\gamma}_{0}\sim\mathcal{N}_{k_{0}}(0,\sigma_{0}^{2}\amsmathbb{I}_{k_{0}})$, $\bm{\gamma}_{1}\sim\mathcal{N}_{k_{1}}(0,\sigma_{1}^{2}\amsmathbb{I}_{k_{1}})$ and all of them are independent. Here $\mathbb{\Sigma}_{m\times m}$ denotes the covariance kernel $\Sigma(s,t)$ evaluated at $S=\{t_{1},t_{2},\ldots,t_{m}\}$. For the moment let us assume $\mathbb{\Sigma}_{m\times m}$ to be known. Of course in reality $\mathbb{\Sigma}_{m\times m}$ will be unknown and we will need to estimate it. We illustrate in Section 2.4 how to estimate $\mathbb{\Sigma}_{m\times m}$ using functional principal component analysis (FPCA). Writing equation (2) in stacked form for $i=1,2,3,\ldots,n$ we have

[TABLE]

where $\amsmathbb{B}=[\amsmathbb{C}_{0}^{T}|\amsmathbb{C}_{0}^{T}|\ldots|\amsmathbb{C}_{0}^{T}]^{T}$, $\amsmathbb{Z}=[\amsmathbb{Z}_{1}^{T}|\amsmathbb{Z}_{2}^{T}|\ldots|\amsmathbb{Z}_{n}^{T}]^{T}$, $\mathbf{Y}=(\mathbf{Y}_{1}^{T},\mathbf{Y}_{2}^{T},\ldots,\mathbf{Y}_{n}^{T})^{T}$ and $\mathbf{\mathcal{E}}=(\bm{\epsilon}_{1}^{T},\bm{\epsilon}_{2}^{T},\ldots,\bm{\epsilon}_{n}^{T})^{T}$, $\mathbf{\mathcal{E}}\sim\mathcal{N}_{N}(0,\mathbb{\Sigma})$ ($N=mn$), where $\mathbb{\Sigma}$=diag {${\mathbb{\Sigma}_{m\times m},\mathbb{\Sigma}_{m\times m},\ldots,\mathbb{\Sigma}_{m\times m}}$}. Note that $\amsmathbb{Z}=\amsmathbb{X}\mathbb{\Sigma}_{1}^{1/2}$ and $\amsmathbb{B}=\mathbf{\mathcal{B}}\mathbb{\Sigma}_{0}^{1/2}$, where $\amsmathbb{X},\mathbf{\mathcal{B}}$ are defined similarly by stacking $\amsmathbb{X}_{i}$ and $\amsmathbb{B}_{0}$ s. In this set up, we are interested in testing $H_{0}:\sigma_{1}^{2}=0$ against the alternative $H_{1}:\sigma_{1}^{2}>0$, which as demonstrated, is equivalent to testing the null hypothesis of no effect of a covariate on the response under functional linear concurrent model. Next we develop a score based test for conducting the test of variance component.

2.3 Testing method

We develop our testing method treating $\amsmathbb{Z}$ as nonrandom (fixed). Namely our test is a conditional test based on observed $\amsmathbb{Z}$ {i.e., observed $X_{i}(t)$}. We show that our conditional testing method has the right levels under null which in turn ensures the unconditional test would also enjoy this property. Marginally $\mathbf{Y}\sim\mathcal{N}(0,\amsmathbb{V})$, which follows from equation (3) with $\amsmathbb{V}=\amsmathbb{V}(\tau_{0},\tau_{1})=\mathbb{\Sigma}+\tau_{0}\amsmathbb{B}\amsmathbb{B}^{T}+\tau_{1}\amsmathbb{Z}\amsmathbb{Z}^{T}$, where $\tau_{0}=\sigma_{0}^{2}$ and $\tau_{1}=\sigma_{1}^{2}$. So the marginal log-likelihood of $\mathbf{Y}$ (upto a constant) is :

[TABLE]

Based on this likelihood, we want to test $H_{0}:\tau_{1}=0$ vs $H_{1}:\tau_{1}>0$.

Let $\bm{\theta}=(\tau_{0},\tau_{1})^{T}$ and $\tilde{\bm{\theta}}$ denote the maximum likelihood estimate (M.L.E) of $\bm{\theta}$ under $H_{0}$. The score function of $\tau_{1}$ is

[TABLE]

where $\amsmathbb{M}=\amsmathbb{Z}\amsmathbb{Z}^{T}$. The information matrix $\amsmathbb{I}(\bm{\theta})$ corresponding to the likelihood in (4) is partitioned as,

[TABLE]

where $I_{11}(\bm{\theta})=tr\{(\amsmathbb{B}^{T}\amsmathbb{V}^{-1}\amsmathbb{B})^{2}\}/2,I_{22}(\bm{\theta})=tr\{(\amsmathbb{Z}^{T}\amsmathbb{V}^{-1}\amsmathbb{Z})^{2}\}/2$, $I_{21}(\bm{\theta})=I_{12}(\bm{\theta})\\ =tr\{(\amsmathbb{B}^{T}\amsmathbb{V}^{-1}\amsmathbb{Z})(\amsmathbb{B}^{T}\amsmathbb{V}^{-1}\amsmathbb{Z})^{T}\}/2$. Then the classical score test statistic is given by,

[TABLE]

As the parameter space is constrained and the null hypothesis is on the boundary of the parameter space, following Molenberghs and Verbeke (2007) we define our one sided score test statistic as

[TABLE]

where $\Lambda(\tilde{\bm{\theta}})=I_{22}(\tilde{\bm{\theta}})-I_{21}(\tilde{\bm{\theta}})^{T}I_{11}^{-1}(\tilde{\bm{\theta}})I_{12}(\tilde{\bm{\theta}}).$ The reason behind using such one sided score test statistic is if the score is negative at the M.L.E under null, then the value of the score gives no evidence in favour of the alternative and therefore the statistic is set to zero in such cases. An illustration for constrained one parameter case is given in Figure 1.

We assume the true covariance matrix $\mathbb{\Sigma}$ to be known for the time being for establishing asymptotic distribution of our test statistic but in reality it is generally unknown, so we will need to estimate it from data by some consistent estimator $\hat{\mathbb{\Sigma}}$ and plug that in for $\mathbb{\Sigma}$ in $T_{S}$. Next, we posit two theorems which establish asymptotic distribution our test statistic.

Theorem 2.1**.**

Suppose the following conditions are true :

1. a)

The null hypothesis is true i.e., $H_{0}:\tau_{1}=0$ holds and $\bm{\theta}_{0}=(\tau_{0}^{*},0)$ is the true value of $\bm{\theta}$, 2. b)

$\mathbb{\Sigma}$* be the true covariance matrix of the residual vector $\mathbf{\mathcal{E}}$ in equation (3).*

*Then $T_{S}(\bm{\theta}_{0},\mathbb{\Sigma})=\frac{S_{\tau_{1}}^{2}({\bm{\theta}_{0}})}{\Lambda({\bm{\theta}_{0}})}I(S_{\tau_{1}}({\bm{\theta}_{0}})\geq 0)\overset{d}{=}(1/2)^{2}\frac{\left(\sum_{\ell=1}^{k_{1}}\lambda_{\ell}x_{\ell}^{2}-\sum_{\ell=1}^{k_{1}}\lambda_{\ell}\right)^{2}}{\Lambda_{n}(\bm{\theta}_{0})}I(\sum_{\ell=1}^{k_{1}}\lambda_{\ell}x_{\ell}^{2}\geq\sum_{\ell=1}^{k_{1}}\lambda_{\ell})$,

where $x_{\ell}\stackrel{{\scriptstyle iid}}{{\sim}}\mathcal{N}(0,1)$ and $\lambda_{\ell}$ are eigenvalues of $\amsmathbb{Z}^{T}\amsmathbb{V}(\theta_{0},\mathbb{\Sigma})^{-1}\amsmathbb{Z}/n$ and $\Lambda_{n}(\bm{\theta}_{0})=\Lambda(\bm{\theta}_{0})/n^{2}=\frac{1}{2}tr\{(\amsmathbb{Z}^{T}\amsmathbb{V}^{-1}\amsmathbb{Z}/n)^{2}\}$ -  $\frac{\left[\frac{1}{2}tr\{(\amsmathbb{B}^{T}\amsmathbb{V}^{-1}\amsmathbb{Z}/n)(\amsmathbb{B}^{T}\amsmathbb{V}^{-1}\amsmathbb{Z}/n)^{T}\}\right]^{2}}{\frac{1}{2}tr\{(\amsmathbb{B}^{T}\amsmathbb{V}^{-1}\amsmathbb{B}/n)^{2}\}}$*

The proof of the Theorem 2.1 is given in Appendix A.

Theorem 2.2**.**

Suppose the following conditions are true :

1. a)

The null hypothesis is true i.e $H_{0}:\tau_{1}=0$ holds and $\bm{\theta}_{0}=(\tau_{0}^{*},0)$ is the true value of $\bm{\theta}$, 2. b)

$\tilde{\bm{\theta}}$* is $\sqrt[]{n}$ consistent estimator of $\bm{\theta}_{0}$ under null and the estimator $\hat{\mathbb{\Sigma}}$ is a consistent estimator of $\mathbb{\Sigma}$ in the sense $||\hat{\mathbb{\Sigma}}^{-1}-\mathbb{\Sigma}^{-1}||_{2}=o_{p}(1)$ (spectral norm).*

Then $T_{S}(\tilde{\bm{\theta}},\hat{\mathbb{\Sigma}})\overset{d}{\rightarrow}T_{S}(\bm{\theta}_{0},\mathbb{\Sigma})$.

The proof is mainly based on application of Slutsky’s theorem and matrix norm inequalities. A detailed proof is given in Appendix S1 of supplemental material. As mentioned in Theorem 2.1 the null distribution of the test statistic is given by $(1/2)^{2}\frac{\left(\sum_{\ell=1}^{k_{1}}\lambda_{\ell}x_{\ell}^{2}-\sum_{\ell=1}^{k_{1}}\lambda_{\ell}\right)^{2}}{\Lambda_{n}(\theta_{0})}I(\sum_{\ell=1}^{k_{1}}\lambda_{\ell}x_{\ell}^{2}\geq\sum_{\ell=1}^{k_{1}}\lambda_{\ell})$. Because $\bm{\theta}_{0}$ and $\lambda_{\ell}$ are unknown in reality, we approximate the null distribution using plug-in estimates of $\tilde{\bm{\theta}}$ and $\hat{\mathbb{\Sigma}}$, i.e, we use the approximate null distribution

[TABLE]

where $\tilde{\lambda_{\ell}}$ are eigenvalues of $\amsmathbb{Z}^{T}\amsmathbb{V}(\tilde{\bm{\theta}},\hat{\mathbb{\Sigma}})^{-1}\amsmathbb{Z}/n$ and $x_{\ell}\stackrel{{\scriptstyle iid}}{{\sim}}\mathcal{N}(0,1)$ for $\ell=1,2,\ldots,k_{1}$. This is justified as it can be shown $\tilde{\lambda_{\ell}}\overset{p}{\rightarrow}\lambda_{\ell}$ and $\Lambda_{n}(\tilde{\bm{\theta}})\overset{p}{\rightarrow}\Lambda_{n}(\bm{\theta}_{0})$, see the proof in Appendix S1 of supplemental material. Our simulations show that we are able to get the correct type I error rates and good power of our test using the above strategy.

Simulation from the null distribution of our test statistic is easy and computationally efficient, as we only need to calculate $k_{1}$ eigenvalues of the matrix $\amsmathbb{Z}^{T}\amsmathbb{V}(\tilde{\bm{\theta}},\hat{\mathbb{\Sigma}})^{-1}\amsmathbb{Z}/n$, and simulate $x_{\ell}\stackrel{{\scriptstyle iid}}{{\sim}}\mathcal{N}(0,1)$ for $\ell=1,2,\ldots,k_{1}$. For calculation of $\amsmathbb{V}(\tilde{\bm{\theta}},\hat{\mathbb{\Sigma}})^{-1}$, we use the Woodbury matrix identity and also the fact $\hat{\mathbb{\Sigma}}$ is block diagonal, which greatly speeds up the calculation. It is well known (Zhang and Lin 2003), that the usual asymptotic ${\chi}^{2}$ distribution of the score test do not work here, so the approximate null distribution in (6) is more appropriate.

As the test statistic is one sided and in particular not continuous at zero, the p-value under null is asymptotically distributed as mixture distribution of degenerate one and $U(0,\alpha)$, where $\alpha=P_{H_{0}}\{S_{\tau_{1}}(\bm{\theta}_{0})\geq 0\}=P(\sum_{\ell=1}^{k_{1}}\lambda_{\ell}x_{\ell}^{2}\geq\sum_{\ell=1}^{k_{1}}\lambda_{\ell})$, see Appendix S2 in supplemental material for a detailed proof. As $k_{1}$ increases it follows by application of CLT, $\alpha\rightarrow\frac{1}{2}$ and the null distribution of our test statistic asymptotically converges to $0.5{\chi}^{2}_{0}$ : $0.5{\chi}^{2}_{1}$. In reality choice of $k_{1}$ will depend on the type of design (dense or sparse) and number of observed time points for each subjects. So such a convergence do not hold (Crainiceanu and Ruppert 2004) for subjects observed only on a finite set of points. Therefore, we use the approximate null distribution (6) to perform our test.

2.4 Estimation of Covariance Matrix

In reality $\mathbb{\Sigma}$ is unknown, and we need a consistent estimator $\hat{\mathbb{\Sigma}}$. In the context of functional data, we want to estimate $\Sigma(\cdot,\cdot)$ completely non-parametrically. If we had the original residuals $\epsilon_{ij}$ available, we could use functional principal component analysis (FPCA), e.g., Yao et al. (2005) or Zhang and Chen (2007) to estimate $\Sigma(s,t)$. The error process $\epsilon(t)$ was defined as $\epsilon(t)=V(t)+w_{t}$. We assume the covariance kernel $G(s,t)$ of the smooth part $V(t)$ is a Mercer kernel (Mercer 1909). Then by Mercer’s theorem $G(s,t)$ must have a spectral decomposition

[TABLE]

where $\lambda_{1}\geq\lambda_{2}\geq\ldots 0$ are the ordered eigenvalues and $\phi_{k}(\cdot)$s are corresponding eigenfunctions. Thus we have $\Sigma(s,t)=\sum_{k=1}^{\infty}\lambda_{k}\phi_{k}(s)\phi_{k}(t)+\sigma^{2}I(s=t)$. Given $\epsilon_{t_{ij}}=V(t_{ij})+w_{ij}$, one could employ FPCA based methods to get $\hat{\phi}_{k}(\cdot)$, $\hat{\lambda}_{k}$s and $\hat{\sigma^{2}}$. Hall et al. (2006) established $L^{2}$ convergence of the FPCA estimates, in particular of the eigenfunctions and eigenvalues, in both the sparse and dense functional data setting under appropriate regularity conditions on the sampling design. Li et al. (2010) established uniform convergence rates for eigenfunctions and eigenvalues under more general framework where the number of observations for each function can be sampled at any rate relative to the sample size. More specifically Li et al. (2010) showed it is possible to get consistent estimators $\hat{\phi}_{k}(\cdot)$, $\hat{\lambda}_{k}$ and $\hat{\sigma^{2}}$ under both sparse and dense functional data settings. So a consistent estimator of $\Sigma(s,t)$ can be formed as

[TABLE]

where $K$ is large enough for the convergence to hold and is typically chosen such that percent of variance explained (PVE) by the selected eigencomponents exceeds some pre-specified value such as $99\%$ or $95\%$. In reality we don’t have the original residuals $\epsilon_{ij}$ and use the full model (1) to obtain residuals $e_{ij}=Y_{i}(t_{j})-\hat{Y_{i}}(t_{j})$. Then treating $e_{ij}$ as our original residuals, we obtain $\hat{\Sigma}(s,t)$ using FPCA. Our simulations show good results using this approach and we are able to maintain the correct levels of the test under both sparse and dense sampling design scenarios.

2.5 Extension to Sparse and Noisy Covariate

In developing our method we assumed that covariates are measured without noise, and data is observed on a regular dense grid of points $S=\{t_{1},t_{2},\ldots,t_{m}\}\subset\mathcal{T}=[0,1]$. Although in reality data might be observed sparsely, and the covariates may be contaminated with measurement error. Our testing method can be extended to these situations in the following ways.

Case 1: Sparse design, no measurement error

We assume response $Y_{i}(t)$ and the covariate $X_{i}(t)$ are observed in $S_{i}=\{t_{i1},t_{i2},\ldots,t_{im_{i}}\}\subset\mathcal{T}=[0,1]$ for each $i=1,2,\ldots,n$ and $\max_{1\leq i\leq n}m_{i}\leq M$, for some fixed $M$. In this case the only difference in our model is that ${\mathbf{Y}_{i}}$ is a $m_{i}\times 1$ dimensional vector and that ${\bm{\epsilon}_{i}}\sim\mathcal{N}_{m_{i}}(0,\mathbb{\Sigma}_{i})$ independently, where $\mathbb{\Sigma}_{i}$ is the covariance kernel $\Sigma(s,t)$ evaluated at $S_{i}=\{t_{i1},t_{i2},\ldots,t_{im_{i}}\}$. So our model described in (2) still holds with $\mathbb{\Sigma}$=diag {${\mathbb{\Sigma}_{1},\mathbb{\Sigma}_{2},\ldots,\mathbb{\Sigma}_{n}}$}. As discussed in Section 2.4, if $\Sigma(s,t)$ can be consistently estimated by $\hat{\Sigma}(s,t)$ then $\hat{\mathbb{\Sigma}_{i}}$ would be a consistent estimator of $\mathbb{\Sigma}_{i}$ and our testing method is still valid.

Case 2: Dense design with measurement error

Suppose now data $Y_{i}(t)$ and covariate $X_{i}(t)$ are observed in a fine regular grid $S=\{t_{1},t_{2},\ldots,t_{m}\}\subset\mathcal{T}=[0,1]$ and the covariate is observed with measurement error. Namely instead of observing $X_{i}(t_{j})$ we observe $U_{ij}=X_{i}(t_{j})+\delta_{ij}$, where $\delta_{ij}$ are i.i.d. mean zero random errors with variance $\kappa^{2}$. There exists several methods to reconstruct the original curve $X_{i}(\cdot)$ from the observed curve with measurement error. Zhang and Chen (2007) proposed to use local polynomial kernel smoothing technique for individual function reconstructions. They showed that under appropriate conditions and suitable choice of bandwidth, the smoothed trajectories, $\hat{X_{i}}(\cdot)$ will estimate the true curves ${X_{i}(\cdot)}$ with negligible error. Thus we can use the reconstructed curves $\hat{X_{i}}(\cdot)$ and the effect of such a substitution is asymptotically negligible.

Case 3 : Sparse design with measurement error

More generally, we consider the case where functional data is observed in irregular and sparse grid of points and covariate is observed with measurement error. Here we have observed response {($Y_{i}(t_{ij}),t_{ij}),j=1,2,\ldots,m_{i}$} and observed covariate {($U(t_{ij}),t_{ij}),j=1,2,\ldots,m_{1i}$}. Again we assume that $U_{ij}=X_{i}(t_{ij})+\delta_{ij}$, where $\delta_{ij}$ are i.i.d. mean zero random errors with variance $\kappa^{2}$. In such a sparse sampling design, it is generally assumed (Kim et al. 2018) although the individual number of observations $m_{i}$ is small, $\bigcup_{i=1}^{n}\bigcup_{j=1}^{m_{i}}{t_{ij}}$ is dense in $\mathcal{T}=[0,1]$. Then reconstructing of the original curves from the observed sparse curves can be done by using FPCA methods of Yao et al. (2005) or Hall et al. (2006). The methods are based on estimating the mean and covariance functions using local linear smoothing, and subsequently estimating the eigenvalues and eigenfunctions from a spectral decomposition of estimated covariance matrix.

As mentioned earlier, Li et al. (2010) proved uniform convergence of the mean, eigenvalues and eigenfunctions for both dense and sparse design under suitable regularity conditions. For prediction of the scores Yao et al. (2005) introduced the PACE method which ensures the estimated scores asymptotically goes to BLUP of the original scores. Then these estimates can be put together using Karhunen-Lo$\grave{e}$ve expansion to get estimates $\hat{X_{i}}(\cdot)$ of the true curve $X_{i}(\cdot)$. We again use the PVE criterion to select the number of PC. So for sparse data observed on irregular grid and observed with measurement error, we employ FPCA to get $\hat{X_{i}}(t)$ and then use {$Y_{i}(t_{ij}),\hat{X_{i}}(t_{ij}),j=1,2,\ldots,m_{i}\}_{i=1}^{n}$ as our original data to perform our proposed test. Our simulations show that we are able to maintain correct type I error rate and obtain satisfactory power of our testing method using this strategy.

2.6 Extension to multiple covariates

Now we illustrate how our testing method can be applied in multiple covariate setting.

In this case the likelihood is as in (2.4) given by

[TABLE]

with the only difference being $\amsmathbb{V}=\amsmathbb{V}(\tau_{0},\tau_{1},\ldots,\tau_{p})=\mathbb{\Sigma}+\tau_{0}\amsmathbb{B}\amsmathbb{B}^{T}+\tau_{1}\amsmathbb{Z}_{1}\amsmathbb{Z}_{1}^{T}+\ldots+\tau_{p}\amsmathbb{Z}_{p}\amsmathbb{Z}_{p}^{T}$. So testing for the effect of $X_{p}(\cdot)$ here similarly reduces to testing for the variance component $\tau_{p}$, namely, we test $H_{0}:\tau_{p}=0$ vs $H_{1}:\tau_{p}>0$. The score function $S_{\tau_{p}}(\tau_{0},\tau_{1},\ldots,\tau_{p})$ and the information matrix $\amsmathbb{I}(\bm{\theta})$ {$\bm{\theta}=(\tau_{0},\tau_{1},\ldots,\tau_{p})^{T}$} are also modified accordingly taking into account the additional variance components, e.g., the information matrix $\amsmathbb{I}(\bm{\theta})$ now has to be partitioned as

[TABLE]

Then our testing method can applied as in Section 2.3 and the generalization of Theorem (2.1) and Theorem (2.2) to this multicovariate setting remains valid. We have given examples of application of our testing method in multicovariate setting, both in our simulation study in Section 3 and real data application in Section 4, where we have considered two covariate scenarios.

3 Simulation Study

3.1 Study design

In this section, we investigate the performance of our testing method via simulation study. We evaluate our test in terms of type I error rate and power. We also compare our method to the existing bootstrapped F-test proposed by Kim et al. (2018). To this end, we consider the following two scenarios.

Scenario A (single covariate): We generate data from the model,

[TABLE]

where $\beta_{0}(t)=1+2t+t^{2}$ and $\beta_{1}(t)=t/8$. The original covariate $X_{i}(\cdot)$ are i.i.d. copies of $X(\cdot)$, where $X(t)=a+b\hskip 5.69054pt\sqrt[]{2}sin(\pi t)+c\hskip 5.69054pt\sqrt[]{2}cos(\pi t)$, where $a\sim\mathcal{N}(0,1)$, $b\sim\mathcal{N}(0,.85^{2})$ and $c\sim\mathcal{N}(0,.70^{2})$ and they are independent. As discussed in Section 2, we assume that we observe $X_{i}(t)$ with measurement error, i.e., we observe $U_{i}(t)=X_{i}(t)+\delta$, where $\delta\sim\mathcal{N}(0,.6^{2})$. The error process $\epsilon_{i}(t)$ is generated as

[TABLE]

where $\xi_{i1}\stackrel{{\scriptstyle iid}}{{\sim}}\mathcal{N}(0,2)$ and $\xi_{i2}\stackrel{{\scriptstyle iid}}{{\sim}}\mathcal{N}(0,0.75^{2})$. We consider the following sampling designs:

• •

Dense design: Functional data are observed in $S$ for each subject, where $S$ is the set of $m=81$ equidistant time points in $\mathcal{T}=[0,1].$

• •

Sparse design: The response $Y_{i}(t)$ and noisy covariate $U_{i}(t)$ both are observed in random $m_{Y_{i}}$ and $m_{U_{i}}$ points in $S$ where $m_{Y_{i}}\stackrel{{\scriptstyle iid}}{{\sim}}Uniform\{20,21,\ldots,31\}$ and also $m_{U_{i}}\stackrel{{\scriptstyle iid}}{{\sim}}Uniform\{20,21,\ldots,31\}.$

We consider two sample sizes, $n=100$ and $300$ for comparison with the bootstrapped F test method. An additional simulation is done in this same set up for dense data with $m=20$, $n=39$ to illustrate the performance of the testing method for finite sample sizes, as considered in the gait analysis in later section of this paper.

Scenario B (Multiple Covariates):

We generate data from the model,

[TABLE]

where $\beta_{0}(t)=1+2t+t^{2}$, $\beta_{1}(t)=t/8$ and $\beta_{2}(t)=sin(\pi t)$. The original covariates $X_{ik}(\cdot)$ are i.i.d. copies of $X_{k}(\cdot)$, where $X_{k}(t)=a_{k}+b_{k}\hskip 5.69054pt\sqrt[]{2}sin(\pi t)+c_{k}\hskip 5.69054pt\sqrt[]{2}cos(\pi t)$, where $a_{k}\sim\mathcal{N}(0,(2^{-.5(k-1)})^{2})$, $b_{k}\sim\mathcal{N}(0,(0.85\times 2^{-.5(k-1)})^{2})$, and $c_{k}\sim\mathcal{N}(0,(0.70\times 2^{-.5(k-1)})^{2})$, and they are independent. The covariates $X_{k}(\cdot)$ are same as considered in Kim et al. (2018). We observe $X_{ik}(t)$ with measurement error, i.e., we observe $U_{ik}(t)=X_{ik}(t)+\delta_{k}$, where $\delta_{k}\sim\mathcal{N}(0,.6^{2})$. The error process $\epsilon_{i}(t)$ are generated as in Scenario A described above. Similar sparse and dense design settings and sample sizes $n\in\{100,300\}$ are considered. Here the main goal is to test $H_{0}:\beta_{2}(t)=0$ against the alternative $H_{1}:\beta_{2}(t)\neq 0$. For each of the scenarios we use 1000 generated data sets to asses type I error and power. To model the regression functions, 7 cubic B-splines with equally spaced knots are used for both simulation scenarios.

3.2 Simulation Results

Scenario A:

We first assess type I error of the test. We use nominal levels of $\alpha=5\%$ for $n=100$ and $n=300$. The results are displayed in Table 1. The estimated standard errors are also given.

                                                [Table 1 about here]

We observe that our test maintains nominal type I error rate. For dense design, the estimated type I error lies below the nominal level of $.05$. For sparse design, As sample size increases, the size performance of the test improves, which is expected as the proposed test is a large sample one and for the asymptotic convergence to the null distribution to hold we need larger sample size. The nominal level lies within two standard error limit of estimated type I error in this case.

Next, we study the power performance of our test for a fixed nominal level $\alpha=5\%$. To this end, we generate data in the simulation set up mentioned earlier with $\beta_{1}(t)=dt/8$. Then $d=0$ corresponds to the null hypothesis and $d>0$ plays the role of effect size and captures the departure from null hypothesis. We compare the power of our test to the bootstrapped-F test for $n\in\{100,300\}$, and for both dense and sparse sampling designs. For comparison of power with the bootstrapped method, we use the results from the simulation study conducted by Kim et al. (2018). Results of our study are displayed in Figure 2.

                                                [Figure 1 about here]

We note that across the sparse and dense design scenarios, the proposed score based test produces higher power than the bootstrapped-F test method. This is expected as our method is a likelihood based method and directly uses value of the score. We also observe that as the sample size increases or $d$ increases the power of our test converges to one across all the settings faster than the bootstrapped-F test. Similar results for the additional simulation study is provided in the Figure S1 of supplemental material.

Scenario B:

The type I error rates for this scenario are also displayed in Table 1. Nominal levels of $\alpha=5\%$ for $n=100$ and $n=300$ are considered. Again we observe the proposed test maintaining nominal type I error rate, and improvement in size performance with increasing sample size, particularly for sparse design. A similar evaluation of power performance is done as in Scenario A. Namely we generate data as in simulation Scenario B with $\beta_{2}^{*}(t)=d\beta_{2}(t)$, where $d$ is a constant. The power curve of our testing method for both dense and sparse settings for $n\in\{100,300\}$ are displayed in Figure 3. We observe the power converging to one, as sample size or $d$ increases, across all the sampling designs.

                                                [Figure 2 about here]

Our simulation results illustrate the proposed testing method is able to maintain the nominal type I error rate as well as capture the departure from null hypothesis successfully, even when data is observed sparsely and there are multiple covariates observed with measurement error.

We have used $7$ cubic B-splines basis functions to model the regression functions in our simulations. The choice of the number of basis in reality depends on the type of design (dense or sparse) and number of observed time points for each subjects. To illustrate the effect of number of basis on the performance of the test an additional simulation is carried out as in Scenario A for dense design and sample size $n=300$, with number of basis $k_{1}=7,10,12,14,16$. The result is illustrated in Figure S2 of supplemental material. The estimated type I error of the test maintain the nominal level for all the cases. Noticeably, we observe a marginal decrease in the power with a larger number of basis functions, particularly for smaller effect size, which might be attributed to more random effects in the model and the coefficient functions getting more rough with the increasing number of basis. Analogous comparison for Scenario B is provided in Figure S3 of supplemental material where similar trend can be observed.

4 Real Data Applications

Through simulations we have shown our proposed testing method is able to identify significant covariates under both dense and sparse sampling design, even when the covariates are observed with measurement error. Next we present two real data applications of our testing method to demonstrate its usefulness in identifying significant time varying covariates in practical problems. We first consider the study of gait deficiency which is a typical case of dense data with small measurement error, subsequently we also apply our method to a study of dietary calcium absorption, where the data is sparse and measurement error is relatively higher.

4.1 Gait Data

In this study, the goal is to understand how the joints in hip and knee interact during a gait cycle (Theologis 2009), (Ramsay and Silverman 2005). Here, there are longitudinal measurements of hip and knee angles taken on 39 children on 20 equispaced evaluation points in $\mathcal{T}=[0,1]$. Figure 4 displays the observed individual trajectories of the hip and knee angles.

                                                [Figure 3 about here]

The study of gait is important as it helps to identify issues causing pain, and also implement and evaluate treatments to correct abnormalities. As discussed earlier, one natural question to ask here is whether the knee angles (response) are at all associated with the hip angles (covariate). In our terminology, here $Y_{i}(t)$ is the knee angle. We assume the hip angles $X_{i}(t)$ are observed with measurement error, i.e, we observe $U_{i}(t_{ij})=X_{i}(t_{ij})+\delta_{ij}$, where $\delta_{ij}$ are assumed to be white noise. We use the functional linear concurrent model (1) in this paper to model the time varying effect of hip angles on knee angles. The FLCM can be see as an extension of classical regression model allowing the covariates to have dynamic effects. In the functional linear concurrent modeling setup, we are therefore interested in testing $H_{0}:\beta_{1}(t)=0$ against the alternative $H_{1}:\beta_{1}(t)\neq 0$, where $\beta_{1}(t)$ denotes the linear concurrent effect of hip angles on knee angles. We use the proposed score test method with 7 cubic B-splines to model the time varying effect $\beta_{1}(t)$. The p-value of the proposed test is calculated to be $<10^{-4}$. So we reject the null hypothesis and conclude that knee angle at any fixed time point is associated with the hip angle at the same time point. Our findings match with that of Kim et al. (2018), who found the association to be significant using a bootstrapped-F test.

The effect of the number of basis function on the test is investigated by varying the number of cubic B-splines basis functions between $5$ to $10$, in all cases the p-value of the proposed test was calculated to be $<10^{-2}$, further confirming existence of significant concurrent association between knee and hip angles.

As the sample size ($n=39$) for this data is small and our testing method is an asymptotic one, we further evaluate the performance of our proposed method using a simulation study that captures the feature of the gait data. This also enables us to see the power performance of our method. This is done in a similar way as that of Kim et al. (2018). We use a model that mimic the feature of the gait data, generate a large simulated data set and assess the power performance of our method on the simulated data. In particular, We generate the covariate $X_{i}(t)$ from a process with the mean and covariance functions that equal their estimated counterparts from the data using FPCA. We also estimate the parameters $\beta_{0}(t)$, $\beta_{1}(t)$ from the fit of our full model and estimate $\Sigma(s,t)$ by doing FPCA on the residuals obtained from the full model fit. Subsequently we generate observations using the model $Y_{i}(t)=\hat{\beta_{0}}(t)+d\{\hat{X_{i}}(t)\hat{\beta_{1}}(t)\}+\epsilon_{i}(t)$, where $d$ is a constant, and $\epsilon_{i}(\cdot){\sim}\mathcal{N}(0,\hat{\Sigma}(\cdot,\cdot))$. We simulate $n=300$ response curves $Y_{i}(t)$ from the above set up. In the above set up $d$ again plays the role of a parameter, which controls departure from null hypothesis. We perform a power analysis simulating 1000 such data sets from above scenario for various $d$. Figure 5 shows the power curve obtained from our analysis using the proposed testing method with 7 cubic B-spline basis functions.

                                                [Figure 4 about here]

For $d=0$, the type I error is $.062$, and nominal level $\alpha=.05$ is within its 2-standard errors limit. As $d$ increases the power gradually increases and ultimately goes to one when $d=0.1$, which ensures power is one at $d=1$. Similar power curves with different number of basis functions are also provided in Figure 5. For a small number basis functions (5) we observe slightly inflated type I error, for larger number of basis (9), a loss in power is noticed for small effect sizes. Again the best results are obtained using a moderate number of basis functions. The results further confirm our conclusion that knee angle made during a gait cycle is concurrently associated with hip angle.

4.2 Calcium Absorption Data

We consider the dietary calcium absorption study given in Davis (2002). In this study the subjects are a group of 188 patients. We have data on calcium absorption, dietary calcium intake and BMI of these patients at irregular intervals between 35 to 64 years of their ages. The number of repeated measurements for each patient is between 1 to 4. Figure S4 in supplemental material shows the individual curves of patients’ calcium absorption and calcium intake along their ages. We notice that this is a typical case of sparse data with relatively high measurement error. In this study, we are interested in finding whether calcium intake has any effect on calcium absorption in presence of BMI. We assume the covariates BMI $X_{i1}(t)$ and calcium intake $X_{i2}(t)$ are observed with measurement error. So our observed data is $U_{i}(t_{ij})=X_{i1}(t_{ij})+\delta_{ij}$ and $V_{i}(t_{ij})=X_{i2}(t_{ij})+\nu_{ij}$. Our response variable is calcium absorption $Y_{i}(t)$. Following the study in Kim et al. (2018) we assume a functional linear concurrent regression model, $Y_{i}(t)=\beta_{0}(t)+X_{i1}(t)\beta_{1}(t)+X_{i2}(t)\beta_{2}(t)+\epsilon_{i}(t)$, and want to test $H_{0}:\beta_{2}(t)=0$ against the alternative $H_{1}:\beta_{2}(t)\neq 0.$

First we use FPCA as discussed in Section 2.5 on noisy covariates to get their smooth counterparts at time points where we have $Y_{i}(t)$ available and then apply our one sided score test method for multiple covariates described as in Section 2.6. We use 7 cubic B-Splines to model $\beta_{0}(t)$, $\beta_{1}(t)$ and $\beta_{2}(t)$. The p-value of our test is calculated to be $<10^{-5}$. Thus we conclude in presence of BMI, calcium intake of patients has a significant effect on calcium absorption which again matches with the findings in Kim et al. (2018) and other studies of dietary calcium absorption.

5 Discussion and Future Work

In this article, we have proposed a likelihood based method for testing of hypothesis in functional linear concurrent regression. We have formulated the problem as a test for variance component and have used a one sided score test approach. We have established the asymptotic null distribution of our test statistic under some standard assumptions. Through simulations we have shown our proposed method maintains the nominal type I error rate and also yield higher power compared to the existing bootstrapped-F test, even when data is observed sparsely and with measurement error. We have successfully applied our method in finding significance of covariates in two real data application, namely the gait study and calcium absorption study. We note our method is a general one and can be applied for testing in longitudinal data setting too, where we have more flexibility in assuming parametric form of error covariance structure and estimating it consistently from the data.

We have considered approximating the true coefficient function using finite cubic B-spline basis expansion. We have used cubic B-spline basis with equispaced knots for computational tractability. The placement and choice of optimal number of knots is itself a challenge and open problem in functional data and there is no single consensus (Ramsay and Silverman 2005; Vsevolozhskaya et al. 2014). Even if the true function lies in infinite dimensional space we have shown the proposed test performs well with a moderate number of basis functions. The choice of number of basis functions is subjective and depends on the type of design and number of observed time points for each subject. We have demonstrated using a large number of basis functions leads to loss of power of the test and therefore using a moderate number of basis functions is recommended.

We have considered Gaussian distribution for the functional error in our model, if the distribution is non Gaussian, one can still use the score test statistic proposed in this paper and employ a subject level bootstrap for performing the test. Similarly it would be of interest to explore how the score test approach can be extended to generalized functional concurrent model such as logistic regression. In developing our test, we have used a random effects formulation of the problem arising from directly penalizing the coefficient functions. It is also plausible to use penalty on $r$-th ($r>0$) derivative of the coefficient functions. In this case the main challenge is to handle the singularity of the resulting covariance matrices, which can be addressed by using mixed effects model and subsequently testing for variance component using variants of likelihood ratio tests (Crainiceanu and Ruppert 2004). Another important work for future could be to prove the consistency results of the FPCA approximation methods used in this article. Further we would like to extend our testing method to nonparametric functional concurrent regression and more general function on function regression models and these are possible areas for future research that could explored based on our work.

Software

Software in the form of R code, together with the data set and complete documentation is available at GitHub (https://github.com/rahulfrodo/FLCM_Score).

Acknowledgements

We would like to thank Dr. Janet Kim and Dr. Ana-Maria Staicu for sharing the results of their method (Kim et al. 2018) with us.

Supplemental Material

Appendix S1, Appendix S2 referenced in Section 2.3, and Figures S1$-$S4 are available online with Supplemental Material.

Appendix A Proof of Theorem 2.1:

Suppose that the conditions (a) and (b) of Theorem 1 hold, i.e., the null hypothesis $H_{0}:\tau_{1}=0$ holds, $\bm{\theta}_{0}=(\tau_{0}^{*},0)$ is the true value of $\bm{\theta}$ and $\mathbb{\Sigma}$ is the true covariance matrix of residual vector $\mathbf{\mathcal{E}}$. Denote $\amsmathbb{V}=\amsmathbb{V}(\bm{\theta}_{0},\mathbb{\Sigma})$, $\mathbf{Y}\sim\mathcal{N}(0,\amsmathbb{V})$ under null. Then $\amsmathbb{W}=\amsmathbb{V}^{-1/2}\mathbf{Y}\sim\mathcal{N}(0,\amsmathbb{I})$. Now $S_{\tau_{1}}(\bm{\theta}_{0})=-1/2\{tr(\amsmathbb{Z}^{T}\amsmathbb{V}^{-1}\amsmathbb{Z})-(\amsmathbb{V}^{-1/2}\mathbf{Y})^{T}\amsmathbb{V}^{-1/2}\amsmathbb{Z}\amsmathbb{Z}^{T}\amsmathbb{V}^{-1/2}(\amsmathbb{V}^{-1/2}\mathbf{Y})\}$. Hence

[TABLE]

Now we use the spectral decomposition $\amsmathbb{V}^{-1/2}\amsmathbb{Z}\amsmathbb{Z}^{T}\amsmathbb{V}^{-1/2}/n=\amsmathbb{U}\amsmathbb{D}\amsmathbb{U}^{T}$. Then we have

[TABLE]

This follows from the fact $\amsmathbb{X}=\amsmathbb{U}^{T}\amsmathbb{W}\sim\mathcal{N}(0,\amsmathbb{I})$ as $\amsmathbb{U}$ is an orthogonal matrix and nonzero eigenvalues of $\amsmathbb{V}^{-1/2}\amsmathbb{Z}\amsmathbb{Z}^{T}\amsmathbb{V}^{-1/2}/n$ and $\amsmathbb{Z}^{T}\amsmathbb{V}^{-1}\amsmathbb{Z}/n$ are same. Therefore we have shown $\frac{S_{\tau_{1}}^{2}(\bm{\theta}_{0})}{\Lambda(\bm{\theta}_{0})}\overset{d}{=}\frac{\{1/2(\sum_{\ell=1}^{k_{1}}\lambda_{\ell}x_{\ell}^{2}-\sum_{\ell=1}^{k_{1}}\lambda_{\ell})\}^{2}}{\Lambda_{n}(\bm{\theta}_{0})}$ under null, which completes the proof of our Theorem 1, namely:

[TABLE]

where $x_{\ell}\stackrel{{\scriptstyle iid}}{{\sim}}\mathcal{N}(0,1)$ and $\lambda_{\ell}$ are eigenvalues of $\amsmathbb{Z}^{T}\amsmathbb{V}(\bm{\theta}_{0},\mathbb{\Sigma})^{-1}\amsmathbb{Z}/n$.