Linear models for classification¶
In [1]:
from scipy import interpolate
import scipy as sp
import sklearn
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LogisticRegression
from sklearn.inspection import DecisionBoundaryDisplay
import warnings
warnings.simplefilter("ignore")
In [2]:
data_heart=pd.read_csv('data/heart_deasease.csv')
X_heart=data_heart.copy()
X_heart.drop(['chd','row.names','famhist', 'obesity'],axis=1,inplace=True)
y_heart=data_heart.chd
data_vowel=pd.read_csv('data/vowel.csv')
X_vowel=data_vowel.copy()
X_vowel.drop(['row.names','y'],axis=1,inplace=True)
y_vowel=data_vowel.y
Logistic regression¶
In [3]:
X_heart.columns
Out[3]:
Index(['sbp', 'tobacco', 'ldl', 'adiposity', 'typea', 'alcohol', 'age'], dtype='str')
One dimensional logistic regression¶
In [26]:
lr=LogisticRegression(C=1)
columns=['sbp']
X_predict=X_heart[columns]
#columns=['sbp', 'tobacco', 'ldl', 'adiposity', 'typea', 'alcohol', 'age']
X=X_heart[columns]
y=y_heart
lr.fit(X,y)
plt.plot(X_heart[columns],lr.predict_proba(X_predict)[:,1],label='Predicted Probability')
plt.scatter(X_heart[columns],y_heart,color='red',marker=3,label='Observed')
plt.title('P(chd | {})'.format(columns))
plt.ylabel('Predicted Probability')
plt.xlabel(columns)
plt.legend()
lr.score(X_predict, y)
Out[26]:
0.6666666666666666
In [45]:
data=pd.read_csv('data/lrdata.csv',sep='\t')
data
XM=data[['X1','X2' ]]
yM=data.G
lr.fit(XM,yM)
b0=lr.intercept_[0]
b1=lr.coef_[0][0]
classifier = LogisticRegression().fit(XM, yM)
disp = DecisionBoundaryDisplay.from_estimator(
lr, XM, #response_method="predict",
alpha=0.5,
)
disp.ax_.scatter(XM.X1, XM.X2, edgecolor="k",c=(yM=='y').astype(int))
plt.plot()
Out[45]:
[]
In [ ]:
## Fitted parameters
(lr.intercept_, lr.coef_)
Out[ ]:
0 3.0 1 4.0 2 2.0 3 2.0 4 3.0 5 3.5 Name: X1, dtype: float64
Multivariate Logistic Regression¶
In [7]:
X=X_heart
y=y_heart
X_predict=X_heart
lr_full=LogisticRegression()
lr_full.fit(X,y)
lr_full.predict_proba(X_predict)
lr_full.score(X_predict, y)
Out[7]:
0.7164502164502164
In [8]:
X_heart.columns
Out[8]:
Index(['sbp', 'tobacco', 'ldl', 'adiposity', 'typea', 'alcohol', 'age'], dtype='str')
In [9]:
lr_full.predict_proba([[142. , 0.0, 3.38, 16.2 , 59. , 2.62, 23. ]])
Out[9]:
array([[0.89456668, 0.10543332]])
Partial Dependence Display¶
In [10]:
from sklearn.inspection import PartialDependenceDisplay
features = [0, 1,2, (0, 2)]
PartialDependenceDisplay.from_estimator(lr_full, X, features)
Out[10]:
<sklearn.inspection._plot.partial_dependence.PartialDependenceDisplay at 0x22794a93a10>
Penalized Logistic Regression¶
In [11]:
lr_full=LogisticRegression(penalty='l2')
lr_full.fit(X,y)
coef=pd.Series(lr_full.coef_[0],index=lr_full.feature_names_in_)
coef
Out[11]:
sbp 0.005275 tobacco 0.073312 ldl 0.190723 adiposity -0.011366 typea 0.038391 alcohol 0.001692 age 0.055438 dtype: float64
Vowel Data¶
- more categories, more fun
In [12]:
from sklearn import preprocessing
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
scaler = preprocessing.StandardScaler().fit(X_vowel)
X_scaled = scaler.transform(X_vowel)
lr_vowel=LogisticRegression(penalty='l2')
lr_vowel.fit(X_scaled, y_vowel)
pipe = make_pipeline(StandardScaler(), LogisticRegression())
pipe.fit(X_vowel, y_vowel);
In [13]:
from mlxtend.plotting import plot_decision_regions
fill=dict(pd.Series(X_scaled.mean(axis=0)[2:],index=range(2,10)))
plot_decision_regions(X_scaled, y_vowel.values,
clf=lr_vowel, filler_feature_values=fill,legend=2)
Out[13]:
<Axes: >
In [14]:
X_scaled.mean(axis=0)[2:]
Out[14]:
array([ 0.00000000e+00, 5.38289951e-17, -2.69144976e-17, 2.69144976e-17,
-1.34572488e-17, -8.07434927e-17, 0.00000000e+00, -2.69144976e-17])
Nearest Neighbours¶
In [15]:
from sklearn.neighbors import KNeighborsClassifier
from IPython.display import display, HTML
clf = KNeighborsClassifier(n_neighbors=15)
clf.fit(X_vowel, y_vowel) ################### fit
plot_decision_regions(X_scaled, y_vowel.values,
clf=clf, filler_feature_values=fill,legend=2)
plt.title('Train data Knn with K={}'.format(2))
Out[15]:
Text(0.5, 1.0, 'Train data Knn with K=2')
Linear Discriminant Analysis¶
In [16]:
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
#sklearn.lda.LDA(solver='svd', shrinkage=None, priors=None, n_components=None, store_covariance=False, tol=0.0001)
lda=LinearDiscriminantAnalysis()
lda.fit(X_vowel, y_vowel)
plot_decision_regions(X_scaled, y_vowel.values,
clf=lda, filler_feature_values=fill,legend=2)
plt.title('Linear Discriminant Analysis');
Sklearn just splines, GAM only in statsmodels¶
- spline interpolation example
In [17]:
def f(x):
x_points = [ 0, 1, 2, 3, 4, 5]
y_points = [12,14,22,39,27,15]
tck = sp.interpolate.make_interp_spline(x_points, y_points,k=3)
return interpolate.splev(x, tck)
print(f(1.25))
14.718750000000004
GAM - regression¶
In [18]:
import statsmodels.api as sm
from statsmodels.gam.api import GLMGam, BSplines
from statsmodels.gam.generalized_additive_model import LogitGam
heart=pd.concat([X, y], axis=1)
# create spline basis for weight and hp
x_spline = heart[['sbp','tobacco']]
bs = BSplines(x_spline, df=[12, 10], degree=[3, 3])
# penalization weight
alpha = np.array([21833888.8, 6460.38479])
gam_bs = GLMGam.from_formula('chd ~ sbp + tobacco', data=heart, smoother=bs, alpha=alpha)
res_bs = gam_bs.fit()
print(res_bs.summary())
Generalized Linear Model Regression Results
==============================================================================
Dep. Variable: chd No. Observations: 462
Model: GLMGam Df Residuals: 458.47
Model Family: Gaussian Df Model: 2.53
Link Function: Identity Scale: 0.20251
Method: PIRLS Log-Likelihood: -284.87
Date: Tue, 03 Mar 2026 Deviance: 92.842
Time: 20:56:13 Pearson chi2: 92.8
No. Iterations: 3 Pseudo R-squ. (CS): 0.1180
Covariance Type: nonrobust
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
Intercept -0.1898 8.33e+04 -2.28e-06 1.000 -1.63e+05 1.63e+05
sbp 0.0030 821.411 3.71e-06 1.000 -1609.933 1609.939
tobacco 0.0171 9.98e+04 1.71e-07 1.000 -1.96e+05 1.96e+05
sbp_s0 9.371e-05 4654.662 2.01e-08 1.000 -9122.970 9122.970
sbp_s1 0.0002 1.1e+04 1.92e-08 1.000 -2.15e+04 2.15e+04
sbp_s2 0.0004 1.83e+04 2.31e-08 1.000 -3.6e+04 3.6e+04
sbp_s3 0.0006 2.22e+04 2.52e-08 1.000 -4.35e+04 4.35e+04
sbp_s4 0.0007 2.55e+04 2.74e-08 1.000 -4.99e+04 4.99e+04
sbp_s5 0.0009 2.93e+04 3.02e-08 1.000 -5.74e+04 5.74e+04
sbp_s6 0.0012 3.45e+04 3.44e-08 1.000 -6.77e+04 6.77e+04
sbp_s7 0.0017 4.22e+04 4.09e-08 1.000 -8.27e+04 8.27e+04
sbp_s8 0.0036 6.3e+04 5.77e-08 1.000 -1.23e+05 1.23e+05
sbp_s9 0.0058 8.13e+04 7.14e-08 1.000 -1.59e+05 1.59e+05
sbp_s10 0.0075 9.61e+04 7.83e-08 1.000 -1.88e+05 1.88e+05
tobacco_s0 -0.0008 7548.636 -1.04e-07 1.000 -1.48e+04 1.48e+04
tobacco_s1 0.0010 6150.103 1.59e-07 1.000 -1.21e+04 1.21e+04
tobacco_s2 0.0086 4.42e+04 1.94e-07 1.000 -8.66e+04 8.66e+04
tobacco_s3 0.0273 1.44e+05 1.9e-07 1.000 -2.82e+05 2.82e+05
tobacco_s4 0.0526 2.89e+05 1.82e-07 1.000 -5.66e+05 5.66e+05
tobacco_s5 0.0851 5.09e+05 1.67e-07 1.000 -9.97e+05 9.97e+05
tobacco_s6 0.2021 1.45e+06 1.4e-07 1.000 -2.84e+06 2.84e+06
tobacco_s7 0.1872 2.33e+06 8.03e-08 1.000 -4.57e+06 4.57e+06
tobacco_s8 0.2077 3.11e+06 6.68e-08 1.000 -6.09e+06 6.09e+06
==============================================================================
In [19]:
res_bs.plot_partial(0, cpr=True);
Poisson¶
In [20]:
data_heart=pd.read_csv('data/heart_deasease.csv')
data_heart.drop(['row.names','famhist', 'obesity'],axis=1,inplace=True)
data_heart.columns
#'sbp', 'tobacco', 'ldl', 'adiposity', 'typea', 'alcohol', 'age'
x_spline = data_heart[['ldl','age', 'adiposity','typea']]
bs = BSplines(x_spline, df=[5,5,5,5], degree=[3, 3,3,3])
gam_bs = GLMGam.from_formula('chd ~ ldl + age + adiposity + typea', data=data_heart,
smoother=bs)#, alpha=alpha)
res_bs = gam_bs.fit()
#alpha = np.array((10000000.0, 10000000.0, 215.44346900318845, 46415.88833612782))
gam_bs_poiss = GLMGam.from_formula('chd ~ ldl + age + adiposity + typea', data=data_heart,
smoother=bs, #alpha=alpha,
family=sm.genmod.families.family.Poisson())
#statsmodels.gam.generalized_additive_model.LogitGam
res_bs_poiss = gam_bs_poiss.fit()
#print(res_bs_poiss.summary())
#gam_bs.select_penweight()[0]
#gam_bs.select_penweight_kfold()[0]
In [21]:
for i in range(3):
res_bs.plot_partial(i, cpr=True)
LogitGam¶
In [22]:
from patsy import dmatrices
alpha = np.array((1.0, 1.0, 1.0))#, 215.4446900318845, 46415.88833612782))
#gam_bs = LogitGam.from_formula('chd ~ ldl + age + adiposity + typea', data=data_heart)#, smoother=bs,alpha=alpha)
y,X=dmatrices('chd ~ ldl + age + adiposity', data=data_heart,return_type='dataframe')
x_spline2 = data_heart[['ldl','age','adiposity']]
bs2 = BSplines(x_spline2, df=[8,8,8], degree=[3, 3,3])
log_m=LogitGam(y, smoother=bs2,alpha=alpha)
#res_lm=log_m.fit()
fit=log_m.fit_regularized()
Optimization terminated successfully (Exit mode 0)
Current function value: 0.5426514025693134
Iterations: 174
Function evaluations: 174
Gradient evaluations: 174
In [23]:
for i in range(2):
res_bs.plot_partial(i, cpr=True)