Kernel regression results in diverse outputs I am using kernel regression to build a prediction model. For the same, I am using np package. It is working fine, but I observed in multiple runs on the same data, it produces different results. Why it results in diverse outputs on the same data? Is there any way to select the best run of the model? Here is the minimal R code:

library(np)
   bw.all = npregbw(formula=power ~ temperature  
                                  + prevday1 + prevday2
                                  + prev_instant1 + prev_instant2 
                                  + prev_2_hour,
                regtype="ll",bwmethod="cv.aic", data=new_tr_dat)
   model.np &lt;- npreg(bws=bw.all)
   summary(model.np) 


I am using following data for my experiments:

        power temperature   prevday1   prevday2 prev_instant1 prev_instant2 prev_2_hour
1   220.59680          38         NA         NA      648.3621     1392.2186    848.7299
2   584.06867          38  220.59680         NA     1012.6853      250.1150    434.7129
3   206.39849          40  584.06867  220.59680      169.9380      105.5796    127.7294
4   177.05559          39  206.39849  584.06867      167.6312      229.3927    249.9871
5   165.71996          41  177.05559  206.39849      214.8291      248.5378    247.0262
6   184.02724          44  165.71996  177.05559      256.9970      314.3742    485.5184
7   187.70557          43  184.02724  165.71996      125.6160      213.9993    174.0830
8   916.78484          43  187.70557  184.02724      668.2840      217.3451    423.8285
9   185.98017          42  916.78484  187.70557      295.7329      331.6580   1227.0293
10  490.42294          42  185.98017  916.78484      241.6590      249.0523    255.3110
11  703.92694          39  490.42294  185.98017      806.5259     1515.1619   1140.4415
12 2038.91747          37  703.92694  490.42294      232.5541      582.5105    632.7118
13  208.66049          26 2038.91747  703.92694      210.5353      217.5053    221.3938
14  281.89860          37  208.66049 2038.91747      796.4336      256.4664    603.0781
15  425.72868          32  281.89860  208.66049      250.6069      187.1751    260.0573
16   86.77193          36  425.72868  281.89860      174.1249      179.6437    164.4359
17  218.06322          39   86.77193  425.72868      223.6548      316.2230    322.8536
18  258.89159          43  218.06322   86.77193      233.4561      372.5123    256.8588
19 1436.19980          40  258.89159  218.06322     1266.2630     1387.2287    791.7056
20  261.68520          42 1436.19980  258.89159      278.3378      230.5614    262.0084
21  225.34517          44  261.68520 1436.19980      211.3332      147.6705    196.8328
22  852.68835          44  225.34517  261.68520     1271.5826     1233.7158    991.7835
23 1729.79826          44  852.68835  225.34517      945.6528      298.0929    412.2199
24  464.58053          43 1729.79826  852.68835      182.6507      184.3031    203.5395
25  902.30950          45  464.58053 1729.79826      308.1398     1743.3495    642.4563
26  428.18792          45  902.30950  464.58053      205.1806      697.9208   1434.5425
27 1508.74739          43  428.18792  902.30950     1371.0550     2165.7173   1918.5236
28  355.01704          42 1508.74739  428.18792     1750.3907     1740.4654   1022.5056
29 3248.62618          43  355.01704 1508.74739      686.8528      360.0539    660.6378
30 1949.63937          44 3248.62618  355.01704      258.4627      217.2683    232.3818
31  725.25368          40 1949.63937 3248.62618     1406.3282     1714.6412   1375.2824
32  261.31252          32  725.25368 1949.63937      553.0443      275.6697    409.9598
You have a very small number of observations (32?) and a non-trivial number of predictors. 

It is known that the cross-validation function possesses multiple local minima/maxima. If you increase the number of multistarts to, say, nmulti=100 (add this option to your call to npregbw()), you ought to see that the same results occur on each invocation of the optimization process. 

Note you can do this all in the call to npreg() and skip the bandwidth call for convenience (npreg() will call npregbw() automatically but accept the arguments intended for npregbw()). 

Also, you will get the same results if you restart R each time you run the routine (seeds are set automatically to ensure this).

model &lt;- npreg(power ~ temperature  
                     + prevday1 
                     + prevday2
                     + prev_instant1 
                     + prev_instant2 
                     + prev_2_hour,
                regtype="ll",
                bwmethod="cv.aic", 
                nmulti=100,
                data=new_tr_dat)


To see whether things are stable with respect to the number of multistarts, look at the value of the cross-validation function and summary provided. Also, you can look at partial regression plots along with resampled variability bounds via plot().

summary(model$bws)
plot(model,common.scale=FALSE,plot.errors.method="bootstrap")


Also, I note that your predictor `temperature' is discrete, so you might consider using

+ ordered(temperature)


(i.e. use a discrete support kernel). Doing so reveals that there is little signal in this model, but the same holds for a simple parametric model (adjusted r-squared is negative).

model.lm &lt;- lm(model$bws$formula)
summary(model.lm)


Hope this helps!