Fractals with point impact in functional linear regression

This paper investigates a linear regression model involving a scalar response $Y$ and a predictor given by the value of the trajectory of a continuous stochastic process $X=\{X(t)$, $t\in[0,1]\}$ at some unknown time point. Specifically, we consider the point impact linear regression model

and focus on the time point $\theta\in(0,1)$ as the target parameter of interest. The intercept $\alpha$ and the slope $\beta$ are scalars, and the error $\varepsilon$ is taken to be independent of $X$, having zero mean and finite variance $\sigma^{2}$. The complete trajectory of $X$ is assumed to be observed (at least on a fine enough grid that it makes no difference in terms of accuracy), even though the model itself only involves the value of $X$ at $\theta$, which represents a “sensitive” time point in terms of the relationship to the response. The main aim of the paper is to show that the precision of estimation of $\theta$ is driven by fractal behavior in $X$, and to develop valid inferential procedures that adapt to a broad range of such behavior. Our model could easily be extended in various ways, for example, to allow multiple sensitive time points or further covariates, but, for simplicity, we restrict attention to (1).

Our motivation for developing this type of model arises from genome-wide expression studies that measure the activity of numerous genes simultaneously. In these studies, it is of interest to locate genes showing activity that is associated with clinical outcomes. Emilsson et al. emilsson , for example, studied gene expression levels at over 24,000 loci in samples of adipose tissue to identify genes correlated with body mass index and other obesity-related outcomes. Gruvberger-Saal et al. gruvberger used gene expression profiles from the tumors of breast cancer patients to predict estrogen receptor protein concentration, an important prognostic marker for breast tumors; see also buness . In such studies, the gene expression profile across a chromosome can be regarded a functional predictor, and a gene associated with the clinical outcome is identified by its base pair position $\theta$ along the chromosome; see Figure 1. Our aim here is to develop a method of estimating a confidence interval for $\theta$, leading to the identification of chromosomal regions that are potentially useful for diagnosis and therapy. Although there is extensive statistical literature on gene expression data, it is almost exclusively concerned with multiple testing procedures for detecting differentially expressed genes; see, for example, dula , salas .

Gene expression profiles (as in Figure 1) clearly display fractal behavior, that is, self-similarity over a range of scales. Indeed, fractals often arise when spatiotemporal patterns at higher levels emerge from localized interactions and selection processes acting at lower levels, as with gene expression activity. Moreover, the recent discovery l-a that chromosomes are folded as “fractal globules,” which can easily unfold during gene activation, also helps explain the fractal appearance of gene expression profiles.

A basic stochastic model for fractal phenomena is provided by fractional Brownian motion (fBm) (see manness ), in which the so-called Hurst exponent $H\in[0,1]$ calibrates the scaling of the self-similarity and provides a natural measure of trajectory roughness. It featured prominently in the pioneering work of Benoît Mandelbrot, who stated (man , page 256) that fBm provides “the most manageable mathematical environment I can think of (for representing fractals).” For background on fBm from a statistical modeling point of view, see gn .

The key issue to be considered in this paper is how to construct a confidence interval for the true sensitive time point $\theta_{0}$ based on its least squares estimator $\hat{\theta}_{n}$, obtained by fitting model (1) from a sample of size $n$,

We show that, when $X$ is fBm, both the rate of convergence $r_{n}$ and limiting distribution of $\hat{\theta}_{n}$ depend on $H$. In addition, we construct bootstrap confidence intervals for $\theta_{0}$ that do not require knowledge of $H$. This facilitates applications (e.g., to gene expression data) in which the type of fractal behavior is not known in advance; the trajectory in Figure 1 has an estimated Hurst exponent of about 0.1, but it would be very difficult to estimate precisely using data in a small neighborhood of $\hat{\theta}_{n}$, so a bootstrap approach becomes crucial. We emphasize that nothing about the distribution of $X$ is used in the construction of the estimators or the bootstrap confidence intervals; the fBm assumption will only be utilized to study the large sample properties of these procedures. Moreover, our main results will make essential use of the fBm assumption only locally, that is, in a small neighborhood of $\theta_{0}$.

The point impact model (1) can be regarded as a simple working model that provides interpretable information about the influence of $X$ at a specific location (e.g., a genetic locus). Such information cannot be extracted using the standard functional linear regression model rs given by

where $f$ is a continuous function and $\alpha$ is an intercept, because the influence of $X(t)$ is spread continuously across $[0,1]$ and point-impact effects are excluded. In the gene expression context, if only a few genes are predictive of $Y$, then a model of the form (1) would be more suitable than (3), which does not allow $f$ to have infinite spikes. In general, however, a continuum of locations is likely to be involved (as well as point-impacts), so it is of interest to study the behavior of $\hat{\theta}_{n}$ in misspecified settings in which the data arise from combinations of (1) and (3).

