Nonparametric tests for semiparametric regression models

Semiparametric regression models have received considerable attention over the last decades, because of their flexibility and their good finite sample performances. Here we propose an innovative nonparametric test for the linear part of the models, based on random sign-flipping of an appropriate transformation of the residuals, that exploits a spectral decomposition of the residualizing matrix associated with the nonparametric part of the model. The test can be applied to a vast class of extensively used semiparametric regression models with roughness penalties, with nonparametric components defined over one-dimensional, as well as over multi-dimensional domains, including, for instance, models based on univariate or multivariate splines. We prove the good asymptotic properties of the proposed test. Moreover, by means of extensive simulation studies, we show the superiority of the proposed test with respect to current parametric alternatives, demonstrating its excellent control of the Type I error, accompanied by a good power, even in challenging data scenarios, where instead current parametric alternatives fail.


Introduction
Semiparametric regression models have a long history in statistics (see, e.g., the textbooks Green and Silverman 1994;Bickel et al. 1998;Ruppert et al. 2003, and references therein).Because of their flexibility and versatility, they have been the object of an extensive and still very active literature.In this work, we propose an efficient (conditional) resampling-based test (Pesarin 2001;Hemerik and Goeman 2018b;Chung and Romano 2013) for the linear component in partially linear and semiparametric regression models with roughness penalties.The test can be applied to a vast class of extensively used models, with nonparametric components defined over one-dimensional, as well as over multi-dimensional domains, including manifold domains.This embraces, for instance, the highly popular semiparametric regression models based on splines (see, e.g., Heckman 1986;Yu and Ruppert 2002;Wand and Ormerod 2008;Wang 2019, and references therein), on thin-plate splines (see, e.g., Wood 2003), and on spherical splines (Wahba 1981), as well as semiparametric models based on recent smoothing techniques over two-dimensional (possibly irregularly shaped or curved) domains, such as soap film smoothing (Wood et al. 2008), bivariate-splines over triangulations (Lai and Schumaker 2007;Baramidze et al. 2006;Lai et al. 2009;Guillas and Lai 2010;Lai and Wang 2013;Wang et al. 2020), and Spatial Regression with Partial Differential Equation regularization (SR-PDE) (see, e.g., Sangalli et al. 2013;Azzimonti et al. 2015;Ettinger et al. 2016;Wilhelm et al. 2016;Sangalli 2021).
Various classical approaches are available to make inference in the context of semiparametric regressions, and different strategies have been proposed to cope with the bias induced by the roughness penalty.Some possibilities include undersmoothing approaches developed for nonparametric models [see, e.g., the review in Hall and Horowitz (2013)], Bayesian approaches (Wahba 1983;Nychka 1988;Marra and Wood 2012) and various corrections of Wald-type test statistics, such as the sandwich estimators in Gray (1994) and Yu and Ruppert (2002), and the Speckman's version in Speckman (1988) and Holland (2017).These approaches might nonetheless have poor performances in the finite sample scenario, due to the effects of the roughness penalty (see, e.g., Maas and Hox 2004;Freedman 2006).In particular, as also evidenced by the simulation studies reported in this work, such tests have a poor control of Type I error.
Here we propose an innovative test for the linear part of semiparametric regression models, based on conditional resampling of a transformation of the residuals.This test, unlike other proposals, allows to overcome the problem of dependence in the residuals that is particularly strong in semiparametric models.Some approaches proposed in the context of classical regression models, such as those in Huh and Jhun (2001) and Kherad-Pajouh and Renaud (2010), derive transformed residuals from spectral decomposition of the residualizing matrix that projects into the residual space.In the setting considered by these authors, the conditional distribution of the test statistic can be defined on the basis of permutations (see, e.g., Pesarin 2001;Chung and Romano 2013;Pauly et al. 2015;Winkler et al. 2014), rotations (Solari et al. 2014) or sign-flips (Hemerik et al. 2020) of such transformed residuals.These approaches are nonetheless not valid in this context, since the residualizing matrix is not idempotent in the case of penalized regression models.Because of this, the transformed residuals are not spherical (i.e., they are not homoscedastic and independent) and the standard permutation, rotation or sign-flip procedures become invalid in our context.To overcome this problem, we here study a conditional sign-flip procedure, named eigen sign-flip test, that preserves the finite sample covariance structure of the residuals, hence ensuring asymptotically exactness of the derived test.This idea has been explored in Ferraccioli et al. (2021), restricted to a specific case of SR-PDE model.The current work addresses instead the broad spectrum of highly popular semiparametric regression models mentioned above.Moreover, we study in detail the asymptotic properties of the test.In particular, we prove the asymptotic exactness of the test and derive similar results for interval hypothesis and confidence intervals.Some of the obtained results leverage on the asymptotic properties of the estimator of the nonparametric part of the model.Such properties in turn depend on conditions that are model-specific, since they depend, for instance, on the dimension and geometry of the domain over which the nonparametric term is defined, on the roughness term being considered, on the type of basis, etc.In the present work, we hence define assumptions that are general enough to cover a variety of semiparametric regression models, and refer the reader to other works for the appropriate specifications of such assumptions for the specific model being considered (e.g., to Claeskens et al. (2009) for univariate penalized splines estimators, to Holland (2017) for multivariate penalized splines estimators, to Xiao (2019) for general penalized splines and to Arnone et al. (2021) for SR-PDE.) The paper is organized as follows.In Sect. 2 we briefly review the semiparametric penalized regression framework, outlining the forms of the associated discrete estimators.In Sect. 3 we recall some classical parametric approaches for inference on the linear part of a semiparametric regression model and summarize the properties of the score test statistic in this context.In Sect. 4 we present the eigen-sign flip test and describe its theoretical and asymptotic properties.In Sect. 5 we compare our proposal to more classical parametric approaches in extensive simulation studies.In Sect.6 we present an application to the study of human development in Nigeria.Finally, some discussions and possible directions for future research are outlined in Sect.7.

