COMPUTERS

Trend

AND

BIOMEDICAL

RESEARCH

10,373-381(1977)

and Homogeneity

Analyses of Proportions Table Data

DONALDG.THOMAS,*NORMANBRESLOW,~' *National

Cancer

Institute,

and Life

ANDJOHN J. GART*

Bethesda, Maryland 20014 and *Department Washington, Seattle, Washington 98195

ofBiostatktics.

Unioersify

of

Received August 18, 1976

A FORTRANprogram for computing various statistical tests relating proportions to dose is described. Life table curves are computed and printer plots may be obtained. Exact and approximate tests of linear trend for proportions are given. Finally two methods of testing for linear trend while adjusting for the time element are programmed. In both cases, tests of homogeneity and departures from trend are also given. The program is available from authors.

In someanimal experiments (such as in carcinogenesisbioassay) and clinical trials the data consist not only of counts of the number of events (animals with tumor or deaths of patients) but also of the time until the event occurs. An added complication is the question of censoredobservations, e.g., deathsbefore the end of the experiment or patients having incomplete follow-up. The usual feature of interest is the relation of theseresponsesto doselevel of the treatment. TABLE

Group

Dose

i

4

0

0

1 2

1.5 2.0

Initial group size (niO) 9 10 10

1

Times to event (I& and censored times (tkt) 73+, 74+. 75+, 76,76,76+,99, X66,246+ 43+,44+, 45+,67,68+, 136, 136, 150, 150, 150 41+, 41+, 47,47+, 47+, 58, 58,58, loo+, 117

Consider the fictitious data of Table 1, which will be considered the result of a carcinogenesisexperiment. The times to event may be considered the time in days until an event occurs (the tumor is observed) and those times with a plus sign are the times at which animals died without tumor in the course of the experiment. In order to take into account this censoring it is necessaryto restructure the data into the life table form given in Table 2. Here the number of events per animal at risk for the Copyright 0 1971 by Academic Press, Inc. All rights of reproduction in any form reserved. Printed in Great Britain.

373

ISSN 0010-4809

374

THOMAS,

BRESLOW,

AND

TABLE

GAR?

2

NOT,4l’ION

Distinct

Doses

times

k

‘A

I

1,

I,

1,

K

lh

_~~

d,,

Totals

.‘.

I,lk’ILi

4

.u,,ln,,,

t71,.1t1.,

.\‘,I ‘11:/

mi/n.l

-%ln,, s,]Jnoh

Crude proportions Effective proportions

x,./n,,,, .u,,.ln,,,

s,.in, x,./n,,

.x, * i n,*

m,ln.,

.\-, ./n ,O .l-, h,,

tn./n.,, min.,

various groups is given at each distinct event time. In Table 3 the explicit structure of the data of Table 1 is given in this form. This program implements the methods for analyzing both the unadjusted TABLE

RATIOS

OF NUMBERS

OF EVENTS DISTINCT

(xii) TO THE NUMBER TIMES OF EVENTS (t,) Group

Distinct

AT RISK

(nlk)

AT THE

numbers

times

k

1,

1 2

41 58

3 4

61 76 99 117 136 150 166

5 6 I 8 9”

Crude proportions Effective proportions (xi.lni,) 0 As there

3

(xi./nio)

are individuals

0 0 o/9 o/9 019 216

1 1.5 o/7

2 2.0

Totals m,ln.,

l/8 315

Ii24 312 I l/18 2113 l/IO l/8 2,fl 3!5 Ii2

o/2

o/7 117 o/5 o/5 o/5 215

o/2 112

3:3 o/o

O/2 o/2 l/l o/o o/o o/o

419

6110

5/10

15129

419

617

518

15124

l/3

o/2

o/2

at risk in only one group

at k = 9. the omission of this and r,), but will not affect the values of the other test statistics (see (IO)). This also applies to pairwise comparisons involving group 2 at k = 7 and k = 8.

