廣義線性模型（GLM）是一種靈活的方法，可以將不同的預測變數組合在一起以預測響應變數的值。 許多常見的統計程式是GLM的型別-t檢驗，多元迴歸，ANOVA，ANCOVA等。此外，許多更復雜的統計程式是GLM的擴充套件。 因此，瞭解GLM的工作原理非常有用。接上一篇的文章做一個詳細的程式碼解析：

內容簡介：

廣義線性模型原理

簡單的線性迴歸

設計矩陣

Statsmodels

Statsmodels formula inputs

以圖形方式探索模型擬合

分類變數和偽編碼

t-檢驗

One-sample test （mean of

high

different to zero）

超過一級的類別變數

參考資料

廣義線性模型原理

General form of the model

The general form of a GLM can be written as

= （

）+ br/> ： the design matrix， where columns are predictors， rows are observations

is the response variable we‘re trying to predict （a column vector）

is a vector of weights， with the same length as the number of columns in the design matrixbr/> is called the

linear predictor

（） is some transformation of the weighted design matrix （the link function）。 Examples include logistic， logarithmic， Poisson。。。

is random noise of some type。

用個推導圖表示吧：

公式模型理解：

Another way to write

is in vector notation：

1 1+ 2 2+。。。+

are column vectors

is a scalar

So the t

e

rm is a linear weighted sum of predictors （i。e。， it’

；s a l

inear model）。

簡單的線性迴歸

匯入依賴庫

# import necessary modules：

import

numpy

as

np

import

pandas

as

pd

import

patsy

import

seaborn

as

sns

import

statsmodels。api

as

sm

import

statsmodels。formula。api

as

smf

import

statsmodels。stats

as

sm_stats

import

matplotlib。pyplot

as

plt

%

matplotlib

inline

sns

。

set_style

（

“white”

）

取一段資料，原始檔可以文末的地步的github連結上獲取全部工程程式碼：

First， we import the data， drop missing cases， and create our log bodyweight variable。

# import the dataset csv file into a Pandas dataframe object：

sleep = pd。read_csv（‘sleep。txt’， sep=‘\t’）

sleep = sleep。dropna（axis=0） # drop on axis 0 （i。e。 rows）

sleep［‘log_BodyWt’］ = np。log（sleep［‘BodyWt’］）

sleep。info（）

g = sns。jointplot（“log_BodyWt”， “TotalSleep”， sleep）

我們可以看到，體重的對數似乎與TotalSleep近似線性相關（儘管低體重的方差更高）。

給定一種動物的體重，該動物平均要睡幾個小時？

為了回答這個問題，我們可以使用線性迴歸。 我們只有一個預測變數（體重的對數）來預測我們感興趣的反應變數（總睡眠時間）。

從中學數學可以知道，直線的方程式可以寫成 = + 。 也就是說，我們想估計兩個引數：作為function的函式的斜率和截距（當 為零時c的值）。

設計矩陣

設計矩陣columns的列中有預測變數。請記住，我們將為設計矩陣的每一列估計一個引數（ 權重）。

設計矩陣不一定只是我們資料框中列的副本。有時，需要對變數進行轉換以使其在模型中具有特定含義（如去心的否個變數）。我們可以使用名為Patsy的模組根據模型的規範自動建立設計矩陣。

Patsy和S樣式的公式

Patsy程式包使我們能夠基於古典S統計語言（現在最常在R中使用）所使用的模型，指定我們要適合高階公式語言的模型。

Patsy允許我們指定我們的模型包含我們要使用的變數的資料框並返回設計矩陣 （還可以選擇返回 ）。

公式的一般形式為

y〜x

。

〜字元左側的所有內容都是我們要預測的內容（響應變數）

右邊的一切都是預測因素

我們可以在資料框中使用變數的名稱來使程式碼更具可讀性

為了更具體一點，讓我們看一下上面的線性迴歸示例。

# this command returns both X and y （hence “dmatrices” plural）。 See “dmatrix” if you just want X。