Semiparametric regression
Let y i ∈ R be the value of the variable of interest observed in correspondence of covariates x i ∈ R q and of p i ∈ ⊆ R d , d ≥ 1.We consider the semiparametric model where β ∈ R q is the vector of regression parameters, f is a real-valued smooth function on , and i are i.i.d.random errors with E( i ) = 0 and E( 2 i ) = σ 2 .The interest is to estimate both the linear coefficients β and the nonparametric component f .However, the estimation of (β, f ) in model (1) via maximum likelihood is usually inappropriate or infeasible, due to the infinite-dimensionality of the nonparametric component f .To avoid this problem, some type of roughness penalty can be imposed, in order to reduce the space of possible solutions.In general, the resulting penalized likelihood estimators for β and f are the solution of the minimization problem argmin where P(•) is some type of roughness penalty.Depending on the assumptions on the domain ⊆ R d , on the dimension d, and on the required smoothness of the function f , various proposals for P( f ) have been considered in the literature, and different discretization procedures have been adopted to reduce the infinite-dimensional estimation problem (2) to a finite dimensional one.For instance, for d = 1 and an interval of the real line, model ( 1)-( 2) can involve the classical and extensively used O'Sullivan splines (O'Sullivan 1986;Heckman 1986;Yu and Ruppert 2002;Wand and Ormerod 2008), whose penalty is the integrated squared derivative of some order, and can, for instance, rely on B-spline bases.When is the real plane, it is possible to use thin-plate splines (see, e.g., Duchon 1977;Wahba 1990;Wood 2003), which involve as penalty the so-called thin-plate energy.Moreover, various recent techniques target twodimensional bounded planar domains ⊂ R 2 , including: soap-film smoothing (Wood et al. 2008) that considers a penalty involving the Laplacian of f ; bivariate-splines over triangulations (Lai and Schumaker 2007;Guillas and Lai 2010;Lai and Wang 2013), whose regularizing term may include high-order derivatives; SR-PDE (Sangalli et al. 2013;Azzimonti et al. 2015), where the regularizing term can involve general second-order partial differential equations, and the estimation problem is discretized via finite element bases (Sangalli et al. 2013;Azzimonti et al. 2015) or advanced spline bases (Wilhelm et al. 2016).Some of these techniques also permit the constructions of semiparametric models over spherical domains (Wahba 1981;Baramidze et al. 2006;Lai et al. 2009) and general surface domains (Ettinger et al. 2016;Wilhelm et al. 2016).

