【实验二】基于回归分析的大学综合得分预测

目的

根据 CWUR 所提供的世界各地知名大学各方面的排名（师资、科研等），一方面通过数据可视化的方式观察不同大学的特点，另一方面希望构建机器学习模型（线性回归）预测一所大学的综合得分。数据来源

方法

使用线性回归模型，尝试使用梯度下降法，最小二乘法实现回归预测，学习使用sklearn便捷的实现线性回归，并做预测结果与模型系数进行评价分析

过程

一、导入工具包

import pandas as pd
import numpy as np
from sklearn.model_selection import train_test_split, cross_validate
import torch
import torch.utils.data as Data
from sklearn import linear_model
from sklearn.metrics import mean_squared_error, r2_score

除去上次实验已经用到过的一些工具包，这次新导入了torch，sklearn中的linear_model，mean_squared_error和r2_score
torch的引入是因为我之前熟悉过相关操作，所以使用其进行矩阵运算
linear_model是sklearn中集成的线性回归工具模块，可以简单快捷的进行运算
mean_squared_error和r2_score则是用来便捷计算RMSE与R^2，具体用法会在下文提及

二、导入数据

data_df = pd.read_csv('./cwurData.csv')  # 读入 csv 文件为 pandas 的 DataFrame
#print(data_df.head(3).T)  # 观察前几列并转置方便观察
data_df = data_df.dropna()  # 舍去包含 NaN 的 row
feature_cols = ['quality_of_faculty', 'publications', 'citations', 'alumni_employment',
                'influence', 'quality_of_education', 'broad_impact', 'patents']
X = data_df[feature_cols]
Y = data_df['score']
all_y = Y.values
all_x = X.values
x_train, x_test, y_train, y_test = train_test_split(all_x, all_y, test_size=0.2)#将数据分为训练集与测试集

上述代码实现了数据导入，特征选择与对训练集与数据集的切分
其中特征选择部分运用的索引方式之前没有学习过，本次进行重点记忆

三、训练模型

本部分为线性回归核心部分，以下分别采用梯度下降，最小二乘，sklearn三种方式实现

1.梯度下降法

for i in range(1,num_epoch+1):
        y_pre = torch.mm(feature,w_pre)
        loss = y_pre - lable.view([1600,1])
        loss = loss/len(feature)
        loss = torch.mm(torch.t(feature),loss)
        w_pre = w_pre - lr*loss

此处代码为核心部分，不包含对数据的进一步处理方式，完整代码详见文件“最小二乘法.ipynb”
梯度下降法利用了均方误差单峰的特点，每一次循环都使参数朝着梯度方向下降一段距离。
但此时如果学习率设置的过大，超过了某一限制，最后结果就会出现错误。经过询问助教与上网查询的方法，了解到这叫做overshoot minimum现象，这种现象会使得参数震荡偏离最低点，最终导致错误出现。
为解决这一问题，我们需要控制学习率的大小，但是这是另一个问题就又出现了：如果直接对原数据进行学习，学习率需要设置的非常非常小，数量级在小数点后六位，最后训练出的模型RMSE数值也达到了12左右，毫无疑问需要继续进行优化。
这个问题最终在群内交流时得到了解决办法：归一化。

from sklearn import preprocessing
min_max_scaler = preprocessing.MinMaxScaler()
all_x = min_max_scaler.fit_transform(all_x)

这是归一化的实现方式，借助了sklearn中的preprocessing，将所有的x数据映射到了0-1中。
经过归一化，学习率设置到0.3即可让程序顺利运行，并且RMSE也降到了4左右，R^2在0.6左右，与后面的sklearn实现相差无几

2.最小二乘法

w_pre = torch.mm(torch.mm(torch.inverse(torch.mm(torch.t(feature),feature)),torch.t(feature)),lable.view(1600,1))

此处为最小二乘法计算部分，完整代码详见“最小二乘法.ipynb”
最小二乘法的实现较为简单，直接使用公式计算即可，公式为(((X^T)X)-1)(X^T)Y
该方法得到的模型较为优秀，RMSE达到了2左右，但R^2仍然在0.6左右

3.sklearn实现

reg = linear_model.LinearRegression()
reg.fit(x_train,y_train)#训练
y_pre=reg.predict(x_test)#预测
rmse = mean_squared_error(y_test,y_pre,squared = False)
print(rmse)
print(r2_score(y_test,y_pre))#评估

不得不说sklearn的却是方便快捷，短短几行就可以完成训练，预测与评估。完整代码见“sklearn实现.ipynb”
首先第一行是定义一个线性模型，接着第二行实现了模型的训练，第三行便完成了预测。
第四行的RMSE计算中有一个参数需要注意：squared。当这个参数为Ture时，函数计算的是MSE，为False时则计算RMSE。
最后一行便完成了R^2的计算。
最后的结果是RMSE大概在4左右，而R^2在0.6左右。

四、分析模型系数

经过查阅有关资料，线性回归模型的解释方法主要有三种：数值特征的解释，分类特征的解释与特征重要性。由于本实验特征均为数值特征，特征重要性的计算方法我尚未完全掌握，因此直接采用第一种方法。
该方法对于特征的解释是：当所有其他特征保持不变时，特征x增加一个单位，则预测结果y增加b。只是我们就可以根据b来判断特征的重要程度。由于本实验较为简单，这种方法等同于直接比较特征值。
经过计较，quality_of_education一项对于结果单影响远大于其他特征，而其他特征中publications，citations较为重要，broad_impact与patents对结果有微弱影响，剩余特征几乎不对结果造成影响。

结果

经过调整，该实验的R^2在各种实现方式中均维持在0.6左右，RMSE在最小二乘法中表现更优秀，维持在2左右，而梯度下降与sklearn辅助实现则在4左右。
特征分析详见过程第四步