y， X = patsy。dmatrices（‘TotalSleep ~ log_BodyWt’， sleep）

X = patsy。dmatrix（‘log_BodyWt’， sleep）

y

X

X：

請注意，X的第二列只是我們資料框中的對數體重值。但是，Patsy添加了另一列：截距，其中每個值均為1。這意味著什麼？

在進行迴歸時，請記住，我們將為每列估計 權重，並且每行的線性預測變數的值為 1 1+ 2 2 + 。。。 +。。。。因此，在

設計矩陣中增加一列表示我們將 1（截距值）的值新增到資料集中的每個觀察值

。這個值是我們方程式中 = x+ 的截距項（ ），y=mx +1。

因此，在我們的迴歸分析中，我們對每次觀察的總睡眠時間的預測為：

̂ = 1 + 2logBodyWt

1和 2是我們的模型係數，將從資料中估算出來。現在，讓我們將它們視為兩個數字，我們將從上面的散點圖中透過眼睛進行估算。

當日志體重為零時，TotalSleep大約為8。3小時。因此，讓我們設定 1 = 8。3。斜率很難估計：我猜是-4。讓我們將這些權重應用於我們的設計矩陣以生成一些預測：

betas = np。array（［8。3， -4。］） # the 。 after the number makes it a floating-point number rather than an integer。

yhat = np。inner（X， betas）

print（yhat）

請注意，取矩陣和行向量的內積（使用命令np。inner）與將每一列乘以其相應的beta權重然後求和（即

1 1 + 2 2

）相同。

我們對真實值的預測有多相似？

# this plot uses straight matplotlib commands， not seaborn：

plt。scatter（yhat， y）

plt。plot（（-30， 40）， （-30， 40）， c=“。2”， ls=“——”）；

plt。xlabel（‘yhat （predicted）’）

plt。ylabel（‘y （real）’）

plt。show（）

從散點圖（帶有虛線），可以看到我認為斜率是合理的。 讓我們計算錯誤：

SSE（差方和）

errors = y - yhat # absolute errors for each data point

squared_errors = errors ** 2 # squared errors

sse = squared_errors。sum（） # sum of squared errors

sse

Statsmodels

import statsmodels。api as sm

現在我們可以使用Statsmodels包從資料中估計最佳係數。 我們將使用statsmodels的GLM函式對其進行擬合。

fit = sm。GLM（y， X）。fit（） # fit the design matrices we specified above using Patsy。。。

# let‘s look at the fit summary：

print（fit。summary（））

上表顯示了擬合簡單線性迴歸的結果。這裡有很多資訊。讓我們來看看其中的一些。

請注意，我們已將“Gaussian”模型與“identity”連結函式配合使用。這意味著我們正在進行普通的線性迴歸，如上所述，這是廣義線性模型的一種情況。identity意味著 （）（參見上文）只是identity（即它不會以任何方式修改線性預測變數））。

我們有42個觀測量（No。 Observations）（即睡眠中所有完整的病例）。剩餘自由度（Df Residuals）等於觀測值減去2，因為我們從資料中估計了2個引數（截距和斜率）。

底部的表格向我們顯示了有關我們估計的各個引數的資訊。截距的

係數值

（Intercept）為11。4，非常接近我的猜測8。3。對數體重（log_BodyWt）的斜率是-0。95，所以我的斜率估計很垃圾。

z得分是係數除以其標準誤差（coefficient divided by its standard error），而P值正在測試該標準化係數得分是否不同於零。這兩個P值在這裡都非常重要，因此在這種簡單情況下，我們可以推斷出這兩個引數都有助於預測睡眠總量。最後兩列為我們提供了

係數估計值的置信區間

的邊界。

和我們估計的模型對比：

yhat_2 = fit。predict（）

# note that internally， the fit。predict（） function does something similar

# to our inner product between X and beta， above。

plt。scatter（yhat_2， y）

