source: trunk/yat/regression/Cox.cc @ 4198

// $Id: Cox.cc 4198 2022-08-19 06:26:14Z peter $

/*
  Copyright (C) 2022 Peter Johansson

  This file is part of the yat library, https://dev.thep.lu.se/yat

  The yat library is free software; you can redistribute it and/or
  modify it under the terms of the GNU General Public License as
  published by the Free Software Foundation; either version 3 of the
  License, or (at your option) any later version.

  The yat library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with yat. If not, see <https://www.gnu.org/licenses/>.
*/

#include <config.h>

#include "Cox.h"

#include "detail/Cox.h"

#include <yat/utility/VectorBase.h>

#include <gsl/gsl_cdf.h>

#include <boost/math/tools/roots.hpp>

#include <algorithm>
#include <memory>

namespace theplu {
namespace yat {
namespace regression {

  class Cox::Impl : public cox::Implementation<double>
  {
  public:
    using cox::Implementation<double>::add;
    void add(const yat::utility::VectorBase& x,
             const yat::utility::VectorBase& time,
             const std::vector<char>& event)
    {
      assert(x.size() == time.size());
      assert(x.size() == event.size());
      for (size_t i=0; i<x.size(); ++i)
        add(x(i), time(i), event[i]);
    }

    double b(void) const { return beta_ ; }
    double hazard_ratio(void) const
    { return exp(beta_); }

    double hazard_ratio_lower_CI(double alpha) const
    { return exp(beta_ - hazard_ratio_CI(alpha)); }

    double hazard_ratio_upper_CI(double alpha) const
    { return exp(beta_ + hazard_ratio_CI(alpha)); }

    double p(void) const
    { return 2 * gsl_cdf_ugaussian_Q(std::abs(z())); }

    void train(void);

    double z(void) const { return beta_ / beta_std_error_; }

  private:
    double hazard_ratio_CI(double alpha) const
    {
      double z = gsl_cdf_ugaussian_Qinv(0.5 * (1.0 - alpha));
      return z * beta_std_error_;
    }
    double beta_;
    double beta_std_error_;

    class logL
    {
    public:
      logL(const std::vector<TimePoint>& times);
      std::pair<double, double> operator()(double beta) const;

      double hessian(double beta) const;
    private:
      const std::vector<TimePoint>& times_;
    };
  };


  void Cox::Impl::train(void)
  {
    if (data_.empty())
      return;
    beta_ = 0;

    prepare_times();
    logL func(times_);
    using boost::math::tools::newton_raphson_iterate;
    beta_ = newton_raphson_iterate(func, beta_, -1e42, 1e42, 30);
    if (std::isnan(beta_))
      throw std::runtime_error("beta is NaN");
    // Calculate 2nd deriviate at beta_;
    double hessian = func.hessian(beta_);
    beta_std_error_ = 1.0 / std::sqrt(hessian);
  }


  Cox::Impl::logL::logL(const std::vector<TimePoint>& times)
    : times_(times)
  {
  }


  std::pair<double, double> Cox::Impl::logL::operator()(double beta) const
  {
    // Using Efron's method:
    //
    // sort data wrt time and denote unique time t_i such that t_i <
    // t_j iff i < j
    // Let denote H_j the indices of events at time t_j, i.e.,
    // Y_i = t_j and event_i = true; n_j = |H_j|
    //
    // log partial likelihood
    // logL =
    //  sum_j (sum_i x_i*beta - sum_k log{sum_i theta_i - k/n_j sum_i theta_i})
    // where j runs over all unique times
    // i runs over all events in H_j
    // k runs over all events in H_j
    // 1st i sum runs : Y_i >= t_j
    // 2nd i sum runs over H_j
    //
    // and the derivative is
    // deriv =
    //  sum_j (sum_i x_i - sum_k {[sum_i theta_i*x_i - k/n_j sum_i theta_i*x_i] / [sum_i theta_i - k/n_j sum_i theta_i]})
    // 1st i sum runs : Y_i >= t_j
    // 2nd i sum runs over H_j
    // 3rd i sum runs : Y_i >= t_j
    // 4th i sum runs over H_j


    /*
      We handle tied ties using Efron's method. Let denote H_j the
      indices of events at time t_j

      logL =
      \sum_j (sum_i(theta_i) - \sum_k(log(sum_i(theta_i) - k/m_j sum_i(theta_i))

      where theta_i = beta * x_i
      where j runs over all unique time points
      1st i runs over H_j, i.e., all events at time t_j.
      k runs over H_j
      2nd i runs over all i: Y_i > t_j
      3rd i runs over H_j
    */

    double logL = 0;
    double deriv = 0;

    double theta_Q = 0;
    double thetaX_Q = 0;
    for (auto time = times_.rbegin(); time!=times_.rend(); ++time) {
      double sum_event_theta = 0;
      double sum_event_thetaX = 0;
      for (auto it = time->events_begin(); it!=time->events_end(); ++it) {
        double theta = it->theta(beta);
        sum_event_theta += theta;
        sum_event_thetaX += theta * it->x;
      }
      theta_Q += sum_event_theta;
      thetaX_Q += sum_event_thetaX;

      for (auto it = time->censored_begin(); it!=time->censored_end(); ++it) {
        double theta = it->theta(beta);
        theta_Q += theta;
        thetaX_Q += theta * it->x;
      }


      // loop over events at time point t
      for (auto it = time->events_begin(); it!=time->events_end(); ++it) {
        const size_t k = it - time->events_begin();
        double r = static_cast<double>(k) / time->size();

        logL += it->x * beta;
        assert(theta_Q > 0.0);
        assert(theta_Q > r * sum_event_theta);
        logL -= std::log(theta_Q - r * sum_event_theta);

        deriv += it->x;
        deriv -= (thetaX_Q - r * sum_event_thetaX) /
          (theta_Q - r * sum_event_theta);
      }

    }

    assert(!std::isnan(logL));
    assert(!std::isnan(deriv));
    return std::make_pair(-logL, -deriv);
  }


