Scalar-on-function local linear regression and beyond

TL;DR

This paper advances nonparametric scalar-on-function regression by developing a local linear approach that outperforms local constant methods and enables consistent estimation of functional derivatives, with demonstrated effectiveness on simulated and real data.

Contribution

It introduces a local linear regression method for scalar-on-function problems using a projection approach, improving estimation accuracy and enabling derivative estimation.

Findings

Local linear regression outperforms local constant methods.

The estimator of the functional derivative is consistent.

Good finite sample properties demonstrated on simulated and real data.

Abstract

Regressing a scalar response on a random function is nowadays a common situation. In the nonparametric setting, this paper paves the way for making the local linear regression based on a projection approach a prominent method for solving this regression problem. Our asymptotic results demonstrate that the functional local linear regression outperforms its functional local constant counterpart. Beyond the estimation of the regression operator itself, the local linear regression is also a useful tool for predicting the functional derivative of the regression operator, a promising mathematical object on its own. The local linear estimator of the functional derivative is shown to be consistent. On simulated datasets we illustrate good finite sample properties of both proposed methods. On a real data example of a single-functional index model we indicate how the functional derivative of the regression operator provides an original and fast, widely applicable estimating method.

Tables

Table 1. Table 1: Number of times, out of 100 100 100 , that the dimension is correctly selected, and running times (in s.).

$n$	$\rho =0.05$	$\rho =0.1$	$\rho =0.2$	$\rho =0.4$	Timing
100	27	20	21	11	0.60
150	46	48	30	11	0.69
200	56	59	46	13	0.51
250	79	67	54	12	0.89
300	85	86	56	20	1.03
350	86	87	67	20	1.62
400	96	93	71	23	1.92
450	94	90	83	18	2.39
500	98	97	86	29	3.29

Table 2. Table 2: Average and standard deviation (in brackets) of O ​ R ​ M ​ S ​ E ​ P 𝑂 𝑅 𝑀 𝑆 𝐸 𝑃 ORMSEP with J = 4 𝐽 4 J=4 and n ​ s ​ r = 0.05 𝑛 𝑠 𝑟 0.05 nsr=0.05 .

	$n$	$\rho =0.05$	$\rho =0.1$	$\rho =0.2$	$\rho =0.4$
$ORMSE{P}_{reg}$	$100$	0.388  (0.058)	0.391  (0.055)	0.438  (0.060)	0.593  (0.065)
	$150$	0.271  (0.042)	0.288  (0.039)	0.354  (0.056)	0.518  (0.049)
	$200$	0.216  (0.032)	0.231  (0.030)	0.294  (0.033)	0.459  (0.040)
	$250$	0.176  (0.026)	0.195  (0.022)	0.260  (0.029)	0.431  (0.043)
	$300$	0.151  (0.018)	0.172  (0.022)	0.233  (0.022)	0.397  (0.036)
	$350$	0.135  (0.022)	0.151  (0.019)	0.218  (0.025)	0.369  (0.033)
	$400$	0.117  (0.014)	0.137  (0.017)	0.201  (0.022)	0.350  (0.030)
	$450$	0.107  (0.014)	0.125  (0.014)	0.184  (0.018)	0.333  (0.027)
	$500$	0.098  (0.011)	0.116  (0.012)	0.174  (0.018)	0.323  (0.029)
$ORMSE{P}_{deriv}$	$100$	0.374  (0.181)	0.399  (0.171)	0.420  (0.162)	0.604  (0.119)
	$150$	0.256  (0.155)	0.256  (0.150)	0.348  (0.166)	0.528  (0.116)
	$200$	0.207  (0.145)	0.202  (0.130)	0.257  (0.126)	0.476  (0.107)
	$250$	0.134  (0.128)	0.170  (0.129)	0.220  (0.119)	0.446  (0.095)
	$300$	0.104  (0.095)	0.113  (0.093)	0.205  (0.115)	0.404  (0.112)
	$350$	0.104  (0.126)	0.103  (0.083)	0.173  (0.107)	0.378  (0.101)
	$400$	0.064  (0.053)	0.086  (0.097)	0.159  (0.104)	0.352  (0.095)
	$450$	0.070  (0.096)	0.084  (0.077)	0.124  (0.088)	0.353  (0.095)
	$500$	0.052  (0.038)	0.061  (0.041)	0.116  (0.083)	0.326  (0.111)

Table 3. Table 3: Model (M3) with n ​ s ​ r = 0.05 𝑛 𝑠 𝑟 0.05 nsr=0.05 and ρ = 0.05 𝜌 0.05 \rho=0.05 .

		$a=0$	$a=0.25$	$a=0.5$	$a=0.75$	$a=1$
Reg.	L	0.001  (0.000)	0.014  (0.001)	0.111  (0.009)	0.531  (0.032)	1.008  (0.010)
	LC	0.048  (0.006)	0.051  (0.007)	0.070  (0.011)	0.159  (0.017)	0.265  (0.027)
	LL	0.026  (0.021)	0.026  (0.018)	0.035  (0.020)	0.061  (0.011)	0.099  (0.012)
	MY	0.028  (0.029)	0.033  (0.021)	0.071  (0.022)	0.165  (0.027)	0.273  (0.071)
Deriv.	L	1.852  (0.024)	1.058  (0.079)	1.016  (0.025)	1.011  (0.016)	1.004  (0.011)
	LL	0.102  (0.089)	8.679  (8.054)	0.989  (1.169)	0.092  (0.157)	0.055  (0.048)
	MY	0.397  (0.126)	3.266  (2.804)	0.573  (0.433)	0.257  (0.127)	0.216  (0.069)

Table 4. Table 4: Model (M3) with n ​ s ​ r = 0.4 𝑛 𝑠 𝑟 0.4 nsr=0.4 and ρ = 0.4 𝜌 0.4 \rho=0.4 .

		$a=0$	$a=0.25$	$a=0.5$	$a=0.75$	$a=1$
Reg.	L	0.008  (0.004)	0.021  (0.004)	0.119  (0.011)	0.542  (0.037)	1.014  (0.020)
	LC	0.153  (0.023)	0.158  (0.024)	0.206  (0.028)	0.403  (0.036)	0.629  (0.051)
	LL	0.086  (0.063)	0.081  (0.060)	0.120  (0.051)	0.252  (0.040)	0.393  (0.039)
	MY	0.055  (0.044)	0.058  (0.025)	0.103  (0.027)	0.245  (0.077)	0.369  (0.079)
Deriv.	L	1.851  (0.028)	1.201  (0.090)	1.034  (0.017)	1.012  (0.007)	1.005  (0.006)
	LL	0.849  (0.138)	7.561  (7.474)	1.419  (0.937)	0.572  (0.172)	0.461  (0.074)
	MY	0.565  (0.204)	7.702  (10.505)	1.022  (0.527)	3.782  (33.740)	0.307  (0.111)

Table 5. Table 5: Stability of the bootstrap bandwidth selection.

(a) $n=100$ and $nsr=0.4$	(b) $n=500$ and $nsr=0.05$
$B$ $h_{deriv}$ $ORMSE{P}_{deriv}$ 50 33.070  (7.395) 0.218  (0.050) 100 35.710  (9.401) 0.239  (0.075) 500 34.820  (6.967) 0.225  (0.048) 1000 34.310  (6.872) 0.217  (0.046)	$B$ $h_{deriv}$ $ORMSE{P}_{deriv}$ 50 33.680  (4.537) 0.043  (0.003) 100 33.800  (4.837) 0.043  (0.003) 500 32.720  (3.893) 0.044  (0.004) 1000 33.760  (4.209) 0.044  (0.003)

Table 6. Table 6: Estimated bandwidths.

$n$	$nsr$	$h_{reg}$	$h_{deriv}$	$h_{deriv}^{oracle}$
100	0.05	11.150  (1.755)	15.060  (5.007)	18.610  (2.357)
	0.2	14.010  (2.560)	26.560  (6.191)	22.990  (2.834)
	0.4	17.910  (4.147)	35.700  (6.882)	26.490  (3.680)
150	0.05	11.330  (1.551)	17.240  (4.490)	21.810  (2.485)
	0.2	14.170  (2.396)	30.810  (5.144)	28.110  (2.998)
	0.4	17.160  (3.369)	37.500  (6.308)	33.540  (3.471)
200	0.05	11.560  (1.321)	19.460  (4.988)	24.800  (2.361)
	0.2	14.960  (2.247)	34.840  (5.158)	33.400  (3.260)
	0.4	18.600  (2.878)	41.500  (5.402)	39.440  (4.246)
250	0.05	11.460  (1.234)	21.840  (5.502)	27.620  (2.440)
	0.2	15.760  (2.151)	37.260  (5.010)	38.080  (3.617)
	0.4	18.920  (3.080)	45.280  (5.760)	45.320  (4.387)
300	0.05	11.870  (1.368)	24.200  (4.767)	30.590  (2.575)
	0.2	16.490  (1.957)	40.880  (4.326)	41.810  (3.786)
	0.4	19.820  (2.823)	49.220  (4.683)	51.650  (4.961)
350	0.05	12.140  (1.524)	26.360  (5.262)	32.930  (2.786)
	0.2	16.220  (2.116)	42.980  (5.067)	46.910  (3.728)
	0.4	20.870  (2.707)	52.730  (4.750)	56.090  (5.071)
400	0.05	12.360  (1.150)	28.080  (4.618)	35.320  (3.011)
	0.2	17.360  (2.067)	46.680  (4.759)	49.640  (3.509)
	0.4	21.360  (3.227)	55.960  (5.895)	60.640  (4.394)
450	0.05	13.000  (0.000)	31.240  (5.190)	37.600  (2.566)
	0.2	17.760  (1.944)	48.320  (4.479)	53.400  (4.566)
	0.4	22.920  (3.240)	59.880  (6.125)	65.000  (5.714)
500	0.05	13.520  (1.306)	33.840  (4.720)	40.200  (3.028)
	0.2	18.120  (2.006)	51.120  (4.959)	56.800  (4.158)
	0.4	23.280  (3.358)	61.800  (5.222)	69.920  (6.016)

Table 7. Table 7: Number of times, out of 100 100 100 , that the dimension is correctly selected with J = 2 𝐽 2 J=2 and n ​ s ​ r = 0.05 𝑛 𝑠 𝑟 0.05 nsr=0.05 .

$n$	$\rho =0.05$	$\rho =0.1$	$\rho =0.2$	$\rho =0.4$
100	99	100	100	99
150	100	98	100	99
200	100	100	100	100
250	100	100	100	99
300	100	100	100	100
350	100	100	100	99
400	100	100	100	100
450	100	100	100	100
500	100	99	100	100

Table 8. Table 8: Average and standard deviation (in brackets) of O ​ R ​ M ​ S ​ E ​ P 𝑂 𝑅 𝑀 𝑆 𝐸 𝑃 ORMSEP with J = 2 𝐽 2 J=2 and n ​ s ​ r = 0.05 𝑛 𝑠 𝑟 0.05 nsr=0.05 .

	$n$	$\rho =0.05$	$\rho =0.1$	$\rho =0.2$	$\rho =0.4$
$ORMSE{P}_{reg}$	$100$	0.031  (0.008)	0.038  (0.009)	0.059  (0.012)	0.116  (0.021)
	$150$	0.021  (0.006)	0.029  (0.006)	0.040  (0.006)	0.080  (0.014)
	$200$	0.017  (0.004)	0.023  (0.004)	0.035  (0.006)	0.066  (0.010)
	$250$	0.014  (0.003)	0.020  (0.005)	0.029  (0.005)	0.056  (0.008)
	$300$	0.012  (0.003)	0.017  (0.003)	0.026  (0.004)	0.047  (0.006)
	$350$	0.011  (0.003)	0.015  (0.003)	0.023  (0.004)	0.043  (0.006)
	$400$	0.010  (0.002)	0.014  (0.003)	0.021  (0.004)	0.038  (0.005)
	$450$	0.009  (0.002)	0.013  (0.002)	0.020  (0.003)	0.035  (0.005)
	$500$	0.009  (0.002)	0.012  (0.002)	0.019  (0.003)	0.033  (0.004)
$ORMSE{P}_{deriv}$	$100$	0.033  (0.011)	0.037  (0.008)	0.050  (0.010)	0.108  (0.023)
	$150$	0.024  (0.005)	0.028  (0.006)	0.037  (0.007)	0.080  (0.017)
	$200$	0.020  (0.004)	0.023  (0.004)	0.031  (0.005)	0.067  (0.013)
	$250$	0.018  (0.003)	0.020  (0.003)	0.027  (0.004)	0.063  (0.015)
	$300$	0.016  (0.003)	0.018  (0.003)	0.024  (0.003)	0.052  (0.013)
	$350$	0.015  (0.004)	0.016  (0.003)	0.022  (0.003)	0.050  (0.014)
	$400$	0.014  (0.002)	0.016  (0.003)	0.021  (0.003)	0.045  (0.010)
	$450$	0.013  (0.003)	0.015  (0.003)	0.020  (0.002)	0.042  (0.011)
	$500$	0.012  (0.002)	0.014  (0.002)	0.019  (0.002)	0.039  (0.009)

Table 9. Table 9: Number of times, out of 100 100 100 , that the dimension is correctly selected with J = 3 𝐽 3 J=3 and n ​ s ​ r = 0.05 𝑛 𝑠 𝑟 0.05 nsr=0.05 .

$n$	$\rho =0.05$	$\rho =0.1$	$\rho =0.2$	$\rho =0.4$
100	96	100	91	82
150	100	98	100	90
200	100	100	100	91
250	99	100	100	99
300	99	100	100	100
350	100	100	100	100
400	100	100	100	100
450	100	100	100	100
500	100	100	100	100

Table 10. Table 10: Average and standard deviation (in brackets) of O ​ R ​ M ​ S ​ E ​ P 𝑂 𝑅 𝑀 𝑆 𝐸 𝑃 ORMSEP with J = 3 𝐽 3 J=3 and n ​ s ​ r = 0.05 𝑛 𝑠 𝑟 0.05 nsr=0.05 .

	$n$	$\rho =0.05$	$\rho =0.1$	$\rho =0.2$	$\rho =0.4$
$ORMSE{P}_{reg}$	$100$	0.137  (0.038)	0.142  (0.027)	0.208  (0.043)	0.369  (0.056)
	$150$	0.086  (0.015)	0.100  (0.018)	0.147  (0.020)	0.285  (0.044)
	$200$	0.064  (0.009)	0.077  (0.012)	0.117  (0.015)	0.242  (0.036)
	$250$	0.054  (0.014)	0.066  (0.010)	0.102  (0.013)	0.207  (0.023)
	$300$	0.046  (0.007)	0.058  (0.008)	0.091  (0.010)	0.181  (0.021)
	$350$	0.040  (0.007)	0.053  (0.007)	0.084  (0.010)	0.169  (0.020)
	$400$	0.036  (0.005)	0.048  (0.007)	0.078  (0.009)	0.157  (0.016)
	$450$	0.033  (0.004)	0.044  (0.005)	0.073  (0.009)	0.143  (0.015)
	$500$	0.030  (0.004)	0.041  (0.005)	0.068  (0.008)	0.134  (0.013)
$ORMSE{P}_{deriv}$	$100$	0.097  (0.102)	0.087  (0.014)	0.153  (0.129)	0.300  (0.121)
	$150$	0.056  (0.009)	0.075  (0.074)	0.086  (0.013)	0.218  (0.097)
	$200$	0.045  (0.006)	0.052  (0.009)	0.074  (0.010)	0.185  (0.085)
	$250$	0.046  (0.073)	0.044  (0.006)	0.063  (0.008)	0.139  (0.033)
	$300$	0.042  (0.072)	0.040  (0.004)	0.057  (0.007)	0.125  (0.017)
	$350$	0.032  (0.003)	0.036  (0.003)	0.051  (0.007)	0.112  (0.013)
	$400$	0.029  (0.003)	0.033  (0.003)	0.048  (0.005)	0.106  (0.014)
	$450$	0.027  (0.002)	0.031  (0.003)	0.044  (0.005)	0.099  (0.015)
	$500$	0.025  (0.002)	0.029  (0.002)	0.041  (0.003)	0.093  (0.011)

Table 11. Table 11: Number of times, out of 100 100 100 , that the dimension is correctly selected with J = 2 𝐽 2 J=2 and n ​ s ​ r = 0.4 𝑛 𝑠 𝑟 0.4 nsr=0.4 .

$n$	$\rho =0.05$	$\rho =0.1$	$\rho =0.2$	$\rho =0.4$
100	91	97	98	91
150	96	95	93	98
200	94	97	97	98
250	99	99	100	98
300	99	100	98	99
350	99	97	100	99
400	98	99	100	100
450	99	100	100	100
500	97	100	100	99

Table 12. Table 12: Average and standard deviation (in brackets) of O ​ R ​ M ​ S ​ E ​ P 𝑂 𝑅 𝑀 𝑆 𝐸 𝑃 ORMSEP with J = 2 𝐽 2 J=2 and n ​ s ​ r = 0.4 𝑛 𝑠 𝑟 0.4 nsr=0.4 .

	$n$	$\rho =0.05$	$\rho =0.1$	$\rho =0.2$	$\rho =0.4$
$ORMSE{P}_{reg}$	$100$	0.107  (0.043)	0.111  (0.041)	0.143  (0.040)	0.235  (0.070)
	$150$	0.076  (0.026)	0.085  (0.030)	0.114  (0.040)	0.172  (0.042)
	$200$	0.063  (0.024)	0.067  (0.023)	0.095  (0.024)	0.143  (0.030)
	$250$	0.049  (0.017)	0.056  (0.017)	0.079  (0.021)	0.125  (0.023)
	$300$	0.042  (0.013)	0.050  (0.013)	0.072  (0.015)	0.110  (0.021)
	$350$	0.039  (0.012)	0.047  (0.013)	0.064  (0.014)	0.103  (0.019)
	$400$	0.036  (0.013)	0.045  (0.012)	0.058  (0.011)	0.095  (0.016)
	$450$	0.032  (0.008)	0.041  (0.009)	0.055  (0.011)	0.085  (0.013)
	$500$	0.030  (0.009)	0.038  (0.009)	0.052  (0.011)	0.082  (0.016)
$ORMSE{P}_{deriv}$	$100$	0.116  (0.117)	0.094  (0.072)	0.111  (0.070)	0.197  (0.089)
	$150$	0.072  (0.061)	0.077  (0.072)	0.094  (0.091)	0.146  (0.068)
	$200$	0.071  (0.086)	0.058  (0.056)	0.064  (0.017)	0.121  (0.027)
	$250$	0.045  (0.014)	0.050  (0.052)	0.055  (0.013)	0.101  (0.021)
	$300$	0.040  (0.011)	0.041  (0.009)	0.049  (0.011)	0.091  (0.016)
	$350$	0.038  (0.012)	0.044  (0.054)	0.044  (0.008)	0.081  (0.015)
	$400$	0.035  (0.015)	0.036  (0.009)	0.042  (0.008)	0.077  (0.014)
	$450$	0.034  (0.014)	0.032  (0.007)	0.040  (0.008)	0.071  (0.012)
	$500$	0.035  (0.054)	0.031  (0.007)	0.038  (0.006)	0.066  (0.013)

Table 13. Table 13: Number of times, out of 100 100 100 , that the dimension is correctly selected with J = 3 𝐽 3 J=3 and n ​ s ​ r = 0.4 𝑛 𝑠 𝑟 0.4 nsr=0.4 .

$n$	$\rho =0.05$	$\rho =0.1$	$\rho =0.2$	$\rho =0.4$
100	55	58	60	37
150	67	71	65	51
200	69	76	71	57
250	81	75	71	68
300	86	88	87	67
350	89	89	85	73
400	87	90	78	78
450	94	93	85	81
500	98	94	92	82

Table 14. Table 14: Average and standard deviation (in brackets) of O ​ R ​ M ​ S ​ E ​ P 𝑂 𝑅 𝑀 𝑆 𝐸 𝑃 ORMSEP with J = 3 𝐽 3 J=3 and n ​ s ​ r = 0.4 𝑛 𝑠 𝑟 0.4 nsr=0.4 .

	$n$	$\rho =0.05$	$\rho =0.1$	$\rho =0.2$	$\rho =0.4$
$ORMSE{P}_{reg}$	$100$	0.287  (0.085)	0.310  (0.076)	0.344  (0.075)	0.503  (0.077)
	$150$	0.202  (0.054)	0.223  (0.052)	0.276  (0.054)	0.426  (0.062)
	$200$	0.168  (0.042)	0.189  (0.040)	0.235  (0.050)	0.368  (0.048)
	$250$	0.141  (0.041)	0.163  (0.040)	0.205  (0.037)	0.324  (0.047)
	$300$	0.120  (0.025)	0.134  (0.030)	0.182  (0.032)	0.298  (0.043)
	$350$	0.108  (0.024)	0.123  (0.023)	0.167  (0.027)	0.274  (0.040)
	$400$	0.099  (0.023)	0.113  (0.023)	0.157  (0.028)	0.252  (0.030)
	$450$	0.089  (0.018)	0.104  (0.019)	0.147  (0.024)	0.239  (0.032)
	$500$	0.083  (0.015)	0.098  (0.018)	0.138  (0.018)	0.232  (0.029)
$ORMSE{P}_{deriv}$	$100$	0.321  (0.227)	0.300  (0.204)	0.315  (0.179)	0.513  (0.142)
	$150$	0.236  (0.212)	0.207  (0.157)	0.267  (0.187)	0.422  (0.160)
	$200$	0.229  (0.234)	0.179  (0.164)	0.217  (0.171)	0.351  (0.140)
	$250$	0.159  (0.186)	0.182  (0.190)	0.200  (0.161)	0.287  (0.133)
	$300$	0.130  (0.151)	0.121  (0.136)	0.133  (0.099)	0.272  (0.134)
	$350$	0.103  (0.112)	0.106  (0.114)	0.130  (0.107)	0.244  (0.131)
	$400$	0.121  (0.170)	0.102  (0.123)	0.154  (0.138)	0.217  (0.116)
	$450$	0.085  (0.123)	0.083  (0.082)	0.120  (0.110)	0.205  (0.123)
	$500$	0.063  (0.073)	0.076  (0.076)	0.093  (0.078)	0.192  (0.112)

Table 15. Table 15: Number of times, out of 100 100 100 , that the dimension is correctly selected with J = 4 𝐽 4 J=4 and n ​ s ​ r = 0.4 𝑛 𝑠 𝑟 0.4 nsr=0.4 .

$n$	$\rho =0.05$	$\rho =0.1$	$\rho =0.2$	$\rho =0.4$
100	6	3	6	9
150	8	7	5	3
200	7	6	7	3
250	9	8	8	4
300	15	12	7	3
350	17	10	11	1
400	15	7	5	2
450	24	12	10	0
500	19	16	12	1

Table 16. Table 16: Average and standard deviation (in brackets) of O ​ R ​ M ​ S ​ E ​ P 𝑂 𝑅 𝑀 𝑆 𝐸 𝑃 ORMSEP with J = 4 𝐽 4 J=4 and n ​ s ​ r = 0.4 𝑛 𝑠 𝑟 0.4 nsr=0.4 .

	$n$	$\rho =0.05$	$\rho =0.1$	$\rho =0.2$	$\rho =0.4$
$ORMSE{P}_{reg}$	$100$	0.516  (0.101)	0.502  (0.075)	0.549  (0.081)	0.701  (0.087)
	$150$	0.394  (0.066)	0.416  (0.071)	0.451  (0.059)	0.608  (0.064)
	$200$	0.337  (0.057)	0.345  (0.054)	0.405  (0.051)	0.550  (0.055)
	$250$	0.298  (0.051)	0.310  (0.048)	0.351  (0.045)	0.512  (0.052)
	$300$	0.258  (0.037)	0.270  (0.044)	0.329  (0.048)	0.480  (0.052)
	$350$	0.239  (0.031)	0.250  (0.043)	0.305  (0.041)	0.453  (0.046)
	$400$	0.211  (0.032)	0.229  (0.032)	0.285  (0.033)	0.428  (0.037)
	$450$	0.199  (0.029)	0.212  (0.024)	0.271  (0.031)	0.409  (0.038)
	$500$	0.190  (0.032)	0.205  (0.024)	0.259  (0.032)	0.396  (0.033)
$ORMSE{P}_{deriv}$	$100$	0.561  (0.176)	0.577  (0.156)	0.559  (0.135)	0.707  (0.114)
	$150$	0.479  (0.180)	0.505  (0.169)	0.508  (0.137)	0.646  (0.111)
	$200$	0.451  (0.183)	0.446  (0.152)	0.478  (0.161)	0.602  (0.104)
	$250$	0.435  (0.207)	0.413  (0.158)	0.401  (0.122)	0.542  (0.100)
	$300$	0.354  (0.151)	0.363  (0.143)	0.403  (0.133)	0.529  (0.107)
	$350$	0.333  (0.157)	0.353  (0.143)	0.362  (0.117)	0.505  (0.091)
	$400$	0.350  (0.178)	0.367  (0.145)	0.368  (0.104)	0.479  (0.085)
	$450$	0.294  (0.158)	0.318  (0.114)	0.356  (0.117)	0.484  (0.087)
	$500$	0.319  (0.171)	0.295  (0.117)	0.319  (0.089)	0.462  (0.081)

Table 17. Table 17: Model (M3) with n ​ s ​ r = 0.1 𝑛 𝑠 𝑟 0.1 nsr=0.1 and ρ = 0.1 𝜌 0.1 \rho=0.1 .

		$a=0$	$a=0.25$	$a=0.5$	$a=0.75$	$a=1$
Reg.	L	0.001  (0.001)	0.015  (0.001)	0.113  (0.009)	0.532  (0.031)	1.007  (0.011)
	LC	0.061  (0.008)	0.064  (0.008)	0.087  (0.011)	0.187  (0.021)	0.292  (0.032)
	LL	0.035  (0.028)	0.034  (0.024)	0.048  (0.021)	0.091  (0.020)	0.142  (0.018)
	MY	0.026  (0.022)	0.035  (0.024)	0.724  (6.494)	0.207  (0.291)	0.271  (0.051)
Deriv.	L	1.847  (0.029)	1.093  (0.096)	1.026  (0.031)	1.014  (0.016)	1.004  (0.007)
	LL	0.228  (0.213)	8.999  (8.559)	1.191  (1.158)	0.210  (0.243)	0.124  (0.116)
	MY	0.407  (0.116)	5.040  (6.388)	0.766  (0.577)	0.287  (0.110)	0.221  (0.062)

Table 18. Table 18: Model (M3) with n ​ s ​ r = 0.2 𝑛 𝑠 𝑟 0.2 nsr=0.2 and ρ = 0.2 𝜌 0.2 \rho=0.2 .

		$a=0$	$a=0.25$	$a=0.5$	$a=0.75$	$a=1$
Reg.	L	0.003  (0.002)	0.017  (0.002)	0.113  (0.010)	0.534  (0.032)	1.013  (0.018)
	LC	0.085  (0.014)	0.092  (0.013)	0.125  (0.017)	0.260  (0.031)	0.402  (0.040)
	LL	0.036  (0.036)	0.050  (0.036)	0.071  (0.030)	0.141  (0.020)	0.220  (0.023)
	MY	0.062  (0.210)	0.051  (0.028)	0.161  (0.769)	0.191  (0.039)	0.653  (3.516)
Deriv.	L	1.852  (0.031)	1.133  (0.115)	1.024  (0.021)	1.011  (0.011)	1.005  (0.008)
	LL	0.514  (0.200)	9.041  (8.499)	1.203  (1.066)	0.333  (0.193)	0.279  (0.109)
	MY	0.602  (0.893)	12.163  (24.987)	1.125  (0.857)	0.505  (1.284)	0.326  (0.479)

Table 19. Table 19: Model (M3) with n ​ s ​ r = 0.05 𝑛 𝑠 𝑟 0.05 nsr=0.05 and ρ = 0.05 𝜌 0.05 \rho=0.05 .

		$a=0$	$a=0.25$	$a=0.5$	$a=0.75$	$a=1$
Reg.	L	0.003  (0.002)	0.017  (0.003)	0.119  (0.010)	0.567  (0.041)	1.029  (0.043)
	LC	0.134  (0.030)	0.138  (0.028)	0.188  (0.040)	0.376  (0.053)	0.593  (0.088)
	LL	0.071  (0.061)	0.084  (0.061)	0.110  (0.062)	0.239  (0.045)	0.382  (0.064)
	MY	0.446  (3.538)	65.192  (651.124)	0.236  (0.593)	0.814  (4.602)	1.609  (11.202)
Deriv.	L	1.865  (0.101)	1.229  (0.354)	1.069  (0.133)	1.045  (0.084)	1.012  (0.049)
	LL	0.760  (0.376)	10.463  (8.701)	1.424  (1.113)	0.472  (0.292)	0.413  (0.204)
	MY	0.680  (0.668)	50.611  (277.149)	5.677  (27.481)	3.090  (11.623)	0.369  (0.359)

Table 20. Table 20: Model (M3) with n ​ s ​ r = 0.1 𝑛 𝑠 𝑟 0.1 nsr=0.1 and ρ = 0.1 𝜌 0.1 \rho=0.1 .

		$a=0$	$a=0.25$	$a=0.5$	$a=0.75$	$a=1$
Reg.	L	0.007  (0.005)	0.021  (0.006)	0.123  (0.014)	0.570  (0.049)	1.038  (0.055)
	LC	0.151  (0.031)	0.148  (0.034)	0.198  (0.045)	0.416  (0.064)	0.637  (0.089)
	LL	0.083  (0.070)	0.095  (0.065)	0.120  (0.066)	0.257  (0.051)	0.411  (0.064)
	MY	0.079  (0.080)	0.135  (0.488)	14.058  (139.119)	0.622  (2.011)	3.253  (24.778)
Deriv.	L	1.869  (0.098)	1.404  (0.554)	1.098  (0.151)	1.053  (0.086)	1.014  (0.050)
	LL	0.770  (0.329)	10.960  (8.586)	1.397  (1.064)	0.533  (0.271)	0.418  (0.177)
	MY	2.244  (7.015)	19007.294  (188516.083)	24.681  (201.610)	23.785  (126.218)	1.217  (3.787)

Table 21. Table 21: Model (M3) with n ​ s ​ r = 0.2 𝑛 𝑠 𝑟 0.2 nsr=0.2 and ρ = 0.2 𝜌 0.2 \rho=0.2 .

		$a=0$	$a=0.25$	$a=0.5$	$a=0.75$	$a=1$
Reg.	L	0.014  (0.008)	0.030  (0.012)	0.131  (0.016)	0.573  (0.046)	1.045  (0.055)
	LC	0.182  (0.038)	0.197  (0.044)	0.252  (0.052)	0.470  (0.079)	0.754  (0.090)
	LL	0.105  (0.075)	0.111  (0.078)	0.180  (0.071)	0.323  (0.066)	0.496  (0.068)
	MY	0.449  (2.173)	0.218  (0.674)	0.221  (0.206)	0.617  (1.084)	5.052  (43.219)
Deriv.	L	1.877  (0.077)	1.551  (0.691)	1.077  (0.076)	1.035  (0.043)	1.017  (0.040)
	LL	0.935  (0.276)	8.854  (8.173)	1.804  (0.958)	0.639  (0.244)	0.507  (0.149)
	MY	116.568  (1067.638)	429.089  (3417.084)	14.088  (42.670)	10.561  (81.778)	62.822  (575.144)

Table 22. Table 22: Model (M3) with n ​ s ​ r = 0.4 𝑛 𝑠 𝑟 0.4 nsr=0.4 and ρ = 0.4 𝜌 0.4 \rho=0.4 .

		$a=0$	$a=0.25$	$a=0.5$	$a=0.75$	$a=1$
Reg.	L	0.042  (0.025)	0.059  (0.021)	0.161  (0.028)	0.601  (0.060)	1.058  (0.075)
	LC	0.303  (0.065)	0.303  (0.066)	0.359  (0.073)	0.632  (0.083)	0.939  (0.101)
	LL	0.191  (0.114)	0.184  (0.109)	0.254  (0.110)	0.451  (0.075)	0.705  (0.081)
	MY	0.235  (0.191)	0.215  (0.167)	0.296  (0.229)	0.544  (0.200)	1.255  (3.361)
Deriv.	L	1.869  (0.049)	2.023  (0.495)	1.146  (0.075)	1.044  (0.028)	1.012  (0.023)
	LL	1.331  (0.213)	10.177  (7.576)	1.902  (0.871)	0.823  (0.177)	0.721  (0.116)
	MY	11.547  (84.817)	309.921  (1678.421)	2772.892  (27642.691)	10.038  (37.057)	6.362  (47.083)

Equations

$$m(x+u)=m(x)+⟨m_x^{\prime },u⟩+\frac{1}{2}⟨m_{\zeta }^{\prime \prime }u,u⟩,$$

$$SWSE(a;\beta ):=\sum _{i=1}^n{\left(Y_i-a-⟨\beta ,X_i-x⟩\right)}^{2}K\left(h^{-1}\Vert X_i-x\Vert \right),$$

$$SWS{E}_J(a;b_1,\dots ,b_J):=\sum _{i=1}^n{\left(Y_i-a-\sum _{j=1}^J{b}_j⟨{\varphi }_j,X_i-x⟩\right)}^{2}K\left(h^{-1}\Vert X_i-x\Vert \right),$$

$$\boldsymbol{\Phi}\coloneqq\left[\begin{array}[]{cccc}1&\langle\phi_{1},\,X_{1}-x\rangle&\cdots&\langle\phi_{J},\,X_{1}-x\rangle\\ \vdots&\vdots&\ddots&\vdots\\ 1&\langle\phi_{1},\,X_{n}-x\rangle&\cdots&\langle\phi_{J},\,X_{n}-x\rangle\end{array}\right].$$

$$m(X_i)=m(x)+⟨m_x^{\prime },X_i-x⟩+\frac{1}{2}⟨m_{\zeta }^{\prime \prime }(X_i-x),X_i-x⟩\text{a.s.},$$

$$\left[\begin{array}[]{c}m(X_{1})\\ \vdots\\ m(X_{n})\end{array}\right]\,=\,\boldsymbol{\Phi}\,\left[\begin{array}[]{c}m(x)\\ \nabla m_{x}\end{array}\right]+\left[\begin{array}[]{c}\langle\mathcal{P}_{\mathcal{S}_{J}^{\perp}}m_{x}^{\prime},X_{1}-x\rangle\\ \vdots\\ \langle\mathcal{P}_{\mathcal{S}_{J}^{\perp}}m_{x}^{\prime},X_{n}-x\rangle\end{array}\right]+\frac{1}{2}\,\boldsymbol{R}_{\zeta,x}\quad\textrm{a.s.}$$

$$E_{\mathbf{X}}\left\{\hat{m}(x)\right\}$$

$$=m(x)+T_1+\frac{1}{2}{T}_2$$

$$\Vert {\hat{m}}_x^{\prime }-m_x^{\prime }\Vert$$

$$\Vert {\hat{m}}_x^{\prime }-E_{\mathbf{X}}\left({\hat{m}}_x^{\prime }\right)\Vert =O_P\left(h^{-1}\{{\lambda }_J^{2}n{\pi }_x(h){\}}^{-1/2}\right)+O_P\left(h^{-1}\{{\lambda }_{J}n{\pi }_x(h){\}}^{-1/2}\sqrt{J}\right).$$

$$E_{\mathbf{X}}\left({\hat{m}}_x^{\prime }\right)-m_x^{\prime }=-{\mathcal{P}}_{{\mathcal{S}}_J^{\perp }}{m}_x^{\prime }+Q_1+\frac{1}{2}{Q}_2,$$

$$E_{\mathbf{X}}\left(\Vert {\hat{m}}_x^{\prime }-E_{\mathbf{X}}\left({\hat{m}}_x^{\prime }\right){\Vert }^2\right)={\int }_0^1{\mathrm{Var}}_{\mathbf{X}}\left\{{\hat{m}}_x^{\prime }(t)\right\}\mathrm{d}t={\sigma }^2(x)Q\{1+o(1)\},$$

$$E_{\mathbf{X}}\left\{\hat{m}(x)\right\}=m(x)+O_P(h{J}^{-p-\nu })+O_P(h^2).$$

$$n^{-1}\sum _{i=1}^n{\Vert {\widehat{m_{X_i}^{\prime }}}^{boot}-\widehat{m_{X_i,-i}^{\prime }}\Vert }^2,$$

$$ORMSE{P}_{reg}:=\frac{1/500{\sum }_{i=n+1}^{n+500}{\left\{m(X_i)-\hat{m}(X_i)\right\}}^2}{1/500{\sum }_{i=n+1}^{n+500}{\left\{m(X_i)-\overline{m(X)}\right\}}^2}$$

$$ORMSE{P}_{deriv}:=\frac{1/500{\sum }_{i=n+1}^{n+500}{\Vert m_{X_i}^{\prime }-\widehat{m_{X_i}^{\prime }}\Vert }^2}{1/500{\sum }_{i=n+1}^{n+500}{\Vert m_{X_i}^{\prime }-\overline{m_X^{\prime }}\Vert }^2}$$

$$\hat{m}(X_j):=\sum _{i=1}^n{Y}_{i}K\left(\Vert X_j-X_i\Vert /h\right)/\sum _{i=1}^{n}K\left(\Vert X_j-X_i\Vert /h\right)$$

$$\text{mbox}heightat18=g\left(⟨\text{mbox}growthvelocityupto10,\beta ⟩\right)+error,$$

$$A_0$$

$$A_j$$

$$\boldsymbol{\Phi}^{\mathsf{T}}\widetilde{\boldsymbol{K}}\boldsymbol{\Phi}\,=\,\left[\begin{array}[]{c|c}\delta_{0}&\boldsymbol{\delta}^{\mathsf{T}}\\ \hline\cr&\\[-10.00002pt] \boldsymbol{\delta}&\boldsymbol{\Delta}\end{array}\right],$$

$$\left(\boldsymbol{\Phi}^{\mathsf{T}}\widetilde{\boldsymbol{K}}\boldsymbol{\Phi}\right)^{-1}=\left[\begin{array}[]{c|c}\mu&-\mu\,\boldsymbol{\delta}^{\mathsf{T}}\boldsymbol{\Delta}^{-1}\\ \hline\cr&\\[-10.00002pt] -\mu\,\boldsymbol{\Delta}^{-1}\boldsymbol{\delta}&\boldsymbol{\Delta}^{-1}\,+\,\mu\,\boldsymbol{\Delta}^{-1}\boldsymbol{\delta}\,\boldsymbol{\delta}^{\mathsf{T}}\boldsymbol{\Delta}^{-1}\end{array}\right]$$

$${\mathbf{\Delta }}^{-1}=b_{x,0,1}{b}_{x,2,1}^{-1}{h}^{-2}{\mathbf{\Gamma }}^{-1}\left\{1+o_P(1)\right\}.$$

$$T_1$$

$$-\left({\mathbf{\gamma }}^{\mathsf{T}}+{\mathbf{\theta }}^{\mathsf{T}}{O}_P\left(\{n{\pi }_x(h){\}}^{-1/2}\right)\right){\mathbf{\Gamma }}^{-1}{\mathbf{\alpha }}_{\mathbf{x},\mathbf{n}}^{\mathbf{bias}}h\left\{1+o_P(1)\right\}$$

$$\displaystyle\phantom{=}+\,\left(\boldsymbol{\gamma}^{\mathsf{T}}+\boldsymbol{\theta}^{\mathsf{T}}O_{P}\left(\{n\pi_{x}(h)\}^{-1/2}\right)\right)\boldsymbol{\Gamma}^{-1}\sqrt{\boldsymbol{\alpha_{x,n}^{var}}}\,O_{P}\left(h\,\{n\pi_{x}(h)\}^{-1/2}\right)\Big{\}}$$

$$\left({\mathbf{\gamma }}^{\mathsf{T}}+{\mathbf{\theta }}^{\mathsf{T}}{O}_P\left(\{n{\pi }_x(h){\}}^{-1/2}\right)\right){\mathbf{\Gamma }}^{-1}{\mathbf{\alpha }}_{\mathbf{x},\mathbf{n}}^{\mathbf{bias}}$$

$$\left({\mathbf{\gamma }}^{\mathsf{T}}+{\mathbf{\theta }}^{\mathsf{T}}{O}_P\left(\{n{\pi }_x(h){\}}^{-1/2}\right)\right){\mathbf{\Gamma }}^{-1}\sqrt{{\mathbf{\alpha }}_{\mathbf{x},\mathbf{n}}^{\mathbf{var}}}$$

$${\mathbf{\delta }}^{\mathsf{T}}{\mathbf{\Delta }}^{-1}\mathbf{\delta }$$

$$+O_P\left(\{n{\pi }_x(h){\}}^{-1/2}\right){\mathbf{\theta }}^{\mathsf{T}}{\mathbf{\Gamma }}^{-1}\mathbf{\gamma }+O_P\left(\{n{\pi }_x(h){\}}^{-1}\right){\mathbf{\theta }}^{\mathsf{T}}{\mathbf{\Gamma }}^{-1}\mathbf{\theta }.$$

$${\mathbf{\delta }}^{\mathsf{T}}{\mathbf{\Delta }}^{-1}\mathbf{\delta }=O_P\left({\lambda }_J^{-1}\right),\text{mbox}whichresultsin\mu =O_P\left({\lambda }_J\right).$$

$$T_{21}$$

$$-b_{x,1,1}{b}_{x,2,1}^{-1}{b}_{x,3,1}\left({\mathbf{\gamma }}^{\mathsf{T}}+{\mathbf{\theta }}^{\mathsf{T}}{O}_P\left(\{n{\pi }_x(h){\}}^{-1/2}\right)\right){\mathbf{\Gamma }}^{-1}{\mathbf{\beta }}_{\mathbf{x}}^{\mathbf{bias}}{h}^2\left\{1+o_P(1)\right\}$$

$$+\left({\mathbf{\gamma }}^{\mathsf{T}}+{\mathbf{\theta }}^{\mathsf{T}}{O}_P\left(\{n{\pi }_x(h){\}}^{-1/2}\right)\right){\mathbf{\Gamma }}^{-1}\sqrt{{\mathbf{\beta }}_{\mathbf{x}}^{\mathbf{var}}}{O}_P\left(h^2\{n{\pi }_x(h){\}}^{-1/2}\right)\}.$$

Full text

keywords: Asymptotics, functional data, functional derivative of regression operator, functional index model, local linear regression, scalar-on-function regression

1 Introduction

Functional data analysis is a toolbox of statistical techniques for dealing with datasets of random functions [46, 47, 31, 37]. When regressing a scalar response $Y$ on an explanatory random function $X$ using the model $Y=m(X)+error$, with $m$ unknown, common terminology refers to the scalar-on-function regression. In this setting, linear modeling [10, 9, 13] or its generalized versions [33, 41] have been intensively studied in the literature. Although the scalar-on-function linear regression is a powerful tool, its lack of flexibility when nonlinearities occur led the statistical community to develop a nonparametric approach as it was in the multivariate case. For an overview on the functional version of the Nadaraya-Watson kernel estimator see [22]. Since the development of local linear estimation for multivariate data in the 1990’s (see for instance [16, 18, 17, 51, 12, 26] and the monograph [19]), it is well known that local linear regression outperforms the usual Nadaraya-Watson kernel estimator (i.e. the local constant regressor). Therefore, it became perhaps the most popular nonparametric regression technique. Surprisingly, in the framework of functional data, there are only two papers that focus on the estimation of the regression operator in the scalar-on-function local linear regression model. In [5], a projection approach to the problem similar to the study presented here is proposed, but the asymptotics derived in that paper appears to suffer from a lack of rigorousness. Paper [6] is a pure theoretical work providing an alternative estimating procedure by regularizing a non-bounded linear operator. The latter method does not directly relate to the approach taken here, as it does not rely on the use of projections. All in all, the theory and practice of scalar-on-function local linear regression is severely underdeveloped, and thus the method is far from being as popular as it is in the multivariate case.

This paper introduces the local linear regression as an indispensable tool in the setting of scalar-on-function nonparametric regression. It turns out that the functional local linear regression (that is, local linear regression when the regressor is a random function) is not only a smart method of estimating the regression operator. As an exciting by-product we obtain an easy and fast method for estimating the functional derivative $m^{\prime}_{x}$ of the regression operator $m$ at any function $x$. The functional derivative is a linear functional that represents a local linear approximation to the regression operator $m$ around $x$ (for a precise statement see (H(H1)) below). In what follows, we use the Riesz representation theorem, and identify the functional derivative $m^{\prime}_{x}$ with its unique representing function. What makes the estimation of functional derivatives of such great interest? A first motivation is given in the pioneering works [28, 43] where estimating procedures are developed without considering the local linear regression setting. There, it was convincingly demonstrated that the concept of functional derivative greatly facilitates the interpretation of results. As a further step in the pursuit for understanding how one can use the functional derivative in a natural way, let us consider the functional Taylor expansion of the regression operator. For a small positive real $\eta$ and a direction $u$ (i.e. a function $u$ such that $\|u\|=1$), Taylor’s expansion and the Riesz representation theorem allow us to write $m(x+\eta\,u)-m(x)=\eta\langle m_{x}^{\prime},u\rangle+O(\eta^{2})$. A first order approximation of the magnitude of the difference $m(x+\eta\,u)-m(x)$ is therefore the interval $[-\eta\|m^{\prime}_{x}\|,\eta\|m^{\prime}_{x}\|]$ — the smaller $\|m^{\prime}_{x}\|$ is, the less sensitive to small perturbations in $x$ is $m$. In a sense, $\|m^{\prime}_{x}\|$ can be seen as a measure of reliability for the prediction of $m$ at $x$. For another example where the functional derivative appears as a successful tool in interesting statistical problems, consider, for instance, the single-functional index model [4, 1, 11, 34]. That model takes the form $m(x)=\mu+g\left(\langle\beta,x\rangle\right)$, where $\mu$ is an unknown scalar and the scalar response interacts with the functional covariate only through an unknown functional direction $\beta$ combined with an unknown real-valued link function $g$. Extending the average derivative estimation method introduced in [29] to the functional setting, it is easy to show that $\operatorname{E}\left(m^{\prime}_{X}\right)$ is proportional to the functional direction $\beta$. Thus, given a sample $(X_{1},Y_{1}),\ldots,(X_{n},Y_{n})$, as soon as one is able to obtain estimates $\widehat{m_{X_{1}}^{\prime}},\ldots,\widehat{m_{X_{n}}^{\prime}}$ of the functional derivatives $m_{X_{1}}^{\prime},\ldots,m_{X_{n}}^{\prime}$, one can compute $\widehat{\operatorname{E}\left(m^{\prime}_{X}\right)}\coloneqq n^{-1}\sum_{i}\widehat{m_{X_{i}}^{\prime}}$. The quantity $\widehat{\operatorname{E}\left(m^{\prime}_{X}\right)}/\|\widehat{\operatorname{E}\left(m^{\prime}_{X}\right)}\|$ is a reasonable estimator of the functional index $\beta$.

All the examples above emphasize the major role that the functional derivative of the regression operator plays in important aspects of statistics: interpretation, reliability and methodology. This is why we propose to revisit the functional local linear regression by focusing not only on the regression operator, but also on its functional derivative. In this work, a projection approach to functional local linear estimation is adopted. Theoretical properties of the estimator of the regression operator $m$ are stated. Simultaneously, a local linear estimator of the functional derivative $m^{\prime}_{x}$ at $x$ is proposed. For both estimators of the regression operator and its functional derivative, original technical tools are developed in order to study their theoretical behavior. If the implementation of the estimator of the regression operator is straightforward, selection of the smoothing parameter for estimating its functional derivative poses a major challenge. This is why an ad hoc bootstrap procedure is introduced to pilot the bandwidth choice. Ease of implementation as well as nice finite sample properties of our local linear estimators are highlighted in a simulation study. The paper is concluded with a benchmark real dataset that is used to illustrate the important role of the functional derivatives in functional data analysis. Most of the proofs given in the paper are deferred to Appendices A and B. The paper is complemented by extensive supplementary material that includes an efficient R [44] implementation of the newly proposed methods333Available at https://bitbucket.org/StanislavNagy/fllr and some additional examples in Appendix C, and theoretical results in Appendix D.

2 Functional local linear estimation

Let $X$ be an $H$-valued random function where $H$ is the separable Hilbert space of square integrable functions defined on $[0,1]$ equipped with the inner product $\langle\cdot,\cdot\rangle$, and let $\|\cdot\|$ be the associated norm. This work focuses on the relationship between $X$ and a scalar response $Y$ by considering the nonparametric regression model $Y\,=\,m(X)\,+\,\varepsilon$ where $\operatorname{E}(\varepsilon|X)=0$. The regression operator $m$ mapping $H$ into $\mathbb{R}$ is unknown, and is assumed to be smooth enough in a neighborhood $\mathcal{N}_{x}$ of a given $x\in H$.

• (H1)

For any $u\in\mathcal{N}_{x}$, there exists $\zeta=x+t\,u$ with $t\in(0,1)$ such that

[TABLE]

where $m^{\prime}_{x}\in H$, $m^{\prime\prime}_{\zeta}$ is a Hilbert-Schmidt linear operator mapping $H$ into $H$, and $v\mapsto m_{v}^{\prime\prime}$ is Lipschitz in $v\in\mathcal{N}_{x}$.

In other words, one focuses on those regression operators $m$ for which the second order Taylor expansion is valid. Note that condition (H(H1)) is the functional counterpart of what is standardly required in the finite-dimensional local linear regression setting.

Based on an $n$-sample $(X_{i},\,Y_{i})_{i=1,\ldots,n}$ of independent identically distributed (iid) copies of $(X,\,Y)$, our main task is to estimate the regression operator $m$, as well as the functional derivative $m^{\prime}_{x}$. To this end, one extends the sum of weighted squared errors (SWSE) to the functional setting

[TABLE]

where $K(\cdot)$ is a kernel function defined on $[0,1]$ and $h$ is a positive smoothing parameter (a bandwidth). The principle of the local linear criterion is to linearize the regression operator in a neighborhood of $x$. In other words, for any $X_{i}$ close to $x$, one considers that $\operatorname{E}(Y_{i}|X_{i})=a+\langle\beta,X_{i}-x\rangle$. Then, the real $a$ (resp. the square integrable function $\beta$) can be interpreted as the regression operator (resp. the functional derivative) at $x$. Consequently, the estimation of the functional derivative $m^{\prime}_{x}$ will be based on the estimation of the function $\beta$. Instead of focusing on $\beta$ itself, it is computationally advantageous to consider its projection $\sum_{j\leq J}\langle\phi_{j},\beta\rangle\phi_{j}$ onto the $J$-dimensional subspace $\mathcal{S}_{J}$ of $H$ spanned by the orthonormal sequence $\phi_{1},\ldots,\phi_{J}$ that is completed by $\phi_{J+1},\phi_{J+2},\ldots$ in order to get an orthonormal basis of $H$. In this notation the SWSE criterion is approximated by

[TABLE]

where for any $j$, $b_{j}\coloneqq\langle\phi_{j},\beta\rangle$. Then, $\widehat{m}(x)\coloneqq\widehat{a}$ and $\widehat{m^{\prime}_{x}}\coloneqq\sum_{j=1}^{J}\widehat{b}_{j}\phi_{j}$ where $(\widehat{a};\,\widehat{b}_{1},\ldots,\widehat{b}_{J})\,\coloneqq\,\arg\inf_{(a;\,b_{1},\ldots,b_{J})}SWSE_{J}(a;\,b_{1},\ldots,b_{J}).$ By using vector and matrix notations, one is able to express the local linear estimators. Let $\boldsymbol{Y}\coloneqq[Y_{1},\dots,Y_{n}]^{\mathsf{T}}$, $\boldsymbol{K}\coloneqq\mbox{diag}\left\{K\left(h^{-1}\|X_{1}-x\|\right),\dots,K\left(h^{-1}\|X_{n}-x\|\right)\right\}$ the diagonal $n\times n$ matrix and $\boldsymbol{\Phi}$ the following $n\times(J+1)$ matrix

[TABLE]

Define $\widehat{\boldsymbol{b}}\coloneqq\left[\widehat{b}_{1},\dots,\widehat{b}_{J}\right]^{\mathsf{T}}$, $\boldsymbol{\phi}\coloneqq\left[\phi_{1},\dots,\phi_{J}\right]^{\mathsf{T}}$, $\boldsymbol{0}$ the $J\times 1$ null vector, $\boldsymbol{e}$ the $(J+1)$-dimensional vector $\left[1,\boldsymbol{0}^{\mathsf{T}}\right]^{\mathsf{T}}$ and $\boldsymbol{I}$ the $J\times J$ identity matrix. Then, it is easy to see that $\left[\widehat{a}\,|\,\widehat{\boldsymbol{b}}^{\mathsf{T}}\right]^{\mathsf{T}}\,=\,\left(\boldsymbol{\Phi}^{\mathsf{T}}\boldsymbol{K}\boldsymbol{\Phi}\right)^{-1}\boldsymbol{\Phi}^{\mathsf{T}}\boldsymbol{K}\boldsymbol{Y}$, $\widehat{m}(x)\,=\,\boldsymbol{e}^{\mathsf{T}}\,\left(\boldsymbol{\Phi}^{\mathsf{T}}\boldsymbol{K}\boldsymbol{\Phi}\right)^{-1}\boldsymbol{\Phi}^{\mathsf{T}}\boldsymbol{K}\boldsymbol{Y}$ and $\widehat{m^{\prime}_{x}}\,=\,\boldsymbol{\phi}^{\mathsf{T}}\,\left[\boldsymbol{0}|\boldsymbol{I}\right]\left(\boldsymbol{\Phi}^{\mathsf{T}}\boldsymbol{K}\boldsymbol{\Phi}\right)^{-1}\boldsymbol{\Phi}^{\mathsf{T}}\boldsymbol{K}\boldsymbol{Y}$. The local linear approach has the nice property that both $\widehat{m}(x)$ and $\widehat{m^{\prime}_{x}}$ are based on the common terms $\left(\boldsymbol{\Phi}^{\mathsf{T}}\boldsymbol{K}\boldsymbol{\Phi}\right)^{-1}\boldsymbol{\Phi}^{\mathsf{T}}\boldsymbol{K}\boldsymbol{Y}$ and this is why the raw computational cost of $\widehat{m^{\prime}_{x}}$ is not much higher than the one for $\widehat{m}(x)$.

3 Asymptotic study

Let us first focus on the assumptions we need to derive the asymptotic behavior of $\widehat{m}(x)$ and $\widehat{m^{\prime}_{x}}$.

• (H3)

The kernel function $K$ is continuously differentiable on its support $(0,1)$ with $K^{\prime}(s)\leq 0$ for all $s\in(0,1)$ and $K(1)>0$.

• (H5)

For any integers $j_{1},\ldots,j_{M},p_{1},\ldots,p_{M}\geq 0$ with $M\geq 1$, let us define $\gamma_{j_{1},\ldots,j_{M}}^{p_{1},\ldots,p_{M}}(t)\coloneqq\operatorname{E}\left(\langle\phi_{j_{1}},X_{1}-x\rangle^{p_{1}}\cdots\langle\phi_{j_{M}},X_{1}-x\rangle^{p_{M}}|\|X_{1}-x\|^{p_{1}+\cdots+p_{M}}=t\right)$ and let ${\gamma_{j_{1},\ldots,j_{M}}^{p_{1},\ldots,p_{M}}}^{{}^{\prime}}(t)$ be its derivative at $t$. The functions $\gamma_{j_{1}}^{1}$, $\gamma_{j_{1},j_{2}}^{1,1}$, $\ldots$, $\gamma_{j_{1},\ldots,j_{4}}^{1,\ldots,1}$, $\gamma_{j_{1}}^{2}$, $\gamma_{j_{1},j_{2},j_{3}}^{2,1,1}$, $\gamma_{j_{1},\dots,j_{5}}^{2,1,1,1,1}$ and $\gamma_{j_{1},j_{2}}^{2,2}$ are assumed to be continuously differentiable around zero and the smallest eigenvalue $\lambda_{J}$ of the $J\times J$ matrix $\boldsymbol{\Gamma}$, whose $(j,k)$-th element is defined by $\left[\boldsymbol{\Gamma}\right]_{jk}\coloneqq{\gamma_{j,k}^{1,1}}^{{}^{\prime}}(0)$, is strictly positive.

• (H7)

$h=h_{n}$ tends to 0 with $n$, $J=J_{n}$ and $n\,\pi_{x}(h)$ grow to infinity with $n$ so that $h\,J^{1/2}=o(1)$, $h^{-1}\lambda_{J}^{-1}\{n\,\pi_{x}(h)\}^{-1/2}=o(1)$ and $h^{-1}(\lambda_{J}/J)^{-1/2}\{n\,\pi_{x}(h)\}^{-1/2}=o(1)$, where $\pi_{x}(h)\coloneqq P\left(\|X_{1}-x\|<h\right)$.

• (H9)

For all $s$ in $[0,1]$, the ratio $\displaystyle\tau_{x,h}(s)\coloneqq\frac{\pi_{x}(hs)}{\pi_{x}(h)}$ tends to $\tau_{x}(s)$ as $h$ goes to 0.

• (H11)

The conditional variance of the error $\sigma^{2}(x)=\operatorname{Var}(Y|X=x)$ mapping $H$ into $\mathbb{R}$ is a uniformly continuous operator.

The original hypothesis (H(H5)) introduces a particular family of functions. It is shown in Lemma 6 (postponed to Appendix B) that the functions $\gamma_{j_{1},\ldots,j_{M}}^{p_{1},\ldots,p_{M}}$ have interesting properties, in particular the matrix $\boldsymbol{\Gamma}$ is positive semi-definite. Of course, since the size $J$ of the square matrix $\boldsymbol{\Gamma}$ tends to infinity with $n$, it is clear that its smallest eigenvalue $\lambda_{J}$ tends to zero with $n$. Nevertheless, it is not too much restrictive to require $\lambda_{J}$ strictly positive for any fixed $J$. It is worth noting that the differentiability assumptions imposed on $\gamma_{j_{1},\ldots,j_{M}}^{p_{1},\ldots,p_{M}}$ make this function very useful to approximate $\operatorname{E}\left\{\langle\phi_{j_{1}},X_{1}-x\rangle^{p_{1}}\cdots\langle\phi_{j_{M}},X_{1}-x\rangle^{p_{M}}\,K^{q}\left(h^{-1}\|X_{1}-x\|\right)\right\}$ that plays a major role in the asymptotic behavior of the estimators $\widehat{m}_{x}$ and $\widehat{m}^{\prime}_{x}$. The reader will find in Appendix D (see Lemma 12) an interesting result providing a general situation where (H(H5)) is fulfilled. Condition (H(H9)) is a more classical assumption. For standard families of processes an explicit form of the function $\tau_{x}(s)$ is available (for more details, see [20] and references therein).

Let $\boldsymbol{X}$ stand for the sample $X_{1},\ldots,X_{n}$ and let $\operatorname{E}_{\boldsymbol{X}}$ (resp. $\operatorname{Var}_{\boldsymbol{X}}$) be the conditional expectation (resp. variance) with respect to $\boldsymbol{X}$. The asymptotic conditional bias and variance of $\widehat{m}(x)$ are provided in the following theorem.

Theorem 1**.**

Under conditions (H(H1))–(H(H11)),

1. (i)

$\operatorname{E}_{\boldsymbol{X}}\left\{\widehat{m}(x)\right\}\,=\,m(x)\,+\,O_{P}\left(\|\mathcal{P}_{\mathcal{S}_{J}^{\perp}}m^{\prime}_{x}\|\,h\right)\,+\,O_{P}\left(h^{2}\right),$** 2. (ii)

$\operatorname{Var}_{\boldsymbol{X}}\left\{\widehat{m}(x)\right\}\,=\,O_{P}\left(\{n\,\pi_{x}(h)\}^{-1}\right)$,

where $\mathcal{S}_{J}^{\perp}$ is the orthogonal complement to the space $\mathcal{S}_{J}$ in $H$, and $\mathcal{P}_{\mathcal{S}_{J}^{\perp}}\colon H\to\mathcal{S}_{J}^{\perp}$ is the orthogonal projection onto $\mathcal{S}_{J}^{\perp}$.

In the same situation (see for instance [22]), the conditional bias of the estimated functional local constant regression is of order $h$ with the same conditional variance. However, the quantity $\|\mathcal{P}_{\mathcal{S}_{J}^{\perp}}m^{\prime}_{x}\|$ tends to zero when $n$ grows to infinity since $m^{\prime}_{x}$ is a square integrable function. Therefore $\widehat{m}(x)$ outperforms the asymptotic behavior of the kernel estimator of the functional local constant regression. Note that the conditional variance in Theorem 1 involves the small ball probability $\pi_{x}(h)$, which is standard in the functional nonparametric setting [22].

The functional setting involves implicitly the dimension $J$ of the approximating subspace $\mathcal{S}_{J}$ in the rate of convergence. Indeed, the asymptotic behavior of the conditional bias depends on the quantity $\|\mathcal{P}_{\mathcal{S}_{J}^{\perp}}m^{\prime}_{x}\|$, which assesses the approximation error of the functional derivative $m^{\prime}_{x}$ in $\mathcal{S}_{J}$. From a theoretical point of view, $\widehat{m}(x)$ involves $J\times J$ matrices and $J$-dimensional vectors with $J$ converging to infinity. This makes the asymptotic study much harder in comparison with the multivariate (finite-dimensional) setting. The issue of infinite dimension is overcome by deriving accurately, in an element-wise sense, the asymptotic behavior of the matrices and vectors involved in both local linear estimators.

Proof of the conditional bias of $\widehat{m}(x)$.

Here only the main guidelines of the proof are given. Details and technical lemmas are postponed to the appendix.

For any $X_{i}$, the Taylor expansion of $m$ at $x$ can be expressed as

[TABLE]

where $\zeta=x+t(X_{i}-x)$ with $t\in(0,1)$. By using a basis expansion of $m_{x}^{\prime}$ in terms of $\phi_{1},\phi_{2},\ldots$, we get $m(X_{i})=m(x)+\sum_{j\leq J}\langle m_{x}^{\prime},\phi_{j}\rangle\langle\phi_{j},X_{i}-x\rangle+\langle\mathcal{P}_{\mathcal{S}_{J}^{\perp}}m_{x}^{\prime},X_{i}-x\rangle+R_{\zeta,x,i}/2$ almost surely with $R_{\zeta,x,i}\coloneqq\langle m_{\zeta}^{\prime\prime}(X_{i}-x),X_{i}-x\rangle$. Let $\nabla m_{x}\coloneqq\left[\langle m_{x}^{\prime},\phi_{1}\rangle,\ldots,\langle m_{x}^{\prime},\phi_{J}\rangle\right]^{\mathsf{T}}$ be the first $J$ coordinates of the gradient of $m$ at $x$ and $\boldsymbol{R}_{\zeta,x}\coloneqq[R_{\zeta,x,1},\ldots,R_{\zeta,x,n}]^{\mathsf{T}}$. Then we can write

[TABLE]

According to the definition of $\widehat{m}(x)$,

[TABLE]

where $\displaystyle T_{1}\coloneqq\boldsymbol{e}^{\mathsf{T}}\,\left(\boldsymbol{\Phi}^{\mathsf{T}}\boldsymbol{K}\boldsymbol{\Phi}\right)^{-1}\boldsymbol{\Phi}^{\mathsf{T}}\boldsymbol{K}\left[\langle\mathcal{P}_{\mathcal{S}_{J}^{\perp}}m_{x}^{\prime},X_{1}-x\rangle,\dots,\langle\mathcal{P}_{\mathcal{S}_{J}^{\perp}}m_{x}^{\prime},X_{n}-x\rangle\right]^{\mathsf{T}}$ and $\displaystyle T_{2}\coloneqq\boldsymbol{e}^{\mathsf{T}}\,\left(\boldsymbol{\Phi}^{\mathsf{T}}\boldsymbol{K}\boldsymbol{\Phi}\right)^{-1}\boldsymbol{\Phi}^{\mathsf{T}}\boldsymbol{K}\boldsymbol{R}_{\zeta,x}^{\mathsf{T}}$. The next lemma focuses on the asymptotic behavior of $T_{1}$ and $T_{2}$. Its proof can be found in Appendix A.

Lemma 1**.**

As soon as conditions (H(H1))–(H(H9)) are fulfilled,

1. (i)

$T_{1}\,=\,O_{P}\left(\|\mathcal{P}_{\mathcal{S}_{J}^{\perp}}m_{x}^{\prime}\|\,h\right)$, 2. (ii)

$T_{2}\,=\,O_{P}\left(h^{2}\right)$.

Now, it is enough to plug in Lemma 1 with (3) to get the claimed conditional bias. ∎

Proof of the conditional variance of $\widehat{m}(x)$.

By the assumptions, the covariance matrix of $\boldsymbol{Y}$ given $X_{1},\ldots,X_{n}$ is the diagonal matrix $\mbox{diag}\{\sigma^{2}(X_{1}),\ldots,\sigma^{2}(X_{n})\}$ and $\operatorname{Var}_{\boldsymbol{X}}\{\widehat{m}(x)\}=\boldsymbol{e}^{\mathsf{T}}\left(\boldsymbol{\Phi}^{\mathsf{T}}\boldsymbol{K}\boldsymbol{\Phi}\right)^{-1}\boldsymbol{\Phi}^{\mathsf{T}}\boldsymbol{D}\boldsymbol{\Phi}\left(\boldsymbol{\Phi}^{\mathsf{T}}\boldsymbol{K}\boldsymbol{\Phi}\right)^{-1}\boldsymbol{e}$. Here, $\boldsymbol{D}$ is also a diagonal matrix such that $[\boldsymbol{D}]_{ii}\,=\,\sigma^{2}(X_{i})\,[\boldsymbol{K}]_{ii}^{2}$. According to (H(H3)), (H(H7)) and (H(H11)), $[\boldsymbol{D}]_{ii}\,=\,\{\sigma^{2}(x)\,+\,o(1)\}\,[\boldsymbol{K}]_{ii}^{2}$ and $\operatorname{Var}_{\boldsymbol{X}}\{\widehat{m}(x)\}\,=\,\{\sigma^{2}(x)\,+\,o(1)\}\,\boldsymbol{e}^{\mathsf{T}}\,\left(\boldsymbol{\Phi}^{\mathsf{T}}\boldsymbol{K}\boldsymbol{\Phi}\right)^{-1}\,\boldsymbol{\Phi}^{\mathsf{T}}\,\boldsymbol{K}^{2}\,\boldsymbol{\Phi}$ $\left(\boldsymbol{\Phi}^{\mathsf{T}}\boldsymbol{K}\boldsymbol{\Phi}\right)^{-1}\boldsymbol{e}.$ The remainder of the proof is based on the technical Lemma 3, which is postponed to Appendix A. ∎

Once the theoretical properties of $\widehat{m}(x)$ have been given, a natural and interesting issue concerns the asymptotic behavior of the functional derivative $\widehat{m}_{x}^{\prime}$ of $m$ at $x$. The next result details the conditional bias and variance of $\widehat{m}_{x}^{\prime}$.

Theorem 2**.**

As soon as (H(H1))–(H(H11)) are fulfilled, conditionally to $X_{1},\ldots,X_{n}$,

[TABLE]

The first two summands in the formula above correspond to the conditional bias of $\widehat{m}_{x}^{\prime}$, while the remaining terms come from its conditional variance. As $n$ tends to infinity, $h$ tends to zero, the dimension $J$ goes to infinity, and the smallest eigenvalue $\lambda_{J}$ to zero. Therefore, the rate of convergence of $\widehat{m}_{x}^{\prime}$ is slower than that of $\widehat{m}(x)$. Compared to Theorem 1, one $h$ is removed in the conditional bias, and $h^{-1}$ is added to the terms that relate to conditional variance, which corresponds to the standard degradation of the convergence rate observed in the multivariate case. But the functional setting adds specific terms like $J$ and $\lambda_{J}$ which deteriorate the asymptotic behavior. Nevertheless, as pointed out below, the finite sample properties of $\widehat{m}_{x}^{\prime}$ are surprisingly good.

Proof of Theorem 2.

The asymptotic behavior of $\|\widehat{m}_{x}^{\prime}-m_{x}^{\prime}\|$ is a direct by-product of the next results which detail the conditional bias and variance of the functional derivative $\widehat{m}_{x}^{\prime}$.

Lemma 2**.**

As soon as (H(H1))–(H(H11)) are fulfilled,

1. (i)

$\|\operatorname{E}_{\boldsymbol{X}}\left(\widehat{m}_{x}^{\prime}\right)-m_{x}^{\prime}\|\,=\,O_{P}\left(\lambda_{J}^{-1}\|\mathcal{P}_{\mathcal{S}_{J}^{\perp}}m^{\prime}_{x}\|\right)\,+\,O_{P}\left(\lambda_{J}^{-1}h\right),$** 2. (ii)

and conditionally on $X_{1},\ldots,X_{n}$,

[TABLE]

Proof of Lemma 2-$(i)$.

Let us first focus on the conditional bias of $\widehat{m}_{x}^{\prime}$. By the definition of $\widehat{m}_{x}^{\prime}$, $\operatorname{E}_{\boldsymbol{X}}\left(\widehat{m}_{x}^{\prime}\right)\,=\,\boldsymbol{\phi}^{\mathsf{T}}\,[\boldsymbol{0}|\boldsymbol{I}]\,\left(\boldsymbol{\Phi}^{\mathsf{T}}\boldsymbol{K}\boldsymbol{\Phi}\right)^{-1}\boldsymbol{\Phi}^{\mathsf{T}}\boldsymbol{K}\left[m(X_{1}),\dots,m(X_{n})\right]^{\mathsf{T}}$ so that, by (1) and (2),

[TABLE]

for $\displaystyle Q_{1}=\boldsymbol{\phi}^{\mathsf{T}}\,[\boldsymbol{0}|\boldsymbol{I}]\,\left(\boldsymbol{\Phi}^{\mathsf{T}}\boldsymbol{K}\boldsymbol{\Phi}\right)^{-1}\boldsymbol{\Phi}^{\mathsf{T}}\boldsymbol{K}\left[\langle\mathcal{P}_{\mathcal{S}_{J}^{\perp}}m_{x}^{\prime},X_{1}-x\rangle,\dots,\langle\mathcal{P}_{\mathcal{S}_{J}^{\perp}}m_{x}^{\prime},X_{n}-x\rangle\right]^{\mathsf{T}}$, and $\displaystyle Q_{2}=\boldsymbol{\phi}^{\mathsf{T}}\,[\boldsymbol{0}|\boldsymbol{I}]\,\left(\boldsymbol{\Phi}^{\mathsf{T}}\boldsymbol{K}\boldsymbol{\Phi}\right)^{-1}\boldsymbol{\Phi}^{\mathsf{T}}\boldsymbol{K}\boldsymbol{R}_{\zeta,x}^{\mathsf{T}}$. Lemma 4 given in Appendix A together with (4) provide the claimed expression for the conditional bias. ∎

Proof of Lemma 2-$(ii)$.

Following the proof of the conditional variance of $\widehat{m}(x)$

[TABLE]

where $\displaystyle Q=\int_{0}^{1}\,\boldsymbol{\phi}(t)^{\mathsf{T}}\,[\boldsymbol{0}|\boldsymbol{I}]\left(\boldsymbol{\Phi}^{\mathsf{T}}\boldsymbol{K}\boldsymbol{\Phi}\right)^{-1}\,\boldsymbol{\Phi}^{\mathsf{T}}\,\boldsymbol{K}^{2}\,\boldsymbol{\Phi}\,\left(\boldsymbol{\Phi}^{\mathsf{T}}\boldsymbol{K}\boldsymbol{\Phi}\right)^{-1}[\boldsymbol{0}|\boldsymbol{I}]^{\mathsf{T}}\boldsymbol{\phi}(t)\,\mathrm{d}\,t$. The remainder of the proof consists in decomposing this conditional variance according to formula (9) given in Appendix A with by-products of Lemma 3. All details are postponed to the appendix. ∎

The statement of Theorem 2 now follows immediately. ∎

4 Examples of approximating bases

The functional local linear estimator depends on the basis $\phi_{1},\phi_{2},\ldots$ In this section we specify the asymptotic behavior of our estimator when considering particular bases. We start with usual approximating function spaces (cases 1 and 2 below). A more challenging issue consists in replacing the deterministic basis $\phi_{1},\phi_{2},\ldots$ with a data driven one. This important question is investigated in case 3.

Case 1: Orthogonal B-spline basis. B-spline basis is a well-known and useful tool for approximating smooth functions. Consider the set of functions that are polynomials of degree $q$ on each interval $[(t-1)/k,\,t/k]$ for $t=1,\ldots,k$ and are $(q-1)$ times continuously differentiable on $[0,\,1]$. This $(k+q)$-dimensional subspace of $H$ defines the well-known space of splines. One can derive an orthogonal basis of B-splines of that space $\{B_{k,1},\ldots,B_{k,k+q}\}$ (see [15] for an overview on spline functions and [49] for the orthogonalization of B-spline basis functions). In this situation, one sets $J=k+q$ and for any $j\leq J$, $\phi_{j}=B_{k,j}/\|B_{k,j}\|$. In this setting $k=k_{n}$ is a sequence that grows to infinity with $n$. Now, let the functional $m_{x}^{\prime}$ be smooth enough so that its $p$th derivative is a Hölder function:

• (H13)

$|m_{x}^{\prime(p)}(u)-m_{x}^{\prime(p)}(v)|\leq C|u-v|^{\nu}$ with $\nu\in[0,1]$.

From Theorem XII.1 in [15] and (H(H13)), $\|\mathcal{P}_{\mathcal{S}_{J}^{\perp}}m_{x}^{\prime}\|\leq C\,(J-q)^{-(p+\nu)}$. Because $\|P_{\mathcal{S}_{J}}m_{x}^{\prime}-m_{x}^{\prime}\|=\|P_{\mathcal{S}_{J}^{\perp}}m_{x}^{\prime}\|$, the rate of convergence of the conditional bias becomes

[TABLE]

Case 2: Fourier basis. When one suspects some periodic features for the functional $m_{x}^{\prime}$, it can be advantageous to expand the functions by means of the Fourier basis $\phi_{1}(t)=1$, $\phi_{2j}(t)=\sqrt{2}\sin(2\pi jt)$, and $\phi_{2j+1}(t)=\sqrt{2}\cos(2\pi jt)$ for $j=1,2,\ldots$ If one assumes that

• (H15)

$m_{x}^{\prime}$ is a periodic function,

then, according to [57] and assumptions (H(H13)) and (H(H15)), one has $\|\mathcal{P}_{\mathcal{S}_{J}^{\perp}}m_{x}^{\prime}\|\,=\,O\left(J^{-p-\nu}\right)$ which leads to the same rate of convergence (5).

Case 3: Data driven basis. The functional principal components analysis (FPCA) allows to expand a random function $X$ into a basis $\phi_{1},\phi_{2},\ldots$ in an optimal way (see [35, 38, 48, 14] for precursor works and [7, 56, 25, 27] for more recent statistical developments). In this setting, functions $\phi_{j}$ are the eigenfunctions of the covariance operator of $X$ and the eigenanalysis of the empirical covariance operator provides a data driven basis $\widehat{\phi}_{1},\widehat{\phi}_{2},\ldots$ (by convention, we assume that $\langle\phi_{j},\widehat{\phi}_{j}\rangle>0$). We propose to investigate the asymptotic properties of the functional local linear estimator when replacing $\phi_{1},\ldots,\phi_{J}$ with the data driven basis $\widehat{\phi}_{1},\ldots,\widehat{\phi}_{J}$. This results in the new estimator $\widehat{\widehat{m}}(x)\,=\,\boldsymbol{e}^{\mathsf{T}}\,\left(\widehat{\boldsymbol{\Phi}}^{\mathsf{T}}\boldsymbol{K}\widehat{\boldsymbol{\Phi}}\right)^{-1}\widehat{\boldsymbol{\Phi}}^{\mathsf{T}}\boldsymbol{K}\boldsymbol{Y}$ where, for $i=1,\ldots,n$, $[\widehat{\boldsymbol{\Phi}}]_{i1}=1$ and for $j=2,\ldots,J$, $[\widehat{\boldsymbol{\Phi}}]_{ij}=\langle\widehat{\phi}_{j},\,X_{i}-x\rangle$. Thanks to [7] and [10], as soon as $\operatorname{E}\|X\|^{4}<\infty$, one has for any $j=1,2,\ldots$ that $\|\widehat{\phi}_{j}-\phi_{j}\|=O_{P}(a_{J}^{-1}n^{-1/2})$ with $a_{J}\coloneqq\min_{j\leq J}\left\{\rho_{j}-\rho_{j+1},\rho_{j-1}-\rho_{j}\right\}$. Here, $\rho_{j}$ are the eigenvalues of the covariance operator of $X$ (placed in descending order).

Theorem 3**.**

If $\operatorname{E}\|X\|^{4}<\infty$, $a_{J}^{-1}n^{-1/2}=o(1)$ and (H(H1))–(H(H11)) hold,

1. (i)

$\operatorname{E}_{\boldsymbol{X}}\left\{\widehat{\widehat{m}}(x)\right\}=m(x)+O_{P}\left(J^{1/2}\,\|\mathcal{P}_{\mathcal{S}_{J}^{\perp}}m_{x}^{\prime}\|\,h\right)+O_{P}\left(J^{1/2}\,h^{2}\right)$

*                       $+\,O_{P}\left(a_{J}^{-1}\,n^{-1/2}\,J^{1/2}\,h\right)$,* 2. (ii)

$\operatorname{Var}_{\boldsymbol{X}}\left\{\widehat{\widehat{m}}(x)\right\}\,=\,O_{P}\left(\{n\,\pi_{x}(h)\}^{-1}\right)$.

Consideration of the data driven basis degrades slightly the conditional bias by introducing in both original terms the quantity $J^{1/2}$ and by adding a third term $O_{P}\left(a_{J}^{-1}\,n^{-1/2}\,J^{1/2}\,h\right)$. However, the conditional variance is not sensitive to the introduction of the data driven basis.

5 Implementation

In this section we discuss the practical aspects, and assess the finite sample performance of our local linear estimators for functional data. Beyond standard issues such as the choice of the tuning parameters in the estimation of the regression operator, a novel heuristic is developed for selection of the bandwidth for the functional derivative. In a comparative simulation study, finite sample performance of our estimator is compared with its competitors available in the literature. We conclude this section with a real data example and an application of the functional local linear estimator. On a benchmark growth dataset we demonstrate a strong link between the considered scalar-on-function local linear regression, and the important single-functional index model, widely considered in the literature.

5.1 Selection of tuning parameters

Methodology. According to the definition of $\widehat{m}(x)$ and $\widehat{m^{\prime}_{x}}$, two parameters have to be selected: the dimension $J$ of the approximating subspace and the bandwidth $h$. Theoretical results of Theorems 1 and 2 emphasize different asymptotic behaviors for the estimators of regression and functional derivative. This is why the optimal parameters for estimating the regression operator do not match necessarily those designed for the functional derivative.

Choosing optimal parameters for the functional derivative is quite challenging because no quantity is directly available to compare with. The following ad hoc methodology is proposed. Firstly, optimal parameters $h_{reg}$ and $J_{reg}$ for the estimator of the regression operator $\widehat{m}$ are selected. The R implementation of the estimating procedure that can be found in the online supplementary material allows for the use of two standard criteria — (leave-one-out) cross-validation [3, 52, 53] and an adaptation of the corrected Akaike information criterion [2, 32] to functional data; a general overview of these criteria can be found in [30]. Secondly, the parameters $h_{reg}$ and $J_{reg}$ are used to build, for each random function $X_{i}$ in the sample, a pilot estimator $\widehat{m_{X_{i}}^{\prime}}^{boot}$ of the functional derivative at $X_{i}$ by means of a wild bootstrap procedure [55, 40, 21]. Then, the optimal parameters $h_{deriv}$ and $J_{deriv}$ are those minimizing the mean squared error

[TABLE]

where $\widehat{m^{\prime}_{{X_{i}},-i}}$ is the estimator of the functional derivative at $X_{i}$ from a dataset with the $i$th observation removed. An important methodological point consists in translating the bandwidths into the number of nearest neighbors. Given a function $x$ in $H$ and an integer $k$, a $k$-nearest neighbors bandwidth at $x$ corresponds to the smallest bandwidth $h$ such that the number of functional regressors $X_{i}$ belonging to ball of radius $h$ centered at $x$ is equal to $k$. An advantage of such a local approach is that it provides more flexible estimators while reducing a continuous set of candidates to a discrete one. From now on, all bandwidths in our results are expressed in terms of nearest neighbors.

Focus on bandwidth choice. In order to assess the quality of the bandwidth choice for both estimators, a first model (M1) is simulated so that the regression operator as well as the functional derivatives can be expressed analytically.

Simulated model (M1). Let $X_{1},\dots,X_{n},X_{n+1},\dots,X_{n+500}$ be iid copies of a functional predictor $X$. The estimators are based on the training set $X_{1},\ldots,X_{n}$; the 500 remaining functions $X_{n+1},\ldots,X_{n+500}$ are used to assess the quality of the estimators. The model takes the form $Y\coloneqq m(X)+\varepsilon$, where $\varepsilon$ is an independent, centered and normally distributed error term with variance $\sigma_{\varepsilon}^{2}$. Let $\phi_{1},\ldots,\phi_{4}$ be the first four elements of the Fourier basis; the random function $X$ is equal to the linear combination $\sum_{j=1}^{4}U_{j}\,\phi_{j}$ where $U_{j}$ are iid uniform random variables on $[-1,1]$. The regression operator is given by $m(X)\coloneqq\sum_{j=1}^{4}\exp(-U_{j}^{2})$. From the expression for $m$, one can derive its functional derivative at $x$ that takes the form $m_{x}^{\prime}(t)=-2\sum_{j=1}^{4}U_{j}\,\exp(-U_{j}^{2})\,\phi_{j}(t)$. Here, the size of the approximating subspace is $J=4$. We consider $n=100,150,200,\ldots,500$. At last, the noise-to-signal ratio ($nsr$) is also controlled by setting $\sigma_{\varepsilon}^{2}\coloneqq nsr\times\mbox{Var}\{m(X)\}$, with $nsr=0.05,0.2,0.4$. Because we focus on the bandwidth selection, in the present model in the estimation procedure we deliberately use the true approximating subspace (i.e. $J=4$ and the first four Fourier basis elements $\phi_{1},\ldots,\phi_{4}$).

Regression operator estimation. The first step is to assess the quality of the bandwidth selection for estimating the regression operator. To this end, the optimal bandwidth $h_{reg}$ is the one minimizing the CV (cross-validation) or AICC (corrected Akaike information criterion). Because in (M1) the true regression operator $m$ is known, one can compute the oracle relative mean squared error of prediction

[TABLE]

where $\overline{m(X)}\coloneqq 1/500\sum_{i=n+1}^{n+500}m(X_{i})$. Repeating the simulation scheme 100 times in various situations, Figure 1 assesses the quality of the local linear estimator of the regression operator. The consistency of the estimator appears clearly even for a large noise-to-signal ratio. Concerning the bandwidth selection, both methods CV and AICC provide similar results, although CV seems to outperform slightly AICC.

Figure 2 displays the true values $m(X_{i})$ against their predictions $\widehat{m}(X_{i})$ when cross-validation is used. Even in the worst case (small learning sample size and high noise-to-signal ratio), the local linear estimator provides reliable results.

Functional derivative estimation. The selection of the bandwidth for estimating the functional derivative is much more challenging because there is no standard criterion to minimize. An original bandwidth selection based on the wild bootstrap procedure is proposed. It aims to build a pilot estimator of the functional derivative:

1. (i)

use $h_{reg}$ for estimating the model error $\widehat{\varepsilon}_{i}\coloneqq Y_{i}-\widehat{m}(X_{i})$ for $i=1,\ldots,n$, 2. (ii)

given iid centered random variables $V_{1},\ldots,V_{n}$ independent of $\widehat{\varepsilon}_{i}$ such that their first moments equal $1$, compute the bootstrapped errors $\varepsilon_{i}^{(b)}\coloneqq\widehat{\varepsilon}_{i}\times V_{i}$, 3. (iii)

derive a bootstrapped sample \mathcal{S}^{(b)}\coloneqq\big{(}X_{i},\ Y_{i}^{(b)}\coloneqq\widehat{m}(X_{i})+\varepsilon_{i}^{(b)}\big{)}_{i=1,\ldots,n} and compute the bootstrapped estimator $\widehat{m^{\prime}_{X_{i}}}^{(b)}$ from $\mathcal{S}^{(b)}$.

Repeat steps (i)–(iii) independently $B$ times and denote $\widehat{m^{\prime}_{X_{i}}}^{boot}\coloneqq 1/B\sum_{b=1}^{B}\widehat{m^{\prime}_{X_{i}}}^{(b)}$ what we name the pilot estimator of the functional derivative at $X_{i}$. The optimal bandwidth $h_{deriv}$ is defined as the one minimizing (6). In order to assess the relevance of this bandwidth choice, a simulation study with model (M1) is conducted with $B=100$. Similarly as in the study of the estimator of the regression operator, Figure 3(a) displays the means of the oracle relative mean squared error

[TABLE]

where $\overline{m^{\prime}_{X}}\coloneqq 1/500\sum_{i=n+1}^{n+500}m^{\prime}_{X_{i}}$. Firstly, selecting $h_{reg}$ with CV or AICC has no significant impact on the prediction quality of the estimator of functional derivatives. Secondly, this plot demonstrates the consistency of the local linear estimator of the functional derivative. Its rate of convergence seems to be slightly slower than the one observed for the local linear estimator of the regression operator, as supported by the asymptotic results. Figure 3(b) reproduces the same plot but $\widehat{m_{X_{i}}^{\prime}}$ is built with $h_{reg}$ (the bandwidth used for estimating the regression operator with cross-validation) instead of $h_{deriv}$ (the specific bandwidth computed for estimating the functional derivative). In this situation, it is the AICC method used for selecting $h_{reg}$ which provides better predictions. Nevertheless, consistency is far less obvious than in part (a). Figures 3(a) and 3(b) illustrate the importance of selecting a specific bandwidth for the estimation of functional derivatives.

Figure 4 compares a sample of true functional derivatives $m_{X_{1}}^{\prime},\ldots,m_{X_{n}}^{\prime}$ with corresponding predictions $\widehat{m_{X_{1}}^{\prime}},\ldots,\widehat{m_{X_{n}}^{\prime}}$ (a) in the worst case ($n=100$ and $nsr=0.4$) and (b) in the most favorable situation ($n=500$ and $nsr=0.05$). Even in the worst situation, the predictions remain adequate.

Another way to assess the bandwidth choice for estimating the functional derivative is to compare the selected bandwidth $h_{deriv}$ itself with the oracle one $h_{deriv}^{oracle}$ minimizing the oracle relative mean squared error $ORMSE_{deriv}$ from (8) computed from the random sample functions $X_{1},\dots,X_{n}$. Given a learning sample size $n$ and a noise-to-signal ratio ($nsr$), 100 simulated datasets are drawn from (M1). Our local linear estimating procedure provides 100 triples of bandwidths $h_{reg}$, $h_{deriv}$ and $h_{deriv}^{oracle}$. Setting $n=100,150,200,\ldots,500$ and $nsr=0.05,0.2,0.4$, Figure 5 displays simultaneously $h_{deriv}^{oracle}$ versus $h_{deriv}$ and $h_{deriv}^{oracle}$ versus $h_{reg}$. For each type of bandwidth, the mean (solid circle for $h_{deriv}$ and square for $h_{reg}$) and standard deviation (whiskers) over the 100 runs are displayed for each of the $9\times 3$ pairs $(n,\,nsr)$. The bandwidth $h_{deriv}$ works quite well even if it underestimates slightly the oracle version $h_{deriv}^{oracle}$. On the other hand, $h_{reg}$ fails drastically, especially for larger sample sizes $n$. Summary statistics of the raw data are provided in Table 6 in Appendix C.

It is worth noting that the bootstrap bandwidth selection introduces additional randomness into the local linear estimation of the functional derivatives. However, Table 5 given in Appendix C indicates clearly that our procedure gives very stable results.

Conclusion. To summarize this section devoted to bandwidth selection, we may conclude that:

(i) cross-validation is a useful method for determining $h_{reg}$, the bandwidth for estimating the local linear regression operator, and (ii) the bootstrap procedure which provides the bandwidth for estimating the functional derivative works well, even for small learning samples size and high noise-to-signal ratios.

Automatic choice of the approximating subspace. So far we focused only on the bandwidth selection. To make our method fully automatic, one has to determine also the approximating subspace spanned by $\phi_{1},\ldots,\phi_{J}$ and its dimension $J$. As explained in Section 4 (case 3), functional principal component analysis is a very useful tool for expanding a random function onto the eigenfunctions of the covariance operator. Let $\phi_{1},\ldots,\phi_{J}$ be the first $J$ eigenfunctions of the covariance operator of the functional predictor $X$ associated to the $J$ largest eigenvalues, and let $\widehat{\phi}_{1},\ldots,\widehat{\phi}_{J}$ be their estimates from the empirical covariance operator. To make estimating procedure fully automatic, we may proceed in three steps:

(i) compute the first $J$ eigenfunctions $\widehat{\phi}_{1},\ldots,\widehat{\phi}_{J}$ of the empirical covariance operator with $J$ large enough, (ii) carry out the local linear estimator of the regression operator with $h_{reg}$ and $J_{opt}$ obtained by minimizing the CV criterion with respect to $J_{opt}\in\left\{0,1,\dots,J\right\}$ and $h_{reg}$, (iii) determine $h_{deriv}$ using the bootstrap procedure by minimizing (6) and compute the corresponding local linear estimator of the functional derivative.

Robustness of the selection procedure for $J_{opt}$. To make the choice of the dimension $J$ more challenging, we add structural perturbation to the functional predictors $X$. Given the first eight Fourier basis elements $\phi_{1},\ldots,\phi_{4},\phi_{5},\ldots,\phi_{8}$, set $X\coloneqq\sum_{j=1}^{4}U_{j}\,\phi_{j}+\eta$ where $\eta\coloneqq\sum_{j=5}^{8}V_{j}\,\phi_{j}$ with $U_{j}$ (resp. $V_{j}$) iid uniform random variables defined on $[-1,\,1]$ (resp. $[-b,\,b]$). The second part $\eta$ provides a structural noise that is controlled by the ratio $\rho\coloneqq\operatorname{E}\left(\|\eta\|^{2}\right)/\operatorname{E}\left(\|X\|^{2}\right)=b^{2}/(1+b^{2})$. We refer to this model as (M2). Given any $\rho\in(0,1)$, one can always find a corresponding bound $b$ for simulating the functional predictors. Table 1 shows how our estimating procedure is robust according to 4 structural noise ratios ($\rho=0.05,\,0.1,\,0.2,\,0.4$) and different learning sample sizes ($n=100,150,\ldots,500$). In each situation, model (M2) is simulated 100 times (with $nsr=0.05$) resulting in 100 estimates $J_{opt}$. The larger $n$ is, the more often the dimension is correctly detected. The results degrade as the ratio $\rho$ increases. The detection is almost perfect for large sample sizes and ratios $\rho$ up to 0.2. Table 2 gives respectively the corresponding $ORMSEP_{reg}$ and $ORMSEP_{deriv}$, each time averaged over 100 runs with standard deviations in brackets. This additional table confirms that our estimating procedure is not too sensitive to perturbations. Complete results of this simulation study which include various choices of the true dimension $J$, $nsr$, and $\rho$, can be found in Appendix C.

Running time. One could think that the introduction of the bootstrap procedure cumulates with the selection of two bandwidths and one dimension, and requires a computation that is quite intensive. Nevertheless, the running time of our R procedure is surprisingly short — at most about 3 seconds (according to the previous simulation scheme with a processor Intel Core i7 2.7 GHz with 16 GB RAM) are necessary to carry out the estimation/prediction for both the regression operator and the functional derivative, including FPCA for the basis expansion $\widehat{\phi}_{1},\ldots,\widehat{\phi}_{J_{opt}}$ with automatic computation of $J_{opt}$, and automatic bandwidths ($h_{reg}$ and $h_{deriv}$) selection, see the last column of Table 1.

5.2 A comparative study

We now conduct a simulation study, in which the finite sample performance of the local linear estimator is compared to its competitors from the literature. In the simulated datasets, we extend models (M1) and (M2) to consider a whole spectrum of scenarios, from a linear to a nonlinear additive one.

Simulated model (M3). A perturbed functional predictor $X\coloneqq\sum_{j=1}^{4}U_{j}\,\phi_{j}+\eta$ is built according to the scheme given in (M2). We now consider $Y\coloneqq m_{a}(X)+\varepsilon$ where $m_{a}(X)\coloneqq(1-a)\langle\beta,X\rangle+a\sum_{j=1}^{4}\exp(-U_{j}^{2})$ and $\beta\coloneqq\sum_{j=1}^{4}\phi_{j}$. Note that the choice $a=0$ corresponds to a standard functional linear model, whereas $a=1$ represents the nonlinear regression model (M2). For all choices of $a$, the model (M3) is additive [42, 43], which means that the regression operator can be expressed as a sum of components where each component is a function that depends only on a single principal score of the regressor $X$. Additivity allows direct computation of the functional derivative of $m_{a}$. The derivative takes the form $m_{a,x}^{\prime}(t)=(1-a)\beta(t)-2a\sum_{j=1}^{4}U_{j}\exp(-U_{j}^{2})\,\phi_{j}(t)$.

For any $X_{j}$ in the testing sample, the predictive performance of our local linear estimators $m(X_{j})$ and $m^{\prime}_{X_{j}}$ is compared with:

• (L)

Functional linear regression estimator: the standard linear regression model applied to the projections of all the involved (centered) functional data into the first $J$ basis functions [50]. The estimator of $m(x)$ is the intercept estimated by this model. A sensible estimator of (the Riesz representation of) the functional derivative can be obtained as $\boldsymbol{\phi}(t)^{\mathsf{T}}\widehat{\boldsymbol{b}_{L}}$, where $\widehat{\boldsymbol{b}_{L}}$ is the estimate of the non-intercept terms in the linear model with $J$ regressors and the intercept. Expansion into the eigenbasis estimated from the random sample functions is considered.

• (LC)

Functional local constant Nadaraya-Watson kernel estimator:

[TABLE]

for $K$ a kernel function, and $h$ a bandwidth [22]. The local constant estimator does not allow direct estimation of the functional derivative $m^{\prime}_{X_{j}}$.

• (LL)

Functional local linear regression estimator: the estimators of $m(X_{j})$ and $m^{\prime}_{X_{j}}$ proposed in this paper (with expansions into the eigenbasis given by the empirical covariance operator of the random sample curves). All parameters (bandwidths and approximating subspace dimension) are automatically selected.

• (MY)

Müller-Yao functional additive model estimator: the (centered) functional data are first projected into the univariate spaces given by their first $J$ estimated eigenfunctions $\widehat{\phi}_{1},\dots,\widehat{\phi}_{J}$ to obtain their principal component scores. For all $j=1,\dots,J$, local polynomial estimates ${\widehat{f_{j}}}$ and ${\widehat{f_{j}}}^{\prime}$ of the regression and its derivative, respectively, in the model of (centered) responses against the univariate scores are obtained. The additive regression operator $m$ is estimated by the sum of the functional values ${\widehat{f_{j}}}$ evaluated at the principal scores of $x$, plus the average of the responses $Y_{i}$. The final estimator of the functional derivative $m_{x}$ is the sum of functions $\widehat{\phi}_{j}$ weighted by the corresponding estimated derivatives ${\widehat{f_{j}}}^{\prime}$ evaluated at the principal scores of $x$. According to the guidelines in [42, 43], $J$ is chosen so that the first $J$ estimated eigenfunctions explain $90~{}\%$ of the variability in the data.

Note that both the kernel estimator (LC) and the linear regression estimator (L) are special cases of the local linear estimator (LL) — for $J=0$ we recover the kernel estimator, and for a kernel $K(t)$ continuous at $t=0$ from the right, the local linear smoother approaches the standard functional linear regression estimator as the bandwidth $h$ tends to infinity.

For all competitors, the asymmetric Epanechnikov kernel $K(t)=0.5(1-t^{2})$ for $t\in[0,\,1]$ is used. The bandwidths, as well as the dimension $J$ in the functional (local) linear regression, are chosen by a leave-one-out cross-validation procedure.

Assessing performances. The learning and testing sample sizes are set to 500. Two perturbations are considered as in Section 5.1: the noise-to-signal ratio $nsr$ of the regression model and the structural perturbation $\rho$ acting on regressors. Parameters $(nsr,\,\rho)$ are set to $(0.05,\,0.05)$ and $(0.4,\,0.4)$, corresponding to a low/high perturbation level and $a\,=\,0,\,0.25,\,0.5$, $0.75,\,1$ successively. Results for other combinations of perturbation levels are provided in Appendix C. 100 runs are performed in each case. To assess the prediction quality, we use the $ORMSEP$ criteria from (7) and (8), except for $a=0$ and the derivative, where the denominator of $ORMSEP_{deriv}$ is null since $m^{\prime}_{0,x}\equiv\beta$ for any $x$. In the latter case, we report only the numerator of $ORMSEP_{deriv}$ from (8). Mean and standard deviation (in brackets) can be found in Tables 3 and 4, each table corresponding to a particular noise level.

From the results of the simulation study we conclude the following:

(i) the functional linear estimator is the best method when the model is linear, or close to linear. In the situation when the model is strongly nonlinear, the estimator fails as expected. (ii) The local linear estimator of the regression operator convincingly outperforms the local constant estimator in all considered scenarios. (iii) For higher values of $a$ which correspond to models far from linear, (MY) performs worse than the local linear estimator, for both regression operator and functional derivatives. This corroborates the good finite sample properties of the local linear estimator observed before, as all models considered in (M3) in the simulation study satisfy the additivity condition, under which (MY) was designed. Note that for (LL) additivity of the regression operator is not required. (iv) Another practical issue regarding the behavior of the estimators is their numerical stability. In the complete results of this simulation study, given in Appendix C, we observed that (MY) tends to be numerically unstable, especially for small learning sample sizes. The instabilities occur mostly when the principal scores of a predictor lie outside the range of the scores of the data, in which case either the functional values $f_{j}$, or the derivatives $f_{j}^{\prime}$ have to be extrapolated. Remarkably, (LL) does not appear to suffer from such drawbacks.

5.3 Benchmark growth data analysis

Berkeley growth data trace back to the pioneering work [54] and were reconsidered in [23, 36, 45, 24]. Recently, in [28] an interesting analysis of this dataset involving estimated functional derivatives was performed. To better understand the growth mechanism, the relationship between the growth velocity profile up to 10 years of age (the functional predictor $X$) and the adult height (observed at 18 years, scalar response $Y$) of the boys (39 individuals) was investigated. Here, we consider the same problem using the local linear methodology. Our approach allows to estimate the functional derivatives corresponding to individual regressors $X$ directly, which greatly facilitates the interpretation of the results. In Figure 6(a) we see the growth velocity profiles obtained via standard univariate local linear regression and in Figure 6(b) the estimated functional derivatives $\widehat{m^{\prime}_{X_{1}}},\ldots,\widehat{m^{\prime}_{X_{39}}}$ are displayed. Focusing on the estimated functional derivatives, a sharp increase at around 6 years is observed for all boys. A possible interpretation is that the growth velocity profile prior to the age of 6 years has little impact on the adult height of an individual. Compared to the previous analyses of the growth dataset, this finding appears to be original.

To better understand the shape of the estimated functional derivatives, we propose to model the relationship between the adult height at 18 and the growth velocity up to 10 using a single-functional index model

[TABLE]

where the link function $g$ and the functional direction $\beta$ are unknown. The single-functional index model is well suited for the studied problem, as all the estimated functional derivatives share a common shape. Therefore, the average functional derivative is a good representative of the collection of the estimated derivatives. Using the average derivative estimation method described in the introduction, we can estimate the functional parameter $\beta$ and the link function $g$ as follows: 1) $\widehat{\operatorname{E}m^{\prime}_{X}}\coloneqq 1/39\sum_{i=1}^{39}\widehat{m^{\prime}_{X_{i}}}$ and $\widehat{\beta}=\widehat{\operatorname{E}m^{\prime}_{X}}/\|\widehat{\operatorname{E}m^{\prime}_{X}}\|$, 2) based on the sample $(Z_{1},Y_{1}),\ldots,(Z_{39},Y_{39})$ where $Z_{i}\coloneqq\langle X_{i},\widehat{\beta}\rangle$, one gets an estimator $\widehat{g}$ of the link function $g$ by any standard univariate nonparametric regression method.

Figures 7(a) and (b) display respectively $\widehat{g}$ and $\widehat{\beta}$. The shape of the estimated functional index $\widehat{\beta}$ reflects the significant jump at around 6 years. The results are quite positive; the estimated heights are strongly correlated with the observed adult heights (Pearson’s correlation $\simeq 0.82$). Since $\widehat{g}$ is a positive and non-decreasing function, the growth velocity after 6 years of age plays a major role in the prediction of the adult height. In order to confirm this interpretation, the single-functional index model was estimated again, but by restricting the growth velocity to the range of 6–10 years of age. As expected, the quality of the estimation is comparable; similar correlation ($\simeq 0.82$) between the estimated adult heights and the observed ones is obtained. For additional details we refer to Appendix C.

Acknowledgments

The work of S. Nagy was supported by the grant 19-16097Y of the Czech Science Foundation, and by the PRIMUS/17/SCI/3 project of Charles University.

\appendixpage

Throughout the appendices, we use several conventions. Sums indicated by $\sum_{i}$ are always meant with $i$ from $1$ to $n$, and sums indicated by $\sum_{j\geq J}$ mean sums with $j$ from $J$ to infinity, for $J$ given; $\boldsymbol{1}$ stands for a column vector of ones of appropriate dimension; $1_{[0,u]}(t)$ is the indicator of $t\in[0,u]$, i.e. $1$ if $t\in[0,u]$ and [math] otherwise. The $(j,k)$-th element of a matrix $\boldsymbol{\Delta}$ can be denoted either by $\left[\boldsymbol{\Delta}\right]_{jk}$, or equivalently by $\Delta_{jk}$.

Appendix A Details of proofs

Proof of Lemma 1-$(i)$.

We know that $T_{1}=\boldsymbol{e}^{\mathsf{T}}\,\left(\boldsymbol{\Phi}^{\mathsf{T}}\widetilde{\boldsymbol{K}}\boldsymbol{\Phi}\right)^{-1}\left[A_{0},A_{1},\dots,A_{J}\right]^{\mathsf{T}}$ with

[TABLE]

for $i=1,\ldots,n$, $K_{i}\coloneqq K\left(h^{-1}\|X_{i}-x\|\right)$ and $\widetilde{\boldsymbol{K}}\coloneqq(n\operatorname{E}K_{1})^{-1}\boldsymbol{K}$. Set $\delta_{0}\coloneqq(n\,\operatorname{E}K_{1})^{-1}\sum_{i}K_{i}$, $\delta_{j}\coloneqq(n\,\operatorname{E}K_{1})^{-1}\sum_{i}\langle\phi_{j},X_{i}-x\rangle K_{i}$ for $j=1,\ldots,J$ and let $\boldsymbol{\Delta}$ be the $J\times J$ matrix whose the $(j,k)$-th entry is equal to $\Delta_{jk}\coloneqq(n\,\operatorname{E}K_{1})^{-1}\sum_{i}\langle\phi_{j},X_{i}-x\rangle\langle\phi_{k},X_{i}-x\rangle K_{i}$ for $j,k=1,\ldots,J$. By using elementary linear algebra,

[TABLE]

where $\boldsymbol{\delta}\coloneqq[\delta_{1},\dots,\delta_{J}]^{\mathsf{T}}$. According to standard results with respect to the inverse of a 2$\times$2 block matrix (see for instance [39]),

[TABLE]

