hide long namesshow long names
hide short namesshow short names
Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

NAG Toolbox: nag_surviv_coxmodel (g12ba)

Purpose

nag_surviv_coxmodel (g12ba) returns parameter estimates and other statistics that are associated with the Cox proportional hazards model for fixed covariates.

Syntax

[dev, b, se, sc, cov, res, nd, tp, sur, ifail] = g12ba(offset, ns, z, isz, t, ic, omega, isi, b, ndmax, tol, maxit, iprint, 'n', n, 'm', m, 'ip', ip)
[dev, b, se, sc, cov, res, nd, tp, sur, ifail] = nag_surviv_coxmodel(offset, ns, z, isz, t, ic, omega, isi, b, ndmax, tol, maxit, iprint, 'n', n, 'm', m, 'ip', ip)

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 titi, for i = 1,2,,ni=1,2,,n, be the failure time or censored time for the iith observation with the vector of pp covariates zizi. It is assumed that censoring and failure mechanisms are independent. The hazard function, λ(t,z)λ(t,z), is the probability that an individual with covariates zz fails at time tt given that the individual survived up to time tt. In the Cox proportional hazards model (see Cox (1972)) λ(t,z)λ(t,z) is of the form:
λ(t,z) = λ0(t)exp(zTβ + ω)
λ(t,z)=λ0(t)exp(zTβ+ω)
where λ0λ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 < nnd<n distinct failure times, t(1) < < t(nd)t(1)<<t(nd) such that didi individuals fail at t(i)t(i), it follows that the marginal likelihood for ββ is well approximated (see Kalbfleisch and Prentice (1980)) by:
nd
L = (exp(siTβ + ωi))/([lR(t(i))exp(zlTβ + ωl)]di)
i = 1
L=i=1ndexp(siTβ+ωi) [lR(t(i))exp(zlTβ+ωl)]di
(1)
where sisi is the sum of the covariates of individuals observed to fail at t(i)t(i) and R(t(i))R(t(i)) is the set of individuals at risk just prior to t(i)t(i), that is, it is all individuals that fail or are censored at time t(i)t(i) along with all individuals that survive beyond time t(i)t(i). 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:
nd
Uj(β) = (lnL)/(βj) = [sjidiαji(β)] = 0
i = 1
Uj(β)= lnL βj =i=1nd[sji-diαji(β)]=0
(2)
for j = 1,2,,pj=1,2,,p, where sjisji is the jjth element in the vector sisi and
αji(β) = (lR(t(i))zjlexp(zlTβ + ωl))/(lR(t(i))exp(zlTβ + ωl)).
αji(β)=lR(t(i))zjlexp(zlTβ+ωl) lR(t(i))exp(zlTβ+ωl) .
Similarly,
nd
Ihj(β) = (2lnL)/(βhβj) = diγhji
i = 1
Ihj(β)=- 2lnL βhβj =i=1nddiγhji
(3)
where
γhji = (lR(t(i)) zhlzjlexp(zlTβ + ωl))/(lR(t(i))exp(zlTβ + ωl))αhi(β)αji(β),   h,j = 1,,p.
γhji=lR(t(i)) zhlzjlexp(zlTβ+ωl) lR(t(i))exp(zlTβ+ωl) -αhi(β)αji(β),   h,j= 1,,p.
Uj(β)Uj(β) is the jjth component of a score vector and Ihj(β)Ihj(β) is the (h,j)(h,j) element of the observed information matrix I(β)I(β) whose inverse I(β)1 = [Ihj(β)]1I(β)-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βj's will be infinite.
If λ0(t)λ0(t) varies across νν strata, where the number of individuals in the kkth stratum is nknk, for k = 1,2,,νk=1,2,,ν with n = k = 1νnkn=k=1νnk, then rather than maximizing (1) to obtain β̂β^, the following marginal likelihood is maximized:
ν
L = Lk,
k = 1
L=k=1νLk,
(4)
where LkLk is the contribution to likelihood for the nknk observations in the kkth 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λ0.
The base-line survivor function associated with a failure time t(i)t(i), is estimated as exp((t(i)))exp(-H^(t(i))), where
(t(i)) = t(j)t(i) ((di)/(lR(t(j))exp(zlTβ̂ + ωl))) ,
H^(t(i))=t(j)t(i) (dilR(t(j))exp(zlTβ^+ωl) ) ,
(5)
where didi is the number of failures at time t(i)t(i). The residual for the llth observation is computed as:
r(tl) = (tl)exp(zlTβ̂ + ωl)
r(tl)= H^(tl)exp(zlTβ^+ωl)
where (tl) = (t(i)),t(i)tl < t(i + 1)H^(tl)=H^(t(i)),t(i)tl<t(i+1). The deviance is defined as 2 × -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χ2-distribution; or, the asymptotic normality of the parameter estimates can be used to form zz 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 zz 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