思考

本实验最难以理解的其实是归一化对模型的优化程度与overshoot minimum现象的出现，虽然问题得以解决，但这两者之间的关系与具体形成原因尚没有完全明确。而这两者又对超参数的设置有很大的影响，因此目前仍然需要对该问题进一步分析研究。

结尾附代码

#最小二乘
import pandas as pd
import numpy as np
from sklearn.model_selection import train_test_split, cross_validate
import torch

data_df = pd.read_csv('./cwurData.csv')  # 读入 csv 文件为 pandas 的 DataFrame
#print(data_df.head(3).T)  # 观察前几列并转置方便观察
data_df = data_df.dropna()  # 舍去包含 NaN 的 row
feature_cols = ['quality_of_faculty', 'publications', 'citations', 'alumni_employment',
                'influence', 'quality_of_education', 'broad_impact', 'patents']
X = data_df[feature_cols]
Y = data_df['score']
all_y = Y.values
all_x = X.values
x_train, x_test, y_train, y_test = train_test_split(all_x, all_y, test_size=0.2)#将数据分为训练集与测试集
feature = torch.tensor(x_train,dtype=torch.float)
lable = torch.tensor(y_train,dtype=torch.float)#numpy转tensor
bais = torch.ones(feature.shape[0])
bais = bais.view([1600,1])
feature = torch.cat((feature,bais),1)#为特征增添偏移量
x_t = torch.tensor(x_test,dtype=torch.float)
baiss = torch.ones(x_t.shape[0])
baiss = baiss.view([len(x_t),1])
x_t = torch.cat((x_t,baiss),1)
w_pre = torch.mm(torch.mm(torch.inverse(torch.mm(torch.t(feature),feature)),torch.t(feature)),lable.view(1600,1))#核心，用最小二乘计算w
print(list(w_pre))
loss = torch.tensor(y_test).view([400,1]) - torch.mm(x_t,w_pre)
loss = torch.mm(torch.t(loss),loss)
RMSE = (loss.item()/len(lable))**0.5
print(RMSE)
var = (torch.tensor(y_test)-torch.mean(torch.tensor(y_test))).view([-1,1])
var = torch.mm(torch.t(var),var)
R = 1-loss.item()/var.item()
print(R)

#梯度下降
import pandas as pd
import numpy as np
from sklearn.model_selection import train_test_split, cross_validate
import torch
from sklearn import preprocessing

data_df = pd.read_csv('./cwurData.csv')  # 读入 csv 文件为 pandas 的 DataFrame
#print(data_df.head(3).T)  # 观察前几列并转置方便观察
data_df = data_df.dropna()  # 舍去包含 NaN 的 row
feature_cols = ['quality_of_faculty', 'publications', 'citations', 'alumni_employment',
                'influence', 'quality_of_education', 'broad_impact', 'patents']
X = data_df[feature_cols]
Y = data_df['score']
all_y = Y.values
all_x = X.values
min_max_scaler = preprocessing.MinMaxScaler()
all_x = min_max_scaler.fit_transform(all_x)#数据归一化
x_train, x_test, y_train, y_test = train_test_split(all_x, all_y, test_size=0.2)#将数据分为训练集与测试集
feature = torch.tensor(x_train,dtype=torch.float)
lable = torch.tensor(y_train,dtype=torch.float)#numpy转tensor
bais = torch.ones(feature.shape[0])
bais = bais.view([1600,1])
feature = torch.cat((feature,bais),1)#为特征增添偏移量
x_t = torch.tensor(x_test,dtype=torch.float)
baiss = torch.ones(x_t.shape[0])
baiss = baiss.view([len(x_t),1])
x_t = torch.cat((x_t,baiss),1)
w_pre = torch.rand(len(feature_cols)+1,1)
num_epoch = 500
lr = 0.3
for i in range(1,num_epoch+1):
        y_pre = torch.mm(feature,w_pre)
        loss = y_pre - lable.view([1600,1])
        loss = loss/len(feature)
        loss = torch.mm(torch.t(feature),loss)
        w_pre = w_pre - lr*loss#核心部分，进行梯度下降
y_pre = torch.mm(x_t,w_pre)
loss = torch.tensor(y_test).view([400,1]) - y_pre
loss = torch.mm(torch.t(loss),loss)
rms = (loss.item()/len(y_pre))**0.5
print(rms)#计算RMSE
var = (torch.tensor(y_test)-torch.mean(torch.tensor(y_test))).view([-1,1])
var = torch.mm(torch.t(var),var)
R = 1-loss.item()/var.item()
print(R)#计算R^2

#sklearn实现
import pandas as pd
import numpy as np
from sklearn.model_selection import train_test_split, cross_validate
import torch
import torch.utils.data as Data
from sklearn import linear_model
from sklearn.metrics import mean_squared_error, r2_score

data_df = pd.read_csv('./cwurData.csv')  # 读入 csv 文件为 pandas 的 DataFrame
#print(data_df.head(3).T)  # 观察前几列并转置方便观察
data_df = data_df.dropna()  # 舍去包含 NaN 的 row
feature_cols = ['quality_of_faculty', 'publications', 'citations', 'alumni_employment',
                'influence', 'quality_of_education', 'broad_impact', 'patents']
X = data_df[feature_cols]
Y = data_df['score']
all_y = Y.values
all_x = X.values
x_train, x_test, y_train, y_test = train_test_split(all_x, all_y, test_size=0.2)#将数据分为训练集与测试集
reg = linear_model.LinearRegression()
reg.fit(x_train,y_train)#训练
y_pre=reg.predict(x_test)#预测
rmse = mean_squared_error(y_test,y_pre,squared = False)
print(rmse)
print(r2_score(y_test,y_pre))#评估