¿Cómo calculo r cuadrado usando Python y Numpy?

92

Estoy usando Python y Numpy para calcular un polinomio de mejor ajuste de grado arbitrario. Paso una lista de valores x, valores y y el grado del polinomio que quiero ajustar (lineal, cuadrático, etc.).

Esto funciona, pero también quiero calcular r (coeficiente de correlación) y r-cuadrado (coeficiente de determinación). Estoy comparando mis resultados con la capacidad de línea de tendencia de mejor ajuste de Excel y el valor r cuadrado que calcula. Al usar esto, sé que estoy calculando r cuadrado correctamente para el mejor ajuste lineal (el grado es igual a 1). Sin embargo, mi función no funciona para polinomios con grado mayor que 1.

Excel puede hacer esto. ¿Cómo calculo r cuadrado para polinomios de orden superior usando Numpy?

Esta es mi función:

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
Travis Beale
fuente
1
Nota: usa el grado solo en el cálculo de coeficientes.
Nick Dandoulakis
tydok es correcto. Está calculando la correlación de xey y r-cuadrado para y = p_0 + p_1 * x. Vea mi respuesta a continuación para ver un código que debería funcionar. Si no le importa que le pregunte, ¿cuál es su objetivo final? ¿Está haciendo una selección de modelo (eligiendo qué grado usar)? ¿O algo mas?
leif
@leif: la solicitud se reduce a "hazlo como lo hace Excel". De estas respuestas tengo la sensación de que los usuarios pueden estar leyendo demasiado en el valor r cuadrado cuando se usa una curva de mejor ajuste no lineal. No obstante, no soy un mago matemático y esta es la funcionalidad solicitada.
Travis Beale

Respuestas:

62

De la documentación de numpy.polyfit , se ajusta a la regresión lineal. Específicamente, numpy.polyfit con grado 'd' se ajusta a una regresión lineal con la función media

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

Entonces, solo necesita calcular el R cuadrado para ese ajuste. La página de wikipedia sobre regresión lineal ofrece todos los detalles. Está interesado en R ^ 2, que puede calcular de un par de formas, probablemente la más sencilla

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

Donde uso 'y_bar' para la media de las y, y 'y_ihat' para ser el valor de ajuste para cada punto.

No estoy muy familiarizado con numpy (normalmente trabajo en R), por lo que probablemente haya una forma más ordenada de calcular su R cuadrado, pero lo siguiente debería ser correcto

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
leif
fuente
5
Solo quiero señalar que el uso de las funciones de matriz numpy en lugar de la comprensión de la lista será mucho más rápido, por ejemplo, numpy.sum ((yi - ybar) ** 2) y más fácil de leer
Josef
17
Según la página wiki en.wikipedia.org/wiki/Coefficient_of_determination , la definición más general de R ^ 2 es R^2 = 1 - SS_err/SS_tot, R^2 = SS_reg/SS_totsiendo solo un caso especial.
LWZ
137

Una respuesta muy tardía, pero en caso de que alguien necesite una función lista para esto:

scipy.stats.linregress

es decir

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

como en la respuesta de @Adam Marples.

Gökhan Sever
fuente
Es razonable analizar con el coeficiente de correlación y luego hacer el trabajo más grande, la regresión .
象 嘉 道
19
Esta respuesta solo funciona para la regresión lineal, que es la regresión polinomial más simple
tashuhka
9
Precaución: r_value aquí es un coeficiente de correlación de Pearson, no R-cuadrado. r_squared = r_value ** 2
Vladimir Lukin
52

De yanl (otra biblioteca más) sklearn.metricstiene una r2_scorefunción;

from sklearn.metrics import r2_score

