Spline regression is one of the most used tools for nonparametric and semiparametric regression and has been extensively developed along the last decades. The main ideia is to select some knots on the function domain to define a partition and then fit a fixed degree polynomial on each segment in such a way that these polynomials smoothly match at the internal knots through continuity conditions imposed on their derivatives.
The result is typically a smooth curve fitting, which accurately estimates a large variety of smooth functions.
Further, spline regression offers a great flexibility in curve estimation problems according to convenient choices of degree and internal knots and, as generally occurs in nonparametric regression by basis expansion, the infinite estimation problem becomes a finite parametric estimation problem, once the spline function is determined by the estimation of its coefficients. For more details about spline based methods, their theory and applications in Statistics, see Wahba (1990) , Ramsay (2004) , Friedman et al. (2001) , and James et al. (2013) .
One of the main issues in spline regression is the selection of the number and positions of the knots. In some data analyses, wrong specification of the knots can cause a high impact on the spline regression performance.
Excessive number of knots may lead the spline regression to overfit the dataset while a not enough number of knots can lead to underfitting. Moreover, the positions of the knots should be well determined, since high concentration of knots on a specific interval of the domain overfits the data on this region and underfits the data outside it. Then, a good procedure for knots selection should provide optimal number and locations of knots, in order to have a well performed spline curve fitting.
In practice, the selection of knots is done by hand and few studies in the literature provide some systematic method to overcome this difficulty. Friedman and Silverman (1989) , Friedman (1991) and Stone et al. (1997) studied stepwise based methods for a set of candidates of knots. Osborne et al. (1998) proposed an algorithm based on the LASSO estimator to estimate the locations of knots. Ruppert (2002) presented two algorithms for selection of a fixed number of knots on penalized spline regression, called P-splines. The first one is based on the minimization of the generalized cross validation statistic and the second one is based on the Demmler-Reinsch type diagonalization andÜlker and Arslan (2009) proposed an artificial immune system to perform knots selection in curve and surface B-splines estimation by considering the knots locations candidates as antibodies. Bayesian methods are proposed by Denison et al. (1998) who proposed a joint prior distribution for the number and locations of knots, and Biller (2000) , who gave an adaptive bayesian approach for knot selection based on reversible jump Markov chain Monte Carlo in semiparametric generalized linear models.
Although the methods described above have been well succeeded, they work with a fixed number of knots and/or require some prior information regarding possible regions where the knots should be set, which can be seen as a limitation in some cases. In this work we propose an automatic method to estimate both the number of knots and their locations, based on the minimization of a penalized least squares criterion. The penalty plays the role of avoiding a high number of knots and consequently to overfit the data. We test the performance of the method on simulated data and we also apply it to daily reported cases of covid-19 on several countries to obtain a fit on the epidemiological curve. This paper is organized as follows: Section 2 is addressed to the description of the spline model and the proposed method of automatic knots selection. Simulation studies about the performance of the method are shown in Section 3. Applications to real datasets involving daily reported cases of Covid-19 are done in Section 4. Section 5 provides conclusion and suggestions of future works.
We start with a nonparametric regression problem of the form
where x i are scalars, f is an unknown smooth function with domain [a, b] such that f (x i ) = E(y i |x i ) and
{ i } n i=1 are zero mean independent random variables with unknown variance σ 2 . To estimate f , we define a sequence of knots a = t 0 < t 1 < ... < t K < t K+1 = b that defines a partition of [a, b] in intervals [a, t 1 ), [t 1 , t 2 ), ..., [t K , b] and use a m−th order spline regression model, m = p + 1, defined by
where p, K ∈ Z + , β = [β 0 , ..., β p+K ] is the vector of coefficients, (u) p + = u p I(u ≥ 0) and t 1 < ... < t K are fixed internal knots which, from the now on, will be referred to as just knots. In general, a spline regression model is a p degree polynomial on each interval of two consecutive knots and has p − 1 continuous derivatives everywhere (supposing that the knots have no multiplicity). In this work we adopt a set of knots sufficiently spaced, i.e, we assume |t k+1 − t k | > δ for some δ > 0, k = 0, ..., K. Thus, the problem of estimating f becomes a finite parametric estimation problem, whose parameters are the vector of coefficients β, the number of interna knots K and their locations t k , k = 1, ..., K. Let us denote by θ the vector of parameters to be estimated, i.e, θ = [β, K, t 1 , ..., t K ] . We assume a fixed value for p, the degree of the spline function.
The estimation of the parameters should consider its performance regarding a loss function, such as in our case, the the residual sum of squares and also some degree of smoothness for f . We propose to control the smoothness with a roughness penalty defined in terms of the number of intervals established by the knots, i.e, the K knots define K + 1 intervals on the domain of f . Thus, we define the penalized sum of squares
where y = [y 1 , ..., y n ] , λ > 0 is a tuning parameter previously chosen and f is the spline function (2.1). Thus, the estimatorθ is obtained byθ
Considering the definition above, we search for a curve estimate that has good fitting properties with low residual sum of squares and which is also smooth enough according to a reasonable number of knots that is controlled by the penalty on the number of segments of the domain. The tuning parameter λ attributes a weight on the roughness penalty term. High λ values increase the importance of the penalization, i.e, a preference for smoothness is given against accuracy. On the other hand, low λ value favors models with high number of parameters, which can lead to overfitting. The specification of a good value for λ is an issue inherent to any regularized estimator and, as in other similar approaches, could be chosen by a standard cross validation procedure.
Our proposed method for automatic knots selection can be applied to other spline basis structures. A well known spline basis is the B-splines, proposed by De Boor et al. (1978) , which are more adapted for computational implementation due to its compact support. The i-th B-spline of order m and knots a = t 0 <
is a spline of order m defined on the same knots such that is nonzero over at most m consecutive subintervals, i.e, over [t i , t i+m ] and zero outside it and can be defined recursively by
. Then, f is expanded in the B-splines basis instead of (2.1) as Another commonly used spline basis structure is the natural spline. This approach imposes the function f to be linear on its extremes, i.e, in regions lesser than t 1 and greater than t K . Such constraint guarantees stable estimation on the boundaries of the function, see James et al. (2013) . For more about spline basis structures, see De Boor (2001) .
To illustrate the role of the penalization on the curve estimation, we generate data from an underlying function defined by cubic B-splines basis on [0, 100] with three not equally spaced knots at t 1 = 20, t 2 = 45, and t 3 = 80 and signal to noise ratio SNR=3. Figure 2 .2(a) shows the underlying function, its knots locations and the generated data. After data generation, we estimate the underlying function using cubic B-splines expansion considering eleven knots with prespecified locations fixed at [t 1 , t 2 , · · · , t 11 ] = [6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26] . Note that we intentionally set several knots at one specific region of domain, around the higher peak of the curve. The estimated curve is the pink line in the knots at [t 1 , t 2 , · · · , t 11 ] = [6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26] (pink) characterizing locally overfitting around the higher peak and undefitting on remaining domain and estimated curve by the proposed method (blue) with optimal number and locations of knots.
The proposed method was applied in simulated data to evaluate its performance in several scenarios. We considered the cubic B-splines with one knot at t 1 = 50, two knots at t 1 = 25 and t 2 = 75 and three knots at t 1 = 25, t 2 = 50 and t 3 = 75 as underlying functions. These functions are shown in Figure 3 .1. Table 3 .3: Simulation results of the proposed method for data generated from a cubic B-splines with three knots at t 1 = 25, t 2 = 50 and t 3 = 75. 
The counts of daily new cases of the novel coronavirus disease (covid-19) varies considerably throughout the affected countries and regions, mainly due to testing criteria used, delays in notification, and population and authorities response to the outbreak, among others. Even though, the analysis of the disease's daily new infected cases is important to understand at what stage of the epidemic each country is in. Therefore, we can use the method presented above to estimate the dates at which the trend of new cases counts had changed.
This may be helpful not only to define public policy such as when safely relax social distance measures but also to verify the effectiveness of mitigations measures taken. Here we perform estimation and prediction for Covid-19 data using our proposed method, the data is publicly available at https://opendata.ecdc.
europa.eu/covid19/.
In general, during an epidemic the daily number of new infected cases shows several trends that are not related to the disease spreading itself, but to other causes, such as delay in test processing time or lack of testing. For this reason, we considered the 7-days moving average smoother on the time series before applying our method. Further, to avoid knots estimation in the beginning of the time series, we remove the 10% initial dates as possible knots locations.
Considering data available until May 30th, In the European countries considered here, we note that mean daily new cases present a decreasing pattern as the peak occurred in late March, period in which our method selects some knots, indicating a change in trend. On the other hand, the daily new reported cases are in a high plateau in the USA, since the righthand side of the curve is nearly parallel to the horizontal axis of Figure 4 .2(b). Moreover, the first change point in late March indicates a decrease in the slope of the fitted curve, the remaining knots are located near the date the country reaches its peak so far. In Brazil, the data suggest the country have not reached its peak of infection yet, change points merely reveal small changes in the increasing pattern. Finally, Iran and South Korea both show oscillations during the observed period, however at latest dates the new cases are apparently increasing, suggesting a new wave of reported cases. We also perform similar analysis for several other countries and for Brazilian states and municipalities, that are updated daily, and the results and coded algorithms in R can be accessed at the webpage https://www.ime.usp.br/~gpeca/covid-19.
We mention that the principal aim of this section is solely to show a descriptive analysis of the Covid-19 situation in the countries based on public available datasets and the proposed automatic method of knots estimation as a tool to do so. For a public policy guide to deal with this sad pandemic context, several social and biological factors beyond the scope of our study should be considered on a deeper analysis. 
In this paper we introduced a method to estimate the number and position of the knots of a spline regression function. The method is based on the minimization of the sum of squared residuals plus a penalty that depends on the number of knots. We evaluate the performance of the criterion on simulated data, and we showed that our proposed method has a great performance in the simulation studies, even with low SNR values. We also applied the method to perform a descriptive data analysis on covid-19 daily reported cases in several countries. In this analysis we show that the penalized least squares estimation guaranteed quite smooth curve fittings that estimated accurately underlying smooth functions automatically, i.e, without requiring to specify the number of internal knots and their locations, which is necessary in most spline based curve fitting methods available in the literature.
The penalizing constant λ used in our data analyses was set to a fixed value. However as in many other similar approaches, it could be chosen by a cross validation procedure. From a theoretical point of view, one open question that remains is what is the rate for the penalizing constant λ in order to obtain consistency of the estimatorθ. Some results on this direction for related penalized models is presented in Castro et al. (2018) , where a consistency result was proved for a penalizing constant of order n −1/2 . Whether this also holds on this setting is an issue that we would like to address in future work.