MADlib——基于SQL的数据挖掘解决方案(18)——回归之稳健方差

2019-05-25 19:36:25 浏览数 (1)

Robust Variance模块中的函数用于计算线性回归、逻辑回归、多类逻辑回归和Cox比例风险回归的稳健方差(Huber-White估计)。它们可用于计算具有潜在噪声异常值的数据集中数据的差异。此处实现的Huber-White与R模块“sandwich”中的“HC0”三明治操作完全相同。

线性、逻辑和多类逻辑回归的稳健方差接口是相似的。每种回归类型都有自己的训练函数。回归结果保存在一个输出表中,取决于回归类型,只具有很小的差异。

警告:请注意,与其它回归方法的接口不同,Cox比例风险的接口接受由coxph_train()函数产生的输出模型表。

稳健线性回归训练函数

robust_variance_linregr()函数语法如下:

代码语言:javascript复制
robust_variance_linregr( source_table,
                         out_table,
                         dependent_varname,
                         independent_varname,
                         grouping_cols
                       )

source_table:VARCHAR类型,包含训练数据的表的名称。

out_table:VARCHAR类型,包含输出模型生成表的名称。输出表包含以下列:

  • coef:DOUBLE PRECISION[]类型,回归系数向量。
  • std_err:DOUBLE PRECISION[]类型,系数标准误向量。
  • t_stats:DOUBLE PRECISION[]类型,系数的t统计向量。
  • p_values:DOUBLE PRECISION[]类型,系数的p值向量。

还会创建一个名为<out_table>_summary的汇总表,它与linregr_train函数创建的汇总表相同。有关详细信息,请参阅线性回归的文档。

dependent_varname:VARCHAR类型,包含因变量的列的名称。

independent_varname:VARCHAR类型,用于评估自变量的表达式列表。截距变量不是假定的。通常在自变量列表中包含单个常数1项来提供明确的截距项。

grouping_cols(可选):VARCHAR类型,缺省值为NULL。一个表达式列表,用于将输入数据集分组为离散组,每组运行一次​​回归。当此值为NULL时,不使用分组,并生成单个结果模型。

稳健逻辑回归训练函数

robust_variance_logregr()函数语法如下:

代码语言:javascript复制
robust_variance_logregr( source_table,
                         out_table,
                         dependent_varname,
                         independent_varname,
                         grouping_cols,
                         max_iter,
                         optimizer,
                         tolerance,
                         verbose_mode
                       )

source_table:VARCHAR类型,包含训练数据的表的名称。

out_table:VARCHAR类型,包含输出模型生成表的名称。输出表包含以下列:

  • coef:DOUBLE PRECISION[]类型,回归系数向量。
  • std_err:DOUBLE PRECISION[]类型,系数标准误向量。
  • z_stats:DOUBLE PRECISION[]类型,系数的z统计向量。
  • p_values:DOUBLE PRECISION[]类型,系数的p值向量。

还会创建一个名为<out_table>_summary的汇总表,它与logregr_train函数创建的汇总表相同。有关详细信息,请参阅逻辑回归的文档。

dependent_varname:VARCHAR类型,包含因变量的列的名称。

independent_varname:VARCHAR类型,用于评估自变量的表达式列表。截距变量不是假定的。通常在自变量列表中包含单个常数1项来提供明确的截距项。

grouping_cols(可选):VARCHAR类型,缺省值为NULL。一个表达式列表,类似于SQL的“GROUP BY”子句,用于将输入数据集分组为离散组,每组运行一次​​回归。当此值为NULL时,不使用分组,并生成单个结果模型。

max_iter(可选):INTEGER类型,缺省值为20。允许的最大迭代次数。

optimizer:VARCHAR类型,缺省值为‘fista’。优化器名称,‘fista’或‘igd’。

tolerance(可选):DOUBLEPRECISION,缺省值为1e-6。结束迭代的标准。‘fista’和‘igd’优化器都计算两个连续迭代系数之间的平均差值,当差值小于tolerance或迭代次数大于max_iter时,计算停止。

verbose_mode(可选):BOOLEAN类型,缺省值为FALSE。回归拟合是否应该打印任何警告信息。

稳健多类逻辑回归函数

robust_variance_mlogregr()函数语法如下:

