cox_prop_hazards.sql_in

Go to the documentation of this file.

/* ----------------------------------------------------------------------- *//**
 *
 * @file cox_prop_hazards.sql_in
 *
 * @brief SQL functions for cox proportional hazards
 * @date July 2012
 *
 * @sa For a brief introduction to cox regression, see the
 *     module description \ref grp_cox_prop_hazards
 *
 *//* ----------------------------------------------------------------------- */

m4_include(`SQLCommon.m4')

/**
@addtogroup grp_cox_prop_hazards

@about
Proportional-Hazard models enable the comparison of various survival models.
These survival models are functions describing the probability of an one-item
event (prototypically, this event is death) with respect to time.
The interval of time before death occurs is the survival time.
Let T be a random variable representing the survival time,
with a cumulative probability function P(t). Informally, P(t) is
the probability that death has happened before time t.

Generally, applications start with a list of \f$ \boldsymbol n \f$ observations,
each with \f$ \boldsymbol m \f$ covariates and a time of death. From this
\f$ \boldsymbol n \times m \f$ matrix, we would like to derive the correlation
between the covariates and the hazard function. This amounts to finding
the parameters \f$ \boldsymbol \beta \f$ that best fit the model described below.

Let us define:
- \f$ \boldsymbol t \in  \mathbf R^{m} \f$ denote the vector of observed dependent
  variables, with \f$ n \f$ rows.
- \f$ X \in \mathbf R^{m} \f$ denote the design matrix with \f$ m \f$
  columns and \f$ n \f$ rows, containing all observed vectors of independent
  variables \f$ \boldsymbol x_i \f$ as rows.
- \f$ R(t_i) \f$ denote the set of observations still alive at time \f$ t_i \f$

Note that this model <b>does not</b> include a <b>constant term</b>, and the data
cannot contain a column of 1s.

By definition,
\f[
    P[T_k = t_i | \boldsymbol R(t_i)]
    =    \frac{e^{\beta^T x_k} }{ \sum_{j \in R(t_i)} e^{\beta^T x_j}}.
    \,.
\f]

The <b>partial likelihood </b>function can now be generated as the product of
conditional probabilities:
\f[
\mathcal L = \prod_{i = 1}^n
                \left(
                \frac{e^{\beta^T x_i}}{ \sum_{j \in R(t_i)} e^{\beta^T x_j}}
                \right).
\f]

The log-likelihood form of this equation is
\f[
L = \sum_{i = 1}^n
                    \left[  \beta^T x_i
                                    - \log\left(\sum_{j \in R(t_i)} e^{\beta^T x_j }\right)
                    \right].
\f]

Using this score function and Hessian matrix, the partial likelihood can be
maximized using the <b> Newton-Raphson algorithm </b>.<b> Breslow's method </b>
is used to resolved tied times of deaths. The time of death for two records are
considered "equal" if they differ by less than 1.0e-6

The inverse of the Hessian matrix, evaluated at the estimate of
\f$ \boldsymbol \beta \f$, can be used as an <b>approximate variance-covariance
matrix </b> for the estimate, and used to produce approximate
<b>standard errors</b> for the regression coefficients.

\f[
    \mathit{se}(c_i) = \left( (H)^{-1} \right)_{ii}
    \,.
\f]
The Wald z-statistic is
\f[
    z_i = \frac{c_i}{\mathit{se}(c_i)}
    \,.
\f]

The Wald \f$ p \f$-value for coefficient \f$ i \f$ gives the probability (under
the assumptions inherent in the Wald test) of seeing a value at least as extreme
as the one observed, provided that the null hypothesis (\f$ c_i = 0 \f$) is
true. Letting \f$ F \f$ denote the cumulative density function of a standard
normal distribution, the Wald \f$ p \f$-value for coefficient \f$ i \f$ is
therefore
\f[
    p_i = \Pr(|Z| \geq |z_i|) = 2 \cdot (1 - F( |z_i| ))
\f]
where \f$ Z \f$ is a standard normally distributed random variable.


The condition number is computed as \f$ \kappa(H) \f$ during the iteration
immediately <em>preceding</em> convergence (i.e., \f$ A \f$ is computed using
the coefficients of the previous iteration). A large condition number (say, more
than 1000) indicates the presence of significant multicollinearity.


@input

The training data is expected to be of the following form:\n
<pre>{TABLE|VIEW} <em>sourceName</em> (
    <em>inputTable</em> VARCHAR,
    <em>outputTable</em> VARCHAR,
    <em>dependentVariable</em> VARCHAR,
    <em>independentVariable</em> VARCHAR,
    [<em>rightCensoringStatus</em> VARCHAR]
)</pre>
Note: Dependent Variables refer to the time of death. There is no need to
pre-sort the data.

NOTE2:'right_censoring_status' is set to TRUE to if the observation is not censored
and 'FALSE' if the observation is censored. The default value for
'right_censoring_status' is TRUE for all observatoions
@usage


<b> The Full Interface</b>

<pre>
SELECT madlib.\ref cox_prop_hazards(
    <em>'source_table'</em>,            -- name of input table, VARCHAR
    <em>'out_table'</em>,               -- name of output table, VARCHAR
    <em>'dependent_varname'</em>,       -- name of dependent variable, VARCHAR
    <em>'independent_varname'</em>,     -- name of independent variable, VARCHAR
    [<em>'right_censoring_status'</em>, -- name of column with right censoring status, VARCHAR (OPTIONAL, default=True)

);
</pre>

Here the <em>'right_censoring_status'</em> can be the name of a column, which contains
array of boolean values. It can also have a format of string 'dependent_variable < 10',
where <em>x1</em>, <em>x2</em> and <em>x3</em> are all column names.

Here the <em>'independent_varname'</em> can be the name of a column, which contains
array of numeric values. It can also have a format of string 'array[1, x1, x2, x3]',
where <em>x1</em>, <em>x2</em> and <em>x3</em> are all column names.


Output is stored in the <em>out_table</em>:
<pre>
[ coef | std_err | stats | p_values |
+------+---------+-------+----------+
</pre>


-#  For function summary information. Run
@verbatim
sql> select cox_prop_hazards('help');
OR
sql> select cox_prop_hazards();
OR
sql> select cox_prop_hazards('?');
@endverbatim

-#  For function usage information.
@verbatim
sql> select cox_prop_hazards('usage');
@endverbatim


Note: The function cox_prop_hazards_regr has been deprecated but maintained
@examp

-# Create the sample data set:
@verbatim
sql> SELECT * FROM data;
      val   | time | status
------------|--------------
 {0,1.95}   |  35  |  t
 {0,2.20}   |  28  |  t
 {1,1.45}   |  32  |  t
 {1,5.25}   |  31  |  t
 {1,0.38}   |  21  |  t
...
@endverbatim
-# Run the cox regression function:
@verbatim
sql> SELECT * FROM cox_prop_hazards('data', 'result_table', 'val', 'time', 'status');
sql> SELECT * from result_table;
--------------|--------------------------------------------------------------
coef           | {0.881089349817059,-0.0756817768938055}
std_err        | {1.16954914708414,0.338426252282655}
z_stats        | {0.753356711368689,-0.223628410729811}
p_values       | {0.451235588326831,0.823046454908087}

@endverbatim

@literature

A somewhat random selection of nice write-ups, with valuable pointers into
further literature:

[1] John Fox: Cox Proportional-Hazards Regression for Survival Data,
        Appendix to An R and S-PLUS companion to Applied Regression Feb 2012,
    http://cran.r-project.org/doc/contrib/Fox-Companion/appendix-cox-regression.pdf

[2] Stephen J Walters: What is a Cox model?
    http://www.medicine.ox.ac.uk/bandolier/painres/download/whatis/cox_model.pdf


@note Source and column names have to be passed as strings (due to limitations
of the SQL syntax).


@sa File cox_prop_hazards.sql_in (documenting the SQL functions)

@internal
@sa Namespace cox_prop_hazards
    \ref madlib::modules::stats documenting the implementation in C++
@endinternal

*/


DROP TYPE IF EXISTS MADLIB_SCHEMA.intermediate_cox_prop_hazards_result;
CREATE TYPE MADLIB_SCHEMA.intermediate_cox_prop_hazards_result AS (
    x DOUBLE PRECISION[],
    status BOOLEAN,
    exp_coef_x DOUBLE PRECISION,
    x_exp_coef_x DOUBLE PRECISION[],
    x_xTrans_exp_coef_x DOUBLE PRECISION[]
);


CREATE OR REPLACE FUNCTION MADLIB_SCHEMA.intermediate_cox_prop_hazards(
    /*+  x */ DOUBLE PRECISION[],
    /*+  status */ BOOLEAN,
    /*+  coef */ DOUBLE PRECISION[])
RETURNS MADLIB_SCHEMA.intermediate_cox_prop_hazards_result AS
'MODULE_PATHNAME'
LANGUAGE c IMMUTABLE;



DROP TYPE IF EXISTS MADLIB_SCHEMA.cox_prop_hazards_result;
CREATE TYPE MADLIB_SCHEMA.cox_prop_hazards_result AS (
    coef DOUBLE PRECISION[],
    logLikelihood DOUBLE PRECISION,
    std_err DOUBLE PRECISION[],
    z_stats DOUBLE PRECISION[],
    p_values DOUBLE PRECISION[],
    num_iterations INTEGER
);




CREATE OR REPLACE FUNCTION MADLIB_SCHEMA.cox_prop_hazards_step_final(
    state DOUBLE PRECISION[])
RETURNS DOUBLE PRECISION[]
AS 'MODULE_PATHNAME'
LANGUAGE C IMMUTABLE STRICT;

CREATE OR REPLACE FUNCTION MADLIB_SCHEMA.cox_prop_hazards_step_transition(
    /*+  state */ DOUBLE PRECISION[],
    /*+  x */ DOUBLE PRECISION[],
    /*+  y */ DOUBLE PRECISION,
    /*+  status */ BOOLEAN,
    /*+  exp_coef_x */ DOUBLE PRECISION,
    /*+  xexp_coef_x */ DOUBLE PRECISION[],
    /*+  x_xTrans_exp_coef_x */ DOUBLE PRECISION[],
    /*+  previous_state */ DOUBLE PRECISION[])
RETURNS DOUBLE PRECISION[] AS
'MODULE_PATHNAME'
LANGUAGE C IMMUTABLE;


/**
 * @internal
 * @brief Perform one iteration the Newton-Rhapson method.
 */
CREATE
m4_ifdef(`__GREENPLUM__',m4_ifdef(`__HAS_ORDERED_AGGREGATES__',`ORDERED'))
AGGREGATE MADLIB_SCHEMA.cox_prop_hazards_step(

    /*+  x */ DOUBLE PRECISION[],
    /*+  y */ DOUBLE PRECISION,
    /*+  status */ BOOLEAN,
    /*+  exp_coef_x */ DOUBLE PRECISION,
    /*+  xexp_coef_x */ DOUBLE PRECISION[],
    /*+  x_xTrans_exp_coef_x */ DOUBLE PRECISION[],
    /*+ previous_state */ DOUBLE PRECISION[]) (
    STYPE=DOUBLE PRECISION[],
    SFUNC=MADLIB_SCHEMA.cox_prop_hazards_step_transition,
    FINALFUNC=MADLIB_SCHEMA.cox_prop_hazards_step_final,
    INITCOND='{0,0,0,0,0,0,0}'
);



CREATE OR REPLACE FUNCTION MADLIB_SCHEMA.internal_cox_prop_hazards_step_distance(
    /*+ state1 */ DOUBLE PRECISION[],
    /*+ state2 */ DOUBLE PRECISION[])
RETURNS DOUBLE PRECISION AS
'MODULE_PATHNAME'
LANGUAGE c IMMUTABLE STRICT;

CREATE OR REPLACE FUNCTION MADLIB_SCHEMA.internal_cox_prop_hazards_result(
    /*+ state */ DOUBLE PRECISION[])
RETURNS MADLIB_SCHEMA.cox_prop_hazards_result AS
'MODULE_PATHNAME'
LANGUAGE c IMMUTABLE STRICT;


-- We only need to document the last one (unfortunately, in Greenplum we have to
-- use function overloading instead of default arguments).
CREATE FUNCTION MADLIB_SCHEMA.compute_cox_prop_hazards_regr(
    "source" VARCHAR,
    "indepColumn" VARCHAR,
    "depColumn" VARCHAR,
    "status" VARCHAR,
    "maxNumIterations" INTEGER,
    "optimizer" VARCHAR,
    "precision" DOUBLE PRECISION)
RETURNS INTEGER
AS $$PythonFunction(stats, cox_prop_hazards, compute_cox_prop_hazards)$$
LANGUAGE plpythonu VOLATILE;


-- Note: The function cox_prop_hazards_regr has been deprecated but maintained
/**
 * @brief Compute cox-regression coefficients and diagnostic statistics
 *
 * To include an intercept in the model, set one coordinate in the
 * <tt>independentVariables</tt> array to 1.
 *
 * @param source Name of the source relation containing the training data
 * @param indepColumn Name of the independent column
 * @param depColumn Name of the dependent column measuring time of death
 * @param status Name of the column that determines right censoring support
 * @param maxNumIterations The maximum number of iterations
 * @param optimizer The optimizer to use (either
 *        <tt>'newton'</tt>/<tt>'newton'</tt> for the newton method
 * @param precision The difference between log-likelihood values in successive
 *        iterations that should indicate convergence. Note that a non-positive
 *        value here disables the convergence criterion, and execution will only
 *        stop after \c maxNumIterations iterations.
 *
 * @return A composite value:
 *  - <tt>coef FLOAT8[]</tt> - Array of coefficients, \f$ \boldsymbol \beta \f$
 *  - <tt>log_likelihood FLOAT8</tt> - Log-likelihood \f$l(\boldsymbol \beta)\f$
 *  - <tt>std_err FLOAT8[]</tt> - Array of standard errors,
 *    \f$ \mathit{se}(c_1), \dots, \mathit{se}(c_k) \f$
 *  - <tt>z_stats FLOAT8[]</tt> - Array of Wald z-statistics, \f$ \boldsymbol z \f$
 *  - <tt>p_values FLOAT8[]</tt> - Array of Wald p-values, \f$ \boldsymbol p \f$
 *  - <tt>condition_no FLOAT8</tt> - The condition number of matrix
 *    \f$ H \f$ during the iteration immediately <em>preceding</em>
 *    convergence (i.e., \f$ H \f$ is computed using the coefficients of the
 *    previous iteration)
 *  - <tt>num_iterations INTEGER</tt> - The number of iterations before the
 *    algorithm terminated
 *
 * - Get vector of coefficients \f$ \boldsymbol \beta \f$ and all diagnostic
 *  statistics:\n
 *  <pre>SELECT * FROM \ref cox_prop_hazards_regr(
 *    '<em>sourceName</em>', '<em>dependentVariable</em>',
 *      '<em>independentVariables</em>'
 *    [, <em>numberOfIterations</em> [, '<em>optimizer</em>' [, <em>precision</em> ] ] ]
 * );</pre>
 * - Get vector of coefficients \f$ \boldsymbol \beta \f$:\n
 *  <pre>SELECT (\ref cox_prop_hazards_regr('<em>sourceName</em>',
 * '<em>dependentVariable</em>', '<em>independentVariables</em>')).coef;</pre>
 * - Get a subset of the output columns, e.g., only the array of coefficients
 *  \f$ \boldsymbol \beta \f$, the log-likelihood of determination:
 *  <pre>SELECT coef, log_likelihood
 * FROM \ref cox_prop_hazards_regr('<em>sourceName</em>', '<em>dependentVariable</em>',
 * '<em>independentVariables</em>');</pre>
 */
CREATE FUNCTION MADLIB_SCHEMA.cox_prop_hazards_regr(
    "source" VARCHAR,
    "indepColumn" VARCHAR,
    "depColumn" VARCHAR,
    "status" VARCHAR,
    "maxNumIterations" INTEGER /*+ DEFAULT 20 */,
    "optimizer" VARCHAR /*+ DEFAULT 'newton' */,
    "precision" DOUBLE PRECISION /*+ DEFAULT 0.0001 */)
RETURNS MADLIB_SCHEMA.cox_prop_hazards_result AS $$
DECLARE
    theIteration INTEGER;
    fnName VARCHAR;
    theResult MADLIB_SCHEMA.cox_prop_hazards_result;
BEGIN
    theIteration := (
        SELECT MADLIB_SCHEMA.compute_cox_prop_hazards_regr($1, $2, $3, $4, $5, $6, $7)
    );
    IF optimizer = 'newton' THEN
        fnName := 'internal_cox_prop_hazards_result';
    ELSE
        RAISE EXCEPTION 'Unknown optimizer (''%'')', optimizer;
    END IF;
    EXECUTE
        $sql$
        SELECT (result).*
        FROM (
            SELECT
                MADLIB_SCHEMA.$sql$ || fnName || $sql$(_madlib_state) AS result
                FROM _madlib_iterative_alg
                WHERE _madlib_iteration = $sql$ || theIteration || $sql$
            ) subq
        $sql$
        INTO theResult;

    -- The number of iterations are not updated in the C++ code. We do it here.
    IF NOT (theResult IS NULL) THEN
        theResult.num_iterations = theIteration;
    END IF;
    RETURN theResult;
END;
$$ LANGUAGE plpgsql VOLATILE;


-- Note: The function cox_prop_hazards_regr has been deprecated but maintained
CREATE FUNCTION MADLIB_SCHEMA.cox_prop_hazards_regr(
    "source" VARCHAR,
    "indepColumn" VARCHAR,
    "depColumn" VARCHAR,
    "status" VARCHAR)
RETURNS MADLIB_SCHEMA.cox_prop_hazards_result AS
$$SELECT MADLIB_SCHEMA.cox_prop_hazards_regr($1, $2, $3, $4, 20, 'newton', 0.0001);$$
LANGUAGE sql VOLATILE;

-- Note: The function cox_prop_hazards_regr has been deprecated but maintained
CREATE FUNCTION MADLIB_SCHEMA.cox_prop_hazards_regr(
    "source" VARCHAR,
    "indepColumn" VARCHAR,
    "depColumn" VARCHAR,
    "status" VARCHAR,
    "maxNumIterations" INTEGER)
RETURNS MADLIB_SCHEMA.cox_prop_hazards_result AS
$$SELECT MADLIB_SCHEMA.cox_prop_hazards_regr($1, $2, $3, $4, $5, 'newton', 0.0001);$$
LANGUAGE sql VOLATILE;

-- Note: The function cox_prop_hazards_regr has been deprecated but maintained
CREATE FUNCTION MADLIB_SCHEMA.cox_prop_hazards_regr(
    "source" VARCHAR,
    "indepColumn" VARCHAR,
    "depColumn" VARCHAR,
    "status" VARCHAR,
    "maxNumIterations" INTEGER,
    "optimizer" VARCHAR)
RETURNS MADLIB_SCHEMA.cox_prop_hazards_result AS
$$SELECT MADLIB_SCHEMA.cox_prop_hazards_regr($1, $2, $3, $4, $5, $6, 0.0001);$$
LANGUAGE sql VOLATILE;


/**
  * @brief Cox regresison training function
 **/
CREATE FUNCTION MADLIB_SCHEMA.cox_prop_hazards(
     source_table               VARCHAR       -- name of input  table
   , out_table                  VARCHAR       -- name of output table
   , dependent_varname          VARCHAR       -- name of dependent variable
   , independent_varname        VARCHAR       -- name of independent variable
   , status                     VARCHAR       -- censoring status
  )
RETURNS VOID AS $$
DECLARE
insert_string           VARCHAR;
cox_rst                 MADLIB_SCHEMA.cox_prop_hazards_result;
old_msg_level           TEXT;
BEGIN
  EXECUTE 'SELECT setting FROM pg_settings WHERE name=''client_min_messages''' INTO old_msg_level;
  EXECUTE 'SET client_min_messages TO warning';
  IF (source_table IS NULL OR source_table = '') THEN
      RAISE EXCEPTION 'Invalid input table name given.';
  END IF;
  IF (NOT MADLIB_SCHEMA.__table_exists(source_table)) THEN
      RAISE EXCEPTION 'Input table name does not exist.';
  END IF;
  IF (MADLIB_SCHEMA.__table_exists(out_table)) THEN
    RAISE EXCEPTION 'Output table name already exists. Drop the table before calling the function.';
  END IF;
  IF (out_table IS NULL OR out_table = '') THEN
      RAISE EXCEPTION 'Invalid output table name given.';
  END IF;
  IF (dependent_varname IS NULL OR dependent_varname = '') THEN
      RAISE EXCEPTION 'Invalid dependent variable name given.';
  END IF;
  IF (independent_varname IS NULL OR independent_varname = '') THEN
      RAISE EXCEPTION 'Invalid independent variable name given.';
  END IF;
  IF (status IS NULL OR independent_varname = '') THEN
      RAISE EXCEPTION 'Invalid status name given.';
  END IF;
  -- create output table with appropriate column names
  EXECUTE 'DROP TABLE IF EXISTS ' || out_table;
  EXECUTE '
          CREATE TABLE ' || out_table || ' (
              coef DOUBLE PRECISION[],
              std_err DOUBLE PRECISION[],
              z_stats DOUBLE PRECISION[],
              p_values DOUBLE PRECISION[])';
  -- compute cox regression results
  cox_rst := MADLIB_SCHEMA.__internal_get_cox_prop_hazards_result(
              source_table,  dependent_varname, independent_varname, status);
  insert_string := MADLIB_SCHEMA.__internal_get_cox_prop_hazards_insert_string(
                  cox_rst, out_table);
  -- Ensure Infinity and NaN are cast properly
  insert_string := REGEXP_REPLACE(insert_string, 'Infinity',
                                  '''Infinity''::double precision', 'gi');
  insert_string := REGEXP_REPLACE(insert_string, 'NaN',
                                  '''NaN''::double precision', 'gi');
  EXECUTE insert_string;
  EXECUTE 'SET client_min_messages TO '|| old_msg_level;
END;
$$ LANGUAGE plpgsql VOLATILE;



/**
  * @brief Return cox regression output for source data
  *
**/
CREATE FUNCTION MADLIB_SCHEMA.__internal_get_cox_prop_hazards_result(
     source_table         VARCHAR       -- name of input  table
   , dependent_varname    VARCHAR       -- name of dependent variable
   , independent_varname  VARCHAR       -- name of independent variable
   , status               VARCHAR       -- name of status variable
)
RETURNS MADLIB_SCHEMA.cox_prop_hazards_result AS $$
DECLARE
cox_rst                 MADLIB_SCHEMA.cox_prop_hazards_result;
BEGIN
  EXECUTE '
        SELECT (MADLIB_SCHEMA.cox_prop_hazards_regr('
                || '''' || source_table || ''', '
                || '''' || independent_varname    || ''', '
                || '''' || dependent_varname    || ''', '
                || '''' || status || ''')).* '
  INTO cox_rst;
  RETURN cox_rst;
END;
$$ LANGUAGE plpgsql VOLATILE;


/**
  * @brief Return cox regression output for insert string
  *
**/
CREATE FUNCTION MADLIB_SCHEMA.__internal_get_cox_prop_hazards_insert_string(
    cox_rst MADLIB_SCHEMA.cox_prop_hazards_result,
    out_table TEXT
)
RETURNS VARCHAR AS $$
DECLARE
  insert_string VARCHAR;
BEGIN
  insert_string := 'INSERT INTO ' || out_table || ' VALUES (';
  insert_string := insert_string  ||
  CASE
    WHEN (cox_rst).coef is NULL
    THEN '''{}'','
    ELSE 'ARRAY[' || array_to_string((cox_rst).coef, ',')     || '], '
  END             ||
  CASE
    WHEN (cox_rst).std_err is NULL
    THEN '''{}'','
    ELSE 'ARRAY[' || array_to_string((cox_rst).std_err, ',')  || '], '
  END             ||
  CASE
    WHEN (cox_rst).z_stats is NULL
    THEN '''{}'','
    ELSE 'ARRAY[' || array_to_string((cox_rst).z_stats, ',')  || '], '
  END             ||
  CASE
    WHEN (cox_rst).p_values is NULL
    THEN '''{}'','
    ELSE 'ARRAY[' || array_to_string((cox_rst).p_values, ',') || ']) '
  END;
  RETURN insert_string;
END;
$$ LANGUAGE plpgsql VOLATILE;


/**
  * @brief Cox regresison training function
 **/
CREATE FUNCTION MADLIB_SCHEMA.cox_prop_hazards(
     source_table               VARCHAR       -- name of input  table
   , out_table                  VARCHAR       -- name of output table
   , dependent_variable          VARCHAR       -- name of dependent variable
   , independent_variable        VARCHAR       -- name of independent variable
  )
RETURNS VOID AS $$
BEGIN
  EXECUTE 'SELECT MADLIB_SCHEMA.cox_prop_hazards(' ||
          ' ''' || source_table || ''' ' ||
          ' ,''' || out_table || ''' ' ||
          ' ,''' || dependent_variable || ''' ' ||
          ' ,''' || independent_variable || ''' ' ||
          ' , ''TRUE'')';
END;
$$ LANGUAGE plpgsql VOLATILE;


CREATE OR REPLACE FUNCTION MADLIB_SCHEMA.cox_prop_hazards()
RETURNS VARCHAR AS $$
BEGIN
  RETURN MADLIB_SCHEMA.cox_prop_hazards('');
END;
$$ LANGUAGE plpgsql VOLATILE;

CREATE OR REPLACE FUNCTION MADLIB_SCHEMA.cox_prop_hazards(
     usage_string VARCHAR                               -- usage string
)
RETURNS VARCHAR AS $$
DECLARE
insert_string       VARCHAR;
BEGIN
    IF (usage_string = '' OR usage_string = 'help' OR usage_string = '?') THEN
    insert_string := '' ||
    E'Summary \n' ||
    E'-----------------------------------------------------------------------------------------\n' ||
    E' Functionality: Cox proprtional hazards regression (Breslow method)\n' ||
    E' SELECT {schema_madlib}.cox_prop_hazards(''source_table''  \n' ||
    E'                                        ,''output_table''  \n' ||
    E'                                        ,''dependent_variable'' \n' ||
    E'                                        ,''independent_variable'' \n' ||
    E'                                        ,''right_censoring_status''  \n' ||
    E'                                        );\n' ||
    E'For more details on function usage:  \n' ||
    E'SELECT {schema_madlib}.cox_prop_hazards(''usage'') \n' ||
    E'';
  ElSIF (usage_string  = 'usage') THEN
    insert_string := '' ||
    E'Usage\n' ||
    E'-----------------------------------------------------------------------------------------\n' ||
    E' SELECT {schema_madlib}.cox_prop_hazards( \n' ||
    E' ''source_table'',            -- Name of data table          \n' ||
    E' ''output_table'',            -- Name of result table (overwrites if exists) \n' ||
    E' ''dependent_variable'',      -- Name of column for dependent variables\n' ||
    E' ''independent_variable'',    -- Name of column for independent variables\n' ||
    E'                                  (can be any SQL expression Eg: ''*'') \n' ||
    E' [''right_censoring_status'', -- Name of the column containing censoring status \n' ||
    E'                                    0/false : If the observation is censored\n' ||
    E'                                    1/true : otherwise\n' ||
    E'                                Default is 1/true for all observations\n' ||
    E'                                Can also be an SQL expression: ''dependent_variable < 10'') \n' ||
    E' );\n' ||
    E'\n'    ||
    E'Output:\n' ||
    E'-----------------------------------------------------------------------------------------\n' ||
    E' The output table (''output_table'' above) has the following columns\n' ||
    E' ''coef''    DOUBLE PRECISION[],  -- Coefficients of regression \n' ||
    E' ''std_err''  DOUBLE PRECISION[], -- Standard errors\n' ||
    E' ''z_stats''    DOUBLE PRECISION[], -- z-stats of the standard errors\n' ||
    E' ''p_values'' DOUBLE PRECISION[], -- p-values of the standard errors\n' ||
    E'\n' ||
    E'';
  ELSE
    insert_string := 'No such option. Run SELECT {schema_madlib}.cox_prop_hazards()';
  END IF;
  RETURN insert_string;
END;
$$ LANGUAGE plpgsql VOLATILE;