Model.fit() gives the wrong 'rsquared' when 'weights' is given

[ydy1206]

First Time Issue Code

Yes, I read the instructions and I am sure this is a GitHub Issue.

Description

The ModelResult instance gained from Model.fit() has the wrong rsquared attribute. Floating point R^2 statistic, defined for data y and best-fit model f as R^2=1-\sum{(y_i-f_i)^2/(y_i-y_mean)^2}. In model.py, line1482, the code to calculate rsquared attribute is self.rsquared = 1.0 - (self.residual**2).sum()/max(tiny, sstot), but in your code, residual is not best_fit-data but the return value of the objective function when using the best-fit values of the parameters, which is (best_fit-data)*weights. if weights=None in Model.fit(), the two are same, otherwise they are different. This is a small bug, but I don't know if there's other part where you use residual as best_fit-data. Hope you check your code carefully. Thank you for all your contributions.

A Minimal, Complete, and Verifiable example

import lmfit as lt
import numpy as np

# Model function
def func(x, k=1, b=0):
    return k*x+b

# function to calculate R^2
def fit_R2(modelresult):
    y = modelresult.data
    f = modelresult.best_fit
    ym = sum(y)/len(y)
    return 1 - sum((y-f)**2)/sum((y-ym)**2)

# data
x = np.array([1, 2, 3, 4])
y = np.array([1.1, 1.9, 3.05, 3.95])
yerr = np.array([0.03, 0.04, 0.01, 0.02])

# fitting
mod = lt.Model(func)
params = mod.make_params()
result = mod.fit(y, params, x=x, weights=1/yerr)
print(result.fit_report())

# comparing
print('\n')
print('R^2 in ModelResult:  ', result.rsquared)
print('My R^2:  ', fit_R2(result))
print('best_fit-data:  ', result.best_fit - result.data )
print('(best_fit-data)*weights:  ', (result.best_fit - result.data ) * result.weights)
print('residual:  ', result.residual)

Fit report:

[[Model]]
    Model(func)
[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 7
    # data points      = 4
    # variables        = 2
    chi-square         = 26.0483871
    reduced chi-square = 13.0241935
    Akaike info crit   = 11.4946460
    Bayesian info crit = 10.2672347
    R-squared          = -4.51288616
[[Variables]]
    k:  0.96370968 +/- 0.04150367 (4.31%) (init = 1)
    b:  0.13774194 +/- 0.12714875 (92.31%) (init = 0)
[[Correlations]] (unreported correlations are < 0.100)
    C(k, b) = -0.9713


R^2 in ModelResult:   -4.51288615804741
My R^2:   0.993748167968698
best_fit-data:   [ 0.00145161  0.16516129 -0.02112903  0.04258065]
(best_fit-data)*weights:   [ 0.0483871   4.12903226 -2.11290323  2.12903226]
residual:   [ 0.0483871   4.12903226 -2.11290323  2.12903226]

Error message:

There's no error message, but as you see, ModelResult gives the wrong rsquared.

Version information

Python: 3.12.0 (tags/v3.12.0:0fb18b0, Oct 2 2023, 13:03:39) [MSC v.1935 64 bit (AMD64)]

lmfit: 1.2.2, scipy: 1.11.3, numpy: 1.26.0,asteval: 0.9.31, uncertainties: 3.1.7

[newville]

@ydy1206 Thanks -- yes, it looks like the R^2 value is definitely very wrong when weights are used!
That should be an easy fix, almost certainly

lmfit-py/lmfit/model.py

Line 1482 in c415990

self.rsquared = 1.0 - (self.residual**2).sum()/max(tiny, sstot)

should replace self.residual (which, as you say, is weighted) with dat-self.best_fit (probably moving the calculation of self.best_fit that is just a few lines below.

If you are willing and able to make a Pull Request, that would be great. If not, let us know and we'll fix this very soon.

[ydy1206]

Thank you for reply, but I'm sorry I don't know how to make a Pull Request. Actually, I just register my Github account today, and I'm afraid I'll mess it up if I make a Pull Request.

[newville]

No worries, we'll do this. Thanks!