MADlib——基于SQL的数据挖掘解决方案（18）——回归之稳健方差
MADlib——基于SQL的数据挖掘解决方案（18）——回归之稳健方差
Robust Variance模块中的函数用于计算线性回归、逻辑回归、多类逻辑回归和Cox比例风险回归的稳健方差（Huber-White估计）。它们可用于计算具有潜在噪声异常值的数据集中数据的差异。此处实现的Huber-White与R模块“sandwich”中的“HC0”三明治操作完全相同。
线性、逻辑和多类逻辑回归的稳健方差接口是相似的。每种回归类型都有自己的训练函数。回归结果保存在一个输出表中，取决于回归类型，只具有很小的差异。
警告：请注意，与其它回归方法的接口不同，Cox比例风险的接口接受由coxph_train()函数产生的输出模型表。
稳健线性回归训练函数
robust_variance_linregr()函数语法如下：
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()函数语法如下：
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()函数语法如下：
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()函数语法如下：
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. 查看逻辑回归训练函数的联机帮助。
SELECT madlib.robust_variance_logregr();2. 创建训练数据表。
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. 运行逻辑回归训练函数并计算回归的稳健逻辑方差：
DROP TABLE IF EXISTS patients_logregr;
SELECT madlib.robust_variance_logregr( 'patients',
'patients_logregr',
'second_attack',
'ARRAY[1, treatment, trait_anxiety]'
);4. 查看回归结果。
\x on
SELECT * FROM patients_logregr;结果：
-[ 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}另外，结果中的数组可以更简单地输出。
\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比例风险训练函数的联机帮助。
SELECT madlib.robust_variance_coxph();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 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回归函数。
DROP TABLE IF EXISTS sample_cox, sample_cox_summary;
SELECT madlib.coxph_train( 'sample_data',
'sample_cox',
'timedeath',
'ARRAY[grp,wbc]',
'status'
);4. 运行稳健Cox回归函数。
DROP TABLE IF EXISTS sample_robust_cox;
SELECT madlib.robust_variance_coxph( 'sample_cox',
'sample_robust_cox'
);5. 查看稳健Cox回归结果。
\x on
SELECT * FROM sample_robust_cox;结果：
-[ 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稳健估计那样直接使用
的公式，而需要额外的条件。