Asymptotic results for the least squares estimator (2) in the correctly specified setting are presented in Section 2. In Section 3 it is shown that the residual bootstrap is consistent for the distribution of $\hat{\theta}_{n}$, leading to the construction of valid bootstrap confidence intervals without knowing $H$. The nonparametric bootstrap is shown to be inconsistent in the same setting. The effect of misspecification is discussed in Section 4. A two-sample problem version of the point impact model is discussed in Section 5. Some numerical examples are presented in Section 6, where we compare the proposed bootstrap confidence interval with Wald-type confidence intervals (in which $H$ is assumed to be known); an application to gene expression data is also discussed. Concluding remarks appear in Section 7. Proofs are placed in Section 8.

Throughout we take $X$ to be a fBm with Hurst exponent $H$, which, as discussed earlier, controls the roughness of the trajectories. We shall see in this section that the rate of convergence of $\hat{\theta}_{n}$ can be expressed explicitly in terms of $H$.

First we recall some basic properties of fBm. A (standard) fBm with Hurst exponent $H\in(0,1]$ is a Gaussian process $B_{H}=\{B_{H}(t)$, $t\in\mathbb{R}\}$ having continuous sample paths, mean zero and covariance function

By comparing their mean and covariance functions, $B_{H}(at)\stackrel{{\scriptstyle d}}{{=}}a^{H}B_{H}(t)$ as processes, for all $a>0$ (self-similarity). Clearly, $B_{1/2}$ is a two-sided Brownian motion, and $B_{1}$ is a random straight line: $B_{1}(t)=tZ$ where $Z\sim N(0,1)$. The increments are negatively correlated if $H<1/2$, and positively correlated if $H>1/2$. Increasing $H$ results in smoother sample paths.

Suppose $(X_{i},Y_{i}),i=1,\ldots,n$, are i.i.d. copies of $(X,Y)$ satisfying the model (1). The unknown parameter is $\eta=(\alpha,\beta,\theta)\in\Xi=\mathbb{R}^{2}\times[0,1]$, and its true value is denoted $\eta_{0}=(\alpha_{0},\beta_{0},\theta_{0})$. The following conditions are needed:

• [(A3)]

• (A1)

  $X$ is a fBm with Hurst exponent $H\in(0,1)$.

• (A2)

  $0<\theta_{0}<1$ and $\beta_{0}\neq 0$.

• (A3)

  $E|\varepsilon|^{2+\delta}<\infty$ for some $\delta>0$.

The construction of the least squares estimator $\hat{\eta}_{n}=(\hat{\alpha}_{n},\hat{\beta}_{n},\hat{\theta}_{n})$, defined by (2), does not involve any assumptions about the distribution of the trajectories, whereas the asymptotic behavior does. Our first result gives the consistency and asymptotic distribution of $\hat{\eta}_{n}$ under the above assumptions.

In general, a valid Wald-type confidence interval for $\theta_{0}$ would at least need a consistent estimator of the Hurst exponent $H$, which is a nuisance parameter in this problem. Unfortunately, however, accurate estimation of $H$ is difficult and quite often unstable. Bootstrap methods have been widely applied to avoid issues of nuisance parameter estimation, and they work well in problems with $\sqrt{n}$-rates; see bickel81 , singh81 , shao95 and the references therein. In this section we study the consistency properties of two bootstrap methods that arise naturally in our setting. One of these methods leads to a valid confidence interval for $\theta_{0}$ without requiring any knowledge of $H$.

The point impact model cannot capture effects that are spread out over the domain of the trajectory, for example, gene expression profiles for which the effect on a clinical outcome involves complex interactions between numerous genes. Such effects, however, may be represented by a functional linear model, and we now examine how the limiting behavior of $\hat{\theta}_{n}$ changes when the data arise in this way.

In this section we discuss a variation of the point impact regression model in which the response takes just two values (say $\pm 1$). This is of interest, for example, in case-control studies in which gene-expression data are available for a sample of cancer patients and a sample of healthy controls, and the target parameter is the locus of a differentially expressed gene.

Suppose we have two independent samples of trajectories $X$, with $n_{1}$ trajectories from class 1, and $n_{2}$ trajectories from class $-1$, for a total sample size of $n=n_{1}+n_{2}$. It is assumed that $\rho=n_{1}/n_{2}>0$ remains fixed, and the $j$th sample satisfies the model

where $\varepsilon_{ij}$, $i=1,\ldots,n_{j}$ are i.i.d. fBms with Hurst exponent $H\in(0,1)$, and $\mu_{j}(t)$ is an unknown mean function, assumed to be continuous. The treatment effect $\mathbb{M}(t)=\mu_{1}(t)-\mu_{2}(t)$ is taken to have a point impact in the sense of having a unique maximum at $\theta_{0}\in(0,1)$; minima can of course be treated in a similar fashion. The least squares estimator of the sensitive time point now becomes

where $\bar{X}_{j}(\theta)=\sum_{i=1}^{n_{j}}X_{ij}(\theta)/{n_{j}}$ is the sample mean for class $j$. Although a studentized version (normalizing the the difference of the sample means by a standard error estimate) might be preferable in some cases, with small or unbalanced samples, say, to keep the discussion simple, we restrict attention to $\hat{\theta}_{n}$. The empirical criterion function $\mathbb{M}_{n}(\theta)=\bar{X}_{1}(\theta)-\bar{X}_{2}(\theta)$ converges uniformly to $\mathbb{M}(\theta)$ a.s. (by the Glivenko–Cantelli theorem), so $\hat{\theta}_{n}$ is a consistent estimator of $\theta_{0}$.

As before, our objective is to find a confidence interval for $\theta_{0}$ based on $\hat{\theta}_{n}$ under appropriate conditions on the treatment effect. Toward this end, we need an assumption on the degree of smoothness of the treatment effect at $\theta_{0}$ in terms of an exponent $0<S\leq 1$:

as $\theta\to\theta_{0}$, where $c>0$. If $\mathbb{M}$ is twice-differentiable at $\theta_{0}$, then this assumption holds only with $S=1$; for it to hold for some $S<1$, a cusp is needed. When the smoothness of the treatment effect and the fBm match, that is, $S=H$, the rate of convergence of $\hat{\theta}_{n}$ is $n^{1/(2H)}$, as before, and $\hat{\theta}_{n}$ has a nondegenerate limit distribution of the same form as in Theorem 2.1:

The key step in the proof (which is simpler than in the regression case) is given at the end of Section 8.

In this section we report some numerical results based on trajectories from fBm simulations and from gene expression data.

We first consider a correctly specified example as in Section 2 and study the behavior of CIs for the sensitive time-point $\theta_{0}$ using the two bootstrap based methods, and compare them with the $100(1-\alpha)\%$ Wald-type CI

with $H$ assumed to be known. Here $\hat{\sigma}_{n}$ is the sample standard deviation of the residuals, and $z_{H,\alpha}$ is the upper $\alpha$-quantile of $\arg\min_{t\in\mathbb{R}}(B_{H}(t)+|t|^{2H}/2)$. In practice, $H$ needs to be estimated to apply (14). Numerous estimators of $H$ based on a single realization of $X$ have been proposed in the literature b , coeur , although observation at fine time scales is required for such estimators to work well, and it is not clear that direct plug-in would be satisfactory. The quantiles $z_{H,\alpha/2}$ needed to compute the Wald-type CIs were extracted from an extensive simulation of the limit distribution, as no closed form expression is available.

Table 1 reports the estimated coverage probabilities and average lengths of nominal 95% confidence intervals for $\theta_{0}$ calculated using 500 independent samples. The data were generated from the model (1), for $\alpha_{0}=0$, $\beta_{0}=1$, $\theta_{0}=1/2$, $\varepsilon\sim N(0,\sigma^{2})$ where $\sigma=0.3$ and $0.5$, the Hurst exponent $H=0.3,0.5,0.7$ and sample sizes $n=20$ and $40$. To calculate the least squares estimators (2), we restricted $\theta$ to a uniform grid of 101 points in $[0,1]$; the fBm trajectories were generated over the same grid. The fBm simulations were carried out in R, using the function fbmSim from the fArma package, and via the Cholesky method of decomposing the covariance matrix of $X$. Histograms and scatterplots of $\hat{\theta}_{n}$ and $\hat{\beta}_{n}$ for $H=0.3,0.5,0.7$ when $\sigma=0.5$ are displayed in Figure 2.

In practice, $X$ can only be observed at discrete time points, so restricting to a grid is the natural formulation for this example. Indeed, the resolution of the observation times in the neighborhood of $\theta_{0}$ is a limiting factor for the accuracy of $\hat{\theta}_{n}$, so the grid resolution needs to be fine enough for the statistical behavior of $\hat{\theta}_{n}$ to be apparent. For large sample sizes, a very fine grid would be needed in the case of a small Hurst exponent (cf. Theorem 2.1). Indeed, the histogram of $\hat{\theta}_{n}$ in the case $H=0.3$ (the first plot in Figure 2) shows that the resolution of the grid is almost too coarse to see the statistical variation, as the bin centered on $\theta_{0}=1/2$ contains almost 80% of the estimates. This phenomenon is also observed in Table 1 when $n=40$ and $\sigma=H=0.3$. The average length of the CIs is smaller than the resolution of the grid and, thus, we observe an over-coverage. The two histograms of $\hat{\theta}_{n}$ for $H=0.5$ and $H=0.7$, however, show increasing dispersion and become closer to bell-shaped as $H$ increases.