plt。plot（（-30， 40）， （-30， 40）， c=“。2”， ls=“——”）；

plt。xlabel（’yhat_2 （predicted）‘）

plt。ylabel（’y （real）‘）

plt。show（）

同理計算SSE

errors = y - yhat_2 # absolute errors for each data point

squared_errors = errors ** 2 # squared errors

sse2 = squared_errors。sum（） # sum of squared errors

sse2

您可以看到，擬合值（SSE = 52753）比我的猜測（SSE = 303857）更接近資料……這很好。

Statsmodels formula inputs

import statsmodels。formula。api as smf

注意，我們也可以透過使用statsmodels公式輸入

將Patsy公式直接傳遞給Statsmodels

來做同樣的事情：

fit_2 = smf。glm（formula=’TotalSleep ~ log_BodyWt‘， data=sleep）。fit（）

print（fit_2。summary（））

Graphically exploring model fit

擬合模型時，最好檢視

模型擬合的圖形表示形式以及殘差（預測值與觀察值之間的距離）以診斷問題

。 Statsmodels中提供了許多模型擬合圖，至少對於線性模型。

我們查看了模型擬合與上面資料之間關係的一種表示形式。 這是一些常見的情節。

# first we add to the sleep data frame， make a data frame of the stuff we want to plot：

plots_df = sleep［［’TotalSleep‘， ’log_BodyWt‘］］。copy（）

plots_df［’yhat‘］ = fit。predict（） # 預測值

plots_df［’resids‘］ = plots_df［’TotalSleep‘］ - plots_df［’yhat‘］ # 殘差

# note the fit object also contains various residuals。

plots_df。head（）

# plot predictions against real data：

plt。scatter（sleep［’log_BodyWt‘］， sleep［’TotalSleep‘］， c=’k‘）

plt。plot（sleep［’log_BodyWt‘］， fit。predict（）， c=’r‘）

plt。xlabel（’Log Bodyweight‘）

plt。ylabel（’Total Sleep‘）

plt。legend（［’Prediction‘， ’Data‘］）

plt。show（）

殘差圖

將殘差作為我們的預測變數的函式來看，可能會發現建模假設存在一些問題

。 回想一下，我們的模型期望觀察結果是獨立的，並且誤差呈正態分佈。 我們可以用一些圖解來圖形檢查。

# plot residuals against bodyweight：

resids = sleep［’TotalSleep‘］ - fit。predict（）

plt。scatter（plots_df［’log_BodyWt‘］， resids， c=’k‘）

plt。plot（（-8， 10）， （0， 0）， ’r——‘）

plt。xlabel（’Log Bodyweight‘）

plt。ylabel（’Residuals‘）

plt。tight_layout（）

plt。show（）

基本上，我們希望殘差在零線上的散佈看起來是隨機的。 那裡不應該有系統的東西

如果是這樣，我們的模型將丟失資料中的非隨機（即可預測）關係。

模型驗證：解釋殘差圖

如果殘差圖看起來像（b）–（f），請考慮擴充套件模型以捕獲資料中的其他關係。

殘差的分佈也應大致為正態：

plt。hist（resids。values）；

在這種情況下，它看起來有點右偏。 也就是說，我們可能違反了殘差的正態性假設，但可能並沒有大刀闊斧。

Seaborn regression plots

方便的是，Seaborn允許您直接進行其中一些繪圖，而無需首先進行擬合的statsmodels。 這對於快速探索資料中的關係非常有用，並且令人高興的是，這些圖也比預設的Matplotlib圖更漂亮。 有關示例和程式碼，請參見此處。

sns。lmplot（’log_BodyWt‘， ’TotalSleep‘， data=sleep）

sns。despine（offset=10， trim=True）；

sns。residplot（’log_BodyWt‘， ’TotalSleep‘， data=sleep）

sns。despine（offset=10， trim=True）；

在內部，這些函式使用statsmodels擬合與我們上面編寫的迴歸模型相同的模型，並將其包裝為漂亮的探索性圖。

