NAG FL Interface
g12baf (coxmodel)

1 Purpose

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

2 Specification

Fortran Interface
Subroutine g12baf ( offset, n, m, ns, z, ldz, isz, ip, t, ic, omega, isi, dev, b, se, sc, cov, res, nd, tp, sur, ndmax, tol, maxit, iprint, wk, iwk, ifail)
Integer, Intent (In) :: n, m, ns, ldz, isz(m), ip, ic(n), isi(*), ndmax, maxit, iprint
Integer, Intent (Inout) :: ifail
Integer, Intent (Out) :: nd, iwk(2*n)
Real (Kind=nag_wp), Intent (In) :: z(ldz,m), t(n), omega(*), tol
Real (Kind=nag_wp), Intent (Inout) :: b(ip), sur(ndmax,*)
Real (Kind=nag_wp), Intent (Out) :: dev, se(ip), sc(ip), cov(ip*(ip+1)/2), res(n), tp(ndmax), wk(ip*(ip+9)/2+n)
Character (1), Intent (In) :: offset
C Header Interface
#include <nag.h>
void  g12baf_ (const char *offset, const Integer *n, const Integer *m, const Integer *ns, const double z[], const Integer *ldz, const Integer isz[], const Integer *ip, const double t[], const Integer ic[], const double omega[], const Integer isi[], double *dev, double b[], double se[], double sc[], double cov[], double res[], Integer *nd, double tp[], double sur[], const Integer *ndmax, const double *tol, const Integer *maxit, const Integer *iprint, double wk[], Integer iwk[], Integer *ifail, const Charlen length_offset)
The routine may be called by the names g12baf or nagf_surviv_coxmodel.

