https://github.com/cran/survival
Tip revision: 38a8ccc0fb014c2b2bfbb2b5153b419e615b01ba authored by Terry M Therneau on 16 April 2016, 19:26:46 UTC
version 2.39-2
version 2.39-2
Tip revision: 38a8ccc
coxsafe.c
/*
** A very few pathologic cases can cause the Newton Raphson iteration
** path in coxph to generate a horrific argument to exp(). Since all these
** calls to exp result in (essentially) relative risks we choose a
** fixed value of LARGE on biological grounds: any number less than
** 1/(population of the earth) is essentially a zero, that is, an exponent
** outside the range of +-23.
** A sensible numeric limit would be log(.Machine$double.xmax) which is
** about 700, perhaps divided by 2 or log(n) to keep a few more bits.
** However, passing this down the R calling chain to the c-routine is a lot
** more hassle than I want to implement for this very rare case.
**
** Actually, the argument does not have to get large enough to have any
** single exponential overflow. In (start, stop] data we keep a running
** sum of scores exp(x[i]*beta), which involves both adding subjects in and
** subtracting them out. An outlier x value that enters and then leaves can
** erase all the digits of accuracy. Most machines have about 16 digits of
** accuracy and exp(21) uses up about 9 of them, leaving enough that the
** routine doesn't fall on it's face. (A user data set with outlier that
** got exp(54) and a overlarge first beta on the first iteration led to this
** paragraph.) When beta-hat is infinite and x well behaved, the loglik
** usually converges before xbeta gets to 15, so this protection should not
** harm the iteration path of even edge cases; only fix those that truely
** go astray.
**
** The truncation turns out not to be necessary for small values, since a risk
** score of exp(-50) or exp(-1000) or 0 all give essentially the same effect.
** We only cut these off enough to avoid underflow.
*/
#define LARGE 22
#define SMALL -200
double coxsafe(double x) {
if (x< SMALL) return(SMALL);
if (x> LARGE) return(LARGE);
return (x);
}