分類變數和偽編碼

一些感興趣的變數將是

分類變數

，而不是度量標準變數/連續變數。 例如：

eg：

性別通常被認為是二分法（男性，女性）

一個人的宗教可能是以下幾種類別中的一種（基督教，穆斯林，佛教，非信徒。。。）

“種類”是我們上面的資料集中的一個示例

在設計矩陣中如何處理分類變數？

將類別變數轉換為一個或多個設計矩陣列，這些列具有有關變數級別的特定含義。 我們在這裡探討的最常用的方法稱為處理或偽編碼。 您可以在此處檢視Patsy內建的其他編碼方案的一些示例。

patsy.Treatment（reference = None）

為了演示的目的，讓我們建立一個類別變數“ danger_cat”，在其中我們考慮具有高/低掠食/暴露危險的動物（

high

or

low

danger of predation / exposure。）。

sleep［’danger_cat‘］ = ’low‘ # set as a string， “low” for all。

sleep。loc［sleep［’Danger‘］ > 3， ’danger_cat‘］ = ’high‘ # set Danger=4 or 5 to “high”

sleep［’danger_cat‘］

# how does total sleep appear to change according to our new categories？

sns。set_context（’notebook‘）

g = sns。boxplot（sleep［’TotalSleep‘］， groupby=sleep［’danger_cat‘］）

sns。despine（ax=g， offset=10）；

# Note how “violinplot” can be a bit misleading here， due to few cases and the kernel density smooth：

# g = sns。violinplot（sleep［’TotalSleep‘］， groupby=sleep［’danger_cat‘］）

# sns。despine（ax=g， offset=10）；

plt。savefig（’blah。pdf‘）

低危群體的睡眠時間似乎比高危群體的睡眠時間更長。

現在想象我們要使用GLM來測試這種現象。 高危群體的睡眠總量是否少於低微人群的睡眠量？

y， X = patsy。dmatrices（’TotalSleep ~ danger_cat‘， sleep）

X

現在，“ danger_cat”的設計矩陣列具有兩個值：0（高危險）和1（低危險）。 這對預測意味著什麼？ 讓我們考慮一下將這些結果乘以happens會發生什麼。

如果“ danger_cat” = 1，

̂ = 1 ∗ 1 + 2 ∗ 1

當“ danger_cat” = 0時，

̂ = 1 ∗ 1 + 2 ∗ 0

即，

截距係數 1表示高危險組的平均值（danger_cat = 0），

2係數表示高危險組的平均值與低危險組的平均值之差(低-高)。

讓我們看看這些係數是多少：

fit = smf。glm（formula=’TotalSleep ~ danger_cat‘， data=sleep）。fit（）

print（fit。summary（））

攔截係數為7。29小時睡眠。 從上方的方框圖中可以看出，這對應於高危險組的平均值。

高低組之間的差異（danger_cat的係數）為5.02小時

。 也就是說，低危險組的平均值為7。29 + 5。02 = 12。31小時（再次，請根據方框圖檢查它是否看起來正確）。

從z分數中我們可以看到，給定樣本的方差（std err），截距顯著不同於零（即高危險組的睡眠時間明顯超過零小時），而高危險組和低危險組之間的差異 這也是很重要的（即，低危人群的睡眠時間比高危人群的睡眠時間長得多）。

我們剛剛檢驗了兩個假設。 這些測試使您想起什麼？

t檢驗

測試“ danger_cat”變數的係數是否等於零（幾乎）與獨立樣本t檢驗相同！ 這並不完全相同，因為glm使用了z檢驗-如果我們用普通的最小二乘［smf。ols］擬合模型，則在這種情況下將完全相同。

1。 獨立樣本T檢驗僅用於分析兩組資料，比如性別，高分組低分組，實驗組對照組等。如超過兩組比較則使用方差分析

2。 獨立樣本T檢驗用於分析定類資料（X）與定量資料（Y）之間的差異情況。如果X、Y均為定類資料，則使用卡方分析

為了向我們證明這些是等效的，讓我們進行適當的測試：

high_array = sleep［’TotalSleep‘］［sleep［’danger_cat‘］==’high‘］

low_array = sleep［’TotalSleep‘］［sleep［’danger_cat‘］==’low‘］

# how does this compare to an independent groups z-test？

t， p = sm_stats。weightstats。ztest（high_array， low_array）

’Independent samples z-test： z = {t}， p = {p}‘。format（t=np。round（t， decimals=3）， p=p）

’Model values： z = {t}， p = {p}‘。format（

t=np。round（fit。tvalues［1］， decimals=3），

p=fit。pvalues［1］）

fit_ols = smf。ols（formula=’TotalSleep ~ danger_cat‘， data=sleep）。fit（）

print（fit_ols。summary（））

One-sample test （mean of

high

different to zero）

截距迴歸係數的檢驗在概念上與一樣本t檢驗相同：我們詢問“高”組的均值是否不同於零。 但是，實際上，由於迴歸的標準誤差是使用整個設計矩陣計算的，因此截距的檢驗統計量和置信區間與一樣本t檢驗不同。 鬆散地，截距的標準誤差取決於資料集中的總方差（即迴歸中包括的所有變數），而在單獨的一樣本測試中則不然。

# using statsmodels， must first compute “descrstats”：

d1 = sm_stats。weightstats。DescrStatsW（high_array）

t， p = d1。ztest_mean（0）

’Independent z-test： z = {t}， p = {p}‘。format（

t=np。round（t， decimals=3）， p=p）

’Model values： z = {t}， p = {p}‘。format（

t=np。round（fit。tvalues［0］， decimals=3），

p=fit。pvalues［0］）

# cell to check “high” alone

high_only = sleep。copy（）

high_only = high_only［high_only［’danger_cat‘］==’high‘］

fit_high = smf。glm（formula=’TotalSleep ~ 1‘， data=high_only）。fit（）

print（fit_high。summary（））

Categorical variable with more than one level

上面我們瞭解到，“虛擬編碼”是指在設計矩陣中將類別變數轉換為0和1時的情況。 這對於具有兩個以上級別的分類變數如何起作用？

讓我們嘗試將“壽命”變數（連續的，以年為單位的壽命）視為分類變數而不是連續變數。 我將建立一個新變數，將壽命分為三類：短，中，長。

為此，我們將使用Pandas的qcut函式，該函式根據分位數將變數切成組。 這將建立一個新的Pandas系列，其型別是分類的。 雖然類別在將來會非常有用，但現在它們與其他軟體包（例如statsmodels）的配合不好。 因此，我們將類別資料型別轉換為普通的舊物件。 請參閱文件以獲取類別。

Let’s work through what this means for the Beta weights：

For rows where 'life_cat' is 'long'

， the value of the two X columns is

［0， 0］

。 Therefore， the prediction for TotalSleep in these rows is：

̂ = 1∗1+ 2∗0+ 3∗0

̂ = 1

The first beta weight （the intercept） is therefore the mean of the ‘long’ group。

For rows where 'life_cat' is 'short'

， the value of the two X columns is

［1， 0］

。 Therefore， the prediction for TotalSleep in these rows is：

̂ = 1∗1+ 2∗1+ 3∗0

̂ = 1+ 2

The second beta weight is therefore the average

difference

between the ‘long’ group and the ‘short’ group。

Finally，

for rows where 'life_cat' is 'med'

， the value of the two X columns is

［0， 1］

。 Therefore， the prediction for TotalSleep in these rows is：

̂ = 1∗1+ 2∗0+ 3∗1

̂ = 1+ 3

The third beta weight is therefore the average

difference

between the ‘long’ group and the ‘med’ group。

Let‘s see what these coefficients turn out to be， by fitting the model：

GLM（General Linear Module ）廣義線性模型：連線函式解析