test/cox.cc

Code
Comments
Other
Rev Date Author Line
4198 19 Aug 22 peter 1 // $Id$
4198 19 Aug 22 peter 2
4198 19 Aug 22 peter 3 /*
4198 19 Aug 22 peter 4   Copyright (C) 2022 Peter Johansson
4198 19 Aug 22 peter 5
4198 19 Aug 22 peter 6   This file is part of the yat library, https://dev.thep.lu.se/yat
4198 19 Aug 22 peter 7
4198 19 Aug 22 peter 8   The yat library is free software; you can redistribute it and/or
4198 19 Aug 22 peter 9   modify it under the terms of the GNU General Public License as
4198 19 Aug 22 peter 10   published by the Free Software Foundation; either version 3 of the
4198 19 Aug 22 peter 11   License, or (at your option) any later version.
4198 19 Aug 22 peter 12
4198 19 Aug 22 peter 13   The yat library is distributed in the hope that it will be useful,
4198 19 Aug 22 peter 14   but WITHOUT ANY WARRANTY; without even the implied warranty of
4198 19 Aug 22 peter 15   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
4198 19 Aug 22 peter 16   General Public License for more details.
4198 19 Aug 22 peter 17
4198 19 Aug 22 peter 18   You should have received a copy of the GNU General Public License
4198 19 Aug 22 peter 19   along with yat. If not, see <https://www.gnu.org/licenses/>.
4198 19 Aug 22 peter 20 */
4198 19 Aug 22 peter 21
4198 19 Aug 22 peter 22 #include <config.h>
4198 19 Aug 22 peter 23
4198 19 Aug 22 peter 24 #include "Suite.h"
4198 19 Aug 22 peter 25
4198 19 Aug 22 peter 26 #include <yat/regression/Cox.h>
4198 19 Aug 22 peter 27 #include <yat/regression/MultiCox.h>
4198 19 Aug 22 peter 28
4198 19 Aug 22 peter 29 #include <yat/utility/Matrix.h>
4198 19 Aug 22 peter 30 #include <yat/utility/split.h>
4198 19 Aug 22 peter 31 #include <yat/utility/utility.h>
4198 19 Aug 22 peter 32 #include <yat/utility/Vector.h>
4198 19 Aug 22 peter 33
4198 19 Aug 22 peter 34 using namespace theplu;
4198 19 Aug 22 peter 35 using namespace yat;
4198 19 Aug 22 peter 36 using regression::Cox;
4198 19 Aug 22 peter 37 using regression::MultiCox;
4198 19 Aug 22 peter 38
4198 19 Aug 22 peter 39 void test1(yat::test::Suite& suite);
4198 19 Aug 22 peter 40
4198 19 Aug 22 peter 41 int main(int argc, char* argv[])
4198 19 Aug 22 peter 42 {
4198 19 Aug 22 peter 43   yat::test::Suite suite(argc, argv);
4198 19 Aug 22 peter 44   try {
4198 19 Aug 22 peter 45     test1(suite);
4198 19 Aug 22 peter 46   }
4198 19 Aug 22 peter 47   catch (std::runtime_error& e) {
4198 19 Aug 22 peter 48     std::cerr << "error: " << e.what() << "\n";
4198 19 Aug 22 peter 49     return EXIT_FAILURE;
4198 19 Aug 22 peter 50   }
4252 18 Nov 22 peter 51   return suite.return_value();
4198 19 Aug 22 peter 52 }
4198 19 Aug 22 peter 53
4198 19 Aug 22 peter 54
4198 19 Aug 22 peter 55 void parse(const std::vector<std::string>& data,
4198 19 Aug 22 peter 56            size_t offset, int value, yat::utility::Vector& time,
4198 19 Aug 22 peter 57            yat::utility::VectorMutable& x, std::vector<char>& y)
4198 19 Aug 22 peter 58 {
4198 19 Aug 22 peter 59   for (size_t i=0; i<data.size(); ++i) {
4198 19 Aug 22 peter 60     size_t j = offset + i;
4198 19 Aug 22 peter 61     x(j) = value;
4198 19 Aug 22 peter 62     if (*data[i].rbegin() == '*') {
4198 19 Aug 22 peter 63       time(j)=yat::utility::convert<double>(data[i].substr(0,data[i].size()-1));
4198 19 Aug 22 peter 64       y[j] = 0;
4198 19 Aug 22 peter 65     }
4198 19 Aug 22 peter 66     else {
4198 19 Aug 22 peter 67       time(j)=yat::utility::convert<double>(data[i]);
4198 19 Aug 22 peter 68       y[j] = 1;
4198 19 Aug 22 peter 69     }
4198 19 Aug 22 peter 70   }
4198 19 Aug 22 peter 71 }
4198 19 Aug 22 peter 72
4198 19 Aug 22 peter 73
4198 19 Aug 22 peter 74 bool equal(test::Suite& suite, const std::string& message,
4198 19 Aug 22 peter 75            double a, double b, double margin)
4198 19 Aug 22 peter 76 {
4198 19 Aug 22 peter 77   suite.out() << message << "\n";
4198 19 Aug 22 peter 78   return suite.add(suite.equal_fix(a, b, margin));
4198 19 Aug 22 peter 79 }
4198 19 Aug 22 peter 80
4198 19 Aug 22 peter 81
4198 19 Aug 22 peter 82 void test1(yat::test::Suite& suite)
4198 19 Aug 22 peter 83 {
4198 19 Aug 22 peter 84   // Test example taken from
4198 19 Aug 22 peter 85   // https://www.statsdirect.com/help/survival_analysis/cox_regression.htm
4198 19 Aug 22 peter 86
4198 19 Aug 22 peter 87   std::vector<std::string> stage3;
4198 19 Aug 22 peter 88   yat::utility::split(stage3, "6,19,32,42,42,43*,94,126*,169*,207,211*,227*,253,255*,270*,310*,316*,335*,346*", ',');
4198 19 Aug 22 peter 89   size_t n1 = stage3.size();
4198 19 Aug 22 peter 90
4198 19 Aug 22 peter 91   std::vector<std::string> stage4;
4198 19 Aug 22 peter 92   yat::utility::split(stage4, "4,6,10,11,11,11,13,17,20,20,21,22,24,24,29,30,30,31,33,34,35,39,40,41*,43*,45,46,50,56,61*,61*,63,68,82,85,88,89,90,93,104,110,134,137,160*,169,171,173,175,184,201,222,235*,247*,260*,284*,290*,291*,302*,304*,341*,345*", ',');
4198 19 Aug 22 peter 93   size_t n2 = stage4.size();
4198 19 Aug 22 peter 94   size_t n = n1 + n2;
4198 19 Aug 22 peter 95
4198 19 Aug 22 peter 96   yat::utility::Vector time(n);
4198 19 Aug 22 peter 97   yat::utility::Vector x(n);
4198 19 Aug 22 peter 98   std::vector<char> event(n);
4198 19 Aug 22 peter 99
4198 19 Aug 22 peter 100   parse(stage3, 0, 5, time, x, event);
4198 19 Aug 22 peter 101   parse(stage4, n1, 6, time, x, event);
4198 19 Aug 22 peter 102
4198 19 Aug 22 peter 103   Cox cox;
4198 19 Aug 22 peter 104   cox.add(x, time, event);
4198 19 Aug 22 peter 105   cox.train();
4198 19 Aug 22 peter 106
4198 19 Aug 22 peter 107   /*
4198 19 Aug 22 peter 108     Test case is taken from
4198 19 Aug 22 peter 109     https://www.statsdirect.com/help/survival_analysis/cox_regression.htm
4198 19 Aug 22 peter 110     and modified to refelect that we use Efron's method for ties
4198 19 Aug 22 peter 111     rather than Breslow's as below.
4198 19 Aug 22 peter 112
4198 19 Aug 22 peter 113     80 subjects with 54 events
4198 19 Aug 22 peter 114     Deviance (likelihood ratio) chi-square = 7.634383 df = 1 P = 0.0057
4198 19 Aug 22 peter 115     Stage group b1 = 0.96102 z = 2.492043 P = 0.0127
4198 19 Aug 22 peter 116     Cox regression - hazard ratios
4198 19 Aug 22 peter 117     ===
4198 19 Aug 22 peter 118     Parameter   Hazard ratio   95% CI
4198 19 Aug 22 peter 119     Stage group   2.614362   1.227756 to 5.566976
4198 19 Aug 22 peter 120     ===
4198 19 Aug 22 peter 121     Parameter   Coefficient   Standard Error
4198 19 Aug 22 peter 122     Stage group   0.96102   0.385636
4198 19 Aug 22 peter 123     ===
4198 19 Aug 22 peter 124     Cox regression - model analysis
4198 19 Aug 22 peter 125     Log likelihood with no covariates = -207.554801
4198 19 Aug 22 peter 126     Log likelihood with all model covariates = -203.737609
4198 19 Aug 22 peter 127     Deviance (likelihood ratio) chi-square = 7.634383 df = 1 P = 0.0057
4198 19 Aug 22 peter 128   */
4198 19 Aug 22 peter 129
4198 19 Aug 22 peter 130   suite.out() << "b: " << cox.b() << "\n";
4198 19 Aug 22 peter 131   suite.out() << "z: " << cox.z() << "\n";
4198 19 Aug 22 peter 132   suite.out() << "P: " << cox.p() << "\n";
4198 19 Aug 22 peter 133
4198 19 Aug 22 peter 134   suite.out() << "HR: " << cox.hazard_ratio() << "\n";
4198 19 Aug 22 peter 135   suite.out() << "HR CI: " << cox.hazard_ratio_lower_CI() << " to "
4198 19 Aug 22 peter 136             << cox.hazard_ratio_upper_CI() << "\n";
4198 19 Aug 22 peter 137
4198 19 Aug 22 peter 138   equal(suite, "testing b", cox.b(), 0.96132, 1e-5);
4198 19 Aug 22 peter 139   equal(suite, "testing z", cox.z(), 2.492922, 1e-6);
4198 19 Aug 22 peter 140   equal(suite, "testing p", cox.p(), 0.0127, 1e-4);
4198 19 Aug 22 peter 141   equal(suite, "testing HR", cox.hazard_ratio(), 2.61515, 1e-6);
4198 19 Aug 22 peter 142   equal(suite, "testing HR lower CI", cox.hazard_ratio_lower_CI(),
4198 19 Aug 22 peter 143         1.228163, 1e-6);
4198 19 Aug 22 peter 144   equal(suite, "testing HR upper_CI", cox.hazard_ratio_upper_CI(),
4198 19 Aug 22 peter 145         5.568489,5e-6);
4198 19 Aug 22 peter 146
4198 19 Aug 22 peter 147   suite.out() << "\n===\nTesting MultiCox\n";
4198 19 Aug 22 peter 148   MultiCox multicox;
4198 19 Aug 22 peter 149   for (size_t i=0; i<x.size(); ++i)
4198 19 Aug 22 peter 150     multicox.add(yat::utility::Vector(1, x(i)), time(i), event[i]);
4198 19 Aug 22 peter 151   multicox.train();
4198 19 Aug 22 peter 152
4198 19 Aug 22 peter 153   equal(suite, "testing b", multicox.b(0), 0.96132, 1e-5);
4198 19 Aug 22 peter 154   equal(suite, "testing z", multicox.z(0), 2.492922, 1e-6);
4198 19 Aug 22 peter 155   equal(suite, "testing p", multicox.p(0), 0.0127, 1e-4);
4198 19 Aug 22 peter 156   equal(suite, "testing HR", multicox.hazard_ratio(0), 2.61515, 1e-6);
4198 19 Aug 22 peter 157   equal(suite, "testing HR lower CI", multicox.hazard_ratio_lower_CI(0),
4198 19 Aug 22 peter 158         1.228163, 5e-6);
4198 19 Aug 22 peter 159   equal(suite, "testing HR upper_CI", multicox.hazard_ratio_upper_CI(0),
4198 19 Aug 22 peter 160         5.568489,5e-6);
4198 19 Aug 22 peter 161 }

Generated on Sat Oct 14 02:20:49 2023 (UTC) by svndigest 0.11.2pre (dev build 1672)