  double Cox::Impl::logL::hessian(double beta) const
  {
    // The second derivative of the log-likelihood evaluated at the
    // maximum likelihood estimates (MLE) is the observed Fisher
    // information

    double hessian = 0;

    double sum_theta = 0;
    double sum_thetaX = 0;
    double sum_thetaXX = 0;
    for (const auto& t : times_) {
      for (auto it = t.begin; it!=t.end; ++it) {
        double theta = it->theta(beta);
        sum_theta += theta;
        sum_thetaX += theta * it->x;
        sum_thetaXX += theta * it->x * it->x;
      }
    }

    // loop over unique time points
    for (const auto& t : times_) {

      // sum over all events in H_j
      double part_sum_theta = 0;
      double part_sum_thetaX = 0;
      double part_sum_thetaXX = 0;
      for (auto it = t.events_begin(); it!=t.events_end(); ++it) {
        double theta = it->theta(beta);
        part_sum_theta += theta;
        part_sum_thetaX += theta * it->x;
        part_sum_thetaXX += theta * it->x * it->x;
      }

      // loop over events at time point t
      for (auto it = t.events_begin(); it!=t.events_end(); ++it) {
        const size_t k = it - t.events_begin();
        double r = static_cast<double>(k) / t.size();

        double S_thetaXX = sum_thetaXX - r * part_sum_thetaXX;
        double S_thetaX  = sum_thetaX  - r * part_sum_thetaX;
        double S_theta   = sum_theta   - r * part_sum_theta;
        hessian += S_thetaXX/S_theta;
        hessian -= std::pow(S_thetaX/S_theta, 2);
      }

      // update the cumulative sums
      for (auto it = t.begin; it!=t.end; ++it) {
        double theta = it->theta(beta);
        sum_theta -= theta;
        sum_thetaX -= theta * it->x;
        sum_thetaXX -= theta * it->x * it->x;
      }
    }
    return hessian;
  }


  // class Cox

  Cox::Cox(void)
    : pimpl_(new Impl)
  {
  }


  Cox::Cox(const Cox& other)
    : pimpl_(new Impl(*other.pimpl_))
  {
  }


  Cox::Cox(Cox&& other)
  {
    std::swap(pimpl_, other.pimpl_);
  }


  Cox::~Cox(void)
  {
  }


  Cox& Cox::operator=(const Cox& other)
  {
    assert(other.pimpl_);
    pimpl_.reset(new Impl(*other.pimpl_));
    return *this;
  }


  Cox& Cox::operator=(Cox&& other)
  {
    std::swap(pimpl_, other.pimpl_);
    return *this;
  }


  void Cox::Cox::add(double x, double time, bool event)
  {
    pimpl_->add(x, time, event);
  }


  void Cox::add(const yat::utility::VectorBase& x,
                const yat::utility::VectorBase& time,
                const std::vector<char>& event)
  {
    pimpl_->add(x, time, event);
  }


  double Cox::b(void) const
  {
    return pimpl_->b();
  }


  void Cox::clear(void)
  {
    pimpl_->clear();
  }


  double Cox::hazard_ratio(void) const
  {
    return pimpl_->hazard_ratio();
  }


  double Cox::hazard_ratio_lower_CI(double alpha) const
  {
    return pimpl_->hazard_ratio_lower_CI(alpha);
  }


  double Cox::hazard_ratio_upper_CI(double alpha) const
  {
    return pimpl_->hazard_ratio_upper_CI(alpha);
  }


  double Cox::p(void) const
  {
    return pimpl_->p();
  }


  void Cox::train(void)
  {
    pimpl_->train();
  }


  double Cox::z(void) const
  {
    return pimpl_->z();
  }

}}}