@David L Post

Thank you. I used Python to calculate R2.

R2 can be negative (it is reported by other researchers):

https://stats.stackexchange.com/questions/12900/when-is-r-squared-negative#:~:text=R2%20is%20negative%20only,Y%20intercept%20must%20equal%201500.

Thanks for the explanations. Unfortunately, I am not ale to change dataset since it is from a client.

Mohammad -

Whatever 50L and 2400L refer to, I don't see that you have any evidence to treat them separately. As David said, the range is too narrow for the 2400L.

Usually the predicted values would be on the x-axis, and it is the range of predicted-y which determines heteroscedasticity. I haven't used Gaussian processes regression, so I don't know if a smooth curve should fit it, but this looks odd. I would expect bars at each predicted value, in the observed value direction, and that you could see then if the observed values fell in there. I would also expect the predicted values to follow a straight line, or a curve, or perhaps have some discontinuities as in splines, but these points seem to be all treated separately. It's hard to picture sideways as it is. I see that this is a Bayesian technique, but I wonder if it is supposed to look this way on a graph. It appears that a slightly polynomial linear or nonlinear regression could also be used, and I wonder if results would be similar.

Back to that narrow range for the 2400L: If the change in predicted-y is small, a regression of any sort should not be very meaningful, as the observed values should also be about the same, appearing somewhat randomly based on the prediction interval. Sigma is a big part of all of this, a relatively large sigma giving you a relatively large R-squared, but you have to have a range of points to estimate R-squared. Still, for a short set of predicted values, you see a substantial sigma here, about like the others.

I don't think there is anything you can really say about the 2400L other than they seem consistent with the 50L in the short range in which you have them.

You might find the following to be of interest: https://data.library.virginia.edu/is-r-squared-useless/

They note that R-squared is not a measure of fit. I think that is because a model may be good for means at each predicted-y-value, but sigma may be high and give you a small R-squared.

Cheers - Jim

In regression with one or more predictors, if the value of y is expected to be zero when all predictors are zero, then there should be no intercept term. But I don't know what is happening here.

I got predicted y using Gaussian processes regressor (non-linear model). The plots I attached show predicted y vs observed y. The line shown on the plots is only for visualization purpose. It is not a regression line.

Something that is surprising for me is that why R2 is low while RMSE is low. I expected that lower RMSE means higher R2.

Considering what you have said, Mohammad -

If model errors are also very low, then a big portion of your overall (small) error could still be from the residuals. A graphical residual analysis set at a scale you can read would show this.

Because you are not graphing one variable against another, but "observed" versus "predicted," the relationship should be that line, or else your model is not model-unbiased, with expected residual zero. You should be putting "observed" (i.e., y) on the y-axis and "predicted" (i.e., predicted-y) on the x-axis. Are those bars around points for prediction intervals? They should be around predicted values on that line, not observed ones. So they should not jump around like that. In each case, the "observed" value would fall somewhere along the bar, or above or below it, along the same line as the bar, for a given predicted-yi.

As I said, I am not familiar with "Gaussian processes regression," but I can't see how this could be read differently.

When I mentioned curves and splines before, I made a mistake. I was thinking of how your variables may relate, but since your graphs are all observed and predicted only, points should all be about that line. Predicted-y may be different under different models, but one would expect to have the expected value of the estimated residuals for that sample to be zero for each and every alternative model. How well this holds up under cross-validation will differ. Some models will do better for the cases outside your sample than others will.

PS - Reviewing the question, apparently the bars are not prediction intervals, but just sigma (RMSE).  Ideally, for larger predicted-y, sigma should be larger, contrary to what is often assumed.  (See https://www.researchgate.net/project/OLS-Regression-Should-Not-Be-a-Default-for-WLS-Regression, and the various updates there.)  Ken Brewer explained why and to what extent this should occur.  See https://www.researchgate.net/publication/320853387_Essential_Heteroscedasticity.  When the ideal 'formula' for predicted-y is not found, this may not hold. I do not know about any special features of the regression used here.

Looking at your last graph, at RMSE (sigma), and R-squared, sigma is about the same for each of the two cases, but R-squared is better for the case where there were more data.  OK.  Sigma (standard deviation) is not reduced with a larger sample size, but the standard errors of the model coefficients are reduced with a larger sample size.  So, that is why you stated "Reason of low R2 is because of lower No. of [observations in the sample]."  (This assumes that the same model is used in either case, or a very similar one.) The smaller sample size made the regression coefficient standard errors larger, which would make a given sigma smaller in comparison. 

However, you cannot count on this all the time.  Typically with a small sample size, you will get a bad estimate of sigma (not to mention that ideally sigma is bigger for bigger predicted-y).  With a small sample, nothing is reliable.  With a very small sample, you could easily get an R-squared of close to one if just a couple of points are used and they 'line-up' by chance.  You could obtain a wide range of results. 

Here sigma looks close to unchanged with larger predicted-y. However, the axes are reversed, so I'm looking at the wrong one.  You might get about the same results, but you should correct that.  It is confusing. 

Please change your axes and recalculate everything.

Hmmm. About the smaller sample size and the lower R-squared. I said "The smaller sample size made the regression coefficient standard errors larger, which would make a given sigma smaller in comparison." But I think I'm mistaken in thinking that makes R-squared smaller, but instead, just more volatile. (A quick check on the internet and remembering back a little indicates I was wrong about R-squared.)

Thank you for comprehensive explanations.