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

// $Id: Cox.cc 4252 2022-11-18 02:54:04Z 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/Steffenson.h>
#include <yat/utility/VectorBase.h>

#include <gsl/gsl_cdf.h>

#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 Score
    {
    public:
      Score(const std::vector<TimePoint>& times);
      double operator()(double beta) const;
      double derivative(double beta) const;
    private:
      const std::vector<TimePoint>& times_;
    };
  };


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

    prepare_times();
    // score is derivative of logL
    Score score(times_);
    for (double b=0; b<1.1; b+=0.1)
      score(b);
    utility::Steffenson solver;
    beta_ = solver(score, 0.0,
                   utility::RootFinderDerivative::Delta(0.0, 1e-5));
    if (std::isnan(beta_))
      throw std::runtime_error("beta is NaN");
    // 2nd derivative of logL is 1st derivative of score
    double hessian = score.derivative(beta_);
    beta_std_error_ = 1.0 / std::sqrt(-hessian);
  }


  /*
    Without ties the log-likelihood is
    logL = sum_i (x_i * beta - log sum_j x_j * beta)
    where i runs over all events and j runs over all data points j
    such that t_j >= t_i

    Setting theta = x * beta we have
    logL = sum_i (theta_i - log sum_j theta_j) =
         = sum_i (theta_i - log theta_Q_i) =

    theta = beta * x -> dtheta/dbeta = x

    We handle ties using Efron's method. Let denote m_i number of data
    points at t_i, H_i the indices of events at time t_i.

    logL = sum_i (theta_H_i - sum_k^m_i-1 log(theta_Q_i - k/m_i theta_H_i))

    where theta_H_i is a sum of theta running over H_i.

    The derivative (wrt beta)
    l' = sum_i(x_H_i - sum_k^m_i-1 (theta*x_Q_i - k/m_j theta+x_H_i) / (theta_Q_i - k/m_j theta_H_i))

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


  double Cox::Impl::Score::operator()(double beta) const
  {
    double score = 0;
    // variables with suffix _Q denote sums running over data points
    // (including events and censored data points) at current time and
    // future
    double theta_Q = 0;
    double thetaX_Q = 0;

    for (auto time = times_.rbegin(); time!=times_.rend(); ++time) {
      // variables with suffix _H denote sums running over events at
      // the current time.
      double theta_H = 0;
      double thetaX_H = 0;
      for (auto it = time->events_begin(); it!=time->events_end(); ++it) {
        double theta = it->theta(beta);
        theta_H += theta;
        thetaX_H += theta * it->x;
      }
      theta_Q += theta_H;
      thetaX_Q += thetaX_H;

      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) {
        score += it->x;

        const size_t k = it - time->events_begin();
        double r = static_cast<double>(k) / time->size();

        assert(theta_Q > r * theta_H);

        score -= (thetaX_Q - r * thetaX_H) / (theta_Q - r * theta_H);
      }
    }

    assert(!std::isnan(score));
    return score;
  }


  double Cox::Impl::Score::derivative(double beta) const
  {
    double deriv = 0;
    // variables with suffix _Q denote sums running over data points
    // (including events and censored data points) at current time and
    // future
    double theta_Q = 0;
    double thetaX_Q = 0;
    double thetaXX_Q = 0;
    double XX_Q = 0;

    for (auto time = times_.rbegin(); time!=times_.rend(); ++time) {
      // variables with suffix _H denote sums running over events at
      // the current time.
      double theta_H = 0;
      double thetaX_H = 0;
      double thetaXX_H = 0;
      double XX_H = 0;
      for (auto it = time->events_begin(); it!=time->events_end(); ++it) {
        double theta = it->theta(beta);
        theta_H += theta;
        thetaX_H += theta * it->x;
        thetaXX_H += theta * it->x * it->x;
        XX_H += it->x * it->x;
      }
      theta_Q += theta_H;
      thetaX_Q += thetaX_H;
      thetaXX_Q += thetaXX_H;
      XX_Q += XX_H;

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

      // f = g/h
      // g = - (thetaX_Q - r*thetaX_H)
      // g'= - (thetaXX_Q - r*thetaXX_H)
      //
      // h = theta_Q - r*theta_H
      // h'= thetaX_Q - r*thetaX_H =

      // f' = (g/h)' = (g'h - gh')/h^2 =

      // 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();
        double g = - (thetaX_Q - r*thetaX_H);
        double dg = - (thetaXX_Q - r * thetaXX_H);

        double h = (theta_Q - r*theta_H);
        double dh= (thetaX_Q - r*thetaX_H);
        deriv += (dg*h - g*dh) / (h*h);
      }
    }

    assert(!std::isnan(deriv));
    return deriv;
  }

  // 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();
  }

}}}