Discrete estimators
The estimation of model (1) usually involves the representation of the nonparametric component f through some type of basis expansion, depending on the penalization being considered.Let ∈ R n × R K be the matrix of the evaluations of the K basis functions ψ 1 , . . ., ψ K at the n data locations p 1 , . . ., p n , that is, Then, we write ( f (p 1 ), . . ., f (p n )) = γ for some vector of coefficients γ ∈ R K .Moreover, let P denote the K × K positive semidefinite matrix representing the discretization of the penalty P(•).Finally, set y = (y 1 , . . ., y n ) and denote by X ∈ R n ×R q the design matrix, whose i-th row is given by x i .The estimation problem (2) is therefore discretized as The solution to (3) is uniquely determined by the normal equations the explicit form of the estimators for β and γ is, respectively, or equivalently 3 Inference on Ǐn semiparametric regression, a natural question is whether the covariates X have an effect on the variable of interest.We are thus interested in the system of hypotheses A standard approach to verify (9) is to use a Wald-type test (see, e.g., Schervish 2012), based on the asymptotic distribution of β.The study of the asymptotic distribution of β, in semiparametric regression models, has been tackled by a number of works.See, for instance, Heckman (1986); Yu and Ruppert (2002); Li and Ruppert (2008); Holland (2017); Xiao (2019); Yu et al. (2019); Wang et al. (2020) for semiparametric models based on univariate and bivariate splines.
The parametric Wald-type test may nonetheless have poor performances in small sample scenarios, due to the overestimation of the variance of the test statistic, induced by the penalization.A number of corrections to Wald-type test have been proposed to avoid this issue, such as the sandwich estimators in Gray (1994) and Yu and Ruppert (2002) and the Speckman's version in Speckman (1988) and Holland (2017).Nonetheless, these approaches can only partially solve the problem, and may lead to a poor control of the Type I error, especially when a strong temporal/spatial structure in the covariates is present, as indicated by the simulations carried on in Sect. 5.
In the Sect. 4 we introduce an innovative nonparametric alternative for testing on β.Such proposal is based on the score statistic.For this reason, in the remainder of this section we review the properties of the score statistic in the context of penalized semiparametric regression.The proposed method does not rely on the estimation of the Fisher information matrix to define the null distribution, which is implicitly recovered by an appropriate nonparametric resampling procedure, as described in Sect. 4.

