Ciencia de datos aplicada: Su majestad, la Regresión Lineal

Luis Eduardo Pino Villarreal MD, MSc, MBA

Abordaremos hoy el primer modelo que todo científico de datos debe(ría) conocer, me refiero a la regresión lineal (RL).

La RL es uno de los modelos que intenta identificar una relación entre dos o más variables y utilizar dicha relación para predecir valores de una variable con base en los valores dados de otra variable. Como su nombre lo indica, en RL se asume que las relaciones entre variables pueden modelarse a través de la ecuación lineal (y = a + bx) y existe en ella una variable que se usa en la predicción a la cual llamaremos independiente, explicatoria o regresora y una variable predicha a la que llamaremos dependiente, objetivo o de respuesta. En RL entonces asumimos que las variables independientes están relacionadas linealmente con la variable dependiente según la siguiente ecuación:

𝑦 = 𝑐 + 𝑚𝑥

En ML esta ecuación la representamos mejor mediante:

𝑦 = 𝑤0 + 𝑤1𝑥

En donde: w0 es el intercepto en el eje Y, w1 es la pendiente de la línea, x es una variable independiente y y la variable independiente u objetivo.

Por ejemplo, si tenemos las cifras de tensión arterial (TA) de un paciente junto con datos de su edad, podemos usar esta información para predecir los niveles de TA de otro paciente de acuerdo con su edad, este es un problema modelado como de regresión en el cual la edad es una variable independiente y la TA la variable objetivo, entendiendo su relación así:

TA = w0 + w1 x (edad)

En el mundo real no existen problemas con una sola variable independiente y por tanto tendremos que utilizar otras variantes de RL como la regresión múltiple,  pero para fines de explicación básica estos ejemplos funcionan bien.

El mejor gráfico relacional entre estos tipos de variables para el mundo de la RL es el scatterplot que nos permite visualizar la posible fuerza de relación entre dos variables (ver nuestra columna sobre visualización de datos), pero esencialmente nos sirve para intuir si nuestro problema podrá o no resolverse a través de un algoritmo lineal. Es importante entonces conocer cuáles son las medidas que nos permiten establecer la fortaleza de asociación entre variables, es decir las medidas de correlación.

Medidas de correlación y Cálculo de Línea

La más común de todas es el coeficiente de correlación que es una medida para determinar la fuerza de asociación entre dos variables. Dicho coeficiente va de -1 (fuerte asociación negativa) a 1 (fuerte asociación positiva), pasando por el cero (0) que indica ausencia de relación entre variables. En ML las matrices de correlación nos pueden dar una idea inicial de este coeficiente:

Ejemplo

En esta matriz podemos ver correlaciones positivas altas como la de las variables Power y Engine y otras negativas como Engine y Mileage (ejemplo tomado de un dataset de carros usados disponible en Kaggle).

Ahora bien, si tenemos un conjunto de datos tabulares con dos variables posiblemente correlacionadas, por ejemplo:

Horas de estudio por díaCalificación final (/10)
67
55
42
79
86
23

Tendremos la siguiente gráfica en un plano cartesiano (en este caso la variable dependiente es la calificación final (X) y la variable independiente es horas de estudio por día (Y))

Si observamos la gráfica podríamos pensar en trazar una línea que uniera estas variables estableciendo su relación lineal, la pregunta es cómo puede calcularse la “mejor línea” para cubrir dicha correlación. Regresemos entonces a nuestra ecuación base:

y= w0 + w1 * x

Recordemos que W0 es el intercepto y X1 es la pendiente. Si trazamos una línea cualquiera es posible que algunos de los puntos del plano cartesiano NO sean cubiertos por ella y excepto que la línea pase directamente por ese punto, habra siempre un error o residual que será la “distancia” entre el punto y la línea, este puede medirse así:

La predicción (yˆi) para cualquier x (xi en este caso) es:

yˆi = w0 + w1*xi

Entonces el valor de los residuales (e) para cualquier punto i es la diferencia:

ei = yi – yˆi

Teniendo en cuenta lo que es un residual tenemos entonces que la mejor línea será aquella que minimice la suma de los errores cuadrados (es decir los ei2)

Ei2 = ( y1 – yˆ1)2 + ( y2 – yˆ2)2 + (yn – yˆn)2

Los cálculos secuenciales para encontrar esta línea se denominan relaciones lineasles cuadráticas y pueden verse así:

