儘管眼下十分艱難、可日後這段經歷說不定就會開花結果.
「碎碎念念」：大家好、我是風一；數據建模時、依據不同的場景、數據結構等等；往往會需要構建不同的數據模型、因此，就需要做模型對比，而比較模型時，評估性能是最重要的一環。在挑選最佳模型時， ANOVA(線性模型)、偏差觀察法(GLM) 和 交叉驗證法都是測量模型誤差和衡量模型性能的常用方法、下面開始模型診斷的學習。
一、簡介創建模型是持續性活動。當向模型中添加或刪除變量時，需要設法比較模型，並需要統一的方法來衡量模型的性能。二、殘差(單個模型)模型的殘差指實際觀測值與模型估計值之差。import pandas as pd
housing = pd.read_csv("housing_renamed.csv")housing.head()
neighborhoodtypeunitsyear_builtsq_ftincomeincome_per_sq_ftexpenseexpense_per_sq_ftnet_incomevaluevalue_per_sq_ftboro0FINANCIALR9-CONDOMINIUM421920.036500133261536.513420059.379906107300000200.00Manhattan1FINANCIALR4-CONDOMINIUM781985.0126420663325752.47176229513.94487096230690000242.76Manhattan2FINANCIALRR-CONDOMINIUM500NaN5541741731000031.2435430006.391376700090970000164.15Manhattan3FINANCIALR4-CONDOMINIUM2821930.02490761177631347.28278467011.18899164367556006271.23Manhattan4TRIBECAR4-CONDOMINIUM2391985.02194951000458245.58278319712.68722138554320996247.48Manhattan首先擬合一個多元線性回歸模型(含三個協變量、也就是自變量)import statsmodels
import statsmodels.api as sm
import statsmodels.formula.api as smf# 因變量在 ~ 左邊、 自變量在 ~ 右邊
house1 = smf.glm('value_per_sq_ft ~ units + sq_ft + boro',
data = housing).fit()house1.summary()
Generalized Linear Model Regression ResultsDep. Variable:value_per_sq_ftNo. Observations:2626Model:GLMDf Residuals:2619Model Family:GaussianDf Model:6Link Function:identityScale:1879.5Method:IRLSLog-Likelihood:-13621.Date:Sun, 18 Jul 2021Deviance:4.9224e+06Time:06:19:39Pearson chi2:4.92e+06No. Iterations:3
Covariance Type:nonrobust
coefstd errzP>|z|[0.0250.975]Intercept43.29095.3308.1220.00032.84553.737boro[T.Brooklyn]34.56215.5356.2440.00023.71445.411boro[T.Manhattan]130.99245.38524.3270.000120.439141.546boro[T.Queens]32.99375.6635.8270.00021.89544.092boro[T.Staten Island]-3.63039.993-0.3630.716-23.21615.956units-0.18810.022-8.5110.000-0.231-0.145sq_ft0.00022.09e-0510.0790.0000.0000.000繪製模型的殘差圖import seaborn as sns
import matplotlib.pyplot as pltfig, ax = plt.subplots()
# 繪製數據和線性回歸模型擬合 (官方文檔：https://www.cntofu.com/book/172/docs/28.md)
ax = sns.regplot(x = house1.fittedvalues, # 擬合值
y = house1.resid_deviance, # 殘留偏差
fit_reg = False) # 如果為 True，則估計並繪製與 x 和 y 變量相關的回歸模型
plt.show()
該殘差圖效果不佳，因為其中包含了明顯的群集和組。下面通過 boro 變量為圖著色,以不同顏色表示紐約各個區res_df = pd.DataFrame({
'fittedvalues': house1.fittedvalues,
'resid_deviance': house1.resid_deviance,
'boro': housing['boro']
})
# 在FacetGrid對象上繪製數據和回歸模型(官方文檔：https://www.bookstack.cn/read/seaborn-0.9/docs-27.md)
fig = sns.lmplot(x = 'fittedvalues', y = 'resid_deviance',
data = res_df, hue = 'boro', # hue: 定義數據子集的變量，將在網格中的不同構面上繪製
fit_reg = False)
plt.show()
Q-Q 圖Q-Q 圖用於判斷數據是否符合某個參考分布。許多模型都是假設數據是正態分布的，因此 Q-Q 圖常用於檢驗數據是否來自正態分布from scipy import statsresid = house1.resid_deviance.copy()
resid_std = stats.zscore(resid) # 計算樣本中每個值的 z 分數，相對於樣本平均值和標準偏差
fig = statsmodels.graphics.gofplots.qqplot(resid, line = 'r')
plt.show()
繪製殘存直方圖，觀察數據是否呈正態分布# fig, ax = plt.subplots()
ax = sns.displot(resid_std)
plt.show()
三、比較多個模型前面介紹了如何評估單個模型，有時需要比較多個模型，以便從中選出 "最佳" 模型。3.1.比較線性模型有些模型使用 "+" 運算符向模型中添加協變量，而有些模型則使用 "" 運算符 進行添加。當在模型中指定交互時應使用 "" 操作符。這意味著交互變量的行為不是彼此獨立的，它們的值相互影響，而不是簡單地相加。# 擬合五個模型
f1 = 'value_per_sq_ft ~ units + sq_ft + boro'
f2 = 'value_per_sq_ft ~ units * sq_ft + boro'
f3 = 'value_per_sq_ft ~ units + sq_ft * boro + type'
f4 = 'value_per_sq_ft ~ units + sq_ft * boro + sq_ft * type'
f5 = 'value_per_sq_ft ~ boro + type'
house1 = smf.ols(f1, data = housing).fit()
house2 = smf.ols(f2, data = housing).fit()
house3 = smf.ols(f3, data = housing).fit()
house4 = smf.ols(f4, data = housing).fit()
house5 = smf.ols(f5, data = housing).fit()
針對所有模型、獲取所有係數和相關模型house1.paramsIntercept 43.290863
boro[T.Brooklyn] 34.562150
boro[T.Manhattan] 130.992363
boro[T.Queens] 32.993674
boro[T.Staten Island] -3.630251
units -0.188115
sq_ft 0.000210
dtype: float64mod_results = pd.concat([house1.params, house2.params, house3.params, house4.params, house5.params],
axis = 1).rename(columns = lambda x: 'house' + str(x + 1)).reset_index().rename(columns = {'index': 'param'})mod_results.head(5)
paramhouse1house2house3house4house50Intercept43.29086344.45183840.23269242.92515838.3558791boro[T.Brooklyn]34.56215032.31498728.43926927.59763928.2309042boro[T.Manhattan]130.992363127.098918115.711927114.434920126.8014463boro[T.Queens]32.99367429.79931828.55323327.06885225.1621084boro[T.Staten Island]-3.630251-7.542656-15.640961-15.941330-16.950990'''
'寬數據'轉換為 '長數據'
melt函數有如下幾個參數：
id_vars: 該參數是一個容器(列表、元組或ndarray),所表示的變量會保持原樣。
value_vars: 指定想"融合"(或轉換為行)的列。它默認會"融合"未在 id_vars 參數指定的所有列。
var_name: 該字符串用於指定 value_vars 融合後的新列名。默認為 variable。
value_name: 該字符串為新列名,代表 var_name 的值。默認為 value。
'''
mod_results = mod_results.melt(
id_vars = 'param',
var_name = 'model'
)mod_results.head()
parammodelvalue0Intercepthouse143.2908631boro[T.Brooklyn]house134.5621502boro[T.Manhattan]house1130.9923633boro[T.Queens]house132.9936744boro[T.Staten Island]house1-3.630251
係數圖形化、直觀展現模型是如何評估彼此相關參數的fig, ax = plt.subplots()
ax = sns.pointplot(
x = 'value', y = 'param', hue = 'model',
data = mod_results,
dodge = True, join = False
)
plt.tight_layout()
plt.show()
這樣就得到了線性模型。下面可以使用方差分析(ANOVA)方法來比較它們。方差分析會給出殘差平方和(RSS), 可以通過它來評估模型的性能(殘差平方和越小，模型的擬合效果越好)。
方差分析(Analysis of Variance，簡稱ANOVA)，又稱「變異數分析」，為數據分析中常見的統計模型，主要為探討連續型（Continuous）因變量（Dependent variable）與類別型自變量（Independent variable）的關係。當自變量的因子等於或超過三個類別時，檢驗各類別平均值是否相等，採用方差分析。model_names = ['house1', 'house2', 'house3', 'house4', 'house5']
house_anova = statsmodels.stats.anova.anova_lm(house1, house2, house3, house4, house5)
house_anova.index = model_names
house_anova.reset_index()
indexdf_residssrdf_diffss_diffFPr(>F)0house12619.04.922389e+060.0NaNNaNNaN1house22618.04.884872e+061.037517.43760520.0390497.912333e-062house32612.04.619926e+066.0264945.53999423.5857282.754431e-273house42609.04.576671e+063.043255.4411927.7012894.025581e-054house52618.04.901463e+06-9.0-324791.84790719.275539NaN評估模型的另一種方法是使用赤池信息量準則(AIC) 和 貝葉斯信息準則 (BIC)。這些方法均引入了與模型參數個數相關的懲罰項。因此，應儘量在模型的性能和簡潔性之間做好平衡(簡單優於複雜)。3.2.比較 GLM可以採用同樣的方法評估和診斷GLM。不過,方差分析只評估模型的偏差。def anova_deviance_table(*models):
return pd.DataFrame({
'df_residuals': [i.df_resid for i in models],
'resid_stddev': [i.deviance for i in models],
'df': [i.df_model for i in models],
'deviance': [i.deviance for i in models]
})f1 = 'value_per_sq_ft ~ units + sq_ft + boro'
f2 = 'value_per_sq_ft ~ units * sq_ft + boro'
f3 = 'value_per_sq_ft ~ units + sq_ft * boro + type'
f4 = 'value_per_sq_ft ~ units + sq_ft * boro + sq_ft * type'
f5 = 'value_per_sq_ft ~ boro + type'
glm1 = smf.glm(f1, data = housing).fit()
glm2 = smf.glm(f2, data = housing).fit()
glm3 = smf.glm(f3, data = housing).fit()
glm4 = smf.glm(f4, data = housing).fit()
glm5 = smf.glm(f5, data = housing).fit()
glm_anova = anova_deviance_table(glm1, glm2, glm3, glm4, glm5)
glm_anova
df_residualsresid_stddevdfdeviance026194.922389e+0664.922389e+06126184.884872e+0674.884872e+06226124.619926e+06134.619926e+06326094.576671e+06164.576671e+06426184.901463e+0674.901463e+06四、k 折交叉驗證交叉驗證是另外一種模型比較方法。它可以解釋模型在新數據上的表現。這種方法會把數據分成 k 個部分，把其中一個部分用作測試集，把其餘部分用作訓練集來擬合模型(往往是測 2 訓 8)。模型擬合好之後，使用測試集進行測試，並計算誤差。不斷重複這個過程, 直到 k 個部分都測試過。模型的最終誤差是所有模型的平均值。進行交叉驗證有多種方法。剛才介紹的方法稱為 "k 折交叉驗證"。而 "留一交叉驗證"方法每次只留下一個樣本用作測試集，其他樣本都用作訓練集。from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegressionhousing.columnsIndex(['neighborhood', 'type', 'units', 'year_built', 'sq_ft', 'income',
'income_per_sq_ft', 'expense', 'expense_per_sq_ft', 'net_income',
'value', 'value_per_sq_ft', 'boro'],
dtype='object')# 獲取訓練數據和測試數據
X_train, X_test, y_train, y_test = train_test_split(
pd.get_dummies(housing[['units', 'sq_ft', 'boro']], drop_first = True), # 刪除 k 個分類級別中的 k-1 個虛擬變量第一級
housing['value_per_sq_ft'],
test_size = 0.20,
random_state = 42
)lr = LinearRegression().fit(X_train, y_train)
# 算出一個分數，衡量模型在測試數據上的表現
lr.score(X_test, y_test)0.6137125285030869
由於 sklearn 高度依賴於 NumPy ndarray，所以 patsy 庫允許指定一個公式，比如 statsmodels 中的公式 API, 並且會返回合適的 NumPy 數組，以便在 sklearn 中使用。from patsy import dmatricesy, X = dmatrices('value_per_sq_ft ~ units + sq_ft + boro', housing, return_type = "dataframe")
X_train, X_test, y_train, y_test = train_test_split(
X, y, test_size = 0.20, random_state = 42
)
lr = LinearRegression().fit(X_train, y_train)
lr.score(X_test, y_test)0.613712528503469
為了執行 k 折交叉驗證，需要先從 sklearn 庫導入該函數from sklearn.model_selection import KFold, cross_val_score# 重新獲取住房數據集
housing = pd.read_csv("housing_renamed.csv")
依據數據的行數，指定需要多少折；如果數據中包含的觀測值不太多，可以選擇一個較小的 k (例如 2)。k 值通常介於 5 到 10 之間。但請注意, k 值越大，所需要的計算時間就越長，需根據實際情況選擇合適的 k。kf = KFold(n_splits = 5)
y, X = dmatrices('value_per_sq_ft ~ units + sq_ft + boro', housing)
# 每一折上訓練並測試模型
coefs = []
scores = []
for train, test in kf.split(X):
X_train, X_test = X[train], X[test]
y_train, y_test = y[train],y[test]
lr = LinearRegression().fit(X_train, y_train)
coefs.append(pd.DataFrame(lr.coef_))
scores.append(lr.score(X_test, y_test))coefs_df = pd.concat(coefs)
coefs_df.columns = X.design_info.column_names
coefs_df
Interceptboro[T.Brooklyn]boro[T.Manhattan]boro[T.Queens]boro[T.Staten Island]unitssq_ft00.033.369037129.90401132.103100-4.381085-0.2058900.00022000.032.889925116.95738531.295956-4.919232-0.1461800.00015500.030.975560141.85932732.043449-4.379916-0.1796710.00019400.041.449196130.77901333.050968-3.430209-0.2079040.00023200.0-38.51191556.069855-17.5579390.000000-0.1458290.000202# 使用 apply 和 np.mean 函數查看所有折的平均係數
import numpy as np
coefs_df.apply(np.mean)Intercept 0.000000
boro[T.Brooklyn] 20.034361
boro[T.Manhattan] 115.113918
boro[T.Queens] 22.187107
boro[T.Staten Island] -3.422088
units -0.177095
sq_ft 0.000201
dtype: float64
查看每個模型的分數。每個模型都有默認的計分方法。例如 LinearRegression 使用 R^2(決定係數)回歸評分函數。scores[0.027314162909376938,
-0.553836221218619,
-0.15636371688031647,
-0.3234202061860394,
-1.6929655586236945]# 使用 cross_val_scores(交叉驗證分數) 計算分數
model = LinearRegression()
scores = cross_val_score(model, X, y, cv = 5)
scoresarray([ 0.02731416, -0.55383622, -0.15636372, -0.32342021, -1.69296556])綜上可見、模型4：模型4：f4 = 'value_per_sq_ft ~ units + sq_ft * boro + sq_ft * type' 性能比較好。
特此申明：以上文章來源於書籍《Python 數據分析活用 Pandas 庫》的閱讀理解以及練習筆記。
往期文章：
Python數據建模：簡單線性回歸
Python數據建模：廣義線性模型