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稳健估计那样直接使用
的公式,而需要额外的条件。