Statistical inference on partially shape-constrained function-on-scalar linear regression models

Let $Y_{i}$ be the functional response coupled with a vector covariate $\mathbf{X}_{i}=(X_{i,1},\ldots,X_{i,p})^{\top}$ as

independently for $i=1,\ldots,n$, where $\boldsymbol{\beta}(t)=\big{(}\beta_{1}(t),\ldots,\beta_{p}(t)\big{)}^{\top}$ is a vector coefficient function and $\varepsilon_{i}(t)$ is a mean-zero stochastic process that models the regression error uncorrelated with $\mathbf{X}_{i}$, i.e., $\mathbb{E}(\varepsilon_{i}(t)|\mathbf{X}_{i})=0$. We assume that $X_{i,1},\ldots,X_{i,p}$ are linearly independent with positive probability, which allows the design matrix to include the intercept, e.g., $X_{i,1}=1$ for all $i=1,\ldots,n$, depending on the inferential interest. Despite the generality of the functional model (2.1), it is practically infeasible to observe the functional outcomes as the infinite-dimensional objects, and we assume that discrete evaluations $\mathbf{Y}_{i}=(Y_{i,1},\ldots,Y_{i,L_{i}})^{\top}$ with $Y_{i,\ell}=Y_{i}(T_{i,\ell})$ are only accessible for $\ell=1,\ldots,L_{i}$, where $T_{i,1},\ldots,T_{i,L_{i}}$ is a random sample of $T$ with a probability density function $\pi$ on $[0,1]$. Here, $L_{i}$ is an independent random integer that may or may not depend on the sample size $n$. More detailed conditions on $L_{i}$ required for kernel and regression spline methods are provided in the Supplementary Material S.1 and S.2. We denote the finite random sample as $\mathcal{X}_{n}=\{(\mathbf{Y}_{i},\mathbf{X}_{i}):i=1,\ldots,n\}$.

Suppose we test the functional shape of $\beta_{j}$ on a pre-specified sub-interval $\mathcal{I}_{j}\subset[0,1]$ for $j\in J$ with $J\subset\{1,\ldots,p\}$. The null hypothesis “$H_{0}:\textrm{$H_{0,j}$ is true for all $j\in J$}$” consists of individual hypotheses

for $j\in J$. For example, if we hypothesize that $\beta_{j}$ is monotone increasing on $\mathcal{I}_{j}$, we read (2.2) as “$H_{0,j}:\textrm{$\beta_{j}$ is monotone increasing on $\mathcal{I}_{j}$}$.” Similarly, if one tests a partial convexity hypothesis, (2.2) can also be written as “$H_{0,j}:\textrm{$\beta_{j}$ is convex on $\mathcal{I}_{j}$}$.” In this paper, monotone and convex hypotheses will mainly be exemplified, yet one can extend our framework to general shape constraints, similarly as covered by many others in the literature [17, 18, 27, 28, 20, 25].

We then propose to reject the null hypothesis $H_{0}$ if the test statistic

is significantly large, where $\widehat{\boldsymbol{\beta}}(t)=\big{(}\hat{\beta}_{1}(t),\ldots,\hat{\beta}_{p}(t)\big{)}^{\top}$ and $\widetilde{\boldsymbol{\beta}}(t)=\big{(}\tilde{\beta}_{1}(t),\ldots,\tilde{\beta}_{p}(t)\big{)}^{\top}$ denote the estimates of coefficient function $\boldsymbol{\beta}(t)=\big{(}\beta_{1}(t),\ldots,\beta_{p}(t)\big{)}^{\top}$ under the null hypothesis $H_{0}$ and its general alternative $H_{1}$, respectively. While this $D_{n}$ is partly motivated by [22], a goodness-of-fit based test statistic for validating functional constraints on $\beta_{j}(t)$ via linear operator, similar types of the $L^{2}$-based test statistics were frequently employed in the literature for testing the nullity of functional difference for general scope [30, 36, 35]. Indeed, our proposed method corresponds to the global shape constraints if one sets $\mathcal{I}_{j}=[0,1]$. The following Sections 2.2 and 2.3 elaborate on how to obtain $\widehat{\boldsymbol{\beta}}(t)$ under kernel-smoothing and spline-smoothing approaches, respectively. Based on those estimators, we present a bootstrap test procedure in Section 2.4.

