Como faço para calcular r-quadrado usando Python e Numpy?

Estou usando Python e Numpy para calcular um polinômio de melhor ajuste de grau arbitrário. Eu passo uma lista de valores x, valores y e o grau do polinômio que desejo ajustar (linear, quadrático, etc.).

Isso funciona, mas também quero calcular r (coeficiente de correlação) e r-quadrado (coeficiente de determinação). Estou comparando meus resultados com a capacidade de linha de tendência de melhor ajuste do Excel e o valor de r quadrado que ele calcula. Usando isso, sei que estou calculando r-quadrado corretamente para o melhor ajuste linear (grau igual a 1). No entanto, minha função não funciona para polinômios com grau maior que 1.

O Excel é capaz de fazer isso. Como calculo r-quadrado para polinômios de ordem superior usando Numpy?

Esta é minha função:

import numpy

# Polynomial Regression
def polyfit(x, y, degree):
    results = {}

    coeffs = numpy.polyfit(x, y, degree)
     # Polynomial Coefficients
    results['polynomial'] = coeffs.tolist()

    correlation = numpy.corrcoef(x, y)[0,1]

     # r
    results['correlation'] = correlation
     # r-squared
    results['determination'] = correlation**2

    return results

A partir da documentação numpy.polyfit , é uma regressão linear adequada. Especificamente, numpy.polyfit com o grau 'd' ajusta uma regressão linear com a função média

E (y | x) = p_d * x ** d + p_ {d-1} * x ** (d-1) + ... + p_1 * x + p_0

Portanto, você só precisa calcular o R-quadrado desse ajuste. A página da Wikipedia sobre regressão linear fornece detalhes completos. Você está interessado em R ^ 2, que você pode calcular de algumas maneiras, provavelmente a mais fácil

SST = Sum(i=1..n) (y_i - y_bar)^2
SSReg = Sum(i=1..n) (y_ihat - y_bar)^2
Rsquared = SSReg/SST

Onde eu uso 'y_bar' para a média dos y, e 'y_ihat' para ser o valor de ajuste para cada ponto.

Não estou muito familiarizado com o numpy (geralmente trabalho com R), então provavelmente há uma maneira mais organizada de calcular seu R ao quadrado, mas o seguinte deve estar correto

import numpy

# Polynomial Regression
def polyfit(x, y, degree):
    results = {}

    coeffs = numpy.polyfit(x, y, degree)

     # Polynomial Coefficients
    results['polynomial'] = coeffs.tolist()

    # r-squared
    p = numpy.poly1d(coeffs)
    # fit values, and mean
    yhat = p(x)                         # or [p(z) for z in x]
    ybar = numpy.sum(y)/len(y)          # or sum(y)/len(y)
    ssreg = numpy.sum((yhat-ybar)**2)   # or sum([ (yihat - ybar)**2 for yihat in yhat])
    sstot = numpy.sum((y - ybar)**2)    # or sum([ (yi - ybar)**2 for yi in y])
    results['determination'] = ssreg / sstot

    return results

Uma resposta muito tardia, mas apenas no caso de alguém precisar de uma função pronta para isso:

scipy.stats.linregress

ie

slope, intercept, r_value, p_value, std_err = scipy.stats.linregress(x, y)

como na resposta de @Adam Marples.

From yanl (yet-another-library) sklearn.metricstem uma r2_scorefunção;

from sklearn.metrics import r2_score

coefficient_of_dermination = r2_score(y, p(x))

Tenho usado isso com sucesso, em que xey são semelhantes a uma matriz.

def rsquared(x, y):
    """ Return R^2 where x and y are array-like."""

    slope, intercept, r_value, p_value, std_err = scipy.stats.linregress(x, y)
    return r_value**2

Eu originalmente postei os benchmarks abaixo com o propósito de recomendar numpy.corrcoef, tolamente não percebendo que a pergunta original já usa corrcoefe, na verdade, estava perguntando sobre ajustes polinomiais de ordem superior. Eu adicionei uma solução real para a questão polinomial r-quadrado usando modelos de estatísticas e deixei os benchmarks originais, que embora fora do tópico, são potencialmente úteis para alguém.

statsmodelstem a capacidade de calcular o r^2de um ajuste polinomial diretamente, aqui estão 2 métodos ...

import statsmodels.api as sm
import statsmodels.formula.api as smf

# Construct the columns for the different powers of x
def get_r2_statsmodels(x, y, k=1):
    xpoly = np.column_stack([x**i for i in range(k+1)])    
    return sm.OLS(y, xpoly).fit().rsquared

# Use the formula API and construct a formula describing the polynomial
def get_r2_statsmodels_formula(x, y, k=1):
    formula = 'y ~ 1 + ' + ' + '.join('I(x**{})'.format(i) for i in range(1, k+1))
    data = {'x': x, 'y': y}
    return smf.ols(formula, data).fit().rsquared # or rsquared_adj

