NAG Library Routine Document
G12BAF
1 Purpose
G12BAF returns parameter estimates and other statistics that are associated with the Cox proportional hazards model for fixed covariates.
2 Specification
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 
N, M, NS, LDZ, ISZ(M), IP, IC(N), ISI(*), ND, NDMAX, MAXIT, IPRINT, IWK(2*N), IFAIL 
REAL (KIND=nag_wp) 
Z(LDZ,M), T(N), OMEGA(*), DEV, B(IP), SE(IP), SC(IP), COV(IP*(IP+1)/2), RES(N), TP(NDMAX), SUR(NDMAX,*), TOL, WK(IP*(IP+9)/2+N) 
CHARACTER(1) 
OFFSET 

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 rightcensored, that is the exact time to failure is not known, only that it is greater than a known time.
Let
${t}_{\mathit{i}}$, for
$\mathit{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 \left(t,z\right)$, 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 \left(t,z\right)$ is of the form:
where
${\lambda}_{0}$ is the baseline 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:
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:
for
$j=1,2,\dots ,p$, where
${s}_{ji}$ is the
$j$th element in the vector
${s}_{i}$ and
Similarly,
where
${U}_{j}\left(\beta \right)$ is the
$j$th component of a score vector and
${I}_{hj}\left(\beta \right)$ is the
$\left(h,j\right)$ 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 variancecovariance matrix of
$\mathbf{\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
$\mathit{k}$th stratum is
${n}_{\mathit{k}}$, for
$\mathit{k}=1,2,\dots ,\nu $ with
$n={\displaystyle \sum _{k=1}^{\nu}}{n}_{k}$, then rather than maximizing
(1) to obtain
$\hat{\beta}$, the following marginal likelihood is maximized:
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 baseline hazard function
${\lambda}_{0}$.
The baseline survivor function associated with a failure time
${t}_{\left(i\right)}$, is estimated as
$\mathrm{exp}\left(\hat{H}\left({t}_{\left(i\right)}\right)\right)$, where
where
${d}_{i}$ is the number of failures at time
${t}_{\left(i\right)}$. The residual for the
$l$th observation is computed as:
where
$\hat{H}\left({t}_{l}\right)=\hat{H}\left({t}_{\left(i\right)}\right),{t}_{\left(i\right)}\le {t}_{l}<{t}_{\left(i+1\right)}$. The deviance is defined as
$2\times \text{}$(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.
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 Parameters
 1: OFFSET – CHARACTER(1)Input
On entry: indicates if an offset is to be used.
 ${\mathbf{OFFSET}}=\text{'Y'}$
 An offset must be included in OMEGA.
 ${\mathbf{OFFSET}}=\text{'N'}$
 No offset is included in the model.
Constraint:
${\mathbf{OFFSET}}=\text{'Y'}$ or $\text{'N'}$.
 2: N – INTEGERInput
On entry: $n$, the number of data points.
Constraint:
${\mathbf{N}}\ge 2$.
 3: M – INTEGERInput
On entry: the number of covariates in array
Z.
Constraint:
${\mathbf{M}}\ge 1$.
 4: NS – INTEGERInput
On entry: the number of strata. If
${\mathbf{NS}}>0$ then the stratum for each observation must be supplied in
ISI.
Constraint:
${\mathbf{NS}}\ge 0$.
 5: Z(LDZ,M) – REAL (KIND=nag_wp) arrayInput
On entry: the
$i$th row must contain the covariates which are associated with the
$i$th failure time given in
T.
 6: LDZ – INTEGERInput
On entry: the first dimension of the array
Z as declared in the (sub)program from which G12BAF is called.
Constraint:
${\mathbf{LDZ}}\ge {\mathbf{N}}$.
 7: ISZ(M) – INTEGER arrayInput
On entry: indicates which subset of covariates is to be included in the model.
 ${\mathbf{ISZ}}\left(j\right)\ge 1$
 The $j$th covariate is included in the model.
 ${\mathbf{ISZ}}\left(j\right)=0$
 The $j$th covariate is excluded from the model and not referenced.
Constraint:
${\mathbf{ISZ}}\left(j\right)\ge 0$ and at least one and at most
${n}_{0}1$ elements of
ISZ must be nonzero where
${n}_{0}$ is the number of observations excluding any with zero value of
ISI.
 8: IP – INTEGERInput
On entry: the number of covariates included in the model as indicated by
ISZ.
Constraints:
 ${\mathbf{IP}}\ge 1$;
 ${\mathbf{IP}}=\text{ number of nonzero values of}{\mathbf{ISZ}}$.
 9: T(N) – REAL (KIND=nag_wp) arrayInput
On entry: the vector of $n$ failure censoring times.
 10: IC(N) – INTEGER arrayInput
On entry: the status of the individual at time
$t$ given in
T.
 ${\mathbf{IC}}\left(i\right)=0$
 The $i$th individual has failed at time ${\mathbf{T}}\left(i\right)$.
 ${\mathbf{IC}}\left(i\right)=1$
 The $i$th individual has been censored at time ${\mathbf{T}}\left(i\right)$.
Constraint:
${\mathbf{IC}}\left(\mathit{i}\right)=0$ or $1$, for $\mathit{i}=1,2,\dots ,{\mathbf{N}}$.
 11: OMEGA($*$) – REAL (KIND=nag_wp) arrayInput

Note: the dimension of the array
OMEGA
must be at least
${\mathbf{N}}$ if
${\mathbf{OFFSET}}=\text{'Y'}$, and at least
$1$ otherwise.
On entry: if
${\mathbf{OFFSET}}=\text{'Y'}$, the offset,
${\omega}_{\mathit{i}}$, for
$\mathit{i}=1,2,\dots ,{\mathbf{N}}$. Otherwise
OMEGA is not referenced.
 12: ISI($*$) – INTEGER arrayInput

Note: the dimension of the array
ISI
must be at least
${\mathbf{N}}$ if
${\mathbf{NS}}>0$, and at least
$1$ otherwise.
On entry: if
${\mathbf{NS}}>0$, the stratum indicators which also allow data points to be excluded from the analysis.
If
${\mathbf{NS}}=0$,
ISI is not referenced.
 ${\mathbf{ISI}}\left(i\right)=k$
 The $i$th data point is in the $k$th stratum, where $k=1,2,\dots ,{\mathbf{NS}}$.
 ${\mathbf{ISI}}\left(i\right)=0$
 The $i$th data point is omitted from the analysis.
Constraint:
if
${\mathbf{NS}}>0$,
$0\le {\mathbf{ISI}}\left(\mathit{i}\right)\le {\mathbf{NS}}$ and more than
IP values of
${\mathbf{ISI}}\left(\mathit{i}\right)>0$, for
$\mathit{i}=1,2,\dots ,{\mathbf{N}}$.
 13: DEV – REAL (KIND=nag_wp)Output
On exit: the deviance, that is $2\times \text{}$(maximized log marginal likelihood).
 14: B(IP) – REAL (KIND=nag_wp) arrayInput/Output
On entry: initial estimates of the covariate coefficient parameters
$\beta $.
${\mathbf{B}}\left(j\right)$ 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
${\mathbf{B}}\left(j\right)$ may be used. For other suggestions see
Section 8.
On exit:
${\mathbf{B}}\left(j\right)$ 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.
 15: SE(IP) – REAL (KIND=nag_wp) arrayOutput
On exit: ${\mathbf{SE}}\left(\mathit{j}\right)$ is the asymptotic standard error of the estimate contained in ${\mathbf{B}}\left(\mathit{j}\right)$ and score function in ${\mathbf{SC}}\left(\mathit{j}\right)$, for $\mathit{j}=1,2,\dots ,{\mathbf{IP}}$.
 16: SC(IP) – REAL (KIND=nag_wp) arrayOutput
On exit: ${\mathbf{SC}}\left(j\right)$ is the value of the score function, ${U}_{j}\left(\beta \right)$, for the estimate contained in ${\mathbf{B}}\left(j\right)$.
 17: COV(${\mathbf{IP}}\times \left({\mathbf{IP}}+1\right)/2$) – REAL (KIND=nag_wp) arrayOutput
On exit: the variancecovariance matrix of the parameter estimates in
B stored in packed form by column, i.e., the covariance between the parameter estimates given in
${\mathbf{B}}\left(i\right)$ and
${\mathbf{B}}\left(j\right)$,
$j\ge i$, is stored in
${\mathbf{COV}}\left(j\left(j1\right)/2+i\right)$.
 18: RES(N) – REAL (KIND=nag_wp) arrayOutput
On exit: the residuals,
$r\left({t}_{\mathit{l}}\right)$, for $\mathit{l}=1,2,\dots ,{\mathbf{N}}$.
 19: ND – INTEGEROutput
On exit: the number of distinct failure times.
 20: TP(NDMAX) – REAL (KIND=nag_wp) arrayOutput
On exit: ${\mathbf{TP}}\left(\mathit{i}\right)$ contains the $\mathit{i}$th distinct failure time, for $\mathit{i}=1,2,\dots ,{\mathbf{ND}}$.
 21: SUR(NDMAX,$*$) – REAL (KIND=nag_wp) arrayOutput

Note: the second dimension of the array
SUR
must be at least
$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{NS}},1\right)$.
On exit: if
${\mathbf{NS}}=0$,
${\mathbf{SUR}}\left(i,1\right)$ contains the estimated survival function for the
$i$th distinct failure time.
If ${\mathbf{NS}}>0$, ${\mathbf{SUR}}\left(i,k\right)$ contains the estimated survival function for the $i$th distinct failure time in the $k$th stratum.
 22: NDMAX – INTEGERInput
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:
${\mathbf{NDMAX}}\ge \text{the number of distinct failure times. This is returned in}{\mathbf{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 ${\mathbf{TOL}}\times \left(1.0+\mathrm{CurrentDeviance}\right)$. This corresponds approximately to an absolute precision if the deviance is small and a relative precision if the deviance is large.
Constraint:
${\mathbf{TOL}}\ge 10\times \mathit{machineprecision}$.
 24: MAXIT – INTEGERInput
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, variancecovariance matrix and the survival function are computed for the input value of
$\beta $ in
B but
$\beta $ is not updated.
Constraint:
${\mathbf{MAXIT}}\ge 0$.
 25: IPRINT – INTEGERInput
On entry: indicates if the printing of information on the iterations is required.
 ${\mathbf{IPRINT}}\le 0$
 No printing.
 ${\mathbf{IPRINT}}\ge 1$
 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: WK(${\mathbf{IP}}\times \left({\mathbf{IP}}+9\right)/2+{\mathbf{N}}$) – REAL (KIND=nag_wp) arrayWorkspace
 27: IWK($2\times {\mathbf{N}}$) – INTEGER arrayWorkspace
 28: IFAIL – INTEGERInput/Output

On entry:
IFAIL must be set to
$0$,
$1\text{ or}1$. If you are unfamiliar with this parameter you should refer to
Section 3.3 in the Essential Introduction for details.
For environments where it might be inappropriate to halt program execution when an error is detected, the value
$1\text{ or}1$ is recommended. If the output of error messages is undesirable, then the value
$1$ is recommended. Otherwise, if you are not familiar with this parameter, the recommended value is
$0$.
When the value $\mathbf{1}\text{ or}\mathbf{1}$ is used it is essential to test the value of IFAIL on exit.
On exit:
${\mathbf{IFAIL}}={\mathbf{0}}$ unless the routine detects an error or a warning has been flagged (see
Section 6).
6 Error Indicators and Warnings
If on entry
${\mathbf{IFAIL}}={\mathbf{0}}$ or
${{\mathbf{1}}}$, explanatory error messages are output on the current error message unit (as defined by
X04AAF).
Errors or warnings detected by the routine:
 ${\mathbf{IFAIL}}=1$
On entry,  ${\mathbf{OFFSET}}\ne \text{'Y'}$ or $\text{'N'}$, 
or  ${\mathbf{M}}<1$, 
or  ${\mathbf{N}}<2$, 
or  ${\mathbf{NS}}<0$, 
or  ${\mathbf{IP}}<1$, 
or  ${\mathbf{LDZ}}<{\mathbf{N}}$, 
or  ${\mathbf{TOL}}<10\times \mathit{machineprecision}$, 
or  ${\mathbf{MAXIT}}<0$. 
 ${\mathbf{IFAIL}}=2$
On entry,  ${\mathbf{ISZ}}\left(i\right)<0$ for some $i$, 
or  the value of IP is incompatible with ISZ, 
or  ${\mathbf{IC}}\left(i\right)\ne 1$ or $0$. 
or  ${\mathbf{ISI}}\left(i\right)<0$ or ${\mathbf{ISI}}\left(i\right)>{\mathbf{NS}}$, 
or  number of values of ${\mathbf{ISZ}}\left(i\right)>0$ is greater than or equal to ${n}_{0}$, the number of observations excluding any with ${\mathbf{ISI}}\left(i\right)=0$, 
or  all observations are censored, i.e., ${\mathbf{IC}}\left(i\right)=1$ for all $i$, 
or  NDMAX is too small. 
 ${\mathbf{IFAIL}}=3$
The matrix of second partial derivatives is singular. Try different starting values or include fewer covariates.
 ${\mathbf{IFAIL}}=4$
Overflow has been detected. Try using different starting values.
 ${\mathbf{IFAIL}}=5$
Convergence has not been achieved in
MAXIT iterations. The progress toward convergence can be examined by using a nonzero value of
IPRINT. Any nonconvergence may be due to a linear combination of covariates being monotonic with time.
Full results are returned.
 ${\mathbf{IFAIL}}=6$
In the current iteration $10$ step halvings have been performed without decreasing the deviance from the previous iteration. Convergence is assumed.
7 Accuracy
The accuracy is specified by
TOL.
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
$\mathrm{exp}\left({\beta}^{\mathrm{T}}{z}_{i}\right)$ or there may be nonconvergence. Reasonable estimates can often be obtained by fitting an exponential model using
G02GCF.
9 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.
9.1 Program Text
Program Text (g12bafe.f90)
9.2 Program Data
Program Data (g12bafe.d)
9.3 Program Results
Program Results (g12bafe.r)