We assume that the vector coefficient $\boldsymbol{\beta}(t)$ is as smooth as it allows local linear approximation $\boldsymbol{\beta}(u)\approx\boldsymbol{\beta}(t)+\dot{\boldsymbol{\beta}}(t)(u-t)$ for $u$ near $t$, where $\dot{\boldsymbol{\beta}}(t)=\big{(}\beta_{1}^{\prime}(t),\ldots,\beta_{p}^{\prime}(t)\big{)}^{\top}$ denotes the gradient of $\boldsymbol{\beta}(t)$. The local linear kernel estimator of $\mathbb{B}(t)=\big{[}\boldsymbol{\beta}(t),\dot{\boldsymbol{\beta}}(t)\big{]}\in\mathbb{R}^{p\times 2}$ is given by the unconstrained minimizer $\widetilde{\mathbb{B}}(t)$ of the following objective function.

where $K_{h}(u)=K(u/h)/h$ with a probability density function $K$ and a bandwidth $h>0$, $\mathbf{Z}_{i,\ell}(t)=\big{(}1,T_{i,\ell}-t\big{)}^{\top}$ is the local smoothing design, and $\mathbf{X}_{i}^{\top}\mathbb{B}(t)\mathbf{Z}_{i,\ell}(t)=\mathbb{W}_{i,\ell}(t)^{\top}\mathrm{vec}\big{(}\mathbb{B}(t)\big{)}$ with $\mathrm{vec}\big{(}\mathbb{B}(t)\big{)}=\big{(}\boldsymbol{\beta}(t)^{\top},\dot{\boldsymbol{\beta}}(t)^{\top}\big{)}^{\top}\in\mathbb{R}^{2p}$ of $\mathbb{B}(t)$ and the Kronecker product $\mathbb{W}_{i,\ell}(t)=\mathbf{Z}_{i,\ell}(t)\otimes\mathbf{X}_{i}$ of the covarite design for local linear smoothing.

The implication of conditions $(K0)$–$(K4)$ are provided in Supplementary Material S.1. It is worth mentioning that Proposition 2.2.1 is valid even if only a few longitudinal observations are available for some subjects. For example, it is common in observational studies that dense observations of functional responses for some subjects may not be available regardless of the sample size $n$, i.e., the minimum number of longitudinal observations $\lambda_{n}=\min_{1\leq i\leq n}L_{i}$ is bounded. Therefore, the technical condition on $\lambda_{n}$ in $(K4)$ is not a restriction because we do not require $\lambda_{n}\to\infty$ as $n\to\infty$.

For the shape-constrained estimation of $\widehat{\boldsymbol{\beta}}$ under $H_{0}$, we extend the kernel-weighted least squares for nonparametric regression with scalar responses [34], where side conditions are empirically examined over a fine grid $\mathcal{G}_{M}=\{t_{m}\in[0,1]:m=0,1,\ldots,M\}$ with $0=t_{0}<t_{1}<\cdots<t_{M-1}<t_{M}=1$ so that $\hat{\beta}_{j}$ satisfies the null hypothesis (2.2) over $\mathcal{G}_{j}=\mathcal{G}|_{\mathcal{I}_{j}}=\{t_{j_{m}}\in\mathcal{I}_{j}:m=0,1,\ldots,M_{j}\}$ for each $j\in J$. Specifically, if we consider “$H_{0,j}:\textrm{$\beta_{j}$ is monotone increasing on $\mathcal{I}_{j}$}$,” we propose to substitute

Similarly, we empirically validate the partial hypothesis “$H_{0,j}:\textrm{$\beta_{k}$ is convex on $\mathcal{I}_{j}$}$” as