Parameters

Compulsory Input Parameters

1:     offset – string (length ≥ 1)
Indicates if an offset is to be used.
offset = 'Y'offset='Y'
An offset must be included in omega.
offset = 'N'offset='N'
No offset is included in the model.
Constraint: offset = 'Y'offset='Y' or 'N''N'.
2:     ns – int64int32nag_int scalar
The number of strata. If ns > 0ns>0 then the stratum for each observation must be supplied in isi.
Constraint: ns0ns0.
3:     z(ldz,m) – double array
ldz, the first dimension of the array, must satisfy the constraint ldznldzn.
The iith row must contain the covariates which are associated with the iith failure time given in t.
4:     isz(m) – int64int32nag_int array
m, the dimension of the array, must satisfy the constraint m1m1.
Indicates which subset of covariates is to be included in the model.
isz(j)1iszj1
The jjth covariate is included in the model.
isz(j) = 0iszj=0
The jjth covariate is excluded from the model and not referenced.
Constraint: isz(j)0iszj0 and at least one and at most n01n0-1 elements of isz must be nonzero where n0n0 is the number of observations excluding any with zero value of isi.
5:     t(n) – double array
n, the dimension of the array, must satisfy the constraint n2n2.
The vector of nn failure censoring times.
6:     ic(n) – int64int32nag_int array
n, the dimension of the array, must satisfy the constraint n2n2.
The status of the individual at time tt given in t.
ic(i) = 0ici=0
The iith individual has failed at time t(i)ti.
ic(i) = 1ici=1
The iith individual has been censored at time t(i)ti.
Constraint: ic(i) = 0ici=0 or 11, for i = 1,2,,ni=1,2,,n.
7:     omega( : :) – double array
Note: the dimension of the array omega must be at least nn if offset = 'Y'offset='Y', and at least 11 otherwise.
If offset = 'Y'offset='Y', the offset, ωiωi, for i = 1,2,,ni=1,2,,n. Otherwise omega is not referenced.
8:     isi( : :) – int64int32nag_int array
Note: the dimension of the array isi must be at least nn if ns > 0ns>0, and at least 11 otherwise.
If ns > 0ns>0, the stratum indicators which also allow data points to be excluded from the analysis.
If ns = 0ns=0, isi is not referenced.
isi(i) = kisii=k
The iith data point is in the kkth stratum, where k = 1,2,,nsk=1,2,,ns.
isi(i) = 0isii=0
The iith data point is omitted from the analysis.
Constraint: if ns > 0ns>0, 0isi(i)ns0isiins and more than ip values of isi(i) > 0isii>0, for i = 1,2,,ni=1,2,,n.
9:     b(ip) – double array
ip, the dimension of the array, must satisfy the constraint
  • ip1ip1
  • ip = ​ number of nonzero values of ​iszip=​ number of nonzero values of ​isz
  • .
    Initial estimates of the covariate coefficient parameters ββ. b(j)bj must contain the initial estimate of the coefficient of the covariate in z corresponding to the jjth nonzero value of isz.
    10:   ndmax – int64int32nag_int scalar
    The dimension of the array tp and the first dimension of the array sur as declared in the (sub)program from which nag_surviv_coxmodel (g12ba) is called.
    Constraint: ndmaxthe number of distinct failure times. This is returned in ​ndndmaxthe number of distinct failure times. This is returned in ​nd.
    11:   tol – double scalar
    Indicates the accuracy required for the estimation. Convergence is assumed when the decrease in deviance is less than tol × (1.0 + CurrentDeviance)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 precisiontol10×machine precision.
    12:   maxit – int64int32nag_int scalar
    The maximum number of iterations to be used for computing the estimates. If maxit is set to 00 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: maxit0maxit0.
    13:   iprint – int64int32nag_int scalar
    Indicates if the printing of information on the iterations is required.
    iprint0iprint0
    No printing.
    iprint1iprint1
    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 nag_file_set_unit_advisory (x04ab)).

    Optional Input Parameters

    1:     n – int64int32nag_int scalar
    Default: The dimension of the arrays t, ic and the first dimension of the array z. (An error is raised if these dimensions are not equal.)
    nn, the number of data points.
    Constraint: n2n2.
    2:     m – int64int32nag_int scalar
    Default: The dimension of the array isz and the second dimension of the array z. (An error is raised if these dimensions are not equal.)
    The number of covariates in array z.
    Constraint: m1m1.
    3:     ip – int64int32nag_int scalar
    Default: The dimension of the array b.
    The number of covariates included in the model as indicated by isz.
    Constraints:
    • ip1ip1;
    • ip = ​ number of nonzero values of ​iszip=​ number of nonzero values of ​isz.

    Input Parameters Omitted from the MATLAB Interface

    ldz wk iwk

    Output Parameters

    1:     dev – double scalar
    The deviance, that is 2 × -2×(maximized log marginal likelihood).
    2:     b(ip) – double array
    b(j)bj contains the estimate β̂iβ^i, the coefficient of the covariate stored in the iith column of z where ii is the jjth nonzero value in the array isz.
    3:     se(ip) – double array
    se(j)sej is the asymptotic standard error of the estimate contained in b(j)bj and score function in sc(j)scj, for j = 1,2,,ipj=1,2,,ip.
    4:     sc(ip) – double array
    sc(j)scj is the value of the score function, Uj(β)Uj(β), for the estimate contained in b(j)bj.
    5:     cov(ip × (ip + 1) / 2ip×(ip+1)/2) – double array
    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)bi and b(j)bj, jiji, is stored in cov(j(j1) / 2 + i)cov(j(j-1)/2+i).
    6:     res(n) – double array
    The residuals, r(tl)r(tl), for l = 1,2,,nl=1,2,,n.
    7:     nd – int64int32nag_int scalar
    The number of distinct failure times.
    8:     tp(ndmax) – double array
    tp(i)tpi contains the iith distinct failure time, for i = 1,2,,ndi=1,2,,nd.
    9:     sur(ndmax, : :) – double array
    The second dimension of the array will be max (ns,1)max(ns,1)
    If ns = 0ns=0, sur(i,1)suri1 contains the estimated survival function for the iith distinct failure time.
    If ns > 0ns>0, sur(i,k)surik contains the estimated survival function for the iith distinct failure time in the kkth stratum.
    10:   ifail – int64int32nag_int scalar
    ifail = 0ifail=0 unless the function detects an error (see [Error Indicators and Warnings]).

    Error Indicators and Warnings

    Errors or warnings detected by the function:

    Cases prefixed with W are classified as warnings and do not generate an error of type NAG:error_n. See nag_issue_warnings.

      ifail = 1ifail=1
    On entry,offset'Y'offset'Y' or 'N''N',
    orm < 1m<1,
    orn < 2n<2,
    orns < 0ns<0,
    orip < 1ip<1,
    orldz < nldz<n,
    ortol < 10 × machine precisiontol<10×machine precision,
    ormaxit < 0maxit<0.
      ifail = 2ifail=2
    On entry,isz(i) < 0iszi<0 for some ii,
    orthe value of ip is incompatible with isz,
    oric(i)1ici1 or 00.
    orisi(i) < 0isii<0 or isi(i) > nsisii>ns,
    ornumber of values of isz(i) > 0iszi>0 is greater than or equal to n0n0, the number of observations excluding any with isi(i) = 0isii=0,
    orall observations are censored, i.e., ic(i) = 1ici=1 for all ii,
    orndmax is too small.
      ifail = 3ifail=3
    The matrix of second partial derivatives is singular. Try different starting values or include fewer covariates.
      ifail = 4ifail=4
    Overflow has been detected. Try using different starting values.
    W ifail = 5ifail=5
    Convergence has not been achieved in maxit 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.
    Full results are returned.
    W ifail = 6ifail=6
    In the current iteration 1010 step halvings have been performed without decreasing the deviance from the previous iteration. Convergence is assumed.

    Accuracy

    The accuracy is specified by tol.

    Further Comments

    nag_surviv_coxmodel (g12ba) 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)exp(βTzi) or there may be non-convergence. Reasonable estimates can often be obtained by fitting an exponential model using nag_correg_glm_poisson (g02gc).

    Example

    function nag_surviv_coxmodel_example
    offset = 'No-offset';
    ns = int64(0);
    z = [0;
         0;
         0;
         0;
         0;
         0;
         0;
         0;
         0;
         0;
         0;
         0;
         0;
         0;
         0;
         0;
         0;
         0;
         0;
         0;
         0;
         1;
         1;
         1;
         1;
         1;
         1;
         1;
         1;
         1;
         1;
         1;
         1;
         1;
         1;
         1;
         1;
         1;
         1;
         1;
         1;
         1];
    isz = [int64(1)];
    t = [1;
         1;
         2;
         2;
         3;
         4;
         4;
         5;
         5;
         8;
         8;
         8;
         8;
         11;
         11;
         12;
         12;
         15;
         17;
         22;
         23;
         6;
         6;
         6;
         7;
         10;
         13;
         16;
         22;
         23;
         6;
         9;
         10;
         11;
         17;
         19;
         20;
         25;
         32;
         32;
         34;
         35];
    ic = [int64(0);0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;1;1;1;1;1;1;1;1;1;1;1;1];
    omega = [0];
    isi = [int64(0)];
    b = [-1.526606964258046];
    ndmax = int64(42);
    tol = 5e-05;
    maxit = int64(20);
    iprint = int64(0);
    [dev, bOut, se, sc, covar, res, nd, tp, sur, ifail] = ...
        nag_surviv_coxmodel(offset, ns, z, isz, t, ic, omega, isi, b, ndmax, tol, maxit, iprint)
    
     
    
    dev =
    
      172.7592
    
    
    bOut =
    
       -1.5091
    
    
    se =
    
        0.4096
    
    
    sc =
    
       3.3775e-04
    
    
    covar =
    
        0.1677
    
    
    res =
    
        0.0780
        0.0780
        0.1626
        0.1626
        0.2088
        0.3057
        0.3057
        0.4130
        0.4130
        0.9141
        0.9141
        0.9141
        0.9141
        1.1864
        1.1864
        1.4175
        1.4175
        1.7233
        2.0993
        2.6630
        3.5226
        0.1312
        0.1312
        0.1312
        0.1452
        0.2216
        0.3466
        0.4217
        0.5888
        0.7789
        0.1312
        0.2021
        0.2216
        0.2623
        0.4642
        0.4642
        0.4642
        0.7789
        0.7789
        0.7789
        0.7789
        0.7789
    
    
    nd =
    
                       17
    
    
    tp =
    
         1
         2
         3
         4
         5
         6
         7
         8
        10
        11
        12
        13
        15
        16
        17
        22
        23
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
    
    
    sur =
    
        0.9640
        0.9264
        0.9065
        0.8661
        0.8235
        0.7566
        0.7343
        0.6506
        0.6241
        0.5724
        0.5135
        0.4784
        0.4447
        0.4078
        0.3727
        0.2859
        0.1908
             0
             0
             0
             0
             0
             0
             0
             0
             0
             0
             0
             0
             0
             0
             0
             0
             0
             0
             0
             0
             0
             0
             0
             0
             0
    
    
    ifail =
    
                        0
    
    
    
    function g12ba_example
    offset = 'No-offset';
    ns = int64(0);
    z = [0;
         0;
         0;
         0;
         0;
         0;
         0;
         0;
         0;
         0;
         0;
         0;
         0;
         0;
         0;
         0;
         0;
         0;
         0;
         0;
         0;
         1;
         1;
         1;
         1;
         1;
         1;
         1;
         1;
         1;
         1;
         1;
         1;
         1;
         1;
         1;
         1;
         1;
         1;
         1;
         1;
         1];
    isz = [int64(1)];
    t = [1;
         1;
         2;
         2;
         3;
         4;
         4;
         5;
         5;
         8;
         8;
         8;
         8;
         11;
         11;
         12;
         12;
         15;
         17;
         22;
         23;
         6;
         6;
         6;
         7;
         10;
         13;
         16;
         22;
         23;
         6;
         9;
         10;
         11;
         17;
         19;
         20;
         25;
         32;
         32;
         34;
         35];
    ic = [int64(0);0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;1;1;1;1;1;1;1;1;1;1;1;1];
    omega = [0];
    isi = [int64(0)];
    b = [-1.526606964258046];
    ndmax = int64(42);
    tol = 5e-05;
    maxit = int64(20);
    iprint = int64(0);
    [dev, bOut, se, sc, covar, res, nd, tp, sur, ifail] = ...
        g12ba(offset, ns, z, isz, t, ic, omega, isi, b, ndmax, tol, maxit, iprint)
    
     
    
    dev =
    
      172.7592
    
    
    bOut =
    
       -1.5091
    
    
    se =
    
        0.4096
    
    
    sc =
    
       3.3775e-04
    
    
    covar =
    
        0.1677
    
    
    res =
    
        0.0780
        0.0780
        0.1626
        0.1626
        0.2088
        0.3057
        0.3057
        0.4130
        0.4130
        0.9141
        0.9141
        0.9141
        0.9141
        1.1864
        1.1864
        1.4175
        1.4175
        1.7233
        2.0993
        2.6630
        3.5226
        0.1312
        0.1312
        0.1312
        0.1452
        0.2216
        0.3466
        0.4217
        0.5888
        0.7789
        0.1312
        0.2021
        0.2216
        0.2623
        0.4642
        0.4642
        0.4642
        0.7789
        0.7789
        0.7789
        0.7789
        0.7789
    
    
    nd =
    
                       17
    
    
    tp =
    
         1
         2
         3
         4
         5
         6
         7
         8
        10
        11
        12
        13
        15
        16
        17
        22
        23
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
    
    
    sur =
    
        0.9640
        0.9264
        0.9065
        0.8661
        0.8235
        0.7566
        0.7343
        0.6506
        0.6241
        0.5724
        0.5135
        0.4784
        0.4447
        0.4078
        0.3727
        0.2859
        0.1908
             0
             0
             0
             0
             0
             0
             0
             0
             0
             0
             0
             0
             0
             0
             0
             0
             0
             0
             0
             0
             0
             0
             0
             0
             0
    
    
    ifail =
    
                        0
    
    
    

    PDF version (NAG web site, 64-bit version, 64-bit version)
    Chapter Contents
    Chapter Introduction
    NAG Toolbox

    © The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2013