Para aproveitar ainda mais statsmodels, deve-se também olhar o resumo do modelo ajustado, que pode ser impresso ou exibido como uma rica tabela HTML no bloco de notas Jupyter / IPython. O objeto de resultados fornece acesso a muitas métricas estatísticas úteis, além de rsquared.

model = sm.OLS(y, xpoly)
results = model.fit()
results.summary()

Abaixo está minha resposta original, onde comparei vários métodos de regressão linear r ^ 2 ...

A função corrcoef usada na pergunta calcula o coeficiente de correlação ,,r apenas para uma única regressão linear, portanto, não aborda a questão de r^2para ajustes polinomiais de ordem superior. No entanto, pelo que vale a pena, descobri que, para a regressão linear, é de fato o método de cálculo mais rápido e direto r.

def get_r2_numpy_corrcoef(x, y):
    return np.corrcoef(x, y)[0, 1]**2

Estes foram os meus resultados timeit comparando um monte de métodos para 1000 pontos aleatórios (x, y):

• Python puro ( rcálculo direto )
  • 1000 loops, melhor de 3: 1,59 ms por loop
• Numpy polyfit (aplicável a ajustes polinomiais de n-ésimo grau)
  • 1000 loops, melhor de 3: 326 µs por loop
• Numpy Manual ( rcálculo direto )
  • 10.000 loops, melhor de 3: 62,1 µs por loop
• Numpy corrcoef ( rcálculo direto )
  • 10.000 loops, melhor de 3: 56,6 µs por loop
• Scipy (regressão linear com rsaída)
  • 1000 loops, melhor de 3: 676 µs por loop
• Modelos de estatísticas (podem fazer polinômios de n-ésimo grau e muitos outros ajustes)
  • 1000 loops, melhor de 3: 422 µs por loop

O método corrcoef supera o cálculo de r ^ 2 "manualmente" usando métodos numpy. É> 5X mais rápido do que o método polyfit e ~ 12X mais rápido do que o scipy.linregress. Apenas para reforçar o que o numpy está fazendo por você, ele é 28 vezes mais rápido do que o python puro. Não sou muito versado em coisas como numba e pypy, então outra pessoa teria que preencher essas lacunas, mas acho que isso é muito convincente para mim que corrcoefé a melhor ferramenta para calcular ruma regressão linear simples.

Aqui está o meu código de benchmarking. Copiei e colei de um Notebook Jupyter (difícil não chamá-lo de Notebook IPython ...), então peço desculpas se alguma coisa quebrou no caminho. O comando mágico% timeit requer IPython.

import numpy as np
from scipy import stats
import statsmodels.api as sm
import math

n=1000
x = np.random.rand(1000)*10
x.sort()
y = 10 * x + (5+np.random.randn(1000)*10-5)

x_list = list(x)
y_list = list(y)

def get_r2_numpy(x, y):
    slope, intercept = np.polyfit(x, y, 1)
    r_squared = 1 - (sum((y - (slope * x + intercept))**2) / ((len(y) - 1) * np.var(y, ddof=1)))
    return r_squared
    
def get_r2_scipy(x, y):
    _, _, r_value, _, _ = stats.linregress(x, y)
    return r_value**2
    
def get_r2_statsmodels(x, y):
    return sm.OLS(y, sm.add_constant(x)).fit().rsquared
    
def get_r2_python(x_list, y_list):
    n = len(x_list)
    x_bar = sum(x_list)/n
    y_bar = sum(y_list)/n
    x_std = math.sqrt(sum([(xi-x_bar)**2 for xi in x_list])/(n-1))
    y_std = math.sqrt(sum([(yi-y_bar)**2 for yi in y_list])/(n-1))
    zx = [(xi-x_bar)/x_std for xi in x_list]
    zy = [(yi-y_bar)/y_std for yi in y_list]
    r = sum(zxi*zyi for zxi, zyi in zip(zx, zy))/(n-1)
    return r**2
    
def get_r2_numpy_manual(x, y):
    zx = (x-np.mean(x))/np.std(x, ddof=1)
    zy = (y-np.mean(y))/np.std(y, ddof=1)
    r = np.sum(zx*zy)/(len(x)-1)
    return r**2
    
def get_r2_numpy_corrcoef(x, y):
    return np.corrcoef(x, y)[0, 1]**2
    
print('Python')
%timeit get_r2_python(x_list, y_list)
print('Numpy polyfit')
%timeit get_r2_numpy(x, y)
print('Numpy Manual')
%timeit get_r2_numpy_manual(x, y)
print('Numpy corrcoef')
%timeit get_r2_numpy_corrcoef(x, y)
print('Scipy')
%timeit get_r2_scipy(x, y)
print('Statsmodels')
%timeit get_r2_statsmodels(x, y)

R-quadrado é uma estatística que se aplica apenas à regressão linear.

Essencialmente, ele mede quanta variação em seus dados pode ser explicada pela regressão linear.

Assim, você calcula a "Soma Total dos Quadrados", que é o desvio quadrático total de cada uma das variáveis ​​de resultado de sua média. . .

\ sum_ {i} (y_ {i} - y_bar) ^ 2

onde y_bar é a média dos y's.