To define the shape-constrained kernel-weighted least squares for the function-on-scalar regression model (2.1) under the empirical null hypothesis on the grid $\mathcal{G}_{M}$, we note that

where $\tilde{e}_{i,\ell}=Y_{i,\ell}-\mathbb{W}_{i,\ell}^{\top}\mathrm{vec}\big{(}\widetilde{\mathbb{B}}(t)\big{)}$ is the unconstrained residual of fitting $Y_{i,\ell}$ and $\mathcal{Q}_{n}\big{(}\mathbb{B}(t),\widetilde{\mathbb{B}}(t)\big{)}=\mathrm{vec}\big{(}\mathbb{B}(t)-\widetilde{\mathbb{B}}(t)\big{)}^{\top}\widehat{\Psi}(t)\mathrm{vec}\big{(}\mathbb{B}(t)-\widetilde{\mathbb{B}}(t)\big{)}$ with $\widehat{\Psi}(t)=n^{-1}\sum_{i=1}^{n}L_{i}^{-1}\sum_{\ell=1}^{L_{i}}\mathbb{W}_{i,\ell}(t)\mathbb{W}_{i,\ell}(t)^{\top}K_{h}(T_{i,\ell}-t)$. To get (2.8), we used the fact that the unconstrained estimator $\widetilde{\mathbb{B}}(t)$ solves the estimating equation $\partial\mathcal{L}_{n}\big{(}\mathbb{B}(t)\big{)}/\partial\mathrm{vec}\big{(}\mathbb{B}(t)\big{)}=0_{2p}$, or equivalently $n^{-1}\sum_{i=1}^{n}L_{i}^{-1}\sum_{\ell=1}^{L_{i}}\mathbb{W}_{i,\ell}(t)\tilde{e}_{i,\ell}K_{h}(T_{i,\ell}-t)=0_{2p}$, where $0_{2p}=(0,\ldots,0)^{\top}$ denotes the zero vector in $\mathbb{R}^{2p}$. Noting that the minimization of $\mathcal{L}_{n}\big{(}\mathbb{B}(t)\big{)}$ only depends on the first term of (2.8) given the unconstrained estimator $\widetilde{B}(t)$, we propose the constrained optimization problem as follows.

where $0_{k}\in\mathbb{R}^{k}$ and $O_{k\times\ell}\in\mathbb{R}^{k\times\ell}$ are the $k$-vector and $(k,\ell)$-matrix of zeros, respectively, and $A_{j}\in\mathbb{R}^{M_{j}\times 2p(M+1)}$ is the linear constraint matrix associated with the individual hypothesis $H_{0,j}$. For the monotone increasing hypothesis (2.6), we set $A_{j}=A_{\mathrm{mono}}(\mathcal{I}_{j})\otimes\big{[}\mathrm{e}_{j}^{\top},0_{p}^{\top}\big{]}$, where $\textrm{e}_{j}\in\mathbb{R}^{p}$ is the unit vector with $1$ only at its $j$-th coordinate and $A_{\mathrm{mono}}(\mathcal{I}_{j})\in\mathbb{R}^{M_{j}\times(M_{j}+1)}$ is defined as

Similarly, for the convexity hypothesis (2.7), we set we set $A_{j}=A_{\mathrm{conv},0}(\mathcal{I}_{j})\otimes\big{[}\mathrm{e}_{j}^{\top},0_{p}^{\top}\big{]}-A_{\mathrm{conv},1}(\mathcal{I}_{j})\otimes\big{[}0_{p}^{\top},\mathrm{e}_{j}^{\top}\big{]}$, where $A_{\mathrm{conv},0}(\mathcal{I}_{k}),A_{\mathrm{conv},1}(\mathcal{I}_{j})\in\mathbb{R}^{M_{j}\times(M_{j}+1)}$ are defined as