Properties of the score statistic in penalized semiparametric regression
We first study the distributional properties of the score statistic, which constitute the base of the nonparametric test defined in Sect. 4. Using the normal equation (4), we can define the classical score test statistic Since γ is unknown, we can use the plug-in γ .Substituting γ in expression (10), we define the test statistic T as with r = y − X β 0 .We make the following assumption: (A1) For n large enough, the matrix is positive definite.
Assumption (A1) is quite general; its specification depends on the basis considered.
In particular, this specification usually involves conditions on the nodes of the basis and their position with respect to the design points p 1 , . . ., p n .More specifically, it involves the type of basis, the rate at which the number of bases K grows with n, the minimum distance between the nodes, and the density of the design points inside the domain.For instance, in the case of univariate penalized splines estimators, (A1) follows from Assumptions 1-3 in Claeskens et al. (2009).In the case of multivariate penalized splines estimators, it follows from Assumptions 1-2 in Holland (2017).In the case of SR-PDE, it follows from Assumptions 3-5 in Arnone et al. (2021).
Here we consider the case of fixed designs, thus implicitly conditioning on the sample points and the covariates.Similar results can be obtained in the random design scenario, by introducing further assumptions on the distribution of the design points and covariates (e.g., that the covariates are realizations of continuous processes on ).
Under (A1), we can consider the Demmler and Reinsch (1975) decomposition where U is the matrix of eigenvectors, and ρ is the corresponding vector of eigenvalues {ρ k } K k=1 (see Eubank 1999, for details).Let us also denote A = ( ) −1/2 U .Note that this matrix is semi-orthogonal, i.e., A A = I K and A A = ( ) −1 .Following Demmler and Reinsch (1975), we can rewrite the matrix in (6) as Using this decomposition, we can now study the behavior of the bias of the test statistic T , in terms of the eigenvalues ρ k .Lemma 3.1 Assume (A1) and let X = A X and γ = A γ .Let also xi be the q-dimensional vectors corresponding to the rows of X , and γ i be the elements of the vector γ .Under the null hypothesis (9), the bias b λ of T is where the inequality is considered element-wise.
Proof Denote by the n-dimensional vector of i.i.d.residuals.Under the null hypothesis, we have since the term E( ) is zero by assumption.Using the decomposition in (3.1), it follows that Substituting ( 13) in b λ , we obtain Using the notation X = A X and γ = A γ , the bias can therefore be rewritten as where xi are the q-dimensional vectors corresponding to the rows of X , and γ i the elements of the vector γ .Equation ( 14) highlights that the bias is a sum of K contributions, weighted by the eigenvalues ρ k , and moderated by λ.Since the function x/(1 + x) < x, for x > 0, we can bound the bias as follows The expression ( 14) highlights how the bias depends on the chosen penalization through the eigenvalues ρ k .We finally make the following assumption.
(A2) The smoothing parameter λ = λ n is chosen so that λ K i=1 ρ i = o(1).Thanks to Lemma 3.1, assumption (A2) implies the asymptotic unbiasedness of score statistic T , since f is a continuous function on the bounded domain and the covariates are realizations of a continuous process on .This is a standard assumption when studying the asymptotic properties of semiparametric and nonparametric penalized regression models.Likewise for Assumption (A1), also Assumption (A2) needs to be specified depending on the penalty and basis considered.Indeed, Assumptions (A1)-(A2) are intentionally left quite general to embrace various semiparametric models; moreover, the precise rates of convergence are not of direct interest in this work.Theorem 1 in Claeskens et al. (2009) gives, for instance, the appropriate rates for λ in the case of univariate penalized splines estimators, Theorem 3 in Holland (2017) gives it for multivariate penalized spline estimators, while Lemma 3 in Arnone et al. (2021) gives it for SR-PDE estimators.
We can now state the main result for the asymptotic distribution of the test statistic T .
Theorem 3.2 Let ν = σ 2 X 2 X .Under the assumptions (A1)-(A2), the test statistic T in (11) is asymptotically normal under the null hypothesis (9), with Proof We know that where the notation [X ] i is used to indicate the i-th column of the q ×n matrix X .Under assumption (A2), it follows from (12) that the bias b λ is asymptotically zero.The expected value E(T ) is therefore asymptotically zero.For the variance, under the null hypothesis we have Substituting the expression of from equation ( 13) in the previous expression, we obtain Using the notation X = A X and completing the square in the second term, we hence get where xi are the q-dimensional vectors corresponding to the rows of X .Note that the first term does not depend on λ.As for the second term, since x 2 /(1 + x) 2 < x 2 , for x > 0, we have where the maximum is taken element-wise.Therefore, for n large enough (since the covariates are realizations of a continuous process on ), assumption (A2) implies that the second term in (15) vanishes faster than the first term.Concerning the first term in (15), with a similar argument it is easy to check that the matrix A A is idempotent with rank K .Thus, it admits the spectral decomposition A A = U diag(1, . . ., 1, 0, . . ., 0)U , with the first K non-null eigenvalues equal to 1.The term X (I n − A A )X is therefore the sum of n − K components with bounded variance, since the covariates are realizations of a continuous process on , thus the Feller condition is satisfied.It follows from the central limit theorem [see, e.g., Van der Vaart (2000)] that the test statistic T is also asymptotically normal.

Eigen sign-flip test for the linear component in penalized semiparametric regression models
In the classical linear regression case, under the standard assumption of i.i.d.random noise, the score statistics can also be viewed as a sum of n contributions that have asymptotically zero mean, under the null hypothesis H 0 (9).This information can be used to derive the null distribution of the test statistic, without the need of a direct estimation of the Fisher information.In the context of semiparametric regression, instead, a first naive attempt to derive the distribution of the test statistic can be made by random permutations (or sign-flips) of the contributions of the score (Winkler et al. 2014;Hemerik et al. 2020).This approach, attempted in Ferraccioli (2020) for a simple type of SR-PDE model (Sangalli et al. 2013), might nonetheless be not optimal in the semiparametric regression setting.The reason for this lies in the fact that naive permutation does not account for the correlation between residuals, nor for the bias of the estimates, which is inherent to semiparametric models.To solve this issue, always considering a special case of SR-PDE model, Ferraccioli et al. ( 2021) defines a new test statistic, that leverages on the spectral decomposition of the matrix , leading to the definition of the eigen sign-flip test.
We here defined the eigen sign-flip test on β for a general forms of penalized semiparametric regression models.We study the properties of the test, proving its asymptotic distribution.A thorough discussion on the nature of the proposed test in given in Sect.4.2.
Definition 1 (Eigen sign-flip test) Let us consider the singular value decomposition = V DV .Set = diag(π 1 , . . ., π n ), where π = (π 1 , . . ., π n ) is a random vector uniformly distributed in {−1, 1} n .Let us also define the n-dimensional vectors X = D 1/2 V X and r = D 1/2 V r = D 1/2 V (y−X β 0 ).The eigen sign-flip statistics is defined as Note that the observed statistic T = T I corresponds to the case where As standard in permutational approaches, the component-wise p-values are thus computed as the rank of T I with respect to a sample of M sign-flips π, divided by M (see, e.g., Pesarin 2001).

Asymptotic properties of the eigen sign-flip test
We now study the asymptotic properties of the test statistic T in Definition 1.We first show that the asymptotic distribution of the test statistic T is the same as T I .We then show that the eigen sign-flip test is asymptotically exact.
Theorem 4.1 Let ν = σ 2 X 2 X .Under the assumptions (A1)-( A2), for any given , the distribution of T is asymptotically normal, with Proof For the expected value, under the null hypothesis we have Following the same reasoning of the proof of Theorem 3.2, but with the quantity V V X in place of X , we can show that the expected value of T is asymptotically zero.
As for the variance, under the null hypothesis we have It follows from the central limit theorem (Van der Vaart 2000), that the test statistic T is also asymptotically normal.
Remark 1 Note that the bias in the mean of the test statistic is intrinsic in the regularization approach, and cannot be avoided in the finite sample scenario.Because of this bias, we are only able to reach asymptotically exact results.
Remark 2 Note also that the matrix is defined so that it commutes with D. This is necessary to ensure that the variance of the test statistic is invariant under the action of .
We now introduce some notation before establishing the main result, that constitutes the pivot point to prove the asymptotic control of the probability of Type I error.For the sake of simplicity of exposition, we consider the results for a single covariate case in the remainder of this section and in Sect.4.3.In Sect.4.4, we outline the procedure for the general multivariate case.Let α ∈ [0, 1).For any a ∈ R, let a be the smallest integer which is larger than or equal to a and let a be the largest integer which is at most a.We consider all the possible w = 2 n sign-flips 1 , . . ., w , where ) .This will prove the asymptotic control of the Type I error through Theorem 15.2.1 in Lehmann and Romano (2008) and Theorem 1 in Hemerik and Goeman (2018a).Let X be the diagonal n × n matrix with elements ( X1 , . . ., Xn ).The test statistic in Definition 1 can be rewritten as where 1 n is the n-dimensional unit vector.The test statistic T can hence be viewed as sum of n contributions, where each element of Xr is sign-flipped through .Similarly, the variance of T can be written as see also Theorem 4.1.
To evaluate the joint distribution of the test statistics T, let us now define as the 2 n × n matrix collecting all the w = 2 n vectors of sign-flip row-wise.Therefore, we can write T = Xr, and • T = Xr.The joint distribution of T is multivariate normal with variance Var(T) = σ 2 n −1 XD X and asymptotically zero mean.We now have to show that • T follows the same asymptotic multivariate normal distribution.First note that the transformation does not affect the expected value that remains asymptotically zero.Furthermore, for the variance we have Thanks to Theorem 15.2.1 in Lehmann and Romano (2008) and Theorem 1 in Hemerik and Goeman (2018a), this yields the null invariance Remark 3 The previous result is still valid in the case when w = 2 n , i.e., when not every element of {1, −1} n is used once (see, e.g., Hemerik and Goeman 2018a).For computational reasons, it is in fact common practice to sample uniformly from {−1, 1} n , with or without replacement.The same holds also for the results in the next section.

On nature of the eigen sign-flip test
A few comments on this approach may be useful to understand its nature.
In the simpler context of linear regression models, it is possible to define a residualizing matrix that projects into the residual space.This is an orthogonal projection matrix; as such, it is idempotent, and all its eigenvalues are zero or one.Thanks to this, for classical linear regression models, Kherad-Pajouh and Renaud (2010) propose to pre-multiply the residuals by the semiorthogonal matrix defined by the eigenvectors corresponding to the non-null eigenvalues.This pre-multiplication transforms the n residuals into pseudo-residuals, reducing their cardinality.The number of these new pseudo-residuals is equal to the rank of the residualizing projection matrix (i.e., the number of non zeros eigenvalues, usually equal to the number of covariates).Being the remaining eigenvalues of residualizing matrix equal to one, the resulting pseudo-residuals are now independent and homoscedastic (i.e., spherical).In particular, Kherad-Pajouh and Renaud (2010) suggest the use of a permutation approach, while Solari et al. ( 2014) extend it to the more general framework of rotations matrices.
Unfortunately, within the semiparametric regression framework, the residualizing matrix is not a projection matrix and is not idempotent; therefore, its eigenvalues do not take values in {0, 1}.A multiplication by these eigenvalues would act as a scaling factor for the residuals, making them independent, but not homoscedastic.For this reason, defining as a permutation or a rotation matrix would not be a valid solution.Defining instead as sign-flipping matrix ensures the commutative property D = D , as highlighted in Remark 2. This property is indeed crucial, as it guarantees that the variance of the test statistics is constant over , as proved in Theorem 4.1.It is also worth to emphasize that the test guarantees only asymptotic exactness since the penalization induces a bias in the estimation of the mean-which vanishes with increasing n-while, for fixed n, the variance remains constant for all the test statistics defining the null distribution.On the contrary, the standard parametric Wald test needs asymptotic results for both the mean and the variance, such as those obtained in Sect.3.1.Similar considerations could be drawn for a naive sign-flip score test that does not decomposes the matrix .In this case, for finite samples, the signflipped test statistic would not have variance equal to observed test statistic.This would lead to performances that are comparable to the parametric Wald test, as shown in Ferraccioli (2020).This difference between the eigen sign-flip test and the other competitors is crucial in providing an adequate control of Type I error, as shown by the simulations in Sect. 5.

Interval hypotheses, two-sided tests and confidence intervals
So far we have defined the eigen sign-flip test for point-wise null hypothesis.The most common situation in practice is nonetheless to define interval null hypotheses.As for standard approaches, we need to prove that the p-value computed under any point-wise null hypothesis H 0 : β = β 0 − (∀ > 0) has also Type I error probability bounded by α.For convenience, let us define the test statistic as a function of the tested coefficient, that is . We now give results for interval hypothesis and two-sided tests, and a third result on confidence intervals.

Proof
We have Note that last inequality holds since Note that (I − ) is a diagonal matrix with non-negative diagonal entries, thus it is positive semi-definite for all .Therefore, any quadratic form of it is nonnegative.
Proof Theorem 4.2 proves that P This, together with the fact that ) → 0, proves the corollary.
As consequence of the two lemmas above, we can also derive confidence intervals for the parameter β.

Testing a subset of the covariates
We now with the case where we have multiple covariates, and we are interested in testing a subset of the covariates.Specifically, assume X ∈ R n ×R q represents the set of covariates of interest, with associated vector of coefficients β, and Z ∈ R n ×R p the set of covariates associated with the vector of nuisance coefficients ζ .The minimization problem in (2) then becomes argmin We might be interested in testing for any value of ζ and γ .Let us define * = [Z | ] the n × ( p + K ) matrix composed by the covariates associated with the nuisance parameters and the bases for the nonparametric part of the model, with coefficients θ = (ζ , γ ).We can then rewrite equation ( 5) as where O is a matrix of zeros.Definition 1 of the eigen sign-flip test remains valid also in this case, with the only modification of the matrix in (6), where is replaced by * .Moreover, the following corollary provides the extension of Theorem 4.2 to the case where β is a vector.

Corollary 4.2.4 Consider the test that rejects H
), where ϕ(•) is any nonparametric combining function (Section 6.2 Pesarin 2001).Then, under the null hypothesis, the test is asymptotically exact and the rejection probability )) is at most α.Proof In order to extend the proof of Theorem 4.2 to the multivariate framework, we need to rely on the Nonparametric Combination of dependent test statistics, as defined, e.g., in Section 6.2 of Pesarin (2001).First of all, recall that the test statistic T is a vector itself.Moreover, Theorem 3.2 proves the asymptotic multivariate normality of T and Theorem 4.1 shows that the sign-flipped vectors of test statistics T n I , . . ., T n w share the same distribution.Therefore, the matrix T = (T n I , . . ., T n w ) is equal in distribution to • T.More precisely, T is the 2 n -dimensional vector of test statistics T n (•) , i.e., each row of T is a sampling from the multivariate test statistics T .We can therefore use any nonparametric combining function Pesarin (2001) to obtain a pvalue.
Among the most commonly used nonparametric combining functions, defined, e.g., in Pesarin (2001), are the max-T, sum-T or Mahalanobis distance.As an illustrative example, a p-value based on the min-p combining function (Westfall and Young 1993) rejects the multivariate hypothesis if the maximum value of T is larger than the 1 − α quantile of the distribution of the maxima computed over the w elements of (T n I , . . ., T n w ).

Simulation studies
In this section we present two simulation studies, to investigate the finite sample performances of the proposed test.Simulation 1, in Sect.5.1, considers a semiparametric model based on classical univariate splines (as, for instance, in Heckman 1986;Wand and Ormerod 2008).Simulation 2, in Sect.5.2, considers instead a semiparametric model based on SR-PDE (Sangalli et al. 2013).In these different settings, we compare the performances of three different tests: • Wald: a classical Wald-type test based on the asymptotic distribution of β; • Speck: a similar Wald-type test based on the asymptotic distribution of the Speckman version of the estimator (Speckman 1988), as derived in Holland (2017); • ESF: the Eigen sign-flip score test introduced in Definition 1.
The results show the performances of the tests over 1000 simulation repetitions.
The covariates and the true f are standardized, before computing the response variable y, so that their relative contributions to the response are comparable.We consider both β 0 = 0 and other 10 different values of β 0 , from 0.01 to 0.1, to check both the Type I error and the power of the test.Finally, we add i.i.d.normal random errors 1 , . . ., n , with zero mean and standard deviation 0.1.For each test case, the generation of the covariates and noise is repeated 1000 times.
The model is estimated using cubic B-spline bases, with 200 equispaced internal nodes on , using the implementation in Wand and Ormerod (2008).The smoothing parameter is chosen via cross-validation.The tests are performed with nominal value 0.05.For the proposed eigen sign-flip test, we consider 1000 random sign-flips.Table 1 shows the control of Type I error, and Fig. 1 shows the power functions for the three competing tests.The table and figure immediately highlight that the most challenging scenarios are cases (b) and (d), where the covariates have been generated with a dependence structure, sampling from a Gaussian process.The classic parametric test (Wald) shows an extremely poor control of the Type I error in these two cases, with an observed proportion of Type I error of over 26%, when the nominal value of the test is 5%.behavior is possibly due to the poor estimation of the variance induced by the regularized estimates.The Speckman variant appears more robust, partly correcting for the misspecified variance.Nonetheless, this test is significantly underconservative in cases (b) and (d), with a proportion of Type I error of almost 9%, while it is over-conservative in cases (a) and (c), where it returns a proportion of Type I error of about 2.3%.The proposed eigen sign-flip score test, on the contrary, maintains an extremely good control of the Type I error, under all scenarios, and it is never underconservative.Also in the challenging cases (b) and (d), at the cost of a slightly loss of power, it manages to keep a proportion of Type I error very close (and just slightly inferior) to the nominal value of the test.
We also considered the case of multiple covariates, following the simulation scheme detailed above, but including simultaneously all four covariates (a), (b), (c) and (d) in the data generation, and testing one parameter at a time, considering the other parameters as nuisance, as detailed in Sect.4.4.The same considerations as those detailed for the simulation in Fig. 1 can be drown (results non included for sake of space).

Simulation 2
In Simulation 2, we simulate from model (1), with = [0, 1] × [0, 1] and p 1 , . . ., p n randomly sampled from a uniform distribution on , with n = 225.For the nonparametric component of the model, we consider the test function 2 from the function gamSim in the R package mgcv (Wood 2015(Wood , 2017)), defined as We consider q = 1 covariate, and we generate x 1 , . . ., x n according to four different stochastic processes: 2 ) 2 with added a Matern random field with ν = 1, σ = 2 and scale 0.1.The covariates and the true f are standardized, before computing the response variable y, so that their relative contributions to the response are comparable.We consider both The model is estimated using SR-PDE, with linear finite elements on a mesh having 225 nodes on a regular lattice over , implemented using the package fdaPDE.The smoothing parameter is chosen via cross-validation.The tests are performed with nominal value 0.05.For the proposed eigen sign-flip test, we consider 1000 random sign-flips.
The results are presented in Table 2 and Fig. 2. The classic parametric test (Wald) has poor performances and very low control of the Type I error in all the scenarios, with proportion of Type I error of about 15% and higher.The Speckman variant is always more robust than the Wald, but it is often severely underconservative, with observed proportion of Type I error of about 10%.The proposed eigen sign-flip, on the contrary, at a loss of some power, permits an extremely good control of the Type I error, even in the more challenging scenarios, where the covariate has a strong spatial structure.

Study of human development in Nigeria
In this section we apply the proposed methodology to the analysis of human development in Nigeria.In particular, we are interested in better understanding the difference in socioeconomic and health conditions in the various states of the country.Unfortunately, data at national and subnational level are often poor or not publicly available.This lack of information and of public domain surveys hamper the efforts to identify and develop targeted interventions in troubled areas (Jerven 2013).An alternative to traditional data consists in using other sources of openly accessible data, such as data from social media, mobile phone networks, or satellites.In particular, a popular recent approach leverages on satellite images of luminosity at night to estimate economic activity (Chen and Nordhaus 2011;Jean et al. 2016).These images highlight urban areas, which typically offer better provisions of basic services such as electricity, water and public health, as well as more job opportunities, with respect to rural areas.
Here we use open satellite data (NASA Worldview Snapshots), together with demographic data, to predict human development.Specifically, as a response variable, we consider the Human Development Index (HDI) (available at https://globaldatalab.org/ shdi), an aggregated index that takes into account multiple dimensions at the household  Since the HDI, the response variable, and one of the covariate, the population density, are available at state level, we also aggregate the other three covariates at state level, considering their areal means.We then apply SR-PDE, considering the data located at the capitals of each state.We use a mesh with 320 nodes and select the smoothing parameter through generalized cross-validation (λ n = 0.1).We hence perform significance tests on each covariate, one at a time, considering the other parameters as nuisance, as described in Sect.4.4, using the eigen sign-flip procedure with 5000 random sign-flips.
Nightlight results significant ( p < 0.005), with a positive impact on human development (the estimated coefficients is 0.29).The finding on nightlight is in line with other recent research studies (Chen and Nordhaus 2011;Jean et al. 2016).The presence of urban areas, in fact, plays a huge role in the overall wealth of the population.This of course does not imply a causality effect, since increased wealth has itself an impact on the development of urban areas.Nightlight is nonetheless a good indicator of the socioeconomic status at local level, that does not require any official statistics, as previously discussed.Short-wave infrared seems to be slightly significant (0.05 < p < 0.1), with a negative impact on human development (the estimated coefficients is −0.016).The result might suggest that the presence of deserted areas with large amounts of bare soil lead to a decrease in human development.The more advanced states are indeed close to the ocean, in the southern part of the country.The northern part instead, that is mostly deserted, is not very populated.It is also worth noting that the aggregation at state level averages localized features, such as the presence of rivers, lakes or small vegetation, possibly reducing important information.The third satellite covariate, near infrared, does not appear significant ( p > 0.1).The same apply for population density ( p > 0.1).This is possibly due to the fact that the distribution is highly skewed, with most of the population residing in the state of Lagos, in the southwest of the country (see panel d in Fig. 3).Panel f in Fig. 3 shows the predicted MPI values, highlighting the very high explicative power of the model.

Discussion
This paper describes a strongly innovative and highly promising inferential approach in the context of semiparametric regression with roughness penalties.The paper focuses on tests for the linear part of the models.On the other hand, similar ideas can be used to develop tests and confidence bands on the nonlinear part of the models.Moreover, the described approach could be extended to deal with semiparametric regression with spatiotemporal components [see, e.g., Ugarte et al. (2009Ugarte et al. ( , 2010)); Aguilera-Morillo et al. (2017); Marra et al. (2012); Augustin et al. (2013); Bernardi et al. (2018)], further broadening the spectrum of potential models that could benefit from our proposal.These developments will be objects of dedicated future studies.We are confident this inferential approach will become popular and will prove to be highly valuable in the varied contexts where semiparametric regression is used.ment.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material.If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Fig. 2
Fig. 2 Simulation 2. Power of the Wald test (green dotted), of its Speckman variant (cyan dashed) and of the proposed eigen sign-flip (red solid)

Fig. 3
Fig. 3 Panels a-d show the covariates used: nightlight luminosity (a), Short-wave infrared (b), Near infrared (c) satellite images, and population density (d).Panel e shows the observed HDI.Panel f shows the HDI predicted by the SR-PDE model.Imagery from the Worldview Snapshots application (https://wvs.earthdata.nasa.gov For a given value of the test statistic T n I , we hence consider all the associated sign-flipped values T n I , T n 2 , . .., T n w , where we use the superscript n to highlight the sample size.We denote by T n(1) ≤ . . .≤ T n (w) the corresponding sorted value.Finally, we write T n [1−α] = T n ( 1−α w) .

Table 2
0 = 0 and other 10 different values of β 0 , from to 0.1, to check both the Type I error and the power of the test.Finally, we add i.i.d.normal random errors 1 , . . ., n , with zero mean and standard deviation 0.1.For each test case, the generation of the covariates and noise is repeated 1000 times. β