Los dos métodos mas utilizados para calcular esta línea son el de mínimos cuadrados y el de gradiente descendente (atención a este último que será muy importante en modelos posteriores de ML).

Dado que no es objeto de la columna profundizar en estos métodos remito al lector a las referencias recomendadas al final de nuestra serie. Debo mencionar que los métodos a ejecutar difieren en su contexto matemático y en sus medidas de desempeño, en el caso de los mínimos cuadrados aplicaremos el R2 y el R2 ajustado mientras en el gradiente descendente utilizaremos funciones de costo. Los dos modelos incluyen cálculos de derivadas y por tanto es importante antes de profundizar en ellos desempolvar un libro de cálculo básico.

Cómo mencionamos previamente, en medicina no existen problemas que se ajusten a una regresión lineal simple, usualmente tenemos múltiples variables que influencia un desenlace. En este caso nos asistimos de los modelos de regresión lineal múltiple cuya ecuación es:

y= w0 + w1*x1 + w2*x2 + w3*x3…+wn*xn

En este caso y es la variable dependiente y x las diversas variables dependientes. Los métodos para calcular los coeficientes y la mejor línea son los descritos previamente.

Medidas de desempeño

Un elemento crucial en el desarrollo de los modelos de ML son las métricas de desempeño. Para el caso de las regresiones las medidas de desempeño que se utilizan son:

  1. R2 y R2 ajustado

El R2 es esencialmente la diferencia en la propagación o el “spread” entre los residuales iniciales y los encontrados post regresión. En ML el R2 es una medida interna de la muestra (in sample).

Puede expresarse como

R2 = 1 – ∑ ei2 / ∑ (y1 – yˆ)2

En general R2 expresa la fracción de la varianza en y que es explicada por la regresión y se expresa de 0 a 1, en donde 1 expresa una correlación máxima.

El R2 ajustado lo utilizamos especialmente en regresiones múltiples dado que su valor se aumenta con cada variable añadida.

  • Error absoluto promedio (MAE)

Es la media de la diferencia entre los variables actuales (y) y predichos (yˆ) de la variable dependiente.

  • Raiz cuadrada del error cuadrático promedio (RMCE)

Esta ecuación cuadrática se utiliza para estimar la desviación estandar de los residuales y se usa la potencia para disminuir el sesgo que produce el n-1, la ecuación se expresa:

Ventajas y Limitaciones de la Regresión Lineal

VentajasLimitaciones
– Es un modelo de amplio uso y muy antiguo .
– Simple y elegante  
– Computacionalmente muy eficiente.  
– Facil interpretación de sus salidas (outputs)
– Muy simple para la complejidad del mundo real
– Asume relaciones lineales que son infrecuentes  
– Los outliers pueden afectar su desempeño  
– Asume independencia entre atributos, es dificl que pueda gestionar la multicolinealidad

Acompañanos a ver un ejemplo de regresión lineal (Anexo 1): file:///C:/Users/57310/Desktop/Ejemplo_RegLineal.html

En nuestra próxima columna hablaremos sobre regresión logística, aquella en la cual la variable dependiente no es un número sino una o varias categorías y por tanto difiere de la regresión lineal simple y múltiple en algunas características y métricas.

Anexo 1. Ejemplo de regresión líneal.

Regresión Lineal by AIpocrates

A continuación un ejemplo básico de Regresión Lineal en una base de datos de carros

Primero importamos las librerías con las que trabajaremos

In [1]:

import numpy as np   # importamos Numpy para funciones matemáticas 
from sklearn.linear_model import LinearRegression # importanos el modelo de regresión desde Sklearn
import pandas as pd    #Importamos Pandas para manejo de dataframes y math
import matplotlib.pyplot as plt # importamos nuestro visor de gráficas
%matplotlib inline # Este código nos permite trabajar las gráficas en el mismo archivo
import seaborn as sns # nos sirve como apoyo en la visualización de datos estadísticos
from sklearn.model_selection import train_test_split # Libreria para fraccionar los datos en entrenamiento y prueba
C:\Users\Rajet singh\Anaconda3\lib\site-packages\statsmodels\tools\_testing.py:19: FutureWarning: pandas.util.testing is deprecated. Use the functions in the public API at pandas.testing instead.
  import pandas.util.testing as tm

Cargando y revisando la base de datos

In [2]:

cData = pd.read_csv("auto-mpg.csv")  # importamos la base de datos en formato CSV
cData.shape # vemos las características de la base (en este caso tiene 398 registros y 9 variables)