with $\mu\coloneqq\left(\delta_{0}-\boldsymbol{\delta}^{\mathsf{T}}\boldsymbol{\Delta}^{-1}\boldsymbol{\delta}\right)^{-1}$. Then, $T_{1}\,=\,\mu\left(A_{0}\,-\,\boldsymbol{\delta}^{\mathsf{T}}\boldsymbol{\Delta}^{-1}\left[A_{1},\dots,A_{J}\right]^{\mathsf{T}}\right)$. Before going on, let us focus on $\boldsymbol{\Delta}^{-1}$, the inverse of $\boldsymbol{\Delta}$. From Lemma 7, one has $\boldsymbol{\Delta}=\boldsymbol{\Delta}_{1}+\boldsymbol{\Delta}_{2}$ with $\displaystyle\boldsymbol{\Delta}_{1}\coloneqq b_{x,0,1}^{-1}\,b_{x,2,1}\,h^{2}\,\boldsymbol{\Gamma}\,\left\{1+o(1)\right\}$ and $\displaystyle\boldsymbol{\Delta}_{2}\coloneqq O_{P}\left(h^{2}\{n\pi_{x}(h)\}^{-1/2}\right)\boldsymbol{\Lambda}$. As soon as $\boldsymbol{\Delta}_{1}^{-1}$ is invertible, $\boldsymbol{\Delta}^{-1}=\boldsymbol{\Delta}_{1}^{-1}(\boldsymbol{I}+\boldsymbol{\Delta}_{2}\boldsymbol{\Delta}_{1}^{-1})^{-1}$. Let $\|\cdot\|_{F}$ stand for the Frobenius matrix norm and recall that $\lambda_{J}$ is the smallest eigenvalue of the $J\times J$ matrix $\boldsymbol{\Gamma}$. We have that $\|\boldsymbol{\Delta}_{1}^{-1}\|_{F}=O_{P}\left(\lambda_{J}^{-1}\,J^{1/2}\,h^{-2}\right)$, and thanks to Lemma 6, $\|\boldsymbol{\Delta}_{2}\|_{F}=O_{P}\left(h^{2}\,\{n\pi_{x}(h)\}^{-1/2}\right)$ so that $\|\boldsymbol{\Delta}_{2}\boldsymbol{\Delta}_{1}^{-1}\|_{F}=O_{P}\left(\lambda_{J}^{-1}\,J^{1/2}\,\{n\pi_{x}(h)\}^{-1/2}\right)$. According to (H(H7)), for $n$ large enough, the Frobenius norm of $\boldsymbol{\Delta}_{2}\boldsymbol{\Delta}_{1}^{-1}$ is smaller than 1. Then $\boldsymbol{\Delta}^{-1}=\boldsymbol{\Delta}_{1}^{-1}$ $\left\{\boldsymbol{I}+\sum_{k\geq 1}(-1)^{k}\left(\boldsymbol{\Delta}_{2}\boldsymbol{\Delta}_{1}^{-1}\right)^{k}\right\}$ which results in

[TABLE]

Lemmas 8 and 9 allow us to write

[TABLE]

where we denote $\boldsymbol{\gamma}\coloneqq[{\gamma_{1}^{1}}^{{}^{\prime}}(0),\dots,{\gamma_{J}^{1}}^{{}^{\prime}}(0)]^{\mathsf{T}}$, $\boldsymbol{\theta}\coloneqq\left[\sqrt{{\gamma_{1}^{2}}^{{}^{\prime}}(0)},\dots,\sqrt{{\gamma_{J}^{2}}^{{}^{\prime}}(0)}\right]^{\mathsf{T}}$, $\boldsymbol{\alpha_{x,n}^{bias}}\coloneqq\left[\alpha_{1,x,n}^{bias},\dots,\alpha_{J,x,n}^{bias}\right]^{\mathsf{T}}$ and $\sqrt{\boldsymbol{\alpha_{x,n}^{var}}}\coloneqq\left[\sqrt{\alpha_{1,x,n}^{var}},\dots,\sqrt{\alpha_{J,x,n}^{var}}\right]^{\mathsf{T}}$. From Lemma 9 we know that $\|\boldsymbol{\alpha_{x,n}^{bias}}\|_{2}\leq\|\mathcal{P}_{\mathcal{S}_{J}^{\perp}}m_{x}^{\prime}\|$ and $\|\sqrt{\boldsymbol{\alpha_{x,n}^{var}}}\|_{2}\leq\|\mathcal{P}_{\mathcal{S}_{J}^{\perp}}m_{x}^{\prime}\|$. Here $\|.\|_{2}$ stands for the Euclidean vector norm. Based on Lemmas 6, 9, and 10,

[TABLE]

where $\lambda_{J}$ is the smallest eigenvalue of the $J\times J$ matrix $\boldsymbol{\Gamma}$. Lemma 6 indicates that the Frobenius norm of $\boldsymbol{\Gamma}$ is finite and thus $\lambda_{J}$ converges to 0 with $n$. Let us now focus on the term $\mu$. Firstly, one has $\operatorname{E}\delta_{0}=1$ and $\operatorname{Var}(\delta_{0})=n^{-1}(\operatorname{E}K_{1})^{-2}\operatorname{Var}(K_{1})$. Thanks to Corollary 1, it is easy to see that $\operatorname{Var}(\delta_{0})=O\left(1/\{n\,\pi_{x}(h)\}\right)$ which leads to $\delta_{0}=1+O_{P}\left(1/\sqrt{n\,\pi_{x}(h)}\right)$. Secondly, Lemma 8 and (10) imply that

[TABLE]

Lemmas 6 and 10 give $\boldsymbol{\theta}^{\mathsf{T}}\boldsymbol{\Gamma}^{-1}\boldsymbol{\gamma}\,/\,\boldsymbol{\gamma}^{\mathsf{T}}\boldsymbol{\Gamma}^{-1}\boldsymbol{\gamma}=O(\lambda_{J}^{-1})=\boldsymbol{\theta}^{\mathsf{T}}\boldsymbol{\Gamma}^{-1}\boldsymbol{\theta}\,/\,\boldsymbol{\gamma}^{\mathsf{T}}\boldsymbol{\Gamma}^{-1}\boldsymbol{\gamma}$. Then (H(H7)) leads to $\{n\,\pi_{x}(h)\}^{-1/2}\boldsymbol{\theta}^{\mathsf{T}}\boldsymbol{\Gamma}^{-1}\boldsymbol{\gamma}=o_{P}(\boldsymbol{\gamma}^{\mathsf{T}}\boldsymbol{\Gamma}^{-1}\boldsymbol{\gamma})=\{n\,\pi_{x}(h)\}^{-1}\boldsymbol{\theta}^{\mathsf{T}}\boldsymbol{\Gamma}^{-1}\boldsymbol{\theta}$ and finally, $\boldsymbol{\delta}^{\mathsf{T}}\boldsymbol{\Delta}^{-1}\boldsymbol{\delta}\,=\,b_{x,0,1}^{-1}\,b_{x,1,1}^{2}\,b_{x,2,1}^{-1}\,\boldsymbol{\gamma}^{\mathsf{T}}\boldsymbol{\Gamma}^{-1}\boldsymbol{\gamma}\,\left\{1+o_{P}(1)\right\}$. Use again Lemma 10 to get

[TABLE]

This last result combined with (11) and (12) gives $T_{1}\,=\,O_{P}\left(h\,\|\mathcal{P}_{\mathcal{S}_{J}^{\perp}}m_{x}^{\prime}\|\right)$. ∎

Proof of Lemma 1-$(ii)$.

We have that $T_{2}=T_{21}+T_{22}$ where

• $\bullet$

$T_{21}\coloneqq\boldsymbol{e}^{\mathsf{T}}\,\left(\boldsymbol{\Phi}^{\mathsf{T}}\widetilde{\boldsymbol{K}}\boldsymbol{\Phi}\right)^{-1}\left[B_{0},B_{1},\dots,B_{J}\right]^{\mathsf{T}}$ with $B_{0}\coloneqq(n\operatorname{E}K_{1})^{-1}\sum_{i}K_{i}R_{x,x,i}$, and for $j=1,\ldots,J$, $B_{j}\coloneqq(n\operatorname{E}K_{1})^{-1}\sum_{i}K_{i}R_{x,x,i}\langle\phi_{j},X_{i}-x\rangle$, where $R_{x,x,i}$ is a term involved in (2) with $\zeta$ replaced by $x$;

• $\bullet$

$T_{22}\coloneqq\boldsymbol{e}^{\mathsf{T}}\,\left(\boldsymbol{\Phi}^{\mathsf{T}}\widetilde{\boldsymbol{K}}\boldsymbol{\Phi}\right)^{-1}\left[C_{0},C_{1},\dots,C_{J}\right]^{\mathsf{T}}$ where $C_{0}\coloneqq(n\operatorname{E}K_{1})^{-1}\sum_{i}K_{i}\left(R_{\zeta,x,i}\right.$ $-\left.R_{x,x,i}\right)$, and $C_{j}\coloneqq(n\operatorname{E}K_{1})^{-1}\sum_{i}K_{i}\left(R_{\zeta,x,i}-R_{x,x,i}\right)\langle\phi_{j},X_{i}-x\rangle$ for $j=1,\ldots,J$.

About $T_{21}$. Lemmas 8 and (11) allow to write

[TABLE]

Based on Lemmas 6, 10 and 11,

[TABLE]

Now, (13), (14), and (15) result in $T_{21}\,=\,O_{P}\left(h^{2}\right)$.

About $T_{22}$. Thanks to (H(H1)), it is easy to show that $C_{0}=O_{P}(h^{3})$ and $C_{j}=O_{P}(h^{4})$ for $j=1,\ldots,J$. Consequently,

[TABLE]

According to (13) and Lemmas 7, 8 and 10, $T_{22}=O_{P}\left(h^{3}\sqrt{J}\right)$.

Back to $T_{2}$. Because $h\,\sqrt{J}=o(1)$ thanks to (H(H7)), $T_{22}=o_{P}(h^{2})$ and $T_{2}=O_{P}(h^{2})$. ∎

Let us now focus on the conditional variance of $\widehat{m}(x)$. With the notations introduced in the proof of Lemma 1-$(i)$, the next results provide a decomposition of the conditional variance with the asymptotic behavior of each term.

Lemma 3**.**

As soon as conditions (H(H1))–(H(H11)) are fulfilled,

1. (i)

$\operatorname{Var}_{\boldsymbol{X}}\left\{\widehat{m}(x)\right\}=\{\sigma^{2}(x)\,+\,o(1)\}\,\mu^{2}\,\left\{\widetilde{\delta}_{0}-\boldsymbol{\delta}^{\mathsf{T}}\boldsymbol{\Delta}^{-1}\widetilde{\boldsymbol{\delta}}-\widetilde{\boldsymbol{\delta}}^{\mathsf{T}}\boldsymbol{\Delta}^{-1}\boldsymbol{\delta}+\boldsymbol{\delta}^{\mathsf{T}}\boldsymbol{\Delta}^{-1}\widetilde{\boldsymbol{\Delta}}\boldsymbol{\Delta}^{-1}\boldsymbol{\delta}\right\},$

where $\widetilde{\delta}_{0}\coloneqq(n\,\operatorname{E}K_{1})^{-2}\,\sum_{i}K_{i}^{2}$, $\widetilde{\delta}_{j}\coloneqq(n\,\operatorname{E}K_{1})^{-2}\sum_{i}\langle\phi_{j},X_{i}-x\rangle\,K_{i}^{2}$ for $j=1,\ldots,J$, $\widetilde{\boldsymbol{\delta}}\coloneqq[\widetilde{\delta}_{1},\dots,\widetilde{\delta}_{J}]^{\mathsf{T}}$, $\widetilde{\Delta}_{jk}\coloneqq(n\,\operatorname{E}K_{1})^{-2}\sum_{i}\langle\phi_{j},X_{i}-x\rangle\,\langle\phi_{k},X_{i}-x\rangle\,K_{i}^{2}$ for $j,k=1,\ldots,J$, and $\widetilde{\boldsymbol{\Delta}}$ is the $J\times J$ matrix such that $[\widetilde{\boldsymbol{\Delta}}]_{jk}\coloneqq\widetilde{\Delta}_{jk}$, 2. (ii)

$\widetilde{\delta}_{0}=b_{x,0,1}^{-2}\,b_{x,0,2}\,\{n\,\pi_{x}(h)\}^{-1}\{1+o_{P}(1)\},$* $\boldsymbol{\delta}^{\mathsf{T}}\boldsymbol{\Delta}^{-1}\widetilde{\boldsymbol{\delta}}=O_{P}\left(\,\{\lambda_{J}\,n\,\pi_{x}(h)\}^{-1}\right),$ and $\boldsymbol{\delta}^{\mathsf{T}}\boldsymbol{\Delta}^{-1}\widetilde{\boldsymbol{\Delta}}\boldsymbol{\Delta}^{-1}\boldsymbol{\delta}=O_{P}\left(\{\lambda_{J}\,n\,\pi_{x}(h)\}^{-1}\right)$.*

Proof of Lemma 3-$(i)$.

The main term of $\operatorname{Var}_{\boldsymbol{X}}\left\{\widehat{m}(x)\right\}$ is given by

[TABLE]

where \boldsymbol{\Phi}^{\mathsf{T}}\widetilde{\boldsymbol{K}}^{2}\boldsymbol{\Phi}\,=\,\left[\begin{array}[]{c|c}\widetilde{\delta}_{0}&\widetilde{\boldsymbol{\delta}}^{\mathsf{T}}\\ \hline\cr&\\[-10.00002pt] \widetilde{\boldsymbol{\delta}}&\widetilde{\boldsymbol{\Delta}}\end{array}\right]. Formula (9) with elementary linear algebra result in the claimed decomposition. ∎

Proof of Lemma 3-$(ii)$.

By Corollary 1, $\operatorname{E}K_{1}^{2}=b_{x,0,2}\,\pi_{x}(h)\,\{1+o(1)\}$ and in the same vein of Lemmas 7 and 8, it is easy to see that

[TABLE]

As a by-product, one has

[TABLE]

and

[TABLE]

These two last equations with (10) lead to

[TABLE]

and

[TABLE]

Combination of Lemmas 6-$(ii)$ and 10 gives $\boldsymbol{\gamma}^{\mathsf{T}}\boldsymbol{\Gamma}^{-1}\boldsymbol{\Lambda}\boldsymbol{\Gamma}^{-1}\boldsymbol{\gamma}=O(\lambda_{J}^{-2})$, $\boldsymbol{\gamma}^{\mathsf{T}}\boldsymbol{\Gamma}^{-1}\boldsymbol{\gamma}=O(\lambda_{J}^{-1})$ and the claimed assertions hold. Finally, the rate of convergence of the conditional variance is derived from Lemma 3 and the fact that $\mu=O_{P}(\lambda_{J})$. ∎

Lemma 4**.**

As soon as conditions (H(H1))–(H(H9)) are fulfilled,

1. (i)

$\|Q_{1}\|\,=\,O_{P}\left(\lambda_{J}^{-1}\|\mathcal{P}_{\mathcal{S}_{J}^{\perp}}m_{x}^{\prime}\|\right)$, 2. (ii)

$\|Q_{2}\|\,=\,O_{P}\left(\lambda_{J}^{-1}\,h\right)$.

Proof of Lemma 4-$(i)$.

With the notations introduced in the proof Lemma 1-$(i)$ and (9), $Q_{1}=-\mu\,A_{0}\,\boldsymbol{\phi}^{\mathsf{T}}\boldsymbol{\Delta}^{-1}\boldsymbol{\delta}+\boldsymbol{\phi}^{\mathsf{T}}(\boldsymbol{\Delta}^{-1}+\mu\,\boldsymbol{\Delta}^{-1}\,\boldsymbol{\delta}\,\boldsymbol{\delta}^{\mathsf{T}}\,\boldsymbol{\Delta}^{-1})\left[A_{1},\dots,A_{J}\right]^{\mathsf{T}}$. According to Lemma 8, Lemma 9 and (10), $\mu\,A_{0}\,\boldsymbol{\phi}^{\mathsf{T}}\boldsymbol{\Delta}^{-1}\boldsymbol{\delta}\,=\,O_{P}\left(\lambda_{J}\|\mathcal{P}_{\mathcal{S}_{J}^{\perp}}m_{x}^{\prime}\|\,\boldsymbol{\phi}^{\mathsf{T}}\boldsymbol{\Gamma}^{-1}\boldsymbol{\gamma}\right)$, $\boldsymbol{\phi}^{\mathsf{T}}\boldsymbol{\Delta}^{-1}\left[A_{1},\dots,A_{J}\right]^{\mathsf{T}}=O_{P}\left(\boldsymbol{\phi}^{\mathsf{T}}\boldsymbol{\Gamma}^{-1}\boldsymbol{\alpha_{x,n}^{bias}}\right)$ and $\mu\,\boldsymbol{\phi}^{\mathsf{T}}\boldsymbol{\Delta}^{-1}\,\boldsymbol{\delta}\,\boldsymbol{\delta}^{\mathsf{T}}\,\boldsymbol{\Delta}^{-1}\left[A_{1},\dots,A_{J}\right]^{\mathsf{T}}=$ $O_{P}\left(\|\mathcal{P}_{\mathcal{S}_{J}^{\perp}}m_{x}^{\prime}\|\,\boldsymbol{\phi}^{\mathsf{T}}\boldsymbol{\Gamma}^{-1}\boldsymbol{\gamma}\right)$. Set $\boldsymbol{u}\coloneqq\lambda_{J}\boldsymbol{\Gamma}^{-1}\boldsymbol{\gamma}$. For any $J$, $\|\boldsymbol{\phi}^{\mathsf{T}}\boldsymbol{\Gamma}^{-1}\boldsymbol{\gamma}\|=\lambda_{J}^{-1}\|\boldsymbol{\phi}^{\mathsf{T}}\boldsymbol{u}\|=\lambda_{J}^{-1}\|\boldsymbol{u}\|_{2}\leq\lambda_{J}^{-1}$. Consequently $\|\boldsymbol{\phi}^{\mathsf{T}}\boldsymbol{\Gamma}^{-1}\boldsymbol{\gamma}\|=O_{P}\left(\lambda_{J}^{-1}\right)$. Now, define the vector $\boldsymbol{v}\coloneqq\lambda_{J}\|\mathcal{P}_{\mathcal{S}_{J}^{\perp}}m_{x}^{\prime}\|^{-1}\boldsymbol{\Gamma}^{-1}\boldsymbol{\alpha_{x,n}^{bias}}$ with $\|\boldsymbol{v}\|_{2}\leq 1$. From the definition of $\boldsymbol{v}$ we have $\|\boldsymbol{\phi}^{\mathsf{T}}\boldsymbol{\Gamma}^{-1}\boldsymbol{\alpha_{x,n}^{bias}}\|=\lambda_{J}^{-1}\|\mathcal{P}_{\mathcal{S}_{J}^{\perp}}m_{x}^{\prime}\|\,\|\boldsymbol{\phi}^{\mathsf{T}}\boldsymbol{v}\|$. Because of the orthonormality of the basis, $\|\boldsymbol{\phi}^{\mathsf{T}}\boldsymbol{v}\|=\|\boldsymbol{v}\|_{2}$, $\|\boldsymbol{\phi}^{\mathsf{T}}\boldsymbol{\Gamma}^{-1}\boldsymbol{\alpha_{x,n}^{bias}}\|=\lambda_{J}^{-1}\|\mathcal{P}_{\mathcal{S}_{J}^{\perp}}m_{x}^{\prime}\|\,\|\boldsymbol{v}\|_{2}\leq\lambda_{J}^{-1}\|\mathcal{P}_{\mathcal{S}_{J}^{\perp}}m_{x}^{\prime}\|$ which gives $\|\boldsymbol{\phi}^{\mathsf{T}}\boldsymbol{\Gamma}^{-1}\boldsymbol{\alpha_{x,n}^{bias}}\|=O_{P}\left(\lambda_{J}^{-1}\|\mathcal{P}_{\mathcal{S}_{J}^{\perp}}m_{x}^{\prime}\|\right)$ and the claimed rate of convergence holds. ∎

Proof of Lemma 4-$(ii)$.

Similarly to the notations and the proof of Lemma 1-$(ii)$, $Q_{2}=Q_{21}+Q_{22}$ with

[TABLE]

