代码收藏家技术教程 2024-12-11

机器学习与python-陈强第五章逻辑回归课后习题答案

题目

5.1 使用 Hastie etal.(2009)的南非心脏病数据 SAheart.csv 进行逻辑回归”。其中响应变量为chd(是否有冠心病，即coronary heart disease)。特征向量包括sbp(收缩压，即 systolic blood pressure)、tobacco(累计抽烟量)、ldl(低密度脂蛋白胆固醇,即 low densitylipoprotein cholesterol)、因子变量 famhist(是否有家族心脏病史)、typea(A 型行为,即type-Abehavior)、alcoho1(当前饮酒量)、age(发病时年龄)，以及两个关于肥胖程度的数值型度量adiposity与obesity。
(1)计算样本中有冠心病的比例。
(2)由于数据包含分类变量，使用命令“x=pd.get_dummies(X)”将数据矩阵X中的分类变量设为虚拟变量。使用“random_state=0”，预留100个观测值作为测试集。
(3)在训练集中，使用skleamn模块的LogisticRegression类(设参数“c=1e10”与“fit intercept=False”)，将chd 对其余变量进行逻辑回归。
(4)展示此逻辑回归的回归系数。
(5)计算测试集的预测概率，并展示前5个预测概率。
(6)在测试集中进行预测，计算准确率、错误率、灵敏度、特异度与召回率。
(7)画 ROC 曲线。
(8)计算 AUC。
(9)计算kappa值。

详细代码

1.引入库

代码如下（示例）：

import pandas as pd
import numpy as np
from sklearn.model_selection import train_test_split
import statsmodels.api as sm
from patsy import dmatrices
import matplotlib.pyplot as plt
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import confusion_matrix
from sklearn.metrics import classification_report
from sklearn.metrics import cohen_kappa_score
from sklearn.metrics import roc_curve, RocCurveDisplay

2.读入数据

代码如下（示例）：

SA = pd.read_csv('SAheart.csv')
# 查看数据
SA.head()
X = pd.DataFrame(SA)

# tips
提示：数据文件要和代码文件放在同一文件夹里面

一：计算样本中冠心病的比例

total_count = len(X) #总样本量
chd_count= np.sum(SA.chd) #患病数量
f=chd_count/total_count
print('冠心病比例为',f)

结果：
冠心病比例为 0.3463203463203463

二：将数据矩阵X中的分类变量设为虚拟变量；预留100个观测值作为测试集

X = pd.DataFrame(SA)
X_xuni=pd.get_dummies(X)
#虚拟变量会生成两列，保留一列即可（得冠心病）
X_xuni.drop(labels=['famhist_Absent'],axis=1,inplace=True)

# tips
random_state=0 无论你何时运行代码，只要使用相同的输入数据和相同
的 random_state 值，你都会得到相同的随机数序列或结果。

#分层抽样 生成训练集和测试集（保证两边所含类别比例相同)
train, test =  train_test_split(X_xuni, test_size=100, stratify=X_xuni.chd, random_state=0)

#测试集与训练集
ytest =test.iloc[:,8] #索引是从 0 开始计数的，这意味着第一个元素的位置是 0
ytrain=train.iloc[:,8]
xtest=test.iloc[:, [0,1,2,3,4,5,6,7,9]]
xtrain=train.iloc[:, [0,1,2,3,4,5,6,7,9]]
# tips
#索引是从 0 开始计数的，这意味着第一个元素的位置是 0

结果：
X_xuni为全是数值型变量得数据表

三：chd 对其余变量进行逻辑回归

其中，参数“c=1e10”为惩罚力度的倒数，故惩罚力度为10-0，非常接近于0(可视为没有惩罚)。使用 skleam与 statsmodels 所估计的回归系数非常接近，但截距项差别较大。尽管此截距项没有实质性影响，若要在skleam 模块得到与 statsmodels 一样的截距项，可输入参数“fit intercept=False”:

model =  LogisticRegression(C=1e10, fit_intercept=False,max_iter=1000, solver='lbfgs')
model.fit(xtrain, ytrain)
# tips
若按书上的代码 LogisticRegression(C=1e10, fit_intercept=False)会报错，也许是因为迭代次数不够，进而增加一些参数
不过加不加参数，回归系数都是一样的

四：展示此逻辑回归的回归系数

具体系数的实际意义是否合理，请读者自行判断

model.coef_
# 使用zip组合数据并使用格式化字符串输出
for name, age in zip(X_xuni.columns, model.coef_[0]):  # 注意这里使用model.coef_[0]，因为coef_是一个二维数组
    print(f'{name:10} {age:8.6f}')  # 调整了字段宽度，使得浮点数显示更加合理

结果：

五：计算测试集的预测概率，并展示前5个预测概率

获取标签类别，为了知道输出的每一列代表什么。0-没得冠心病，1-得冠心病。

