Pasos para realizar una regresión lineal múltiple con Python
1). Importar los siguientes módulos o paquetes:
Éstos son necesarios para realizar la regresión lineal múltiple.
In [1]: import pandas as pd ## Este proporciona una estructura similiar a los data.
frame
import [Link] as sm ## Este proporciona funciones para la estimación de muc
hos modelos estadísticos
import [Link] as smf ## Permite ajustar modelos estadísticos utilizando fórm
ulas de estilo R
2). Leer la base de datos:
Para esto, se usará la función read_csv del paquete pandas.
In [2]: file = '[Link]
ion/Regresi%C3%B3n%20lineal%20m%C3%BAltiple/[Link]'
dt = pd.read_csv(file)
3). Estructura y primeras observaciones de la base de datos:
En este paso, se debe observar el conjunto de datos importados, con el fin de familiarizarse con las características esenciales de
los datos.
In [3]: [Link]() # Sirve para mostrar la parte inicial de los datos
# [Link](15) muestra las 15 primeras observaciones
# (Existen más funciones importantes y complementarias para ésto)
Out[3]:
Obs y x1 x2
0 1 16.68 7 560
1 2 11.50 3 220
2 3 12.03 3 340
3 4 14.88 4 80
4 5 13.75 6 150
4). Diagrama de dispersión:
Éste gráfico ayuda a analizar la forma en la que variables están relacionadas. (Es opcional, pero útil).
In [5]: from matplotlib import pyplot # Permite la generación de gráficos
from mpl_toolkits.mplot3d import Axes3D # Permite agregar eje tridimensionales
import random # Permiten obtener de distintos modos números aleato
rios
fig = [Link](figsize=(8, 6)) # Ajustes del gráfico
ax = Axes3D(fig)
x1 = dt["x1"] # Datos eje X
x2 = dt["x2"] # Datos eje Y
y = dt["y"] # Datos eje Z (Var. Respuesta)
[Link](x1, x2, y, marker='*', c='r')
ax.set_xlabel('Cantidad de cajas') # Etiqueta del eje X
ax.set_ylabel('Distancia (metros)') # Etiqueta del eje Y
ax.set_zlabel('Tiempo (minutos)'); # Etiqueta del eje Z (Var. Respuesta)
5). Ajustar el modelo de regresión mediante fórmulas:
Para ésto, se usa [Link].
In [6]: mod = [Link]('y ~ x1 + x2', data=dt).fit() # Ajusta el modelo usando el registro natural de
uno de los regresores
6). Resumen del modelo ajustado:
In [7]: print([Link]())
OLS Regression Results
==============================================================================
Dep. Variable: y R-squared: 0.960
Model: OLS Adj. R-squared: 0.956
Method: Least Squares F-statistic: 261.2
Date: Sun, 02 Jun 2019 Prob (F-statistic): 4.69e-16
Time: 15:38:16 Log-Likelihood: -63.415
No. Observations: 25 AIC: 132.8
Df Residuals: 22 BIC: 136.5
Df Model: 2
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
Intercept 2.3412 1.097 2.135 0.044 0.067 4.616
x1 1.6159 0.171 9.464 0.000 1.262 1.970
x2 0.0144 0.004 3.981 0.001 0.007 0.022
==============================================================================
Omnibus: 0.421 Durbin-Watson: 1.170
Prob(Omnibus): 0.810 Jarque-Bera (JB): 0.010
Skew: 0.032 Prob(JB): 0.995
Kurtosis: 3.073 Cond. No. 873.
==============================================================================
Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
7). Extraer los β estimados:
In [8]: [Link]
Out[8]: Intercept 2.341231
x1 1.615907
x2 0.014385
dtype: float64
8). Extraer el σ 2 estimado:
Para esto se usa el atributo [Link] sobre el objeto mod.
In [9]: mod.mse_resid
Out[9]: 10.624167155479673
8). Agregar el plano de regresión:
In [10]: from mpl_toolkits.mplot3d import Axes3D
import [Link] as plt
import numpy as np
import pandas as pd
import [Link] as sm
from matplotlib import cm
csv = pd.read_csv('[Link] index_col=0)
model = [Link](formula='sales ~ TV + radio', data = csv)
fit = [Link]()
[Link]()
xx, yy = [Link](csv['TV'], csv['radio'])
fig = [Link]()
ax = fig.add_subplot(111, projection='3d')
x_surf = [Link](0, 350, 20) # generar una malla
y_surf = [Link](0, 60, 4)
x_surf, y_surf = [Link](x_surf, y_surf)
exog = [Link]({'TV': x_surf.ravel(), 'radio': y_surf.ravel()})
out = [Link](exog = exog)
exog = [Link]({'TV':[Link](),'radio':[Link]()})
out = [Link](exog=exog)
#ax.plot_surface(xx, yy, [Link]([Link]), color='None')
print([Link])
[Link](csv['TV'], csv['radio'], csv['sales'],
c='blue',
marker='o',
alpha=1)
ax.set_xlabel('TV')
ax.set_ylabel('Radio')
ax.set_zlabel('Sales')
[Link]()
(200, 200)
In [12]: import [Link] as plt
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
from sklearn import datasets, linear_model
diabetes = datasets.load_diabetes()
indices = (0, 1)
X_train = [Link][:-20, indices]
X_test = [Link][-20:, indices]
y_train = [Link][:-20]
y_test = [Link][-20:]
ols = linear_model.LinearRegression()
[Link](X_train, y_train)
# #############################################################################
# Traza la figura
def plot_figs(fig_num, elev, azim, X_train, clf):
fig = [Link](fig_num, figsize=(4, 3))
[Link]()
ax = Axes3D(fig, elev=elev, azim=azim)
[Link](X_train[:, 0], X_train[:, 1], y_train, c='k', marker='+')
ax.plot_surface([Link]([[-.1, -.1], [.15, .15]]),
[Link]([[-.1, .15], [-.1, .15]]),
[Link]([Link]([[-.1, -.1, .15, .15],
[-.1, .15, -.1, .15]]).T
).reshape((2, 2)),
alpha=.5)
ax.set_xlabel('X_1')
ax.set_ylabel('X_2')
ax.set_zlabel('Y')
ax.w_xaxis.set_ticklabels([])
ax.w_yaxis.set_ticklabels([])
ax.w_zaxis.set_ticklabels([])
# Generar las tres diferentes figuras desde distintas perspectivas
elev = 43.5
azim = -110
plot_figs(1, elev, azim, X_train, ols)
elev = -.5
azim = 0
plot_figs(2, elev, azim, X_train, ols)
elev = -.5
azim = 90
plot_figs(3, elev, azim, X_train, ols)
[Link]()
Fuente: la base de datos utilizada son del ejemplo 2.1 del libro de Montgomery, Peck and Vining (2003). Además, dicha base se
encuentra alojada en el Github del profesor Freddy Hernández Barajas, [Link] ([Link]
Documento creado por Yuberth Anderson Saavedra Coneo, durante el curso Análisis de Regresión de la Universidad
Nacional de Colombia sede Medellín.