and $A_{\mathrm{conv},1}(\mathcal{I}_{j})=\big{[}I_{M_{j}},0_{M_{j}}\big{]}$. This illustrates that the constrained estimates $\widehat{\mathbb{B}}(t_{0}),\ldots,\widehat{\mathbb{B}}(t_{M})$ can be obtained by the standard quadratic programming with $\sum_{j\in J}M_{j}$ linear constraints.

Theorem 2.2.2 indicates that the shape-constrained kernel least square estimator achieves the same rate of convergence as the unconstrained estimator under the null hypothesis. Moreover, suppose the empirical distribution of $\mathcal{G}_{M}$ converges to the uniform distribution on $[0,1]$. Then, it can be verified that $\int_{0}^{1}\big{\|}\widehat{\boldsymbol{\beta}}(t)-\boldsymbol{\beta}_{0}(t)\big{\|}\,\mathrm{d}t=O_{P}(n^{-2/5}\sqrt{\log n})$ similarly to Remark 1.

We now consider fitting the functional regression model by means of spline functions for smooth, flexible, and parsimonious estimation of coefficient functions. Let $\phi_{k}(t)$, $k=1,\ldots,d$ denote spline basis functions defined over $[0,1]$ for a given sequence of knots, then we express $\beta_{j}(t)=\sum_{k=1}^{d}c_{jk}\phi_{k}(t)$, where $c_{jk}$ denotes the basis coefficients. Under the shape constraints on $\beta_{j}(t)$ over $\mathcal{I}_{j}$, we can still employ the regression spline approach based on chosen basis functions, such as $I$-splines [23], $C$-splines [18] or general $B$-splines, with corresponding constraints assigned on a subset of basis coefficients $c_{jk}$. While we focus on partially constrained estimation based on $I$- or $C$- splines in this section, estimation via $B$-splines is also discussed in the Supplementary Material S.6.

[23] first introduced $I$-splines for curve estimation under the monotonicity restriction over the entire domain. As illustrated in the left panel of Figure 2, $I$-spline basis functions are piece-wise quadratic, and at each of the knots indicated by dotted vertical lines, there is exactly one basis function with a non-zero slope. Thus, a monotone increasing (decreasing) spline function can be constructed as a non-negative (non-positive) linear combination of $I$-spline basis functions. For readers interested in generating process of $I$-splines, we refer [23] on how piece-wise quadratic $I$-splines are formed from nonnegative $M$-spline family. Then, if the specific aim is on assigning, for example, monotone increasing restriction over a subset of the domain $\mathcal{I}_{j}$, non-negative coefficients condition should apply to a subset of $c_{jk}$ for which corresponding $I$-spline $\phi_{k}(t)$ display positive slopes within the range of $\mathcal{I}_{j}$.

Next, in terms of convexity constraints, [18] proposed $C$-splines by integrating $I$-splines. As shown in the right panel of Figure 2, $C$-splines form convex and piece-wise cubic functions with non-zero second derivatives at each knot, implying that a linear combination of $C$-splines with non-negative (non-positive) basis coefficients ensures the convexity (concavity) of the function estimator. See Section 2 of [18] for details. Similarly, for the convexity constraints over sub-interval(s) $\mathcal{I}_{j}$, we can impose non-negative restrictions to basis coefficients corresponding to splines showing positive second derivates features within $\mathcal{I}_{j}$.

We now apply the shape-restricted spline function technique to the functional regression model framework to estimate functional coefficients under (2.2). Although methodological or theoretical developments are general to other shape constraints, such as convexity or monotone convexity, we elaborate on the estimation under the monotonicity restriction below and refer to remarks for other restrictions. Let $\phi_{k}^{j}(t)$, $k=1,\ldots,d_{j}$, denote piece-wise quadratic $I$-spline basis functions under a given knot sequence, employed to fit $\beta_{j}(t)$, $j=1,\ldots,p$. Especially for partially constrained coefficients $j\in J$, the knot sequence should include lower and upper boundaries of $\mathcal{I}_{j}$. For the case of multiple disjoint sub-intervals form of $\mathcal{I}_{j}$, all boundaries should be a part of the knot sequence. We then estimate the functional regression coefficients by minimizing the following objective function with respect to non-negative conditions assigning on a subset of basis coefficients $c_{jk}$ depending on the set $J$ and $A_{j}$ as,

where $\phi_{0}^{j}(t)$ indicates the intercept function, i.e., $\phi_{0}^{j}(t)=1$, so that $c_{j0}$ represents the location of monotonicity functions of $\beta_{j}(t)$, the dimensionality of splines $d_{j}$ is determined by the number of knots chosen to approximate $\beta_{j}$, the set $J$ from (2.2) includes indices which coefficients are constrained in $H_{0}$, and the set $A_{j}$ is defined as $\big{\{}k\in\{1,\ldots,d_{j}\}:{\phi_{k}^{j}}{{}^{\prime}}(t)>0,~{}\textrm{for}~{}t\in\mathcal{I}_{j}\big{\}}$.

To obtain shape-restricted $\hat{\beta}_{j}(t)$, we first write our regression problem associated with (2.17) as below. We elaborate for the case of regular evaluation grids, $T_{1},\ldots,T_{L}$, for $i=1,\ldots,n$, and $d_{j}=d$ for representation simplicity, but having different $T_{i,\ell}$ or $d_{j}$ does not affect methodological developments. Let $\textrm{vec}\big{(}[\boldsymbol{c}_{j}]_{j=1}^{p}\big{)}=(\boldsymbol{c}_{1}^{\top},\ldots,\boldsymbol{c}_{p}^{\top})^{\top}$ for $\boldsymbol{c}_{j}=(c_{j0},c_{j1},\ldots,c_{jd})^{\top}$, then

under $c_{jk}\geq 0$ for certain $(j,k)$ specified in (2.17), where $\Phi$ denotes a $L\times(d+1)$ basis matrix having its $\ell$-th row consisting of $\big{(}1,\phi_{1}(T_{\ell}),$ $\ldots,$ $\phi_{d}(T_{\ell})\big{)}$ for $\ell=1,\ldots,L$. Then, the coefficient estimation is associated with finding the projection of $nL$-dimensional $(\mathbf{Y}_{1}^{\top}\ldots,\mathbf{Y}_{n}^{\top})^{\top}$ onto the space spanned by columns of $(\mathbf{X}\otimes\Phi)$ with corresponding constraints on $\boldsymbol{c}_{j}$, where $\mathbf{X}$ denotes the $(n\times p)$ design matrix. By letting $\boldsymbol{\sigma}^{h}\in\mathbb{R}^{nL}$, for $h=1,\ldots,p(d+1)$, denote the set of column vectors on $(\mathbf{X}\otimes\Phi)$ and $\mathcal{L}(\cdot)$ represent the linear space spanned by a given set of vectors, we separate $\boldsymbol{\sigma}^{h}$, $h\in H_{1}$, a set of column vectors associated with basis coefficients with no non-negativity constraint, and define $\mathcal{V}=\mathcal{L}\big{(}\{\boldsymbol{\sigma}^{h};~{}h\in H_{1}\}\big{)}$ indicating the linear space spanned by them. We then further obtain a set of generating vectors, denoted as $\boldsymbol{\delta}^{h}$ for $h\in H_{2}$ satisfying $H_{2}=\{1,\ldots,p(d+1)\}\backslash H_{1}$, that are orthogonal to $\mathcal{V}$. That is, $\boldsymbol{\delta}^{h}=\boldsymbol{\sigma}^{h}-\Pi(\boldsymbol{\sigma}^{h}|\mathcal{V})$, for $h\in H_{2}$, where $\Pi$ indicates a projection operator, such as the Gram-Schmidt process. By extending [18], we can characterize the constrained set as

Then the projection of $(\mathbf{Y}_{1}^{\top}\ldots,\mathbf{Y}_{n}^{\top})^{\top}$ to $\mathcal{C}$ can be fulfilled by projecting it onto the set of nonnegative linear combinations of $\boldsymbol{\delta}^{h}$, $h\in H_{2}$, called as constraint cone, and onto the $\mathcal{V}$, separately. We refer to Section 2 of [18] for readers interested in more comprehensive descriptions under a nonparametric regression setting. Owing to the fact that the constraint set $\mathcal{C}$ is a closed convex polyhedral cone in $\mathbb{R}^{nL}$, the projection is unique, and we employ the cone projection algorithm of [19] to find a solution, available as the function qprog or coneA in the R package coneproj. We note that the unconstrained estimates $\tilde{\beta}_{j}(t)$ can be calculated by the least square estimates of $\boldsymbol{c}_{j}$ under (2.18).

Next, we investigate rates of convergence for estimated functional coefficients under the constraints. Let $K_{jn}$ denote the number of knots to fit $\beta_{j}(t)$ growing with $n$, $q$ is the order of the spline, and $\|a\|_{L_{2}}$ be the $L_{2}$ norm of a square-integrable function $a(t)$. Under Conditions $(S1)$–$(S6)$, deferred in Supplementary Material S.2, and $\lim_{n\rightarrow\infty}K_{n}\log K_{n}/n=0$ for $K_{n}=\text{max}_{1\leq j\leq p}K_{jn}$, Theorem 2 of [14] and Theorem 6.25 of [26] imply the consistency of unconstrained $\tilde{\boldsymbol{\beta}}(t)=\big{(}\tilde{\beta}_{1}(t),\ldots,\tilde{\beta}_{p}(t)\big{)}^{\top}$ written as, $\|\tilde{\beta}_{j}-\beta_{j}\|^{2}_{L_{2}}=O_{p}(K_{n}n^{-1}+K_{n}^{-2q})$, $j=1,\ldots,p.$ Next, we show the consistency of the constrained estimators by adding the condition (S7) specified in Supplementary Material.

We employ a resampling method to assess the statistical significance of the observed $D_{n}$ given a random sample $\mathcal{X}_{n}$. We note that resampling functional observations , say $\widetilde{\mathbf{Y}}_{i}$ $=$ $\big{(}\widetilde{Y}_{i}(T_{i,1}),$ $\ldots,$ $\widetilde{Y}_{i}(T_{i,L_{i}})\big{)}^{\top}$, is not straightforward because $\widetilde{Y}_{i}(T_{i,1}),\ldots,\widetilde{Y}_{i}(T_{i,L_{i}})$ should inherit the joint distribution of ${Y}_{i}(T_{i,1}),\ldots,{Y}_{i}(T_{i,L_{i}})$ with within-subject heteroscedasticity. In literature, the wild bootstrap method [33, 13, 16] has been extensively investigated that can effectively handle the heteroscedasticity, and there have been many advances for handling dependent data [3, 7, 29, 11]. In this study, we use the wild bootstrap procedure as described in Algorithm 1, which can be regarded as the wild bootstrap for clustered data [3, 7].

Once $B$ bootstrap samples $\mathcal{X}_{n}^{(1)},\ldots,\mathcal{X}_{n}^{(B)}$ are independently generated, we calculate the shape-constrained and unconstrained estimates $\widehat{\boldsymbol{\beta}}^{(b)}$ and $\widetilde{\boldsymbol{\beta}}^{(b)}$, and calculate the test statistic $D_{n}^{(b)}$, for each $b=1,\ldots,B$. These bootstrap test statistics are compared to the original test statistic $D_{n}$ to determine the bootstrap $p$-value, $p_{n}=B^{-1}\sum_{b=1}^{B}\mathbb{I}(D_{n}^{(b)}>D_{n})$. Then, we reject the null hypothesis $H_{0}$ if $p_{n}<\alpha$, and retain $H_{0}$ otherwise at the significance level $\alpha\in(0,1)$. Our numerical study (Section 3) demonstrates that the proposed test procedure fulfills consistency in the sense that the type I error rate meets the significance level as specified and that the power of the test increases as sample size increases.