coefficient_of_dermination = r2_score(y, p(x))
danodonovan
fuente
1
(Atención: "El valor predeterminado corresponde a 'variance_weighted', este comportamiento está obsoleto desde la versión 0.17 y se cambiará a 'uniform_average' a partir de 0.19")
Franck Dernoncourt
4
r2_score en sklearn podría ser un valor negativo, que no es el caso normal.
Qinqing Liu
1
¿Por qué es r2_score([1,2,3],[4,5,7])= -16?
cz
22

He estado usando esto con éxito, donde xey son como una 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
Adam Marples
fuente
20

Originalmente publiqué los puntos de referencia a continuación con el propósito de recomendar numpy.corrcoef, tontamente sin darme cuenta de que la pregunta original ya se usa corrcoefy, de hecho, preguntaba acerca de los ajustes polinomiales de orden superior. Agregué una solución real a la pregunta del polinomio r-cuadrado usando modelos de estadísticas, y dejé los puntos de referencia originales, que aunque están fuera del tema, son potencialmente útiles para alguien.


statsmodelstiene la capacidad de calcular el r^2ajuste de un polinomio directamente, aquí hay 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 aprovechar aún más statsmodels, también se debe mirar el resumen del modelo ajustado, que puede imprimirse o mostrarse como una tabla HTML enriquecida en el cuaderno Jupyter / IPython. El objeto de resultados proporciona acceso a muchas métricas estadísticas útiles además de rsquared.

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

A continuación se muestra mi respuesta original donde comparé varios métodos de regresión lineal r ^ 2 ...

La función corrcoef utilizada en la Pregunta calcula el coeficiente de correlación r, solo para una regresión lineal simple, por lo que no aborda la cuestión de r^2los ajustes de polinomios de orden superior. Sin embargo, por lo que vale, he llegado a encontrar que para la regresión lineal, es de hecho el método de cálculo más rápido y directo r.

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

Estos fueron los resultados de mi tiempo al comparar varios métodos para 1000 puntos aleatorios (x, y):

  • Python puro ( rcálculo directo )
    • 1000 bucles, lo mejor de 3: 1,59 ms por bucle
  • Numpy polyfit (aplicable a ajustes polinomiales de n-ésimo grado)
    • 1000 bucles, lo mejor de 3: 326 µs por bucle
  • Numpy Manual ( rcálculo directo )
    • 10000 bucles, lo mejor de 3: 62,1 µs por bucle
  • Numpy corrcoef ( rcálculo directo )
    • 10000 bucles, lo mejor de 3: 56,6 µs por bucle
  • Scipy (regresión lineal con rcomo salida)
    • 1000 bucles, lo mejor de 3: 676 µs por bucle
  • Statsmodels (puede hacer polinomios de n-ésimo grado y muchos otros ajustes)
    • 1000 bucles, lo mejor de 3: 422 µs por bucle

El método corrcoef supera por poco el cálculo de r ^ 2 "manualmente" utilizando métodos numerosos. Es> 5 veces más rápido que el método polyfit y ~ 12 veces más rápido que el scipy.linregress. Solo para reforzar lo que Numpy está haciendo por ti, es 28 veces más rápido que Python puro. No estoy muy versado en cosas como numba y pypy, por lo que alguien más tendría que llenar esos vacíos, pero creo que esto me convence bastante de que corrcoefes la mejor herramienta para calcular runa regresión lineal simple.

Aquí está mi código de evaluación comparativa. Copié y pegué de un Jupyter Notebook (es difícil no llamarlo un IPython Notebook ...), así que me disculpo si algo se rompió en el camino. El comando% timeit magic requiere 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)
flutefreak7
fuente
1
Estás comparando 3 métodos ajustando una pendiente y regresión con 3 métodos sin ajustar una pendiente.
Josef
Sí, sabía mucho ... pero ahora me siento tonto por no leer la pregunta original y ver que ya usa corrcoef y se dirige específicamente a r ^ 2 para polinomios de orden superior ... ahora me siento tonto por publicar mis puntos de referencia que eran para un propósito diferente. Ups ...
flutefreak7
1
Actualicé mi respuesta con una solución a la pregunta original usando statsmodelsy me disculpé por la evaluación comparativa innecesaria de los métodos de regresión lineal r ^ 2, que mantuve como información interesante, pero fuera del tema.
flutefreak7
Todavía encuentro interesante el punto de referencia porque no esperaba que el linregress de scipy fuera más lento que los modelos de estadísticas que hacen un trabajo más genérico.
Josef
1
Tenga en cuenta que np.column_stack([x**i for i in range(k+1)])se puede vectorizar en numpy con x[:,None]**np.arange(k+1)o usando las funciones vander de numpy que tienen el orden inverso en las columnas.
Josef
5

