naginterfaces.library.surviv.coxmodel

naginterfaces.library.surviv.coxmodel(ns, z, isz, t, ic, isi, b, ndmax, omega=None, tol=1.1102230246251565e-14, maxit=1000, iprint=0, io_manager=None)

coxmodel returns parameter estimates and other statistics that are associated with the Cox proportional hazards model for fixed covariates.

For full information please refer to the NAG Library document for g12ba

https://support.nag.com/numeric/nl/nagdoc_30/flhtml/g12/g12baf.html

Parameters

nsint

The number of strata. If $ns>0$ then the stratum for each observation must be supplied in $isi$.

zfloat, array-like, shape $(n,m)$

The $i$th row must contain the covariates which are associated with the $i$th failure time given in $t$.

iszint, array-like, shape $\left(m\right)$

Indicates which subset of covariates is to be included in the model.

$isz[j-1]\ge 1$

The $j$th covariate is included in the model.

$isz[j-1]=0$

The $j$th covariate is excluded from the model and not referenced.

tfloat, array-like, shape $\left(n\right)$

The vector of $n$ failure censoring times.

icint, array-like, shape $\left(n\right)$

The status of the individual at time $t$ given in $t$.

$ic[i-1]=0$

The $i$th individual has failed at time $t[i-1]$.

$ic[i-1]=1$

The $i$th individual has been censored at time $t[i-1]$.

isiint, array-like, shape $(:)$

Note: the required length for this argument is determined as follows: if $ns>0$: $n$; otherwise: $1$.

If $ns>0$, the stratum indicators which also allow data points to be excluded from the analysis.

If $ns=0$, $isi$ is not referenced.

$isi[i-1]=k$

The $i$th data point is in the $k$th stratum, where $k=1,2,\dots ,ns$.

$isi[i-1]=0$

The $i$th data point is omitted from the analysis.

bfloat, array-like, shape $\left(\text{ip}\right)$

Initial estimates of the covariate coefficient parameters $\beta$. $b[j-1]$ must contain the initial estimate of the coefficient of the covariate in $z$ corresponding to the $j$th nonzero value of $isz$.

Suggested value: in many cases an initial value of zero for $b[j-1]$ may be used. For other suggestions see Further Comments.

ndmaxint

The dimension of the array $tp$.

The first dimension of the array $sur$.

omegaNone or float, array-like, shape $\left(n\right)$, optional

If $omega\text{is not}None$, the offset, ${\omega }_{\text{i}}$, for $\text{i}=1,2,\dots ,n$. Otherwise $omega$ is not referenced.

tolfloat, optional

Indicates the accuracy required for the estimation. Convergence is assumed when the decrease in deviance is less than $tol\times (1.0+CurrentDeviance)$. This corresponds approximately to an absolute precision if the deviance is small and a relative precision if the deviance is large.

maxitint, optional

The maximum number of iterations to be used for computing the estimates. If $maxit$ is set to $0$ then the standard errors, score functions, variance-covariance matrix and the survival function are computed for the input value of $\beta$ in $b$ but $\beta$ is not updated.

iprintint, optional

Indicates if the printing of information on the iterations is required.

$iprint\le 0$

No printing.

$iprint\ge 1$

The deviance and the current estimates are printed every $iprint$ iterations. When printing occurs the output is directed to the file object associated with the advisory I/O unit (see FileObjManager).

io_managerFileObjManager, optional

Manager for I/O in this routine.

Returns

devfloat

The deviance, that is $-2\times$(maximized log marginal likelihood).

bfloat, ndarray, shape $\left(\text{ip}\right)$

$b[j-1]$ contains the estimate ${\hat{\beta }}_i$, the coefficient of the covariate stored in the $i$th column of $z$ where $i$ is the $j$th nonzero value in the array $isz$.

sefloat, ndarray, shape $\left(\text{ip}\right)$

$se[\text{j}-1]$ is the asymptotic standard error of the estimate contained in $b[\text{j}-1]$ and score function in $sc[\text{j}-1]$, for $\text{j}=1,2,\dots ,\text{ip}$.

scfloat, ndarray, shape $\left(\text{ip}\right)$

$sc[j-1]$ is the value of the score function, $U_j\left(\beta \right)$, for the estimate contained in $b[j-1]$.

covfloat, ndarray, shape $(\text{ip}\times (\text{ip}+1)/2)$

The variance-covariance matrix of the parameter estimates in $b$ stored in packed form by column, i.e., the covariance between the parameter estimates given in $b[i-1]$ and $b[j-1]$, $j\ge i$, is stored in $cov[j(j-1)/2+i-1]$.

resfloat, ndarray, shape $\left(n\right)$