Então, você calcula a "soma dos quadrados da regressão", que é o quanto seus valores FITTED diferem da média

\ sum_ {i} (yHat_ {i} - y_bar) ^ 2

e encontre a proporção desses dois.

Agora, tudo que você teria que fazer para um ajuste polinomial é conectar o y_hat desse modelo, mas não é preciso chamar isso de r ao quadrado.

Aqui está um link que encontrei que fala um pouco sobre isso.

O artigo da Wikipedia sobre r-quadrados sugere que ele pode ser usado para ajuste de modelo geral, em vez de apenas regressão linear.

Aqui está uma função para calcular o r-quadrado ponderado com Python e Numpy (a maior parte do código vem de sklearn):

from __future__ import division 
import numpy as np

def compute_r2_weighted(y_true, y_pred, weight):
    sse = (weight * (y_true - y_pred) ** 2).sum(axis=0, dtype=np.float64)
    tse = (weight * (y_true - np.average(
        y_true, axis=0, weights=weight)) ** 2).sum(axis=0, dtype=np.float64)
    r2_score = 1 - (sse / tse)
    return r2_score, sse, tse

Exemplo:

from __future__ import print_function, division 
import sklearn.metrics 

def compute_r2_weighted(y_true, y_pred, weight):
    sse = (weight * (y_true - y_pred) ** 2).sum(axis=0, dtype=np.float64)
    tse = (weight * (y_true - np.average(
        y_true, axis=0, weights=weight)) ** 2).sum(axis=0, dtype=np.float64)
    r2_score = 1 - (sse / tse)
    return r2_score, sse, tse    

def compute_r2(y_true, y_predicted):
    sse = sum((y_true - y_predicted)**2)
    tse = (len(y_true) - 1) * np.var(y_true, ddof=1)
    r2_score = 1 - (sse / tse)
    return r2_score, sse, tse

def main():
    '''
    Demonstrate the use of compute_r2_weighted() and checks the results against sklearn
    '''        
    y_true = [3, -0.5, 2, 7]
    y_pred = [2.5, 0.0, 2, 8]
    weight = [1, 5, 1, 2]
    r2_score = sklearn.metrics.r2_score(y_true, y_pred)
    print('r2_score: {0}'.format(r2_score))  
    r2_score,_,_ = compute_r2(np.array(y_true), np.array(y_pred))
    print('r2_score: {0}'.format(r2_score))
    r2_score = sklearn.metrics.r2_score(y_true, y_pred,weight)
    print('r2_score weighted: {0}'.format(r2_score))
    r2_score,_,_ = compute_r2_weighted(np.array(y_true), np.array(y_pred), np.array(weight))
    print('r2_score weighted: {0}'.format(r2_score))

if __name__ == "__main__":
    main()
    #cProfile.run('main()') # if you want to do some profiling

saídas:

r2_score: 0.9486081370449679
r2_score: 0.9486081370449679
r2_score weighted: 0.9573170731707317
r2_score weighted: 0.9573170731707317

Isso corresponde à fórmula ( espelho ):

com f_i é o valor previsto do ajuste, y_ {av} é a média dos dados observados y_i é o valor dos dados observados. w_i é a ponderação aplicada a cada ponto de dados, geralmente w_i = 1. SSE é a soma dos quadrados devido ao erro e SST é a soma total dos quadrados.

Se estiver interessado, o código em R: https://gist.github.com/dhimmel/588d64a73fa4fef02c8f ( espelho )

Aqui está uma função python muito simples para calcular R ^ 2 a partir dos valores reais e previstos, assumindo y e y_hat são séries pandas:

def r_squared(y, y_hat):
    y_bar = y.mean()
    ss_tot = ((y-y_bar)**2).sum()
    ss_res = ((y-y_hat)**2).sum()
    return 1 - (ss_res/ss_tot)

Da fonte scipy.stats.linregress. Eles usam o método da soma média dos quadrados.

import numpy as np

x = np.array(x)
y = np.array(y)

# average sum of squares:
ssxm, ssxym, ssyxm, ssym = np.cov(x, y, bias=1).flat

r_num = ssxym
r_den = np.sqrt(ssxm * ssym)
r = r_num / r_den

if r_den == 0.0:
    r = 0.0
else:
    r = r_num / r_den

    if r > 1.0:
        r = 1.0
    elif r < -1.0:
        r = -1.0

Você pode executar este código diretamente, ele encontrará o polinômio e o valor R. Você pode colocar um comentário abaixo se precisar de mais explicações.

from scipy.stats import linregress
import numpy as np

x = np.array([1,2,3,4,5,6])
y = np.array([2,3,5,6,7,8])

p3 = np.polyfit(x,y,3) # 3rd degree polynomial, you can change it to any degree you want
xp = np.linspace(1,6,6)  # 6 means the length of the line
poly_arr = np.polyval(p3,xp)

poly_list = [round(num, 3) for num in list(poly_arr)]
slope, intercept, r_value, p_value, std_err = linregress(x, poly_list)
print(r_value**2)