Recall that, for simplicity, we pretend as if we know $H$, which should be an advantage, yet the residual bootstrap has better performance based on the results in Table 1. We see that usually the Wald-type CIs have coverage less than the nominal 95%, whereas the inconsistent nonparametric bootstrap method over-covers with observed coverage probability close to 1. Accordingly, the average lengths of the Wald-type CIs are the smallest, whereas those obtained from the nonparametric bootstrap method are the widest. The behavior of CIs obtained from the nonparametric bootstrap method also illustrates the inconsistency of this procedure. A similar phenomenon was also observed in LM06 in connection with estimators that converge at $n^{1/3}$-rate.

Despite the asymptotic independence of $\hat{\theta}_{n}$ and $\hat{\beta}_{n}$, considerable correlation is apparent in the scatterplots in Figure 2, with increasing negative correlation as $H$ increases; note, however, that when $H=1$ there is a lack of identifiability of $\theta$ and $\beta$, so the trend in the correlation as $H$ approaches 1 is to be expected in small samples.

Next we consider a partially misspecified example, in which the data are now generated from (9) by setting $f(t)=1/2$ and $\theta=\theta_{*}=1/2$, but the other ingredients are unchanged from the correctly specified example. The plots in Figure 2 correspond to those in Figure 3. The effect of misspecification on $\hat{\theta}_{n}$ is a slight increase in dispersion but no change in mean; the effect on $\hat{\beta}_{n}$ is a substantial shift in mean along with a slight increase in dispersion.

In this paper we have developed a point impact functional linear regression model for use with trajectories as predictors of a continuous scalar response. It is expected that the proposed approach will be useful when there are sensitive time points at which the trajectory has an effect on the response. We have derived the rates of convergence and the explicit limiting distributions of the least squares estimator of such a parameter in terms of the Hurst exponent for fBm trajectories. We also established the validity of the residual bootstrap method for obtaining CIs for sensitive time points, avoiding the need to estimate the Hurst exponent. In addition, we have developed some results in the misspecified case in which the data are generated partially or completely from a standard functional linear model, and in the two-sample setting.

Although for simplicity of presentation we have assumed that the trajectories are fBm, it is clear from the proofs that it is only local properties of the trajectories in the neighborhood of the sensitive time point that drive the theory, and thus the validity of the confidence intervals. The consistency of the least squares estimator is of course needed, but this could be established under much weaker assumptions (namely, uniform convergence of the empirical criterion function and the existence of a well-separated minimum of the limiting criterion function; cf. vw , page 287).

Exploiting the fractal behavior of the trajectories plays a crucial role in developing confidence intervals based on the least squares estimator of the sensitive time point, in contrast to standard functional linear regression where it is customary to smooth the predictor trajectories prior to fitting the model (rs , Chapter 15). Our approach does not require any preprocessing of the trajectories involving a choice of smoothing parameters, nor any estimation of nuisance parameters (namely, the Hurst exponent). On the other hand, functional linear regression is designed with prediction in mind, rather than interpretability, so in a sense the two approaches are complimentary. The tendency of functional linear regression to over-smooth a point impact (see lmck for detailed discussion) is also due to the use of a roughness penalty on the regression function; the smoothing parameter is usually chosen by cross-validation, a criterion that optimizes for predictive performance.

Our model naturally extends to allow multiple sensitive time points, and any model selection procedure having the oracle property (such as the adaptive lasso) could be used to estimate the number of those sensitive time points. The bootstrap procedure for the (unregularized) least squares estimator can then be adapted to provide individual CIs around each time point, although developing theoretical justification would be challenging. Other challenging problems would be to develop bootstrap procedures that are suitable for the two-sample problem and for the misspecified model settings.

It should be feasible to carry through much of our program for certain types of diffusion processes driven by fBm, and also for processes having jumps. In the case of piecewise constant trajectories that have a single jump, the theory specializes to an existing type of change-point analysis kqs . Other possibilities include Lévy processes (which have stationary independent increments) and multi-parameter fBm. It should also be possible to develop versions of our results in the setting of censored survival data (e.g., Cox regression). Lindquist and McKeague lmck recently studied point impact generalized linear regression models in the case that $X$ is standard Brownian motion and we expect that our approach can be extended to such models as well.

To avoid measurability problems and for simplicity of notation, we will always use outer expectation/probability, and denote them by $E$ and $P$; $E^{*}$ and $P^{*}$ will denote bootstrap conditional expectation/probability given the data.

We begin with the proof of Theorem 2.1. The strategy is to establish (a) consistency, (b) the rate of convergence, (c) the weak convergence of a suitably localized version of the criterion function, and (d) apply the argmax (or argmin) continuous mapping theorem.