The residuals, $r\left(t_{\text{l}}\right)$, for $\text{l}=1,2,\dots ,n$.

ndint

The number of distinct failure times.

tpfloat, ndarray, shape $\left(ndmax\right)$

$tp[\text{i}-1]$ contains the $\text{i}$th distinct failure time, for $\text{i}=1,2,\dots ,nd$.

surfloat, ndarray, shape $(ndmax,max(ns,1))$

If $ns=0$, $sur[i-1,0]$ contains the estimated survival function for the $i$th distinct failure time.

If $ns>0$, $sur[i-1,k-1]$ contains the estimated survival function for the $i$th distinct failure time in the $k$th stratum.

Raises

NagValueError

(errno $1$)

On entry, $\text{ip}=⟨value⟩$.

Constraint: $\text{ip}\ge 1$.

(errno $1$)

On entry, $tol=⟨value⟩$.

Constraint: $tol\ge 10\times \text{machine precision}$.

(errno $1$)

On entry, $maxit=⟨value⟩$.

Constraint: $maxit\ge 0$.

(errno $1$)

On entry, $ns=⟨value⟩$.

Constraint: $ns\ge 0$.

(errno $1$)

On entry, $n=⟨value⟩$.

Constraint: $n\ge 2$.

(errno $1$)

On entry, $m=⟨value⟩$.

Constraint: $m\ge 1$.

(errno $2$)

On entry, there are not $\text{ip}$ values of $isz>0$.

(errno $2$)

On entry, $i=⟨value⟩$ and $isz[i-1]=⟨value⟩$.

Constraint: $isz[i-1]\ge 0$.

(errno $2$)

On entry too few observations included in model.

(errno $2$)

On entry, $i=⟨value⟩$, $isi[i-1]=⟨value⟩$ and $ns=⟨value⟩$.

Constraint: $0\le isi[i-1]\le ns$.

(errno $2$)

On entry, $i=⟨value⟩$ and $ic[i-1]=⟨value⟩$.

Constraint: $ic[i-1]=0$ or $1$.

(errno $2$)

On entry, $ndmax=⟨value⟩$ and minimum value for $ndmax=⟨value⟩$.

Constraint: $ndmax\ge \text{number}$ of distinct failure times.

(errno $2$)

All observations are censored.

(errno $3$)

The matrix of second partial derivative is singular.

(errno $4$)

Overflow has been detected in the calculations.

Warns

NagAlgorithmicWarning

(errno $5$)

Convergence not achieved in $⟨value⟩$ iterations.

(errno $6$)

Too many step halvings required.

Notes

In the NAG Library the traditional C interface for this routine uses a different algorithmic base. Please contact NAG if you have any questions about compatibility.

The proportional hazard model relates the time to an event, usually death or failure, to a number of explanatory variables known as covariates. Some of the observations may be right-censored, that is the exact time to failure is not known, only that it is greater than a known time.

Let $t_{\text{i}}$, for $\text{i}=1,2,\dots ,n$, be the failure time or censored time for the $i$th observation with the vector of $p$ covariates $z_i$. It is assumed that censoring and failure mechanisms are independent. The hazard function, $\lambda (t,z)$, is the probability that an individual with covariates $z$ fails at time $t$ given that the individual survived up to time $t$. In the Cox proportional hazards model (see Cox (1972)) $\lambda (t,z)$ is of the form:

$$\lambda (t,z)={\lambda }_0\left(t\right)exp(z^T\beta +\omega )$$

where ${\lambda }_0$ is the base-line hazard function, an unspecified function of time, $\beta$ is a vector of unknown parameters and $\omega$ is a known offset.

Assuming there are ties in the failure times giving $n_d<n$ distinct failure times, $t_{\left(1\right)}<\cdots <t_{\left(n_d\right)}$ such that $d_i$ individuals fail at $t_{\left(i\right)}$, it follows that the marginal likelihood for $\beta$ is well approximated (see Kalbfleisch and Prentice (1980)) by:

$$L=\prod _{i=1}^{n_d}\frac{exp(s_i^T\beta +{\omega }_i)}{{\left[{\sum }_{l\in R\left(t_{\left(i\right)}\right)}exp(z_l^T\beta +{\omega }_l)\right]}^{d_i}}$$

where $s_i$ is the sum of the covariates of individuals observed to fail at $t_{\left(i\right)}$ and $R\left(t_{\left(i\right)}\right)$ is the set of individuals at risk just prior to $t_{\left(i\right)}$, that is, it is all individuals that fail or are censored at time $t_{\left(i\right)}$ along with all individuals that survive beyond time $t_{\left(i\right)}$. The maximum likelihood estimates (MLEs) of $\beta$, given by $\hat{\beta }$, are obtained by maximizing (1) using a Newton–Raphson iteration technique that includes step halving and utilizes the first and second partial derivatives of (1) which are given by equations (2) and (3) below:

$$U_j\left(\beta \right)=\frac{\partial ln\left(L\right)}{\partial {\beta }_j}=\sum _{i=1}^{n_d}[s_{ji}-d_i{\alpha }_{ji}\left(\beta \right)]=0$$

for $j=1,2,\dots ,p$, where $s_{ji}$ is the $j$th element in the vector $s_i$ and

$${\alpha }_{ji}\left(\beta \right)=\frac{{\sum }_{l\in R\left(t_{\left(i\right)}\right)}{z}_{jl}exp(z_l^T\beta +{\omega }_l)}{{\sum }_{l\in R\left(t_{\left(i\right)}\right)}exp(z_l^T\beta +{\omega }_l)}\text{.}$$

Similarly,

$$I_{hj}\left(\beta \right)=-\frac{{\partial }^{2}ln\left(L\right)}{\partial {\beta }_h\partial {\beta }_j}=\sum _{i=1}^{n_d}{d}_i{\gamma }_{hji}$$

where

$${\gamma }_{hji}=\frac{{\sum }_{l\in R\left(t_{\left(i\right)}\right)}{z}_{hl}{z}_{jl}exp(z_l^T\beta +{\omega }_l)}{{\sum }_{l\in R\left(t_{\left(i\right)}\right)}exp(z_l^T\beta +{\omega }_l)}-{\alpha }_{hi}\left(\beta \right){\alpha }_{ji}\left(\beta \right)\text{,}h,j=1,\dots ,p\text{.}$$

$U_j\left(\beta \right)$ is the $j$th component of a score vector and $I_{hj}\left(\beta \right)$ is the $(h,j)$ element of the observed information matrix $I\left(\beta \right)$ whose inverse $I{\left(\beta \right)}^{-1}={\left[I_{hj}\left(\beta \right)\right]}^{-1}$ gives the variance-covariance matrix of $\beta$.

It should be noted that if a covariate or a linear combination of covariates is monotonically increasing or decreasing with time then one or more of the ${\beta }_j$’s will be infinite.

If ${\lambda }_0\left(t\right)$ varies across $\nu$ strata, where the number of individuals in the $\text{k}$th stratum is $n_{\text{k}}$, for $\text{k}=1,2,\dots ,\nu$ with $n={\sum }_{k=1}^{\nu }{n}_k$, then rather than maximizing (1) to obtain $\hat{\beta }$, the following marginal likelihood is maximized:

$$L=\prod _{k=1}^{\nu }{L}_k\text{,}$$

where $L_k$ is the contribution to likelihood for the $n_k$ observations in the $k$th stratum treated as a single sample in (1). When strata are included the covariate coefficients are constant across strata but there is a different base-line hazard function ${\lambda }_0$.

The base-line survivor function associated with a failure time $t_{\left(i\right)}$, is estimated as $exp(-\hat{H}\left(t_{\left(i\right)}\right))$, where

$$\hat{H}\left(t_{\left(i\right)}\right)=\sum _{t_{\left(j\right)}\le t_{\left(i\right)}}\left(\frac{d_i}{{\sum }_{l\in R\left(t_{\left(j\right)}\right)}exp(z_l^T\hat{\beta }+{\omega }_l)}\right)\text{,}$$

where $d_i$ is the number of failures at time $t_{\left(i\right)}$. The residual for the $l$th observation is computed as:

$$r\left(t_l\right)=\hat{H}\left(t_l\right)exp(z_l^T\hat{\beta }+{\omega }_l)$$

where $\hat{H}\left(t_l\right)=\hat{H}\left(t_{\left(i\right)}\right),t_{\left(i\right)}\le t_l<t_{(i+1)}$. The deviance is defined as $-2\times$(logarithm of marginal likelihood). There are two ways to test whether individual covariates are significant: the differences between the deviances of nested models can be compared with the appropriate ${\chi }^2$-distribution; or, the asymptotic normality of the parameter estimates can be used to form $z$ tests by dividing the estimates by their standard errors or the score function for the model under the null hypothesis can be used to form $z$ tests.

References

Cox, D R, 1972, Regression models in life tables (with discussion), J. Roy. Statist. Soc. Ser. B (34), 187–220

Gross, A J and Clark, V A, 1975, Survival Distributions: Reliability Applications in the Biomedical Sciences, Wiley

Kalbfleisch, J D and Prentice, R L, 1980, The Statistical Analysis of Failure Time Data, Wiley

See also

naginterfaces.library.examples.surviv.coxmodel_ex.main()