代码语言:javascript复制
robust_variance_mlogregr( source_table,
                          out_table,
                          dependent_varname,
                          independent_varname,
                          ref_category,
                          grouping_cols,
                          optimizer_params,
                          verbose_mode
                        )

source_table:VARCHAR类型,包含培训数据的表格的名称。

out_table:VARCHAR类型,存储回归模型的表的名称,具有以下列:

  • category:类别。
  • ref_category:用于建模的参考类别。
  • coef:回归系数向量。
  • std_err:系数标准误向量。
  • z_stats:系数的z统计向量。
  • p_values:系数的p值向量。

还会创建一个名为<out_table>_summary的汇总表,它与mlogregr_train函数创建的汇总表相同。有关详细信息,请参阅多类逻辑回归的文档。

dependent_varname:VARCHAR类型,包含因变量的列的名称。

independent_varname:VARCHAR类型,用于评估自变量的表达式列表。截距变量不是假定的。通常在自变量列表中包含单个常数1项来提供明确的截距项。independent_varname可以是包含数值数组的列的名称,也可以是格式为“ARRAY[1,x1,x2,x3]”的字符串,其中x1,x2和x3是列名。

ref_category(可选):INTEGER类型,参考类别,缺省值为0。

grouping_cols(可选):VARCHAR类型,缺省值为NULL。当前未实现,忽略任何非NULL值。一个表达式列表,类似于SQL的“GROUP BY”子句,用于将输入数据集分组为离散组,每组运行一次​​回归。当此值为NULL时,不使用分组,并生成单个结果模型。

optimizer_params(可选):TEXT类型,缺省值为NULL,使用优化器参数的缺省值:max_iter=20, optimizer='newton', tolerance=1e-4。它应该是一个包含由逗号分隔的‘key = value’对的字符串。

verbose_mode(可选):BOOLEAN类型,缺省值为FALSE。当前未实现。如果为TRUE,则打印回归拟合警告消息。

Cox比例风险稳健方差函数

robust_variance_coxph()函数语法如下:

代码语言:javascript复制
robust_variance_coxph(model_table, output_table)

参数:

model_table:TEXT类型,模型表的名称,与coxph_train()函数的“output_table”参数相同。

output_table:TEXT类型,保存输出的表的名称,有以下几列:

  • coef:FLOAT8[]类型,回归系数向量。
  • loglikelihood:FLOAT8类型,极大似然估计(Maximum Likelihood Estimate,MLE)的对数似然值。
  • std_err:FLOAT8[]类型,系数标准误向量。
  • robust_se:FLOAT8[]类型,系数的稳健标准误向量。
  • robust_z:FLOAT8[]类型,系数的稳健z统计向量。
  • robust_p:FLOAT8[]类型,系数的稳健p值向量。
  • hessian:FLOAT8[][]类型,黑塞矩阵(Hessian Matrix)。

示例

逻辑回归示例

1. 查看逻辑回归训练函数的联机帮助。

代码语言:javascript复制
SELECT madlib.robust_variance_logregr();

2. 创建训练数据表。

代码语言:javascript复制
DROP TABLE IF EXISTS patients;
CREATE TABLE patients (id INTEGER NOT NULL, second_attack INTEGER,
    treatment INTEGER, trait_anxiety INTEGER);
COPY patients FROM STDIN WITH DELIMITER '|';
  1 |             1 |         1 |            70
  3 |             1 |         1 |            50
  5 |             1 |         0 |            40
  7 |             1 |         0 |            75
  9 |             1 |         0 |            70
 11 |             0 |         1 |            65
 13 |             0 |         1 |            45
 15 |             0 |         1 |            40
 17 |             0 |         0 |            55
 19 |             0 |         0 |            50
  2 |             1 |         1 |            80
  4 |             1 |         0 |            60
  6 |             1 |         0 |            65
  8 |             1 |         0 |            80
 10 |             1 |         0 |            60
 12 |             0 |         1 |            50
 14 |             0 |         1 |            35
 16 |             0 |         1 |            50
 18 |             0 |         0 |            45
 20 |             0 |         0 |            60
.

3. 运行逻辑回归训练函数并计算回归的稳健逻辑方差:

代码语言:javascript复制
DROP TABLE IF EXISTS patients_logregr;
SELECT madlib.robust_variance_logregr( 'patients',
                                       'patients_logregr',
                                       'second_attack',
                                       'ARRAY[1, treatment, trait_anxiety]'
                                     );

4. 查看回归结果。

代码语言:javascript复制
x on
SELECT * FROM patients_logregr;

结果:

代码语言:javascript复制
-[ RECORD 1 ]-------------------------------------------------------
coef     | {-6.36346994178,-1.02410605239,0.119044916669}
std_err  | {3.4587206233261,1.171619257826,0.0534328864184024}
z_stats  | {-1.83983346294692,-0.874094587938307,2.22793348157963}
p_values | {0.0657926909731544,0.382066744588117,0.0258849510749644}

另外,结果中的数组可以更简单地输出。

代码语言:javascript复制
x off
SELECT unnest(array['intercept', 'treatment', 'trait_anxiety' ]) as attribute,
       unnest(coef) as coefficient,
       unnest(std_err) as standard_error,
       unnest(z_stats) as z_stat,
       unnest(p_values) as pvalue
FROM patients_logregr;

Cox比例风险示例

1. 查看稳健Cox比例风险训练函数的联机帮助。

代码语言:javascript复制
SELECT madlib.robust_variance_coxph();

2. 创建一个输入数据集。

代码语言:javascript复制
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 DELIMITER '|';
  0 |   0 | 1.45 |        35 | t
  1 |   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
.

3. 运行Cox回归函数。

代码语言:javascript复制
DROP TABLE IF EXISTS sample_cox, sample_cox_summary;
SELECT madlib.coxph_train( 'sample_data',
                           'sample_cox',
                           'timedeath',
                           'ARRAY[grp,wbc]',
                           'status'
                         );

4. 运行稳健Cox回归函数。

代码语言:javascript复制
DROP TABLE IF EXISTS sample_robust_cox;
SELECT madlib.robust_variance_coxph( 'sample_cox',
                           'sample_robust_cox'
                         );

5. 查看稳健Cox回归结果。

代码语言:javascript复制
x on
SELECT * FROM sample_robust_cox;

结果:

代码语言:javascript复制
-[ RECORD 1 ]- ----------------------------------------------------------------------------
coef          | {2.54407073265254,1.67172094779487}
loglikelihood | -37.8532498733
std_err       | {0.677180599294897,0.387195514577534}
robust_se     | {0.62109558107367,0.27477352143822}
robust_z      | {4.09610180812216,6.08399579058692}
robust_p      | {4.2016521207968e-05,1.17223683102586e-09}
hessian       | {{2.78043065745617,-2.25848560642414},{-2.25848560642414,8.50472838284472}}

技术背景

在进行回归分析时,我们有时对计算系数

的方差感兴趣。虽然内置的回归函数提供方差估计,但是我们可能更需要稳健方差估计。

稳健方差计算可以表示为三明治形式,即下面的形式:

是矩阵。

矩阵,也被称为面包,是相对直接的,可以计算为:

其中

是黑塞矩阵(Hessian Matrix)。

矩阵有几个变化,每个具有不同的稳健性属性。这里实现的是Huber-White三明治运算,它采取的形式为:

上述计算稳健方差的方法(Huber-White估计)用于线性回归、逻辑回归和多项式逻辑回归。在计算具有潜在噪声异常值的数据集中数据的差异时是很有用。此处实现的Huber-White等同于R模块“sandwich”中的“HC0”三明治操作。

在计算多类逻辑回归的稳健方差时,它使用默认参考类别零,并且回归系数被包括在输出表中。输出中的回归系数与多类逻辑回归函数的顺序相同。对于K个因变量(1,...,K)和J个类别(0,...,J-1)的问题,令

表示因变量k和类别j的系数。输出是

。该顺序与函数marginal_mlogregr的多类回归边际效应计算不一致。这是故意为之,因为所有多类回归(稳健、聚类、...)的接口将被移动到匹配边际效应使用的接口中。

Cox比例风险的稳健方差更复杂,因为系数是通过最大化部分对数似然来训练的。因此,不能像Huber-White稳健估计那样直接使用

的公式,而需要额外的条件。

0 人点赞