time stratum will improve the conservative approximations (7’,

PROPORTIONS

AND LIFE TABLE

DATA

315

proportions and the life table adjusted data. Additional features include associated graphical output which aids in the interpretation of the results. Options are included for testing two kinds of unadjustedproportions: (1) the crude proportions, xi.lni,,, the number of events per number of individuals initially treated; and (2) the effective proportion, x,../nil, the number of events per number at risk when the first event occurs among any of the groups. Another option permits pairwise comparison of all the groups by all the statistical tests. It is also possibleto pool groups easily.

METHODS

1. Life Table Curves and Standard Errors Life table curves are computed for each of the groups by using time to event and adjusting for the censored observations. This is called the “Kaplan--Meier Survival Curve P(t)” in the output and is graphed as the “Proportion Surviving without Known Tumor.” A secondlife table analysis takes all the time points and treats them as deaths in the life table formulas. This is called the “Standard Survival Curve, S(t)” in the output and is graphed as the “Proportion Surviving.” The product limit estimatesof the proportion of individuals in the ith group who survive (say, without known tumor) past t, are calculated recursively from Kaplan and Meier (I.?),

where sik is the number in the ith group who survive past t, and Pi(&) = 1.. Its approximate variance is similarly estimated as

vIPi(tk)l = V{pi(tk-l)\ sfklnfk+ {P(t312XJ(nikSJ9

k=l,2

,..., K,

where V{16&) } = 0 and theseformulas are sequentially applied to tied observations. 2. Unadjusted Analysis of Proportions The logistic model (4) provides an elegant framework for developing an exact regressiontest of crude proportions. Let the unconditional probability, Pi, of an event in group i be functionally related to the dose, d, through the relation: Pi =

ea+p, 1 + ea+pd,’

i = 0, 1, 2, . . ., r,

(2.1)

where a and /3 are unknown parameters. We shall be interested in testing the hypothesis of homogeneity of the Pi usually against the alternative hypothesis of a positive linear trend in the logit scale, H,:p= 0 Hi:/3 < 0

or or

In {Pi/(1 -Pi)} = a, In {Pi/(1 -Pi)} = a+ /Id,.

376

THOMAS,

BRESLOW,

AND

CART

If the events are assumed to be identically and independently distributed within the groups (i.e., the time element dependence is ignored) then the likelihood can be written as a product of r + 1 binomial distributions, exp(ccx..+plxi.di) L=

! I& 11 + exp (a + pdi)}nid

(2.2)

When the effective proportion rather than the crude proportion is to be analyzed, lzi, replaces rzi,,, the initial number on test, in (2.2). This shows x.. and xxi.di to be a minimal set of sufficient statistics. An optimal test of H, can be based on the conditional distribution with the observed total number of events x.. fixed. This is the r-dimensional hypergeometric distribution, f(X.lx..

; n,) = I&

where the capital X. indicates the random variable corresponding is (t&j,, nio9. . ., n,J. The P value for the exact conditional test is P = lf(X. R

(2.3)

(~:)A3

to the xI.‘s and n,,

Ix. .; n,).

(2.4)

where R is the set of Y + 1 fold partitions {X.) ofx.. such that 1 Xi. dj 2 \‘ xi. d, the observed value. When the alternative is p < 0, the sense of this inequality is reversed. When r = 1, the test is equivalent to the Fisher-Irwin exact test. The corresponding asymptotic test for positive trend is based on the normal deviate:

z = z‘ Xi. d; - (X. * 1 niodi)/t~.,- 412 IV(~xi.diIx..;p=O)]“* ’

(2.5)

where

V(EXi. diIx..;p=o) I

=

x..(n.,-x..) Y n, (d, -- a,*, L n., (n.,- 1) ;

where d = 1 niodi/n., and A/2 is the continuity correction. If the doses are equally spaced then A is the interval between adjacent doses. When the doses are not equally spaced, no constant correction is adequate for all outcomes. In this case the program gives a “lower bound” where A is taken to be the minimum interval between adjacent doses and an “upper bound” where A is taken to be the maximum difference between adjacent values of the d;s. The numerator of (2.5), apart from the continuity correction, and V(li Xi. d,lx. .; /I = O), the square of the denominator, are also printed out. These can be used to combine similar experiments or clinical trials (for instance, over sex, laboratories, or hospitals) into a single normal deviate test (1.5).

PROPORTIONS

AND

LIFE

TABLE

311

DATA

One can also test whether (2.1) is an adequate model. This is most conveniently done by finding the homogeneity x2 with r degrees of freedom: n., (n., - 1) Xr2

= x..

1

(n.,-x..)

XT

X.. 2

*iO

n-0

‘-c

i

. >

The chi-square statistic for linear trend is the square of z, omitting correction, lCXi.diXT2Z

(X-a

the continuity

~niodi)/n.o)2

i v(~xrdiIx..;p=o)

.

The chi-square for departure from trend is (see (I)) x’,- I = xr2 - x;. When effective proportions above equations.

are to be used ni, is substituted for n,, throughout the

3. Adjusted Analyses Under the null hypotheses of equality among the r + 1 survival distributions, the vector X, = (Xlk, Xzk, . . ., X,,)’ considered conditionally on x.~ = mk and nk = (* ,k, . . .’ nrk)’ fixed will have an r-dimensional hypergeometric distribution,f(X,lx.,; nk, (see (2.3)). It follows that its mean, e, and covariance matrix, V,, have components i=l,2 ,..., r, eik = E(X,k) = nik(m!/n.k)y and nik(6ij - n,~kh’k~mk(n.k - mk) II VkII 0 = COV(xi@ x/k) = -3 n.k(n.k - 1) respectively, where 6, = 1 or 0 as i = j or i #j. Furthermore, following the partial likelihood argument of Cox (7), these K conditional distributions may be regarded formally as statistically independent. The adjusted analyses consist of various weighted combinations of the information in the K x (r + 1) x 2 tables (e.g., Tables 2 and 3). Cox (6) (see also (16)) suggests using the summary statistic, 0 - E, where 0 = C, xk and E = E,kek represent the vectors of observed and “expected” numbers’ in each of the last r groups. This has covariance matrix2 V = C,V, so that T, = (0 - E)’ V-’ (0 -E) I If a large proportion of the individuals in the study are followed to death, the “expected” can exceed the total initial number “at risk.” * Note that the computer output gives the full (r + 1) x (r + 1) covariance matrices.

numbers

318

THOMAS,

BRESLOW.

AND

GART

will have an approximate xrZ distribution in large samples. This is referred to as Cox’s test in the output. When r = 1, T, is the “log rank” statistic of Peto and Peto (18). This test has been shown to have certain optimal properties for testing survival curves having hazard functions which are multiples of one another. A conservative approximation to this statistic, not requiring matrix inversion, is T2 = 1 (Oi - E;)‘/Di, i-0

where

The conservatism should be minimal unless the proportions at risk niJn., vary markedly with time (8). If there are no ties among observed times of death, i.e., m, =-. 1 for all k, then Dj = E, and Tz reduces to the familiar y (O-E)2/E criterion (see (19)). The test for trend of the adjusted data was derived by Tarone (20) and is formally identical to the analogous contingency table test of Mantel (15). The appropriate one degree of freedom chi-square statistic is then3

T,=

i di (Oi - EJ i=o

; ;’ V.Ii{/d.d. -i=O j-.0

The departure from trend may be tested by T, = T, - T, which has an asymptotic chi-square distribution with r - 1 degrees of freedom. Utilizing concepts introduced by Cox (.5), Tarone (20) derived an exact test for trend under some further assumptions in the structure of the model. Introduction of a continuity correction will permit a better approximation to the exact test. The corrected asymptotic normal deviate is

z=

i di (Oi - Ei) - A/2

i=O

I,

1/2

Bounds on the continuity correction A/2 are handled as in (2.5). The numerator of the z statistic apart from the continuity correction and its variance, d’V d, are given so that they may be used in combining results. 3 Note that in the linear are used.

trend tests, the full r + 1 vectors

and (r + 1) Y (r + 1) covariance

matrices

PROPORTIONS

AND

LIFE

TABLE

379

DATA

The d’s entered into the trend statistic need not necessarily be actual doses, but may be coefficients for a linear contrast of interest. For instance if d, = -2, d, = d, = 1, then T3 will test in some sense whether the mean response rates of groups 1 and 2 differ from that of group 0. The generalized Wilcoxon or Kruskal-Wallis statistic, introduced for survival data by Gehan (II) and Breslow (2), may also be calculated from the basic data. This uses another scoring procedure which gives more weight to early deaths. Specifically, if

wik = n.pik - mknik is the score associated with the Ith group at time t,, and Wi = &IV,, is its total score then the r-dimensional vector4 W = (IV,, . . . , IV,)’ has, under the null hypothesis, mean 0 and covariance matrix v,=

; nZ,V,=

pyq.

k=l

Consequently T5 = W’ V,’ W may also be referred to tables of xr2 to test the global null hypothesis of no difference in the r + 1 survival distributions. Using the same arguments as previously (8), we have a conservative approximation to this statistic which requires no matrix inversion : T, = x WJGi, i=O

where

Gi = 5 {n.kmk(n.k - mk) nikl(n.k - 1)). k=l

Provided the proportions nik/n.k stay reasonably constant over time, T6 will be only slightly smaller than T5. A single d.f. x,* variate for trend based on W is

while may be referred to x:- 1 tables for testing departure from trend. In contrast to the situations with O-E, the tests based on the IV scores are not derivable by likelihood methods from any known model for the data and hence lack optimality properties. Nevertheless they are known to be more powerful than the OE procedures for testing against certain alternatives where the ratios of hazard functions are nonconstant (14). As already noted, W gives relatively more weight to early differences in the death rates for the r + 1 groups. 4 This definition of W treats ties between suggested in (2) rather than in (II).

censored

and uncensored

observations

in the manner

380

THOMAS,

BRESLOW,

AND

GART

4. Available Options One option is deleting early deaths, that is, using ni, in the unadjusted analyses. When these n’s and the T’S are large (say, n . 0 or n. , > 100 or I > 3) the exact analyses may be omitted, thus saving computer time with little loss of accuracy. There are various options for printing out either the full life table survival curves or only their median and final points. Several options for plotting the curves are available. Other manipulations of the data are possible. Groups may easily be pooled. An option is available for making all pairwise comparisons of the groups. This option produces the chi-square tests for all the unadjusted and adjusted analyses.

5. Discussion of the Methods Several authors have discussed the problems in applying to animal experimentation the methodology programmed in this paper. Hoe1 and Walburg (12) and Peto (17) have emphasized that uncensored observations should be entered by cause of death. Gart (IO) has questioned the feasibility of this in some feeding experiments. It should be noted that many of the methods of multiple contingency tables are formally equivalent to Cox’s formulation of life table regression problems (see (3, 9)). In particular Breslow pointed out how the ratio of hazard rates (i to j) may be estimated by Oi Ei/qi E, when this ratio is not too far from 1. In the formally equivalent 2 x 2 x k contingency table problem, the ratio of the hazard rates is estimated by the common odds ratio. Thomas (21) gives a comprehensive program for doing such related analyses.

REFERENCES 1. ARMITAGE, P. Tests for linear trends in proportions and frequencies. Biometrics II, 375-386 (1955). 2. BRESLOW, N. A generalized Kruskal-Wallis test for comparing K samples subject to unequal patterns of censorship. Biometriku 57,579-594 (1970). 3. BRESLOW, N. Analysis of survival data under the proportional hazards model. Inf. Stat. Rev. 43, 45-57 (1975). 4. Cox, D. R. The regression analysis of binary sequences (with discussion). .I. Roy. Staf. Sot. 20, 2 15-242 (1958). 5. Cox. D. R. The analysis of exponentially distributed life-times with two types of failure. J. Roy. Stat. Sot. Ser. B. 21,411-421 (1959). 6. Cox, D. R. Regression models and life tables (with discussion). J. Roy Stat. Sot. Ser. B. 34, 187.. 220 (1972). 7. Cox, D. R. Partial likelihood. Biomefrika 62,269-276 (1975). 8. CROWLEY, J. AND BRESLOW, N. Remarks on the conservatism of \(O-f?)‘/E in survival data. Biometrics 31,957-961 (1975). 9. GART, J. J. Contribution to the discussion on the paper by D. R. Cox, Regression Models and Life Tables. J. Roy. Statist. Sot. Ser. B. 34,212-213 (1972). 10. GART. J. J. Letter to the editor. Br. J. Cancer 31,696-697 (1975). Il. GEHAN, E. A. A generalized Wilcoxon test for comparing K samples subject to unequal patterns of censorship. Biometrika 52,203-223 (1965).

PROPORTIONS

AND

LIFE

TABLE

DATA

381

12. HOEL, D. G. AND WALBURG, H. E. Statistical analysis of survival experiments. J. Nat. Cancer Inst. 49,361-372 (1972). 13. KAPLAN, E. L. AND MEIER, P. Non-parametric estimation from incomplete observations. J. Am. Stat. Assoc. 53,457-481 (1958). 14. LEE, E. T., DESU, M. M., AND GEHAN, E. A. A Monte Carlo study of the power of some twosample tests. Biometrika 62,425-432 (1975). 15. MANTEL, N. Chi-square tests with one degree of freedom; extensions of the Mantel-Haenszel procedure. .I. Am. Stat. Assoc. 58,690-700 (1963). 16. MANTEL, N. AND HAENSZEL, W. Statistical aspects of the analysis of data from retrospective studies of disease. J. Nat. Cancer Inst. 22,719-748 (1959). 17. PETO, R. Guidelines on the analysis of tumor rates and death rates in experimental animals (editorial). Brit. J. Cancer 29, 101-105 (1974). 18. PETO, R. AND PETO, .I. Asympototically efficient rank invariant test procedures (with discussion). J. Roy. Stat. Sot. Ser. A. 135,185-206 (1972). 19. PETO, R. AND PIKE, M. C. Conservatism of the approximation ~(&E)*/E in the logrank test for survival data or tumor incidence data. Biometrics 29,579-584 (1973). 20. TARONE, R. E. Tests for trend in life table analysis. Biometrika 62,679-682 (1975). 21. THOMAS, D. G. Exact and asymptotic methods for the combination of 2 x 2 tables. Comput. Biomed. Res. 8,423-446 (1975).