By following the proof of Lemma 4-$(i)$ and replacing Lemma 9 with Lemma 11, $\mu\,B_{0}\,\|\boldsymbol{\phi}^{\mathsf{T}}\boldsymbol{\Delta}^{-1}\boldsymbol{\delta}\|=O_{P}\left(h\right)$, $\|\boldsymbol{\phi}^{\mathsf{T}}\boldsymbol{\Delta}^{-1}\left[B_{1},\dots,B_{J}\right]^{\mathsf{T}}\|=O_{P}\left(\lambda_{J}^{-1}\,h\right)$ and $\mu\,\|\boldsymbol{\phi}^{\mathsf{T}}\boldsymbol{\Delta}^{-1}\,\boldsymbol{\delta}\,\boldsymbol{\delta}^{\mathsf{T}}\,\boldsymbol{\Delta}^{-1}\left[B_{1},\dots,B_{J}\right]^{\mathsf{T}}\|=O_{P}\left(\lambda_{J}^{-1}\,h\right)$ so that $\|Q_{21}\|=O_{P}\left(\lambda_{J}^{-1}\,h\right)$. For studying $Q_{22}$, recall that $C_{0}=O_{P}(h^{3})$ and $C_{j}=O_{P}(h^{4})$ for $j=1,\ldots,J$ (see again the proof of Lemma 1-$(ii)$). Then, we obtain that $\mu\,C_{0}\,\|\boldsymbol{\phi}^{\mathsf{T}}\boldsymbol{\Delta}^{-1}\boldsymbol{\delta}\|=O_{P}\left(h^{2}\right)$, $\|\boldsymbol{\phi}^{\mathsf{T}}\boldsymbol{\Delta}^{-1}\left[C_{1},\dots,C_{J}\right]^{\mathsf{T}}\|=O_{P}\left(\lambda_{J}^{-1}\,\sqrt{J}\,h^{2}\right)$, $\mu\,\|\boldsymbol{\phi}^{\mathsf{T}}\boldsymbol{\Delta}^{-1}\,\boldsymbol{\delta}\,\boldsymbol{\delta}^{\mathsf{T}}\,\boldsymbol{\Delta}^{-1}\left[C_{1},\dots,C_{J}\right]^{\mathsf{T}}\|=$ $O_{P}\left(\lambda_{J}^{-1}\,\sqrt{J}\,h^{2}\right)$ and $\|Q_{22}\|=o_{P}\left(\lambda_{J}^{-1}\,h\right)$ from (H(H7)). This is enough to get the claimed result. ∎

Details of the proof of Lemma 2-$(ii)$.

Set $\boldsymbol{M}_{0}\coloneqq\boldsymbol{\Delta}^{-1}\boldsymbol{\delta}\boldsymbol{\delta}^{\mathsf{T}}\boldsymbol{\Delta}^{-1}$. Thanks to (9), elementary linear algebra and notations introduced in Lemma 3, the following decomposition of the conditional variance

[TABLE]

holds with

[TABLE]

According to (9), (13), (16), (17), (18), and (19), Lemmas 3, 7, 8, and 10, with the fact that $\phi_{1},\phi_{2},\ldots$ is an orthonormal basis,

[TABLE]

and $\int_{0}^{1}\boldsymbol{\phi}(t)^{\mathsf{T}}\boldsymbol{M}_{8}\boldsymbol{\phi}(t)\,\mathrm{d}\,t=O_{P}\left(\lambda_{J}^{-1}u_{n}\right)$, where $u_{n}\coloneqq h^{-2}\{n\pi_{x}(h)\}^{-1}$. In conclusion, $\operatorname{E}_{\boldsymbol{X}}\left(\|\widehat{m}_{x}^{\prime}-\operatorname{E}_{\boldsymbol{X}}\left(\widehat{m}_{x}^{\prime}\right)\|^{2}\right)\,=\,O_{P}\left(\lambda_{J}^{-2}u_{n}\right)+O_{P}\left(J\lambda_{J}^{-1}u_{n}\right)$ which is exactly the claimed rate of convergence. ∎

Proof of Theorem 3.

Let us start with the conditional bias. Similarly to (3), $\operatorname{E}_{\boldsymbol{X}}\left\{\widehat{\widehat{m}}(x)\right\}=m(x)+\widehat{T}_{1}+\widehat{T}_{2}+\widehat{T}_{3}$ where $\widehat{T}_{1}$ (resp. $\widehat{T}_{2}$) corresponds to $T_{1}$ (resp. $T_{2}$) when replacing $\boldsymbol{\Phi}$ with $\widehat{\boldsymbol{\Phi}}$, and $\widehat{T}_{3}\coloneqq\boldsymbol{e}^{\mathsf{T}}\,\left(\widehat{\boldsymbol{\Phi}}^{\mathsf{T}}\boldsymbol{K}\widehat{\boldsymbol{\Phi}}\right)^{-1}\widehat{\boldsymbol{\Phi}}^{\mathsf{T}}\boldsymbol{K}(\boldsymbol{\Phi}-\widehat{\boldsymbol{\Phi}})\left[m(x)|\nabla{m_{x}}^{\mathsf{T}}\right]^{\mathsf{T}}$.

About $\widehat{T}_{1}$. Similarly to the proof of Lemma 1-$(i)$, $\widehat{T}_{1}=\widehat{\mu}\left(A_{0}-\widehat{\boldsymbol{\delta}}^{\mathsf{T}}\widehat{\boldsymbol{\Delta}}^{-1}[\widehat{A}_{1}\ldots\widehat{A}_{J}]^{\mathsf{T}}\right)$ where $\widehat{\boldsymbol{\delta}}$ (resp. $\widehat{\boldsymbol{\Delta}}$, $\widehat{A}_{1},\ldots,\widehat{A}_{J}$) are defined as $\boldsymbol{\delta}$ (resp. $\boldsymbol{\Delta}$, $A_{1},\ldots,A_{J}$) but $\phi_{1},\ldots,\phi_{J}$ are replaced with $\widehat{\phi}_{1},\ldots,\widehat{\phi}_{J}$ and where $\widehat{\mu}\coloneqq\left(\delta_{0}-\widehat{\boldsymbol{\delta}}^{\mathsf{T}}\widehat{\boldsymbol{\Delta}}^{-1}\widehat{\boldsymbol{\delta}}\right)^{-1}$. Because $A_{0}=O_{P}\left(\|\mathcal{P}_{\mathcal{S}_{J}^{\perp}}m_{x}^{\prime}\|\,h\right)$ and $\widehat{A}_{j}=O_{P}\left(\|\mathcal{P}_{\mathcal{S}_{J}^{\perp}}m_{x}^{\prime}\|\,h^{2}\right)$ for $j=1,\ldots,J$,

[TABLE]

According to the definitions of $\widehat{\boldsymbol{\delta}}$ and $\widehat{\boldsymbol{\Delta}}$, it is easy to see that $\widehat{\boldsymbol{\delta}}=\boldsymbol{\delta}\left(1+U_{n}\right)$ and $\widehat{\boldsymbol{\Delta}}^{-1}=\boldsymbol{\Delta}^{-1}\left(1+U_{n}\right)$ with $U_{n}\coloneqq O_{P}\left(a_{J}^{-1}\,n^{-1/2}\right)$ being the rate of convergence of $\max_{j\leq J}\|\phi_{j}-\widehat{\phi}_{j}\|$. On one side, $\widehat{\boldsymbol{\delta}}^{\mathsf{T}}\widehat{\boldsymbol{\Delta}}^{-1}\widehat{\boldsymbol{\delta}}=\boldsymbol{\delta}^{\mathsf{T}}\boldsymbol{\Delta}^{-1}\boldsymbol{\delta}\left(1+U_{n}\right)$, and on the other side, $\widehat{\boldsymbol{\delta}}^{\mathsf{T}}\widehat{\boldsymbol{\Delta}}^{-1}\boldsymbol{1}=\boldsymbol{\delta}^{\mathsf{T}}\boldsymbol{\Delta}^{-1}\boldsymbol{1}\left(1+U_{n}\right)$. According to (13), $\boldsymbol{\delta}^{\mathsf{T}}\boldsymbol{\Delta}^{-1}\boldsymbol{\delta}=O_{P}\left(\lambda_{J}^{-1}\right)$ and Lemmas 7, 8 and 10 result in $\widehat{\boldsymbol{\delta}}^{\mathsf{T}}\widehat{\boldsymbol{\Delta}}^{-1}\boldsymbol{1}=O_{P}\left(h^{-1}\,\lambda_{J}^{-1}\,J^{1/2}\right)$. Finally, $\widehat{\mu}=O_{P}\left(\lambda_{J}\right)$, and with (20)

[TABLE]

About $\widehat{T}_{2}$. Following the guidelines in the proof of Lemma 1-$(ii)$, $\widehat{T}_{2}=\widehat{T}_{21}+\widehat{T}_{22}$ where $\widehat{T}_{21}\coloneqq\boldsymbol{e}^{\mathsf{T}}\,\left(\widehat{\boldsymbol{\Phi}}^{\mathsf{T}}\widetilde{\boldsymbol{K}}\widehat{\boldsymbol{\Phi}}\right)^{-1}\left[B_{0},\widehat{B}_{1},\dots,\widehat{B}_{J}\right]^{\mathsf{T}}$ and $\widehat{T}_{22}\coloneqq\boldsymbol{e}^{\mathsf{T}}\,\left(\widehat{\boldsymbol{\Phi}}^{\mathsf{T}}\widetilde{\boldsymbol{K}}\widehat{\boldsymbol{\Phi}}\right)^{-1}$ $\left[C_{0},\widehat{C}_{1},\dots,\widehat{C}_{J}\right]^{\mathsf{T}}$ with $\widehat{B}_{1},\ldots,\widehat{B}_{J}$ (resp. $\widehat{C}_{1},\ldots,\widehat{C}_{J}$) defined as $B_{1},\ldots,B_{J}$ (resp. $C_{1},\ldots,C_{J}$) but $\phi_{1},\ldots,\phi_{J}$ are replaced with $\widehat{\phi}_{1},\ldots,\widehat{\phi}_{J}$. One has $B_{0}=O_{P}(h^{2})$, $\widehat{B}_{j}=O_{P}(h^{3})$ for $j=1,\ldots,J$, and with (H(H1)), $C_{0}=O_{P}(h^{3})$ and $\widehat{B}_{j}=O_{P}(h^{4})$ for $j=1,\ldots,J$. Consequently, $\widehat{T}_{22}$ is negligible with respect to $\widehat{T}_{21}$. Similarly to $\widehat{T}_{1}$, $\widehat{T}_{21}=\widehat{\mu}\left(1-h\,\widehat{\boldsymbol{\delta}}^{\mathsf{T}}\widehat{\boldsymbol{\Delta}}^{-1}\boldsymbol{1}\right)O_{P}\left(h^{2}\right)$ so that

[TABLE]

About $\widehat{T}_{3}$. By the definition of $\widehat{T}_{3}$ and since $\max_{j\leq J}\|\phi_{j}-\widehat{\phi}_{j}\|=O_{P}\left(a_{J}^{-1}\,n^{-1/2}\right)$,

[TABLE]

Since $\widehat{\mu}\left(\delta_{0}-\widehat{\boldsymbol{\delta}}^{\mathsf{T}}\widehat{\boldsymbol{\Delta}}^{-1}\widehat{\boldsymbol{\delta}}\right)=1$, we have that

[TABLE]

Now it is enough to combine the decomposition of $\operatorname{E}_{\boldsymbol{X}}\left\{\widehat{\widehat{m}}(x)\right\}$ with (21), (22), and (23) to get the claimed conditional bias.

Let us now focus on the conditional variance. Similarly to Lemma 3-$(i)$,

[TABLE]

where $\widehat{\widetilde{\boldsymbol{\delta}}}$ (resp. $\widehat{\widetilde{\boldsymbol{\Delta}}}$) is defined as $\widetilde{\boldsymbol{\delta}}$ (resp. $\widetilde{\boldsymbol{\Delta}}$) but $\phi_{1},\ldots,\phi_{J}$ are replaced with $\widehat{\phi}_{1},\ldots,\widehat{\phi}_{J}$. From the definition of $\widehat{\widetilde{\boldsymbol{\delta}}}$ and $\widehat{\widetilde{\boldsymbol{\Delta}}}$, it is easy to state that $\widehat{\widetilde{\boldsymbol{\delta}}}=\widetilde{\boldsymbol{\delta}}\left(1+U_{n}\right)$ and $\widehat{\widetilde{\boldsymbol{\Delta}}}=\widetilde{\boldsymbol{\Delta}}\left(1+U_{n}\right)$ which results in $\operatorname{Var}_{\boldsymbol{X}}\left\{\widehat{\widehat{m}}(x)\right\}=\operatorname{Var}_{\boldsymbol{X}}\left\{\widehat{m}(x)\right\}\left(1+U_{n}\right)+$ $\{\sigma^{2}(x)\,+\,o(1)\}\,\widehat{\mu}^{2}\,\left(-\boldsymbol{\delta}^{\mathsf{T}}\boldsymbol{\Delta}^{-1}\widetilde{\boldsymbol{\delta}}-\widetilde{\boldsymbol{\delta}}^{\mathsf{T}}\boldsymbol{\Delta}^{-1}\boldsymbol{\delta}+\boldsymbol{\delta}^{\mathsf{T}}\boldsymbol{\Delta}^{-1}\widetilde{\boldsymbol{\Delta}}\boldsymbol{\Delta}^{-1}\boldsymbol{\delta}\right)U_{n}$. Let us now use Lemma 3-$(ii)$ with the fact that $\widehat{\mu}=O_{P}(\lambda_{J})$. This leads to $\operatorname{Var}_{\boldsymbol{X}}\left\{\widehat{\widehat{m}}(x)\right\}=O_{P}\left(\{n\,\pi_{x}(h)\}^{-1}\right)+O_{P}\left(a_{J}^{-1}\,n^{-1/2}\,\lambda_{J}\{n\,\pi_{x}(h)\}^{-1}\right)$. Since $a_{J}^{-1}\,n^{-1/2}=o(1)$, the claimed conditional variance holds. ∎

Appendix B Technical results

Lemma 5**.**

Let $p\geq 0$ and $q>0$. As soon as (H(H3)) and (H(H9)) are fulfilled,

[TABLE]

where $b_{x,p,q}\coloneqq K^{q}(1)-\int_{0}^{1}\left\{u^{p}K^{q}(u)\right\}^{\prime}\tau_{x}(u)\,\mathrm{d}\,u$.

A useful by-product of this lemma is the following corollary.

Corollary 1**.**

Under (H(H3)) and (H(H9)), we have that $\operatorname{E}K_{1}^{q}\,=\,b_{x,0,q}\,\pi_{x}(h)\,\left\{1+o(1)\right\}$ for any $q>0$.

Proof of Lemma 5..

In order to shorten the notations, set $Z_{h}=\|X_{1}-x\|/h$. From the differentiability of $K$ it comes that $t^{p}K^{q}(t)=K^{q}(1)-\int_{t}^{1}\left\{u^{p}K^{q}(u)\right\}^{\prime}\,\mathrm{d}\,u$ and

[TABLE]

Finally, $\int_{0}^{1}t^{p}K^{q}(t)\,\mathrm{d}\,P^{Z_{h}}(t)=\pi_{x}(h)\left\{K^{q}(1)-\int_{0}^{1}\left\{u^{p}K^{q}(u)\right\}^{\prime}\,\left\{\pi_{x}(hu)/\pi_{x}(h)\right\}\,\mathrm{d}\,u\right\}$ and thanks to (H(H9)) the claimed result holds. ∎

Lemma 6**.**

Under (H(H3))–(H(H9)) one has

1. (i)

$\operatorname{E}\left\{K_{1}^{a}\,\gamma_{j_{1},\ldots,j_{K}}^{p_{1},\ldots,p_{K}}\left(\|X_{1}-x\|^{p_{+}}\right)\right\}\,=\,b_{x,p_{+},a}\,{\gamma_{j_{1},\ldots,j_{K}}^{p_{1},\ldots,p_{K}}}^{{}^{\prime}}(0)\,h^{p_{+}}\,\pi_{x}(h)\,\left\{1+o(1)\right\}$* for any $a\geq 0$, where $p_{+}=p_{1}+\cdots+p_{K}$,* 2. (ii)

$\displaystyle\sum_{j\geq 1}\left\{{\gamma_{j}^{1}}^{{}^{\prime}}(0)\right\}^{2}\leq 1$, $\displaystyle\sum_{j_{1},j_{2}\geq 1}\left\{{\gamma_{j_{1},j_{2}}^{1,1}}^{{}^{\prime}}(0)\right\}^{2}\leq 1$, $\displaystyle\sum_{j_{1},j_{2},j_{3}\geq 1}\left\{{\gamma_{j_{1},j_{2},j_{3}}^{1,1,1}}^{{}^{\prime}}(0)\right\}^{2}\leq 1$, $\displaystyle\sum_{j_{1},\ldots,j_{4}\geq 1}\left\{{\gamma_{j_{1},\ldots,j_{4}}^{1,\ldots,1}}^{{}^{\prime}}(0)\right\}^{2}\leq 1$, $\displaystyle\sum_{j\geq 1}{\gamma_{j}^{2}}^{{}^{\prime}}(0)\leq 1$, and $\displaystyle\sum_{j_{1},j_{2}\geq 1}{\gamma_{j_{1},j_{2}}^{2,2}}^{{}^{\prime}}(0)\leq 1$, 3. (iii)

the $J\times J$ matrix $\boldsymbol{\Gamma}$ is positive semi-definite.

An interesting by-product of this lemma indicates that the Frobrenius norm of $\boldsymbol{\Gamma}$ is bounded from above by 1, and so is its largest eigenvalue.

Proof of Lemma 6-$(i)$..

We can write

[TABLE]

where $\eta(\epsilon)={{\gamma_{j_{1},\ldots,j_{K}}^{p_{1},\ldots,p_{K}}}^{{}^{\prime}}(\epsilon)}/{{\gamma_{j_{1},\ldots,j_{K}}^{p_{1},\ldots,p_{K}}}^{{}^{\prime}}(0)}-1$ with $0<\epsilon<h^{p_{+}}t^{p_{+}}$. Because $\gamma_{j_{1},\ldots,j_{K}}^{p_{1},\ldots,p_{K}}$ is continuously differentiable, $\sup_{\epsilon}|\eta(\epsilon)|=o(1)$, which implies

[TABLE]

Use of Lemma 5 gives in the claimed result. ∎

Proof of Lemma 6-$(ii)$..

Let us first remark that

[TABLE]

by using the definition of $\gamma_{j}^{1}(t)$. Because $\sum_{j\geq 1}\langle\phi_{j},X_{1}-x\rangle^{2}=\|X_{1}-x\|^{2}$, $\sum_{j\geq 1}\left\{\gamma_{j}^{1}(t)\right\}^{2}\leq\operatorname{E}\left(\|X_{1}-x\|^{2}\,|\,\|X_{1}-x\|=t\right)=t^{2}$. Now, by the definition of the derivative,

[TABLE]

So, for any $d\geq 1$, $\left|S_{d}(t)-S_{\infty}(t)\right|\leq\operatorname{E}\left(\zeta_{d}\,|\,\|X_{1}-x\|=t\right)$ where we write $S_{d}(t)\coloneqq\sum_{j=1}^{d}\left\{\gamma_{j}^{1}(t)\right\}^{2}$ and $\zeta_{d}\coloneqq\sum_{j>d}\,\langle\phi_{j},X_{1}-x\rangle^{2}$. It is clear that $\left\{\zeta_{d}\right\}_{d}$ is a non-increasing sequence of random variables that converges almost surely to zero with $d\to\infty$. The monotone convergence theorem implies that for any $t$, the sequence $\operatorname{E}\left(\zeta_{d}\,|\,\|X_{1}-x\|=t\right)$ is also a non-increasing sequence of random variables converging almost surely to zero with $d$. Consequently, $\sum_{j\geq 1}t^{-2}\left\{\gamma_{j}^{1}(t)\right\}^{2}$ converges uniformly on $(0,h)$. By remarking that for any $j$, $\lim_{t\rightarrow 0}t^{-2}\left\{\gamma_{j}^{1}(t)\right\}^{2}=\left\{{\gamma_{j}^{1}}^{{}^{\prime}}(0)\right\}^{2}$, one can exchange the limit with the sum in (24)

[TABLE]

The expression $\sum_{j>d}\langle\phi_{j},X_{1}-x\rangle^{2}=\|\mathcal{P}_{\mathcal{S}_{d}^{\perp}}(X_{1}-x)\|^{2}$ converges almost surely to 0 when $d$ tends to infinity, where $\operatorname{E}\left(\|\mathcal{P}_{\mathcal{S}_{d}^{\perp}}(X_{1}-x)\|^{2}\,|\,\|X_{1}-x\|=t\right)\leq t^{2}$. Thanks to [8] (see Corollary 2 - INT IV.37), one can put the infinite summation into the expectation

[TABLE]

which corresponds to the claimed first assertion. The other ones can be obtained by using similar arguments. ∎

Proof of Lemma 6-$(iii)$..

For any $J$-dimensional vector $\boldsymbol{u}$ and $t\geq 0$,

[TABLE]

∎

Lemma 7**.**

Under (H(H3))–(H(H9)), one has

[TABLE]

where $\boldsymbol{\Lambda}$ is the $J\times J$ matrix such that $\left[\boldsymbol{\Lambda}\right]_{j,k}\coloneqq\sqrt{{\gamma_{j,k}^{2,2}}^{{}^{\prime}}(0)}$.

Proof of Lemma 7..

From the definition of $\Delta_{jk}$,

[TABLE]

Use of Lemma 6 and Corollary 1 results in

[TABLE]

One can also derive the asymptotic behavior of $\operatorname{Var}\left(\Delta_{jk}\right)$

[TABLE]

where the last equality comes from Lemma 6 and Corollary 1. Equations (25) and (26) result in

[TABLE]

which is the element-wise version of the claimed result. ∎

Lemma 8**.**

As soon as (H(H3))–(H(H9)) are fulfilled,

[TABLE]