3 Description

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 ti, for i=1,2,,n, be the failure time or censored time for the ith observation with the vector of p covariates zi. It is assumed that censoring and failure mechanisms are independent. The hazard function, λ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)) λt,z is of the form:
λt,z=λ0texpzTβ+ω  
where λ0 is the base-line hazard function, an unspecified function of time, β is a vector of unknown parameters and ω is a known offset.
Assuming there are ties in the failure times giving nd<n distinct failure times, t1<<tnd such that di individuals fail at ti, it follows that the marginal likelihood for β is well approximated (see Kalbfleisch and Prentice (1980)) by:
L=i=1ndexpsiTβ+ωi lRtiexpzlTβ+ωldi (1)
where si is the sum of the covariates of individuals observed to fail at ti and Rti is the set of individuals at risk just prior to ti, that is, it is all individuals that fail or are censored at time ti along with all individuals that survive beyond time ti. The maximum likelihood estimates (MLEs) of β, given by β^, 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:
Ujβ= lnL βj =i=1ndsji-diαjiβ=0 (2)
for j=1,2,,p, where sji is the jth element in the vector si and
αjiβ=lRtizjlexpzlTβ+ωl lRtiexpzlTβ+ωl .  
Similarly,
Ihjβ=- 2lnL βhβj =i=1nddiγhji (3)
where
γhji=lRti zhlzjlexpzlTβ+ωl lRtiexpzlTβ+ωl -αhiβαjiβ,   h,j= 1,,p.  
Ujβ is the jth component of a score vector and Ihjβ is the h,j element of the observed information matrix Iβ whose inverse Iβ-1=Ihjβ -1 gives the variance-covariance matrix of β.
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 βj's will be infinite.
If λ0t varies across ν strata, where the number of individuals in the kth stratum is nk, for k=1,2,,ν with n=k=1νnk, then rather than maximizing (1) to obtain β^, the following marginal likelihood is maximized:
L=k=1νLk, (4)
where Lk is the contribution to likelihood for the nk observations in the kth 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 λ0.
The base-line survivor function associated with a failure time ti, is estimated as exp-H^ti, where
H^ti=tjti dilRtjexpzlTβ^+ωl , (5)
where di is the number of failures at time ti. The residual for the lth observation is computed as:
rtl= H^tlexpzlTβ^+ωl  
where H^tl=H^ti,titl<ti+1. The deviance is defined as -2×(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 χ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.

4 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

5 Arguments

1: offset Character(1) Input
On entry: indicates if an offset is to be used.
offset='Y'
An offset must be included in omega.
offset='N'
No offset is included in the model.
Constraint: offset='Y' or 'N'.
2: n Integer Input
On entry: n, the number of data points.
Constraint: n2.
3: m Integer Input
On entry: the number of covariates in array z.
Constraint: m1.
4: ns Integer Input
On entry: the number of strata. If ns>0 then the stratum for each observation must be supplied in isi.
Constraint: ns0.
5: zldzm Real (Kind=nag_wp) array Input
On entry: the ith row must contain the covariates which are associated with the ith failure time given in t.
6: ldz Integer Input
On entry: the first dimension of the array z as declared in the (sub)program from which g12baf is called.
Constraint: ldzn.
7: iszm Integer array Input
On entry: indicates which subset of covariates is to be included in the model.
iszj1
The jth covariate is included in the model.
iszj=0
The jth covariate is excluded from the model and not referenced.
Constraint: iszj0 and at least one and at most n0-1 elements of isz must be nonzero where n0 is the number of observations excluding any with zero value of isi.
8: ip Integer Input
On entry: the number of covariates included in the model as indicated by isz.
Constraints:
  • ip1;
  • ip=​ number of nonzero values of ​isz.
9: tn Real (Kind=nag_wp) array Input
On entry: the vector of n failure censoring times.
10: icn Integer array Input
On entry: the status of the individual at time t given in t.
ici=0
The ith individual has failed at time ti.
ici=1
The ith individual has been censored at time ti.
Constraint: ici=0 or 1, for i=1,2,,n.
11: omega* Real (Kind=nag_wp) array Input
Note: the dimension of the array omega must be at least n if offset='Y', and at least 1 otherwise.
On entry: if offset='Y', the offset, ωi, for i=1,2,,n. Otherwise omega is not referenced.
12: isi* Integer array Input
Note: the dimension of the array isi must be at least n if ns>0, and at least 1 otherwise.
On entry: if ns>0, the stratum indicators which also allow data points to be excluded from the analysis.
If ns=0, isi is not referenced.
isii=k
The ith data point is in the kth stratum, where k=1,2,,ns.
isii=0
The ith data point is omitted from the analysis.
Constraint: if ns>0, 0isiins and more than ip values of isii>0, for i=1,2,,n.
13: dev Real (Kind=nag_wp) Output
On exit: the deviance, that is -2×(maximized log marginal likelihood).
14: bip Real (Kind=nag_wp) array Input/Output
On entry: initial estimates of the covariate coefficient parameters β. bj must contain the initial estimate of the coefficient of the covariate in z corresponding to the jth nonzero value of isz.
Suggested value: in many cases an initial value of zero for bj may be used. For other suggestions see Section 9.
On exit: bj contains the estimate β^i, the coefficient of the covariate stored in the ith column of z where i is the jth nonzero value in the array isz.
15: seip Real (Kind=nag_wp) array Output
On exit: sej is the asymptotic standard error of the estimate contained in bj and score function in scj, for j=1,2,,ip.
16: scip Real (Kind=nag_wp) array Output
On exit: scj is the value of the score function, Ujβ, for the estimate contained in bj.
17: covip×ip+1/2 Real (Kind=nag_wp) array Output
On exit: 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 bi and bj, ji, is stored in covjj-1/2+i .
18: resn Real (Kind=nag_wp) array Output
On exit: the residuals, rtl, for l=1,2,,n.
19: nd Integer Output
On exit: the number of distinct failure times.
20: tpndmax Real (Kind=nag_wp) array Output
On exit: tpi contains the ith distinct failure time, for i=1,2,,nd.
21: surndmax* Real (Kind=nag_wp) array Output
Note: the second dimension of the array sur must be at least maxns,1.
On exit: if ns=0, suri1 contains the estimated survival function for the ith distinct failure time.
If ns>0, surik contains the estimated survival function for the ith distinct failure time in the kth stratum.
22: ndmax Integer Input
On entry: the dimension of the array tp and the first dimension of the array sur as declared in the (sub)program from which g12baf is called.
Constraint: ndmaxthe number of distinct failure times. This is returned in ​nd.
23: tol Real (Kind=nag_wp) Input
On entry: indicates the accuracy required for the estimation. Convergence is assumed when the decrease in deviance is less than tol×1.0+CurrentDeviance. This corresponds approximately to an absolute precision if the deviance is small and a relative precision if the deviance is large.
Constraint: tol10×machine precision.
24: maxit Integer Input
On entry: 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 β in b but β is not updated.
Constraint: maxit0.
25: iprint Integer Input
On entry: indicates if the printing of information on the iterations is required.
iprint0
No printing.
iprint1
The deviance and the current estimates are printed every iprint iterations. When printing occurs the output is directed to the current advisory message unit (see x04abf).
26: wkip×ip+9/2+n Real (Kind=nag_wp) array Workspace
27: iwk2×n Integer array Workspace
28: ifail Integer Input/Output
On entry: ifail must be set to 0, -1 or 1 to set behaviour on detection of an error; these values have no effect when no error is detected.
A value of 0 causes the printing of an error message and program execution will be halted; otherwise program execution continues. A value of -1 means that an error message is printed while a value of 1 means that it is not.
If halting is not appropriate, the value -1 or 1 is recommended. If message printing is undesirable, then the value 1 is recommended. Otherwise, the value 0 is recommended. When the value -1 or 1 is used it is essential to test the value of ifail on exit.
On exit: ifail=0 unless the routine detects an error or a warning has been flagged (see Section 6).

6 Error Indicators and Warnings

If on entry ifail=0 or -1, explanatory error messages are output on the current error message unit (as defined by x04aaf).
Errors or warnings detected by the routine:
ifail=1
On entry, ip=value.
Constraint: ip1.
On entry, ldz=value.
Constraint: ldzn.
On entry, m=value.
Constraint: m1.
On entry, maxit=value.
Constraint: maxit0.
On entry, n=value.
Constraint: n2.
On entry, ns=value.
Constraint: ns0.
On entry, offset=value.
Constraint: offset='Y' or 'N'.
On entry, tol=value.
Constraint: tol10×machine precision.
ifail=2
All observations are censored.
On entry, i=value, isii=value and ns=value.
Constraint: 0isiins.
On entry, i=value and ici=value.
Constraint: ici=0 or 1.
On entry, i=value and iszi=value.
Constraint: iszi0.
On entry, ndmax=value and minimum value for ndmax=value.
Constraint: ndmaxnumber of distinct failure times.
On entry, there are not ip values of isz>0.
On entry too few observations included in model.
ifail=3
The matrix of second partial derivative is singular. Try different starting values or include fewer covariates.
ifail=4
Overflow has been detected in the calculations. Try using different starting values.
ifail=5
Convergence not achieved in value iterations. The progress toward convergence can be examined by using a nonzero value of iprint. Any non-convergence may be due to a linear combination of covariates being monotonic with time. The full results are returned.
ifail=6
Too many step halvings required. In the current iteration 10 step halvings have been preformed without decreasing the deviance from the previous iteration. The process is assumed to have converged.
ifail=-99
An unexpected error has been triggered by this routine. Please contact NAG.
See Section 7 in the Introduction to the NAG Library FL Interface for further information.
ifail=-399
Your licence key may have expired or may not have been installed correctly.
See Section 8 in the Introduction to the NAG Library FL Interface for further information.
ifail=-999
Dynamic memory allocation failed.
See Section 9 in the Introduction to the NAG Library FL Interface for further information.

7 Accuracy

The accuracy is specified by tol.

8 Parallelism and Performance

g12baf is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.
g12baf makes calls to BLAS and/or LAPACK routines, which may be threaded within the vendor library used by this implementation. Consult the documentation for the vendor library for further information.
Please consult the X06 Chapter Introduction for information on how to control and interrogate the OpenMP environment used within this routine. Please also consult the Users' Note for your implementation for any additional implementation-specific information.

9 Further Comments

g12baf uses mean centering which involves subtracting the means from the covariables prior to computation of any statistics. This helps to minimize the effect of outlying observations and accelerates convergence.
If the initial estimates are poor then there may be a problem with overflow in calculating expβTzi or there may be non-convergence. Reasonable estimates can often be obtained by fitting an exponential model using g02gcf.

10 Example

The data are the remission times for two groups of leukemia patients (see page 242 of Gross and Clark (1975)). A dummy variable indicates which group they come from. An initial estimate is computed using the exponential model and then the Cox proportional hazard model is fitted and parameter estimates and the survival function are printed.

10.1 Program Text

Program Text (g12bafe.f90)

10.2 Program Data

Program Data (g12bafe.d)

10.3 Program Results

Program Results (g12bafe.r)