model.score(xtest, ytest) #预测准确率

# 获取类别标签
class_labels = model.classes_

# 打印类别标签及其对应的列索引
print("Class labels and their corresponding column indices:")
for i, label in enumerate(class_labels):
    print(f"Column {i}: {label}")

# 引入更具描述性的变量名
predicted_probabilities = model.predict_proba(xtest)

# 输出预测概率
print("Predicted probabilities:")
print(predicted_probabilities[:5])

结果：

# 直接输出是否有冠心病
pred = model.predict(xtest)
pred[:5]

结果：

六：在测试集中进行预测，计算准确率、错误率、灵敏度、特异度与召回率

ytest真实数据，pred预测数据

print(classification_report(ytest, pred, target_names=['Not chd', 'chd']))

结果：

精确度（Precision）:
对于“无冠心病”类别，精确度为0.79，意味着在模型预测为“无冠心病”的样本中，有79%是真正的“无冠心病”样本。
对于“冠心病”类别，精确度为0.67，意味着在模型预测为“冠心病”的样本中，有76%是真正的“冠心病”样本。

召回率（Recall）:
对于“无冠心病”类别，召回率为0.85，意味着在所有真正的“无冠心病”样本中，有85%被模型正确预测为“无冠心病”。
对于“冠心病”类别，召回率为0.57，意味着在所有真正的“冠心病”样本中，只有57%被模型正确预测为“冠心病”。

F1 分数（F1-score）:
对于“无冠心病”类别，F1分数为0.81，是精确度和召回率的调和平均数，表明模型在这个类别上的性能相对较好。
对于“冠心病”类别，F1分数为0.62，相对较低，表明模型在这个类别上的性能较差，尤其是在召回率方面。

支持度（Support）:
“无冠心病”类别的支持度为65，意味着数据集中有65个样本属于这个类别。
“无、冠心病”类别的支持度为35，意味着数据集中有35个样本属于这个类别。
准确率（Accuracy）:
整体的准确率为0.75，意味着在所有样本中，有75%被模型正确预测。

七：画 ROC 曲线

位于 ROC 曲线之下的蓝色区域即为 AUC,一般介于 0.5与1之间。
如果 AUC为1，则意味着模型对于所有正例与反例的预测都是正确的，这一般是无法达到的理想状态。
如果 ROC 曲线与 45°线(从原点到(1,1)的对角线)重合，则意味着该模型的预测结果无异于随机猜测。比如，样本中正例与负例各占一半，而通过从[0,1]区间的均匀分布随机抽样来预测概率。此时，AUC为0.5。AUC小于0.5的情形十分罕见，这意味着模型的预测结果还不如随机猜测。

# 计算 ROC 曲线的数据
yscores = model.predict_proba(xtest)[:,1]# 索引从0开始
fpr, tpr, _ = roc_curve(ytest, yscores)

# 使用 RocCurveDisplay 绘制 ROC 曲线
roc_display = RocCurveDisplay(fpr=fpr, tpr=tpr).plot()

x = np.linspace(0, 1, 100)
plt.plot(x, x, 'k--', linewidth=1)
plt.title('ROC Curve (Test Set)')
# 显示图例
plt.legend()
# 显示图形
plt.show()

结果：
有点梯度上升，不过总体模型效果还是较好的。

八：计算AUC

from sklearn.metrics import roc_auc_score
auc = roc_auc_score(ytest, yscores)
print('AUC:', auc)

结果：
AUC值为0.798，表示该分类模型在区分正样本和负样本方面具有一定的能力。这个值越接近于1，模型的性能通常越好。
通常，AUC值在0.5到0.7之间被认为是低效的，0.7到0.9之间被认为是中等的，而0.9以上则被认为是非常优秀的。因此，0.798可以被视为一个中等偏上的性能表现。

九：计算kappa值

cohen_kappa_score(ytest, pred)

结果：
Kappa值（κ）是一种衡量两个评判者（或分类模型）在评判（或分类）同一类事物时一致性的指标。其取值范围通常为0到1，其中1表示完全一致，0表示偶然一致或完全不一致（有时也包括-1表示完全不一致，但在此我们主要关注0到1的范围）。

对于Kappa=0.4318，一致性程度：Kappa值0.4318落在0.41到0.60的区间内，这通常被认为是“中等的一致性”（Moderate Agreement）。
当Kappa值达到中等一致性时，说明两个评判者或分类模型在大多数情况下能够给出相同或相近的评判结果或分类标签。
在实际应用中，这可能意味着该分类模型或评判标准在特定领域或数据集上具有一定的可靠性和准确性，但仍有改进的空间。

以上代码或内容是我根据当前的理解和实践编写的，虽然我已经尽力确保它们的准确性和完整性，但难免会有疏漏或不足之处。欢迎各位提出宝贵的批评和指正。
感谢您的关注和帮助，期待与您共同探讨和进步！