where $\boldsymbol{\theta}$ is the $J$-dimensional vector such that $\left[\boldsymbol{\theta}\right]_{j}\coloneqq\sqrt{{\gamma_{j}^{2}}^{{}^{\prime}}(0)}$.

Proof of Lemma 8..

This proof is shortened since it follows the same lines as the previous one. From the definition of $\delta_{j}$, $\operatorname{E}\delta_{j}\,=\,(\operatorname{E}K_{1})^{-1}\,\operatorname{E}\left\{K_{1}\,\gamma_{j}^{1}\left(\|X-x\|\right)\right\}$ and use of Lemma 6 and Corollary 1 results in

[TABLE]

Let us now focus on the variance of the $\delta_{j}$.

[TABLE]

the last equality resulting from Lemma 6 and Corollary 1. By combining this last equation with (27),

[TABLE]

which is the element-wise version of the claimed result. ∎

Lemma 9**.**

Under (H(H3))–(H(H9)), one has

1. (i)

$\displaystyle A_{0}\,=\,b_{x,0,1}^{-1}\,b_{x,1,1}\,\alpha_{0,x,n}^{bias}\,h\,\left\{1+o(1)\right\}\,+\,O_{P}\left(\|\mathcal{P}_{\mathcal{S}_{J}^{\perp}}m_{x}^{\prime}\|\,h\,\{n\,\pi_{x}(h)\}^{-1/2}\right)$*, *

where the sequence $\alpha_{0,x,n}^{bias}$ is upper bounded by $\|\mathcal{P}_{\mathcal{S}_{J}^{\perp}}m_{x}^{\prime}\|$, 2. (ii)

$\displaystyle A_{j}\,=\,b_{x,0,1}^{-1}\,b_{x,2,1}\,\alpha_{j,x,n}^{bias}\,h^{2}\,\left\{1+o(1)\right\}\,+\,O_{P}\left(h^{2}\sqrt{\alpha_{j,x,n}^{var}}\{n\,\pi_{x}(h)\}^{-1/2}\right)$* for $j=1,\ldots,J$, where $\sum_{j=1}^{J}\left\{\alpha_{j,x,n}^{bias}\right\}^{2}\leq\|\mathcal{P}_{\mathcal{S}_{J}^{\perp}}m_{x}^{\prime}\|^{2}$ and $\sum_{j=1}^{J}\alpha_{j,x,n}^{var}\leq\|\mathcal{P}_{\mathcal{S}_{J}^{\perp}}m_{x}^{\prime}\|^{2}$.*

Proof of Lemma 9-$(i)$..

We can write

[TABLE]

where $\alpha_{0,x,n}^{bias}=\sum_{j>J}\langle\phi_{j},m_{x}^{\prime}\rangle\,{\gamma_{j}^{1}}^{{}^{\prime}}(0)$, the last equality coming from Lemma 6-$(i)$. Moreover, the Cauchy-Schwartz inequality and Lemma 6-$(ii)$ imply that $|\alpha_{0,x,n}^{bias}|\leq\|\mathcal{P}_{\mathcal{S}_{J}^{\perp}}m_{x}^{\prime}\|$. In the same way, one has

[TABLE]

where $\alpha_{0,x,n}^{var}=\sum_{j>J}\sum_{k>J}\langle\phi_{j},m_{x}^{\prime}\rangle\langle\phi_{k},m_{x}^{\prime}\rangle\,{\gamma_{j,k}^{1,1}}^{{}^{\prime}}(0)$, the last equality using again Lemma 6-$(i)$. Moreover,

[TABLE]

which results in $\alpha_{0,x,n}^{var}\leq\|\mathcal{P}_{\mathcal{S}_{J}^{\perp}}m_{x}^{\prime}\|^{2}$ (use again Lemma 6-$(ii)$). Then $\operatorname{Var}(A_{0})\,=\,O\left(\|\mathcal{P}_{\mathcal{S}_{J}^{\perp}}m_{x}^{\prime}\|^{2}\,h^{2}\{n\,\pi_{x}(h)\}^{-1}\right)$ and the claimed result holds. ∎

Proof of Lemma 9-$(ii)$..

We have

[TABLE]

where $\alpha_{j,x,n}^{bias}=\sum_{k>J}\langle\phi_{k},m_{x}^{\prime}\rangle\,{\gamma_{j,k}^{1,1}}^{{}^{\prime}}(0)$, the last equality coming from the use of Lemma 6-$(i)$. In addition,

[TABLE]

To derive the asymptotic behavior of the variance of $A_{j}$, we follow similar arguments

[TABLE]

where $\alpha_{j,x,n}^{var}=\sum_{k>J}\sum_{\ell>J}\langle\phi_{k},m_{x}^{\prime}\rangle\langle\phi_{\ell},m_{x}^{\prime}\rangle\,{\gamma_{x,j,k,\ell}^{2,1,1}}^{{}^{\prime}}(0)$. Now, by involving arguments similar to those used in Lemma 6-$(ii)$ one is able to show $\sum_{j=1}^{J}\alpha_{j,x,n}^{var}\leq\|\mathcal{P}_{\mathcal{S}_{J}^{\perp}}m_{x}^{\prime}\|^{2}$. Just combine (28) and (29) to get the claimed result. ∎

Lemma 10**.**

For any $\boldsymbol{u},\boldsymbol{v}\in\mathbb{R}^{J}$, $\left|\boldsymbol{u}^{\mathsf{T}}\,\boldsymbol{\Gamma}^{-1}\,\boldsymbol{v}\right|\,\leq\,\lambda_{J}^{-1}\,\|\boldsymbol{u}\|_{2}\,\|\boldsymbol{v}\|_{2}$, where $\lambda_{J}$ is the smallest eigenvalue of $\boldsymbol{\Gamma}$.

Proof.

Let us remark that $\lambda_{J}>0$ according to (H(H5)). This result involves the Cauchy-Schwartz inequality and the Rayleigh quotient of the inverse of the $J\times J$ matrix $\boldsymbol{\Gamma}$:

[TABLE]

where $R_{1}$ and $R_{2}$ are the Rayleigh quotients of $\boldsymbol{\Gamma}^{-1}$. Let $\lambda_{J}$ be the smallest eigenvalue of $\boldsymbol{\Gamma}$. Then $\lambda_{J}^{-1}$ is the greatest eigenvalue of $\boldsymbol{\Gamma}^{-1}$ with $R_{1}\leq\lambda_{J}^{-1}$ and $R_{2}\leq\lambda_{J}^{-1}$, and the claimed result holds. ∎

Lemma 11**.**

Under (H(H3))–(H(H9)), one has

1. (i)

$\displaystyle B_{0}\,=\,b_{x,0,1}^{-1}b_{x,2,1}\beta_{0,x}^{bias}\,h^{2}\left\{1+o(1)\right\}\,+\,O_{P}\left(h^{2}\,\{n\,\pi_{x}(h)\}^{-1/2}\right)$* with $\beta_{0,x}^{bias}=O(1)$,* 2. (ii)

$\displaystyle B_{j}\,=\,b_{x,0,1}^{-1}\,b_{x,3,1}\,\beta_{j,x}^{bias}\,h^{3}\,\left\{1+o(1)\right\}\,+\,\sqrt{\beta_{j,x}^{var}}\,O_{P}\left(h^{3}\,\{n\,\pi_{x}(h)\}^{-1/2}\right)$* for all $j=1,\ldots,J$, where $\sum_{j=1}^{J}\left\{\beta_{j,x}^{bias}\right\}^{2}=O(1)$ and $\sum_{j=1}^{J}\beta_{j,x}^{var}=O(1)$.*

Proof of Lemma 11-$(i)$..

A standard expansion of the linear operator $m_{x}^{\prime\prime}$ results in

[TABLE]

where $\beta_{0,x}^{bias}\coloneqq\sum_{j\geq 1}\sum_{k\geq 1}\langle m_{x}^{\prime\prime}\phi_{j},\phi_{k}\rangle\,{\gamma_{j,k}^{1,1}}^{{}^{\prime}}(0)$, the last equality by Lemma 6-$(i)$. Moreover, the Cauchy-Schwartz inequality combined with (H(H1)) and Lemma 6-$(ii)$ imply that $\beta_{0,x}^{bias}=O(1)$. In the same way, one has

[TABLE]

where $\beta_{0,x}^{var}\coloneqq\sum_{j_{1},\ldots,j_{4}\geq 1}\langle m_{x}^{\prime\prime}\phi_{j_{1}},\phi_{j_{2}}\rangle\langle m_{x}^{\prime\prime}\phi_{j_{3}},\phi_{j_{4}}\rangle\,{\gamma_{j_{1},\ldots,j_{4}}^{1,1,1,1}}^{{}^{\prime}}(0)$, the last equality using again Lemma 6-$(i)$. Moreover,

[TABLE]

which gives that $\beta_{0,x}^{var}$ is a finite quantity (use again (H(H1)) and Lemma 6-$(ii)$). Then $\operatorname{Var}(B_{0})\,=\,O\left(h^{4}\{n\,\pi_{x}(h)\}^{-1}\right)$ and the claimed result holds. ∎

Proof of Lemma 11-$(ii)$..

Similarly,

[TABLE]

where $\beta_{j,x}^{bias}\coloneqq\sum_{k\geq 1}\sum_{\ell\geq 1}\langle m_{x}^{\prime\prime}\phi_{k},\phi_{\ell}\rangle\,{\gamma_{j,k,\ell}^{1,1,1}}^{{}^{\prime}}(0)$. Moreover,

[TABLE]

where $\left\|\cdot\right\|_{HS}$ denotes the Hilbert-Schmidt norm of an operator. So, $\sum_{j=1}^{J}\left\{\beta_{j,x}^{bias}\right\}^{2}=O(1)$. Let us now focus on the variance of $B_{j}$

[TABLE]

where $\beta_{j,x}^{var}\coloneqq\sum_{k,\ell,p,q\geq 1}\langle m_{x}^{\prime\prime}\phi_{k},\phi_{\ell}\rangle\langle m_{x}^{\prime\prime}\phi_{p},\phi_{q}\rangle\,{\gamma_{j,k,\ell,p,q}^{2,1,1,1,1}}^{{}^{\prime}}(0)$, the last equality using again Lemma 6-$(i)$. Using similar arguments as those involved to show the second assertion of Lemma 6-$(ii)$,

[TABLE]

Consequently,

[TABLE]

where, of course, $\gamma_{k,\ell,p,q}^{2,2,2,2}\left(t^{4/3}\right)$ is the same as

[TABLE]

∎

Appendix C Practical aspects

C.1 fllr R package

An efficient, fully documented R implementation of all the considered estimating procedures, including the automated selection of all their parameters, is freely available as a part of R package fllr. The package can be downloaded from https://bitbucket.org/StanislavNagy/fllr. Using the procedures from fllr and the source codes accompanying the present manuscript, all the results from the main paper, and the simulation studies presented in the supplementary material, can be replicated in full.

C.2 Bandwidth selection

Functional derivative: bootstrap bandwidth selection. The bootstrap procedure introduces additional randomness into our local linear estimation method. A natural question raises: is the proposed bandwidth selection stable? In other words, if the functional derivatives are estimated several times on the same dataset, are the resulting bandwidths similar? A related issue concerns the number $B$ of bootstrap repetitions set by the user: is the procedure sensitive to this parameter, and how to choose it? To address both these concerns, one dataset is simulated according to (M1). For different $B$ and learning sample sizes $n$, our estimating algorithm is launched 100 times. Table 5 displays the mean and standard deviation (in brackets) of $h_{deriv}$ and $ORMSEP_{deriv}$ in (a) the worst case, and (b) the most favorable case in the simulation study from the main document.

The variability of the selected bandwidth is smaller in the most favorable case as expected. The bootstrap bandwidth selection remains stable in both cases despite the additional randomness introduced. In both situations the variability of $ORMSEP_{deriv}$ is close to 0. The default value ($B=100$) used in our procedure seems to be large enough to ensure stability as well as accuracy for predictions.

Raw results. Table 6 gives the mean and the standard deviation (in brackets) of the computed bandwidths over 100 runs in 27 situations (9 learning sample sizes $\times$ 3 noise-to-signal ratios) and corresponds to Figure 5 in the main document.

C.3 Robustness to the model complexity

We now provide several additional tables of results to assess the robustness of our estimating procedure with respect to the complexity of the simulated model (see Section 5.1 in the main paper). Similarly to model (M2), $Y\coloneqq\sum_{j=1}^{J}\exp(-U_{j}^{2})+\varepsilon$ with $X\coloneqq\sum_{j=1}^{J}U_{j}\,\phi_{j}\,+\,\eta$ and where $\eta\coloneqq\sum_{j=J+1}^{2J}V_{j}\,\phi_{j}(t)$ is a structural perturbation acting on the functional predictor $X$. Variables $U_{j}$ (resp. $V_{j}$) are iid uniform on $[-1,\,1]$ (resp. $[-b,\,b]$). Again, the structural perturbation is controlled by the ratio $\rho\coloneqq b^{2}/(1+b^{2})$. The noise-to-signal ratio ($nsr$) of the regression model is set to 0.05 and 0.4. Tables 7–16 display, for $J$ respectively set to 2, 3 and 4 ($J=4$ only when $nsr=0.4$), i) the number of times, out of 100, that the dimension is correctly selected (cf. Table 1 in the main document), and ii) the corresponding $ORMSEP_{reg}$ and $ORMSEP_{deriv}$ averaged over 100 runs with standard deviation in brackets (cf. Table 2 in the main document).

C.4 Complement to the study from Section 5.2

Tables 17 and 18 supplement results given in Tables 3 and 4 for two additional perturbation levels and $500$ learning functions, and Tables 19–22 provide analogous results for samples of $100$ learning functions. The conclusions from the main paper all hold true. Especially in the situation when the learning sample size is low, the method (MY) seems to be numerically unstable. Remarkably, the local linear estimation method does not appear to suffer from such drawbacks.

C.4.1 Results for 500 learning functions

C.4.2 Results for 100 learning functions

C.5 Complement to the data analysis from Section 5.3

As explained in Section 5.3, restricting the growth velocity profiles from ages 1–10 to 6–10 does not degrade the quality of estimation. Figure 8 displays the observed responses versus their estimates when considering the whole growth velocity profile (1–10), or the restricted one (6–10). To quantify the performance of the estimating procedure, the empirical (Pearson’s) correlation coefficient between the observations and their estimates is computed in each situation. When the regression model involves the whole trajectory of the growth velocity, the correlation equals 0.823; in the other case where estimates are based on ages 6–10, one gets 0.819. The accuracy of the estimating procedures are almost the same, which confirms that the behavior of the growth velocity profile under 6 years of age does not influence the adult height at 18.

Appendix D Theoretical complement on (H(H5))

In this section, we investigate hypothesis (H(H5)) that requires regularity of the family of functions

[TABLE]

The next lemma provides a general condition on the functional predictor $X$ in order to fulfill (H(H5)).

Lemma 12**.**

Suppose that for some $x\in H$ the random vector

[TABLE]

is absolutely continuous with a density in $\mathbb{R}^{J+1}$ that is positive at the origin, continuous at the origin in its first $J$ coordinates, and continuous at the origin from the right in its last coordinate. Then condition (H(H5)) is satisfied for all functions in $H$.

Proof.

For $x\in H$ fixed and $j=1,\dots,J$ denote $Z_{j}=\langle\phi_{j},X_{1}-x\rangle$. We want to establish that the matrix

[TABLE]

is positive definite. That is equivalent with the fact that for any $\boldsymbol{u}=\left(u_{1},\dots,u_{J}\right)^{\mathsf{T}}\in\mathbb{R}^{J}$, $\left\|\boldsymbol{u}\right\|=1$

[TABLE]

For any univariate random variable $Z$ with variance and $z\in\mathbb{R}$ we know that

[TABLE]

with equality if and only if $z=\operatorname{E}Z$. Use the conditional version of this inequality to obtain

[TABLE]

Therefore, it suffices to show that the conditional variance of no projection of the vector $\left(Z_{1},\dots,Z_{J}\right)^{\mathsf{T}}$ into a line spanned by a unit vector is of order $o(t)$ with $t\to 0$.

We assume that the random vector $\boldsymbol{Z}=\left(Z_{1},\dots,Z_{J},\sqrt{\sum_{i=J+1}^{\infty}Z_{i}^{2}}\right)^{\mathsf{T}}$ is absolutely continuous in $\mathbb{R}^{J+1}$. For an independent Rademacher random variable $R$, i.e. $P\left(R=1\right)=P\left(R=-1\right)=1/2$, define $\widetilde{\boldsymbol{Z}}=\left(Z_{1},\dots,Z_{J},R\sqrt{\sum_{i=J+1}^{\infty}Z_{i}^{2}}\right)^{\mathsf{T}}$. This random vector is absolutely continuous, with density $f_{J}$ positive and continuous at the origin. It differs from the original random vector $\boldsymbol{Z}$ only in its last coordinate, and $\boldsymbol{Z}^{\mathsf{T}}\boldsymbol{Z}$ has the same distribution as $\widetilde{\boldsymbol{Z}}^{\mathsf{T}}\widetilde{\boldsymbol{Z}}$. The conditional density of $\widetilde{\boldsymbol{Z}}$ given $\sum_{i=1}^{\infty}Z_{i}^{2}=\widetilde{\boldsymbol{Z}}^{\mathsf{T}}\widetilde{\boldsymbol{Z}}=t$ takes the form

[TABLE]

where $1\left[t=\boldsymbol{z}^{\mathsf{T}}\boldsymbol{z}\right]$ is $1$ if $t=\boldsymbol{z}^{\mathsf{T}}\boldsymbol{z}$, [math] otherwise. The integral in (32), and in analogous expressions below, is taken with respect to the Hausdorff measure on an appropriate sphere in $\mathbb{R}^{J+1}$. By our assumptions, $f_{J}$ is positive and continuous in the neighborhood of the origin. Then, for $t$ small enough, $c_{J,t}=\inf_{\{\boldsymbol{z}^{\mathsf{T}}\boldsymbol{z}=t\}}f_{J}(\boldsymbol{z})$ must be positive. Using (30) again, we can therefore write

[TABLE]

Here, $g_{t}(\boldsymbol{z})=\left(\int_{\{\boldsymbol{v}^{\mathsf{T}}\boldsymbol{v}=t\}}1\,\mathrm{d}\,\boldsymbol{v}\right)^{-1}$ is the reciprocal of the Hausdorff measure of a sphere, and $\boldsymbol{U}=\left(U_{1},\dots,U_{J+1}\right)^{\mathsf{T}}$ is a random vector distributed uniformly on the unit sphere in $\mathbb{R}^{J+1}$. The second inequality is from (30). The first equality on the last line in the formula above follows from the spherical symmetry of the vector $\boldsymbol{U}$ — any projection of $\boldsymbol{U}$ onto a line has the same distribution as $U_{1}$. The final equality follows from $\operatorname{Var}\boldsymbol{U}=\boldsymbol{I}/(J+1)$, for $\boldsymbol{I}$ the $(J+1)\times(J+1)$ identity matrix.

From the continuity of $f_{J}$ around the origin it follows that

[TABLE]

Therefore, we obtain

[TABLE]

and also the desired

[TABLE]

Note that the uniformity in $\boldsymbol{u}$ is assured by the spherical symmetry of $\boldsymbol{U}$ utilized above.

Finally, to see that, for instance, function ${\gamma_{j,k}^{1,1}}$ is continuously differentiable at the origin, note that from (32) we get

[TABLE]

for $z_{j}$ and $z_{k}$ the $j$-th and $k$-th elements of the vector $\boldsymbol{z}\in\mathbb{R}^{J+1}$, respectively. The assertion follows from the continuity of $f_{J}$, and both integrands above in the formula for $\gamma_{j,k}^{1,1}(t)$, around the origin and the fundamental theorem of calculus. ∎

For a non-degenerate Gaussian process we know that each $Z_{i}$ has a univariate, non-degenerate normal distribution. The distributions of $Z_{i}$ are generally correlated, with

[TABLE]

for $B$ the covariance operator of the process $X$. But, for the special case of the $\left\{\phi_{i}\right\}$ being the eigenbasis of $B$, $\left\{Z_{i}\right\}$ is a sequence of independent random variables with non-degenerate normal distributions. Therefore, in this case the elements of the random vector $\widetilde{\boldsymbol{Z}}$ are independent, and absolutely continuous random variables. Absolute continuity of the last element $\sqrt{\sum_{i=J+1}^{\infty}\langle\phi_{i},X_{1}-x\rangle^{2}}$ follows from the orthonormality of the basis $\left\{\phi_{i}\right\}$, independence and absolute continuity of the terms of the sequence $\left\{Z_{i}\right\}=\left\{\langle\phi_{i},X_{1}-x\rangle\right\}$, and absolute continuity of $\left\|X_{1}-x\right\|$. All the marginal densities of vector $\widetilde{\boldsymbol{Z}}$ are positive and continuous at the origin. Therefore, for any non-degenerate Gaussian process with $\left\{\phi_{i}\right\}$ its collection of eigenfunctions, the assumptions of Lemma 12 are satisfied, and $X$ satisfies (H(H5)).