A method for determining groups in nonparametric regression curves: application to prefrontal cortex neural activity analysis

http://www.aimspress.com/journal/mbe MBE, 19(7): 6435–6454. DOI: 10.3934/mbe.2022302 Received: 03 March 2022 Revised: 21 March 2022 Accepted: 28 March 2022 Published: 24 April 2022 Research article A method for determining groups in nonparametric regression curves: Application to prefrontal cortex neural activity analysis Javier Roca-Pardi˜ nas1, Celestino Ord´ o˜ nez2,*and Lu´ ıs Meira Machado3 1Department of Statistics and Operational Research, Vigo University, Vigo 36310, Spain 2Department of Mining Exploitation and Prospecting, Geomatics and Computer Graphics Lab, Oviedo University, Mieres 33600, Spain 3Center of Mathematics, Minho University, Braga 4704-553, Portugal *Correspondence: Email: [email protected]; Tel: +34985458027. Abstract: Generalized additive models provide a flexible and easily-interpretable method for uncovering a nonlinear relationship between response and covariates. In many situations, the effect of a continuous covariate on the response varies across groups defined by the levels of a categorical variable. When confronted with a considerable number of groups defined by the levels of the categorical variable and a factor-by-curve interaction is detected in the model, it then becomes important to compare these regression curves. When the null hypothesis of equality of curves is rejected, leading to the clear conclusion that at least one curve is different, we may assume that individuals can be grouped into a number of classes whose members all share the same regression function. We propose a method that allows determining such groups with an automatic selection of their number by means of bootstrapping. The validity and behavior of the proposed method were evaluated through simulation studies. The applicability of the proposed method is illustrated using real data from an experimental study in neurology. Keywords: clustering of regression curves; factor-by-curve interaction; generalized additive model; multiple regression curves; nonlinear regression; number of groups 1. Introduction One of the main goals of statistical modeling is to understand the effect of some explanatory variables, Xi,i=1, ..., p, on a dependent variable, Y, also known as the response. This type of dependence is often modeled using a generalized linear regression model (GLM) [1] that imposes a linear relationship between explanatory and dependent variables without specifying in advance the function that links 6436 them. Generalized additive models (GAMs) [2, 3] offer an extension of the GLM through the incorporation of nonlinear forms for the explanatory variables. GAM has been widely used as an effective technique for conducting nonlinear regression analysis in various fields such as economics, finance, the environment, medicine, and biology. In a wide range of these applications, problems can be observed which it might be useful to compare two or more regression curves. This often occurs when it is necessary to check if the effect of a continuous covariate on the response varies across groups defined by levels of a categorical variable. Testing the hypothesis of equality of the regression functions is a topic of statistical inference that has been widely researched in the literature. One important review on this topic may be found in [4] (see their Section 7). Examples of methods for clustering regression curves using L2distance and other metrics are in [5–7]. For nonparametric models, most approaches rely on smoothing techniques, such as splines or Nadaraya-Watson [8], in the construction of the nonparametric estimators of regression functions. The general idea is to use these nonparametric estimators directly, either by contrasting all individual estimators with the pooled estimator, or by performing pairwise comparisons. Other authors considered the use of empirical processes to avoid the necessity of selecting a smoothing parameter required in the construction of the nonparametric estimators of regression functions [9–14]. The idea behind these tests is to use the distribution of the residuals and to compare them via Kolmogorov-Smirnov or Cram´ er-von Mises type statistics. ANOVA-type test statistics have also been used to compare regression functions (see for example [15] and [16]). Reference [17] illustrates how the SiZer exploratory tool is capable of comparing multiple curves based on the residuals. In [18] a kernel-based nonparametric approach is considered while [19] proposed a testing procedure for single index models. A second, though related question is how to determine groups among a series of regression curves when the null hypothesis of equality of the regression functions is rejected. Though the aforementioned methods can be used to compare regression curves, to the best of our knowledge, those methods cannot be used to determine groups among a series of regression curves, for example, among groups defined by the levels of a categorical variable. Some of them are not recommended when confronted with a considerable number of curves. If the null hypothesis of equality of curves is rejected, then this leads to the clear conclusion that at least one regression curve is different. However, these methods cannot be used to ascertain whether individuals can be grouped into a reduced number of classes whose members all share the same regression function or if all the regression curves are different from each other. One na¨ ıve approach would be to perform pairwise comparisons. However, this approach would lead to a large number of comparisons (e.g., 7 groups would lead to 21 pairwise comparisons). This could be done but without the possibility of determining groups with similar regression curves. A statistical procedure to estimate the unknown group structure was proposed in [20] and later in [21]. The first paper of these authors was about financial time series whereas the second deals with environmental statistics. More recently, the comparison of survival curves was addressed [22]. Similarly, in this paper, we propose an approach that allows determining groups that share the same regression function with an automatic selection of their number. The proposed method can be used, for instance, to establish groups with similar monotonic trends in the regression curves (e.g., over the levels of categorical variables). The proposed methodology will be shown in the framework of the generalized additive model with a binary outcome and a logit link function. The generalization of the proposed methods to other link functions is possible. The remainder of the paper is organized as follows: The following section provides the notation Mathematical Biosciences and Engineering Volume 19, Issue 7, 6435–6454. 6437 and the methodological background. Then, the performance of the proposed methods is investigated through simulations, and their usage is illustrated through the analysis of a real data set from neuronal activity in the prefrontal cortex of monkeys during a discrimination task. Finally, the last section contains a discussion and the main conclusions of the work. 2. Materials and methods 2.1. Mathematical model Let Ydenote the binary (0/1) target response and X=(X1,...,Xp) a set of pcontinuous covariates. In this regression context, the logistic generalized additive model expresses the conditional probability P(X)=P(Y=1|X) as log P(X) 1−P(X)=α+ p X j=1 fj(Xj) (2.1) where αis a constant and each partial function fjrepresents the effect of the covariate Xj. These unknown smooth partial functions are modeled from the data without specifying in advance any parametric structure. Therefore, the GAM in (2.1) combines flexibility with interpretability, in which each of the additive components describes the influence of each covariate separately. In practice, the effect of a continuous covariate Xjcan vary across the levels 1,...,Kof a categorical factor F. For simplicity of notation, here we consider that only the effect of the last covariate, Xp, depends on the levels of F, so the pure GAM in (2.1) can incorporate this type of factor-by-curve relationship as follows: logit(F,X)=log P(Y=1|F,X) 1−P(Y=1|F,X)=α+ p−1 X j=1 fj(Xj)+                  g1(Xp) if F=1 g2(Xp) if F=2 . . .. . .. . . gK(Xp) if F=K (2.2) where g1,...,gKrepresent the specific effect of Xpfor each of the possible Klevels established by the factor F. 2.2. Estimating partitions The study of the partial functions g1,...,gKcan be useful in the comparison of two or more groups, which is an important problem associated with statistical inference. Interest centers on the null hypothesis H0:g1=. . . =gK, namely, that the effect of Xpdoes not depend on the levels of the factor F. When the equality of the Kcurves is rejected, leading to the clear conclusion that at least one curve is different, it can be interesting to ascertain whether groups can be performed or, by contrast, that all these curves are different from each other. More clearly, if H0is rejected, it could be interesting to determine: a) if any pair of curves are different from each other, that is gi,gj, for all i,j∈ {1,...,K}, Mathematical Biosciences and Engineering Volume 19, Issue 7, 6435–6454. 6438 or b) if the indices in {1,...,K}can be grouped in a partition of J(J<K) groups I=(I1,...,IJ) with J<K, so that gi=gjfor each i,j∈Ik, for some k. Note that Imust be a partition of {1,...,K}and, therefore, must satisfy I1∪. . . ∪IJ={1,...,K}and Ii∩Ij=∅for all i,j∈I. In some situations, the possible partition I=(I1,...,IJ), can be established in advance, and the interest is to validate such partition by testing the null model logit(F,X)=α+ p−1 X j=1 fj(Xj)+                 g1(Xp) if F∈I1 g2(Xp) if F∈I2 . . .. . .. . . gK(Xp) if F∈IJ (2.3) versus the general model given in (2.2). However, in general, the partitions are unknown and will have to be estimated from the data. Before explaining the procedure for determining the groups, let us introduce some notation. Given a partition I=(I1, . . . IJ), we define a function GI:{1,...,K} → {1,...,J}so GI(k)=jif k∈Ij. Now, given a sample {Fi,Xi,Yi}n i=1and using the estimates ˆgk(k=1,...,K) of the model in (2.2), the estimated partition ˆ I=(ˆ I1, . . . ˆ IJ) can be obtained by minimizing the following Cram´ er-von Mises type distance min I=(I1,...,IJ) K X k=1Zxˆgk(x)−ˆ CGI(k)(x)2dx (2.4) Alternatively, the estimated partition can be obtained using the Kolmogorov-Smirnov type distance min I=(I1,...,IJ) K X k=1Zxˆgk(x)−ˆ CGI(k)(x)dx (2.5) In both cases the estimated centroids for j=1,...,Jare defined as ˆ Cj(x)=PK k=1ˆgk(x)·I{GI(k)=j} PK k=1I{GI(k)=j}(2.6) The above minimization problems can have a high computational cost for large values of K, because they require the evaluation of all the different combinations of the Kcurves into Jgroups. To solve the quadratic minimization problem in (2.4) and L1-norm minimization problem in (2.5) we propose the use of the k-means or k-medians algorithms. In both cases, the partial functions have to be estimated in a common grid of Xp,x• 1< . . . < x• N, of size N, leading to a matrix             ˆg1(x• 1) . . . ˆgK(x• 1) . . .. . . ˆg1(x• N) . . . ˆgK(x• N)            N×K This matrix will be the input of both heuristic methods, k-means and k-medians, and from these, the estimated partition ˆ I=(ˆ I1,...,ˆ IJ) is obtained. Mathematical Biosciences and Engineering Volume 19, Issue 7, 6435–6454. 6439 2.3. Determining the number of groups The procedure to obtain the estimation partition ˆ I=(ˆ I1,...,ˆ IJ) for a given J(J<K) has been explained in the previous section. This section focuses on determining the number of groups for the regression functions. With this goal in mind, one possible approach consists in fitting the model in (2.3) for each value of j=1,...,J, and selecting the number of groups Kthat minimizes the Bayesian Information Criterion (BIC). BIC is one of the most widely used tools in statistical model selection. Its popularity is derived from its computational simplicity and effective performance in many modeling approaches. Our simulations on curve clustering demonstrate that the BIC criterion can be considered a valid option, performing quite satisfactorily choosing the correct number of clusters in simulation scenarios with a large number of curves and complex data sets. Another possibility would be to test, for a given J, the null hypothesis H0(J) that at least one partition Iof length Jexists, so that the model in (2.3) is fulfilled, against the alternative that for any J-partition at least a group Ijexists in which gm,gnfor some m,n∈Ij. In order to test H0(J), for a given sample {Fi,Xi,Yi}n i=1, we consider the following two statistics DCM = K X k=1 n X i=1ˆgk(Xip)−ˆg◦ Gˆ I(k)(Xip)2(2.7) DKS = K X k=1 n X i=1ˆgk(Xip)−ˆg◦ Gˆ I(k)(Xip)(2.8) where the estimated partition ˆ I=(ˆ I1,...,ˆ IJ) is obtained by solving the minimization problems in (2.4) or (2.5), and ˆg◦ jfunctions are obtained from model (2.3) using the same partition. Alternatively, we can use the deviance as an appropriate measure of discrepancy between observed and fitted values and consider the following likelihood ratio test DLR =Pn i=1Devi(Yi,ˆ Pi)−Pn i=1Devi(Yi,ˆ P◦ i) Pn i=1Devi(Yi,ˆ Pi)(2.9) which compares the deviance of the null model (2.3) with the deviance of the general model in (2.2), where ˆ P◦ iand ˆ Piare the estimates values of the true probabilities Pi=P(Yi=1|Fi,Xi) under the null and general model, respectively. The individual deviance Devi(Yi,ˆ Pi) is defined as Devi(Yi,ˆ Pi)= −2Yilog ˆ Pi+(1 −Yi) log(1 −ˆ Pi). Note that if H0(J) is verified, the value of the test statistic Dshould be close to zero. The decision rule based on each of the three statistics, D, consists in rejecting the null hypothesis if Dis greater than (1 −α)-percentile obtained under the null hypothesis. Nevertheless, the theory for determining such percentiles is not closed. To approximate the distributions of the test statistics, resampling methods such as the bootstrap method can be applied [23–26]. The binary bootstrap used in our paper is a particular case of the bootstrap techniques suggested by [27] and [28] for inference in nonparametric models with response belonging to the binary family. The binary bootstrap involves the following steps: Step 1. Using the original sample data {Fi,Xi,Yi}n i=1, compute the test statistics Das explained above. Then, obtain for i=1,...,nthe estimated ˆ P◦ iof the true probabilities Pi=P(Yi=1|Fi,Xi) obtained from the null model in (2.3). Mathematical Biosciences and Engineering Volume 19, Issue 7, 6435–6454. 6440 Step 2. For b=1,...,B, generate the bootstrap sample {Fi,Xi,Y• i}n i=1, with Y• i∈Bernoulli(ˆ P◦ i), and compute the bootstrap statistics D•,b. Note that the estimation of the partition is required in this Step. Finally, the decision rule, based on each test statistic D, consists in rejecting the null hypothesis if D>ˆ D1−α, where ˆ D1−αis the empirical (1−α)-percentile of the values D•1,...,D•Bobtained previously. The bootstrap-based test introduced here can be very useful to automatically determine the number of groups J. Specifically, the rule proposed is as follows: the procedure begins by testing H0(1). If this hypothesis is true, we decide that J=1 (all regression curves are equal). When the previous hypothesis is rejected, it will be necessary to test H0(2). If this new hypothesis is not rejected, we decide that J=2. When H0(2) is rejected, it will be required to test H0(3) and so on until a certain H0(J) is not rejected. Finally, note that we are aware that this approach could deal with the problem of multiple hypothesis testing where a set of Jp-values corresponding to the Jnull hypotheses, H0(1),H0(2),...,H0(J), are given. Even though several methods have been proposed to deal with this problem (see e.g. [29] for an introduction to this area), there is still a challenge because there is no information about the minimum number of tests needed to apply these techniques. Accordingly, and considering that our rule finishes with a low number of them, we have not considered this problem. 2.4. Simulation study This section reports the results of a simulation study conducted to evaluate the practical performance of the proposed methodology. The response Ywas generated under model (2.3) with log P(Y=1|F,X) 1−P(Y=1|F,X)=                      g1(X)=1.5Xif F∈I1={1,2,3} g2(X)=2X2−3 if F∈I2={4,5} g3(X)=2 sin(2X)−2 if F∈I3={6,7} g4(X)=2 cos(2X) if F∈I4={8,9} g5(X)=2 cos(2X)+aexp(0.5X) if F∈I5={10,11,12,13} g6(X)=4−2Xif F∈I6={14,15} (2.10) The categorical covariate Xwas drawn from an uniform distribution U[−2,2], and its factors, F, were generated by taking a random value in {1,...,15}with associated probabilities p=(p1,...,p15)= n/Pnk,k=1,...,15, where: n=(n1,...,n15)=(1.0,2.0,1.0,1.5,1.0,2.0,1.5,2.0,2.0,1.0,1.0,1.0,1.5,1.5,1.0) In this way, and as can be seen in Table 1, we have considered unequal sample sizes for each (k= 1,...,15) curve. We explore the validity of the proposed tests assuming that we aim to test if the K=15 regression curves can be grouped in five (J=5) groups. To study the size and power of the proposed tests, different values were considered for a, ranging from 0 to 2. It should be noted that g5(X)=g4(X)+ aexp(0.5X), so a value of a=0 corresponds to a null hypothesis of five groups, while a value of a>0 corresponds to a null hypothesis of six groups. 2.5. Application to neural activity in the prefrontal cortex during a discrimination task The data analyzed in this section come from laboratory experiments of the extra-cellular single unit activity in the prefrontal cortex of a monkey. These experiments were conducted in the Laboratory of Mathematical Biosciences and Engineering Volume 19, Issue 7, 6435–6454. 6441 Table 1. Sample size nkfor each regression curve k(k=1,...,15) for different total sample sizes n=Pnk. n n1n2n3n4n5n6n7n8n9n10 n11 n12 n13 n14 n15 300 14 28 14 21 14 28 21 28 28 14 14 14 21 21 14 500 23 47 23 35 23 47 35 47 47 23 23 23 35 35 23 1000 47 95 47 71 47 95 71 95 95 47 47 47 71 71 47 2000 95 190 95 142 95 190 142 190 190 95 95 95 142 142 95 4000 190 380 190 285 190 380 285 380 380 190 190 190 285 285 190 Neurophysiology of the University of Vigo (Spain). The monkey was trained to discriminate between different stimuli. The stimuli consisted of stationary bright line segments presented on a monitor screen in front of the monkey that changed their orientation over time. A trial was initiated when the monkey pressed a lever key with its right hand and then two stimuli (reference and test), each of 500 ms duration, were presented in sequence with a fixed inter-stimulus interval (ISI: 1000 ms). At the end of the second stimulus, the subject released the key in a 1200 ms time window and pressed one of the two switches (left or right), indicating whether the orientation of the second stimulus was clockwise (right) or counter-clockwise (left) to the reference stimulus. While the monkey worked on the task, its extracellular unit activity was recorded. The monkey was rewarded for correct discrimination. To gather enough data and to account for the cell response variability, the neuron was recorded over a number of T=80 trials. For each of the trials, we consider the angle (trial specific), TA, corresponding to the test stimulus. The reference angle is 90◦. The following TAs were considered in the experiment: •T A ∈ {78◦,102◦}: test stimuli more separated from the reference, and therefore very easy to discriminate •T A ∈ {81◦, 99◦}: test stimuli is easy to discriminate •T A ∈ {84◦, 96◦}: test stimuli is difficult to discriminate •T A ∈ {87◦,93◦}: test stimuli closest to the reference, therefore more difficult to discriminate In this experiment, the outcome of interest is the neuronal activity in the interval t∈[tmin,tmax]= [−500,3000] in ms. At each instant tand trial j=1,...,T, this outcome may then be represented by a temporal binary sequence, Yj t, where Yj t=1 if there is a spike in [t,t+1) ms and 0 otherwise. The explanatory variables are the time twhen a spike was detected and the interval of time ∆tsince the last spike. The last one is an adjustment variable that takes into account that the neuron might fire more easily if it has been fired recently. Accordingly, for each trial, the data set consists of the following information {T A,(t,∆t),Yt}tmax t=tmin We use generalized additive modeling to assess whether the association between neural activity and decision-making depends on the difficulty of the discrimination task related to the angle of the test stimulus. We consider the following GAM: Mathematical Biosciences and Engineering Volume 19, Issue 7, 6435–6454. 6442 log P(T A,t,∆t) 1−P(T A,t,∆t)!=α+f(∆t)+                                      g1(t) if T A =78◦ g2(t) if T A =102◦ g3(t) if T A =81◦ g4(t) if T A =99◦ g5(t) if T A =84◦ g6(t) if T A =96◦ g7(t) if T A =87◦ g8(t) if T A =93◦ (2.11) where P(T A,t,∆t)=P(Yt=1|F,t,∆t), αis a fixed parameter, and gkfor k=1,...,8 are the specific time functions associated to each of the test angles considered in the study. 3. Results 3.1. Results on simulation data set Type I error rates and power values of the proposed tests as a function of aare shown in Figure 1. They were calculated from 1000 simulation runs using B=400 bootstrap samples at each repetition. We compare the results of the three test statistics under sample sizes n=1000 and 2000, at the significance levels of 0.05 and 0.10. As can be seen in Figure 1, the three curves show the expected behavior pattern, with an increase in the power as aincreases, and an improvement in it as the sample size grows. n=1000 a Power(%) DCM DKS DLR 0.00 0.50 1.00 1.50 2.00 0 20 40 60 80 100 n=1000 a Power(%) 0.00 0.50 1.00 1.50 2.00 0 20 40 60 80 100 n=2000 a Power(%) 0.00 0.50 1.00 1.50 2.00 0 20 40 60 80 100 n=2000 a Power(%) 0.00 0.50 1.00 1.50 2.00 0 20 40 60 80 100 Figure 1. Rejection probabilities (power of the hypothesis test) of the proposed test statistics as a function of a, for sample size n=1000 and n=2000 (upper and lower panels, respectively) at a 0.05 and 0.10 (left and right panels, respectively) significant levels (red line). Mathematical Biosciences and Engineering Volume 19, Issue 7, 6435–6454. 6443 The test statistic based on the likelihood ratio test (labeled as DLR) is that with the highest values of power for n=1000, but the differences between the three test statistics are minimal for n=2000. The Cram´ er-von Mises–type statistic test (labeled as DCM) revealed a poorer behavior with lower rejection probabilities for a sample size of 1000. Values for type I error, for the three test statistics and for different significance levels and different sample sizes (n=1000,2000, and 3000), are reported in Table 2. The three test statistics work satisfactorily according to the type I error, coming quite close to the nominal level regardless of the sample size. Table 2. Estimated type I error (in %) and nominal level percentages (1, 5, 10, 15, and 20) for different sample sizes. level DCM DKS DLR n=1000 1 1.4 1.3 1.0 5 5.5 5.3 5.6 10 9.7 9.8 11.5 15 14.5 15.3 16.5 20 18.7 20.0 21.1 n=2000 1 0.9 0.7 0.8 5 4.2 3.8 5.7 10 8.4 8.0 9.4 15 12.8 12.4 15.6 20 16.5 18.0 21.8 n=3000 1 0.9 0.7 1.2 5 4.2 4.0 6.5 10 8.8 7.6 11.1 15 13.0 12.1 16.2 20 16.2 17.0 22.3 In Table 3, we report results that can be used to evaluate the accuracy of the bootstrap-based algorithm introduced in Section 2.3. Again, different values of awere considered, ranging from 0 to 2. Recall that the value a=0 corresponds to the null hypothesis, which assumes that the fifteen regression functions can be classified into five groups, while, when a,0, the number of groups is six. Note that to select the correct number Jof groups of regression functions, the bootstrap-based algorithm must first reject the first null hypothesis, H0(1), then reject the second hypothesis, H0(2), and so on until it accepts H0(5) if a=0, or until H0(6) when a>0. Results shown in Table 3 for the three test statistics, using a nominal level of 5%, display the number of times that the procedure selects the number of groups J. Results in bold denote the correct classifications according to model (2.3), revealing the high accuracy of the proposed bootstrap-based algorithm for a=0, showing the correct number of groups in percentages quite close to the nominal level. As the value of aincreases, so does the percentage of cases in which the proposed method suggests six groups. When comparing the three test statistics, it can be observed that all have a similar performance with a small advantage for the test statistic based on the likelihood ratio test (labeled as DLR). Mathematical Biosciences and Engineering Volume 19, Issue 7, 6435–6454. 6450 than 90◦(clockwise rotation) were previously reported in a study [26]. Time(ms) Firing rate (spikes/s) −500 0 500 1000 1500 2000 2500 3000 0 5 10 15 20 25 30 Time (ms) Firing rate (spikes/s) TestReference ISI −500 0 500 1000 1500 2000 2500 3000 10 15 20 25 30 35 Figure 5. Representation of neural activity over time. Observed spikes pool counts of spikes within 10 ms intervals (left) and the corresponding smoothing version (right plot). Time (ms) Firing rate (spikes/s) 1500 2000 2500 3000 10 20 30 40 50 60 78º 81º 84º 87º 93º 96º 99º 102º Figure 6. Fitted regression curves for the eight different experimental conditions. Results for the bootstrap method using the three test statistics defined in Section 3 are shown in Table 7. As can be seen, the three test statistics indicate that there are two groups. DKD, and slightly less so DCM, are close to tipping the balance in favor of three groups. We have also applied the BIC criterion to determine the number of groups. The results of this cluster solution also lead to the same conclusion, namely, J=2 groups. These findings are in agreement with those reported in our simulations that indicate that BIC is more conservative, but capable of detecting Mathematical Biosciences and Engineering Volume 19, Issue 7, 6435–6454. 6451 important differences. Table 7. Probability values for testing the null hypothesis H0(J). Results based on the bootstrap method with different test statistics. JDCM DKS DLR 1<0.001 <0.001 <0.001 2<0.06 <0.05 <0.29 3<0.19 <0.13 <0.53 4<0.09 <0.76 <0.48 Clusters for a fixed number of groups between J=2 and J=5 are shown in Figure 7 to facilitate the comprehension of the problem, although the results obtained indicate that there are only two groups. Firing rate (spikes/s) 12345 78º 81º 84º 87º 93º 96º 99º 102º 2 groups Firing rate (spikes/s) 78º 81º 84º 87º 93º 96º 99º 102º 3 groups Time (ms) Firing rate (spikes/s) 12345 1500 2000 2500 3000 78º 81º 84º 87º 93º 96º 99º 102º 4 groups Time (ms) Firing rate (spikes/s) 1500 2000 2500 3000 78º 81º 84º 87º 93º 96º 99º 102º 5 groups Figure 7. Estimated regression curves according to the groups to which they belong. The curves assigned into two (J =2), until five groups (J =5) are represented in each of the four panels. The statistics used to estimate the number of groups found only two. 4. Conclusions In this work we propose a bootstrap-based method to determine the number of groups in generalized additive models when the effect of a continuous covariate on the response varies across groups defined by levels of a categorical variable. The simulation analysis confirms the capacity of the three proposed statistical tests to reproduce the theoretical values of type I and type II errors, even in the presence of a Mathematical Biosciences and Engineering Volume 19, Issue 7, 6435–6454. 6452 large number of curves. However, the power turned out to be higher for the test based on the likelihood test ratio. As expected, the percentage of misclassifications decreases with the sample size, n, being almost null for n=4000. The simulated data were also used to compare our method with BIC, a well-known approach that has been largely used previously in clustering analysis to estimate the number of clusters in an artificial data set. The number of detected groups was the same for both methods, but the simulation study showed that BIC has a lower statistical power than the bootstrap method, although they are close for larger sample sizes. Accordingly, BIC could be used to obtain an initial estimation of the number of clusters that could be tuned using our method. The application of the proposed method to an experiment to determine the activity of a neuron of a monkey subject to visual stimuli led to the same conclusions when we compared our approach with BIC. In both cases the neuron only reacts to two of the eight stimuli, corresponding to the most evident differences with respect to a reference state. Regarding the three statistics tested, they provide the same result, although the evidence for two clusters is stronger for the statistic based on the likelihood ratio test. Finally, it can be said that although the proposed method was designed to detect groups of regression curves in generalized additive models with a binary response, it can be extended without much difficulty to determine groups in models with another kind of response. Acknowledgments This work was partially supported by project 2017/00001/006/001/097: Ayudas para el mantenimiento de actividades de investigaci´on de institutos universitarios de investigaci´on y grupos de investigaci´on de la Universidad de Oviedo para el ejercicio 2021. Lu´ ıs Meira-Machado acknowledges financial support from Portuguese Funds through FCT - ”Fundac¸˜ ao para a Ciˆ encia e a Tecnologia”, within the projects UIDB/00013/2020, UIDP/00013/2020. Javier Roca-Pardi˜ nas acknowledges financial support from Grant PID2020-118101GB-I00, Ministerio de Ciencia e Innovaci´ on (MCIN/AEI /10.13039/501100011033). Conflict of interest The authors declare there is no conflict of interest. References 1. P. McCullagh, J. Nelder, Generalized Linear Models, 2nd edition, Chapman and Hall/CRC, Boca Raton, 1989. https://doi.org/10.1201/9780203753736 2. T. J. Hastie, R. J. Tibshirani, Generalized Additive Models, Chapman & Hall/CRC, New York, 1990. 3. S. Wood, Generalized Additive Models: An Introduction with R, Chapman & Hall/CRC, 2006. https://doi.org/10.1201/9781420010404 4. W. Gonz´ alez-Manteiga, R. M. Crujeiras, An updated review of Goodness-of-Fit tests for regression models, Test,22 (2013), 361–411. https://doi.org/10.1007/s11749-013-0327-5 Mathematical Biosciences and Engineering Volume 19, Issue 7, 6435–6454. 6453 5. H. Dette, A. Munk, Testing heterocedasticity in nonparametric regression, J. R. Stat. Soc. B,60 (1998), 693–708. https://doi.org/10.1111/1467-9868.00149 6. H. Dette, N. Neumeyer, Nonparametric analysis of covariance, Ann. Stat.,29 (2001), 1361–1400. https://doi.org/10.1214/aos/1013203458 7. L. Garc´ ıa-Escudero, A. Gordaliza, A proposal for robust curve clustering, J. Classif.,22 (2005), 185–201. https://doi.org/10.1007/s00357-005-0013-8 8. E. A. Nadaraya, On estimating regression, Theory Probab. Its Appl.,9(1964), 141–142. https://doi.org/10.1137/1109020 9. M. A. Delgado, Testing the equality of nonparametric regression curves, Stat. Probab. Lett.,17 (1993), 199–204. https://doi.org/10.1016/0167-7152(93)90167-H 10. K. B. Kulasekera, Comparison of regression curves using quasi-residuals, J. Am. Stat. Assoc.,90 (1995), 1085–1093. https://doi.org/10.1080/01621459.1995.10476611 11. K. B. Kulasekera, J. Wang, Smoothing parameter selection for power optimality in testing of regression curves, J. Am. Stat. Assoc.,92 (1997), 500–511. https://doi.org/10.1080/01621459.1997.10474003 12. K. B. Kulasekera, J. Wang, Bandwidth selection for power optimality in a test of equality of regression curves, Stat. Probab. Lett.,37 (1998), 287–293. https://doi.org/10.1016/S01677152(97)84155-7 13. N. Neumeyer, H. Dette, Nonparametric comparison of regression curves: An empirical process approach, Ann. Stat.,31 (2003), 31880–31920. 14. J. C. Pardo-Fern´ andez, I. Keilegom, W. Gonz´ alez-Manteiga, Testing for the equality of k regression curves, Stat. Sin.,17 (2007), 1115–1137. 15. S. G. Young, A. W. Bowman, Non-parametric analysis of covariance, Biometrics,51 (1995), 920– 931. https://doi.org/10.2307/2532993 16. J. C. Pardo-Fern´ andez, M. D. Jim´ enez-Gamero, A. Ghouch, A non-parametric ANOVA-type test for regression curves based on characteristic functions. Scand. J. Stat.,42 (2015), 197–213. https://doi.org/10.1111/sjos.12102 17. C. Park, K. Kang, Sizer analysis for the comparison of regression curves, Comput. Stat. Data. Anal.,52 (2008), 3954–3970. https://doi.org/10.1016/j.csda.2008.01.006 18. C. Park, J. Hannig, K. Kang, Nonparametric comparison of multiple regression curves in scale-space, J. Comput. Graphical Stat.,23 (2014), 657–677. https://doi.org/10.1080/10618600.2013.822816 19. W. Lin, K. B. Kulasekera, Testing the equality of linear single-index models, J. Multivar. Anal., 101 (2010), 1156–1167. 20. M. Vogt, O. Linton, Classification of non-parametric regression functions in longitudinal data models, J. R. Stat. Soc. Ser. B Stat. Methodol.,79 (2017), 5–27. https://doi.org/10.1111/rssb.12155 21. M. Vogt, O. Linton, Multiscale clustering of nonparametric regression curves, J. Econometrics, 216 (2020), 305–325. Mathematical Biosciences and Engineering Volume 19, Issue 7, 6435–6454. 6454 22. N. M. Villanueva, M. Sestelo, L. Meira-Machado, A method for determining groups in multiple survival curves, Stat. Med.,38 (2019), 866–877. https://doi.org/10.1002/sim.8016 23. P. Hall, J. D. Hart. Bootstrap test for difference between means in nonparametric regression, J. Am. Stat. Assoc.,85 (412), 1039–1049. 24. M. C. Rodr´ ıguez-Campos, W. Gonz´ alez-Manteiga, R. Cao, Testing the hypothesis of a generalized linear regression model using nonparametric regression estimation, J. Stat. Plan. Infer.,67 (1998), 99–122. https://doi.org/10.1016/S0378-3758(97)00098-0 25. J. Roca-Pardi˜ nas, C. Cadarso-Su´ arez, V. N´ acher, C. Acu˜ na, Bootstrap-based methods for testing factor-by-curve interactions in Generalized Additive Models: assessing prefrontal cortex neural activity related to decision-making, Stat. Med.,25(2006), 2483–2501. https://doi.org/10.1002/sim.2415 26. C. Cadarso-Su´ arez, J. Roca-Pardi˜ nas, G. Molenberghs, F. Faes, V. N´ acher, S. Ojeda, et al., Flexible modelling of neuron firing rates across different experimental conditions. An application to neural activity in the prefrontal cortex during a discrimination task, J. R. Stat. Soc. Ser. C,55 (2006), 431–447. 27. S. Sperlich, D. Tjøstheim, L.Yang, Nonparametric estimation and testing of interaction in additive models, Econom. Theory,18 (2002), 197–251. https://doi.org/10.1017/S0266466602182016 28. L. Yang, S. Sperlich, W. H¨ ardle, Derivative estimation and testing in generalized additive models, J. Stat. Plan. Infer.,115 (2003), 521–542. https://doi.org/10.1016/S0378-3758(02)00163-5 29. S. Dudoit, M. J. Van Der Laan, Multiple Testing Procedures with Applications to Genomics, Springer, Springer Series in Statistics, New York, 2007. ©2022 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0) Mathematical Biosciences and Engineering Volume 19, Issue 7, 6435–6454.