Out[2]:

(398, 9)

In [3]:

# Las 9 variables: 
# Numero
# MPG (miles per gallon), 
# cylinders, 
# engine displacement (cu. inches), 
# horsepower,
# vehicle weight (lbs.), 
# time to accelerate from O to 60 mph (sec.),
# model year (modulo 100), and 
# origin of car (1. American, 2. European,3. Japanese).

# veamos las primeras cinco columnas con este código:
cData.head()

Out[3]:

mpgcylindersdisplacementhorsepowerweightaccelerationmodel yearorigincar name
018.08307.0130350412.0701chevrolet chevelle malibu
115.08350.0165369311.5701buick skylark 320
218.08318.0150343611.0701plymouth satellite
316.08304.0150343312.0701amc rebel sst
417.08302.0140344910.5701ford torino
La idea de este ejercicio será evaluar qué variable o variables se relacionan con el consumo de millas por galón (mpg) de este grupo de carros.

In [4]:

# vamos a quitar la variable car name ya que no contribuye al ejercicio
cData = cData.drop('car name', axis=1)
# Reemplazemos la variable origin con los nombres respectivos
cData['origin'] = cData['origin'].replace({1: 'america', 2: 'europe', 3: 'asia'})
cData.head()

Out[4]:

mpgcylindersdisplacementhorsepowerweightaccelerationmodel yearorigin
018.08307.0130350412.070america
115.08350.0165369311.570america
218.08318.0150343611.070america
316.08304.0150343312.070america
417.08302.0140344910.570america

Creación de variables Dummy

Creamos variables binarias para facilitar el análisis de los datos

In [5]:

cData = pd.get_dummies(cData, columns=['origin'])
cData.head()

Out[5]:

mpgcylindersdisplacementhorsepowerweightaccelerationmodel yearorigin_americaorigin_asiaorigin_europe
018.08307.0130350412.070100
115.08350.0165369311.570100
218.08318.0150343611.070100
316.08304.0150343312.070100
417.08302.0140344910.570100

Tratamiento de los datos perdidos

In [6]:

#Un breve resumen de las columnas:
cData.describe()

Out[6]:

mpgcylindersdisplacementweightaccelerationmodel yearorigin_americaorigin_asiaorigin_europe
count398.000000398.000000398.000000398.000000398.000000398.000000398.000000398.000000398.000000
mean23.5145735.454774193.4258792970.42462315.56809076.0100500.6256280.1984920.175879
std7.8159841.701004104.269838846.8417742.7576893.6976270.4845690.3993670.381197
min9.0000003.00000068.0000001613.0000008.00000070.0000000.0000000.0000000.000000
25%17.5000004.000000104.2500002223.75000013.82500073.0000000.0000000.0000000.000000
50%23.0000004.000000148.5000002803.50000015.50000076.0000001.0000000.0000000.000000
75%29.0000008.000000262.0000003608.00000017.17500079.0000001.0000000.0000000.000000
max46.6000008.000000455.0000005140.00000024.80000082.0000001.0000001.0000001.000000

In [9]:

# Reemplazamos los valores perdidos con NaN (not a number)
cData = cData.replace('?', np.nan)
cData[hpIsDigit['horsepower'] == False] 

Out[9]:

mpgcylindersdisplacementhorsepowerweightaccelerationmodel yearorigin_americaorigin_asiaorigin_europe
3225.0498.0NaN204619.071100
12621.06200.0NaN287517.074100
33040.9485.0NaN183517.380001
33623.64140.0NaN290514.380100
35434.54100.0NaN232015.881001
37423.04151.0NaN303520.582100

In [10]:

# Los datos perdidos los vamos a tratar reemplazandolos por una medida de tendencia central (mediana) que es para cada variable:
cData.median()

Out[10]:

mpg                 23.0
cylinders            4.0
displacement       148.5
horsepower          93.5
weight            2803.5
acceleration        15.5
model year          76.0
origin_america       1.0
origin_asia          0.0
origin_europe        0.0
dtype: float64

In [11]:

# Para reemplazarlos creamos una función lambda:

medianFiller = lambda x: x.fillna(x.median())
cData = cData.apply(medianFiller,axis=0)

cData['horsepower'] = cData['horsepower'].astype('float64')  # converting the hp column from object / string type to float

Análisis Bivariado

Vamos a ver el análisis gráfico bivariado usando pairplot

In [12]:

cData_attr = cData.iloc[:, 0:7]
sns.pairplot(cData_attr, diag_kind='kde')   # Esto nos permite ver densidades y las posibles relaciones entre dos variables

Out[12]:

<seaborn.axisgrid.PairGrid at 0x2aba2117108>

Nuestra variable objetivo es mpg la cual parece tener algun grado de correlación con horsepower y weight. Vamos a construir el modelo para verificarlo.

Hacemos split de los datos antes de ajustar el modelo de RL

In [13]:

## Primero las independientes en donde retiramos a mpg y a origin_europe (esta no la queremos en el modelo)
X = cData.drop(['mpg','origin_europe'], axis=1)
# y en Y nuestra variable objetivo o dependiente
y = cData[['mpg']]

In [14]:

# Rompemos los datos en sets de entrenamiento y de prueba en relación 70:30 con estado random 1 para mantener esta estructura de azar en las iteraciones

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.30, random_state=1)

Ajustando el modelo de RL

In [15]:

# elegimos el modelo LR desde la sklearn y lo ajustaremos por el intercepto (mejor linea) sin normalizacion
regression_model = LinearRegression()
regression_model.fit(X_train, y_train)

Out[15]:

LinearRegression(copy_X=True, fit_intercept=True, n_jobs=None, normalize=False)

Ahora calcularemos los coeficientes de regresion de acuerdo con el intercepto:

In [16]:

for idx, col_name in enumerate(X_train.columns):
    print("The coefficient for {} is {}".format(col_name, regression_model.coef_[0][idx]))
The coefficient for cylinders is -0.3948079661648244
The coefficient for displacement is 0.028945510765487174
The coefficient for horsepower is -0.02175220772354676
The coefficient for weight is -0.007352032065147351
The coefficient for acceleration is 0.061919366007618326
The coefficient for model year is 0.8369338917644993
The coefficient for origin_america is -3.0012830009185123
The coefficient for origin_asia is -0.6060179643247369
El año del modelo tiene una fuerte correlación positiva con mpg, La mayoria de variables restantes tienen una correlación negativa siendo el origen en América la de mayor correlación negativa. Ahora calcularemos el intercepto de todo el modelo

In [17]:

intercept = regression_model.intercept_[0]
print("The intercept for our model is {}".format(intercept))
The intercept for our model is -18.283451116372067

Ahora calcularemos una métrica de desempeño que es el R2:

In [18]:

regression_model.score(X_train, y_train)

Out[18]:

0.8141025501610559

el R2 es de 0.81 que es elevado lo cual quiere decir que el modelo acertará en la correlación 82% de las veces, ahora calcularemos este R2 por fuera de la muestra:

In [19]:

#out of sample score (R^2)

regression_model.score(X_test, y_test)

Out[19]:

0.8433135132808833

Al ser positivo el R2 fuera de muestra nuestro modelo tiene menor error cuadrático promedio lo que quiere decir que es reproducible.

Adicionemos términos de interacción

Esto nos permitirá hacer ajustes polinomiales para optimizar variables y mejorar (posiblemente) el rendimiento del modelo

In [20]:

from sklearn.preprocessing import PolynomialFeatures 
from sklearn import linear_model

poly = PolynomialFeatures(degree=2, interaction_only=True)
X_train2 = poly.fit_transform(X_train)
X_test2 = poly.fit_transform(X_test)

poly_clf = linear_model.LinearRegression()

poly_clf.fit(X_train2, y_train)

y_pred = poly_clf.predict(X_test2)

# veamos ahora nuestro nuevo R2 en muestra
print(poly_clf.score(X_train2, y_train))
0.9015975294607972

In [21]:

# y el R2 fuera de muestra
print(poly_clf.score(X_test2, y_test))
0.8647441061208313

Nuestros R2 mejoraron ampliamente, pero esto ha sido a costa de aumentar variables (de 8 a 37):

In [22]:

print(X_train.shape)
print(X_train2.shape)
(278, 8)
(278, 37)

Este ha sido un pequeño ejemplo de como puede trabajarse un modelo de RL en python versión 3.3. Gracias por su lectura

Luis Eduardo Pino Villarreal MD, MSc. MBA
Fundador de AIpocrates
CEO OxLER

luispino@oxler.me http://www.aipocrates.org http://www.oxler.info

Deja una respuesta

Introduce tus datos o haz clic en un icono para iniciar sesión:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Salir /  Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Salir /  Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Salir /  Cambiar )

Conectando a %s

A %d blogueros les gusta esto: