MADlib——基於SQL的資料探勘解決方案（17）——迴歸之Cox比例風險迴歸

阿新 • • 發佈：2018-12-31

一、Cox比例風險迴歸簡介

Cox比例風險迴歸模型（Cox’s proportional hazards regression model），簡稱Cox迴歸模型，由英國統計學家D.R.Cox於1972年提出，主要用於腫瘤和其它慢性病的預後分析，也可用於佇列研究的病因探索。

1. 基本概念

生存函式：又稱累計生存率，簡稱生存率，常用S(t)表示。代表被觀察物件的生存時間大於t時刻的概率，實際中t時刻的取值估算公式為：S(t) ≈ 生存時間長於t的患者數 / 患者總數
死亡概率函式：簡稱為死亡概率，常用F(t)表示。代表一個被觀察物件從開始觀察到時間t為止的死亡概率，它與生存函式的關係為：F(t) = 1 - S(t)。

死亡密度函式：簡稱為密度函式，用f(t)表示。代表所有被觀察物件在t時刻的瞬時死亡率。
風險函式：即生存時間已達到t的一群被觀察物件在t時刻的瞬時死亡率，用h(t)表示。代表已存活到時間t的每個觀察物件從t到t+∆t這一非常小的區間內死亡的概率極限，它與生存函式、死亡密度函式的關係為：h(t) = f(t) / S(t)。

2. Cox迴歸模型結構

Cox迴歸模型不直接考察生存函式與協變數（影響因素）的關係，而是用風險函式作為因變數。設有n名病人（i=1,2,...,n），第i名病人的生存時間為，同時該病人具有一組協變數，則模型為：

：在時間t處的風險函式。

：基準風險函式，為所有協變數取零時t時刻的風險函式，即沒有協變數下的風險函式。這是模型中的非引數部分，因此Cox迴歸是一種半引數分析方法。

：協變數。

：根據觀察值估算出的迴歸係數。

：預後指數。大於0時表示該病人對應的危險度大於平均水平；等於0時為達到平均水平；小於0時表示該病人的危險度小於平均水平。

迴歸係數時，協變數的取值越大，風險函式的值越大，表示病人死亡的風險越大。
迴歸係數時，表示協變數對風險函式沒有影響。
迴歸係數時，協變數的取值越大，風險函式的值越小，表示病人死亡的風險越小。

3. Cox迴歸模型的前提條件

Cox迴歸模型必須滿足比例風險假設（Proportional Hazards Assumption，PHA）：

（1）任何兩個個體的風險函式之比，即風險比（HazardRatio，HR）保持一個恆定的比例，而與時間t無關。

（2）模型中協變數的效應不隨時間改變而改變。

檢查某協變數是否滿足PHA，最簡單的方法是觀察該變數分組的生存曲線。若生存曲線交叉，表示不滿足PHA，此時可採用分層比例風險模型。

4. 引數估算與假設檢驗

Cox迴歸的引數估計同邏輯迴歸分析一樣採用最大似然估計法。其基本思想是先建立偏似然函式或對數偏似然函式，求偏似然函式或對數偏似然函式達到極大時引數的取值，即為引數的最大似然估計值。

假設檢驗的方法有時協變數法、線性相關檢驗法、加權殘差Score法等。這三種檢驗法有較高的準確率，且三種方法的檢驗效能相近。MADlib的Cox模型PHA檢驗函式使用線性相關檢驗法實現。

5. Cox模型的注意事項

研究的協變數在被研究物件中的分佈要適中，否則會給迴歸引數的估計帶來困難。
Cox模型應用較靈活，被觀察物件進入研究佇列的早晚、時間長短可以不一致，但如果研究的變數隨時間而變化，可以採用時依協變數模型進行分析。
樣本含量不宜過小，一般應在40例以上，經驗上的做法是樣本含量為協變數個數的5-20倍。
儘管可以分析刪失資料，但在觀察時，要儘量避免觀察物件的失訪，過多失訪容易造成結果的偏倚。
Cox模型對異常值較為敏感，所以在進行模型擬合時要注意擬合優度的檢驗。

二、MADlib中Cox比例風險迴歸相關函式

1. 訓練函式

（1）語法

coxph_train( source_table,  
             output_table,  
             dependent_variable,  
             independent_variable,  
             right_censoring_status,  
             strata,  
             optimizer_params  
           )

（2）引數

引數名稱	資料型別	描述
source_table	VARCHAR	包含訓練資料的源表名。
output_table	VARCHAR	儲存模型的輸出表名，主輸出表列和概要輸出表列分別如表2、表3所示。
dependent_variable	VARCHAR	因變數名稱，指死亡時間，不需要對資料進行預排序。
independent_variable	VARCHAR	自變數名稱陣列。
right_censoring_status（可選）	VARCHAR	預設值為TRUE。表示右刪失狀態的字串，如果無刪失則為TRUE，否則為FALSE。該引數可以包含是右刪失狀態的列名，或者是一個可以對每個觀察值進行評估的布林表示式，如‘column_name < 10’。
strata（可選）	VARCHAR	預設值為NULL，不做任何分層，可以是逗號分隔的列名。
optimizer_params（可選）	VARCHAR	預設值為NULL，此時使用預設的優化引數：max_iter=100, optimizer=newton, tolerance=1e-8, array_agg_size=10000000, sample_size=1000000。應該是一個逗號分隔的‘key=value’字串，引數含義如下： l max_iter：最大迭代數。如果迭代次數達到此數值則停止計算，這通常意味著無法收斂。 l optimizer：優化方法，當前僅支援‘newton’。 l tolerance：停止條件。當連續兩次迭代的對數似然值之差小於此引數，計算已經收斂並停止。 l array_agg_size：為了加速計算，將原始資料表切分成多個數據片，每片資料聚合成一個大行。計算處理時，整個大行被匯入記憶體以提高運算速度。此引數控制一個大行中包含多少資料，引數值越大速度越快，但由於PostgreSQL資料庫的限制，一個大行的大小不能超過1G。 l sample_size：為了將資料切分成大致相等的片，先進行資料取樣，然後利用取樣資料找出斷點。該引數值越大，產生的斷點越精確。

表1 coxph_train函式引數說明

列名	資料型別	描述
Coef	FLOAT8[]	迴歸係數向量。
loglikelihood	FLOAT8	極大似然估計的對數似然值。
std_err	FLOAT8[]	迴歸係數標準差向量。
stats	FLOAT8[]	迴歸係數統計向量。
p_values	FLOAT8[]	迴歸係數p值向量。
hessian	FLOAT8[]	Hessian矩陣。
num_iterations	INTEGER	優化器執行的迭代次數。

表2 coxph_train函式主輸出表列說明

訓練函式在產生輸出表的同時，還會建立一個名為<output_table>_summary的概要表，具有以下列：

列名	資料型別	描述
method	VARCHAR	'coxph'，描述模型的字串。
source_table	VARCHAR	源表名。
out_table	VARCHAR	模型表名。
dependent_varname	VARCHAR	因變量表達式。
independent_varname	VARCHAR	自變量表達式。
right_censoring_status	VARCHAR	右刪失狀態。
Strata	VARCHAR	分層列。
num_processed	INTEGER	處理行數。
num_missing_rows_skipped	INTEGER	跳過行數。

表3 coxph_train函式概要輸出表列說明

2. 比例風險假設檢驗函式

cox_zph()函式檢驗Cox迴歸的比例風險假設，它通過計算coxph_train()輸出模型中殘差與時間的相關性驗證比例風險假設。

（1）語法

cox_zph(cox_model_table, output_table)

（2）引數

cox_model_table：TEXT型別，包含Cox模型的表名，應該是coxph_train()函式的輸出表。
output_table：TEXT型別，儲存測試結果的表名，主輸出表列和殘差輸出表列分別如表4、表5所示。

列名	資料型別	描述
covariate	TEXT	協變數。
rho	FLOAT8[]	生存時間與Schoenfeld殘差相關係數向量。
chi_square	FLOAT8[]	相關分析卡方檢驗統計量。
p_value	FLOAT8[]	卡方統計雙尾p值。

表4 cox_zph函式輸出表列說明

檢驗函式還建立一個儲存殘差值的輸出表，名為<output_table>_residual，具有以下列：

列名	資料型別	描述
<dep_column_name>	FLOAT8	原始源表中存在的時間值（因變數）。
residual	FLOAT8[]	原始協變數與coxph_train模型中的期望協變數之差。
scaled_residual	FLOAT8[]	由係數的方差來衡量的殘差值。

表5 cox_zph函式殘差輸出表列說明

3. 預測函式

（1）語法

coxph_predict(model_table,  
              source_table,  
              id_col_name,  
              output_table,  
              pred_type,  
              reference)

（2）引數

引數名稱	資料型別	描述
model_table	TEXT	包含cox模型的表名。
source_table	TEXT	包含預測資料的表名。
id_col_name	TEXT	id列名。
output_table	TEXT	儲存預測結果的輸出表名，輸出表具有以下列： l id：TEXT型別，id列。 l predicted_result：FLOAT8型別，基於預測型別引數值的預測結果。
pred_type（可選）	TEXT	預測型別。預設為‘linear_predictors’，可以是‘linear_predictors’、‘risk’或‘terms’。 l linear_predictors：計算自變數和係數的點積，也即預後指數。 l risk：預後指數預測結果的指數值。 l terms：每個協變數與它們相應的係數值的乘積，不求和。
reference（可選）	TEXT	用於中心預測。可以為‘strata’、‘overall’，預設值為‘strata’。

表6 coxph_predict函式引數說明

注：Cox迴歸模型的因變數是風險函式，因此與其它模型的預測函式不同，它不直接返回生存時間的預測值。

三、示例

1. 檢視聯機幫助

select madlib.coxph_train();  
-- 訓練函式用法  
select madlib.coxph_train('usage');  
-- 示例  
select madlib.coxph_train('example');

2. 建立測試資料表並裝載原始資料

drop table if exists sample_data;  
create table sample_data (  
    id integer not null,    -- 邏輯主鍵  
    grp double precision,   -- 分組（治療方法）  
    wbc double precision,   -- 白細胞值  
    timedeath integer,      -- 生存天數  
    status boolean          -- 結局  
);  
copy sample_data from stdin with delimiter '|';  
  1 |   0 | 1.45 |        35 | t  
  2 |   0 | 1.47 |        34 | t  
  3 |   0 |  2.2 |        32 | t  
  4 |   0 | 1.78 |        25 | t  
  5 |   0 | 2.57 |        23 | t  
  6 |   0 | 2.32 |        22 | t  
  7 |   0 | 2.01 |        20 | t  
  8 |   0 | 2.05 |        19 | t  
  9 |   0 | 2.16 |        17 | t  
 10 |   0 |  3.6 |        16 | t  
 11 |   1 |  2.3 |        15 | t  
 12 |   0 | 2.88 |        13 | t  
 13 |   1 |  1.5 |        12 | t  
 14 |   0 |  2.6 |        11 | t  
 15 |   0 |  2.7 |        10 | t  
 16 |   0 |  2.8 |         9 | t  
 17 |   1 | 2.32 |         8 | t  
 18 |   0 | 4.43 |         7 | t  
 19 |   0 | 2.31 |         6 | t  
 20 |   1 | 3.49 |         5 | t  
 21 |   1 | 2.42 |         4 | t  
 22 |   1 | 4.01 |         3 | t  
 23 |   1 | 4.91 |         2 | t  
 24 |   1 |    5 |         1 | t  
\.

說明：在這個假設的生存分析案例中，將24名患者分為兩組（如模擬兩種治療方法）進行觀察。協變數有兩個，分組與白細胞值，樣本量是協變數個數的12倍。因變數為生存天數。所有患者結局已知，不存在刪失情況。

3. 訓練模型

drop table if exists sample_cox, sample_cox_summary;  
select madlib.coxph_train( 'sample_data',  
                           'sample_cox',  
                           'timedeath',  
                           'array[grp,wbc]',  
                           'status'  
                         );

4. 查看回歸結果

\x on  
select * from sample_cox;

結果：

-[ RECORD 1 ]--+----------------------------------------------------------------------------
coef           | {2.54407073265254,1.67172094779487}
loglikelihood  | -37.8532498733
std_err        | {0.677180599294897,0.387195514577534}
z_stats        | {3.7568570855419,4.31751114064138}
p_values       | {0.000172060691513886,1.5779844638453e-05}
hessian        | {{2.78043065745617,-2.25848560642414},{-2.25848560642414,8.50472838284472}}
num_iterations | 5

經過5次迭代後求得兩個協變數的迴歸係數。z檢驗得到的p值很小，說明樣本間的差異由抽樣誤差所致的概率極小，具有顯著統計意義。

5. 檢驗所得模型的PHA

drop table if exists sample_zph_output, sample_zph_output_residual;  
select madlib.cox_zph( 'sample_cox',  
                       'sample_zph_output'  
                     );  
  
select * from sample_zph_output;

結果：

-[ RECORD 1 ]-----------------------------------------  
covariate  | array[grp,wbc]  
rho        | {0.00237308357328641,0.037560056884043}  
chi_square | {0.000100675718191977,0.0232317400546175}  
p_value    | {0.991994376850758,0.878855984657948}

madlib.cox_zph函式使用卡方線性相關法檢驗比例風險假設是否成立。它的原理很簡單：如果資料滿足比例風險假設，則Schoenfeld殘差和生存時間的秩次之間無相關性。計算步驟按以下3步進行：①用未刪失資料計算每個協變數Schoenfeld殘差；②將未刪失的生存時間排序，並以新變數（協變數殘差）記錄秩次1、2、3...，如出現相同生存時間（結點），則以平均秩次記錄。③進行假設檢驗，檢驗Schoenfeld殘差和生存時間秩次間的相關性（原假設為H0：r=0，無相關性）。若原假設被拒絕，則表明資料不滿足比例風險假設。

從本例的檢驗p值結果看，協變數對應的雙尾p值接近於1，說明應該接受原假設，模型滿足比例風險假設。

6. 用模型進行預測

本例使用源資料表演示預測。

（1）條目預測

\x off  
drop table if exists sample_pred_terms;  
select madlib.coxph_predict('sample_cox',  
                            'sample_data',  
                            'id',  
                            'sample_pred_terms',  
                            'terms');  
  
select id, terms, sum(c) linear_predictors  
  from (select id, terms, unnest(terms) c  
          from sample_pred_terms) t  
 group by id, terms  
 order by linear_predictors;

結果：

 id |                  terms                   | linear_predictors    
----+------------------------------------------+--------------------  
  1 | {-0.848023577550848,-2.12308560369949}   |  -2.97110918125034  
  2 | {-0.848023577550848,-2.08965118474359}   |  -2.93767476229444  
  4 | {-0.848023577550848,-1.57141769092718}   |  -2.41944126847803  
  7 | {-0.848023577550848,-1.18692187293436}   |  -2.03494545048521  
  8 | {-0.848023577550848,-1.12005303502256}   |  -1.96807661257341  
  9 | {-0.848023577550848,-0.936163730765128}  |  -1.78418730831598  
  3 | {-0.848023577550848,-0.869294892853333}  |  -1.71731847040418  
 19 | {-0.848023577550848,-0.685405588595898}  |  -1.53342916614675  
  6 | {-0.848023577550848,-0.668688379117949}  |   -1.5167119566688  
  5 | {-0.848023577550848,-0.250758142169231}  |  -1.09878171972008  
 14 | {-0.848023577550848,-0.200606513735385}  |  -1.04863009128623  
 15 | {-0.848023577550848,-0.0334344189558975} | -0.881457996506746  
 16 | {-0.848023577550848,0.133737675823589}   | -0.714285901727259  
 12 | {-0.848023577550848,0.267475351647179}   | -0.580548225903669  
 13 | {1.6960471551017,-2.03949955630974}      | -0.343452401208047  
 10 | {-0.848023577550848,1.47111443405949}    |  0.623090856508639  
 11 | {1.6960471551017,-0.702122798073847}     |   0.99392435702785  
 17 | {1.6960471551017,-0.668688379117949}     |   1.02735877598375  
 21 | {1.6960471551017,-0.501516284338462}     |   1.19453087076323  
 18 | {-0.848023577550848,2.85864282072923}    |   2.01061924317838  
 20 | {1.6960471551017,1.28722512980205}       |   2.98327228490375  
 22 | {1.6960471551017,2.15652002265538}       |   3.85256717775708  
 23 | {1.6960471551017,3.66106887567077}       |   5.35711603077247  
 24 | {1.6960471551017,3.81152376097231}       |     5.507570916074  
(24 rows)

條目預測的輸出是一個協變數陣列，每個元素值是對應協變數與迴歸係數的點積。

（2）預後指數預測

\x off  
drop table if exists sample_pred;  
select madlib.coxph_predict('sample_cox',  
                            'sample_data',  
                            'id',  
                            'sample_pred');  
  
select id, linear_predictors, exp(linear_predictors) risk  
  from sample_pred order by linear_predictors;

結果：

 id | linear_predictors  |        risk          
----+--------------------+--------------------  
  1 |  -2.97110918125034 | 0.0512464372091503  
  2 |  -2.93767476229444 |  0.052988797150351  
  4 |  -2.41944126847803 | 0.0889713146524469  
  7 |  -2.03494545048521 |  0.130687611255372  
  8 |  -1.96807661257341 |  0.139725343898993  
  9 |  -1.78418730831598 |  0.167933483703554  
  3 |  -1.71731847040418 |  0.179546963459175  
 19 |  -1.53342916614675 |  0.215794402223054  
  6 |   -1.5167119566688 |  0.219432204682558  
  5 |  -1.09878171972008 |   0.33327686110022  
 14 |  -1.04863009128623 |  0.350417460388247  
 15 | -0.881457996506746 |  0.414178600294289  
 16 | -0.714285901727259 |  0.489541567796517  
 12 | -0.580548225903669 |  0.559591499901242  
 13 | -0.343452401208048 |  0.709317242984782  
 10 |  0.623090856508639 |   1.86468261037507  
 11 |   0.99392435702785 |   2.70181658768287  
 17 |   1.02735877598375 |   2.79367735395697  
 21 |   1.19453087076323 |   3.30200833843654  
 18 |   2.01061924317838 |   7.46794038691976  
 20 |   2.98327228490375 |   19.7523463121038  
 22 |   3.85256717775708 |   47.1138577624205  
 23 |   5.35711603077247 |   212.112338046552  
 24 |     5.507570916074 |   246.551504798816  
(24 rows)

可以看到，1的風險度最小，24的風險度最大，並且預後指數預測的結果，就是對條目預測結果陣列的元素求和而得。

（3）風險預測

drop table if exists sample_pred_risk;  
select madlib.coxph_predict('sample_cox',  
                            'sample_data',  
                            'id',  
                            'sample_pred_risk',  
                            'risk');  
  
select * from sample_pred_risk order by risk;

結果：

 id |        risk          
----+--------------------  
  1 | 0.0512464372091503  
  2 |  0.052988797150351  
  4 | 0.0889713146524469  
  7 |  0.130687611255372  
  8 |  0.139725343898993  
  9 |  0.167933483703554  
  3 |  0.179546963459175  
 19 |  0.215794402223054  
  6 |  0.219432204682558  
  5 |   0.33327686110022  
 14 |  0.350417460388247  
 15 |  0.414178600294289  
 16 |  0.489541567796517  
 12 |  0.559591499901242  
 13 |  0.709317242984782  
 10 |   1.86468261037507  
 11 |   2.70181658768287  
 17 |   2.79367735395697  
 21 |   3.30200833843654  
 18 |   7.46794038691976  
 20 |   19.7523463121038  
 22 |   47.1138577624205  
 23 |   212.112338046552  
 24 |   246.551504798816  
(24 rows)

可以看到，風險預測就是預後指數結果的指數值。

MADlib——基於SQL的資料探勘解決方案（17）——迴歸之Cox比例風險迴歸

一、Cox比例風險迴歸簡介

1. 基本概念

2. Cox迴歸模型結構

3. Cox迴歸模型的前提條件

4. 引數估算與假設檢驗

5. Cox模型的注意事項

二、MADlib中Cox比例風險迴歸相關函式

1. 訓練函式

2. 比例風險假設檢驗函式

3. 預測函式

三、示例

1. 檢視聯機幫助

2. 建立測試資料表並裝載原始資料

3. 訓練模型

4. 查看回歸結果

5. 檢驗所得模型的PHA

6. 用模型進行預測

MADlib——基於SQL的資料探勘解決方案（17）——迴歸之Cox比例風險迴歸

MADlib——基於SQL的資料探勘解決方案（21）——分類之KNN

MADlib——基於SQL的資料探勘解決方案（8）——資料探索之描述性統計

MADlib——基於SQL的資料探勘解決方案（26）——聚類之k-means方法

Python資料探勘學習筆記（12）淘寶圖片爬蟲實戰

Python資料探勘學習筆記（1）文字挖掘入門

java.lang.IllegalArgumentException java.sql.Date.getHours的解決方案（SSH）

大資料常見錯誤解決方案（轉載）

資料探勘乾貨總結（六）--推薦演算法之CF

資料探勘乾貨總結（一）--NLP基礎

資料探勘文字分類（七）特徵提取

資料探勘文字分類（八）訓練和測試

資料探勘乾貨總結（五）--推薦演算法之CB

資料探勘文字分類（二）蒐集中文語料庫與ICTCLAS分詞

資料探勘乾貨總結（二）--NLP進階-詳解Jieba分詞工具

HAWQ + MADlib 玩轉資料探勘之（七）——關聯規則方法之Apriori演算法

HAWQ + MADlib 玩轉資料探勘之（六）——主成分分析與主成分投影

收藏 | 大資料應用及其解決方案（完整版）

社交網站的資料探勘與分析（高清版）PDF

機器學習&資料探勘筆記_16（常見面試之機器學習演算法思想簡單梳理）

MADlib——基於SQL的資料探勘解決方案（17）——迴歸之Cox比例風險迴歸

一、Cox比例風險迴歸簡介

1. 基本概念

2. Cox迴歸模型結構

3. Cox迴歸模型的前提條件

4. 引數估算與假設檢驗

5. Cox模型的注意事項

二、MADlib中Cox比例風險迴歸相關函式

1. 訓練函式

2. 比例風險假設檢驗函式

3. 預測函式

三、示例

1. 檢視聯機幫助

2. 建立測試資料表並裝載原始資料

3. 訓練模型

4. 查看回歸結果

5. 檢驗所得模型的PHA

6. 用模型進行預測

相關推薦