R cuadrado es una estadística que solo se aplica a la regresión lineal.

Esencialmente, mide cuánta variación en sus datos se puede explicar mediante la regresión lineal.

Entonces, calcula la "Suma total de cuadrados", que es la desviación al cuadrado total de cada una de las variables de resultado de su media. . .

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

donde y_bar es la media de las y.

Luego, calcula la "suma de cuadrados de regresión", que es cuánto difieren sus valores FITTED de la media.

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

y encuentra la razón de esos dos.

Ahora, todo lo que tendría que hacer para un ajuste polinomial es conectar el y_que es de ese modelo, pero no es exacto llamarlo r cuadrado.

Aquí hay un enlace que encontré que habla un poco.

Baltimark
fuente
Esta parece ser la raíz de mi problema. Entonces, ¿cómo obtiene Excel un valor r cuadrado diferente para un ajuste polinomial frente a una regresión lineal?
Travis Beale
1
¿Solo le está dando a Excel los ajustes de una regresión lineal y los ajustes de un modelo polinomial? Va a calcular el rsq a partir de dos matrices de datos y simplemente asumirá que le está dando los ajustes de un modelo lineal. ¿Qué le estás dando a Excel? ¿Cuál es el comando 'línea de tendencia de mejor ajuste' en Excel?
Baltimark
Es parte de las funciones gráficas de Excel. Puede trazar algunos datos, hacer clic con el botón derecho en ellos y luego elegir entre varios tipos diferentes de líneas de tendencia. Existe la opción de ver la ecuación de la línea, así como un valor de r cuadrado para cada tipo. El valor de r cuadrado también es diferente para cada tipo.
Travis Beale
@Travis Beale: obtendrá un r-cuadrado diferente para cada función media diferente que pruebe (a menos que dos modelos estén anidados y los coeficientes adicionales en el modelo más grande funcionen para ser 0). Entonces, por supuesto, Excel da diferentes valores de r cuadrado. @Baltimark: esta es una regresión lineal, por lo que es r-cuadrado.
leif
5

Aquí hay una función para calcular el r-cuadrado ponderado con Python y Numpy (la mayor parte del código proviene 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

Ejemplo:

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

salidas:

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

Esto corresponde a la fórmula ( espejo ):

ingrese la descripción de la imagen aquí

con f_i es el valor predicho del ajuste, y_ {av} es la media de los datos observados y_i es el valor de los datos observados. w_i es la ponderación aplicada a cada punto de datos, generalmente w_i = 1. SSE es la suma de cuadrados debido al error y SST es la suma total de cuadrados.


Si está interesado, el código en R: https://gist.github.com/dhimmel/588d64a73fa4fef02c8f ( espejo )

Franck Dernoncourt
fuente
2

Aquí hay una función de Python muy simple para calcular R ^ 2 a partir de los valores reales y predichos, asumiendo que y e y_que son series de 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)
Michel Floyd
fuente
0

De la fuente scipy.stats.linregress. Usan el método de suma de cuadrados promedio.

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
Mott la tupla
fuente
0

Puede ejecutar este código directamente, esto le encontrará el polinomio, y le encontrará el valor R, puede poner un comentario abajo si necesita más explicación.

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)
Karam Qusai
fuente